stardis

Perform coupled heat transfer calculations
git clone git://git.meso-star.fr/stardis.git
Log | Files | Refs | README | LICENSE

commit 8f577c6286759d0405e81e3eb9217d5ebd36c155
parent 07105c9726f8146f972f6435bc67e1d1ac27cd0f
Author: Christophe Coustet <christophe.coustet@meso-star.com>
Date:   Tue,  3 Mar 2020 18:01:17 +0100

Major refactoring

Change model description and parsing, now using star-geometry-3d

Stop using data structure without out-of-memory management

Diffstat:
Msrc/main.c | 60+++++++++++++++++++++++++++++++++++++-----------------------
Msrc/stardis-app.c | 645++++++++++++++++++++++++-------------------------------------------------------
Msrc/stardis-app.h | 895++++++++++++++++++++++++++++---------------------------------------------------
Msrc/stardis-compute.c | 659+++++++++++++++++++++++++++++++++++++++++--------------------------------------
Msrc/stardis-compute.h | 14+++++++++++++-
Msrc/stardis-fluid.c | 42++++++++++++++++++++++--------------------
Msrc/stardis-fluid.h | 11+++++++----
Msrc/stardis-intface.c | 96++++++++++++++++++++++++++++++++++++++++----------------------------------------
Msrc/stardis-intface.h | 8+++++---
Msrc/stardis-output.c | 575++++++++++++++++++++++++++++++++++---------------------------------------------
Msrc/stardis-output.h | 29++++++++++++++---------------
Msrc/stardis-parsing.c | 1239+++++++++++++++++++++++++++++++++++++++++++++++--------------------------------
Msrc/stardis-parsing.h | 86++++++++++++++++++++++++++++++++++++++++++++++++++++++++-----------------------
Msrc/stardis-solid.c | 42++++++++++++++++++++++--------------------
Msrc/stardis-solid.h | 10++++++----
15 files changed, 2076 insertions(+), 2335 deletions(-)

diff --git a/src/main.c b/src/main.c @@ -5,31 +5,29 @@ #include "stardis-output.h" #include <rsys/rsys.h> +#include <rsys/mem_allocator.h> #include <stdio.h> -#if !defined(ARCH_64BITS) -#error "Please write hash functions for the current arch" -#endif - int main (int argc, char** argv) { - struct args args = ARGS_DEFAULT; - struct stardis stardis = NULL_STARDIS; + struct args* args = NULL; + struct stardis stardis; + int allocator_initialized = 0, stardis_initialized = 0; int err = 0; + struct mem_allocator allocator; res_T res = RES_OK; - /* Check proper init for ARGS_DEFAULT */ - ASSERT(ARGS_DEFAULT.u.probe[3] == INF); - - res = parse_args(argc, argv, &args); - if (res != RES_OK) goto error; - if (args.just_help) goto exit; + ERR(mem_init_proxy_allocator(&allocator, &mem_default_allocator)); + allocator_initialized = 1; + ERR(init_args(&allocator, &args)); + ERR(parse_args(argc, argv, args)); + if(args->just_help) goto exit; - if (args.mode == UNDEF_MODE) { + if(args->mode == UNDEF_MODE) { fprintf(stderr, "Nothing to do.\n" "One of the following options should be used:"); print_multiple_modes(stderr, EXCLUSIVE_MODES); @@ -38,22 +36,38 @@ main goto exit; } - res = stardis_init(&args, &stardis); - if (res != RES_OK) goto error; + ERR(stardis_init(args, &stardis)); + stardis_initialized = 1; - if (args.mode & DUMP_VTK) { - dump_vtk(&stardis.geometry); + if(args->mode & DUMP_VTK) { + ERR(sg3d_geometry_dump_as_vtk(stardis.geometry.sg3d, stdout)); + /* If it applies dump the compute region */ + if(args->mode & REGION_COMPUTE_MODES) + ERR(dump_compute_region_at_the_end_of_vtk(&stardis, stdout)); + /* If it applies dump holes */ + if(has_holes(&stardis)) + ERR(dump_holes_at_the_end_of_vtk(&stardis, stdout)); + goto exit; /* If dump flag set, then dump and exit */ } - else if (args.mode & COMPUTE_MODES) { - res = stardis_compute(&stardis, args.mode); - if (res != RES_OK) goto error; + else if(args->mode & COMPUTE_MODES) { + ERR(stardis_compute(&stardis, args->mode)); } else { - fprintf(stderr, "Unknown mode: %c.\n", mode_option(args.mode)); + fprintf(stderr, "Unknown mode: %c.\n", mode_option(args->mode)); } exit: - stardis_release(&stardis); - CHK(mem_allocated_size() == 0); + release_args(args); + if(stardis_initialized) stardis_release(&stardis); + if(allocator_initialized) { + if(MEM_ALLOCATED_SIZE(&allocator) != 0) { + char dump[4096] = { '\0' }; + MEM_DUMP(&allocator, dump, sizeof(dump)); + fprintf(stderr, "%s\n", dump); + fprintf(stderr, "\nMemory leaks: %zu Bytes\n", + MEM_ALLOCATED_SIZE(&allocator)); + } + mem_shutdown_proxy_allocator(&allocator); + } return err; error: err = -1; diff --git a/src/stardis-app.c b/src/stardis-app.c @@ -2,16 +2,9 @@ #include "stardis-app.h" #include "stardis-output.h" +#include "stardis-compute.h" -#include <rsys/hash_table.h> -#define HTABLE_NAME descriptions -#define HTABLE_KEY struct description -#define HTABLE_DATA unsigned -#define HTABLE_KEY_FUNCTOR_EQ eq_description -#define HTABLE_KEY_FUNCTOR_HASH hash_description -#define HTABLE_KEY_FUNCTOR_INIT init_description -#define HTABLE_KEY_FUNCTOR_COPY cp_description -#include <rsys/hash_table.h> +#include <rsys/text_reader.h> #include <string.h> @@ -19,414 +12,120 @@ * Local Functions ******************************************************************************/ -static char* -read_line - (char** pb, - FILE* stream) -{ - const int chunk_sz = 32; - char* buf = *pb; - size_t idx; - - ASSERT(pb && stream); - - if(!buf) buf = sa_add(buf, 128); - - do { - /* Read stream until EOL, EOF or up to sa_size(b) */ - if(!fgets(buf, (int)sa_size(buf) + 1, stream)) return NULL; - - /* Ensure that he whole line is read */ - while(!strrchr(buf, '\n') && !feof(stream)) { - /* Add some room to buffer */ - char* more = sa_add(buf, (size_t)chunk_sz); - /* Continue reading chunk_sz more characters until EOL or EOF */ - fgets(more, chunk_sz + 1, stream); - } - /* Search idx of first occurence of comment char or EOL */ - idx = strcspn(buf, "#\n\r"); - /* Remove meaningless text */ - buf[idx] = '\0'; - /* Find first text */ - idx = strspn(buf, " \t"); - } while(buf[idx] == '\0'); /* Skip empty lines */ - *pb = buf; - return buf; -} - -static res_T -read_vertices - (struct sstl_desc* desc, - struct htable_vertex* vertex2id, - unsigned** id2id, - struct geometry* geometry) -{ - res_T res = RES_OK; - unsigned vtx_index = 0; - - ASSERT(desc && vertex2id && id2id && geometry); - - for (vtx_index = 0; vtx_index < desc->vertices_count; ++vtx_index){ - struct vertex vtx = NULL_VERTEX; - unsigned *found_id = NULL; - - vtx.xyz[0] = desc->vertices[3*vtx_index + 0]; - vtx.xyz[1] = desc->vertices[3*vtx_index + 1]; - vtx.xyz[2] = desc->vertices[3*vtx_index + 2]; - found_id = htable_vertex_find(vertex2id, &vtx); - - if (found_id){ - sa_push(*id2id, *found_id); - } else { - size_t sz = htable_vertex_size_get(vertex2id); - unsigned size; - ASSERT(sz < UINT_MAX); - size = (unsigned)sz; - res = htable_vertex_set(vertex2id, &vtx, &size); - if (res != RES_OK) return res; - sa_push(geometry->vertex, vtx); - sa_push(*id2id, size); - } - } - - return res; -} - -static int -indices_to_key - (struct unsigned3* key, - const struct unsigned3* indices) -{ - int i, reversed = 0; - - ASSERT(key && indices); - - FOR_EACH(i, 0, 3) key->data[i] = indices->data[i]; - if (key->data[0] > key->data[1]) { - reversed = !reversed; - SWAP(unsigned, key->data[0], key->data[1]); - } - if (key->data[1] > key->data[2]) { - reversed = !reversed; - SWAP(unsigned, key->data[1], key->data[2]); - } - if (key->data[0] > key->data[1]) { - reversed = !reversed; - SWAP(unsigned, key->data[0], key->data[1]); - } - return reversed; -} - static res_T -read_triangles - (struct sstl_desc* stl_desc, - struct htable_triangle* triangle2id, - const char* stl_filename, - const unsigned* id2id, - const enum content_type file_type, - const unsigned desc_id, - struct geometry* geometry, - struct solve_boundary* boundary) +read_model + (const struct darray_str* model_files, + struct stardis* stardis) { res_T res = RES_OK; - unsigned* packed_indices = NULL; - float* packed_normals = NULL; - unsigned tri_index = 0, drop_cpt = 0; - - ASSERT(stl_desc && triangle2id && stl_filename && id2id && geometry && boundary); - - for (tri_index = 0; tri_index < stl_desc->triangles_count; ++tri_index) { - struct triangle tri = NULL_TRIANGLE; - const unsigned *found_id = NULL; - unsigned i, id; - struct unsigned3 key; - int reversed; - int drop = 0; - enum sdis_side current_side = SDIS_FRONT; - unsigned* side_desc_ptr = NULL; - float coord[3][3]; - - FOR_EACH(i, 0, 3) - tri.indices.data[i] = id2id[stl_desc->indices[3*tri_index + i]]; - reversed = indices_to_key(&key, &tri.indices); - - FOR_EACH(i, 0, 3) { - unsigned j; - f3_set(coord[i], geometry->vertex[key.data[i]].xyz); - FOR_EACH(j, 0, i) { - if (f3_eq(coord[i], coord[j])) { - drop = 1; - drop_cpt++; - } - } - } - if (drop) continue; - FOR_EACH(i, 0, 3) { - sa_push(packed_indices, stl_desc->indices[3 * tri_index + i]); - sa_push(packed_normals, stl_desc->normals[3 * tri_index + i]); - } - - /* Find triangle if already known, or insert it in known list */ - found_id = htable_triangle_find(triangle2id, &key); - if (found_id) { - struct triangle* original_tri; - int original_reversed; - id = *found_id; - original_tri = geometry->triangle + id; - original_reversed = indices_to_key(&key, &original_tri->indices); - current_side = (reversed == original_reversed) ? SDIS_FRONT : SDIS_BACK; - ASSERT(id < sa_size(geometry->triangle)); - } else { - size_t sz; - if (file_type == CONTENT_BOUNDARY) { - /* Cannot create new geometry! */ - fprintf(stderr, - "Triangle %u from file %s cannot compute temperature on this surface\n" - "Not part of the model\n", - tri_index, stl_filename); - return RES_BAD_ARG; - } - sz = htable_triangle_size_get(triangle2id); - ASSERT(sz == sa_size(geometry->triangle)); - ASSERT(sz < UINT_MAX); - id = (unsigned)sz; - /* Store tri as described by stl_desc regardless of key order */ - res = htable_triangle_set(triangle2id, &key, &id); - if (res != RES_OK) goto error; - sa_push(geometry->triangle, tri); - } - - switch (file_type) { - case CONTENT_CONNECTION: - side_desc_ptr = &geometry->triangle[id].connection_description; - break; - case CONTENT_GEOMETRY: - /* Is the triangle in triangle list reversed wrt stl_desc? */ - side_desc_ptr = (current_side == SDIS_FRONT) ? - &geometry->triangle[id].front_description - : &geometry->triangle[id].back_description; - break; - case CONTENT_BOUNDARY: - sa_push(boundary->primitives, id); - sa_push(boundary->sides, current_side); - continue; /* Next triangle */ - break; - default: - FATAL("error:" STR(__FILE__) ":" STR(__LINE__)": Invalid type.\n"); + const struct str* files = NULL; + size_t i; + FILE* f = NULL; + struct txtrdr* txtrdr = NULL; + + ASSERT(model_files && stardis); + files = darray_str_cdata_get(model_files); + FOR_EACH(i, 0, darray_str_size_get(model_files)) { + const char* name = str_cget(files + i); + f = fopen(name, "r"); + if(!f) { + fprintf(stderr, "Cannot open model file '%s'\n", name); + res = RES_IO_ERR; + goto error; } - if (*side_desc_ptr != UINT_MAX && *side_desc_ptr != desc_id) { - /* Already described with a different description! */ - fprintf(stderr, "Triangle %u %s from file %s with 2 different descriptions\n", - tri_index, - (file_type == CONTENT_CONNECTION - ? "connection" : (current_side == SDIS_FRONT ? "front" : "back")), - stl_filename); - print_trg_as_obj(stderr, (const struct vertex*)stl_desc->vertices, - stl_desc->indices + (3 * tri_index)); - return RES_BAD_ARG; + txtrdr_stream(stardis->allocator, f, name, '#', &txtrdr); + for(;;) { + char* line; + ERR(txtrdr_read_line(txtrdr)); + line = txtrdr_get_line(txtrdr); + if(!line) break; + ERR(process_model_line(line, stardis)); } - /* Everithing is OK: store description */ - *side_desc_ptr = desc_id; + txtrdr_ref_put(txtrdr); + txtrdr = NULL; + fclose(f); + f = NULL; } - ASSERT(drop_cpt == 0 || packed_indices != NULL); - if (drop_cpt) { - fprintf(stderr, "File '%s' includes %u degenerate triangles (dropped)\n", - stl_filename, drop_cpt); - /* Keep only valid triangles */ - ASSERT(sa_size(packed_indices) % 3 == 0); - ASSERT(sa_size(packed_indices) == sa_size(packed_normals)); - ASSERT(sa_size(packed_indices) + 3 * drop_cpt == sa_size(stl_desc->indices)); - /* Don't sa_release stl_desc->indices or stl_desc->normals as they have - * their own release mechanism */ - stl_desc->indices = packed_indices; - stl_desc->normals = packed_normals; - stl_desc->triangles_count = sa_size(stl_desc->indices) / 3; - } else { - sa_release(packed_indices); - sa_release(packed_normals); - packed_indices = NULL; - packed_normals = NULL; - } - -end: - sa_release(packed_indices); - sa_release(packed_normals); + ASSERT(!f && !txtrdr); +exit: return res; - error: - goto end; +error: + if(f) fclose(f); + if(txtrdr) txtrdr_ref_put(txtrdr); + goto exit; } -static res_T -read_stl - (const unsigned desc_id, - const enum content_type file_type, - const char* stl_filename, - struct htable_vertex* vertex2id, - struct htable_triangle* triangle2id, - struct stardis* stardis) -{ - res_T res = RES_OK; - struct sstl* sstl = NULL; - struct sstl_desc stl_desc; - unsigned* id2id = NULL; - - ASSERT(stl_filename && vertex2id && triangle2id && stardis); - - res = sstl_create(NULL, &stardis->allocator, 0, &sstl); - if (res != RES_OK) goto error; - - res = sstl_load(sstl, stl_filename); - if (res != RES_OK) { - if (res == RES_IO_ERR) - fprintf(stderr, "Cannot read STL file: %s\n", stl_filename); - goto error; - } - SSTL(get_desc(sstl, &stl_desc)); - - res = read_vertices(&stl_desc, vertex2id, &id2id, &stardis->geometry); - if (res != RES_OK) goto error; - res = read_triangles(&stl_desc, triangle2id, stl_filename, id2id, file_type, - desc_id, &stardis->geometry, &stardis->boundary); - if (res != RES_OK) goto error; - -end: - sa_release(id2id); - SSTL(ref_put(sstl)); - return res; -error: - goto end; +#define COUNT_SIDE(Rank) {\ + if(properties[(Rank)] == SG3D_UNSPECIFIED_PROPERTY) u_count++;\ + else {\ + ASSERT(properties[(Rank)] < darray_descriptions_size_get(descriptions));\ + ASSERT(descs[properties[(Rank)]].type == DESC_MAT_SOLID\ + || descs[properties[(Rank)]].type == DESC_MAT_FLUID);\ + if(descs[properties[(Rank)]].type == DESC_MAT_SOLID) s_count++;\ + else f_count++;\ + }\ } static res_T -geometry_analyse - (const char* medium_filename, - const char* bc_filename, - const char* bound_filename, /* Can be NULL */ - struct stardis* stardis) +validate_properties + (const unsigned itri, + const unsigned properties[SG3D_PROP_TYPES_COUNT__], + void* context, + int* properties_conflict_status) { - res_T res = RES_OK; - FILE* input = NULL; - char* line = NULL; - unsigned sz; - struct htable_vertex vertex2id; - struct htable_triangle triangle2id; - struct htable_descriptions descriptions; - int tables_initialised = 0; - unsigned desc_id; - unsigned *p_desc; - - ASSERT(medium_filename && bc_filename && stardis); - - sz = (unsigned)sa_size(stardis->descriptions); - ASSERT(sz == sa_size(stardis->descriptions)); - - stardis->geometry = NULL_GEOMETRY; - - /* parse medium files */ - input = fopen(medium_filename,"r"); - if (!input){ - fprintf(stderr,"Cannot open %s\n", medium_filename); - res = RES_IO_ERR; - goto error; - } - - htable_vertex_init(&stardis->allocator, &vertex2id); - htable_triangle_init(&stardis->allocator, &triangle2id); - htable_descriptions_init(&stardis->allocator, &descriptions); - tables_initialised = 1; - /* loop on media */ - while (read_line(&line, input)){ - char* stl_filename = NULL; - struct description desc; - - init_description(&stardis->allocator, &desc); - res = parse_medium_line(line, &stl_filename, &desc); - if (res != RES_OK) goto error; - - /* Deduplicate media: find if the same description already exists - * (including name) */ - p_desc = htable_descriptions_find(&descriptions, &desc); - if (!p_desc) { - sa_push(stardis->descriptions, desc); - desc_id = sz++; - ASSERT(sz == sa_size(stardis->descriptions)); - /* Register new medium description */ - res = htable_descriptions_set(&descriptions, &desc, &desc_id); - if (res != RES_OK) goto error; - } else { - fprintf(stderr, "Duplicate media description found: deduplicated.\n"); - desc_id = *p_desc; + struct darray_descriptions* descriptions = context; + unsigned u_count, s_count, f_count, i_count; + const struct description* descs, * d = NULL; + + (void)itri; + ASSERT(descriptions && properties_conflict_status); + descs = darray_descriptions_cdata_get(descriptions); + u_count = s_count = f_count = i_count = 0; + + COUNT_SIDE(SG3D_FRONT); + COUNT_SIDE(SG3D_BACK); + if(properties[SG3D_INTFACE] == SG3D_UNSPECIFIED_PROPERTY) + u_count++; + else i_count++; + + ASSERT(u_count + s_count + f_count + i_count == 3); + if(i_count) { + ASSERT(properties[SG3D_INTFACE] < darray_descriptions_size_get(descriptions)); + d = descs + properties[SG3D_INTFACE]; + } + if(f_count > 1) { + /* Cannot have 2 fluids */ + *properties_conflict_status = 1; + } + else if(u_count > 1) { + /* Cannot have undefined interface if 1 side is undefined */ + *properties_conflict_status = 1; + } + else if(u_count == 1) { + if(!(s_count == 2 && !d) + && !(s_count == 1 && d + && (d->type == DESC_BOUND_H_FOR_SOLID || d->type == DESC_BOUND_H_FOR_SOLID))) + { + /* The only valid configurations are Solid / Undef / Solid + * and Solid / Solid_bound / Undef */ + *properties_conflict_status = 1; } - - res = read_stl(sz-1, CONTENT_GEOMETRY, stl_filename, &vertex2id, - &triangle2id, stardis); - if (stl_filename) free(stl_filename); - if (res != RES_OK) goto error; - } - fclose(input); - input = NULL; - - /* parse boundary files */ - input = fopen(bc_filename,"r"); - if (!input){ - fprintf(stderr,"Cannot open %s\n", bc_filename); - res = RES_IO_ERR; - goto error; - } - - /* loop on boundaries */ - while (read_line(&line, input)){ - char* stl_filename = NULL; - struct description desc; - - init_description(&stardis->allocator, &desc); - res = parse_boundary_line(line, &stl_filename, &desc); - if (res != RES_OK) goto error; - - /* Deduplicate media: find if the same description already exists - * (including name) */ - p_desc = htable_descriptions_find(&descriptions, &desc); - if (!p_desc) { - sa_push(stardis->descriptions, desc); - desc_id = sz++; - ASSERT(sz == sa_size(stardis->descriptions)); - /* Register new boundary description */ - res = htable_descriptions_set(&descriptions, &desc, &desc_id); - if (res != RES_OK) goto error; - } else { - fprintf(stderr, "Duplicate boundary description found: deduplicated.\n"); - desc_id = *p_desc; + } else { + ASSERT(u_count == 0 && d); + if(s_count != 1 || f_count != 1 || d->type != DESC_SOLID_FLUID_CONNECT) { + /* The only valid configuration is Solid / Solid_Fluid_bound / Fluid */ + *properties_conflict_status = 1; } - - res = read_stl(desc_id, CONTENT_CONNECTION, stl_filename, &vertex2id, - &triangle2id, stardis); - if (stl_filename) free(stl_filename); - if (res != RES_OK) goto error; } - fclose(input); - input = NULL; - /* Process STL file describing boundary compute */ - if (bound_filename) { - res = read_stl(0, CONTENT_BOUNDARY, bound_filename, &vertex2id, - &triangle2id, stardis); - if (res != RES_OK) goto error; - } - -exit: - if(input) fclose(input); - if (tables_initialised) { - htable_descriptions_release(&descriptions); - htable_vertex_release(&vertex2id); - htable_triangle_release(&triangle2id); - } - sa_release(line); - return res; -error: - goto exit; + /* No conflict found */ + *properties_conflict_status = 0; + return RES_OK; } +#undef COUNT_SIDE + /******************************************************************************* * Public Functions ******************************************************************************/ @@ -437,80 +136,118 @@ stardis_init struct stardis* stardis) { res_T res = RES_OK; - - ASSERT(args && stardis); - - res = mem_init_proxy_allocator(&stardis->allocator, - &mem_default_allocator); - if (res != RES_OK) goto error; - stardis->allocator_initialized = 1; - - res = geometry_analyse(args->medium_filename, - args->bc_filename, args->solve_boundary_filename, stardis); - if (res != RES_OK) goto error; - - stardis->N = args->N; + int geom_initialized = 0; + unsigned count; + int run = 1; + + ASSERT(args && stardis && args->allocator); + + stardis->allocator = args->allocator; + darray_descriptions_init(stardis->allocator, &stardis->descriptions); + stardis->probe[0] = args->probe[0]; + stardis->probe[1] = args->probe[1]; + stardis->probe[2] = args->probe[2]; + stardis->probe[3] = args->probe[3]; + init_camera(&stardis->camera); + str_init(stardis->allocator, &stardis->solve_name); + darray_size_t_init(stardis->allocator, &stardis->compute_region.primitives); + darray_sides_init(stardis->allocator, &stardis->compute_region.sides); + darray_uint_init(stardis->allocator, &stardis->compute_region.err_triangles); + stardis->samples = args->samples; stardis->nthreads = args->nthreads; - stardis->probe[0] = args->u.probe[0]; - stardis->probe[1] = args->u.probe[1]; - stardis->probe[2] = args->u.probe[2]; - stardis->probe[3] = args->u.probe[3]; stardis->scale_factor = args->scale_factor; stardis->radiative_temp[0] = args->radiative_temp[0]; stardis->radiative_temp[1] = args->radiative_temp[1]; stardis->dump_paths = SDIS_HEAT_PATH_NONE; - if (args->dump_paths & DUMP_ERROR) - stardis->dump_paths |= SDIS_HEAT_PATH_FAILED; - if (args->dump_paths & DUMP_SUCCESS) - stardis->dump_paths |= SDIS_HEAT_PATH_SUCCEED; - - if (args->mode & IR_COMPUTE){ - res = parse_camera(args->camera, &stardis->camera); - if (res != RES_OK) goto error; - } - else if (args->mode & MEDIUM_COMPUTE) { - strcpy(stardis->solve_name, args->medium_name); - } - else if (args->mode & (BOUNDARY_COMPUTE | MAP_COMPUTE)) { - strcpy(stardis->solve_name, args->solve_boundary_filename); - } + if(args->dump_paths & DUMP_ERROR) + stardis->dump_paths |= SDIS_HEAT_PATH_FAILURE; + if(args->dump_paths & DUMP_SUCCESS) + stardis->dump_paths |= SDIS_HEAT_PATH_SUCCESS; + stardis->verbose = args->verbose; + ERR(init_geometry(stardis->allocator, stardis->verbose, &stardis->geometry)); + geom_initialized = 1; + + if(args->mode & IR_COMPUTE) + ERR(parse_camera(args->camera, &stardis->camera)); + else if(args->mode & MEDIUM_COMPUTE) + ERR(str_set(&stardis->solve_name, args->medium_name)); + else if(args->mode & (BOUNDARY_COMPUTE | MAP_COMPUTE)) + ERR(str_set(&stardis->solve_name, args->solve_boundary_filename)); + + ERR(read_model(&args->model_files, stardis)); + + /* If computation is on a compute region, read it */ + if(args->mode & REGION_COMPUTE_MODES) { + ASSERT(!str_is_empty(&stardis->solve_name)); + ERR(read_compute_region(stardis)); + } + + /* Check conflicts */ + ERR(sg3d_geometry_get_unique_triangles_with_merge_conflict_count( + stardis->geometry.sg3d, &count)); + if(count) { + fprintf(stderr, + "Merge conflicts found in the the model (%u triangles).\n", count); + run = 0; + res = RES_BAD_ARG; + /* Don't goto error to allow further checking */ + } + ERR(sg3d_geometry_validate_properties(stardis->geometry.sg3d, + validate_properties, &stardis->descriptions)); + ERR(sg3d_geometry_get_unique_triangles_with_properties_conflict_count( + stardis->geometry.sg3d, &count)); + if(count) { + fprintf(stderr, + "Properties conflicts found in the the model (%u triangles).\n", count); + run = 0; + res = RES_BAD_ARG; + /* Don't goto error to allow further checking */ + } + if(args->mode & REGION_COMPUTE_MODES) { + ERR(sg3d_geometry_validate_properties(stardis->geometry.sg3d, + validate_properties, &stardis->descriptions)); + ERR(sg3d_geometry_get_unique_triangles_with_properties_conflict_count( + stardis->geometry.sg3d, &count)); + if(count) { + fprintf(stderr, "Invalid compute region defined by file '%s'.\n", + str_cget(&stardis->solve_name)); + fprintf(stderr, + "The file contains %u triangles not in the model.\n", count); + run = 0; + res = RES_BAD_ARG; + } + } + if(!run) { + fprintf(stderr, "Cannot run the solver.\n"); + goto error; + } exit: return res; error: + if(geom_initialized) release_geometry(&stardis->geometry); + darray_descriptions_release(&stardis->descriptions); + str_release(&stardis->solve_name); + darray_size_t_release(&stardis->compute_region.primitives); + darray_sides_release(&stardis->compute_region.sides); + darray_uint_release(&stardis->compute_region.err_triangles); goto exit; } void stardis_release(struct stardis* stardis) { - size_t memsz, i = 0; + struct sdis_interface** itf; + struct description* dsc; if(!stardis) return; - if (stardis->geometry.interfaces) { - for (i = 0; i < sa_size(stardis->geometry.interfaces); ++i) - SDIS(interface_ref_put(stardis->geometry.interfaces[i])); - sa_release(stardis->geometry.interfaces); - stardis->geometry.interfaces = NULL; - } - sa_release(stardis->geometry.interf_bytrg); - stardis->geometry.interf_bytrg = NULL; - sa_release(stardis->descriptions); + itf = darray_interface_ptrs_data_get(&stardis->geometry.interfaces); + dsc = darray_descriptions_data_get(&stardis->descriptions); + str_release(&stardis->solve_name); + darray_descriptions_release(&stardis->descriptions); release_geometry(&stardis->geometry); - sa_release(stardis->boundary.primitives); - sa_release(stardis->boundary.sides); - - if (stardis->allocator_initialized - && (memsz = MEM_ALLOCATED_SIZE(&stardis->allocator)) != 0) - { - char dump[4096] = { '\0' }; - MEM_DUMP(&stardis->allocator, dump, sizeof(dump)); - fprintf(stderr, "%s\n", dump); - fprintf(stderr, "Memory leaks: %zu Bytes\n", memsz); - } - if (stardis->allocator_initialized) { - mem_shutdown_proxy_allocator(&stardis->allocator); - stardis->allocator_initialized = 0; - } + darray_size_t_release(&stardis->compute_region.primitives); + darray_sides_release(&stardis->compute_region.sides); + darray_uint_release(&stardis->compute_region.err_triangles); } diff --git a/src/stardis-app.h b/src/stardis-app.h @@ -6,21 +6,21 @@ #include "stardis-parsing.h" #include <star/sstl.h> +#include <star/sg3d.h> + #include <rsys/rsys.h> #include <rsys/float3.h> -#include <rsys/stretchy_array.h> -#include <rsys/hash_table.h> +#include <rsys/double3.h> +#include <rsys/dynamic_array_size_t.h> +#include <rsys/dynamic_array.h> +#include <rsys/str.h> + #include <sdis.h> #include <limits.h> -#define FREE_AARRAY(ARRAY) if (ARRAY) {\ - int i = 0; \ - for (i=0; *(ARRAY+i);i++){\ - free(ARRAY[i]);\ - }\ - free(ARRAY);\ -} +/* Utility macros */ +#define ERR(Expr) if((res = (Expr)) != RES_OK) goto error; else (void)0 /* Different types of descriptions */ enum description_type { @@ -36,118 +36,101 @@ enum description_type { DESC_OUTSIDE }; -struct vertex{ - float xyz[3]; -}; -#define NULL_VERTEX__ {{0.0 ,0.0 ,0.0}} -static const struct vertex NULL_VERTEX = NULL_VERTEX__; +#define DARRAY_NAME interface_ptrs +#define DARRAY_DATA struct sdis_interface* +#include <rsys/dynamic_array.h> -static INLINE char -eq_vertex(const struct vertex* a, const struct vertex* b) -{ - return (char)f3_eq(a->xyz, b->xyz); -} - -/* Declare the hash table that map a vertex to its index */ -#define HTABLE_NAME vertex -#define HTABLE_DATA unsigned -#define HTABLE_KEY struct vertex -#define HTABLE_KEY_FUNCTOR_EQ eq_vertex -#include <rsys/hash_table.h> - -struct unsigned3 { unsigned data[3]; }; -struct triangle { - struct unsigned3 indices; - unsigned front_description; - unsigned back_description; - unsigned connection_description; +struct geometry { + struct sg3d_geometry* sg3d; + struct sdis_scene* sdis_scn; + struct darray_interface_ptrs interf_bytrg; + struct darray_interface_ptrs interfaces; }; -#define NULL_TRIANGLE__ {{{0, 0, 0}}, UINT_MAX, UINT_MAX, UINT_MAX } -static const struct triangle NULL_TRIANGLE = NULL_TRIANGLE__; - -static INLINE char -eq_indices(const struct unsigned3* a, const struct unsigned3* b) -{ - return a->data[0] == b->data[0] - && a->data[1] == b->data[1] - && a->data[2] == b->data[2]; -} static INLINE res_T -cp_indices(struct unsigned3* dst, const struct unsigned3* src) -{ - int i; - FOR_EACH(i, 0, 3) dst->data[i] = src->data[i]; - return RES_OK; +init_geometry + (struct mem_allocator* allocator, + const int verbose, + struct geometry* geom) +{ + res_T res = RES_OK; + struct sg3d_device* sg3d_dev = NULL; + struct sg3d_geometry* sg3d = NULL; + ASSERT(allocator && geom); + + ERR(sg3d_device_create(NULL, allocator, verbose, &sg3d_dev)); + ERR(sg3d_geometry_create(sg3d_dev, &sg3d)); + darray_interface_ptrs_init(allocator, &geom->interfaces); + darray_interface_ptrs_init(allocator, &geom->interf_bytrg); + +exit: + geom->sg3d = sg3d; + geom->sdis_scn = NULL; + if(sg3d_dev) SG3D(device_ref_put(sg3d_dev)); + return res; +error: + if(sg3d) SG3D(geometry_ref_put(sg3d)); + sg3d = NULL; + goto exit; } -/* Declare the hash table that maps a triangle id to its index */ -#define HTABLE_NAME triangle -#define HTABLE_DATA unsigned -#define HTABLE_KEY struct unsigned3 -#define HTABLE_KEY_FUNCTOR_EQ eq_indices -#define HTABLE_KEY_FUNCTOR_COPY cp_indices -#include <rsys/hash_table.h> - -struct geometry { - struct vertex* vertex; - struct triangle* triangle; - struct sdis_interface** interf_bytrg; - struct sdis_interface** interfaces; -}; -#define NULL_GEOMETRY__ {NULL, NULL, NULL, NULL } -static const struct geometry NULL_GEOMETRY = NULL_GEOMETRY__; - static INLINE void release_geometry(struct geometry* geom) { size_t i; - for (i = 0; i < sa_size(geom->interfaces); ++i) - SDIS(interface_ref_put(geom->interfaces[i])); - sa_release(geom->vertex); geom->vertex = NULL; - sa_release(geom->triangle); geom->triangle = NULL; - sa_release(geom->interfaces); geom->interfaces = NULL; - sa_release(geom->interf_bytrg); geom->interf_bytrg = NULL; + struct sdis_interface + **intf = darray_interface_ptrs_data_get(&geom->interfaces); + if(geom->sg3d) SG3D(geometry_ref_put(geom->sg3d)); + if(geom->sdis_scn) SDIS(scene_ref_put(geom->sdis_scn)); + for(i = 0; i < darray_interface_ptrs_size_get(&geom->interfaces); ++i) + SDIS(interface_ref_put(intf[i])); + darray_interface_ptrs_release(&geom->interfaces); + darray_interface_ptrs_release(&geom->interf_bytrg); } struct mat_fluid { - char name[32]; + struct str name; unsigned fluid_id; double rho; double cp; double tinit; double imposed_temperature; }; -#define NULL_FLUID__ { "", UINT_MAX, 0, 0, -1, -1 } -static const struct mat_fluid NULL_FLUID = NULL_FLUID__; + +static FINLINE res_T +init_fluid(struct mem_allocator* allocator, struct mat_fluid* dst) +{ + str_init(allocator, &dst->name); + dst->fluid_id = UINT_MAX; + dst->rho = 0; + dst->cp = 0; + dst->tinit = 0; + dst->imposed_temperature = -1; + return RES_OK; +} + +static FINLINE void +release_fluid(struct mat_fluid* fluid) +{ + str_release(&fluid->name); +} static void print_fluid(FILE* stream, const struct mat_fluid* f) { ASSERT(stream && f); - fprintf(stream, "Fluid '%s': cp=%g rho=%g", - f->name, f->cp, f->rho); - if (f->tinit >= 0) fprintf(stream, " Tinit=%f", f->tinit); - if (f->imposed_temperature >= 0) fprintf(stream, " Temp=%f", f->imposed_temperature); + fprintf(stream, "Fluid '%s': rho=%g cp=%g", + str_cget(&f->name), f->rho, f->cp); + if(f->tinit >= 0) fprintf(stream, " Tinit=%g", f->tinit); + if(f->imposed_temperature >= 0) + fprintf(stream, " Temp=%g", f->imposed_temperature); fprintf(stream, "\n"); -}static char - -eq_fluid(const struct mat_fluid* a, const struct mat_fluid* b) -{ - if (strcmp(a->name, b->name) - || a->fluid_id != b->fluid_id - || a->rho != b->rho - || a->cp != b->cp - || a->tinit != b->tinit - || a->imposed_temperature != b->imposed_temperature) - return 0; - return 1; } static FINLINE res_T cp_fluid(struct mat_fluid* dst, const struct mat_fluid* src) { - strcpy(dst->name, src->name); + str_copy(&dst->name, &src->name); dst->fluid_id = src->fluid_id; dst->rho = src->rho; dst->cp = src->cp; @@ -156,49 +139,8 @@ cp_fluid(struct mat_fluid* dst, const struct mat_fluid* src) return RES_OK; } -static size_t -hash_fluid(const struct mat_fluid* key) -{ - /* 64-bits Fowler/Noll/Vo hash function */ - const uint64_t FNV64_PRIME = - (uint64_t)(((uint64_t)1 << 40) + ((uint64_t)1 << 8) + 0xB3); - const uint64_t OFFSET64_BASIS = (uint64_t)14695981039346656037u; - uint64_t hash = OFFSET64_BASIS; - size_t i; - ASSERT(key); - FOR_EACH(i, 0, sizeof(key->name)) { - hash = hash ^ (uint64_t)((unsigned char)((const char*)&key->name)[i]); - hash = hash * FNV64_PRIME; - } - FOR_EACH(i, 0, sizeof(key->fluid_id)) { - hash = hash ^ (uint64_t)((unsigned char)((const char*)&key->fluid_id)[i]); - hash = hash * FNV64_PRIME; - } - FOR_EACH(i, 0, sizeof(key->rho)) { - hash = hash ^ (uint64_t)((unsigned char)((const char*)&key->rho)[i]); - hash = hash * FNV64_PRIME; - } - FOR_EACH(i, 0, sizeof(key->cp)) { - hash = hash ^ (uint64_t)((unsigned char)((const char*)&key->cp)[i]); - hash = hash * FNV64_PRIME; - } - hash = hash ^ (uint64_t)((unsigned char)(key->tinit == 0)); - hash = hash * FNV64_PRIME; - FOR_EACH(i, 0, sizeof(key->tinit)) { - hash = hash ^ (uint64_t)((unsigned char)((const char*)&key->tinit)[i]); - hash = hash * FNV64_PRIME; - } - hash = hash ^ (uint64_t)((unsigned char)(key->imposed_temperature == 0)); - hash = hash * FNV64_PRIME; - FOR_EACH(i, 0, sizeof(key->imposed_temperature)) { - hash = hash ^ (uint64_t)((unsigned char)((const char*)&key->imposed_temperature)[i]); - hash = hash * FNV64_PRIME; - } - return hash; -} - struct mat_solid { - char name[32]; + struct str name; unsigned solid_id; double lambda; double rho; @@ -208,41 +150,46 @@ struct mat_solid { double imposed_temperature; double vpower; }; -#define NULL_SOLID__ { "", UINT_MAX, 0, 0, 0, 0, -1, -1, 0 } -static const struct mat_solid NULL_SOLID = NULL_SOLID__; + +static FINLINE res_T +init_solid(struct mem_allocator* allocator, struct mat_solid* dst) +{ + str_init(allocator, &dst->name); + dst->solid_id = UINT_MAX; + dst->lambda = 0; + dst->rho = 0; + dst->cp = 0; + dst->delta = 0; + dst->tinit = -1; + dst->imposed_temperature = -1; + dst->vpower = 0; + return RES_OK; +} + +static FINLINE void +release_solid(struct mat_solid* solid) +{ + str_release(&solid->name); +} static void print_solid(FILE* stream, const struct mat_solid* s) { ASSERT(stream && s); fprintf(stream, - "Solid '%s': lambda=%g cp=%g rho=%g delta=%g Power=%f", - s->name, s->lambda, s->cp, s->rho, s->delta, s->vpower); - if (s->tinit >= 0) fprintf(stream, " Tinit=%f", s->tinit); - if (s->imposed_temperature >= 0) - fprintf(stream, " Temp=%f", s->imposed_temperature); + "Solid '%s': lambda=%g rho=%g cp=%g delta=%g", + str_cget(&s->name), s->lambda, s->rho, s->cp, s->delta); + if(s->vpower != 0) fprintf(stream, " VPower=%g", s->vpower); + if(s->tinit >= 0) fprintf(stream, " Tinit=%g", s->tinit); + if(s->imposed_temperature >= 0) + fprintf(stream, " Temp=%g", s->imposed_temperature); fprintf(stream, "\n"); } -static char -eq_solid(const struct mat_solid* a, const struct mat_solid* b) -{ - if (strcmp(a->name, b->name) - || a->solid_id != b->solid_id - || a->lambda != b->lambda - || a->rho != b->rho - || a->cp != b->cp - || a->delta != b->delta - || a->tinit != b->tinit - || a->imposed_temperature != b->imposed_temperature) - return 0; - return 1; -} - static FINLINE res_T cp_solid(struct mat_solid* dst, const struct mat_solid* src) { - strcpy(dst->name, src->name); + str_copy(&dst->name, &src->name); dst->solid_id = src->solid_id; dst->lambda = src->lambda; dst->rho = src->rho; @@ -250,75 +197,36 @@ cp_solid(struct mat_solid* dst, const struct mat_solid* src) dst->delta = src->delta; dst->tinit = src->tinit; dst->imposed_temperature = src->imposed_temperature; + dst->vpower = src->vpower; return RES_OK; } -static size_t -hash_solid(const struct mat_solid* key) -{ - /* 64-bits Fowler/Noll/Vo hash function */ - const uint64_t FNV64_PRIME = - (uint64_t)(((uint64_t)1 << 40) + ((uint64_t)1 << 8) + 0xB3); - const uint64_t OFFSET64_BASIS = (uint64_t)14695981039346656037u; - uint64_t hash = OFFSET64_BASIS; - size_t i; - ASSERT(key); - FOR_EACH(i, 0, sizeof(key->name)) { - hash = hash ^ (uint64_t)((unsigned char)((const char*)&key->name)[i]); - hash = hash * FNV64_PRIME; - } - FOR_EACH(i, 0, sizeof(key->solid_id)) { - hash = hash ^ (uint64_t)((unsigned char)((const char*)&key->solid_id)[i]); - hash = hash * FNV64_PRIME; - } - FOR_EACH(i, 0, sizeof(key->lambda)) { - hash = hash ^ (uint64_t)((unsigned char)((const char*)&key->lambda)[i]); - hash = hash * FNV64_PRIME; - } - FOR_EACH(i, 0, sizeof(key->rho)) { - hash = hash ^ (uint64_t)((unsigned char)((const char*)&key->rho)[i]); - hash = hash * FNV64_PRIME; - } - FOR_EACH(i, 0, sizeof(key->cp)) { - hash = hash ^ (uint64_t)((unsigned char)((const char*)&key->cp)[i]); - hash = hash * FNV64_PRIME; - } - FOR_EACH(i, 0, sizeof(key->delta)) { - hash = hash ^ (uint64_t)((unsigned char)((const char*)&key->delta)[i]); - hash = hash * FNV64_PRIME; - } - hash = hash ^ (uint64_t)((unsigned char)(key->tinit == 0)); - hash = hash * FNV64_PRIME; - FOR_EACH(i, 0,sizeof(key->tinit)) { - hash = hash ^ (uint64_t)((unsigned char)((const char*)&key->tinit)[i]); - hash = hash * FNV64_PRIME; - } - hash = hash ^ (uint64_t)((unsigned char)(key->imposed_temperature == 0)); - hash = hash * FNV64_PRIME; - FOR_EACH(i, 0, sizeof(key->imposed_temperature)) { - hash = hash ^ (uint64_t)((unsigned char) - ((const char*)&key->imposed_temperature)[i]); - hash = hash * FNV64_PRIME; - } - FOR_EACH(i, 0, sizeof(key->vpower)) { - hash = hash ^ (uint64_t)((unsigned char)((const char*)&key->vpower)[i]); - hash = hash * FNV64_PRIME; - } - return hash; -} - struct h_boundary { - char name[32]; + struct str name; unsigned mat_id; double emissivity; double specular_fraction; double hc; - double hc_max; double imposed_temperature; - int has_emissivity, has_hc; }; -#define NULL_HBOUND__ { "", UINT_MAX, 0, 0, 0, -1, 0, 0 } -static const struct h_boundary NULL_HBOUND = NULL_HBOUND__; + +static FINLINE res_T +init_h(struct mem_allocator* allocator, struct h_boundary* dst) +{ + str_init(allocator, &dst->name); + dst->mat_id = UINT_MAX; + dst->emissivity = 0; + dst->specular_fraction = 0; + dst->hc = 0; + dst->imposed_temperature = -1; + return RES_OK; +} + +static FINLINE void +release_h_boundary(struct h_boundary* bound) +{ + str_release(&bound->name); +} static void print_h_boundary @@ -329,108 +237,47 @@ print_h_boundary ASSERT(stream && b && (type == DESC_BOUND_H_FOR_SOLID || type == DESC_BOUND_H_FOR_FLUID)); fprintf(stream, - "H boundary %s for %s: emissivity=%g specular_fraction=%g hc=%g hc_max=%g T=%f\n", - b->name, + "H boundary for %s '%s': emissivity=%g specular_fraction=%g hc=%g T=%g\n", (type == DESC_BOUND_H_FOR_SOLID ? "solid" : "fluid"), - (b->has_emissivity ? b->emissivity : 0), - (b->has_emissivity ? b->specular_fraction : 0), - (b->has_hc ? b->hc : 0), - (b->has_hc ? b->hc_max : 0), + str_cget(&b->name), + b->emissivity, b->specular_fraction, b->hc, b->imposed_temperature); } -static char -eq_h_boundary(const struct h_boundary* a, const struct h_boundary* b) -{ - if (a->mat_id != b->mat_id - || strcmp(a->name, b->name) - || a->specular_fraction != b->specular_fraction - || a->imposed_temperature != b->imposed_temperature - || a->has_emissivity != b->has_emissivity - || a->has_hc != b->has_hc) - return 0; - if (a->has_emissivity && a->emissivity != b->emissivity) return 0; - if (a->has_hc && (a->hc != b->hc || a->hc_max != b->hc_max)) return 0; - return 1; -} - static FINLINE res_T cp_h_boundary(struct h_boundary* dst, const struct h_boundary* src) { - strcpy(dst->name, src->name); + str_copy(&dst->name, &src->name); dst->mat_id = src->mat_id; dst->specular_fraction = src->specular_fraction; dst->imposed_temperature = src->imposed_temperature; - dst->has_emissivity = src->has_emissivity; - dst->has_hc = src->has_hc; dst->emissivity = src->emissivity; - dst->hc_max = src->hc_max; dst->hc = src->hc; return RES_OK; } -static size_t -hash_h_boundary(const struct h_boundary* key) -{ - /* 64-bits Fowler/Noll/Vo hash function */ - const uint64_t FNV64_PRIME = - (uint64_t)(((uint64_t)1 << 40) + ((uint64_t)1 << 8) + 0xB3); - const uint64_t OFFSET64_BASIS = (uint64_t)14695981039346656037u; - uint64_t hash = OFFSET64_BASIS; - size_t i; - ASSERT(key); - FOR_EACH(i, 0, sizeof(key->name)) { - hash = hash ^ (uint64_t)((unsigned char)((const char*)&key->name)[i]); - hash = hash * FNV64_PRIME; - } - FOR_EACH(i, 0, sizeof(key->mat_id)) { - hash = hash ^ (uint32_t)((unsigned char)((const char*)&key->mat_id)[i]); - hash = hash * FNV64_PRIME; - } - FOR_EACH(i, 0, sizeof(key->emissivity)) { - hash = hash ^ (uint32_t)((unsigned char)((const char*)&key->emissivity)[i]); - hash = hash * FNV64_PRIME; - } - FOR_EACH(i, 0, sizeof(key->specular_fraction)) { - hash = hash ^ (uint32_t)((unsigned char) - ((const char*)&key->specular_fraction)[i]); - hash = hash * FNV64_PRIME; - } - FOR_EACH(i, 0, sizeof(key->imposed_temperature)) { - hash = hash ^ (uint32_t)((unsigned char) - ((const char*)&key->imposed_temperature)[i]); - hash = hash * FNV64_PRIME; - } - FOR_EACH(i, 0, sizeof(key->has_emissivity)) { - hash = hash ^ (uint32_t)((unsigned char)((const char*)&key->has_emissivity)[i]); - hash = hash * FNV64_PRIME; - } - FOR_EACH(i, 0, sizeof(key->has_hc)) { - hash = hash ^ (uint32_t)((unsigned char)((const char*)&key->has_hc)[i]); - hash = hash * FNV64_PRIME; - } - if (!key->has_hc) return hash; - FOR_EACH(i, 0, sizeof(key->hc)) { - hash = hash ^ (uint32_t)((unsigned char)((const char*)&key->hc)[i]); - hash = hash * FNV64_PRIME; - } - FOR_EACH(i, 0, sizeof(key->hc_max)) { - hash = hash ^ (uint32_t)((unsigned char)((const char*)&key->hc_max)[i]); - hash = hash * FNV64_PRIME; - } - return hash; -} - struct t_boundary { - char name[32]; + struct str name; unsigned mat_id; - double imposed_temperature; double hc; - double hc_max; - int has_hc; + double imposed_temperature; }; -#define NULL_TBOUND__ { "", UINT_MAX, -1, 0, 0, 0 } -static const struct t_boundary NULL_TBOUND = NULL_TBOUND__; + +static FINLINE res_T +init_t(struct mem_allocator* allocator, struct t_boundary* dst) +{ + str_init(allocator, &dst->name); + dst->mat_id = UINT_MAX; + dst->hc = 0; + dst->imposed_temperature = -1; + return RES_OK; +} + +static FINLINE void +release_t_boundary(struct t_boundary* bound) +{ + str_release(&bound->name); +} static void print_t_boundary @@ -441,98 +288,42 @@ print_t_boundary ASSERT(stream && b && (type == DESC_BOUND_T_FOR_SOLID || type == DESC_BOUND_T_FOR_FLUID)); fprintf(stream, - "T boundary %s for %s: hc=%g hc_max=%g T=%f\n", - b->name, + "T boundary for %s' %s': hc=%g T=%f\n", (type == DESC_BOUND_T_FOR_SOLID ? "solid" : "fluid"), - (b->has_hc ? b->hc : 0), - (b->has_hc ? b->hc_max : 0), - b->imposed_temperature); -} - -static char -eq_t_boundary - (const struct t_boundary* a, - const struct t_boundary* b, - const enum description_type type) -{ - ASSERT(type == DESC_BOUND_T_FOR_SOLID || type == DESC_BOUND_T_FOR_FLUID); - /* These fields are meaningful by both types */ - if (strcmp(a->name, b->name)) - return 0; - if (a->mat_id != b->mat_id || a->imposed_temperature != b->imposed_temperature) - return 0; - if (type != DESC_BOUND_T_FOR_FLUID) - return 1; - /* These ones are only relevant for fluids */ - if (a->has_hc != b->has_hc) - return 0; - if (a->has_hc && (a->hc != b->hc || a->hc_max != b->hc_max)) - return 0; - return 1; + str_cget(&b->name), + b->hc, b->imposed_temperature); } static FINLINE res_T cp_t_boundary(struct t_boundary* dst, const struct t_boundary* src) { - strcpy(dst->name, src->name); + str_copy(&dst->name, &src->name); dst->mat_id = src->mat_id; dst->imposed_temperature = src->imposed_temperature; - dst->has_hc = src->has_hc; - dst->hc_max = src->hc_max; dst->hc = src->hc; return RES_OK; } -static size_t -hash_t_boundary -(const struct t_boundary* key, - const enum description_type type) -{ - (void)type; - ASSERT(type == DESC_BOUND_T_FOR_SOLID || type == DESC_BOUND_T_FOR_FLUID); - /* 64-bits Fowler/Noll/Vo hash function */ - const uint64_t FNV64_PRIME = - (uint64_t)(((uint64_t)1 << 40) + ((uint64_t)1 << 8) + 0xB3); - const uint64_t OFFSET64_BASIS = (uint64_t)14695981039346656037u; - uint64_t hash = OFFSET64_BASIS; - size_t i; - ASSERT(key); - FOR_EACH(i, 0, sizeof(key->name)) { - hash = hash ^ (uint64_t)((unsigned char)((const char*)&key->name)[i]); - hash = hash * FNV64_PRIME; - } - FOR_EACH(i, 0, sizeof(key->mat_id)) { - hash = hash ^ (uint32_t)((unsigned char)((const char*)&key->mat_id)[i]); - hash = hash * FNV64_PRIME; - } - FOR_EACH(i, 0, sizeof(key->has_hc)) { - hash = hash ^ (uint32_t)((unsigned char)((const char*)&key->has_hc)[i]); - hash = hash * FNV64_PRIME; - } - FOR_EACH(i, 0, sizeof(key->imposed_temperature)) { - hash = hash ^ (uint32_t)((unsigned char) - ((const char*)&key->imposed_temperature)[i]); - hash = hash * FNV64_PRIME; - } - if (!key->has_hc) return hash; - FOR_EACH(i, 0, sizeof(key->hc)) { - hash = hash ^ (uint32_t)((unsigned char)((const char*)&key->hc)[i]); - hash = hash * FNV64_PRIME; - } - FOR_EACH(i, 0, sizeof(key->hc_max)) { - hash = hash ^ (uint32_t)((unsigned char)((const char*)&key->hc_max)[i]); - hash = hash * FNV64_PRIME; - } - return hash; -} - struct f_boundary { - char name[32]; + struct str name; unsigned mat_id; double imposed_flux; }; -#define NULL_FBOUND__ { "", UINT_MAX, SDIS_FLUX_NONE } -static const struct f_boundary NULL_FBOUND = NULL_FBOUND__; + +static FINLINE res_T +init_f(struct mem_allocator* allocator, struct f_boundary* dst) +{ + str_init(allocator, &dst->name); + dst->mat_id = UINT_MAX; + dst->imposed_flux = -1; + return RES_OK; +} + +static FINLINE void +release_f_boundary(struct f_boundary* bound) +{ + str_release(&bound->name); +} static void print_f_boundary @@ -541,63 +332,42 @@ print_f_boundary { ASSERT(stream && b); fprintf(stream, - "F boundary %s for SOLID: flux=%f\n", - b->name, + "F boundary for SOLID '%s': flux=%f\n", + str_cget(&b->name), b->imposed_flux); } -static char -eq_f_boundary(const struct f_boundary* a, const struct f_boundary* b) -{ - return (a->mat_id == b->mat_id - && 0 == strcmp(a->name, b->name) - && a->imposed_flux == b->imposed_flux); -} - static FINLINE res_T cp_f_boundary(struct f_boundary* dst, const struct f_boundary* src) { - strcpy(dst->name, src->name); + str_copy(&dst->name, &src->name); dst->mat_id = src->mat_id; dst->imposed_flux = src->imposed_flux; return RES_OK; } -static size_t -hash_f_boundary(const struct f_boundary* key) -{ - /* 64-bits Fowler/Noll/Vo hash function */ - const uint64_t FNV64_PRIME = - (uint64_t)(((uint64_t)1 << 40) + ((uint64_t)1 << 8) + 0xB3); - const uint64_t OFFSET64_BASIS = (uint64_t)14695981039346656037u; - uint64_t hash = OFFSET64_BASIS; - size_t i; - ASSERT(key); - FOR_EACH(i, 0, sizeof(key->name)) { - hash = hash ^ (uint64_t)((unsigned char)((const char*)&key->name)[i]); - hash = hash * FNV64_PRIME; - } - FOR_EACH(i, 0, sizeof(key->mat_id)) { - hash = hash ^ (uint32_t)((unsigned char)((const char*)&key->mat_id)[i]); - hash = hash * FNV64_PRIME; - } - FOR_EACH(i, 0, sizeof(key->imposed_flux)) { - hash = hash ^ (uint32_t)((unsigned char)((const char*)&key->imposed_flux)[i]); - hash = hash * FNV64_PRIME; - } - return hash; -} - struct solid_fluid_connect { - char name[32]; + struct str name; double emissivity; double specular_fraction; double hc; - double hc_max; - int has_emissivity, has_hc; }; -#define NULL_SFCONNECT__ { "", 0, 0, 0, 0, 0, 0} -static const struct solid_fluid_connect NULL_SFCONNECT = NULL_SFCONNECT__; + +static FINLINE void +release_sf_connect(struct solid_fluid_connect* connect) +{ + str_release(&connect->name); +} + +static FINLINE res_T +init_sf(struct mem_allocator* allocator, struct solid_fluid_connect* dst) +{ + str_init(allocator, &dst->name); + dst->emissivity = 0; + dst->specular_fraction = 0; + dst->hc = 0; + return RES_OK; +} static void print_sf_connect @@ -605,35 +375,19 @@ print_sf_connect const struct solid_fluid_connect* c) { ASSERT(stream && c); - fprintf(stream, "Solid-Fluid connection %s:", c->name); - if(c->has_emissivity) - fprintf(stream, " emissivity=%f, specular_fraction=%f", + fprintf(stream, "Solid-Fluid connection '%s':", str_cget(&c->name)); + fprintf(stream, " emissivity=%f, specular_fraction=%f", c->emissivity, c->specular_fraction); - if (c->has_hc) fprintf(stream, " hc=%f", c->hc); + fprintf(stream, " hc=%f", c->hc); fprintf(stream, "\n"); } -static char -eq_sf_connect(const struct solid_fluid_connect* a, const struct solid_fluid_connect* b) -{ - return (char)(!strcmp(a->name, b->name) - && a->emissivity == b->emissivity - && a->specular_fraction == b->specular_fraction - && a->hc == b->hc - && a->hc_max == b->hc_max - && a->has_emissivity == b->has_emissivity - && a->has_hc == b->has_hc); -} - static FINLINE res_T cp_sf_connect(struct solid_fluid_connect* dst, const struct solid_fluid_connect* src) { - strcpy(dst->name, src->name); + str_copy(&dst->name, &src->name); dst->specular_fraction = src->specular_fraction; - dst->has_emissivity = src->has_emissivity; - dst->has_hc = src->has_hc; dst->emissivity = src->emissivity; - dst->hc_max = src->hc_max; dst->hc = src->hc; return RES_OK; } @@ -644,56 +398,9 @@ cp_release_sf_connect(struct solid_fluid_connect* dst, struct solid_fluid_connec return cp_sf_connect(dst, src); } -static FINLINE res_T -release_sf_connect(struct solid_fluid_connect* a) -{ - (void)a; - return RES_OK; -} - -static size_t -hash_sf_connect(const struct solid_fluid_connect* key) -{ - /* 64-bits Fowler/Noll/Vo hash function */ - const uint64_t FNV64_PRIME = - (uint64_t)(((uint64_t)1 << 40) + ((uint64_t)1 << 8) + 0xB3); - const uint64_t OFFSET64_BASIS = (uint64_t)14695981039346656037u; - uint64_t hash = OFFSET64_BASIS; - size_t i; - ASSERT(key); - FOR_EACH(i, 0, sizeof(key->name)) { - hash = hash ^ (uint64_t)((unsigned char)((const char*)&key->name)[i]); - hash = hash * FNV64_PRIME; - } - FOR_EACH(i, 0, sizeof(key->emissivity)) { - hash = hash ^ (uint32_t)((unsigned char)((const char*)&key->emissivity)[i]); - hash = hash * FNV64_PRIME; - } - FOR_EACH(i, 0, sizeof(key->specular_fraction)) { - hash = hash ^ (uint32_t)((unsigned char)((const char*)&key->specular_fraction)[i]); - hash = hash * FNV64_PRIME; - } - FOR_EACH(i, 0, sizeof(key->hc)) { - hash = hash ^ (uint32_t)((unsigned char)((const char*)&key->hc)[i]); - hash = hash * FNV64_PRIME; - } - FOR_EACH(i, 0, sizeof(key->hc_max)) { - hash = hash ^ (uint32_t)((unsigned char)((const char*)&key->hc_max)[i]); - hash = hash * FNV64_PRIME; - } - FOR_EACH(i, 0, sizeof(key->has_emissivity)) { - hash = hash ^ (uint32_t)((unsigned char)((const char*)&key->has_emissivity)[i]); - hash = hash * FNV64_PRIME; - } - FOR_EACH(i, 0, sizeof(key->has_hc)) { - hash = hash ^ (uint32_t)((unsigned char)((const char*)&key->has_hc)[i]); - hash = hash * FNV64_PRIME; - } - return hash; -} struct description { enum description_type type; - union d { + union { struct mat_fluid fluid; struct mat_solid solid; struct t_boundary t_boundary; @@ -703,116 +410,133 @@ struct description { } d; }; -static INLINE void -print_description - (FILE* stream, - const struct description* desc) +static FINLINE res_T +init_description(struct mem_allocator* alloc, struct description* desc) +{ + ASSERT(desc); + (void)alloc; + desc->type = DESCRIPTION_TYPE_COUNT__; + return RES_OK; +} + +static FINLINE void +release_description(struct description* desc) { - ASSERT(stream && desc); switch (desc->type) { case DESC_MAT_SOLID: - print_solid(stream, &desc->d.solid); + release_solid(&desc->d.solid); break; case DESC_MAT_FLUID: - print_fluid(stream, &desc->d.fluid); - break; - case DESC_BOUND_T_FOR_SOLID: - case DESC_BOUND_T_FOR_FLUID: - print_t_boundary(stream, &desc->d.t_boundary, desc->type); + release_fluid(&desc->d.fluid); break; case DESC_BOUND_H_FOR_SOLID: case DESC_BOUND_H_FOR_FLUID: - print_h_boundary(stream, &desc->d.h_boundary, desc->type); + release_h_boundary(&desc->d.h_boundary); + break; + case DESC_BOUND_T_FOR_SOLID: + case DESC_BOUND_T_FOR_FLUID: + release_t_boundary(&desc->d.t_boundary); break; case DESC_BOUND_F_FOR_SOLID: - print_f_boundary(stream, &desc->d.f_boundary); + release_f_boundary(&desc->d.f_boundary); break; case DESC_SOLID_FLUID_CONNECT: - print_sf_connect(stream, &desc->d.sf_connect); + release_sf_connect(&desc->d.sf_connect); break; default: FATAL("error:" STR(__FILE__) ":" STR(__LINE__)": Invalid type.\n"); } } -static INLINE char -eq_description(const struct description* a, const struct description* b) +static INLINE void +print_description + (FILE* stream, + const struct description* desc) { - if (a->type != b->type) return 0; - switch (a->type) { + ASSERT(stream && desc); + switch (desc->type) { case DESC_MAT_SOLID: - return eq_solid(&a->d.solid, &b->d.solid); + print_solid(stream, &desc->d.solid); + break; case DESC_MAT_FLUID: - return eq_fluid(&a->d.fluid, &b->d.fluid); - case DESC_BOUND_H_FOR_SOLID: - case DESC_BOUND_H_FOR_FLUID: - return eq_h_boundary(&a->d.h_boundary, &b->d.h_boundary); + print_fluid(stream, &desc->d.fluid); + break; case DESC_BOUND_T_FOR_SOLID: case DESC_BOUND_T_FOR_FLUID: - return eq_t_boundary(&a->d.t_boundary, &b->d.t_boundary, a->type); - case DESC_BOUND_F_FOR_SOLID: - return eq_f_boundary(&a->d.f_boundary, &b->d.f_boundary); - case DESC_SOLID_FLUID_CONNECT: - return eq_sf_connect(&a->d.sf_connect, &b->d.sf_connect); - default: - FATAL("error:" STR(__FILE__) ":" STR(__LINE__)": Invalid type.\n"); - } -} - -static INLINE size_t -hash_description(const struct description* key) -{ - switch (key->type) { - case DESC_MAT_SOLID: - return hash_solid(&key->d.solid); - case DESC_MAT_FLUID: - return hash_fluid(&key->d.fluid); + print_t_boundary(stream, &desc->d.t_boundary, desc->type); + break; case DESC_BOUND_H_FOR_SOLID: case DESC_BOUND_H_FOR_FLUID: - return hash_h_boundary(&key->d.h_boundary); - case DESC_BOUND_T_FOR_SOLID: - case DESC_BOUND_T_FOR_FLUID: - return hash_t_boundary(&key->d.t_boundary, key->type); + print_h_boundary(stream, &desc->d.h_boundary, desc->type); + break; case DESC_BOUND_F_FOR_SOLID: - return hash_f_boundary(&key->d.f_boundary); + print_f_boundary(stream, &desc->d.f_boundary); + break; case DESC_SOLID_FLUID_CONNECT: - return hash_sf_connect(&key->d.sf_connect); + print_sf_connect(stream, &desc->d.sf_connect); + break; default: FATAL("error:" STR(__FILE__) ":" STR(__LINE__)": Invalid type.\n"); } } static FINLINE res_T -init_description(struct mem_allocator* alloc, struct description* key) -{ - (void)alloc; - key->type = DESCRIPTION_TYPE_COUNT__; - return RES_OK; -} - -static FINLINE res_T cp_description(struct description* dst, const struct description* src) { + res_T res = RES_OK; ASSERT(src && dst); - dst->type = src->type; - switch (dst->type) { + switch (src->type) { case DESC_MAT_SOLID: - return cp_solid(&dst->d.solid, &src->d.solid); + if(dst->type == DESCRIPTION_TYPE_COUNT__) { + dst->type = src->type; + ERR(init_solid(src->d.solid.name.allocator, &dst->d.solid)); + } + ERR(cp_solid(&dst->d.solid, &src->d.solid)); + break; case DESC_MAT_FLUID: - return cp_fluid(&dst->d.fluid, &src->d.fluid); + if(dst->type == DESCRIPTION_TYPE_COUNT__) { + dst->type = src->type; + ERR(init_fluid(src->d.fluid.name.allocator, &dst->d.fluid)); + } + ERR(cp_fluid(&dst->d.fluid, &src->d.fluid)); + break; case DESC_BOUND_H_FOR_SOLID: case DESC_BOUND_H_FOR_FLUID: - return cp_h_boundary(&dst->d.h_boundary, &src->d.h_boundary); + if(dst->type == DESCRIPTION_TYPE_COUNT__) { + dst->type = src->type; + ERR(init_h(src->d.h_boundary.name.allocator, &dst->d.h_boundary)); + } + ERR(cp_h_boundary(&dst->d.h_boundary, &src->d.h_boundary)); + break; case DESC_BOUND_T_FOR_SOLID: case DESC_BOUND_T_FOR_FLUID: - return cp_t_boundary(&dst->d.t_boundary, &src->d.t_boundary); + if(dst->type == DESCRIPTION_TYPE_COUNT__) { + dst->type = src->type; + ERR(init_t(src->d.t_boundary.name.allocator, &dst->d.t_boundary)); + } + ERR(cp_t_boundary(&dst->d.t_boundary, &src->d.t_boundary)); + break; case DESC_BOUND_F_FOR_SOLID: - return cp_f_boundary(&dst->d.f_boundary, &src->d.f_boundary); + if(dst->type == DESCRIPTION_TYPE_COUNT__) { + dst->type = src->type; + ERR(init_f(src->d.f_boundary.name.allocator, &dst->d.f_boundary)); + } + ERR(cp_f_boundary(&dst->d.f_boundary, &src->d.f_boundary)); + break; case DESC_SOLID_FLUID_CONNECT: - return cp_sf_connect(&dst->d.sf_connect, &src->d.sf_connect); + if(dst->type == DESCRIPTION_TYPE_COUNT__) { + dst->type = src->type; + ERR(init_sf(src->d.sf_connect.name.allocator, &dst->d.sf_connect)); + } + ERR(cp_sf_connect(&dst->d.sf_connect, &src->d.sf_connect)); + break; default: FATAL("error:" STR(__FILE__) ":" STR(__LINE__)": Invalid type.\n"); } +end: + return res; +error: + goto end; } struct camera { @@ -822,66 +546,69 @@ struct camera { double fov; unsigned spp; unsigned img[2]; - union { /* Trick to allow static initialization with INF */ - uint64_t t; - double time; - } u; + double time; }; -#define NULL_CAMERA__ {\ - {1,1,1},{0,0,0},{0,0,1},30,4,{640,480},\ - { 0x7FF0000000000000 } /* time=INF */\ -} -static const struct camera NULL_CAMERA = NULL_CAMERA__; - -struct solve_boundary { - size_t* primitives; - enum sdis_side* sides; + +static void +init_camera(struct camera* cam) { + d3_splat(cam->pos, 0); + d3_splat(cam->tgt, 0); + d3_splat(cam->up, 0); + cam->fov = 0; + cam->spp = 0; + cam->img[0] = cam->img[1] = 0; + cam->time = INF; +} + +/* Type to store the primitives of a compute region */ +#define DARRAY_NAME sides +#define DARRAY_DATA enum sdis_side +#include <rsys/dynamic_array.h> + +#define DARRAY_NAME descriptions +#define DARRAY_DATA struct description +#define DARRAY_FUNCTOR_INIT init_description +#define DARRAY_FUNCTOR_COPY cp_description +#define DARRAY_FUNCTOR_RELEASE release_description +#include <rsys/dynamic_array.h> + +struct compute_region { + struct darray_size_t primitives; + struct darray_sides sides; + struct darray_uint err_triangles; }; struct stardis { struct geometry geometry; - struct description* descriptions; /*array of materials and boundaries */ + struct darray_descriptions descriptions; /* Materials and boundaries */ double probe[4]; /* x,y,z,t of probe when mode is PROBE_COMPUTE */ struct camera camera; /* camera when mode is IR_COMPUTE */ - char solve_name[64]; /* medium name when mode is MEDIUM_COMPUTE, - boundary file name when BOUNDARY_COMPUTE*/ - struct solve_boundary boundary; /* Boundary to solve when mode - is BOUNDARY_COMPUTE */ - - size_t N; /*number of MC realizations*/ + struct str solve_name; /* medium name when mode is MEDIUM_COMPUTE, + boundary file name when [FLUX_]BOUNDARY_COMPUTE + map file name when MAP_COMPUTE */ + struct compute_region compute_region; /* Compute region when mode is + BOUNDARY_COMPUTE or BOUNDARY_COMPUTE */ + size_t samples; unsigned nthreads; double scale_factor; double radiative_temp[2]; - struct mem_allocator allocator; - int allocator_initialized; + struct mem_allocator* allocator; enum sdis_heat_path_flag dump_paths; -}; -#define NULL_ALLOCATOR__ {NULL} -#define NULL_STARDIS__ {\ - NULL_GEOMETRY__, NULL,\ - {0,0,0,0}, NULL_CAMERA__, "", {NULL,NULL},\ - 0, 0, 0, {300,300}, NULL_ALLOCATOR__, 0, SDIS_HEAT_PATH_NONE} -static const struct stardis NULL_STARDIS = NULL_STARDIS__; - -enum content_type { - CONTENT_GEOMETRY, - CONTENT_CONNECTION, - CONTENT_BOUNDARY, - CONTENT_TYPES_COUNT + int verbose; }; -extern res_T +extern LOCAL_SYM res_T stardis_init (const struct args* args, struct stardis* stardis); -extern res_T +extern LOCAL_SYM res_T stardis_compute (struct stardis* stardis, const int mode); -extern void +extern LOCAL_SYM void stardis_release (struct stardis* stardis); diff --git a/src/stardis-compute.c b/src/stardis-compute.c @@ -9,9 +9,11 @@ #include <sdis.h> #include <sdis_version.h> + +#include <star/sg3d_sdisXd_helper.h> + #include <rsys/double3.h> #include <rsys/double2.h> - #include <rsys/image.h> struct dummies { @@ -29,43 +31,6 @@ static const struct dummies DUMMIES_NULL = DUMMIES_NULL__; * Local Functions ******************************************************************************/ -static void -geometry_get_position - (const size_t ivert, - double pos[3], - void* context) -{ - struct geometry* geom = context; - ASSERT(geom); - pos[0] = geom->vertex[ivert].xyz[0]; - pos[1] = geom->vertex[ivert].xyz[1]; - pos[2] = geom->vertex[ivert].xyz[2]; -} - -static void -geometry_get_indices - (const size_t itri, - size_t ids[3], - void* context) -{ - struct geometry* geom = context; - ASSERT(geom); - ids[0] = geom->triangle[itri].indices.data[0]; - ids[1] = geom->triangle[itri].indices.data[1]; - ids[2] = geom->triangle[itri].indices.data[2]; -} - -static void -geometry_get_interface - (const size_t itri, - struct sdis_interface** interf, - void* context) -{ - struct geometry* geom = context; - ASSERT(geom); - *interf = geom->interf_bytrg[itri]; -} - /* Check probe position and move it to the closest point on an interface * if on_interface is set */ static res_T @@ -80,38 +45,37 @@ select_probe_type res_T res = RES_OK; double min_d = DBL_MAX, min_delta = DBL_MAX, new_pos[3], tmp_uv[2], min_uv[2]; const struct description *min_desc = NULL; - size_t i = 0, found = SIZE_MAX; + size_t found = SIZE_MAX; enum description_type min_type = DESCRIPTION_TYPE_COUNT__; - const struct triangle* triangle; - const struct vertex* vertex; const struct description* descriptions; + unsigned i, tcount; ASSERT(scn && stardis && pos && iprim && uv); - triangle = stardis->geometry.triangle; - vertex = stardis->geometry.vertex; - descriptions = stardis->descriptions; - for (i = 0; i < sa_size(triangle); ++i) { - const struct triangle* tr = triangle + i; - double dp[3], d, projected_pos[3], nn[3]; + descriptions = darray_descriptions_cdata_get(&stardis->descriptions); + ERR(sg3d_geometry_get_unique_triangles_count(stardis->geometry.sg3d, &tcount)); + for(i = 0; i < tcount; ++i) { + unsigned indices[3], properties[3]; + double dp[3], d, projected_pos[3], nn[3], tmp[3]; float v0[3], v1[3], v2[3], e1[3], e2[3], n[3]; const struct description *desc = NULL; enum sdis_side s; - res = sdis_scene_boundary_project_position(scn, i, pos, tmp_uv); - if (res != RES_OK) return res; - res = sdis_scene_get_boundary_position(scn, i, tmp_uv, projected_pos); - if (res != RES_OK) return res; + ERR(sg3d_geometry_get_unique_triangle_vertices(stardis->geometry.sg3d, i, + indices)); + + ERR(sdis_scene_boundary_project_position(scn, i, pos, tmp_uv)); + ERR(sdis_scene_get_boundary_position(scn, i, tmp_uv, projected_pos)); d3_sub(dp, projected_pos, pos); d = d3_len(dp); - if (d == 0) { + if(d == 0) { /* Best possible match */ found = i; fprintf(stderr, "The probe is on the primitive %llu.\n", (long long int)*iprim); break; } - if (d >= min_d) { + if(d >= min_d) { /* No improvement */ continue; } @@ -119,9 +83,12 @@ select_probe_type d2_set(min_uv, tmp_uv); /* Find side */ - f3_set(v0, vertex[tr->indices.data[0]].xyz); - f3_set(v1, vertex[tr->indices.data[1]].xyz); - f3_set(v2, vertex[tr->indices.data[2]].xyz); + ERR(sg3d_geometry_get_unique_vertex(stardis->geometry.sg3d, indices[0], tmp)); + f3_set_d3(v0, tmp); + ERR(sg3d_geometry_get_unique_vertex(stardis->geometry.sg3d, indices[1], tmp)); + f3_set_d3(v1, tmp); + ERR(sg3d_geometry_get_unique_vertex(stardis->geometry.sg3d, indices[2], tmp)); + f3_set_d3(v2, tmp); f3_sub(e1, v1, v0); f3_sub(e2, v2, v0); f3_cross(n, e1, e2); @@ -130,49 +97,47 @@ select_probe_type d3_set_f3(nn, n); s = (d3_dot(nn, dp) < 0) ? SDIS_BACK : SDIS_FRONT; - if (s == SDIS_FRONT) { - if (tr->front_description != UINT_MAX) { - desc = descriptions + tr->front_description; - } - } else { - if (tr->back_description != UINT_MAX) { - desc = descriptions + tr->back_description; - } - } + ERR(sg3d_geometry_get_unique_triangle_properties(stardis->geometry.sg3d, i, + properties)); + + if(s == SDIS_FRONT && properties[SG3D_FRONT] != SG3D_UNSPECIFIED_PROPERTY) + desc = descriptions + properties[SG3D_FRONT]; + if(s == SDIS_BACK && properties[SDIS_BACK] != SG3D_UNSPECIFIED_PROPERTY) + desc = descriptions + properties[SG3D_BACK]; min_desc = desc; - if (!desc) continue; + if(!desc) continue; min_type = desc->type; - if (desc->type == DESC_MAT_SOLID) { + if(desc->type == DESC_MAT_SOLID) { found = i; min_delta = desc->d.solid.delta; d3_set(new_pos, projected_pos); } } - if (!min_desc) { + if(!min_desc) { /* Doesn't handle the case where pos is outside the model * but very close and could be considered on the boundary */ fprintf(stderr, "The probe is outside the model.\n"); return RES_BAD_ARG; } - if (on_interface) { + if(on_interface) { *iprim = found; - if (min_type == DESC_MAT_FLUID) + if(min_type == DESC_MAT_FLUID) fprintf(stderr, "Probe was on a fluid.\n"); fprintf(stderr, "The probe is moved to %g,%g,%g (primitive %llu)\n", SPLIT3(new_pos), (long long int)*iprim); - if (min_type == DESC_MAT_SOLID && min_d > min_delta) + if(min_type == DESC_MAT_SOLID && min_d > min_delta) fprintf(stderr, "WARNING: The probe move is %.1f delta long.\n", min_d / min_delta); d3_set(pos, new_pos); d2_set(uv, min_uv); } else { - if (min_type == DESC_MAT_FLUID) + if(min_type == DESC_MAT_FLUID) fprintf(stderr, "The probe is in a fluid: computing fluid temperature, not using a specific position.\n"); - if (min_type == DESC_MAT_SOLID && min_d < min_delta) { + if(min_type == DESC_MAT_SOLID && min_d < min_delta) { fprintf(stderr, "WARNING: The probe is very close to the primitive %llu (%.1f delta).\n", (long long int)found, min_d / min_delta); @@ -181,76 +146,80 @@ select_probe_type } } +end: return res; +error: + goto end; } static res_T create_holder - (struct counts* counts, + (struct mem_allocator* allocator, + struct counts* counts, struct dummies* dummies, struct description* description, struct sdis_device* dev, - struct sdis_medium*** media, + struct darray_media_ptr* media_ptr, const int is_green) { res_T res = RES_OK; + ASSERT(allocator && counts && dummies && description && dev && media_ptr); + switch (description->type) { case DESC_BOUND_H_FOR_SOLID: counts->hbound_count++; /* Create an external fluid with bound temp as fluid temp */ - res = create_fluid(dev, NULL, 1, 1, is_green, 1, -1, - description->d.h_boundary.imposed_temperature, media, - &description->d.h_boundary.mat_id); - if (res != RES_OK) goto error; + ERR(create_fluid(dev, allocator, NULL, 1, 1, is_green, 1, -1, + description->d.h_boundary.imposed_temperature, media_ptr, + &description->d.h_boundary.mat_id)); break; case DESC_BOUND_H_FOR_FLUID: /* Create an external solid with bound temp as solid temp */ counts->hbound_count++; - res = create_solid(dev, NULL, INF, 1, 1, 1, is_green, 1, -1, - description->d.h_boundary.imposed_temperature, 0, media, - &description->d.h_boundary.mat_id); - if (res != RES_OK) goto error; + ERR(create_solid(dev, allocator, NULL, INF, 1, 1, 1, is_green, 1, -1, + description->d.h_boundary.imposed_temperature, 0, media_ptr, + &description->d.h_boundary.mat_id)); break; case DESC_BOUND_T_FOR_SOLID: counts->tbound_count++; - if (dummies->dummy_fluid) { + if(dummies->dummy_fluid) { /* Reuse external dummy fluid */ description->d.t_boundary.mat_id = dummies->dummy_fluid_id; } else { /* Create dummy fluid */ - res = create_fluid(dev, NULL, 1, 1, is_green, 1, -1, - -1, media, &description->d.t_boundary.mat_id); - if (res != RES_OK) goto error; - dummies->dummy_fluid = sa_last(*media); + size_t sz = darray_media_ptr_size_get(media_ptr); + ERR(create_fluid(dev, allocator, NULL, 1, 1, is_green, 1, -1, + -1, media_ptr, &description->d.t_boundary.mat_id)); + dummies->dummy_fluid = darray_media_ptr_data_get(media_ptr)[sz]; dummies->dummy_fluid_id = description->d.t_boundary.mat_id; } break; case DESC_BOUND_T_FOR_FLUID: counts->tbound_count++; - if (dummies->dummy_solid) { + if(dummies->dummy_solid) { /* Reuse external dummy solid */ description->d.t_boundary.mat_id = dummies->dummy_solid_id; } else { /* Create dummy solid */ - res = create_solid(dev, NULL, 1, 1, 1, 1, is_green, 1, -1, - -1, 0, media, &description->d.t_boundary.mat_id); - if (res != RES_OK) goto error; - dummies->dummy_solid = sa_last(*media); + size_t sz = darray_media_ptr_size_get(media_ptr); + ERR(create_solid(dev, allocator, NULL, 1, 1, 1, 1, is_green, 1, -1, + -1, 0, media_ptr, &description->d.t_boundary.mat_id)); + dummies->dummy_solid = darray_media_ptr_data_get(media_ptr)[sz]; dummies->dummy_solid_id = description->d.t_boundary.mat_id; } break; case DESC_BOUND_F_FOR_SOLID: counts->fbound_count++; - if (dummies->dummy_fluid) { + if(dummies->dummy_fluid) { /* Reuse external dummy fluid */ description->d.f_boundary.mat_id = dummies->dummy_fluid_id; } else { /* Create dummy fluid */ - res = create_fluid(dev, NULL, 1, 1, is_green, 1, -1, -1, - media, &description->d.f_boundary.mat_id); - if (res != RES_OK) goto error; - dummies->dummy_fluid = sa_last(*media); + size_t sz = darray_media_ptr_size_get(media_ptr); + ERR(create_fluid(dev, allocator, NULL, 1, 1, is_green, 1, -1, -1, + media_ptr, &description->d.f_boundary.mat_id)); + dummies->dummy_fluid = darray_media_ptr_data_get(media_ptr)[sz]; dummies->dummy_fluid_id = description->d.f_boundary.mat_id; } break; @@ -263,8 +232,9 @@ create_holder break; case DESC_MAT_SOLID: counts->smed_count++; - res = create_solid(dev, - description->d.solid.name, + ERR(create_solid(dev, + allocator, + &description->d.solid.name, description->d.solid.lambda, description->d.solid.rho, description->d.solid.cp, @@ -274,23 +244,22 @@ create_holder description->d.solid.tinit, description->d.solid.imposed_temperature, description->d.solid.vpower, - media, - &description->d.solid.solid_id); - if (res != RES_OK) goto error; + media_ptr, + &description->d.solid.solid_id)); break; case DESC_MAT_FLUID: counts->fmed_count++; - res = create_fluid(dev, - description->d.fluid.name, + ERR(create_fluid(dev, + allocator, + &description->d.fluid.name, description->d.fluid.rho, description->d.fluid.cp, is_green, 0, description->d.fluid.tinit, description->d.fluid.imposed_temperature, - media, - &description->d.fluid.fluid_id); - if (res != RES_OK) goto error; + media_ptr, + &description->d.fluid.fluid_id)); break; default: @@ -319,41 +288,35 @@ compute_probe ASSERT(scn && counts && stardis && (mode & PROBE_COMPUTE)); d3_set(pos, stardis->probe); - res = select_probe_type(scn, stardis, 0, pos, &iprim, uv); - if (res != RES_OK) goto error; + ERR(select_probe_type(scn, stardis, 0, pos, &iprim, uv)); - if (mode & GREEN_MODE) { + if(mode & GREEN_MODE) { ASSERT(iprim == SIZE_MAX); - res = sdis_solve_probe_green_function(scn, - stardis->N, + ERR(sdis_solve_probe_green_function(scn, + stardis->samples, pos, stardis->scale_factor, stardis->radiative_temp[0], stardis->radiative_temp[1], - &green); - if (res != RES_OK) goto error; - res = dump_green(green, counts, stardis->descriptions, - stardis->radiative_temp); - if (res != RES_OK) goto error; + &green)); + ERR(dump_green(green, counts, &stardis->descriptions, stardis->radiative_temp)); } else { ASSERT(iprim == SIZE_MAX); d2_splat(time, stardis->probe[3]); - res = sdis_solve_probe(scn, - stardis->N, + ERR(sdis_solve_probe(scn, + stardis->samples, pos, time, stardis->scale_factor, stardis->radiative_temp[0], stardis->radiative_temp[1], stardis->dump_paths, - &estimator); - if (res != RES_OK) goto error; - res = print_single_MC_result(mode, estimator, stardis); - if (res != RES_OK) goto error; + &estimator)); + ERR(print_single_MC_result(mode, estimator, stardis)); } end: - if (estimator) SDIS(estimator_ref_put(estimator)); - if (green) SDIS(green_function_ref_put(green)); + if(estimator) SDIS(estimator_ref_put(estimator)); + if(green) SDIS(green_function_ref_put(green)); return res; error: goto end; @@ -374,29 +337,25 @@ compute_probe_on_interface ASSERT(scn && counts && stardis && (mode & PROBE_COMPUTE_ON_INTERFACE)); d3_set(pos, stardis->probe); - res = select_probe_type(scn, stardis, 1, pos, &iprim, uv); - if (res != RES_OK) goto error; + ERR(select_probe_type(scn, stardis, 1, pos, &iprim, uv)); - if (mode & GREEN_MODE) { + if(mode & GREEN_MODE) { ASSERT(iprim == SIZE_MAX); - res = sdis_solve_probe_boundary_green_function(scn, - stardis->N, + ERR(sdis_solve_probe_boundary_green_function(scn, + stardis->samples, iprim, uv, SDIS_FRONT, stardis->scale_factor, stardis->radiative_temp[0], stardis->radiative_temp[1], - &green); - if (res != RES_OK) goto error; - res = dump_green(green, counts, stardis->descriptions, - stardis->radiative_temp); - if (res != RES_OK) goto error; + &green)); + ERR(dump_green(green, counts, &stardis->descriptions, stardis->radiative_temp)); } else { ASSERT(iprim == SIZE_MAX); d2_splat(time, stardis->probe[3]); - res = sdis_solve_probe_boundary(scn, - stardis->N, + ERR(sdis_solve_probe_boundary(scn, + stardis->samples, iprim, uv, time, @@ -405,14 +364,12 @@ compute_probe_on_interface stardis->radiative_temp[0], stardis->radiative_temp[1], stardis->dump_paths, - &estimator); - if (res != RES_OK) goto error; - res = print_single_MC_result(mode, estimator, stardis); - if (res != RES_OK) goto error; + &estimator)); + ERR(print_single_MC_result(mode, estimator, stardis)); } end: - if (estimator) SDIS(estimator_ref_put(estimator)); - if (green) SDIS(green_function_ref_put(green)); + if(estimator) SDIS(estimator_ref_put(estimator)); + if(green) SDIS(green_function_ref_put(green)); return res; error: goto end; @@ -437,21 +394,17 @@ compute_camera width = (size_t)stardis->camera.img[0]; height = (size_t)stardis->camera.img[1]; /* Setup the camera */ - res = sdis_camera_create(dev, &cam); - if (res != RES_OK) goto error; - res = sdis_camera_set_proj_ratio(cam, (double)width / (double)height); - if (res != RES_OK) goto error; - res = sdis_camera_set_fov(cam, MDEG2RAD(stardis->camera.fov)); - if (res != RES_OK) goto error; - res = sdis_camera_look_at(cam, + ERR(sdis_camera_create(dev, &cam)); + ERR(sdis_camera_set_proj_ratio(cam, (double)width / (double)height)); + ERR(sdis_camera_set_fov(cam, MDEG2RAD(stardis->camera.fov))); + ERR(sdis_camera_look_at(cam, stardis->camera.pos, stardis->camera.tgt, - stardis->camera.up); - if (res != RES_OK) goto error; + stardis->camera.up)); /* Launch the simulation */ - time[0] = time[1] = stardis->camera.u.time; - res = sdis_solve_camera(scn, + time[0] = time[1] = stardis->camera.time; + ERR(sdis_solve_camera(scn, cam, time, stardis->scale_factor, @@ -461,16 +414,14 @@ compute_camera height, (size_t)stardis->camera.spp, stardis->dump_paths, - &buf); - if (res != RES_OK) goto error; + &buf)); /* Write the image */ - res = dump_image(buf); - if (res != RES_OK) goto error; + ERR(dump_image(buf)); end: - if (cam) SDIS(camera_ref_put(cam)); - if (buf) SDIS(estimator_buffer_ref_put(buf)); + if(cam) SDIS(camera_ref_put(cam)); + if(buf) SDIS(estimator_buffer_ref_put(buf)); return res; error: goto end; @@ -479,7 +430,7 @@ error: static res_T compute_medium (struct sdis_scene* scn, - struct sdis_medium** media, + struct darray_media_ptr* media_ptr, struct stardis* stardis, const int mode) { @@ -490,62 +441,161 @@ compute_medium struct sdis_estimator* estimator = NULL; struct sdis_green_function* green = NULL; - ASSERT(scn && media && stardis && (mode & MEDIUM_COMPUTE)); + ASSERT(scn && media_ptr && stardis && (mode & MEDIUM_COMPUTE)); /* Find medium */ - sz = sa_size(stardis->descriptions); + sz = darray_descriptions_size_get(&stardis->descriptions); FOR_EACH(ii, 0, sz) { - struct description* desc = stardis->descriptions + ii; - if (desc->type == DESC_MAT_SOLID) { - if (0 == strcmp(stardis->solve_name, desc->d.solid.name)) { - ASSERT(sa_size(media) > ii); - medium = media[ii]; + struct description* + desc = darray_descriptions_data_get(&stardis->descriptions) + ii; + if(desc->type == DESC_MAT_SOLID) { + if(str_eq(&stardis->solve_name, &desc->d.solid.name)) { + ASSERT(darray_media_ptr_size_get(media_ptr) > ii); + medium = darray_media_ptr_data_get(media_ptr)[ii]; break; } } - else if (desc->type == DESC_MAT_FLUID) { - if (0 == strcmp(stardis->solve_name, desc->d.fluid.name)) { - ASSERT(sa_size(media) > ii); - medium = media[ii]; + else if(desc->type == DESC_MAT_FLUID) { + if(str_eq(&stardis->solve_name, &desc->d.fluid.name) ){ + ASSERT(darray_media_ptr_size_get(media_ptr) > ii); + medium = darray_media_ptr_data_get(media_ptr)[ii]; break; } } } - if (medium == NULL) { /* Not found */ + if(medium == NULL) { /* Not found */ fprintf(stderr, "Cannot solve medium %s (unknown medium)\n", - stardis->solve_name); + str_cget(&stardis->solve_name)); res = RES_BAD_ARG; goto error; } - if (mode & GREEN_MODE) { - res = sdis_solve_medium_green_function(scn, - stardis->N, + if(mode & GREEN_MODE) { + ERR(sdis_solve_medium_green_function(scn, + stardis->samples, medium, stardis->scale_factor, stardis->radiative_temp[0], stardis->radiative_temp[1], - &green); - if (res != RES_OK) goto error; + &green)); } else { d2_splat(time, stardis->probe[3]); time[0] = time[1] = stardis->probe[3]; - res = sdis_solve_medium(scn, - stardis->N, + ERR(sdis_solve_medium(scn, + stardis->samples, medium, time, stardis->scale_factor, stardis->radiative_temp[0], stardis->radiative_temp[1], stardis->dump_paths, - &estimator); - if (res != RES_OK) goto error; - res = print_single_MC_result(mode, estimator, stardis); - if (res != RES_OK) goto error; + &estimator)); + ERR(print_single_MC_result(mode, estimator, stardis)); } end: - if (estimator) SDIS(estimator_ref_put(estimator)); - if (green) SDIS(green_function_ref_put(green)); + if(estimator) SDIS(estimator_ref_put(estimator)); + if(green) SDIS(green_function_ref_put(green)); + return res; +error: + goto end; +} + +static res_T +add_compute_region_triangle + (const unsigned unique_id, const unsigned itri, void* context) +{ + res_T res = RES_OK; + const struct add_geom_ctx* ctx = context; + struct compute_region* region; + ASSERT(ctx); (void)itri; + region = ctx->custom; + /* The triangle should be already in the model, but is not + * Keep it in err_triangles to allow dumps */ + ERR(darray_uint_push_back(&region->err_triangles, &unique_id)); + +end: + return res; +error: + goto end; +} + +static res_T +merge_compute_region_triangle + (const unsigned unique_id, + const unsigned itri, + const int reversed_triangle, + unsigned triangle_properties[SG3D_PROP_TYPES_COUNT__], + const unsigned merged_properties[SG3D_PROP_TYPES_COUNT__], + void* context, + int* merge_conflict_status) +{ + res_T res = RES_OK; + const struct add_geom_ctx* ctx = context; + struct compute_region* region; + size_t uid; + enum sdis_side side = reversed_triangle ? SDIS_BACK : SDIS_FRONT; + ASSERT(ctx); + (void)itri; (void)reversed_triangle; (void)triangle_properties; + (void)merged_properties; (void)merge_conflict_status; + region = ctx->custom; + /* Build the compute region */ + uid = unique_id; + ERR(darray_size_t_push_back(&region->primitives, &uid)); + ERR(darray_sides_push_back(&region->sides, &side)); +end: + return res; +error: + goto end; +} + +/* Process an STL file describing a compute region * + * Triangles must be members of the geometry already + * TODO: what if same triangle appears more than once (same side or not)? + */ +res_T +read_compute_region(struct stardis* stardis) +{ + res_T res = RES_OK; + struct sstl* sstl = NULL; + struct add_geom_ctx add_geom_ctx; + struct sg3d_geometry_add_callbacks callbacks = SG3D_ADD_CALLBACKS_NULL__; + int new_triangles_initialized = 0; + const char* file; + + ASSERT(stardis); + + file = str_cget(&stardis->solve_name); + callbacks.get_indices = add_geom_ctx_indices; + callbacks.get_position = add_geom_ctx_position; + callbacks.add_triangle = add_compute_region_triangle; + callbacks.merge_triangle = merge_compute_region_triangle; + + ERR(sstl_create(NULL, stardis->allocator, 0, &sstl)); + res = sstl_load(sstl, file); + if(res != RES_OK) { + fprintf(stderr, "Cannot read STL file: %s\n", file); + goto error; + } + ERR(sstl_get_desc(sstl, &add_geom_ctx.stl_desc)); + ASSERT(add_geom_ctx.stl_desc.vertices_count <= UINT_MAX + && add_geom_ctx.stl_desc.triangles_count <= UINT_MAX); + + new_triangles_initialized = 1; + add_geom_ctx.custom = &stardis->compute_region; + + res = sg3d_geometry_add( + stardis->geometry.sg3d, + (unsigned)add_geom_ctx.stl_desc.vertices_count, + (unsigned)add_geom_ctx.stl_desc.triangles_count, + &callbacks, + &add_geom_ctx); + if(res != RES_OK) { + fprintf(stderr, "Cannot add file content: %s\n", file); + goto error; + } + +end: + if(sstl) SSTL(ref_put(sstl)); return res; error: goto end; @@ -561,54 +611,41 @@ compute_boundary res_T res = RES_OK; double time[2], pos[3]; struct sdis_green_function* green = NULL; - struct sdis_estimator* estimator = NULL; + struct sdis_estimator* estimator = NULL; ASSERT(scn && counts && stardis && (mode & BOUNDARY_COMPUTE)); - /* FIXME: MOVE TO BOUNDARY SETTING */ - ASSERT(sa_size(stardis->boundary.primitives) - == sa_size(stardis->boundary.sides)); - if (sa_size(stardis->boundary.primitives) == 0) { - fprintf(stderr, "Cannot solve boundary %s (empty geometry)\n", - stardis->solve_name); - } - d3_set(pos, stardis->probe); - if (mode & GREEN_MODE) { - res = sdis_solve_boundary_green_function(scn, - stardis->N, - stardis->boundary.primitives, - stardis->boundary.sides, - sa_size(stardis->boundary.primitives), + if(mode & GREEN_MODE) { + ERR(sdis_solve_boundary_green_function(scn, + stardis->samples, + darray_size_t_cdata_get(&stardis->compute_region.primitives), + darray_sides_cdata_get(&stardis->compute_region.sides), + darray_size_t_size_get(&stardis->compute_region.primitives), stardis->scale_factor, stardis->radiative_temp[0], stardis->radiative_temp[1], - &green); - if (res != RES_OK) goto error; - res = dump_green(green, counts, stardis->descriptions, - stardis->radiative_temp); - if (res != RES_OK) goto error; + &green)); + ERR(dump_green(green, counts, &stardis->descriptions, stardis->radiative_temp)); } else { d2_splat(time, stardis->probe[3]); - res = sdis_solve_boundary(scn, - stardis->N, - stardis->boundary.primitives, - stardis->boundary.sides, - sa_size(stardis->boundary.primitives), + ERR(sdis_solve_boundary(scn, + stardis->samples, + darray_size_t_cdata_get(&stardis->compute_region.primitives), + darray_sides_cdata_get(&stardis->compute_region.sides), + darray_size_t_size_get(&stardis->compute_region.primitives), time, stardis->scale_factor, stardis->radiative_temp[0], stardis->radiative_temp[1], stardis->dump_paths, - &estimator); - if (res != RES_OK) goto error; - res = print_single_MC_result(mode, estimator, stardis); - if (res != RES_OK) goto error; + &estimator)); + ERR(print_single_MC_result(mode, estimator, stardis)); } end: - if (estimator) SDIS(estimator_ref_put(estimator)); - if (green) SDIS(green_function_ref_put(green)); + if(estimator) SDIS(estimator_ref_put(estimator)); + if(green) SDIS(green_function_ref_put(green)); return res; error: goto end; @@ -627,29 +664,21 @@ compute_flux_boundary ASSERT(scn && stardis && (mode & FLUX_BOUNDARY_COMPUTE)); - /* FIXME: MOVE TO BOUNDARY SETTING */ - if (sa_size(stardis->boundary.primitives) == 0) { - fprintf(stderr, "Cannot solve boundary %s (empty geometry)\n", - stardis->solve_name); - } - d3_set(pos, stardis->probe); d2_splat(time, stardis->probe[3]); - res = sdis_solve_boundary_flux(scn, - stardis->N, - stardis->boundary.primitives, - sa_size(stardis->boundary.primitives), + ERR(sdis_solve_boundary_flux(scn, + stardis->samples, + darray_size_t_cdata_get(&stardis->compute_region.primitives), + darray_size_t_size_get(&stardis->compute_region.primitives), time, stardis->scale_factor, stardis->radiative_temp[0], stardis->radiative_temp[1], - &estimator); - if (res != RES_OK) goto error; - res = print_single_MC_result(mode, estimator, stardis); - if (res != RES_OK) goto error; + &estimator)); + ERR(print_single_MC_result(mode, estimator, stardis)); end: - if (estimator) SDIS(estimator_ref_put(estimator)); - if (green) SDIS(green_function_ref_put(green)); + if(estimator) SDIS(estimator_ref_put(estimator)); + if(green) SDIS(green_function_ref_put(green)); return res; error: goto end; @@ -664,55 +693,62 @@ compute_map res_T res = RES_OK; double time[2], pos[3]; struct sdis_green_function* green = NULL; - struct sdis_estimator** estimators = NULL; + struct darray_estimators estimators; + int estimators_initialized = 0; size_t p; ASSERT(scn && stardis && (mode & MAP_COMPUTE) && !(mode & GREEN_MODE)); (void)mode; - /* FIXME: MOVE TO BOUNDARY SETTING */ - ASSERT(sa_size(stardis->boundary.primitives) - == sa_size(stardis->boundary.sides)); - if (sa_size(stardis->boundary.primitives) == 0) { - fprintf(stderr, "Cannot solve boundary %s (empty geometry)\n", - stardis->solve_name); - } - + darray_estimators_init(stardis->allocator, &estimators); + estimators_initialized = 1; d3_set(pos, stardis->probe); d2_splat(time, stardis->probe[3]); - sa_add(estimators, sa_size(stardis->boundary.primitives)); + ERR(darray_estimators_resize(&estimators, + darray_estimators_size_get(&estimators) + + darray_size_t_size_get(&stardis->compute_region.primitives))); - FOR_EACH(p, 0, sa_size(stardis->boundary.primitives)) { - res = sdis_solve_boundary(scn, - stardis->N, - stardis->boundary.primitives + p, - stardis->boundary.sides + p, + FOR_EACH(p, 0, darray_size_t_size_get(&stardis->compute_region.primitives)) { + ERR(sdis_solve_boundary(scn, + stardis->samples, + darray_size_t_cdata_get(&stardis->compute_region.primitives) + p, + darray_sides_cdata_get(&stardis->compute_region.sides) + p, 1, time, stardis->scale_factor, stardis->radiative_temp[0], stardis->radiative_temp[1], stardis->dump_paths, - estimators + p); - if (res != RES_OK) goto error; + darray_estimators_data_get(&estimators) + p)); } - dump_map(stardis, estimators); + dump_map(stardis, &estimators); end: - if (estimators) { - FOR_EACH(p, 0, sa_size(stardis->boundary.primitives)) { - SDIS(estimator_ref_put(estimators[p])); + if(estimators_initialized) { + struct sdis_estimator** est = darray_estimators_data_get(&estimators); + FOR_EACH(p, 0, darray_size_t_size_get(&stardis->compute_region.primitives)) { + SDIS(estimator_ref_put(est[p])); } - sa_release(estimators); + darray_estimators_release(&estimators); } - if (green) SDIS(green_function_ref_put(green)); + if(green) SDIS(green_function_ref_put(green)); return res; error: goto end; } +static struct sdis_interface* +geometry_get_interface + (const size_t itri, void* ctx) +{ + struct darray_interface_ptrs* app_interface_data = ctx; + ASSERT(app_interface_data + && itri < darray_interface_ptrs_size_get(app_interface_data)); + return darray_interface_ptrs_cdata_get(app_interface_data)[itri]; +} + /******************************************************************************* * Public Functions ******************************************************************************/ @@ -724,83 +760,80 @@ stardis_compute { res_T res = RES_OK; struct sdis_device* dev = NULL; - struct sdis_medium** media = NULL; - struct sdis_scene* scn = NULL; + struct darray_media_ptr media; + struct sg3d_sdisXd_scene_create_context create_context; struct dummies dummies = DUMMIES_NULL; struct htable_intface htable_interfaces; - int htable_interfaces_initialised = 0; + int struct_initialised = 0; struct counts counts = COUNTS_NULL; - unsigned i = 0; - int verbose = 1; + int try_compute = 0, try_scene = 0; + unsigned i, vcount, tcount; ASSERT(stardis && (mode & COMPUTE_MODES)); - res = sdis_device_create(NULL, &stardis->allocator, stardis->nthreads, verbose, &dev); - if (res != RES_OK) goto error; + ERR(sdis_device_create(NULL, stardis->allocator, stardis->nthreads, + stardis->verbose, &dev)); + + htable_intface_init(stardis->allocator, &htable_interfaces); + darray_media_ptr_init(stardis->allocator, &media); + struct_initialised = 1; /* Create media and property holders found in descriptions */ - for (i = 0; i < sa_size(stardis->descriptions); ++i) { - res = create_holder(&counts, &dummies, stardis->descriptions + i, - dev, &media, mode & GREEN_MODE); - if (res != RES_OK) goto error; + for(i = 0; i < darray_descriptions_size_get(&stardis->descriptions); ++i) { + ERR(create_holder(stardis->allocator, &counts, &dummies, + darray_descriptions_data_get(&stardis->descriptions) + i, dev, + &media, mode & GREEN_MODE)); } /* Create interfaces */ - htable_intface_init(&stardis->allocator, &htable_interfaces); - htable_interfaces_initialised = 1; - for (i = 0; i < sa_size(stardis->geometry.triangle); ++i) { - res = create_intface(dev, i, &htable_interfaces, &stardis->geometry, - stardis->descriptions, media); - if (res != RES_OK) goto error; - } + ERR(sg3d_geometry_get_unique_vertices_count(stardis->geometry.sg3d, &vcount)); + ERR(sg3d_geometry_get_unique_triangles_count(stardis->geometry.sg3d, &tcount)); - /* Create scene */ - res = sdis_scene_create(dev, - sa_size(stardis->geometry.triangle), - geometry_get_indices, geometry_get_interface, - sa_size(stardis->geometry.vertex), - geometry_get_position, &stardis->geometry, &scn); - if(res != RES_OK) { - fprintf(stderr, "%s: could not setup the scene!\n", FUNC_NAME); - goto error; + for(i = 0; i < tcount; ++i) { + ERR(create_intface(dev, i, &htable_interfaces, &stardis->geometry, + &stardis->descriptions, &media)); } - /* If the scene has holes dump them */ - res = create_edge_file_if(scn, &stardis->allocator); - if (res != RES_OK) goto error; + /* Create scene */ + try_scene = 1; + create_context.geometry = stardis->geometry.sg3d; + create_context.app_interface_getter = geometry_get_interface; + create_context.app_interface_data = &stardis->geometry.interf_bytrg; + ERR(sdis_scene_create(dev, tcount, sg3d_sdisXd_geometry_get_indices, + sg3d_sdisXd_geometry_get_interface, vcount, sg3d_sdisXd_geometry_get_position, + &create_context, &stardis->geometry.sdis_scn)); /* Compute */ - if (mode & PROBE_COMPUTE) - res = compute_probe(scn, &counts, stardis, mode); - else if (mode & PROBE_COMPUTE_ON_INTERFACE) - res = compute_probe_on_interface(scn, &counts, stardis, mode); - else if (mode & IR_COMPUTE) - res = compute_camera(dev, scn, stardis, mode); - else if (mode & MEDIUM_COMPUTE) - res = compute_medium(scn, media, stardis, mode); - else if (mode & BOUNDARY_COMPUTE) - res = compute_boundary(scn, &counts, stardis, mode); - else if (mode & FLUX_BOUNDARY_COMPUTE) - res = compute_flux_boundary(scn, stardis, mode); - else if (mode & MAP_COMPUTE) - res = compute_map(scn, stardis, mode); + try_compute = 1; + if(mode & PROBE_COMPUTE) + ERR(compute_probe(stardis->geometry.sdis_scn, &counts, stardis, mode)); + else if(mode & PROBE_COMPUTE_ON_INTERFACE) + ERR(compute_probe_on_interface(stardis->geometry.sdis_scn, + &counts, stardis, mode)); + else if(mode & IR_COMPUTE) + ERR(compute_camera(dev, stardis->geometry.sdis_scn, stardis, mode)); + else if(mode & MEDIUM_COMPUTE) + ERR(compute_medium(stardis->geometry.sdis_scn, &media, stardis, mode)); + else if(mode & BOUNDARY_COMPUTE) + ERR(compute_boundary(stardis->geometry.sdis_scn, &counts, stardis, mode)); + else if(mode & FLUX_BOUNDARY_COMPUTE) + ERR(compute_flux_boundary(stardis->geometry.sdis_scn, stardis, mode)); + else if(mode & MAP_COMPUTE) + ERR(compute_map(stardis->geometry.sdis_scn, stardis, mode)); else FATAL("Unknown mode.\n"); - if(res != RES_OK) { - fprintf(stderr, "%s: computation failed!\n", FUNC_NAME); - goto error; - } end: - if (htable_interfaces_initialised) + if(struct_initialised) { htable_intface_release(&htable_interfaces); - for (i = 0; i < sa_size(media); ++i) - SDIS(medium_ref_put(media[i])); - sa_release(media); - - if (scn) SDIS(scene_ref_put(scn)); - if (dev) SDIS(device_ref_put(dev)); + FOR_EACH(i, 0, darray_media_ptr_size_get(&media)) + SDIS(medium_ref_put(darray_media_ptr_data_get(&media)[i])); + darray_media_ptr_release(&media); + } + if(dev) SDIS(device_ref_put(dev)); return res; error: + if(try_compute) fprintf(stderr, "%s: computation failed!\n", FUNC_NAME); + else if(try_scene) fprintf(stderr, "%s: could not setup the scene!\n", FUNC_NAME); goto end; } diff --git a/src/stardis-compute.h b/src/stardis-compute.h @@ -12,6 +12,15 @@ #include <rsys/hash_table.h> struct stardis; +struct sdis_medium; + +#define DARRAY_NAME media_ptr +#define DARRAY_DATA struct sdis_medium* +#include <rsys/dynamic_array.h> + +#define DARRAY_NAME estimators +#define DARRAY_DATA struct sdis_estimator* +#include <rsys/dynamic_array.h> struct counts { size_t smed_count, fmed_count, tbound_count, hbound_count, @@ -22,9 +31,12 @@ struct counts { } static const struct counts COUNTS_NULL = COUNTS_NULL__; -extern res_T +extern LOCAL_SYM res_T stardis_compute (struct stardis* stardis, const int mode); +extern LOCAL_SYM res_T +read_compute_region(struct stardis* stardis); + #endif \ No newline at end of file diff --git a/src/stardis-fluid.c b/src/stardis-fluid.c @@ -1,8 +1,8 @@ #include "stardis-fluid.h" #include "stardis-compute.h" +#include "stardis-app.h" #include <sdis.h> -#include<rsys/stretchy_array.h> #include <limits.h> @@ -36,23 +36,24 @@ fluid_get_temperature struct sdis_data* data) { const struct fluid* fluid_props = sdis_data_cget(data); - if (fluid_props->is_green || vtx->time > fluid_props->t0) { + if(fluid_props->is_green || vtx->time > fluid_props->t0) { /* Always use temp for Green mode, regardless of time */ - if (fluid_props->imposed_temperature >= 0) + if(fluid_props->imposed_temperature >= 0) return fluid_props->imposed_temperature; else return -1; } /* Time is <= 0: use tinit */ - if (fluid_props->tinit >= 0) return fluid_props->tinit; + if(fluid_props->tinit >= 0) return fluid_props->tinit; /* Must have had tinit defined: error! */ - if (fluid_props->name[0]) { + if(str_is_empty(&fluid_props->name)) { + FATAL("fluid_get_temperature: getting undefined Tinit\n"); + } else { fprintf(stderr, "fluid_get_temperature: getting undefined Tinit (fluid '%s')\n", - fluid_props->name); + str_cget(&fluid_props->name)); ASSERT(0); abort(); - } else - FATAL("fluid_get_temperature: getting undefined Tinit\n"); + } } /******************************************************************************* @@ -62,14 +63,15 @@ fluid_get_temperature res_T create_fluid (struct sdis_device* dev, - const char* name, + struct mem_allocator* allocator, + const struct str* name, const double rho, const double cp, const int is_green, const int is_outside, const double tinit, const double imposed_temperature, - struct sdis_medium*** media_ptr, + struct darray_media_ptr* media_ptr, unsigned* out_id) { res_T res = RES_OK; @@ -83,15 +85,14 @@ create_fluid fluid_shader.calorific_capacity = fluid_get_calorific_capacity; fluid_shader.volumic_mass = fluid_get_volumic_mass; fluid_shader.temperature = fluid_get_temperature; - res = sdis_data_create(dev, sizeof(struct fluid), ALIGNOF(struct fluid), - NULL, &data); - if (res != RES_OK) goto error; + ERR(sdis_data_create(dev, sizeof(struct fluid), ALIGNOF(struct fluid), + NULL, &data)); - sz = sa_size(*media_ptr); + sz = darray_media_ptr_size_get(media_ptr); ASSERT(sz < INT_MAX); fluid_props = sdis_data_get(data); /* Fetch the allocated memory space */ - if (name) strncpy(fluid_props->name, name, sizeof(fluid_props->name)); - else fluid_props->name[0] = '\0'; + str_init(allocator, &fluid_props->name); + if(name) str_copy(&fluid_props->name, name); fluid_props->cp = cp; fluid_props->rho = rho; fluid_props->tinit = tinit; @@ -99,13 +100,14 @@ create_fluid fluid_props->is_green = is_green; fluid_props->is_outside = is_outside; fluid_props->id = (unsigned)sz; - - res = sdis_fluid_create(dev, &fluid_shader, data, sa_add(*media_ptr, 1)); - if (res != RES_OK) goto error; + + ERR(darray_media_ptr_resize(media_ptr, sz+1)); + ERR(sdis_fluid_create(dev, &fluid_shader, data, + darray_media_ptr_data_get(media_ptr) + sz)); *out_id = fluid_props->id; end: - if (data) SDIS(data_ref_put(data)); + if(data) SDIS(data_ref_put(data)); return res; error: goto end; diff --git a/src/stardis-fluid.h b/src/stardis-fluid.h @@ -4,16 +4,18 @@ #define SDIS_FLUID_H #include <sdis.h> + #include <rsys/rsys.h> +#include <rsys/str.h> struct sdis_device; -struct sdis_medium; +struct darray_media_ptr; /******************************************************************************* * Fluid data ******************************************************************************/ struct fluid { - char name[32]; + struct str name; double cp; /* Calorific capacity */ double rho; /* Volumic mass */ /* Compute mode */ @@ -28,14 +30,15 @@ struct fluid { extern res_T create_fluid (struct sdis_device* dev, - const char* name, + struct mem_allocator* allocator, + const struct str* name, const double rho, const double cp, const int is_green, const int is_outside, const double tinit, const double imposed_temperature, - struct sdis_medium*** media_ptr, + struct darray_media_ptr* media_ptr, unsigned* out_id); #endif diff --git a/src/stardis-intface.c b/src/stardis-intface.c @@ -1,9 +1,9 @@ #include "stardis-intface.h" #include "stardis-app.h" +#include "stardis-compute.h" #include "stardis-output.h" #include <sdis.h> -#include <rsys/stretchy_array.h> /******************************************************************************* * Local Functions @@ -69,10 +69,9 @@ create_intface unsigned tr_idx, struct htable_intface* htable_interfaces, struct geometry* geometry, - const struct description* descriptions, - struct sdis_medium* const* media) + const struct darray_descriptions* desc_array, + struct darray_media_ptr* med_array) { - struct triangle *trg; struct int_descs int_descs = INT_DESCS_NULL; struct sdis_interface** p_intface; struct sdis_interface* intface; @@ -80,26 +79,29 @@ create_intface struct sdis_medium* back_med = NULL; struct sdis_interface_side_shader* fluid_side_shader = NULL; struct sdis_data* data = NULL; - unsigned fd, bd, cd; + unsigned fd, bd, cd, properties[3]; int fluid_count = 0, solid_count = 0; int solid_fluid_connection_count = 0, connection_count = 0, boundary_count = 0; + const struct description* descriptions; + struct sdis_medium** media; res_T res = RES_OK; - ASSERT(dev && htable_interfaces && geometry && descriptions && media); + ASSERT(dev && htable_interfaces && geometry && desc_array && med_array); - trg = geometry->triangle + tr_idx; - fd = trg->front_description; - bd = trg->back_description; - cd = trg->connection_description; + SG3D(geometry_get_unique_triangle_properties(geometry->sg3d, tr_idx, + properties)); + + descriptions = darray_descriptions_cdata_get(desc_array); + media = darray_media_ptr_data_get(med_array); /* Create key */ - int_descs.front = fd; - int_descs.back = bd; - int_descs.connect = trg->connection_description; + int_descs.front = fd = properties[SG3D_FRONT]; + int_descs.back = bd = properties[SG3D_BACK]; + int_descs.connect = cd = properties[SG3D_INTFACE]; /* Search if interface already exists, or create it */ p_intface = htable_intface_find(htable_interfaces, &int_descs); - if (p_intface) { + if(p_intface) { intface = *p_intface; } else { /* create new interface and register */ @@ -118,7 +120,7 @@ create_intface interface_props->front_boundary_id = UINT_MAX; interface_props->back_boundary_id = UINT_MAX; - if (front_defined) { + if(front_defined) { switch (descriptions[fd].type) { case DESC_MAT_SOLID: id = descriptions[fd].d.solid.solid_id; @@ -135,7 +137,7 @@ create_intface FATAL("error:" STR(__FILE__) ":" STR(__LINE__)": Invalid type.\n"); } } - if (back_defined) { + if(back_defined) { switch (descriptions[bd].type) { case DESC_MAT_SOLID: id = descriptions[bd].d.solid.solid_id; @@ -152,13 +154,13 @@ create_intface FATAL("error:" STR(__FILE__) ":" STR(__LINE__)": Invalid type.\n"); } } - if (connect_defined) { + if(connect_defined) { const struct description* connect = descriptions + cd; int type_checked = 0; struct sdis_medium* def_medium = NULL; unsigned ext_id; - if (connect->type == DESC_SOLID_FLUID_CONNECT) { - if (solid_count != 1 || fluid_count != 1) { + if(connect->type == DESC_SOLID_FLUID_CONNECT) { + if(solid_count != 1 || fluid_count != 1) { ASSERT(front_defined && back_defined); fprintf(stderr, "Can only define a DESC_SOLID_FLUID_CONNECT between a fluid and a solid\n"); @@ -166,7 +168,7 @@ create_intface goto error; } } else { - if (front_defined == back_defined) { + if(front_defined == back_defined) { fprintf(stderr, "Cannot define a boundary between 2 %s regions\n", front_defined ? "defined" : "undefined"); @@ -174,13 +176,13 @@ create_intface goto error; } def_medium = front_defined ? front_med : back_med; - if (front_defined) + if(front_defined) interface_props->front_boundary_id = cd; else interface_props->back_boundary_id = cd; } switch (connect->type) { case DESC_BOUND_H_FOR_FLUID: - if (sdis_medium_get_type(def_medium) != SDIS_FLUID) { + if(sdis_medium_get_type(def_medium) != SDIS_FLUID) { fprintf(stderr, "Can only define a DESC_BOUND_H_FOR_FLUID boundary for a fluid region\n"); res = RES_BAD_ARG; @@ -189,7 +191,7 @@ create_intface type_checked = 1; /* fall through */ case DESC_BOUND_H_FOR_SOLID: - if (!type_checked + if(!type_checked && sdis_medium_get_type(def_medium) != SDIS_SOLID) { fprintf(stderr, "Can only define a DESC_BOUND_H_FOR_SOLID boundary for a solid region\n"); @@ -197,18 +199,18 @@ create_intface goto error; } ext_id = connect->d.h_boundary.mat_id; /* External material id */ - ASSERT(ext_id < sa_size(media)); + ASSERT(ext_id < darray_media_ptr_size_get(med_array)); ASSERT(sdis_medium_get_type(media[ext_id]) == (connect->type == DESC_BOUND_H_FOR_SOLID ? SDIS_FLUID : SDIS_SOLID)); connection_count++; boundary_count++; - if (front_defined) back_med = media[ext_id]; + if(front_defined) back_med = media[ext_id]; else front_med = media[ext_id]; interface_shader.convection_coef = interface_get_convection_coef; - interface_shader.convection_coef_upper_bound = connect->d.h_boundary.hc_max; + interface_shader.convection_coef_upper_bound = connect->d.h_boundary.hc; interface_props->hc = connect->d.h_boundary.hc; ASSERT(!fluid_side_shader); /* Cause defined medium is solid */ - if (connect->d.h_boundary.has_emissivity) { + if(connect->d.h_boundary.hc > 0) { fluid_side_shader = front_defined ? &interface_shader.back : &interface_shader.front; fluid_side_shader->emissivity = interface_get_emissivity; @@ -218,7 +220,7 @@ create_intface } break; case DESC_BOUND_T_FOR_FLUID: - if (sdis_medium_get_type(def_medium) != SDIS_FLUID) { + if(sdis_medium_get_type(def_medium) != SDIS_FLUID) { fprintf(stderr, "Can only define a DESC_BOUND_T_FOR_FLUID boundary for a fluid region\n"); res = RES_BAD_ARG; @@ -227,7 +229,7 @@ create_intface type_checked = 1; /* fall through */ case DESC_BOUND_T_FOR_SOLID: - if (!type_checked + if(!type_checked && sdis_medium_get_type(def_medium) != SDIS_SOLID) { fprintf(stderr, "Can only define a DESC_BOUND_T_FOR_SOLID boundary for a solid region\n"); @@ -235,11 +237,11 @@ create_intface goto error; } ext_id = connect->d.t_boundary.mat_id; /* External material id */ - ASSERT(ext_id < sa_size(media)); + ASSERT(ext_id < darray_media_ptr_size_get(med_array)); ASSERT(sdis_medium_get_type(media[ext_id]) == SDIS_FLUID); connection_count++; boundary_count++; - if (front_defined) { + if(front_defined) { back_med = media[ext_id]; /* We set the known T inside * TODO: should be outside to allow contact resistances */ @@ -253,15 +255,15 @@ create_intface ASSERT(connect->d.t_boundary.imposed_temperature >= 0); interface_props->imposed_temperature = connect->d.t_boundary.imposed_temperature; - if (connect->d.t_boundary.has_hc) { + if(connect->d.t_boundary.hc > 0) { ASSERT(connect->type == DESC_BOUND_T_FOR_FLUID); interface_shader.convection_coef = interface_get_convection_coef; - interface_shader.convection_coef_upper_bound = connect->d.t_boundary.hc_max; + interface_shader.convection_coef_upper_bound = connect->d.t_boundary.hc; interface_props->hc = connect->d.t_boundary.hc; } break; case DESC_BOUND_F_FOR_SOLID: - if (sdis_medium_get_type(def_medium) != SDIS_SOLID) { + if(sdis_medium_get_type(def_medium) != SDIS_SOLID) { fprintf(stderr, "Can only define a DESC_BOUND_F_FOR_SOLID boundary for a solid region\n"); res = RES_BAD_ARG; @@ -269,7 +271,7 @@ create_intface } connection_count++; boundary_count++; - if (front_defined) { + if(front_defined) { back_med = media[connect->d.f_boundary.mat_id]; interface_shader.front.flux = interface_get_flux; } else { @@ -285,12 +287,12 @@ create_intface connection_count++; solid_fluid_connection_count++; ASSERT(fluid_side_shader); - if (connect->d.sf_connect.has_hc) { + if(connect->d.sf_connect.hc > 0) { interface_shader.convection_coef = interface_get_convection_coef; - interface_shader.convection_coef_upper_bound = connect->d.sf_connect.hc_max; + interface_shader.convection_coef_upper_bound = connect->d.sf_connect.hc; interface_props->hc = connect->d.sf_connect.hc; } - if (connect->d.sf_connect.has_emissivity) { + if(connect->d.sf_connect.emissivity > 0) { fluid_side_shader->emissivity = interface_get_emissivity; fluid_side_shader->specular_fraction = interface_get_alpha; interface_props->emissivity = connect->d.sf_connect.emissivity; @@ -302,7 +304,7 @@ create_intface } } - if ((fluid_count == 2) + if((fluid_count == 2) || (fluid_count + solid_count + connection_count < 2) || (boundary_count ? (fluid_count + solid_count != 1) : (fluid_count + solid_count != 2)) @@ -310,15 +312,14 @@ create_intface { /* Incoherent triangle description */ fprintf(stderr, "Incoherent triangle description (%u)\n", tr_idx); - print_trg_as_obj(stderr, geometry->vertex, trg->indices.data); fprintf(stderr, "Front: "); - if (!front_defined) fprintf(stderr, "undefined\n"); + if(!front_defined) fprintf(stderr, "undefined\n"); else print_description(stderr, descriptions + fd); fprintf(stderr, "Back: "); - if (!back_defined) fprintf(stderr, "undefined\n"); + if(!back_defined) fprintf(stderr, "undefined\n"); else print_description(stderr, descriptions + bd); fprintf(stderr, "Connection: "); - if (!connect_defined) fprintf(stderr, "undefined\n"); + if(!connect_defined) fprintf(stderr, "undefined\n"); else print_description(stderr, descriptions + cd); res = RES_BAD_ARG; goto error; @@ -327,16 +328,15 @@ create_intface res = sdis_interface_create(dev, front_med, back_med, &interface_shader, data, &intface); SDIS(data_ref_put(data)); data = NULL; - if (res != RES_OK) { + if(res != RES_OK) { fprintf(stderr, "Cannot create interface associated to triangle %u\n", tr_idx); goto error; } - sa_push(geometry->interfaces, intface); - res = htable_intface_set(htable_interfaces, &int_descs, &intface); - if (res != RES_OK) goto error; + ERR(darray_interface_ptrs_push_back(&geometry->interfaces, &intface)); + ERR(htable_intface_set(htable_interfaces, &int_descs, &intface)); } - sa_push(geometry->interf_bytrg, intface); + ERR(darray_interface_ptrs_push_back(&geometry->interf_bytrg, &intface)); end: return res; diff --git a/src/stardis-intface.h b/src/stardis-intface.h @@ -12,12 +12,14 @@ struct sdis_device; struct geometry; struct description; struct sdis_medium; +struct darray_descriptions; +struct darray_media_ptr; /******************************************************************************* * Interface data ******************************************************************************/ struct intface { - double hc; /* Convection coefficient */ + double hc; double emissivity; double alpha; /* Imposed compute temperature & flux */ @@ -60,7 +62,7 @@ create_intface unsigned tr_idx, struct htable_intface* htable_interfaces, struct geometry* geometry, - const struct description* descriptions, - struct sdis_medium* const* media); + const struct darray_descriptions* descriptions, + struct darray_media_ptr* media); #endif diff --git a/src/stardis-output.c b/src/stardis-output.c @@ -16,7 +16,7 @@ #include <stdio.h> struct w_ctx { - const struct description* desc; + const struct darray_descriptions* desc; struct htable_weigth weigths; }; @@ -40,7 +40,9 @@ print_power_term struct sdis_data* data = NULL; enum sdis_medium_type type; unsigned id; - ASSERT(mdm && ctx); (void)ctx; + const struct w_ctx* w_ctx = ctx; + + ASSERT(mdm && ctx); data = sdis_medium_get_data(mdm); type = sdis_medium_get_type(mdm); @@ -51,8 +53,9 @@ print_power_term } case SDIS_SOLID: { struct solid* d__ = sdis_data_get(data); + const struct description* descs = darray_descriptions_cdata_get(w_ctx->desc); id = d__->id; - ASSERT(id == ((struct w_ctx*)ctx)->desc[id].d.solid.solid_id); + ASSERT(id == descs[id].d.solid.solid_id); (void)descs; printf("\tS\t%u\t%g", id, power_term); break; } @@ -73,14 +76,16 @@ get_flux_terms struct intface* d__; unsigned id; struct w_ctx* w_ctx = ctx; + const struct description* descs; ASSERT(interf && ctx); data = sdis_interface_get_data(interf); d__ = sdis_data_get(data); id = (side == SDIS_FRONT) ? d__->front_boundary_id : d__->back_boundary_id; + descs = darray_descriptions_cdata_get(w_ctx->desc); - switch (w_ctx->desc[id].type) { + switch (descs[id].type) { case DESC_BOUND_T_FOR_SOLID: case DESC_BOUND_T_FOR_FLUID: case DESC_BOUND_H_FOR_SOLID: @@ -88,18 +93,18 @@ get_flux_terms FATAL("Cannot have a flux term here.\n"); break; case DESC_BOUND_F_FOR_SOLID: { double* w; - ASSERT(id == w_ctx->desc[id].d.f_boundary.mat_id); + ASSERT(id == descs[id].d.f_boundary.mat_id); w = htable_weigth_find(&w_ctx->weigths, &id); - if (w) - *w += flux_term; - else { - res = htable_weigth_set(&w_ctx->weigths, &id, &flux_term); - } + if(w) *w += flux_term; + else ERR(htable_weigth_set(&w_ctx->weigths, &id, &flux_term)); break; } default: FATAL("Unreachable code.\n"); break; } +end: return res; +error: + goto end; } static res_T @@ -123,8 +128,7 @@ dump_path fprintf(stream, "POINTS %zu double\n", vcount); FOR_EACH(i, 0, vcount) { struct sdis_heat_vertex vtx; - res = sdis_heat_path_get_vertex(path, i, &vtx); - if (res != RES_OK) goto error; + ERR(sdis_heat_path_get_vertex(path, i, &vtx)); fprintf(stream, "%g %g %g\n", SPLIT3(vtx.P)); } /* Write the segment of the path */ @@ -138,8 +142,7 @@ dump_path fprintf(stream, "LOOKUP_TABLE vertex_type\n"); FOR_EACH(i, 0, vcount) { struct sdis_heat_vertex vtx; - res = sdis_heat_path_get_vertex(path, i, &vtx); - if (res != RES_OK) goto error; + ERR(sdis_heat_path_get_vertex(path, i, &vtx)); switch (vtx.type) { case SDIS_HEAT_VERTEX_CONDUCTION: fprintf(stream, "0.0\n"); break; case SDIS_HEAT_VERTEX_CONVECTION: fprintf(stream, "0.5\n"); break; @@ -156,8 +159,7 @@ dump_path fprintf(stream, "LOOKUP_TABLE default\n"); FOR_EACH(i, 0, vcount) { struct sdis_heat_vertex vtx; - res = sdis_heat_path_get_vertex(path, i, &vtx); - if (res != RES_OK) goto error; + ERR(sdis_heat_path_get_vertex(path, i, &vtx)); fprintf(stream, "%g\n", vtx.weight); } /* Write the time of the random walk vertices */ @@ -165,19 +167,17 @@ dump_path fprintf(stream, "LOOKUP_TABLE default\n"); FOR_EACH(i, 0, vcount) { struct sdis_heat_vertex vtx; - res = sdis_heat_path_get_vertex(path, i, &vtx); - if (res != RES_OK) goto error; + ERR(sdis_heat_path_get_vertex(path, i, &vtx)); fprintf(stream, "%g\n", IS_INF(vtx.time) ? FLT_MAX : vtx.time); } /* Write path type */ fprintf(stream, "CELL_DATA %u\n", 1u); fprintf(stream, "SCALARS Path_Type float 1\n"); fprintf(stream, "LOOKUP_TABLE path_type\n"); - res = sdis_heat_path_get_status(path, &status); - if (res != RES_OK) goto error; + ERR(sdis_heat_path_get_status(path, &status)); switch (status) { - case SDIS_HEAT_PATH_SUCCEED: fprintf(stream, "0.0\n"); break; - case SDIS_HEAT_PATH_FAILED: fprintf(stream, "1.0\n"); break; + case SDIS_HEAT_PATH_SUCCESS: fprintf(stream, "0.0\n"); break; + case SDIS_HEAT_PATH_FAILURE: fprintf(stream, "1.0\n"); break; default: FATAL("Unreachable code.\n"); break; } fprintf(stream, "LOOKUP_TABLE path_type 2\n"); @@ -207,10 +207,10 @@ print_sample unsigned id; size_t pcount, fcount; struct w_ctx* w_ctx = ctx; + const struct description* descs; CHK(path && ctx); - res = sdis_green_path_get_limit_point(path, &pt); - if (res != RES_OK) goto error; + ERR(sdis_green_path_get_limit_point(path, &pt)); /* For each path, prints: * # end #power_terms #flux_terms power_term_1 ... power_term_n flux_term_1 ... flux_term_n @@ -220,6 +220,7 @@ print_sample * - flux_term_i = flux_id_i factor_i */ + descs = darray_descriptions_cdata_get(w_ctx->desc); switch (pt.type) { case SDIS_FRAGMENT: { struct intface* d__; @@ -228,19 +229,19 @@ print_sample id = (pt.data.itfrag.fragment.side == SDIS_FRONT) ? d__->front_boundary_id : d__->back_boundary_id; ASSERT(id != UINT_MAX); - switch (w_ctx->desc[id].type) { + switch (descs[id].type) { case DESC_BOUND_T_FOR_SOLID: case DESC_BOUND_T_FOR_FLUID: - ASSERT(id == w_ctx->desc[id].d.t_boundary.mat_id); + ASSERT(id == descs[id].d.t_boundary.mat_id); printf("T\t%u", id); break; case DESC_BOUND_H_FOR_SOLID: case DESC_BOUND_H_FOR_FLUID: - ASSERT(id == w_ctx->desc[id].d.h_boundary.mat_id); + ASSERT(id == descs[id].d.h_boundary.mat_id); printf("H\t%u", id); break; case DESC_BOUND_F_FOR_SOLID: - ASSERT(id == w_ctx->desc[id].d.f_boundary.mat_id); + ASSERT(id == descs[id].d.f_boundary.mat_id); printf("X\t%u", id); break; default: FATAL("Unreachable code.\n"); break; @@ -250,15 +251,17 @@ print_sample case SDIS_VERTEX: type = sdis_medium_get_type(pt.data.mdmvert.medium); data = sdis_medium_get_data(pt.data.mdmvert.medium); - if (pt.data.mdmvert.vertex.P[0] == INF) { + if(pt.data.mdmvert.vertex.P[0] == INF) { /* Radiative output */ - printf("R\t%u", (unsigned)sa_size(w_ctx->desc)); + size_t sz = darray_descriptions_size_get(w_ctx->desc); + ASSERT(sz <= UINT_MAX); + printf("R\t%u", (unsigned)sz); } - else if (type == SDIS_FLUID) { + else if(type == SDIS_FLUID) { struct fluid* d__ = sdis_data_get(data); id = d__->id; - ASSERT(id == w_ctx->desc[id].d.fluid.fluid_id); - if (d__->is_outside) + ASSERT(id == descs[id].d.fluid.fluid_id); + if(d__->is_outside) /* If outside the model and in a fluid with known temperature, * its a fluid attached to a H boundary */ printf("H\t%u", id); @@ -268,7 +271,7 @@ print_sample struct solid* d__ = sdis_data_get(data); ASSERT(type == SDIS_SOLID); id = d__->id; - ASSERT(id == w_ctx->desc[id].d.solid.solid_id); + ASSERT(id == descs[id].d.solid.solid_id); ASSERT(!d__->is_outside); /* FIXME: what if in external solid? */ printf("S\t%u", id); } @@ -276,18 +279,15 @@ print_sample default: FATAL("Unreachable code.\n"); break; } - res = sdis_green_function_get_power_terms_count(path, &pcount); - if (res != RES_OK) goto error; + ERR(sdis_green_function_get_power_terms_count(path, &pcount)); htable_weigth_clear(&w_ctx->weigths); - res = sdis_green_path_for_each_flux_term(path, get_flux_terms, w_ctx); - if (res != RES_OK) goto error; + ERR(sdis_green_path_for_each_flux_term(path, get_flux_terms, w_ctx)); fcount = htable_weigth_size_get(&w_ctx->weigths); printf("\t%zu\t%zu", pcount, fcount); - res = sdis_green_path_for_each_power_term(path, print_power_term, ctx); - if (res != RES_OK) goto error; + ERR(sdis_green_path_for_each_power_term(path, print_power_term, ctx)); htable_weigth_begin(&w_ctx->weigths, &it); htable_weigth_end(&w_ctx->weigths, &end); @@ -315,11 +315,10 @@ dump_image size_t ix, iy; CHK(buf != NULL); - res = sdis_estimator_buffer_get_definition(buf, definition); - if (res != RES_OK) goto error; + ERR(sdis_estimator_buffer_get_definition(buf, definition)); temps = mem_alloc(definition[0] * definition[1] * sizeof(double)); - if (temps == NULL) { + if(temps == NULL) { res = RES_MEM_ERR; goto error; } @@ -336,10 +335,8 @@ dump_image FOR_EACH(ix, 0, definition[0]) { const struct sdis_estimator* estimator; struct sdis_mc T; - res = sdis_estimator_buffer_at(buf, ix, iy, &estimator); - if (res != RES_OK) goto error; - res = sdis_estimator_get_temperature(estimator, &T); - if (res != RES_OK) goto error; + ERR(sdis_estimator_buffer_at(buf, ix, iy, &estimator)); + ERR(sdis_estimator_get_temperature(estimator, &T)); fprintf(stdout, "%f\n", T.E); } } @@ -349,10 +346,8 @@ dump_image FOR_EACH(ix, 0, definition[0]) { const struct sdis_estimator* estimator; struct sdis_mc T; - res = sdis_estimator_buffer_at(buf, ix, iy, &estimator); - if (res != RES_OK) goto error; - res = sdis_estimator_get_temperature(estimator, &T); - if (res != RES_OK) goto error; + ERR(sdis_estimator_buffer_at(buf, ix, iy, &estimator)); + ERR(sdis_estimator_get_temperature(estimator, &T)); fprintf(stdout, "%f\n", T.SE); } } @@ -362,10 +357,8 @@ dump_image FOR_EACH(ix, 0, definition[0]) { const struct sdis_estimator* estimator; struct sdis_mc time; - res = sdis_estimator_buffer_at(buf, ix, iy, &estimator); - if (res != RES_OK) goto error; - res = sdis_estimator_get_realisation_time(estimator, &time); - if (res != RES_OK) goto error; + ERR(sdis_estimator_buffer_at(buf, ix, iy, &estimator)); + ERR(sdis_estimator_get_realisation_time(estimator, &time)); fprintf(stdout, "%f\n", time.E); } } @@ -375,10 +368,8 @@ dump_image FOR_EACH(ix, 0, definition[0]) { const struct sdis_estimator* estimator; struct sdis_mc time; - res = sdis_estimator_buffer_at(buf, ix, iy, &estimator); - if (res != RES_OK) goto error; - res = sdis_estimator_get_realisation_time(estimator, &time); - if (res != RES_OK) goto error; + ERR(sdis_estimator_buffer_at(buf, ix, iy, &estimator)); + ERR(sdis_estimator_get_realisation_time(estimator, &time)); fprintf(stdout, "%f\n", time.SE); } } @@ -388,10 +379,8 @@ dump_image FOR_EACH(ix, 0, definition[0]) { const struct sdis_estimator* estimator; size_t nfails; - res = sdis_estimator_buffer_at(buf, ix, iy, &estimator); - if (res != RES_OK) goto error; - res = sdis_estimator_get_failure_count(estimator, &nfails); - if (res != RES_OK) goto error; + ERR(sdis_estimator_buffer_at(buf, ix, iy, &estimator)); + ERR(sdis_estimator_get_failure_count(estimator, &nfails)); fprintf(stdout, "%zu\n", nfails); } } @@ -407,21 +396,24 @@ res_T dump_green (struct sdis_green_function* green, const struct counts* counts, - const struct description* descriptions, + const struct darray_descriptions* descriptions, const double radiative_temps[2]) { res_T res = RES_OK; - size_t ok_count, failed_count; + size_t ok_count, failed_count, sz; struct w_ctx w_ctx; - unsigned i; + int table_initialized = 0; + unsigned i, szd; + const struct description* descs; ASSERT(green && counts && descriptions && radiative_temps); - res = sdis_green_function_get_invalid_paths_count(green, &failed_count); - if (res != RES_OK) goto error; - - res = sdis_green_function_get_paths_count(green, &ok_count); - if (res != RES_OK) goto error; + ERR(sdis_green_function_get_invalid_paths_count(green, &failed_count)); + ERR(sdis_green_function_get_paths_count(green, &ok_count)); + sz = darray_descriptions_size_get(descriptions); + ASSERT(sz <= UINT_MAX); + szd = (unsigned)sz; + descs = darray_descriptions_cdata_get(descriptions); /* Output counts */ printf("---BEGIN GREEN---\n"); @@ -432,76 +424,73 @@ dump_green ok_count, failed_count); /* List Media */ - if (counts->smed_count) { + if(counts->smed_count) { printf("# Solids\n"); printf("# ID Name lambda rho cp power\n"); - FOR_EACH(i, 0, (unsigned)sa_size(descriptions)) { - const struct description* desc = descriptions + i; + FOR_EACH(i, 0, szd) { + const struct description* desc = descs + i; const struct mat_solid* sl; - if (desc->type != DESC_MAT_SOLID) continue; + if(desc->type != DESC_MAT_SOLID) continue; sl = &desc->d.solid; printf("%u\t%s\t%g\t%g\t%g\t%g\n", - i, sl->name, sl->lambda, sl->rho, sl->cp, sl->vpower); + i, str_cget(&sl->name), sl->lambda, sl->rho, sl->cp, sl->vpower); } } - if (counts->fmed_count) { + if(counts->fmed_count) { printf("# Fluids\n"); printf("# ID Name rho cp\n"); - FOR_EACH(i, 0, sa_size(descriptions)) { - const struct description* desc = descriptions + i; + FOR_EACH(i, 0, szd) { + const struct description* desc = descs + i; const struct mat_fluid* fl; fl = &desc->d.fluid; - printf("%u\t%s\t%g\t%g\n", i, fl->name, fl->rho, fl->cp); + printf("%u\t%s\t%g\t%g\n", i, str_cget(&fl->name), fl->rho, fl->cp); } } /* List Boundaries */ - if (counts->tbound_count) { + if(counts->tbound_count) { printf("# T Boundaries\n"); printf("# ID Name temperature\n"); - FOR_EACH(i, 0, sa_size(descriptions)) { - const struct description* desc = descriptions + i; + FOR_EACH(i, 0, szd) { + const struct description* desc = descs + i; const struct t_boundary* bd; - if (desc->type != DESC_BOUND_T_FOR_SOLID + if(desc->type != DESC_BOUND_T_FOR_SOLID && desc->type != DESC_BOUND_T_FOR_FLUID) continue; bd = &desc->d.t_boundary; - printf("%u\t%s\t%g\n", i, bd->name, bd->imposed_temperature); + printf("%u\t%s\t%g\n", i, str_cget(&bd->name), bd->imposed_temperature); } } - if (counts->hbound_count) { + if(counts->hbound_count) { printf("# H Boundaries\n"); printf("# ID Name emissivity specular_fraction hc hc_max T_env\n"); - FOR_EACH(i, 0, sa_size(descriptions)) { - const struct description* desc = descriptions + i; + FOR_EACH(i, 0, szd) { + const struct description* desc = descs + i; const struct h_boundary* bd; - if (desc->type != DESC_BOUND_H_FOR_SOLID + if(desc->type != DESC_BOUND_H_FOR_SOLID && desc->type != DESC_BOUND_H_FOR_FLUID) continue; bd = &desc->d.h_boundary; - printf("%u\t%s\t%g\t%g\t%g\t%g\t%g\n", - i, bd->name, (bd->has_emissivity ? bd->emissivity : 0), - (bd->has_emissivity ? bd->specular_fraction : 0), - (bd->has_hc ? bd->hc : 0), (bd->has_hc ? bd->hc_max : 0), - bd->imposed_temperature); + printf("%u\t%s\t%g\t%g\t%g\t%g\n", + i, str_cget(&bd->name), bd->emissivity, bd->specular_fraction, + bd->hc, bd->imposed_temperature); } } - if (counts->fbound_count) { + if(counts->fbound_count) { printf("# F Boundaries\n"); printf("# ID Name flux\n"); - FOR_EACH(i, 0, sa_size(descriptions)) { - const struct description* desc = descriptions + i; + FOR_EACH(i, 0, szd) { + const struct description* desc = descs + i; const struct f_boundary* bd; - if (desc->type != DESC_BOUND_F_FOR_SOLID) continue; + if(desc->type != DESC_BOUND_F_FOR_SOLID) continue; bd = &desc->d.f_boundary; printf("%u\t%s\t%g\n", - i, bd->name, bd->imposed_flux); + i, str_cget(&bd->name), bd->imposed_flux); } } /* Radiative Temperatures */ printf("# Radiative Temperatures\n"); printf("# ID Rad_Temp Lin_Temp\n"); - printf("%u\t%g\t%g\n", (unsigned)sa_size(descriptions), - radiative_temps[0], radiative_temps[1]); + printf("%u\t%g\t%g\n", szd, radiative_temps[0], radiative_temps[1]); printf("# Samples\n"); printf("# end #power_terms #flux_terms power_term_1 ... power_term_n flux_term_1 ... flux_term_n\n"); @@ -511,146 +500,70 @@ dump_green w_ctx.desc = descriptions; htable_weigth_init(NULL, &w_ctx.weigths); + table_initialized = 1; - res = sdis_green_function_for_each_path(green, print_sample, &w_ctx); - htable_weigth_release(&w_ctx.weigths); - if (res != RES_OK) goto error; - + ERR(sdis_green_function_for_each_path(green, print_sample, &w_ctx)); printf("---END GREEN---\n"); - CHK(sdis_green_function_ref_put(green) == RES_OK); + ERR(sdis_green_function_ref_put(green)); end: + if(table_initialized) htable_weigth_release(&w_ctx.weigths); return res; error: goto end; } +int +has_holes(const struct stardis* stardis) +{ + struct senc3d_scene* senc3d_scn; + unsigned count; + SDIS(scene_get_senc3d_scene(stardis->geometry.sdis_scn, &senc3d_scn)); + SENC3D(scene_get_frontier_segments_count(senc3d_scn, &count)); + return count != 0; +} + res_T -create_edge_file_if - (struct sdis_scene* scn, - struct mem_allocator* allocator) +dump_holes_at_the_end_of_vtk + (const struct stardis* stardis, + FILE* stream) { res_T res = RES_OK; - unsigned i, comp, scount, vcount; - char name[64]; - FILE* obj = NULL; - struct htable_vrtx_rank ranks; - char structs_initialized = 0; + unsigned* trgs = NULL; struct senc3d_scene* senc3d_scn = NULL; - struct darray_uint comps; - struct darray_uint edges; - unsigned cc_count = 0; - unsigned *cc, *ee, uimax = UINT_MAX; - - ASSERT(scn && allocator); - - res = sdis_scene_get_senc3d_scene(scn, &senc3d_scn); - if (res != RES_OK) goto error; - res = senc3d_scene_get_frontier_segments_count(senc3d_scn, &scount); - if (res != RES_OK) goto error; - - if (scount == 0) - /* No edge to exit */ - goto end; - - fprintf(stderr, "Some frontier edges found: create OBJ files\n"); + unsigned tsz, scount, i; + ASSERT(stardis && stream); + ASSERT(stardis); - htable_vrtx_rank_init(allocator, &ranks); - darray_uint_init(allocator, &comps); - darray_uint_init(allocator, &edges); - structs_initialized = 1; + ERR(sg3d_geometry_get_unique_triangles_count(stardis->geometry.sg3d, &tsz)); + /* Store frontier triangles */ + trgs = MEM_CALLOC(stardis->allocator, tsz, sizeof(*trgs)); + if(!trgs) { + res = RES_MEM_ERR; + goto error; + } - /* Exit vertices, with a new numbering scheme */ - vcount = 0; + /* Keep the involved segments (not the vertices) */ + ERR(sdis_scene_get_senc3d_scene(stardis->geometry.sdis_scn, &senc3d_scn)); + ERR(senc3d_scene_get_frontier_segments_count(senc3d_scn, &scount)); FOR_EACH(i, 0, scount) { - unsigned v, edge[2], rk[2]; - res = senc3d_scene_get_frontier_segment(senc3d_scn, i, edge); - if (res != RES_OK) goto error; - FOR_EACH(v, 0, 2) { - unsigned* p_rank; - p_rank = htable_vrtx_rank_find(&ranks, edge + v); - if (p_rank) { - rk[v] = *p_rank; - } else { - /* First appearance: allocate a rank and exit it */ - unsigned rank; - size_t tmp = htable_vrtx_rank_size_get(&ranks); - ASSERT(tmp < UINT_MAX); - rank = (unsigned)tmp; - res = htable_vrtx_rank_set(&ranks, edge + v, &rank); - if (res != RES_OK) goto error; - rk[v] = rank; - darray_uint_push_back(&comps, &uimax); - darray_uint_push_back(&edges, edge + v); - vcount++; - } - /* Manage components */ - cc = darray_uint_data_get(&comps); - if (cc[rk[0]] != UINT_MAX || cc[rk[1]] != UINT_MAX) { - /* If both components known, they must be the same */ - ASSERT(cc[rk[0]] == UINT_MAX || cc[rk[1]] == UINT_MAX - || cc[rk[0]] == cc[rk[1]]); - FOR_EACH(v, 0, 2) { - unsigned w = 1 - v; /* The other vertex */ - if (cc[rk[v]] == UINT_MAX) cc[rk[v]] = cc[rk[w]]; - } - } else { - /* None are in a know component: create a new one */ - cc[rk[0]] = cc[rk[1]] = cc_count++; - } - } + unsigned vrtc[2], trid; + ERR(senc3d_scene_get_frontier_segment(senc3d_scn, i, vrtc, &trid)); + trgs[trid] = 1; } - /* Exit segments by component using the new numbering scheme */ - cc = darray_uint_data_get(&comps); - ee = darray_uint_data_get(&edges); - FOR_EACH(comp, 0, cc_count) { - unsigned v; - snprintf(name, sizeof(name), "frontier_component_%u.obj", comp); - obj = fopen(name, "w"); - if (!obj) goto error; - /* Exit all vertices even unused - * (same numbering scheme for all components) */ - FOR_EACH(v, 0, vcount) { - double coord[3]; - res = senc3d_scene_get_vertex(senc3d_scn, ee[v], coord); - if (res != RES_OK) goto error; - fprintf(obj, "v %f %f %f\n", SPLIT3(coord)); - } - FOR_EACH(i, 0, scount) { - unsigned edge[2], new_ranks[2]; - res = senc3d_scene_get_frontier_segment(senc3d_scn, i, edge); - if (res != RES_OK) goto error; - ASSERT(cc[edge[0]] == cc[edge[1]]); - if (cc[edge[0]] != comp) - /* Segment not in this component */ - continue; - FOR_EACH(v, 0, 2) { - unsigned* p_rank; - p_rank = htable_vrtx_rank_find(&ranks, edge + v); - ASSERT(p_rank && *p_rank < vcount); - new_ranks[v] = *p_rank; - } -#define OBJ_IDX(Rank) ((Rank)+1) - fprintf(obj, "l %u %u\n", OBJ_IDX(new_ranks[0]), OBJ_IDX(new_ranks[1])); -#undef OBJ_IDX - } - fclose(obj); obj = NULL; - } + fprintf(stream, "SCALARS Holes_frontiers int\n"); + fprintf(stream, "LOOKUP_TABLE default\n"); + FOR_EACH(i, 0, tsz) fprintf(stream, "%d\n", trgs[i]); -end: - if (obj) fclose(obj); - if (senc3d_scn) senc3d_scene_ref_put(senc3d_scn); - if (structs_initialized) { - htable_vrtx_rank_release(&ranks); - darray_uint_release(&comps); - darray_uint_release(&edges); - } +exit: + if(trgs) MEM_RM(stardis->allocator, trgs); + if(senc3d_scn) senc3d_scene_ref_put(senc3d_scn); return res; error: - goto end; + goto exit; } res_T @@ -666,10 +579,8 @@ print_single_MC_result ASSERT(estimator && stardis); /* Fetch the estimation data */ - res = sdis_estimator_get_temperature(estimator, &result); - if (res != RES_OK) goto error; - res = sdis_estimator_get_failure_count(estimator, &nfailures); - if (res != RES_OK) goto error; + ERR(sdis_estimator_get_temperature(estimator, &result)); + ERR(sdis_estimator_get_failure_count(estimator, &nfailures)); /* Print the results */ switch (mode & COMPUTE_MODES) { @@ -682,55 +593,50 @@ print_single_MC_result break; case MEDIUM_COMPUTE: printf("Temperature in medium %s at t=%g = %g +/- %g\n", - stardis->solve_name, stardis->probe[3], + str_cget(&stardis->solve_name), stardis->probe[3], result.E, /* Expected value */ result.SE); /* Standard error */ break; case BOUNDARY_COMPUTE: printf("Temperature at boundary %s at t=%g = %g +/- %g\n", - stardis->solve_name, stardis->probe[3], + str_cget(&stardis->solve_name), stardis->probe[3], result.E, /* Expected value */ result.SE); /* Standard error */ break; case FLUX_BOUNDARY_COMPUTE: { enum sdis_estimator_type type; - sdis_estimator_get_type(estimator, &type); - if (res != RES_OK) goto error; + ERR(sdis_estimator_get_type(estimator, &type)); ASSERT(type == SDIS_ESTIMATOR_FLUX); printf("Temperature at boundary %s at t=%g = %g +/- %g\n", - stardis->solve_name, stardis->probe[3], + str_cget(&stardis->solve_name), stardis->probe[3], result.E, /* Expected value */ result.SE); /* Standard error */ - res = sdis_estimator_get_convective_flux(estimator, &result); - if (res != RES_OK) goto error; + ERR(sdis_estimator_get_convective_flux(estimator, &result)); printf("Convective flux at boundary %s at t=%g = %g +/- %g\n", - stardis->solve_name, stardis->probe[3], + str_cget(&stardis->solve_name), stardis->probe[3], result.E, /* Expected value */ result.SE); /* Standard error */ - res = sdis_estimator_get_radiative_flux(estimator, &result); - if (res != RES_OK) goto error; + ERR(sdis_estimator_get_radiative_flux(estimator, &result)); printf("Radiative flux at boundary %s at t=%g = %g +/- %g\n", - stardis->solve_name, stardis->probe[3], + str_cget(&stardis->solve_name), stardis->probe[3], result.E, /* Expected value */ result.SE); /* Standard error */ - res = sdis_estimator_get_total_flux(estimator, &result); - if (res != RES_OK) goto error; + ERR(sdis_estimator_get_total_flux(estimator, &result)); printf("Total flux Flux at boundary %s at t=%g = %g +/- %g\n", - stardis->solve_name, stardis->probe[3], + str_cget(&stardis->solve_name), stardis->probe[3], result.E, /* Expected value */ result.SE); /* Standard error */ break; } default: FATAL("Invalid mode."); } - printf("#failures: %zu/%zu\n", nfailures, stardis->N); - if (nfailures) - fprintf(stderr, "#failures: %zu/%zu\n", nfailures, stardis->N); + printf("#failures: %zu/%zu\n", nfailures, stardis->samples); + if(nfailures) + fprintf(stderr, "#failures: %zu/%zu\n", nfailures, stardis->samples); /* Dump paths according to user settings */ - res = sdis_estimator_for_each_path(estimator, dump_path, stdout); - if (res != RES_OK) goto error; + ERR(sdis_estimator_for_each_path(estimator, dump_path, stdout)); end: return res; error: @@ -738,127 +644,141 @@ error: } void -dump_vtk - (const struct geometry* geometry) -{ - size_t i; - struct vertex* vtx; - struct triangle* tri; - - ASSERT(geometry); - - vtx = geometry->vertex; - tri = geometry->triangle; - - printf("# vtk DataFile Version 2.0\nvtk output\nASCII\nDATASET POLYDATA\n"); - printf("POINTS %zu float\n\n", sa_size(vtx)); - for (i = 0; i < sa_size(vtx); ++i) { - printf("%f %f %f\n", SPLIT3(vtx[i].xyz)); - } - printf("\nPOLYGONS %zu %zu\n", sa_size(tri), 4 * sa_size(tri)); - for (i = 0; i < sa_size(tri); ++i) { - printf("3 %u %u %u\n", SPLIT3(tri[i].indices.data)); - } - printf("\nCELL_DATA %zu\n", sa_size(tri)); - printf("SCALARS front_description long_long_int 1\n"); - printf("LOOKUP_TABLE default\n"); - for (i = 0; i < sa_size(tri); ++i) { - printf("%lld\n", (long long int)tri[i].front_description); - } - printf("SCALARS back_description long_long_int 1\n"); - printf("LOOKUP_TABLE default\n"); - for (i = 0; i < sa_size(tri); ++i) { - printf("%lld\n", (long long int)tri[i].back_description); - } -} - -void dump_map (const struct stardis* stardis, - struct sdis_estimator** estimators) + const struct darray_estimators* estimators) { - size_t i, last_v = 0; - size_t* idx; - struct vertex* vtx; - struct triangle* tri; - - ASSERT(stardis); - - idx = stardis->boundary.primitives; - vtx = stardis->geometry.vertex; - tri = stardis->geometry.triangle; + unsigned i, vcount, tcount, last_v = 0; + const size_t* idx; + size_t sz; + unsigned szp; + const struct sdis_estimator* const* est; + + ASSERT(stardis && estimators); + + est = darray_estimators_cdata_get(estimators); + idx = darray_size_t_cdata_get(&stardis->compute_region.primitives); + sz = darray_size_t_size_get(&stardis->compute_region.primitives); + ASSERT(sz <= UINT_MAX); + szp = (unsigned)sz; + SG3D(geometry_get_unique_vertices_count(stardis->geometry.sg3d, &vcount)); + SG3D(geometry_get_unique_triangles_count(stardis->geometry.sg3d, &tcount)); /* Find last used vertex */ - for (i = 0; i < sa_size(idx); ++i) { - size_t t = idx[i]; - last_v = MMAX(MMAX(last_v, tri[t].indices.data[0]), - MMAX(tri[t].indices.data[1], tri[t].indices.data[2])); + for(i = 0; i < szp; ++i) { + unsigned t; + unsigned indices[3]; + ASSERT(idx[i] <= UINT_MAX); + t = (unsigned)idx[i]; + SG3D(geometry_get_unique_triangle_vertices(stardis->geometry.sg3d, t, + indices)); + last_v = MMAX(MMAX(last_v, indices[0]), MMAX(indices[1], indices[2])); } /* Dump vertices up to last_v, even unused ones, to avoid reindexing */ printf("# vtk DataFile Version 2.0\nvtk output\nASCII\nDATASET POLYDATA\n"); - printf("POINTS %zu float\n\n", last_v + 1); - for (i = 0; i <= last_v; ++i) { - printf("%f %f %f\n", SPLIT3(vtx[i].xyz)); + printf("POINTS %u float\n\n", last_v + 1); + for(i = 0; i <= last_v; ++i) { + double coord[3]; + SG3D(geometry_get_unique_vertex(stardis->geometry.sg3d, i, coord)); + printf("%f %f %f\n", SPLIT3(coord)); } /* Dump only primitives in boundary */ - printf("\nPOLYGONS %zu %zu\n", sa_size(idx), 4 * sa_size(idx)); - for (i = 0; i < sa_size(idx); ++i) { - size_t t = idx[i]; - printf("3 %u %u %u\n", SPLIT3(tri[t].indices.data)); - } - printf("\nCELL_DATA %zu\n", sa_size(idx)); + printf("\nPOLYGONS %u %u\n", szp, 4 * szp); + for(i = 0; i < szp; ++i) { + unsigned t; + unsigned indices[3]; + ASSERT(idx[i] <= UINT_MAX); + t = (unsigned)idx[i]; + SG3D(geometry_get_unique_triangle_vertices(stardis->geometry.sg3d, t, + indices)); + printf("3 %u %u %u\n", SPLIT3(indices)); + } + printf("\nCELL_DATA %u\n", szp); printf("SCALARS temperature_estimate float 1\n"); printf("LOOKUP_TABLE default\n"); - for (i = 0; i < sa_size(idx); ++i) { - const struct sdis_estimator* estimator = estimators[i]; + for(i = 0; i < szp; ++i) { struct sdis_mc T; - SDIS(estimator_get_temperature(estimator, &T)); + SDIS(estimator_get_temperature(est[i], &T)); printf("%f\n", T.E); } printf("SCALARS temperature_std_dev float 1\n"); printf("LOOKUP_TABLE default\n"); - for (i = 0; i < sa_size(idx); ++i) { - const struct sdis_estimator* estimator = estimators[i]; + for(i = 0; i < szp; ++i) { struct sdis_mc T; - SDIS(estimator_get_temperature(estimator, &T)); + SDIS(estimator_get_temperature(est[i], &T)); printf("%f\n", T.SE); } printf("SCALARS temperature_failures_count unsigned_long_long 1\n"); printf("LOOKUP_TABLE default\n"); - for (i = 0; i < sa_size(idx); ++i) { - const struct sdis_estimator* estimator = estimators[i]; + for(i = 0; i < szp; ++i) { size_t nfails; - SDIS(estimator_get_failure_count(estimator, &nfails)); - printf("%zu\n", nfails); + SDIS(estimator_get_failure_count(est[i], &nfails)); + ASSERT(nfails <= UINT_MAX); + printf("%u\n", (unsigned)nfails); } printf("SCALARS computation_time_estimate float 1\n"); printf("LOOKUP_TABLE default\n"); - for (i = 0; i < sa_size(idx); ++i) { - const struct sdis_estimator* estimator = estimators[i]; + for(i = 0; i < szp; ++i) { struct sdis_mc time; - SDIS(estimator_get_realisation_time(estimator, &time)); + SDIS(estimator_get_realisation_time(est[i], &time)); printf("%f\n", time.E); } printf("SCALARS computation_time_std_dev float 1\n"); printf("LOOKUP_TABLE default\n"); - for (i = 0; i < sa_size(idx); ++i) { - const struct sdis_estimator* estimator = estimators[i]; + for(i = 0; i < szp; ++i) { struct sdis_mc time; - SDIS(estimator_get_realisation_time(estimator, &time)); + SDIS(estimator_get_realisation_time(est[i], &time)); printf("%f\n", time.SE); } } -void -print_trg_as_obj - (FILE* stream, - const struct vertex* vertices, - const unsigned* indices) +res_T +dump_compute_region_at_the_end_of_vtk + (const struct stardis* stardis, + FILE* stream) { - ASSERT(stream && vertices && indices); - fprintf(stream, "v %.8f %.8f %.8f\nv %.8f %.8f %.8f\nv %.8f %.8f %.8f\nf 1 2 3\n", - SPLIT3(vertices[indices[0]].xyz), - SPLIT3(vertices[indices[1]].xyz), - SPLIT3(vertices[indices[2]].xyz)); -} + res_T res = RES_OK; + unsigned char* v = NULL; + unsigned tsz, i; + size_t j; + ASSERT(stardis && stream); + ASSERT(darray_size_t_size_get(&stardis->compute_region.primitives) + == darray_sides_size_get(&stardis->compute_region.sides)); + + ERR(sg3d_geometry_get_unique_triangles_count(stardis->geometry.sg3d, &tsz)); + /* For triangles not in compute region v==0 */ + v = MEM_CALLOC(stardis->allocator, tsz, sizeof(*v)); + if(!v) { + res = RES_MEM_ERR; + goto error; + } + + /* For triangles in compute region, v==1 if FRONT or v==2 for BACK */ + FOR_EACH(j, 0, darray_size_t_size_get(&stardis->compute_region.primitives)) { + size_t prim = darray_size_t_cdata_get(&stardis->compute_region.primitives)[j]; + enum sdis_side side = darray_sides_cdata_get(&stardis->compute_region.sides)[j]; + ASSERT(prim <= tsz); + ASSERT(v[(unsigned)prim] == 0); /* Not the same triangle twice */ + v[(unsigned)prim] = (side == SDIS_FRONT) ? 1 : 2; + } + + /* For triangles in compute region with error v==MAX */ + FOR_EACH(j, 0, darray_uint_size_get(&stardis->compute_region.err_triangles)) { + unsigned prim = darray_uint_cdata_get(&stardis->compute_region.err_triangles)[j]; + ASSERT(prim <= tsz); + ASSERT(v[prim] == 0); /* Not the same triangle twice */ + v[(unsigned)prim] = UCHAR_MAX; + } + + fprintf(stream, "SCALARS Compute_region int\n"); + fprintf(stream, "LOOKUP_TABLE default\n"); + FOR_EACH(i, 0, tsz) + fprintf(stream, "%d\n", v[i] == UCHAR_MAX ? INT_MAX : v[i]); + +exit: + if(v) MEM_RM(stardis->allocator, v); + return res; +error: + goto exit; +} +\ No newline at end of file diff --git a/src/stardis-output.h b/src/stardis-output.h @@ -12,12 +12,13 @@ struct sdis_estimator_buffer; struct sdis_green_function; struct counts; struct description; -struct sdis_scene; struct mem_allocator; struct sdis_estimator; struct stardis; struct geometry; struct vertex; +struct darray_descriptions; +struct darray_estimators; extern res_T print_sample @@ -32,13 +33,16 @@ extern res_T dump_green (struct sdis_green_function* green, const struct counts* counts, - const struct description* descriptions, + const struct darray_descriptions* descriptions, const double radiative_temps[2]); +extern int +has_holes(const struct stardis* stardis); + extern res_T -create_edge_file_if - (struct sdis_scene* scn, - struct mem_allocator* allocator); +dump_holes_at_the_end_of_vtk + (const struct stardis* stardis, + FILE* stream); extern res_T print_single_MC_result @@ -47,18 +51,13 @@ print_single_MC_result struct stardis* stardis); extern void -dump_vtk - (const struct geometry* geometry); - -extern void dump_map (const struct stardis* stardis, - struct sdis_estimator** estimators); + const struct darray_estimators* estimators); -extern void -print_trg_as_obj - (FILE* stream, - const struct vertex* vertices, - const unsigned* indices); +extern res_T +dump_compute_region_at_the_end_of_vtk + (const struct stardis* stardis, + FILE* stream); #endif diff --git a/src/stardis-parsing.c b/src/stardis-parsing.c @@ -4,6 +4,7 @@ #include "stardis-version.h" #include <rsys/cstr.h> +#include <rsys/double3.h> #include <sdis_version.h> #include <getopt.h> @@ -22,7 +23,7 @@ _strupr (char* s) { char* ptr; - for (ptr = s; *ptr; ++ptr) { + for(ptr = s; *ptr; ++ptr) { int tmp = toupper(*ptr); ASSERT(tmp == (char)tmp); *ptr = (char)tmp; @@ -49,7 +50,7 @@ split_line /* Count how many elements will be extracted. */ while (*tmp) { - if (a_delim == *tmp) + if(a_delim == *tmp) { count++; last_comma = tmp; @@ -66,7 +67,7 @@ split_line result = malloc(sizeof(char*) * count); - if (result) + if(result) { size_t idx = 0; char* token = strtok(a_str, delim); @@ -92,28 +93,215 @@ static void print_version () { - printf("stardis-app version %i.%i.%i built on stardis library version %i.%i.%i\n", + printf("Stardis version %i.%i.%i built on stardis solver version %i.%i.%i\n", StardisApp_VERSION_MAJOR, StardisApp_VERSION_MINOR, StardisApp_VERSION_PATCH, Stardis_VERSION_MAJOR, Stardis_VERSION_MINOR, Stardis_VERSION_PATCH); } +void +add_geom_ctx_indices + (const unsigned itri, + unsigned ids[3], + void* context) +{ + const struct add_geom_ctx* ctx = context; + const unsigned* trg; + int i; + ASSERT(ids && ctx); + ASSERT(itri < ctx->stl_desc.triangles_count); + trg = ctx->stl_desc.indices + 3 * itri; + for(i = 0; i < 3; i++) ids[i] = trg[i]; +} + +void +add_geom_ctx_properties + (const unsigned itri, + unsigned prop[3], + void* context) +{ + const struct add_geom_ctx* ctx = context; + int i; + ASSERT(prop && ctx); (void)itri; + ASSERT(itri < ctx->stl_desc.triangles_count); + /* Same media for the whole add_geometry set of triangles */ + for(i = 0; i < SG3D_PROP_TYPES_COUNT__; i++) prop[i] = ctx->properties[i]; +} + +void +add_geom_ctx_position + (const unsigned ivert, + double pos[3], + void* context) +{ + const struct add_geom_ctx* ctx = context; + const float* v; + ASSERT(pos && ctx); + ASSERT(ivert < ctx->stl_desc.vertices_count); + v = ctx->stl_desc.vertices + 3 * ivert; + d3_set_f3(pos, v); +} + +static res_T +read_sides_and_files + (struct stardis* stardis, + const int side_is_intface, /* if 1, don't read side */ + const unsigned id, + char** tok_ctx) +{ + char* tk = NULL; + int file_count = 0; + struct sstl* sstl = NULL; + struct add_geom_ctx add_geom_ctx; + struct sg3d_geometry_add_callbacks callbacks = SG3D_ADD_CALLBACKS_NULL__; + res_T res = RES_OK; + + ASSERT(stardis && tok_ctx); + + callbacks.get_indices = add_geom_ctx_indices; + callbacks.get_properties = add_geom_ctx_properties; + callbacks.get_position = add_geom_ctx_position; + + /* At least one side+name, no side without name */ + ERR(sstl_create(NULL, stardis->allocator, 0, &sstl)); + for(;;) { + struct str name; + if(side_is_intface) { + add_geom_ctx.properties[SG3D_FRONT] = SG3D_UNSPECIFIED_PROPERTY; + add_geom_ctx.properties[SG3D_BACK] = SG3D_UNSPECIFIED_PROPERTY; + add_geom_ctx.properties[SG3D_INTFACE] = id; + } else { + tk = strtok_r(NULL, " ", tok_ctx); + if(!tk) { + if(file_count == 0) { + /* At least 1 side */ + fprintf(stderr, "Invalid data (missing token 'side')\n"); + res = RES_BAD_ARG; + goto error; + } + else break; + } + _strupr(tk); + add_geom_ctx.properties[SG3D_INTFACE] = SG3D_UNSPECIFIED_PROPERTY; + if(0 == strcmp(tk, "FRONT")) { + add_geom_ctx.properties[SG3D_FRONT] = id; + add_geom_ctx.properties[SG3D_BACK] = SG3D_UNSPECIFIED_PROPERTY; + } + else if(0 == strcmp(tk, "BACK")) { + add_geom_ctx.properties[SG3D_FRONT] = SG3D_UNSPECIFIED_PROPERTY; + add_geom_ctx.properties[SG3D_BACK] = id; + } + else if(0 == strcmp(tk, "BOTH")) { + add_geom_ctx.properties[SG3D_FRONT] = id; + add_geom_ctx.properties[SG3D_BACK] = id; + } + else { + fprintf(stderr, "Invalid side specifier: %s\n", tk); + res = RES_BAD_ARG; + goto error; + } + } + tk = strtok_r(NULL, " ", tok_ctx); + if(!tk) { + if(!side_is_intface /* Has read a side */ + || !file_count) /* Need at least 1 file */ + { + fprintf(stderr, "Invalid data (missing token 'file name')\n"); + res = RES_BAD_ARG; + goto error; + } + else break; + } + file_count++; + res = sstl_load(sstl, tk); + if(res != RES_OK) { + fprintf(stderr, "Cannot read STL file: %s\n", tk); + goto error; + } + ERR(sstl_get_desc(sstl, &add_geom_ctx.stl_desc)); + ASSERT(add_geom_ctx.stl_desc.vertices_count <= UINT_MAX + && add_geom_ctx.stl_desc.triangles_count <= UINT_MAX); + /* Keep file names for error reports */ + str_init(stardis->allocator, &name); + ERR(str_set(&name, tk)); + str_release(&name); + + res = sg3d_geometry_add( + stardis->geometry.sg3d, + (unsigned)add_geom_ctx.stl_desc.vertices_count, + (unsigned)add_geom_ctx.stl_desc.triangles_count, + &callbacks, + &add_geom_ctx); + + if(res != RES_OK) { + fprintf(stderr, "Cannot add file content: %s\n", tk); + goto error; + } + } + +end: + if(sstl) SSTL(ref_put(sstl)); + return res; +error: + goto end; +} + /******************************************************************************* * Public Functions ******************************************************************************/ +res_T +init_args(struct mem_allocator* alloc, struct args** out_args) +{ + res_T res = RES_OK; + struct args* args = NULL; + ASSERT(alloc && out_args); + + args = calloc(sizeof(struct args), 1); + if(!args) { + res = RES_MEM_ERR; + goto error; + } + + args->allocator = alloc; + darray_str_init(alloc, &args->model_files); + /* Some fields have non-zero default values */ + args->samples = 10000; + args->nthreads = DEFAULT_NTHREADS; + args->probe[3] = INF; /* Probe time */ + args->scale_factor = 1; + args->radiative_temp[0] = args->radiative_temp[1] = 300; + +end: + *out_args = args; + return res; +error: + if(args) release_args(args); + args = NULL; + goto end; +} + +void +release_args(struct args* args) +{ + ASSERT(args && args->allocator); + args->allocator = NULL; + darray_str_release(&args->model_files); +} + char mode_option (const enum stardis_mode m) { int found = 0; char res = '?'; - if (m & PROBE_COMPUTE) { found++; res = 'p'; } - if (m & PROBE_COMPUTE_ON_INTERFACE) { found++; res = 'P'; } - if (m & MEDIUM_COMPUTE) { found++; res = 'm'; } - if (m & BOUNDARY_COMPUTE) { found++; res = 's'; } - if (m & IR_COMPUTE) { found++; res = 'R'; } - if (m & MAP_COMPUTE) { found++; res = 'S'; } + if(m & PROBE_COMPUTE) { found++; res = 'p'; } + if(m & PROBE_COMPUTE_ON_INTERFACE) { found++; res = 'P'; } + if(m & MEDIUM_COMPUTE) { found++; res = 'm'; } + if(m & BOUNDARY_COMPUTE) { found++; res = 's'; } + if(m & IR_COMPUTE) { found++; res = 'R'; } + if(m & MAP_COMPUTE) { found++; res = 'S'; } if (m & DUMP_VTK) { found++; res = 'd'; } + if (m & FLUX_BOUNDARY_COMPUTE) { found++; res = 'F'; } ASSERT(found == 1); return res; } @@ -127,7 +315,7 @@ void print_multiple_modes ASSERT(stream); do { m = BIT(b++); - if (m & modes) { + if(m & modes) { fprintf(stream, (fst ? " -%c" : ", -%c"), mode_option(m)); fst = 0; } @@ -136,82 +324,86 @@ void print_multiple_modes void print_help - (char* prog) + (const char* prog) { + char* name; +#ifdef COMPILER_GCC + name = 1 + strrchr(prog, '/'); +#else + name = 1 + strrchr(prog, '\\'); +#endif ASSERT(prog); print_version(); - printf("\nUsage: \n%s -M MEDIUM.txt -B BOUNDARY.txt\n", prog); - printf("[ -p X,Y,Z,TIME | -m MEDIUM_NAME,TIME | -d \n"); - printf(" | -s SURFACE.vtk,TIME | -S SURFACE.vtk,TIME\n"); - printf(" | -R t=TIME:spp=SPP:fov=FOV:up=XUP,YUP,ZUP:pos=X,Y,Z:tgt=XT,YT,ZT:img=WxH ]\n"); - printf("[ -g] [-f SCALE_FACTOR] [-D {all | error | success}] -d\n"); - printf("[-n NUM_OF_REALIZATIONS] [-t NUM_OF_THREADS]\n"); - printf("[-r AMBIENT_RAD_TEMP:REFERENCE_TEMP]\n"); - printf("\n -h: print this help.\n"); - printf(" -v: print stardis-app version.\n"); - printf("\nMandatory arguments\n"); - printf("-------------------\n"); - printf("\n -M MEDIUM.txt:\n"); - printf(" MEDIUM.txt is a text file which contains the description of the\n"); - printf(" media.\n"); - printf("\n -B BOUNDARY.txt:\n"); - printf(" BOUNDARY.txt is a text file which contains the description of the\n"); - printf(" boundary conditions.\n"); - printf("\nExclusive arguments\n"); - printf("-------------------\n"); + printf("\nUsage: \n"); + printf("------\n"); + printf("\n%s {-M MODEL_FILE_NAME}+\n", name); + printf(" [ -p X,Y,Z,TIME | -P X,Y,Z,TIME | -m MEDIUM_NAME,TIME\n"); + printf(" | -s STL_FILE_NAME,TIME | -S STL_FILE_NAME,TIME | -F STL_FILE_NAME,TIME\n"); + printf(" | -R t=TIME:spp=SPP:fov=FOV:up=XUP,YUP,ZUP:pos=X,Y,Z:tgt=XT,YT,ZT:img=WxH ]\n"); + printf(" [-d] [-g] [-f SCALE_FACTOR] [-D {all|error|success}] [-n REALIZATIONS_COUNT]\n"); + printf(" [-t NUM_OF_THREADS] [-r AMBIENT_RAD_TEMP:REFERENCE_TEMP]\n"); + printf("\n%s -h: print this help.\n", name); + printf("\n%s -v: print stardis-app version.\n", name); + printf("\nMandatory arguments:\n"); + printf("--------------------\n"); + printf("\n -M <MODEL_FILE_NAME>:\n"); + printf(" Add a text file that contains (partial) description of the model.\n"); + printf(" Can include both media and boundary conditions, in any order.\n"); + printf(" Can be used more than once if description is split across different files.\n"); + printf("\nExclusive arguments:\n"); + printf("--------------------\n"); printf("\n -p X,Y,Z:TIME:\n"); printf(" Compute the temperature at the given probe.\n"); printf(" The probe must be in a medium.\n"); - printf(" If some boundary conditions are time-dependant, TIME cannot be INF.\n"); printf("\n -P X,Y,Z:TIME:\n"); printf(" Compute the temperature at the given probe on an interface.\n"); - printf(" The probe must be in a medium and is moved to the closest interface.\n"); - printf(" If some boundary conditions are time-dependant, TIME cannot be INF.\n"); + printf(" The probe is suposed to be on an interface and is adjusted to the\n"); + printf(" closest point of the closest interface.\n"); printf("\n -m MEDIUM_NAME,TIME:\n"); - printf(" Compute the mean temperature in a given medium at a given time.\n"); - printf(" If some boundary conditions are time-dependant, TIME cannot be INF.\n"); - printf("\n -s SURFACE.stl,TIME:\n"); - printf(" Compute the mean temperature on a given 2D shape at a given time.\n"); - printf(" Mean temperature is computed on the front sides of the primitives\n"); - printf(" listed in SURFACE.stl. These primitives are not added to the geometry,\n"); - printf(" but must be already known. The shape doesn't need to be connex.\n"); - printf(" If some boundary conditions are time-dependant, TIME cannot be INF.\n"); - printf("\n -S SURFACE.stl,TIME:\n"); - printf(" Compute by-primitives mean temperature on a given 2D shape at a given time.\n"); - printf(" Mean temperature is computed on the front sides of the primitives\n"); - printf(" listed in SURFACE.stl. These primitives are not added to the geometry,\n"); - printf(" but must be already known. The shape doesn't need to be connex.\n"); - printf(" If some boundary conditions are time-dependant, TIME cannot be INF.\n"); - printf("\n -F SURFACE.stl,TIME:\n"); - printf(" Compute the mean flux on a given 2D shape at a given time.\n"); - printf(" Mean flux is computed on the primitives listed in SURFACE.stl and is\n"); - printf(" accounted positive, at the primitive scale, when going from the front\n"); - printf(" side to the back side. These primitives are not added to the geometry,\n"); - printf(" but must be already known. The shape doesn't need to be connex.\n"); - printf(" If some boundary conditions are time-dependant, TIME cannot be INF.\n"); + printf(" Compute the mean temperature in a given medium at a given time. The medium\n"); + printf(" doesn't need to be connex.\n"); + printf("\n -s STL_FILE_NAME,TIME:\n"); + printf(" Compute the mean temperature on a given 2D region at a given time, the\n"); + printf(" region being defined as the front sides of the triangles in the provided\n"); + printf(" file. These triangles are not added to the geometry, but must be already\n"); + printf(" part of it. The region doesn't need to be connex.\n"); + printf("\n -S STL_FILE_NAME,TIME:\n"); + printf(" Compute the by-triangle mean temperature on a given 2D region at a given\n"); + printf(" time, the region being defined as the front sides of the triangles in the\n"); + printf(" provided file. These triangles are not added to the geometry, but must be\n"); + printf(" already part of it. The region doesn't need to be connex.\n"); + printf("\n -F STL_FILE_NAME,TIME:\n"); + printf(" Compute the mean flux on a given 2D region at a given time, the region\n"); + printf(" being defined as the front sides of the triangles in the provided file.\n"); + printf(" Flux is accounted positive when going from the front side to the back side,\n"); + printf(" at a single-triangle scale. These triangles are not added to the geometry,\n"); + printf(" but must be already part of it. The region doesn't need to be connex.\n"); printf("\n -R t=TIME:spp=SPP:fov=FOV:up=XUP,YUP,ZUP:pos=X,Y,Z:tgt=XT,YT,ZT:img=WxH:\n"); printf(" The infra-red rendering mode.\n"); printf(" t is the rendering time (default: INF).\n"); printf(" spp is the number of samples per pixel (default: 4).\n"); printf(" fov is the field of view (default: 70 degrees).\n"); printf(" up is the up direction (default: 0,0,1).\n"); - printf(" pos is the position of camera (default: 1,1,1).\n"); - printf(" tgt is the target (default: 0,0,0).\n"); - printf(" img is the resolution (default: 640x480).\n"); + printf(" pos is the position of the camera (default: 1,1,1).\n"); + printf(" tgt is the target of the camera (default: 0,0,0).\n"); + printf(" img is the image resolution (default: 640x480).\n"); + printf("\nOptionnal arguments:\n"); + printf("--------------------\n"); printf("\n -d:\n"); - printf(" Dump the geometry in VTK format with medium front/back id and\n"); - printf(" boundary id.\n"); - printf("\nOptionnal arguments\n"); - printf("-------------------\n"); + printf(" Dump the geometry in VTK format with various properties, including possible\n"); + printf(" errors. If this option is used, no computation occurs. Including both this\n"); + printf(" option and one of the options specifying a compute region ("); + print_multiple_modes(stdout, REGION_COMPUTE_MODES); printf(") has\n"); + printf(" the effect to include the region in the dump.\n"); printf("\n -g:\n"); - printf(" Compute the green function.\n"); + printf(" Compute the green function at t=INF.\n"); printf(" Can only be used in conjonction with one these options:"); print_multiple_modes(stdout, GREEN_COMPATIBLE_MODES); printf("\n"); printf(" Provided TIME is silently ignored.\n"); printf("\n -f SCALE_FACTOR.\n"); printf(" Default value: 1.0.\n"); - printf("\n -D { all | error | success }:\n"); - printf(" dump paths in VTK format.\n"); + printf("\n -D { all|error|success }:\n"); + printf(" Dump paths in VTK format.\n"); printf("\n -n NUM_OF_REALIZATIONS:\n"); printf(" Number of Monte-Carlo realizations.\n"); printf(" Default value: 10000.\n"); @@ -223,6 +415,8 @@ print_help printf(" REFERENCE_TEMP is the temperature used for the linearization\n"); printf(" of the radiative transfer.\n"); printf(" Default values: 300 and 300.\n"); + printf("\n -V:\n"); + printf(" Verbose mode.\n"); } res_T @@ -237,14 +431,14 @@ parse_args ASSERT(argv && args); - if (argc == 1) { + if(argc == 1) { print_help(argv[0]); res = RES_BAD_ARG; goto error; } opterr = 0; /* No default error messages */ - while ((opt = getopt(argc, argv, "hvgn:t:B:M:m:p:P:dD:s:S:F:f:r:R:")) != -1) { + while ((opt = getopt(argc, argv, "hvgn:t:B:M:m:p:P:dD:s:S:F:f:r:R:V")) != -1) { switch (opt) { case '?': /* Unreconised option */ res = RES_BAD_ARG; @@ -254,12 +448,12 @@ parse_args case 'h': print_help(argv[0]); args->just_help = 1; - goto exit; + goto end; case 'v': print_version(); args->just_help = 1; - goto exit; + goto end; case 'g': args->mode |= GREEN_MODE; @@ -268,67 +462,71 @@ parse_args case 'n': { unsigned long n; res = cstr_to_ulong(optarg, &n); - if (res != RES_OK || n == 0) { - res = RES_BAD_ARG; + if(res != RES_OK + || n == 0) + { + if(res == RES_OK) res = RES_BAD_ARG; fprintf(stderr, "Invalid argument for option -%c: %s\n", opt, optarg); goto error; } - args->N = n; + args->samples = n; n_used = 1; break; } - case 'B': - args->bc_filename = optarg; - break; - - case 'M': - args->medium_filename = optarg; + case 'M': { + struct str name; + str_init(args->allocator, &name); + str_set(&name, optarg); + ERR(darray_str_push_back(&args->model_files, &name)); + str_release(&name); break; + } case 't': res = cstr_to_uint(optarg, &args->nthreads); - if (res != RES_OK - || args->nthreads <= 0) { - res = RES_BAD_ARG; + if(res != RES_OK + || args->nthreads <= 0) + { + if(res == RES_OK) res = RES_BAD_ARG; fprintf(stderr, "Invalid argument for option -%c: %s\n", opt, optarg); goto error; } break; case 'p': - if (args->mode & EXCLUSIVE_MODES) { + if(args->mode & EXCLUSIVE_MODES) { res = RES_BAD_ARG; fprintf(stderr, "Modes -%c and -%c are exclusive.\n", (char)opt, mode_option(args->mode)); goto error; } args->mode |= PROBE_COMPUTE; - cstr_to_list_double(optarg, ',', args->u.probe, &len, 4); - if (len != 4 + cstr_to_list_double(optarg, ',', args->probe, &len, 4); + if(len != 4 || res != RES_OK - || args->u.probe[3] < 0) + || args->probe[3] < 0) { - res = RES_BAD_ARG; + if(res == RES_OK) res = RES_BAD_ARG; fprintf(stderr, "Invalid argument for option -%c: %s\n", opt, optarg); goto error; } break; case 'P': - if (args->mode & EXCLUSIVE_MODES) { + if(args->mode & EXCLUSIVE_MODES) { res = RES_BAD_ARG; fprintf(stderr, "Modes -%c and -%c are exclusive.\n", (char)opt, mode_option(args->mode)); goto error; } args->mode |= PROBE_COMPUTE_ON_INTERFACE; - cstr_to_list_double(optarg, ',', args->u.probe, &len, 4); - if (len != 4 + cstr_to_list_double(optarg, ',', args->probe, &len, 4); + if(len != 4 || res != RES_OK - || args->u.probe[3] < 0) + || args->probe[3] < 0) { - res = RES_BAD_ARG; + if(res == RES_OK) res = RES_BAD_ARG; fprintf(stderr, "Invalid argument for option -%c: %s\n", opt, optarg); goto error; } @@ -339,7 +537,7 @@ parse_args case 'F': { int err = 0; char *ptr; - if (args->mode & EXCLUSIVE_MODES) { + if(args->mode & EXCLUSIVE_MODES) { res = RES_BAD_ARG; fprintf(stderr, "Modes -%c and -%c are exclusive.\n", (char)opt, mode_option(args->mode)); @@ -358,15 +556,15 @@ parse_args } ptr = strrchr(optarg, ','); /* Single ',' allowed */ - if (!ptr || ptr != strchr(optarg, ',')) + if(!ptr || ptr != strchr(optarg, ',')) err = 1; else { *ptr = '\0'; ptr++; args->solve_boundary_filename = optarg; - err = (RES_OK != cstr_to_double(ptr, args->u.probe + 3) || args->u.probe[3] < 0); + err = (RES_OK != cstr_to_double(ptr, args->probe + 3) || args->probe[3] < 0); } - if (err) { + if(err) { res = RES_BAD_ARG; fprintf(stderr, "Invalid argument for option -%c: %s\n", opt, optarg); goto error; @@ -376,7 +574,7 @@ parse_args case 'm': { char* ptr; - if (args->mode & EXCLUSIVE_MODES) { + if(args->mode & EXCLUSIVE_MODES) { res = RES_BAD_ARG; fprintf(stderr, "Modes -%c and -%c are exclusive.\n", (char)opt, mode_option(args->mode)); @@ -384,10 +582,10 @@ parse_args } args->mode |= MEDIUM_COMPUTE; ptr = strpbrk(optarg, ","); - if (!ptr || ptr == optarg) + if(!ptr || ptr == optarg) res = RES_BAD_ARG; - else res = cstr_to_double(ptr + 1, args->u.probe + 3); - if (res != RES_OK) { + else res = cstr_to_double(ptr + 1, args->probe + 3); + if(res != RES_OK) { fprintf(stderr, "Invalid argument for option -%c: %s\n", opt, optarg); goto error; } @@ -397,23 +595,17 @@ parse_args } case 'd': - if (args->mode & EXCLUSIVE_MODES) { - res = RES_BAD_ARG; - fprintf(stderr, "Modes -%c and -%c are exclusive.\n", - (char)opt, mode_option(args->mode)); - goto error; - } args->mode |= DUMP_VTK; break; case 'D': - if (0 == strcmp(optarg, "all")) { + if(0 == strcmp(optarg, "all")) { args->dump_paths = DUMP_ALL; } - else if (0 == strcmp(optarg, "error")) { + else if(0 == strcmp(optarg, "error")) { args->dump_paths = DUMP_ERROR; } - else if (0 == strcmp(optarg, "success")) { + else if(0 == strcmp(optarg, "success")) { args->dump_paths = DUMP_SUCCESS; } else { res = RES_BAD_ARG; @@ -424,9 +616,10 @@ parse_args case 'f': cstr_to_double(optarg, &args->scale_factor); - if (res != RES_OK - || args->scale_factor <= 0) { - res = RES_BAD_ARG; + if(res != RES_OK + || args->scale_factor <= 0) + { + if(res == RES_OK) res = RES_BAD_ARG; fprintf(stderr, "Invalid argument for option -%c: %s\n", opt, optarg); goto error; } @@ -434,18 +627,19 @@ parse_args case 'r': cstr_to_list_double(optarg, ',', args->radiative_temp, &len, 2); - if (res != RES_OK + if(res != RES_OK || len != 2 || args->radiative_temp[0] <= 0 - || args->radiative_temp[1] <= 0) { - res = RES_BAD_ARG; + || args->radiative_temp[1] <= 0) + { + if(res == RES_OK) res = RES_BAD_ARG; fprintf(stderr, "Invalid argument for option -%c: %s\n", opt, optarg); goto error; } break; case 'R': - if (args->mode & EXCLUSIVE_MODES) { + if(args->mode & EXCLUSIVE_MODES) { res = RES_BAD_ARG; fprintf(stderr, "Modes -%c and -%c are exclusive.\n", (char)opt, mode_option(args->mode)); @@ -454,10 +648,14 @@ parse_args args->mode |= IR_COMPUTE; args->camera = optarg; break; + + case 'V': + args->verbose = 1; + break; } } - if (args->mode & GREEN_MODE + if(args->mode & GREEN_MODE && (args->mode != (args->mode & (GREEN_MODE | GREEN_COMPATIBLE_MODES)))) { fprintf(stderr, "Option -g can only be used in conjunction with:"); @@ -467,30 +665,24 @@ parse_args goto error; } - if (!args->medium_filename) { - fprintf(stderr, "Missing mandatory argument: -M MEDIUM.txt\n"); + if(!darray_str_size_get(&args->model_files)) { + fprintf(stderr, "Missing mandatory argument: -M <model_file_name>\n"); res = RES_BAD_ARG; goto error; } - if (!args->bc_filename) { - fprintf(stderr, "Missing mandatory argument: -B BOUNDARY.txt\n"); - res = RES_BAD_ARG; - goto error; - } - - if (args->mode & IR_COMPUTE && n_used) { + if(args->mode & IR_COMPUTE && n_used) { fprintf(stderr, "The -n option has no effect in rendering mode;" " use rendering's SPP suboption instead.\n"); res = RES_BAD_ARG; goto error; } -exit: +end: return res; error: fprintf(stderr, "Use the -h option to get help.\n"); - goto exit; + goto end; } res_T @@ -506,41 +698,32 @@ parse_camera line = split_line(cam_param, ':'); - if (line) - { + if(line) { int i = 0; - for (i = 0; *(line + i); i++) - { + for(i = 0; *(line + i); i++) { size_t len = 0; opt = split_line(line[i], '='); _strupr(opt[0]); - if (strcmp(opt[0], "T") == 0) { - res = cstr_to_double(opt[1], &cam->u.time); - if (res != RES_OK) goto error; + if(strcmp(opt[0], "T") == 0) { + ERR(cstr_to_double(opt[1], &cam->time)); } - else if (strcmp(opt[0], "FOV") == 0) { - res = cstr_to_double(opt[1], &cam->fov); - if (res != RES_OK) goto error; + else if(strcmp(opt[0], "FOV") == 0) { + ERR(cstr_to_double(opt[1], &cam->fov)); } - else if (strcmp(opt[0], "UP") == 0) { - res = cstr_to_list_double(opt[1], ',', cam->up, &len, 3); - if (res != RES_OK || len != 3) goto error; + else if(strcmp(opt[0], "UP") == 0) { + ERR(cstr_to_list_double(opt[1], ',', cam->up, &len, 3)); } - else if (strcmp(opt[0], "TGT") == 0) { - res = cstr_to_list_double(opt[1], ',', cam->tgt, &len, 3); - if (res != RES_OK || len != 3) goto error; + else if(strcmp(opt[0], "TGT") == 0) { + ERR(cstr_to_list_double(opt[1], ',', cam->tgt, &len, 3)); } - else if (strcmp(opt[0], "POS") == 0) { - res = cstr_to_list_double(opt[1], ',', cam->pos, &len, 3); - if (res != RES_OK || len != 3) goto error; + else if(strcmp(opt[0], "POS") == 0) { + ERR(cstr_to_list_double(opt[1], ',', cam->pos, &len, 3)); } - else if (strcmp(opt[0], "IMG") == 0) { - res = cstr_to_list_uint(opt[1], 'x', cam->img, &len, 2); - if (res != RES_OK || len != 2) goto error; + else if(strcmp(opt[0], "IMG") == 0) { + ERR(cstr_to_list_uint(opt[1], 'x', cam->img, &len, 2)); } - else if (strcmp(opt[0], "SPP") == 0) { - res = cstr_to_uint(opt[1], &cam->spp); - if (res != RES_OK) goto error; + else if(strcmp(opt[0], "SPP") == 0) { + ERR(cstr_to_uint(opt[1], &cam->spp)); } else { fprintf(stderr, "Warning: unexpected option for rendering mode: %s.\n", opt[0]); @@ -550,388 +733,455 @@ parse_camera } } -exit: +end: +#define FREE_AARRAY(ARRAY) if(ARRAY) {\ + int i = 0; \ + for(i=0; *(ARRAY+i);i++){\ + free(ARRAY[i]);\ + }\ + free(ARRAY);\ +} FREE_AARRAY(line) - FREE_AARRAY(opt) + FREE_AARRAY(opt) +#undef FREE_AARRAY - return res; + return res; error: fprintf(stderr, "Use the -h option to get help.\n"); - goto exit; + goto end; } -/* Read medium line; should be one of: - * SOLID medium_name STL_filename lambda rho cp delta "Tinit(x,y,z)" [ "volumic_power(x,y,z,t)" ] - * FLUID medium_name STL_filename rho cp "Tinit(x,y,z)" -*/ -res_T -parse_medium_line - (char* line, - char** stl_filename, - struct description* desc) +/* H_BOUNDARY_FOR_SOLID Name emissivity specular_fraction hc T_env STL_filenames + * H_BOUNDARY_FOR_FLUID Name emissivity specular_fraction hc T_env STL_filenames */ +static res_T +process_h + (struct stardis* stardis, + const enum description_type type, + char** tok_ctx) { char* tk = NULL; + struct description* desc; + size_t sz; res_T res = RES_OK; - ASSERT(line && stl_filename && desc); - -#define CHK_TOK(x, Name) if ((tk = (x)) == NULL) {\ - fprintf(stderr, "Invalid data (missing token " Name ")\n");\ - res = RES_BAD_ARG;\ - goto exit;\ - } (void)0 - CHK_TOK(_strupr(strtok(line, " ")), "medium type"); - if (strcmp(tk, "SOLID") == 0) { - desc->type = DESC_MAT_SOLID; - desc->d.solid = NULL_SOLID; - - CHK_TOK(strtok(NULL, " "), "medium name"); - if (strlen(tk) >= sizeof(desc->d.solid.name)) { - fprintf(stderr, "Medium name is too long: %s\n", tk); - res = RES_BAD_ARG; - goto exit; - } - strncpy(desc->d.solid.name, tk, sizeof(desc->d.solid.name)); + ASSERT(stardis && tok_ctx); - CHK_TOK(strtok(NULL, " "), "file name"); - *stl_filename = malloc(strlen(tk) + 1); - strcpy(*stl_filename, tk); + sz = darray_descriptions_size_get(&stardis->descriptions); + ASSERT(sz <= UINT_MAX); + ERR(darray_descriptions_resize(&stardis->descriptions, sz+1)); + desc = darray_descriptions_data_get(&stardis->descriptions) + sz; + init_h(stardis->allocator, &desc->d.h_boundary); - CHK_TOK(strtok(NULL, " "), "lambda"); - res = cstr_to_double(tk, &desc->d.solid.lambda); - if (res != RES_OK || desc->d.solid.lambda <= 0) { - fprintf(stderr, "Invalid lambda: %s\n", tk); - res = RES_BAD_ARG; - goto exit; - } - CHK_TOK(strtok(NULL, " "), "rho"); - res = cstr_to_double(tk, &desc->d.solid.rho); - if (res != RES_OK || desc->d.solid.rho <= 0) { - fprintf(stderr, "Invalid rho: %s\n", tk); - res = RES_BAD_ARG; - goto exit; - } - CHK_TOK(strtok(NULL, " "), "cp"); - res = cstr_to_double(tk, &desc->d.solid.cp); - if (res != RES_OK || desc->d.solid.cp <= 0) { - fprintf(stderr, "Invalid cp: %s\n", tk); - res = RES_BAD_ARG; - goto exit; - } - CHK_TOK(strtok(NULL, " "), "delta"); - res = cstr_to_double(tk, &desc->d.solid.delta); - if (res != RES_OK || desc->d.solid.delta <= 0) { - fprintf(stderr, "Invalid delta: %s\n", tk); - res = RES_BAD_ARG; - goto exit; - } - CHK_TOK(strtok(NULL, " "), "Tinit"); - res = cstr_to_double(tk, &desc->d.solid.tinit); - if (res != RES_OK || desc->d.solid.tinit < 0) { - fprintf(stderr, "Invalid Tinit: %s\n", tk); - res = RES_BAD_ARG; - goto exit; - } + desc->d.h_boundary.mat_id = (unsigned)sz; + desc->type = type; - /* Volumic Power is optional */ - tk = strtok(NULL, " "); - if (tk) { - res = cstr_to_double(tk, &desc->d.solid.vpower); - if (res != RES_OK) { - fprintf(stderr, "Invalid volumic power: %s\n", tk); - res = RES_BAD_ARG; - goto exit; - } - } - else desc->d.solid.vpower = SDIS_VOLUMIC_POWER_NONE; + CHK_TOK(strtok_r(NULL, " ", tok_ctx), "h boundary name"); + ERR(str_set(&desc->d.h_boundary.name, tk)); + + CHK_TOK(strtok_r(NULL, " ", tok_ctx), "emissivity"); + res = cstr_to_double(tk, &desc->d.h_boundary.emissivity); + if(res != RES_OK + || desc->d.h_boundary.emissivity < 0 + || desc->d.h_boundary.emissivity > 1) + { + fprintf(stderr, "Invalid emissivity: %s\n", tk); + if(res == RES_OK) res = RES_BAD_ARG; + goto end; + } + CHK_TOK(strtok_r(NULL, " ", tok_ctx), "specular fraction"); + res = cstr_to_double(tk, &desc->d.h_boundary.specular_fraction); + if(res != RES_OK + || desc->d.h_boundary.specular_fraction < 0 + || desc->d.h_boundary.specular_fraction > 1) + { + fprintf(stderr, "Invalid specular fraction: %s\n", tk); + if(res == RES_OK) res = RES_BAD_ARG; + goto end; + } + CHK_TOK(strtok_r(NULL, " ", tok_ctx), "hc"); + res = cstr_to_double(tk, &desc->d.h_boundary.hc); + if(res != RES_OK + || desc->d.h_boundary.hc < 0) + { + fprintf(stderr, "Invalid hc: %s\n", tk); + if(res == RES_OK) res = RES_BAD_ARG; + goto end; + } + CHK_TOK(strtok_r(NULL, " ", tok_ctx), "temperature"); + res = cstr_to_double(tk, &desc->d.h_boundary.imposed_temperature); + if(res != RES_OK + || desc->d.h_boundary.imposed_temperature < 0) + { + fprintf(stderr, "Invalid temperature: %s\n", tk); + if(res == RES_OK) res = RES_BAD_ARG; + goto end; } - else if (strcmp(tk, "FLUID") == 0) { - desc->type = DESC_MAT_FLUID; - desc->d.fluid = NULL_FLUID; - CHK_TOK(strtok(NULL, " "), "medium name"); - if (strlen(tk) >= sizeof(desc->d.fluid.name)) { - fprintf(stderr, "Medium name is too long: %s\n", tk); - res = RES_BAD_ARG; - goto exit; - } - strncpy(desc->d.fluid.name, tk, sizeof(desc->d.fluid.name)); + ERR(read_sides_and_files(stardis, 1, desc->d.h_boundary.mat_id, tok_ctx)); - CHK_TOK(strtok(NULL, " "), "file name"); - *stl_filename = malloc(strlen(tk) + 1); - strcpy(*stl_filename, tk); +end: + return res; +error: + goto end; +} - CHK_TOK(strtok(NULL, " "), "rho"); - res = cstr_to_double(tk, &desc->d.fluid.rho); - if (res != RES_OK || desc->d.fluid.rho <= 0) { - fprintf(stderr, "Invalid rho: %s\n", tk); - res = RES_BAD_ARG; - goto exit; - } - CHK_TOK(strtok(NULL, " "), "cp"); - res = cstr_to_double(tk, &desc->d.fluid.cp); - if (res != RES_OK || desc->d.fluid.cp <= 0) { - fprintf(stderr, "Invalid cp: %s\n", tk); - res = RES_BAD_ARG; - goto exit; - } - /* Depending of the number of spaces before '"' strtok should - * be called once or twice */ - CHK_TOK(strtok(NULL, " "), "Tinit"); - res = cstr_to_double(tk, &desc->d.fluid.tinit); - if (res != RES_OK || desc->d.fluid.tinit < 0) { - fprintf(stderr, "Invalid Tinit: %s\n", tk); - res = RES_BAD_ARG; - goto exit; +/* T_BOUNDARY_FOR_SOLID Name T STL_filenames + * T_BOUNDARY_FOR_FLUID Name T hc STL_filenames */ +static res_T +process_t + (struct stardis* stardis, + const enum description_type type, + char** tok_ctx) +{ + char* tk = NULL; + struct description* desc; + size_t sz; + res_T res = RES_OK; + + ASSERT(stardis && tok_ctx); + + sz = darray_descriptions_size_get(&stardis->descriptions); + ASSERT(sz <= UINT_MAX); + ERR(darray_descriptions_resize(&stardis->descriptions, sz + 1)); + desc = darray_descriptions_data_get(&stardis->descriptions) + sz; + init_t(stardis->allocator, &desc->d.t_boundary); + + desc->d.t_boundary.mat_id = (unsigned)sz; + desc->type = type; + + CHK_TOK(strtok_r(NULL, " ", tok_ctx), "temperature boundary name"); + ERR(str_set(&desc->d.t_boundary.name, tk)); + + CHK_TOK(strtok_r(NULL, " ", tok_ctx), "temperature"); + res = cstr_to_double(tk, &desc->d.t_boundary.imposed_temperature); + if(res != RES_OK + || desc->d.t_boundary.imposed_temperature < 0) + { + fprintf(stderr, "Invalid temperature: %s\n", tk); + if(res == RES_OK) res = RES_BAD_ARG; + goto end; + } + if(type == DESC_BOUND_T_FOR_FLUID) { + CHK_TOK(strtok_r(NULL, " ", tok_ctx), "hc"); + res = cstr_to_double(tk, &desc->d.t_boundary.hc); + if(res != RES_OK + || desc->d.t_boundary.hc < 0) + { + fprintf(stderr, "Invalid hc: %s\n", tk); + if(res == RES_OK) res = RES_BAD_ARG; + goto end; } - } else { - res = RES_BAD_ARG; - fprintf(stderr, "Invalid medium type: %s\n", tk); } -#undef CHK_TOK - exit : - return res; + ERR(read_sides_and_files(stardis, 1, desc->d.t_boundary.mat_id, tok_ctx)); + +end: + return res; +error: + goto end; } -/* Read boundary line; should be one of: -# H_BOUNDARY_FOR_SOLID Boundary_name STL_filename emissivity specular_fraction hc hc_max "T_env(x,y,z,t)" -# | H_BOUNDARY_FOR_FLUID Boundary_name STL_filename emissivity specular_fraction hc hc_max "T_env(x,y,z,t)" -# | T_BOUNDARY_FOR_SOLID Boundary_name STL_filename "T_value(x,y,z,t)" -# | T_BOUNDARY_FOR_FLUID Boundary_name STL_filename "T_value(x,y,z,t)" hc hc_max -# | F_BOUNDARY_FOR_SOLID Boundary_name STL_filename "F_value(x,y,z,t)" -# ; no F_BOUNDARY_FOR_FLUID -# | SOLID_FLUID_CONNECTION Connection_name STL_filename emissivity specular_fraction hc hc_max - */ -res_T -parse_boundary_line - (char* line, - char** stl_filename, - struct description* desc) +/* F_BOUNDARY_FOR_SOLID Name F STL_filenames */ +static res_T +process_flx + (struct stardis* stardis, + char** tok_ctx) { char* tk = NULL; - int h_solid = 0, h_fluid = 0; - int t_solid = 0, t_fluid = 0; - int f_solid = 0; - int sf_connect = 0; + struct description* desc; + size_t sz; res_T res = RES_OK; - ASSERT(line && stl_filename && desc); - -#define CHK_TOK(x, Name) if ((tk = (x)) == NULL) {\ - fprintf(stderr, "Invalid data (missing token " Name ")\n");\ - res = RES_BAD_ARG;\ - goto exit;\ - } (void)0 - CHK_TOK(strtok(line, " "), "boundary type"); - h_solid = (0 == strcmp(tk, "H_BOUNDARY_FOR_SOLID")); - if (!h_solid) { - h_fluid = (0 == strcmp(tk, "H_BOUNDARY_FOR_FLUID")); - if (!h_fluid) { - t_solid = (0 == strcmp(tk, "T_BOUNDARY_FOR_SOLID")); - if (!t_solid) { - t_fluid = (0 == strcmp(tk, "T_BOUNDARY_FOR_FLUID")); - if (!t_fluid) { - f_solid = (0 == strcmp(tk, "F_BOUNDARY_FOR_SOLID")); - if (!f_solid) { - sf_connect = (0 == strcmp(tk, "SOLID_FLUID_CONNECTION")); - } - } - } - } + ASSERT(stardis && tok_ctx); + + sz = darray_descriptions_size_get(&stardis->descriptions); + ASSERT(sz <= UINT_MAX); + ERR(darray_descriptions_resize(&stardis->descriptions, sz + 1)); + desc = darray_descriptions_data_get(&stardis->descriptions) + sz; + init_f(stardis->allocator, &desc->d.f_boundary); + + desc->d.f_boundary.mat_id = (unsigned)sz; + desc->type = DESC_BOUND_F_FOR_SOLID; + + CHK_TOK(strtok_r(NULL, " ", tok_ctx), "flux boundary name"); + ERR(str_set(&desc->d.f_boundary.name, tk)); + + CHK_TOK(strtok_r(NULL, " ", tok_ctx), "flux"); + res = cstr_to_double(tk, &desc->d.f_boundary.imposed_flux); + if(res != RES_OK) { + /* Flux can be < 0 */ + fprintf(stderr, "Invalid flux: %s\n", tk); + goto end; } - if (h_solid || h_fluid) { - desc->type = h_solid ? DESC_BOUND_H_FOR_SOLID : DESC_BOUND_H_FOR_FLUID; - desc->d.h_boundary = NULL_HBOUND; + ERR(read_sides_and_files(stardis, 1, desc->d.f_boundary.mat_id, tok_ctx)); - CHK_TOK(strtok(NULL, " "), "boundary name"); - if (strlen(tk) >= sizeof(desc->d.h_boundary.name)) { - fprintf(stderr, "Boundary name is too long: %s\n", tk); - res = RES_BAD_ARG; - goto exit; - } - strncpy(desc->d.h_boundary.name, tk, sizeof(desc->d.h_boundary.name)); - - /* The description is parsed the same way for both types */ - CHK_TOK(strtok(NULL, " "), "file name"); - *stl_filename = malloc(strlen(tk) + 1); - strcpy(*stl_filename, tk); - - CHK_TOK(strtok(NULL, " "), "emissivity"); - res = cstr_to_double(tk, &desc->d.h_boundary.emissivity); - if (res != RES_OK || - desc->d.h_boundary.emissivity < 0 || desc->d.h_boundary.emissivity > 1) { - fprintf(stderr, "Invalid emissivity: %s\n", tk); - res = RES_BAD_ARG; - goto exit; - } - desc->d.h_boundary.has_emissivity = (desc->d.h_boundary.emissivity > 0); - CHK_TOK(strtok(NULL, " "), "specular fraction"); - res = cstr_to_double(tk, &desc->d.h_boundary.specular_fraction); - if (res != RES_OK - || desc->d.h_boundary.specular_fraction < 0.0 - || desc->d.h_boundary.specular_fraction > 1.0) { - fprintf(stderr, "Invalid specular_fraction: %s\n", tk); - res = RES_BAD_ARG; - goto exit; - } - CHK_TOK(strtok(NULL, " "), "hc"); - res = cstr_to_double(tk, &desc->d.h_boundary.hc); - if (res != RES_OK || desc->d.h_boundary.hc < 0) { - fprintf(stderr, "Invalid hc value: %s\n", tk); - res = RES_BAD_ARG; - goto exit; - } - desc->d.h_boundary.has_hc = (desc->d.h_boundary.hc > 0); - CHK_TOK(strtok(NULL, " "), "hc max"); - res = cstr_to_double(tk, &desc->d.h_boundary.hc_max); - if (res != RES_OK - || desc->d.h_boundary.hc_max < desc->d.h_boundary.hc - || desc->d.h_boundary.hc_max < 0) { - fprintf(stderr, "Invalid hc_max value: %s\n", tk); - res = RES_BAD_ARG; - goto exit; - } - CHK_TOK(strtok(NULL, " "), "temperature"); - res = cstr_to_double(tk, &desc->d.h_boundary.imposed_temperature); - if (res != RES_OK || desc->d.h_boundary.imposed_temperature < 0) { - fprintf(stderr, "Invalid temperature value: %s\n", tk); - res = RES_BAD_ARG; - goto exit; - } +end: + return res; +error: + goto end; +} + +/* SOLID_FLUID_CONNECTION Name emissivity specular_fraction hc STL_filenames */ +static res_T +process_sfc + (struct stardis* stardis, + char** tok_ctx) +{ + char* tk = NULL; + struct description* desc; + size_t sz; + unsigned connection_id; + res_T res = RES_OK; + + ASSERT(stardis && tok_ctx); + + sz = darray_descriptions_size_get(&stardis->descriptions); + ASSERT(sz <= UINT_MAX); + ERR(darray_descriptions_resize(&stardis->descriptions, sz + 1)); + desc = darray_descriptions_data_get(&stardis->descriptions) + sz; + init_sf(stardis->allocator, &desc->d.sf_connect); + + connection_id = (unsigned)sz; + desc->type = DESC_SOLID_FLUID_CONNECT; + + CHK_TOK(strtok_r(NULL, " ", tok_ctx), "connection name"); + ERR(str_set(&desc->d.sf_connect.name, tk)); + + CHK_TOK(strtok_r(NULL, " ", tok_ctx), "emissivity"); + res = cstr_to_double(tk, &desc->d.sf_connect.emissivity); + if(res != RES_OK + || desc->d.sf_connect.emissivity < 0 + || desc->d.h_boundary.emissivity > 1) + { + fprintf(stderr, "Invalid emissivity: %s\n", tk); + if(res == RES_OK) res = RES_BAD_ARG; + goto end; + } + CHK_TOK(strtok_r(NULL, " ", tok_ctx), "specular fraction"); + res = cstr_to_double(tk, &desc->d.sf_connect.specular_fraction); + if(res != RES_OK + || desc->d.sf_connect.specular_fraction < 0 + || desc->d.sf_connect.specular_fraction > 1) + { + fprintf(stderr, "Invalid specular fraction: %s\n", tk); + if(res == RES_OK) res = RES_BAD_ARG; + goto end; + } + CHK_TOK(strtok_r(NULL, " ", tok_ctx), "hc"); + res = cstr_to_double(tk, &desc->d.sf_connect.hc); + if(res != RES_OK + || desc->d.sf_connect.hc <= 0) + { + fprintf(stderr, "Invalid hc: %s\n", tk); + if(res == RES_OK) res = RES_BAD_ARG; + goto end; } - else if (t_solid || t_fluid) { - desc->type = t_solid ? DESC_BOUND_T_FOR_SOLID : DESC_BOUND_T_FOR_FLUID; - desc->d.t_boundary = NULL_TBOUND; - CHK_TOK(strtok(NULL, " "), "boundary name"); - if (strlen(tk) >= sizeof(desc->d.t_boundary.name)) { - fprintf(stderr, "Boundary name is too long: %s\n", tk); - res = RES_BAD_ARG; - goto exit; - } - strncpy(desc->d.t_boundary.name, tk, sizeof(desc->d.t_boundary.name)); + ERR(read_sides_and_files(stardis, 1, connection_id, tok_ctx)); - /* This part of the description is parsed the same way for both types */ - CHK_TOK(strtok(NULL, " "), "file name"); - *stl_filename = malloc(strlen(tk) + 1); - strcpy(*stl_filename, tk); +end: + return res; +error: + goto end; +} - CHK_TOK(strtok(NULL, " "), "temperature"); - res = cstr_to_double(tk, &desc->d.t_boundary.imposed_temperature); - if (res != RES_OK || desc->d.t_boundary.imposed_temperature < 0) { - fprintf(stderr, "Invalid temperature value: %s\n", tk); - res = RES_BAD_ARG; - goto exit; - } +/* SOLID Name lambda rho cp delta Tinit volumic_power STL_filenames */ +static res_T +process_solid + (struct stardis* stardis, + char** tok_ctx) +{ + char* tk = NULL; + struct description* desc; + size_t sz; + res_T res = RES_OK; - /* Parse hc + hc_max only if fluid */ - if (t_fluid) { - CHK_TOK(strtok(NULL, " "), "hc"); - res = cstr_to_double(tk, &desc->d.t_boundary.hc); - if (res != RES_OK || desc->d.t_boundary.hc < 0) { - fprintf(stderr, "Invalid hc value: %s\n", tk); - res = RES_BAD_ARG; - goto exit; - } - desc->d.t_boundary.has_hc = (desc->d.t_boundary.hc > 0); - CHK_TOK(strtok(NULL, " "), "hc max"); - res = cstr_to_double(tk, &desc->d.t_boundary.hc_max); - if (res != RES_OK - || desc->d.t_boundary.hc_max < desc->d.t_boundary.hc - || desc->d.t_boundary.hc_max < 0) { - fprintf(stderr, "Invalid hc_max value: %s\n", tk); - res = RES_BAD_ARG; - goto exit; - } - } - else desc->d.t_boundary.has_hc = 0; + ASSERT(stardis && tok_ctx); + + sz = darray_descriptions_size_get(&stardis->descriptions); + ASSERT(sz <= UINT_MAX); + ERR(darray_descriptions_resize(&stardis->descriptions, sz + 1)); + desc = darray_descriptions_data_get(&stardis->descriptions) + sz; + init_solid(stardis->allocator, &desc->d.solid); + + desc->d.solid.solid_id = (unsigned)sz; + desc->type = DESC_MAT_SOLID; + + CHK_TOK(strtok_r(NULL, " ", tok_ctx), "solid name"); + ERR(str_set(&desc->d.solid.name, tk)); + + CHK_TOK(strtok_r(NULL, " ", tok_ctx), "lambda"); + res = cstr_to_double(tk, &desc->d.solid.lambda); + if(res != RES_OK + || desc->d.solid.lambda <= 0) + { + fprintf(stderr, "Invalid lambda: %s\n", tk); + if(res == RES_OK) res = RES_BAD_ARG; + goto end; + } + CHK_TOK(strtok_r(NULL, " ", tok_ctx), "rho"); + res = cstr_to_double(tk, &desc->d.solid.rho); + if(res != RES_OK + || desc->d.solid.rho <= 0) + { + fprintf(stderr, "Invalid rho: %s\n", tk); + if(res == RES_OK) res = RES_BAD_ARG; + goto end; + } + CHK_TOK(strtok_r(NULL, " ", tok_ctx), "cp"); + res = cstr_to_double(tk, &desc->d.solid.cp); + if(res != RES_OK + || desc->d.solid.cp <= 0) + { + fprintf(stderr, "Invalid cp: %s\n", tk); + if(res == RES_OK) res = RES_BAD_ARG; + goto end; + } + CHK_TOK(strtok_r(NULL, " ", tok_ctx), "delta"); + res = cstr_to_double(tk, &desc->d.solid.delta); + if(res != RES_OK + || desc->d.solid.delta <= 0) + { + fprintf(stderr, "Invalid delta: %s\n", tk); + if(res == RES_OK) res = RES_BAD_ARG; + goto end; + } + CHK_TOK(strtok_r(NULL, " ", tok_ctx), "Tinit"); + res = cstr_to_double(tk, &desc->d.solid.tinit); + if(res != RES_OK + || desc->d.solid.tinit < 0) + { + fprintf(stderr, "Invalid Tinit: %s\n", tk); + if(res == RES_OK) res = RES_BAD_ARG; + goto end; + } + CHK_TOK(strtok_r(NULL, " ", tok_ctx), "volumic power"); + res = cstr_to_double(tk, &desc->d.solid.vpower); + if(res != RES_OK) { + /* VPower can be < 0 */ + fprintf(stderr, "Invalid volumic power: %s\n", tk); + goto end; } - else if (f_solid) { - desc->type = DESC_BOUND_F_FOR_SOLID; - desc->d.f_boundary = NULL_FBOUND; - CHK_TOK(strtok(NULL, " "), "boundary name"); - if (strlen(tk) >= sizeof(desc->d.f_boundary.name)) { - fprintf(stderr, "Boundary name is too long: %s\n", tk); - res = RES_BAD_ARG; - goto exit; - } - strncpy(desc->d.f_boundary.name, tk, sizeof(desc->d.f_boundary.name)); - - CHK_TOK(strtok(NULL, " "), "file name"); - *stl_filename = malloc(strlen(tk) + 1); - strcpy(*stl_filename, tk); - - /* Depending of the number of spaces before '"' strtok should - * be called once or twice */ - CHK_TOK(strtok(NULL, " "), "flux"); - res = cstr_to_double(tk, &desc->d.f_boundary.imposed_flux); - if (res != RES_OK) { - fprintf(stderr, "Invalid flux value: %s\n", tk); - res = RES_BAD_ARG; - goto exit; - } + ERR(read_sides_and_files(stardis, 0, desc->d.solid.solid_id, tok_ctx)); + +end: + return res; +error: + goto end; +} + +/* FLUID Name rho cp Tinit STL_filenames */ +static res_T +process_fluid + (struct stardis* stardis, + char** tok_ctx) +{ + char* tk = NULL; + struct description* desc; + size_t sz; + res_T res = RES_OK; + + ASSERT(stardis && tok_ctx); + + sz = darray_descriptions_size_get(&stardis->descriptions); + ASSERT(sz <= UINT_MAX); + ERR(darray_descriptions_resize(&stardis->descriptions, sz + 1)); + desc = darray_descriptions_data_get(&stardis->descriptions) + sz; + init_fluid(stardis->allocator, &desc->d.fluid); + + desc->d.fluid.fluid_id = (unsigned)sz; + desc->type = DESC_MAT_FLUID; + + CHK_TOK(strtok_r(NULL, " ", tok_ctx), "fluid name"); + ERR(str_set(&desc->d.fluid.name, tk)); + + CHK_TOK(strtok_r(NULL, " ", tok_ctx), "rho"); + res = cstr_to_double(tk, &desc->d.fluid.rho); + if(res != RES_OK + || desc->d.fluid.rho <= 0) + { + fprintf(stderr, "Invalid rho: %s\n", tk); + if(res == RES_OK) res = RES_BAD_ARG; + goto end; + } + CHK_TOK(strtok_r(NULL, " ", tok_ctx), "cp"); + res = cstr_to_double(tk, &desc->d.fluid.cp); + if(res != RES_OK + || desc->d.fluid.cp <= 0) + { + fprintf(stderr, "Invalid cp: %s\n", tk); + if(res == RES_OK) res = RES_BAD_ARG; + goto end; + } + CHK_TOK(strtok_r(NULL, " ", tok_ctx), "Tinit"); + res = cstr_to_double(tk, &desc->d.fluid.tinit); + if(res != RES_OK + || desc->d.fluid.tinit < 0) + { + fprintf(stderr, "Invalid Tinit: %s\n", tk); + if(res == RES_OK) res = RES_BAD_ARG; + goto end; } - else if (sf_connect) { - desc->type = DESC_SOLID_FLUID_CONNECT; - desc->d.sf_connect = NULL_SFCONNECT; - CHK_TOK(strtok(NULL, " "), "boundary name"); - if (strlen(tk) >= sizeof(desc->d.sf_connect.name)) { - fprintf(stderr, "Boundary name is too long: %s\n", tk); - res = RES_BAD_ARG; - goto exit; - } - strncpy(desc->d.sf_connect.name, tk, sizeof(desc->d.sf_connect.name)); + ERR(read_sides_and_files(stardis, 0, desc->d.fluid.fluid_id, tok_ctx)); - CHK_TOK(strtok(NULL, " "), "file name"); - *stl_filename = malloc(strlen(tk) + 1); - strcpy(*stl_filename, tk); +end: + return res; +error: + goto end; +} - CHK_TOK(strtok(NULL, " "), "emissivity"); - res = cstr_to_double(tk, &desc->d.sf_connect.emissivity); - if (res != RES_OK || - desc->d.sf_connect.emissivity < 0 || desc->d.sf_connect.emissivity > 1) { - fprintf(stderr, "Invalid emissivity: %s\n", tk); - res = RES_BAD_ARG; - goto exit; - } - desc->d.sf_connect.has_emissivity = (desc->d.sf_connect.emissivity > 0); - CHK_TOK(strtok(NULL, " "), "specular fraction"); - res = cstr_to_double(tk, &desc->d.sf_connect.specular_fraction); - if (res != RES_OK - || desc->d.sf_connect.specular_fraction < 0.0 - || desc->d.sf_connect.specular_fraction > 1.0) { - fprintf(stderr, "Invalid specular_fraction: %s\n", tk); - res = RES_BAD_ARG; - goto exit; - } - CHK_TOK(strtok(NULL, " "), "hc"); - res = cstr_to_double(tk, &desc->d.sf_connect.hc); - if (res != RES_OK || desc->d.sf_connect.hc < 0) { - fprintf(stderr, "Invalid hc value: %s\n", tk); - res = RES_BAD_ARG; - goto exit; - } - desc->d.sf_connect.has_hc = (desc->d.sf_connect.hc > 0); - CHK_TOK(strtok(NULL, " "), "hc max"); - res = cstr_to_double(tk, &desc->d.sf_connect.hc_max); - if (res != RES_OK - || desc->d.sf_connect.hc_max < desc->d.sf_connect.hc - || desc->d.sf_connect.hc_max < 0) { - fprintf(stderr, "Invalid hc_max value: %s\n", tk); - res = RES_BAD_ARG; - goto exit; - } - } else { - fprintf(stderr, "Invalid boundary type: %s", tk); +/* Read medium or boundary line; should be one of: + * SOLID Name lambda rho cp delta Tinit volumic_power STL_filenames + * FLUID Name rho cp Tinit STL_filenames + * H_BOUNDARY_FOR_SOLID Name emissivity specular_fraction hc T_env STL_filenames + * H_BOUNDARY_FOR_FLUID Name emissivity specular_fraction hc T_env STL_filenames + * T_BOUNDARY_FOR_SOLID Name T STL_filenames + * T_BOUNDARY_FOR_FLUID Name T hc STL_filenames + * F_BOUNDARY_FOR_SOLID Name F STL_filenames + * SOLID_FLUID_CONNECTION Name emissivity specular_fraction hc STL_filenames + * + * STL_filenames = { { FRONT | BACK | BOTH } STL_filename }+ + */ +res_T +process_model_line + (char* line, + struct stardis* stardis) +{ + res_T res = RES_OK; + char* tk = NULL, * tok_ctx = NULL; + + ASSERT(line && stardis); + + CHK_TOK(strtok_r(line, " ", &tok_ctx), "model line type"); + _strupr(tk); + + if(0 == strcmp(tk, "H_BOUNDARY_FOR_SOLID")) + ERR(process_h(stardis, DESC_BOUND_H_FOR_SOLID, &tok_ctx)); + else if(0 == strcmp(tk, "H_BOUNDARY_FOR_FLUID")) + ERR(process_h(stardis, DESC_BOUND_H_FOR_FLUID, &tok_ctx)); + else if(0 == strcmp(tk, "T_BOUNDARY_FOR_SOLID")) + ERR(process_t(stardis, DESC_BOUND_T_FOR_SOLID, &tok_ctx)); + else if(0 == strcmp(tk, "T_BOUNDARY_FOR_FLUID")) + ERR(process_t(stardis, DESC_BOUND_T_FOR_FLUID, &tok_ctx)); + else if(0 == strcmp(tk, "F_BOUNDARY_FOR_SOLID")) + ERR(process_flx(stardis, &tok_ctx)); + else if(0 == strcmp(tk, "SOLID_FLUID_CONNECTION")) + ERR(process_sfc(stardis, &tok_ctx)); + else if(0 == strcmp(tk, "SOLID")) + ERR(process_solid(stardis, &tok_ctx)); + else if(0 == strcmp(tk, "FLUID")) + ERR(process_fluid(stardis, &tok_ctx)); + else { + fprintf(stderr, "Unknown description type: %s\n", tk); res = RES_BAD_ARG; - goto exit; + goto error; } -#undef CHK_TOK - exit : - return res; -} +end: + return res; +error: + goto end; +} +\ No newline at end of file diff --git a/src/stardis-parsing.h b/src/stardis-parsing.h @@ -6,9 +6,25 @@ #include <sdis.h> #include <rsys/rsys.h> +#include <rsys/dynamic_array_str.h> + +#include <star/sstl.h> struct camera; struct description; +struct mem_allocator; +struct stardis; + +/* Utility macros */ +#define CHK_TOK(x, Name) if((tk = (x)) == NULL) {\ + fprintf(stderr, "Invalid data (missing token '" Name "')\n");\ + res = RES_BAD_ARG;\ + goto error;\ + } (void)0 + +#ifdef COMPILER_CL +#define strtok_r strtok_s +#endif enum stardis_mode { UNDEF_MODE = 0, @@ -25,11 +41,13 @@ enum stardis_mode { GREEN_COMPATIBLE_MODES = PROBE_COMPUTE | PROBE_COMPUTE_ON_INTERFACE | MEDIUM_COMPUTE | BOUNDARY_COMPUTE, - COMPUTE_MODES = GREEN_COMPATIBLE_MODES | IR_COMPUTE | MAP_COMPUTE | FLUX_BOUNDARY_COMPUTE, + REGION_COMPUTE_MODES = BOUNDARY_COMPUTE | FLUX_BOUNDARY_COMPUTE | MAP_COMPUTE, + + COMPUTE_MODES = GREEN_COMPATIBLE_MODES | IR_COMPUTE | REGION_COMPUTE_MODES, NON_COMPUTE_MODES = UNDEF_MODE | DUMP_VTK | GREEN_MODE, - EXCLUSIVE_MODES = COMPUTE_MODES | DUMP_VTK | FLUX_BOUNDARY_COMPUTE + EXCLUSIVE_MODES = COMPUTE_MODES | FLUX_BOUNDARY_COMPUTE }; STATIC_ASSERT(GREEN_COMPATIBLE_MODES == (COMPUTE_MODES & GREEN_COMPATIBLE_MODES), @@ -43,35 +61,60 @@ enum dump_path_type { }; struct args { - char* medium_filename; - char* bc_filename; + struct mem_allocator* allocator; + struct darray_str model_files; char* medium_name; char* solve_boundary_filename; - size_t N; + size_t samples; unsigned nthreads; - union u { /* Trick to allow static initialization with INF */ - struct pt { double p[3]; uint64_t t; } pt; - double probe[4]; - } u; + double probe[4]; enum stardis_mode mode; double scale_factor; double radiative_temp[2]; char* camera; enum dump_path_type dump_paths; int just_help; + int verbose; }; +/* Same ctx used for both media and interface add (some unused parts) */ +struct add_geom_ctx { + struct sstl_desc stl_desc; + unsigned properties[3]; + void* custom; +}; + +/* Possible callbacks for sg3d_geometry_add calls + * when void* context is a struct add_geom_ctx */ +extern LOCAL_SYM void +add_geom_ctx_indices + (const unsigned itri, + unsigned ids[3], + void* context); + +extern LOCAL_SYM void +add_geom_ctx_properties + (const unsigned itri, + unsigned prop[3], + void* context); + +extern LOCAL_SYM void +add_geom_ctx_position + (const unsigned ivert, + double pos[3], + void* context); + #ifdef NDEBUG #define DEFAULT_NTHREADS SDIS_NTHREADS_DEFAULT #else #define DEFAULT_NTHREADS 1 #endif -#define ARGS_DEFAULT__ {\ - NULL, NULL, NULL, NULL, 10000, DEFAULT_NTHREADS,\ - { { {0.0, 0.0, 0.0}, 0x7FF0000000000000 /* probe[3]=INF */}},\ - UNDEF_MODE, 1.0, {300.0, 300.0}, NULL, DUMP_NONE, 0} -static const struct args ARGS_DEFAULT = ARGS_DEFAULT__; +extern LOCAL_SYM res_T +init_args(struct mem_allocator* mem, struct args** args); + +extern LOCAL_SYM void +release_args(struct args* args); extern char mode_option @@ -84,7 +127,7 @@ print_multiple_modes extern void print_help - (char* prog); + (const char* prog); extern res_T parse_args @@ -97,16 +140,9 @@ parse_camera (char* cam_param, struct camera* cam); -extern res_T -parse_medium_line - (char* line, - char** stl_filename, - struct description* desc); - -extern res_T -parse_boundary_line +extern LOCAL_SYM res_T +process_model_line (char* line, - char** stl_filename, - struct description* desc); + struct stardis* stardis); #endif /*ARGS_H*/ diff --git a/src/stardis-solid.c b/src/stardis-solid.c @@ -1,8 +1,8 @@ #include "stardis-solid.h" #include "stardis-compute.h" +#include "stardis-app.h" #include <sdis.h> -#include<rsys/stretchy_array.h> #include <limits.h> @@ -67,23 +67,24 @@ solid_get_temperature struct sdis_data* data) { const struct solid* solid_props = sdis_data_cget(data); - if (solid_props->is_green || vtx->time > solid_props->t0) { + if(solid_props->is_green || vtx->time > solid_props->t0) { /* Always use temp for Green mode, regardless of time */ - if (solid_props->imposed_temperature >= 0) + if(solid_props->imposed_temperature >= 0) return solid_props->imposed_temperature; else return -1; } /* Time is <= 0: use tinit */ - if (solid_props->tinit >= 0) return solid_props->tinit; + if(solid_props->tinit >= 0) return solid_props->tinit; /* Must have had t_init defined: error! */ - if (solid_props->name[0]) { + if(str_is_empty(&solid_props->name)) { + FATAL("solid_get_temperature: getting undefined Tinit\n"); + } else { fprintf(stderr, "solid_get_temperature: getting undefined Tinit (solid '%s')\n", - solid_props->name); + str_cget(&solid_props->name)); ASSERT(0); abort(); - } else - FATAL("solid_get_temperature: getting undefined Tinit\n"); + } } static double @@ -103,7 +104,8 @@ solid_get_power res_T create_solid (struct sdis_device* dev, - const char* name, + struct mem_allocator* allocator, + const struct str* name, const double lambda, const double rho, const double cp, @@ -113,7 +115,7 @@ create_solid const double tinit, const double imposed_temperature, const double vpower, - struct sdis_medium*** media_ptr, + struct darray_media_ptr* media_ptr, unsigned* out_id) { res_T res = RES_OK; @@ -133,15 +135,14 @@ create_solid solid_shader.delta_boundary = solid_get_delta_boundary; #endif solid_shader.temperature = solid_get_temperature; - res = sdis_data_create(dev, sizeof(struct solid), ALIGNOF(struct solid), - NULL, &data); - if (res != RES_OK) goto error; + ERR(sdis_data_create(dev, sizeof(struct solid), ALIGNOF(struct solid), + NULL, &data)); - sz = sa_size(*media_ptr); + sz = darray_media_ptr_size_get(media_ptr); ASSERT(sz < INT_MAX); solid_props = sdis_data_get(data); /* Fetch the allocated memory space */ - if (name) strncpy(solid_props->name, name, sizeof(solid_props->name)); - else solid_props->name[0] = '\0'; + str_init(allocator, &solid_props->name); + if(name) str_copy(&solid_props->name, name); solid_props->lambda = lambda; solid_props->rho = rho; solid_props->cp = cp; @@ -152,13 +153,14 @@ create_solid solid_props->is_green = is_green; solid_props->is_outside = is_outside; solid_props->id = (unsigned)sz; - if (vpower != 0) solid_shader.volumic_power = solid_get_power; - res = sdis_solid_create(dev, &solid_shader, data, sa_add(*media_ptr, 1)); - if (res != RES_OK) goto error; + if(vpower != 0) solid_shader.volumic_power = solid_get_power; + ERR(darray_media_ptr_resize(media_ptr, sz + 1)); + ERR(sdis_solid_create(dev, &solid_shader, data, + darray_media_ptr_data_get(media_ptr) + sz)); *out_id = solid_props->id; end: - if (data) SDIS(data_ref_put(data)); + if(data) SDIS(data_ref_put(data)); return res; error: goto end; diff --git a/src/stardis-solid.h b/src/stardis-solid.h @@ -4,15 +4,16 @@ #define SDIS_SOLID_H #include <rsys/rsys.h> +#include <rsys/str.h> struct sdis_device; -struct sdis_medium; +struct darray_media_ptr; /******************************************************************************* * Solid data ******************************************************************************/ struct solid { - char name[32]; + struct str name; double cp; /* Calorific capacity */ double lambda; /* Conductivity */ double rho; /* Volumic mass */ @@ -30,7 +31,8 @@ struct solid { extern res_T create_solid (struct sdis_device* dev, - const char* name, + struct mem_allocator* allocator, + const struct str* name, const double lambda, const double rho, const double cp, @@ -40,7 +42,7 @@ create_solid const double tinit, const double imposed_temperature, const double vpower, - struct sdis_medium*** media_ptr, + struct darray_media_ptr* media_ptr, unsigned* out_id); #endif