solpp.c (11565B)
1 /* Copyright (C) 2017, 2018, 2025 |Méso|Star> 2 * 3 * This program is free software: you can redistribute it and/or modify 4 * it under the terms of the GNU General Public License as published by 5 * the Free Software Foundation, either version 3 of the License, or 6 * (at your option) any later version. 7 * 8 * This program is distributed in the hope that it will be useful, 9 * but WITHOUT ANY WARRANTY; without even the implied warranty of 10 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 11 * GNU General Public License for more details. 12 * 13 * You should have received a copy of the GNU General Public License 14 * along with this program. If not, see <http://www.gnu.org/licenses/>. */ 15 16 #include "solpp.h" 17 18 struct mesh { 19 BUF(double) coords; 20 BUF(size_t) ids; 21 BUF(size_t) entities; 22 BUF(size_t) ncells; 23 size_t voffset; 24 }; 25 static const struct mesh MESH_NULL = {BUF_NULL, BUF_NULL, BUF_NULL, BUF_NULL, 0}; 26 27 static inline void 28 mesh_release(struct mesh* msh) 29 { 30 BUF_RELEASE(msh->coords); 31 BUF_RELEASE(msh->ids); 32 BUF_RELEASE(msh->entities); 33 BUF_RELEASE(msh->ncells); 34 memset(msh, 0, sizeof(struct mesh)); 35 } 36 37 static inline void 38 mesh_write_vtk(FILE* output, struct mesh* msh) 39 { 40 size_t i, n; 41 42 fprintf(output, "# vtk DataFile Version 2.0\n"); 43 fprintf(output, "Test\n"); 44 fprintf(output, "ASCII\n"); 45 fprintf(output, "DATASET POLYDATA\n"); 46 47 n = BUF_SZ(msh->coords)/3; 48 fprintf(output, "POINTS %zu double\n", n); 49 FOR_EACH(i, 0, n) { 50 fprintf(output, "%g %g %g\n", 51 BUF_AT(msh->coords, i*3+0), 52 BUF_AT(msh->coords, i*3+1), 53 BUF_AT(msh->coords, i*3+2)); 54 } 55 n = BUF_SZ(msh->ids)/3; 56 fprintf(output, "POLYGONS %zu %zu\n", n, n*4); 57 FOR_EACH(i, 0, n) { 58 fprintf(output, "3 %zu %zu %zu\n", 59 BUF_AT(msh->ids, i*3+0), 60 BUF_AT(msh->ids, i*3+1), 61 BUF_AT(msh->ids, i*3+2)); 62 } 63 } 64 65 static inline void 66 mesh_write_obj(FILE* output, struct mesh* msh) 67 { 68 size_t i, n; 69 70 n = BUF_SZ(msh->coords)/3; 71 FOR_EACH(i, 0, n) { 72 fprintf(output, "v %g %g %g\n", 73 BUF_AT(msh->coords, i*3+0), 74 BUF_AT(msh->coords, i*3+1), 75 BUF_AT(msh->coords, i*3+2)); 76 } 77 78 n = BUF_SZ(msh->ids)/3; 79 FOR_EACH(i, 0, n) { 80 fprintf(output, "f %zu %zu %zu\n", 81 BUF_AT(msh->ids, i*3+0) + 1, 82 BUF_AT(msh->ids, i*3+1) + 1, 83 BUF_AT(msh->ids, i*3+2) + 1); 84 } 85 } 86 87 static inline void 88 mesh_write_prim_data_vtk 89 (FILE* output, const struct mesh* msh, struct simul* simul) 90 { 91 struct prim* prim; 92 struct rcvXprim* rcvXprim; 93 size_t ircv; 94 size_t iprim; 95 size_t icell; 96 size_t n; 97 98 n = BUF_SZ(msh->ids) / 3; 99 fprintf(output, "CELL_DATA %zu\n", n); 100 101 fprintf(output, "FIELD PrimaryData %zu\n", 2 + BUF_SZ(simul->rcvs)*6); 102 fprintf(output, "cos_factor 2 %zu double\n", n); 103 FOR_EACH(iprim, 0, BUF_SZ(msh->entities)) { 104 CHK(prim = find_primary_by_id(simul, BUF_AT(msh->entities, iprim))); 105 FOR_EACH(icell, 0, BUF_AT(msh->ncells, iprim)) { 106 fprintf(output, "%g %g\n", prim->cos_factor.E, prim->cos_factor.SE); 107 } 108 } 109 110 fprintf(output, "shadow_loss 2 %zu double\n", n); 111 FOR_EACH(iprim, 0, BUF_SZ(msh->entities)) { 112 CHK(prim = find_primary_by_id(simul, BUF_AT(msh->entities, iprim))); 113 FOR_EACH(icell, 0, BUF_AT(msh->ncells, iprim)) { 114 fprintf(output, "%g %g\n", prim->shadow_loss.E, prim->shadow_loss.SE); 115 } 116 } 117 118 #define WRITE(Side, Name) { \ 119 fprintf(output, "%s_"STR(Side)"_"STR(Name)" 2 %zu double\n", rcv->name, n);\ 120 FOR_EACH(iprim, 0, BUF_SZ(msh->entities)) { \ 121 CHK(rcvXprim = find_rcvXprim(simul,rcv->id,BUF_AT(msh->entities,iprim)));\ 122 FOR_EACH(icell, 0, BUF_AT(msh->ncells, iprim)) { \ 123 fprintf(output, "%g %g\n", \ 124 rcvXprim->Name[Side].E, \ 125 rcvXprim->Name[Side].SE); \ 126 } \ 127 } \ 128 } (void)0 129 FOR_EACH(ircv, 0, BUF_SZ(simul->rcvs)) { 130 const struct rcv* rcv = &BUF_AT(simul->rcvs, ircv); 131 WRITE(FRONT, in.flux); 132 WRITE(FRONT, in.flux_mat_loss); 133 WRITE(FRONT, in.flux_atm_loss); 134 WRITE(BACK, in.flux); 135 WRITE(BACK, in.flux_mat_loss); 136 WRITE(BACK, in.flux_atm_loss); 137 } 138 #undef WRITE 139 } 140 141 static inline void 142 mesh_write_rcv_data_vtk 143 (FILE* output, const struct mesh* msh, struct simul* simul) 144 { 145 struct rcv* rcv; 146 size_t ircv; 147 size_t icell; 148 size_t n; 149 150 n = BUF_SZ(msh->ids)/3; 151 fprintf(output, "CELL_DATA %zu\n", n); 152 fprintf(output, "FIELD PrimaryData 12\n"); 153 154 #define WRITE(Side, Name) { \ 155 fprintf(output, STR(Side)"_"STR(Name)" 2 %zu double\n", n); \ 156 FOR_EACH(ircv, 0, BUF_SZ(msh->entities)) { \ 157 CHK(rcv = find_receiver_by_id(simul, BUF_AT(msh->entities, ircv))); \ 158 FOR_EACH(icell, 0, BUF_AT(msh->ncells, ircv)) { \ 159 fprintf(output, "%g %g\n", rcv->Name[Side].E, rcv->Name[Side].SE); \ 160 } \ 161 } \ 162 } (void)0 163 WRITE(FRONT, in.flux); 164 WRITE(FRONT, in.flux_mat_loss); 165 WRITE(FRONT, in.flux_atm_loss); 166 WRITE(FRONT, efficiency); 167 WRITE(BACK, in.flux); 168 WRITE(BACK, in.flux_mat_loss); 169 WRITE(BACK, in.flux_atm_loss); 170 WRITE(BACK, efficiency); 171 #undef WRITE 172 173 #define WRITE_MAP(Flux, Side) { \ 174 fprintf(output, STR(Side)"_"STR(Flux)"_map 2 %zu double\n", n); \ 175 FOR_EACH(ircv, 0, BUF_SZ(msh->entities)) { \ 176 CHK(rcv = find_receiver_by_id(simul, BUF_AT(msh->entities, ircv))); \ 177 if(!BUF_SZ(rcv->map[Flux][Side])) { \ 178 FOR_EACH(icell, 0, BUF_AT(msh->ncells,ircv)) fprintf(output,"-1 -1\n");\ 179 } else { \ 180 FOR_EACH(icell, 0, BUF_AT(msh->ncells,ircv)) { \ 181 const struct mc* mc = &BUF_AT(rcv->map[Flux][Side], icell); \ 182 fprintf(output, "%g %g\n", mc->E, mc->SE); \ 183 } \ 184 } \ 185 } \ 186 } (void)0 187 WRITE_MAP(INCOMING, FRONT); 188 WRITE_MAP(INCOMING, BACK); 189 WRITE_MAP(ABSORBED, FRONT); 190 WRITE_MAP(ABSORBED, BACK); 191 #undef WRITE_MAP 192 } 193 194 int 195 main(int argc, char** argv) 196 { 197 buf_char_T buf = BUF_NULL; 198 char* line = NULL; 199 200 FILE* input; 201 FILE* geom; 202 FILE* output; 203 204 if(argc < 3) { 205 fprintf(stderr, "Usage: %s solstice-geometry solstice-simulation\n", argv[0]); 206 return 1; 207 } 208 209 CHK(geom = fopen(argv[1], "r")); 210 CHK(input = fopen(argv[2], "r")); 211 212 CHK(read_line(&buf, geom)); /* Skip the sun direction of the geometry */ 213 while((line = read_line(&buf, input))) { 214 struct simul simul; 215 216 const struct rcv* rcv = NULL; 217 const struct prim* prim = NULL; 218 219 struct mesh msh_rcv = MESH_NULL; 220 struct mesh msh_prim = MESH_NULL; 221 struct mesh msh_misc = MESH_NULL; 222 223 char filename[128]; 224 char prefix[64]; 225 size_t ntris_grp = 0, nverts_grp = 0; 226 size_t nverts_obj = 0, off_obj = 0, off_grp = 0; 227 228 simul_init(&simul); 229 230 CHK(!strncmp(line, "#--- Sun direction:", 19)); 231 CHK(sscanf(line+19, "%lf %lf (%*f %*f %*f)", &simul.azimuth, &simul.elevation)==2); 232 CHK(snprintf(prefix, sizeof(prefix), "%g-%g-", 233 simul.azimuth, simul.elevation) < sizeof(prefix)); 234 read_simulation(&simul, input); 235 236 while((line = read_line(&buf, geom)) && strncmp(line, "#--- Sun", 8)) { 237 if(!strncmp(line, "g ", 2)) { 238 if(prim) { 239 BUF_PUSH(msh_prim.ncells, ntris_grp); 240 msh_prim.voffset += nverts_grp; 241 } 242 if(rcv) { 243 BUF_PUSH(msh_rcv.ncells, ntris_grp); 244 msh_rcv.voffset += nverts_grp; 245 } 246 if(!prim && !rcv) msh_misc.voffset += nverts_grp; 247 prim = find_primary(&simul, line+2); 248 rcv = find_receiver(&simul, line+2); 249 if(prim) BUF_PUSH(msh_prim.entities, prim->id); 250 if(rcv) BUF_PUSH(msh_rcv.entities, rcv->id); 251 ntris_grp = 0; 252 nverts_grp = 0; 253 off_grp = 0; 254 off_obj = nverts_obj; 255 256 } else if(!strncmp(line, "v ", 2)) { 257 double pos[3]; 258 CHK(sscanf(line+2, "%lf %lf %lf", pos+0, pos+1, pos+2) == 3); 259 if(prim) { 260 BUF_PUSH(msh_prim.coords, pos[0]); 261 BUF_PUSH(msh_prim.coords, pos[1]); 262 BUF_PUSH(msh_prim.coords, pos[2]); 263 } 264 if(rcv) { 265 BUF_PUSH(msh_rcv.coords, pos[0]); 266 BUF_PUSH(msh_rcv.coords, pos[1]); 267 BUF_PUSH(msh_rcv.coords, pos[2]); 268 } 269 if(!prim && !rcv) { 270 BUF_PUSH(msh_misc.coords, pos[0]); 271 BUF_PUSH(msh_misc.coords, pos[1]); 272 BUF_PUSH(msh_misc.coords, pos[2]); 273 } 274 ++nverts_grp; 275 ++nverts_obj; 276 } else if(!strncmp(line, "f ", 2)) { 277 size_t tri[3]; 278 CHK(sscanf(line+2, "%zu %zu %zu", tri+0, tri+1, tri+2) == 3); 279 if(prim) { 280 BUF_PUSH(msh_prim.ids, tri[0]-1 + msh_prim.voffset + off_grp - off_obj); 281 BUF_PUSH(msh_prim.ids, tri[1]-1 + msh_prim.voffset + off_grp - off_obj); 282 BUF_PUSH(msh_prim.ids, tri[2]-1 + msh_prim.voffset + off_grp - off_obj); 283 } 284 if(rcv) { 285 BUF_PUSH(msh_rcv.ids, tri[0]-1 + msh_rcv.voffset + off_grp - off_obj); 286 BUF_PUSH(msh_rcv.ids, tri[1]-1 + msh_rcv.voffset + off_grp - off_obj); 287 BUF_PUSH(msh_rcv.ids, tri[2]-1 + msh_rcv.voffset + off_grp - off_obj); 288 } 289 if(!prim && !rcv) { 290 BUF_PUSH(msh_misc.ids, tri[0]-1 + msh_misc.voffset + off_grp - off_obj); 291 BUF_PUSH(msh_misc.ids, tri[1]-1 + msh_misc.voffset + off_grp - off_obj); 292 BUF_PUSH(msh_misc.ids, tri[2]-1 + msh_misc.voffset + off_grp - off_obj); 293 } 294 ++ntris_grp; 295 } else if(!strcmp(line, "---")) { 296 nverts_obj = 0; 297 off_obj = 0; 298 off_grp = nverts_grp; 299 } 300 } 301 if(prim) BUF_PUSH(msh_prim.ncells, ntris_grp); 302 if(rcv) BUF_PUSH(msh_rcv.ncells, ntris_grp); 303 304 CHK(snprintf(filename, sizeof(filename), 305 "%sprimaries.vtk", prefix) < sizeof(prefix)); 306 printf("Writing `%s'\n", filename); 307 CHK(output = fopen(filename, "w")); 308 mesh_write_vtk(output, &msh_prim); 309 mesh_write_prim_data_vtk(output, &msh_prim, &simul); 310 fclose(output); 311 312 CHK(snprintf(filename, sizeof(filename), 313 "%sreceivers.vtk", prefix) < sizeof(prefix)); 314 printf("Writing `%s'\n", filename); 315 CHK(output = fopen(filename, "w")); 316 mesh_write_vtk(output, &msh_rcv); 317 mesh_write_rcv_data_vtk(output, &msh_rcv, &simul); 318 fclose(output); 319 320 CHK(snprintf(filename, sizeof(filename), 321 "%smiscellaneous.obj", prefix) < sizeof(prefix)); 322 printf("Writing `%s'\n", filename); 323 CHK(output = fopen(filename, "w")); 324 mesh_write_obj(output, &msh_misc); 325 fclose(output); 326 327 mesh_release(&msh_prim); 328 mesh_release(&msh_rcv); 329 mesh_release(&msh_misc); 330 simul_release(&simul); 331 } 332 333 fclose(geom); 334 fclose(input); 335 BUF_RELEASE(buf); 336 return 0; 337 } 338