solpp.c (12301B)
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 buf_char_T dont_read_i = BUF_NULL; 200 buf_char_T dont_read_g = BUF_NULL; 201 FILE* input; 202 FILE* geom; 203 FILE* output; 204 205 if(argc < 3) { 206 fprintf(stderr, "Usage: %s solstice-geometry solstice-simulation\n", argv[0]); 207 return 1; 208 } 209 210 CHK(geom = fopen(argv[1], "r")); 211 CHK(input = fopen(argv[2], "r")); 212 213 CHK(read_line(&buf, geom, &dont_read_g)); /* Skip the sun direction of the geometry */ 214 while((line = read_line(&buf, input, &dont_read_i))) { 215 struct simul simul; 216 217 const struct rcv* rcv = NULL; 218 const struct prim* prim = NULL; 219 220 struct mesh msh_rcv = MESH_NULL; 221 struct mesh msh_prim = MESH_NULL; 222 struct mesh msh_misc = MESH_NULL; 223 224 char filename[128]; 225 char prefix[64]; 226 size_t ntris_grp = 0, nverts_grp = 0; 227 size_t nverts_obj = 0, off_obj = 0, off_grp = 0; 228 229 simul_init(&simul); 230 231 if(!strncmp(line, "#--- Sun direction:", 19)) { 232 CHK(sscanf(line+19, "%lf %lf (%*f %*f %*f)", 233 &simul.azimuth, &simul.elevation)==2); 234 CHK(snprintf(prefix, sizeof(prefix), "%g-%g-", 235 simul.azimuth, simul.elevation) 236 < (long)sizeof(prefix)); 237 } else if(!strncmp(line, "#--- Sun location and time:", 27)) { 238 memset(simul.time, 0, sizeof(simul.time)); 239 CHK(sscanf(line+27, "%lf %lf %s (%*f %*f %*f)", 240 &simul.latitude, &simul.longitude, simul.time) == 3); 241 CHK(snprintf(prefix, sizeof(prefix), "%g-%g-%s-", 242 simul.latitude, simul.longitude, simul.time) 243 < (long)sizeof(prefix)); 244 } 245 read_simulation(&simul, input, &dont_read_i); 246 247 while((line = read_line(&buf, geom, &dont_read_g)) && strncmp(line, "#--- Sun", 8)) { 248 if(!strncmp(line, "g ", 2)) { 249 if(prim) { 250 BUF_PUSH(msh_prim.ncells, ntris_grp); 251 msh_prim.voffset += nverts_grp; 252 } 253 if(rcv) { 254 BUF_PUSH(msh_rcv.ncells, ntris_grp); 255 msh_rcv.voffset += nverts_grp; 256 } 257 if(!prim && !rcv) msh_misc.voffset += nverts_grp; 258 prim = find_primary(&simul, line+2); 259 rcv = find_receiver(&simul, line+2); 260 if(prim) { BUF_PUSH(msh_prim.entities, prim->id); } 261 if(rcv) { BUF_PUSH(msh_rcv.entities, rcv->id); } 262 ntris_grp = 0; 263 nverts_grp = 0; 264 off_grp = 0; 265 off_obj = nverts_obj; 266 267 } else if(!strncmp(line, "v ", 2)) { 268 double pos[3]; 269 CHK(sscanf(line+2, "%lf %lf %lf", pos+0, pos+1, pos+2) == 3); 270 if(prim) { 271 BUF_PUSH(msh_prim.coords, pos[0]); 272 BUF_PUSH(msh_prim.coords, pos[1]); 273 BUF_PUSH(msh_prim.coords, pos[2]); 274 } 275 if(rcv) { 276 BUF_PUSH(msh_rcv.coords, pos[0]); 277 BUF_PUSH(msh_rcv.coords, pos[1]); 278 BUF_PUSH(msh_rcv.coords, pos[2]); 279 } 280 if(!prim && !rcv) { 281 BUF_PUSH(msh_misc.coords, pos[0]); 282 BUF_PUSH(msh_misc.coords, pos[1]); 283 BUF_PUSH(msh_misc.coords, pos[2]); 284 } 285 ++nverts_grp; 286 ++nverts_obj; 287 } else if(!strncmp(line, "f ", 2)) { 288 size_t tri[3]; 289 CHK(sscanf(line+2, "%zu %zu %zu", tri+0, tri+1, tri+2) == 3); 290 if(prim) { 291 BUF_PUSH(msh_prim.ids, tri[0]-1 + msh_prim.voffset + off_grp - off_obj); 292 BUF_PUSH(msh_prim.ids, tri[1]-1 + msh_prim.voffset + off_grp - off_obj); 293 BUF_PUSH(msh_prim.ids, tri[2]-1 + msh_prim.voffset + off_grp - off_obj); 294 } 295 if(rcv) { 296 BUF_PUSH(msh_rcv.ids, tri[0]-1 + msh_rcv.voffset + off_grp - off_obj); 297 BUF_PUSH(msh_rcv.ids, tri[1]-1 + msh_rcv.voffset + off_grp - off_obj); 298 BUF_PUSH(msh_rcv.ids, tri[2]-1 + msh_rcv.voffset + off_grp - off_obj); 299 } 300 if(!prim && !rcv) { 301 BUF_PUSH(msh_misc.ids, tri[0]-1 + msh_misc.voffset + off_grp - off_obj); 302 BUF_PUSH(msh_misc.ids, tri[1]-1 + msh_misc.voffset + off_grp - off_obj); 303 BUF_PUSH(msh_misc.ids, tri[2]-1 + msh_misc.voffset + off_grp - off_obj); 304 } 305 ++ntris_grp; 306 } else if(!strcmp(line, "---")) { 307 nverts_obj = 0; 308 off_obj = 0; 309 off_grp = nverts_grp; 310 } 311 } 312 if(prim) { BUF_PUSH(msh_prim.ncells, ntris_grp); } 313 if(rcv) { BUF_PUSH(msh_rcv.ncells, ntris_grp); } 314 315 CHK(snprintf(filename, sizeof(filename), "%sprimaries.vtk", prefix) 316 < (long)sizeof(filename)); 317 printf("Writing `%s'\n", filename); 318 CHK(output = fopen(filename, "w")); 319 mesh_write_vtk(output, &msh_prim); 320 mesh_write_prim_data_vtk(output, &msh_prim, &simul); 321 fclose(output); 322 323 CHK(snprintf(filename, sizeof(filename), "%sreceivers.vtk", prefix) 324 < (long)sizeof(filename)); 325 printf("Writing `%s'\n", filename); 326 CHK(output = fopen(filename, "w")); 327 mesh_write_vtk(output, &msh_rcv); 328 mesh_write_rcv_data_vtk(output, &msh_rcv, &simul); 329 fclose(output); 330 331 CHK(snprintf(filename, sizeof(filename), "%smiscellaneous.obj", prefix) 332 < (long)sizeof(filename)); 333 if(BUF_SZ(msh_misc.coords) == 0) { 334 printf("No miscellanous geometry to write `%s'\n", filename); 335 } else { 336 printf("Writing `%s'\n", filename); 337 CHK(output = fopen(filename, "w")); 338 mesh_write_obj(output, &msh_misc); 339 fclose(output); 340 } 341 342 mesh_release(&msh_prim); 343 mesh_release(&msh_rcv); 344 mesh_release(&msh_misc); 345 simul_release(&simul); 346 } 347 348 fclose(geom); 349 fclose(input); 350 BUF_RELEASE(buf); 351 return 0; 352 }