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