solpp.h (14199B)
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 #define _POSIX_C_SOURCE 200809L /* Support of strdup */ 17 #include <string.h> 18 #include <stdio.h> 19 #include <stdlib.h> 20 21 #ifndef SOLPP_H 22 #define SOLPP_H 23 24 enum side { FRONT, BACK }; 25 enum flux_density { ABSORBED, INCOMING }; 26 struct mc { double E/* Expected value */, SE/* Standard Error */; }; 27 28 /******************************************************************************* 29 * Helper macros 30 ******************************************************************************/ 31 #define STR__(X) #X 32 #define STR(X) STR__(X) 33 34 #define FOR_EACH(Id, Start, End) for((Id)=(Start); (Id)<(End); ++(Id)) 35 36 #define CHK(Cond) { \ 37 if(!(Cond)) { \ 38 fprintf(stderr, "error:%s:%d\n", __FILE__, __LINE__); \ 39 abort(); \ 40 } \ 41 } (void)0 42 43 /******************************************************************************* 44 * Dynamic buffer 45 ******************************************************************************/ 46 #define BUF(Type) struct { Type* mem; size_t ca; size_t sz; } 47 #define BUF_NULL {NULL, 0, 0} 48 #define BUF_RELEASE(B) free((B).mem) 49 #define BUF_RESERVE(B, Sz) { \ 50 if((Sz)>(B).ca) { \ 51 CHK((B).mem = realloc((B).mem, sizeof(*(B).mem)*((B).ca=(Sz)))); \ 52 } \ 53 }(void)0 54 #define BUF_RESIZE(B, Sz) {BUF_RESERVE((B), Sz); (B).sz = Sz;} (void)0 55 #define BUF_PUSH(B, E) { \ 56 if((B).sz >= (B).ca) BUF_RESERVE((B), ((B).ca?(B).ca*2:32)); \ 57 (B).mem[(B).sz++]=(E); \ 58 } (void)0 59 #define BUF_SZ(B) (B).sz 60 #define BUF_MEM(B) (B).mem 61 #define BUF_AT(B, I) (B).mem[I] 62 63 /******************************************************************************* 64 * Helper function 65 ******************************************************************************/ 66 typedef BUF(char) buf_char_T; 67 68 static inline char* 69 read_line(buf_char_T* b, FILE* stream) 70 { 71 char* c; 72 73 if(!BUF_SZ(*b)) BUF_RESIZE(*b, 32); 74 if(!fgets(BUF_MEM(*b), (int)BUF_SZ(*b), stream)) return NULL; 75 76 /* Ensure that the whole line is read */ 77 while(!strrchr(BUF_MEM(*b), '\n') && !feof(stream)) { 78 BUF_RESIZE(*b, BUF_SZ(*b)+32); 79 CHK(fgets 80 (BUF_MEM(*b) + strlen(BUF_MEM(*b)), 81 (int)(BUF_SZ(*b) - strlen(BUF_MEM(*b))), stream)); 82 } 83 84 /* Remove the carriage return */ 85 if((c = strrchr(BUF_MEM(*b), '\n'))) *c = '\0'; 86 87 return BUF_MEM(*b); 88 } 89 90 /******************************************************************************* 91 * Per receiver double sided estimations 92 ******************************************************************************/ 93 struct flux { 94 struct mc flux[2]; 95 struct mc flux_no_mat_loss[2]; 96 struct mc flux_no_atm_loss[2]; 97 struct mc flux_mat_loss[2]; 98 struct mc flux_atm_loss[2]; 99 }; 100 101 struct rcv { 102 struct flux in; 103 struct flux abs; 104 struct mc efficiency[2]; 105 BUF(struct mc) map[2][2]; /* Absorbed/Incoming, Front/Back */ 106 char* name; 107 size_t id; 108 double area; 109 }; 110 111 static inline void rcv_init(struct rcv* r) {memset(r, 0, sizeof(*r));} 112 static inline void rcv_release(struct rcv* r) 113 { 114 free(r->name); 115 BUF_RELEASE(r->map[ABSORBED][FRONT]); 116 BUF_RELEASE(r->map[ABSORBED][BACK]); 117 BUF_RELEASE(r->map[INCOMING][FRONT]); 118 BUF_RELEASE(r->map[INCOMING][BACK]); 119 } 120 121 /******************************************************************************* 122 * Per primary estimations 123 ******************************************************************************/ 124 struct prim { 125 struct mc cos_factor; 126 struct mc shadow_loss; 127 char* name; 128 size_t id; 129 size_t nsamps; 130 double area; 131 }; 132 static inline void prim_init(struct prim* p) {memset(p, 0, sizeof(*p));} 133 static inline void prim_release(struct prim* p) {free(p->name);} 134 135 /******************************************************************************* 136 * Per receiver X primary estimations 137 ******************************************************************************/ 138 struct rcvXprim { 139 struct flux in; 140 struct flux abs; 141 size_t rcv_id; 142 size_t prim_id; 143 }; 144 static inline void rcvXprim_init(struct rcvXprim* r) {memset(r, 0, sizeof(*r));} 145 static inline void rcvXprim_release(struct rcvXprim* r) { /*Do nothing*/ } 146 147 /******************************************************************************* 148 * Overall estimations of a simulation 149 ******************************************************************************/ 150 struct simul { 151 struct mc potential_flux; 152 struct mc absorbed_flux; 153 struct mc cos_factor; 154 struct mc shadow_loss; 155 struct mc missing_loss; 156 struct mc materials_loss; 157 struct mc atmospheric_loss; 158 159 BUF(struct rcv) rcvs; 160 BUF(struct prim) prims; 161 BUF(struct rcvXprim) rcvXprims; 162 163 double azimuth; 164 double elevation; 165 size_t nsamps; 166 }; 167 168 static inline void simul_init(struct simul* s) {memset(s, 0, sizeof(*s));} 169 170 static inline void 171 simul_release(struct simul* s) 172 { 173 size_t i; 174 FOR_EACH(i,0,BUF_SZ(s->rcvs)) rcv_release(&BUF_AT(s->rcvs, i)); 175 FOR_EACH(i,0,BUF_SZ(s->prims)) prim_release(&BUF_AT(s->prims, i)); 176 FOR_EACH(i,0,BUF_SZ(s->rcvXprims)) rcvXprim_release(&BUF_AT(s->rcvXprims, i)); 177 BUF_RELEASE(s->rcvs); 178 BUF_RELEASE(s->prims); 179 BUF_RELEASE(s->rcvXprims); 180 } 181 182 /******************************************************************************* 183 * Look for entities 184 ******************************************************************************/ 185 static inline struct prim* 186 find_primary(struct simul* s, const char* name) 187 { 188 size_t i; 189 FOR_EACH(i, 0, BUF_SZ(s->prims)) 190 if(!strcmp(BUF_AT(s->prims, i).name, name)) return &BUF_AT(s->prims, i); 191 return NULL; 192 } 193 194 static inline struct prim* 195 find_primary_by_id(struct simul* s, const size_t id) 196 { 197 size_t i; 198 FOR_EACH(i, 0, BUF_SZ(s->prims)) 199 if(BUF_AT(s->prims, i).id == id) return &BUF_AT(s->prims, i); 200 return NULL; 201 } 202 203 static inline struct rcv* 204 find_receiver(struct simul* s, const char* name) 205 { 206 size_t i; 207 FOR_EACH(i, 0, BUF_SZ(s->rcvs)) 208 if(!strcmp(BUF_AT(s->rcvs, i).name, name)) return &BUF_AT(s->rcvs, i); 209 return NULL; 210 } 211 212 static inline struct rcv* 213 find_receiver_by_id(struct simul* s, const size_t id) 214 { 215 size_t i; 216 FOR_EACH(i, 0, BUF_SZ(s->rcvs)) 217 if(BUF_AT(s->rcvs, i).id == id) return &BUF_AT(s->rcvs, i); 218 return NULL; 219 } 220 221 static inline struct rcvXprim* 222 find_rcvXprim(struct simul* s, const size_t rcv_id, const size_t prim_id) 223 { 224 size_t i; 225 FOR_EACH(i, 0, BUF_SZ(s->rcvXprims)) { 226 struct rcvXprim* rXp = &BUF_AT(s->rcvXprims, i); 227 if(rXp->rcv_id == rcv_id && rXp->prim_id == prim_id) return rXp; 228 } 229 return NULL; 230 } 231 232 /******************************************************************************* 233 * Read simulation data from a solstice-output 234 ******************************************************************************/ 235 static inline void 236 read_receiver_map_side_data(struct rcv* rcv, const size_t n, FILE* input) 237 { 238 buf_char_T buf = BUF_NULL; 239 char* line = NULL; 240 size_t i; 241 enum side side; 242 enum flux_density flux; 243 const char* str; 244 245 CHK(line = read_line(&buf, input)); 246 str = line + 8; 247 if(!strncmp(str, "Front_faces", 11)) { side = FRONT; str += 12; } 248 else if(!strncmp(str, "Back_faces", 10)) { side = BACK; str += 11; } 249 else { fprintf(stderr, "Unexpected side name\n"); abort(); } 250 251 if(!strncmp(str, "Incoming_flux", 13)) { flux = INCOMING; } 252 else if(!strncmp(str, "Absorbed_flux", 13)) { flux = ABSORBED; } 253 else { fprintf(stderr, "Unexpected flux name\n"); abort(); } 254 255 CHK(read_line(&buf, input)); /* Discard "LOOKUP_TABLE default" line */ 256 257 BUF_RESIZE(rcv->map[flux][side], n); 258 FOR_EACH(i, 0, n) { 259 struct mc* mc = &BUF_AT(rcv->map[flux][side], i); 260 CHK(line = read_line(&buf, input)); 261 CHK(sscanf(line, "%lf %lf", &mc->E, &mc->SE) == 2); 262 } 263 264 BUF_RELEASE(buf); 265 } 266 267 static inline void 268 read_receiver_map(struct simul* simul, FILE* input) 269 { 270 struct rcv* rcv; 271 buf_char_T buf = BUF_NULL; 272 char* line = NULL; 273 size_t i, n; 274 long fp; 275 276 CHK(line = read_line(&buf, input)); 277 CHK(rcv = find_receiver(simul, line)); 278 279 /* Skip header */ 280 CHK(read_line(&buf, input)); 281 CHK(read_line(&buf, input)); 282 /* Skip vertices */ 283 CHK(line = read_line(&buf, input)); 284 CHK(sscanf(line, "POINTS %zu float", &n) == 1); 285 FOR_EACH(i, 0, n) CHK(read_line(&buf, input)); 286 /* Skip polygons */ 287 CHK(line = read_line(&buf, input)); 288 CHK(sscanf(line, "POLYGONS %zu %*u", &n) == 1); 289 FOR_EACH(i, 0, n) CHK(read_line(&buf, input)); 290 /* Read the map data of one side */ 291 CHK(line = read_line(&buf, input)); 292 CHK(sscanf(line, "CELL_DATA %zu", &n) == 1); 293 /* Read map data */ 294 do { 295 read_receiver_map_side_data(rcv, n, input); 296 fp = ftell(input); 297 line = read_line(&buf, input); 298 fseek(input, fp, SEEK_SET); 299 } while(line && !strncmp(line, "SCALARS", 7)); 300 301 BUF_RELEASE(buf); 302 } 303 304 static inline void 305 read_simulation(struct simul* simul, FILE* input) 306 { 307 buf_char_T buf = BUF_NULL; 308 char* line = NULL; 309 char* tk = NULL; 310 size_t nrcvs, nprims; 311 size_t i; 312 313 /* Counters */ 314 CHK(line = read_line(&buf, input)); 315 CHK(sscanf(line, "%*u %zu %zu %zu %*u", &nrcvs, &nprims, &simul->nsamps)==3); 316 317 /* Global results */ 318 #define READ(Name) { \ 319 CHK(line = read_line(&buf, input)); \ 320 CHK(sscanf(line, "%lf %lf", &simul->Name.E, &simul->Name.SE) == 2); \ 321 } (void)0 322 READ(potential_flux); 323 READ(absorbed_flux); 324 READ(cos_factor); 325 READ(shadow_loss); 326 READ(missing_loss); 327 READ(materials_loss); 328 READ(atmospheric_loss); 329 #undef READ 330 331 /* Read per receiver results */ 332 BUF_RESIZE(simul->rcvs, nrcvs); 333 FOR_EACH(i, 0, nrcvs) { 334 struct rcv* rcv = &BUF_AT(simul->rcvs, i); 335 rcv_init(rcv); 336 337 CHK(line = read_line(&buf, input)); 338 CHK(tk = strtok(line, " \t")); 339 CHK(rcv->name = strdup(tk)); 340 341 CHK(tk = strtok(NULL, "")); 342 #define GET(Side, Name) &rcv->Name[Side].E, &rcv->Name[Side].SE 343 CHK(sscanf 344 (tk, 345 "%zu %lf " 346 "%lf %lf %lf %lf %lf %lf %lf %lf %lf %lf " 347 "%lf %lf %lf %lf %lf %lf %lf %lf %lf %lf " 348 "%lf %lf " 349 "%lf %lf %lf %lf %lf %lf %lf %lf %lf %lf " 350 "%lf %lf %lf %lf %lf %lf %lf %lf %lf %lf " 351 "%lf %lf", 352 &rcv->id, &rcv->area, 353 GET(FRONT, in.flux), 354 GET(FRONT, in.flux_no_mat_loss), 355 GET(FRONT, in.flux_no_atm_loss), 356 GET(FRONT, in.flux_mat_loss), 357 GET(FRONT, in.flux_atm_loss), 358 GET(FRONT, abs.flux), 359 GET(FRONT, abs.flux_no_mat_loss), 360 GET(FRONT, abs.flux_no_atm_loss), 361 GET(FRONT, abs.flux_mat_loss), 362 GET(FRONT, abs.flux_atm_loss), 363 GET(FRONT, efficiency), 364 GET(BACK, in.flux), 365 GET(BACK, in.flux_no_mat_loss), 366 GET(BACK, in.flux_no_atm_loss), 367 GET(BACK, in.flux_mat_loss), 368 GET(BACK, in.flux_atm_loss), 369 GET(BACK, abs.flux), 370 GET(BACK, abs.flux_no_mat_loss), 371 GET(BACK, abs.flux_no_atm_loss), 372 GET(BACK, abs.flux_mat_loss), 373 GET(BACK, abs.flux_atm_loss), 374 GET(BACK, efficiency)) == 46); 375 #undef GET 376 } 377 378 /* Read per primary results */ 379 BUF_RESIZE(simul->prims, nprims); 380 FOR_EACH(i, 0, nprims) { 381 struct prim* prim = &BUF_AT(simul->prims, i); 382 prim_init(prim); 383 384 CHK(line = read_line(&buf, input)); 385 CHK(tk = strtok(line, " \t")); 386 CHK(prim->name = strdup(tk)); 387 388 CHK(tk = strtok(NULL, "")); 389 CHK(sscanf(tk, "%zu %lf %zu %lf %lf %lf %lf", 390 &prim->id, &prim->area, &prim->nsamps, 391 &prim->cos_factor.E, &prim->cos_factor.SE, 392 &prim->shadow_loss.E, &prim->shadow_loss.SE) == 7); 393 } 394 395 /* Per receiverXprimary results */ 396 BUF_RESIZE(simul->rcvXprims, nprims*nrcvs); 397 FOR_EACH(i, 0, nprims*nrcvs) { 398 struct rcvXprim* rcvXprim = &BUF_AT(simul->rcvXprims, i); 399 rcvXprim_init(rcvXprim); 400 401 CHK(line = read_line(&buf, input)); 402 #define GET(Side, Name) &rcvXprim->Name[Side].E, &rcvXprim->Name[Side].SE 403 CHK(sscanf 404 (line, 405 "%zu %zu " 406 "%lf %lf %lf %lf %lf %lf %lf %lf %lf %lf " 407 "%lf %lf %lf %lf %lf %lf %lf %lf %lf %lf " 408 "%lf %lf %lf %lf %lf %lf %lf %lf %lf %lf " 409 "%lf %lf %lf %lf %lf %lf %lf %lf %lf %lf", 410 &rcvXprim->rcv_id, &rcvXprim->prim_id, 411 GET(FRONT, in.flux), 412 GET(FRONT, in.flux_no_mat_loss), 413 GET(FRONT, in.flux_no_atm_loss), 414 GET(FRONT, in.flux_mat_loss), 415 GET(FRONT, in.flux_atm_loss), 416 GET(FRONT, abs.flux), 417 GET(FRONT, abs.flux_no_mat_loss), 418 GET(FRONT, abs.flux_no_atm_loss), 419 GET(FRONT, abs.flux_mat_loss), 420 GET(FRONT, abs.flux_atm_loss), 421 GET(BACK, in.flux), 422 GET(BACK, in.flux_no_mat_loss), 423 GET(BACK, in.flux_no_atm_loss), 424 GET(BACK, in.flux_mat_loss), 425 GET(BACK, in.flux_atm_loss), 426 GET(BACK, abs.flux), 427 GET(BACK, abs.flux_no_mat_loss), 428 GET(BACK, abs.flux_no_atm_loss), 429 GET(BACK, abs.flux_mat_loss), 430 GET(BACK, abs.flux_atm_loss)) == 42); 431 #undef GET 432 } 433 434 /* Read receiver maps */ 435 for(;;) { 436 const long fp = ftell(input); 437 line = read_line(&buf, input); 438 439 if(line && !strncmp(line, "# vtk", 5)) { 440 read_receiver_map(simul, input); 441 } else { 442 fseek(input, fp, SEEK_SET); 443 break; 444 } 445 } 446 BUF_RELEASE(buf); 447 } 448 449 #endif /* SOLPP_H */ 450