commit 53a4cdc26d934088a1ec5f43c1c9dd8685e5b234
parent 9d793c5596a1fcdc56422a481ad4d564362eaa36
Author: Christophe Coustet <christophe.coustet@meso-star.com>
Date: Fri, 21 Dec 2018 11:37:14 +0100
For saving purposes; has memleaks
Diffstat:
6 files changed, 1650 insertions(+), 527 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -72,8 +72,15 @@ if(CMAKE_COMPILER_IS_GNUCC)
set(MATH_LIB m)
endif()
+# Default is right to left pow and log is natural log
+# Build tinyExpr without closure support and with function_1 to function_3 support only
+# (these ones cannot be disabled as predefined functions require them).
+set(TE_DEFAULTS
+ "-DTE_POW_FROM_RIGHT -DTE_NAT_LOG -DTE_WITHOUT_CLOSURES -DTE_WITHOUT_FUNCTION_0 -DTE_MAX_FUNCTION_ARITY=3")
+
ADD_LIBRARY(tinyexpr STATIC
${TINYEXPR_SOURCE_DIR}/tinyexpr.c ${TINYEXPR_SOURCE_DIR}/tinyexpr.h)
+set_target_properties(tinyexpr PROPERTIES COMPILE_FLAGS ${TE_DEFAULTS})
if(CMAKE_COMPILER_IS_GNUCC)
set_target_properties(tinyexpr PROPERTIES COMPILE_FLAGS "-std=c99")
endif()
@@ -91,6 +98,7 @@ if(MSVC)
endif()
set_target_properties(sdis-app PROPERTIES
+ COMPILE_FLAGS ${TE_DEFAULTS}
VERSION ${VERSION})
rcmake_copy_runtime_libraries(sdis-app)
diff --git a/src/args.h b/src/args.h
@@ -31,7 +31,7 @@ struct args{
char* camera;
};
#define ARGS_DEFAULT__ {\
- NULL, NULL, 10000, SDIS_NTHREADS_DEFAULT,\
+ NULL, NULL, 10000, 1,\
{ { {0.0, 0.0, 0.0}, 0x7FF0000000000000 /* probe[3]=INF */}},\
PROBE_COMPUTE, 1.0, {300.0, 300.0}, NULL}
static const struct args ARGS_DEFAULT = ARGS_DEFAULT__;
diff --git a/src/main.c b/src/main.c
@@ -20,7 +20,7 @@ int main(int argc, char** argv){
if (res != RES_OK) goto error;
if (args.mode == DUMP_VTK){
- res = dump_vtk(stdout, &stardis.geometry, stardis.geometry.boundary_count);
+ res = dump_vtk(stdout, &stardis.geometry);
if (res != RES_OK) goto error;
goto exit;
}
@@ -30,7 +30,11 @@ int main(int argc, char** argv){
exit:
stardis_release(&stardis);
+ mem_shutdown_proxy_allocator(&stardis.allocator);
if((memsz = mem_allocated_size()) != 0) {
+ char dump[4096] = { '\0' };
+ MEM_DUMP(&stardis.allocator, dump, sizeof(dump));
+ fprintf(stderr, "%s\n", dump);
fprintf(stderr, "Memory leaks: %lu Bytes\n", (unsigned long)memsz);
}
diff --git a/src/stardis-app.c b/src/stardis-app.c
@@ -4,8 +4,14 @@
#include <string.h>
#include "stardis-app.h"
-#define MEDIUM 0
-#define BOUNDARY 1
+
+#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
+#include <rsys/hash_table.h>
/*******************************************************************************
*
@@ -93,125 +99,315 @@ static char** split_line(char* a_str, const char a_delim)
return result;
}
+#ifdef COMPILER_GCC
+static char*
+_strupr(char* s)
+{
+ char* tmp;
+ for (tmp = s; *tmp; ++tmp) *tmp = toupper((unsigned char)*tmp);
+ return s;
+}
+#endif
+
+/* Read medium line; should be one of:
+ * SOLID STL_filename lambda rho cp delta "Tinit(x,y,z)" "volumic_power(x,y,z,t)"
+ * FLUID STL_filename rho cp "Tinit(x,y,z)"
+*/
static res_T
-parse_medium_line(char* line, char** stl_filename, struct material* mat)
+parse_medium_line(char* line, char** stl_filename, struct description* desc)
{
char* tk = NULL;
res_T res = RES_OK;
- tk = strtok(line," ");
- *stl_filename = malloc(strlen(tk) + 1);
- strcpy(*stl_filename,tk);
- tk = strtok(NULL," ");
- res = cstr_to_double(tk, &mat->lambda);
- if (res != RES_OK || mat->lambda <= 0) {
- fprintf(stderr, "Invalid lambda: %s\n", tk);
- res = RES_BAD_ARG;
- goto error;
- }
- tk = strtok(NULL," ");
- res = cstr_to_double(tk, &mat->rho);
- if (res != RES_OK || mat->rho <= 0) {
- fprintf(stderr, "Invalid rho: %s\n", tk);
- res = RES_BAD_ARG;
- goto error;
+#define CHK_TOK(tk, x) if(((tk) = (x)) == NULL) goto error
+ CHK_TOK(tk, _strupr(strtok(line, " ")));
+ if (strcmp(tk, "SOLID") == 0) {
+ desc->type = DESC_MAT_SOLID;
+
+ CHK_TOK(tk, strtok(NULL, " "));
+ *stl_filename = malloc(strlen(tk) + 1);
+ strcpy(*stl_filename, tk);
+
+ CHK_TOK(tk, strtok(NULL, " "));
+ 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 error;
+ }
+ CHK_TOK(tk, strtok(NULL, " "));
+ 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 error;
+ }
+ CHK_TOK(tk, strtok(NULL, " "));
+ 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 error;
+ }
+ CHK_TOK(tk, strtok(NULL, " "));
+ 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 error;
+ }
+ CHK_TOK(tk, strtok(NULL, "\""));
+ desc->d.solid.Tinit = malloc(strlen(tk) + 1);
+ strcpy(desc->d.solid.Tinit, tk);
+ CHK_TOK(tk, strtok(NULL, "\""));
+
+ CHK_TOK(tk, strtok(NULL, "\""));
+ desc->d.solid.power = malloc(strlen(tk) + 1);
+ strcpy(desc->d.solid.power, tk);
}
- tk = strtok(NULL," ");
- res = cstr_to_double(tk, &mat->cp);
- if (res != RES_OK || mat->cp <= 0){
- fprintf(stderr, "Invalid cp: %s\n", tk);
- res = RES_BAD_ARG;
- goto error;
+ else if (strcmp(tk, "FLUID") == 0) {
+ desc->type = DESC_MAT_FLUID;
+
+ CHK_TOK(tk, strtok(NULL, " "));
+ *stl_filename = malloc(strlen(tk) + 1);
+ strcpy(*stl_filename, tk);
+
+ CHK_TOK(tk, strtok(NULL, " "));
+ 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 error;
+ }
+ CHK_TOK(tk, strtok(NULL, " "));
+ 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 error;
+ }
+ CHK_TOK(tk, strtok(NULL, "\""));
+ desc->d.fluid.Tinit = malloc(strlen(tk) + 1);
+ strcpy(desc->d.fluid.Tinit, tk);
}
- tk = strtok(NULL," ");
- res = cstr_to_double(tk, &mat->delta);
- if (res != RES_OK || mat->delta <= 0){
- fprintf(stderr, "Invalid delta: %s\n", tk);
+ else {
res = RES_BAD_ARG;
- goto error;
+ printf("Invalid medium type: %s\n", tk);
}
- tk = strtok(NULL,"\"");
- tk = strtok(NULL,"\"");
- mat->Tinit = malloc(strlen(tk) + 1);
- strcpy(mat->Tinit,tk);
-
- tk = strtok(NULL,"\"");
- tk = strtok(NULL,"\"");
- mat->power = malloc(strlen(tk) + 1);
- strcpy(mat->power,tk);
+#undef GET_TOK
exit:
return res;
error:
goto exit;
}
+/* Read boundary line; should be one of:
+# H_BOUNDARY_FOR_SOLID STL_filename emissivity specular_fraction hc hc_max "T_env(x,y,z,t)"
+# | H_BOUNDARY_FOR_FLUID STL_filename emissivity specular_fraction hc hc_max "T_env(x,y,z,t)"
+# | T_BOUNDARY_FOR_SOLID STL_filename "T_value(x,y,z,t)"
+# | T_BOUNDARY_FOR_FLUID STL_filename "T_value(x,y,z,t)" hc hc_max
+# | F_BOUNDARY_FOR_SOLID STL_filename "F_value(x,y,z,t)"
+# ; no F_BOUNDARY_FOR_FLUID
+ */
static res_T
-parse_boundary_line(char* line, char** stl_filename, struct boundary* bound)
+parse_boundary_line(char* line, char** stl_filename, struct description* desc)
{
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;
res_T res = RES_OK;
- tk = strtok(line," ");
- *stl_filename = malloc(strlen(tk) + 1);
- strcpy(*stl_filename,tk);
-
- tk = strtok(NULL," ");
- res = cstr_to_double(tk, &bound->emissivity);
- if (res != RES_OK || bound->emissivity < 0 || bound->emissivity > 1){
- fprintf(stderr, "Invalid emissivity: %s\n", tk);
- res = RES_BAD_ARG;
- goto error;
+#define CHK_TOK(tk, x) if(((tk) = (x)) == NULL) goto error
+ CHK_TOK(tk, strtok(line, " "));
+ 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"));
+ }
+ }
+ }
+ }
}
- tk = strtok(NULL," ");
- res = cstr_to_double(tk, &bound->specular_fraction);
- if (res != RES_OK ||
- bound->specular_fraction < 0.0 || bound->specular_fraction > 1.0){
- fprintf(stderr, "Invalid specular_fraction: %s\n", tk);
- res = RES_BAD_ARG;
- goto error;
+ if (h_solid || h_fluid) {
+ desc->type = h_solid ? DESC_BOUND_H_FOR_SOLID : DESC_BOUND_H_FOR_FLUID;
+
+ /* The description is parsed the same way for both types */
+ CHK_TOK(tk, strtok(NULL, " "));
+ *stl_filename = malloc(strlen(tk) + 1);
+ strcpy(*stl_filename, tk);
+
+ CHK_TOK(tk, strtok(NULL, " "));
+ 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 error;
+ }
+ desc->d.h_boundary.has_emissivity = (desc->d.h_boundary.emissivity > 0);
+ CHK_TOK(tk, strtok(NULL, " "));
+ 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 error;
+ }
+ CHK_TOK(tk, strtok(NULL, " "));
+ 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 error;
+ }
+ desc->d.h_boundary.has_hc = (desc->d.h_boundary.hc > 0);
+ CHK_TOK(tk, strtok(NULL, " "));
+ 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 error;
+ }
+ CHK_TOK(tk, strtok(NULL, "\""));
+ desc->d.h_boundary.T = malloc(strlen(tk) + 1);
+ strcpy(desc->d.h_boundary.T, tk);
}
-
- tk = strtok(NULL," ");
- if (strcmp(tk,"h")==0 || strcmp(tk,"H")==0 ){
- tk = strtok(NULL," ");
- res = cstr_to_double(tk, &bound->hc);
- if (res != RES_OK || bound->hc < 0) {
- fprintf(stderr, "Invalid value parameter: %s\n", tk);
+ else if (t_solid || t_fluid) {
+ desc->type = t_solid ? DESC_BOUND_T_FOR_SOLID : DESC_BOUND_T_FOR_FLUID;
+
+ /* This part of the description is parsed the same way for both types */
+ CHK_TOK(tk, strtok(NULL, " "));
+ *stl_filename = malloc(strlen(tk) + 1);
+ strcpy(*stl_filename, tk);
+
+ CHK_TOK(tk, strtok(NULL, "\""));
+ desc->d.t_boundary.T = malloc(strlen(tk) + 1);
+ strcpy(desc->d.t_boundary.T, tk);
+
+ /* Parse hc + hc_max only if fluid */
+ if (t_fluid) {
+ CHK_TOK(tk, strtok(NULL, " "));
+ 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 error;
+ }
+ desc->d.t_boundary.has_hc = (desc->d.t_boundary.hc > 0);
+ CHK_TOK(tk, strtok(NULL, " "));
+ 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 error;
+ }
+ }
+ else desc->d.t_boundary.has_hc = 0;
+ }
+ else if (f_solid) {
+ desc->type = DESC_BOUND_F_FOR_SOLID;
+
+ CHK_TOK(tk, strtok(NULL, " "));
+ *stl_filename = malloc(strlen(tk) + 1);
+ strcpy(*stl_filename, tk);
+
+ CHK_TOK(tk, strtok(NULL, "\""));
+ desc->d.f_boundary.flux = malloc(strlen(tk) + 1);
+ strcpy(desc->d.f_boundary.flux, tk);
+
+ /* Optional hc + hc_max */
+ tk = strtok(NULL, " ");
+ if (tk) {
+ res = cstr_to_double(tk, &desc->d.f_boundary.hc);
+ if (res != RES_OK || desc->d.f_boundary.hc < 0) {
+ fprintf(stderr, "Invalid hc value: %s\n", tk);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ desc->d.f_boundary.has_hc = (desc->d.f_boundary.hc > 0);
+ CHK_TOK(tk, strtok(NULL, " "));
+ res = cstr_to_double(tk, &desc->d.f_boundary.hc_max);
+ if (res != RES_OK
+ || desc->d.f_boundary.hc_max < desc->d.f_boundary.hc
+ || desc->d.f_boundary.hc_max < 0) {
+ fprintf(stderr, "Invalid hc_max value: %s\n", tk);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ }
+ }
+ else if (sf_connect) {
+ desc->type = DESC_SOLID_FLUID_CONNECT;
+
+ CHK_TOK(tk, strtok(NULL, " "));
+ *stl_filename = malloc(strlen(tk) + 1);
+ strcpy(*stl_filename, tk);
+
+ CHK_TOK(tk, strtok(NULL, " "));
+ 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 error;
}
-
- tk = strtok(NULL,"\"");
- tk = strtok(NULL,"\"");
- bound->T = malloc(strlen(tk) + 1);
- strcpy(bound->T,tk);
-
- } else if (strcmp(tk,"t")==0 || strcmp(tk,"T")==0 ){
-
- tk = strtok(NULL,"\"");
- tk = strtok(NULL,"\"");
- bound->T = malloc(strlen(tk) + 1);
- strcpy(bound->T,tk);
-
- } else if (strcmp(tk,"f")==0 || strcmp(tk,"F")==0 ){
-
- tk = strtok(NULL,"\"");
- tk = strtok(NULL,"\"");
- bound->flux = malloc(strlen(tk) + 1);
- strcpy(bound->flux,tk);
-
- } else {
- fprintf(stderr,"unknown boundary type: %s\n", tk);
+ desc->d.sf_connect.has_emissivity = (desc->d.sf_connect.emissivity > 0);
+ CHK_TOK(tk, strtok(NULL, " "));
+ 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 error;
+ }
+ CHK_TOK(tk, strtok(NULL, " "));
+ 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 error;
+ }
+ desc->d.sf_connect.has_hc = (desc->d.sf_connect.hc > 0);
+ CHK_TOK(tk, strtok(NULL, " "));
+ 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 error;
+ }
+ }
+ else {
+ printf("Invalid boundary type: %s", tk);
+ res = RES_BAD_ARG;
goto error;
}
+#undef GET_TOK
exit:
return res;
error:
goto exit;
}
+
static res_T
read_vertices
(struct sstl_desc* desc,
@@ -248,132 +444,198 @@ read_vertices
return res;
}
+static int
+indices_to_key(unsigned3 key, const unsigned3 indices)
+{
+ int i, reversed = 0;
+ FOR_EACH(i, 0, 3) key[i] = indices[i];
+ if (key[0] > key[1]) {
+ reversed = !reversed;
+ SWAP(unsigned, key[0], key[1]);
+ }
+ if (key[1] > key[2]) {
+ reversed = !reversed;
+ SWAP(unsigned, key[1], key[2]);
+ }
+ if (key[0] > key[1]) {
+ reversed = !reversed;
+ SWAP(unsigned, key[0], key[1]);
+ }
+ return reversed;
+}
+
static res_T
read_triangles
- (int FLAG,
- struct sstl_desc* desc,
- struct htable_triangle* triangle2id,
- unsigned* id2id, char** cl_type,
- const unsigned* id,
- struct triangle** tri_list)
+ (const struct sstl_desc* stl_desc,
+ struct htable_triangle* triangle2id,
+ const unsigned* id2id,
+ const int is_connection,
+ const unsigned desc_id,
+ struct triangle** tri_list)
{
res_T res = RES_OK;
unsigned tri_index = 0;
- (void)cl_type;
- for (tri_index = 0; tri_index < desc->triangles_count; ++tri_index){
+ for (tri_index = 0; tri_index < stl_desc->triangles_count; ++tri_index) {
struct triangle tri = NULL_TRIANGLE;
- unsigned *found_id = NULL;
-
- tri.sids[0] = tri.indices[0] = id2id[desc->indices[3*tri_index + 0]];
- tri.sids[1] = tri.indices[1] = id2id[desc->indices[3*tri_index + 1]];
- tri.sids[2] = tri.indices[2] = id2id[desc->indices[3*tri_index + 2]];
-
- if(tri.sids[0] > tri.sids[1]) SWAP(unsigned, tri.sids[0], tri.sids[1]);
- if(tri.sids[1] > tri.sids[2]) SWAP(unsigned, tri.sids[1], tri.sids[2]);
- if(tri.sids[0] > tri.sids[1]) SWAP(unsigned, tri.sids[0], tri.sids[1]);
- found_id = htable_triangle_find(triangle2id, &tri);
-
- if (found_id){
- if (FLAG==MEDIUM) (*tri_list)[*found_id].medium_back = *id;
- if (FLAG==BOUNDARY) (*tri_list)[*found_id].bound_id = *id;
+ const unsigned *found_id = NULL;
+ unsigned id;
+ unsigned3 key;
+ int reversed;
+ enum sdis_side current_side = SDIS_FRONT;
+ unsigned* side_desc_ptr = NULL;
+
+ tri.indices[0] = id2id[stl_desc->indices[3*tri_index + 0]];
+ tri.indices[1] = id2id[stl_desc->indices[3*tri_index + 1]];
+ tri.indices[2] = id2id[stl_desc->indices[3*tri_index + 2]];
+ reversed = indices_to_key(key, tri.indices);
+
+ /* 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 = (*tri_list) + id;
+ original_reversed = indices_to_key(key, original_tri->indices);
+ current_side = (reversed == original_reversed) ? SDIS_FRONT : SDIS_BACK;
+ ASSERT(id < sa_size(*tri_list));
} else {
size_t sz;
- unsigned size;
- if (FLAG==MEDIUM) tri.medium_front = *id;
sz = htable_triangle_size_get(triangle2id);
+ ASSERT(sz == sa_size(*tri_list));
ASSERT(sz < UINT_MAX);
- size = (unsigned)sz;
- htable_triangle_set(triangle2id, &tri, &size);
- /*printf("triangle %i %i %i\n added at index %i\n", SPLIT3(tri.indices), size);*/
+ id = (unsigned)sz;
+ /* Store tri as described by stl_desc regardless of key order */
+ htable_triangle_set(triangle2id, &key, &id);
sa_push(*tri_list, tri);
}
+
+ if (is_connection) {
+ side_desc_ptr = &(*tri_list)[id].connection_description;
+ }
+ else {
+ /* Is the triangle in tri_list reversed wrt stl_desc? */
+ side_desc_ptr = (current_side == SDIS_FRONT) ?
+ &(*tri_list)[id].front_description : &(*tri_list)[id].back_description;
+ }
+ if (*side_desc_ptr != UINT_MAX && *side_desc_ptr != desc_id) {
+ /* Already described with a different description! */
+ fprintf(stderr, "Triangle %s with 2 different descriptions\n",
+ is_connection ? "connection" : "side");
+ return RES_BAD_ARG;
+ }
+ /* Everithing is OK: store description */
+ *side_desc_ptr = desc_id;
}
return res;
}
-
static res_T
read_stl
-(char FLAG,
- unsigned id,
- char** stl_filename,
- struct htable_vertex* vertex2id,
- struct htable_triangle* triangle2id,
- struct geometry* geometry)
+ (const unsigned desc_id,
+ const int is_connection,
+ 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 desc;
+ struct sstl_desc stl_desc;
unsigned* id2id = NULL;
- SSTL(create(NULL, NULL, 0, &sstl));
+ if (!stl_filename) return RES_BAD_ARG;
- SSTL(load(sstl, *stl_filename));
- SSTL(get_desc(sstl, &desc));
+ SSTL(create(NULL, &stardis->allocator, 0, &sstl));
- read_vertices(&desc, vertex2id, &id2id, &geometry->vertex);
- read_triangles(FLAG, &desc, triangle2id, id2id, NULL, &id, &geometry->triangle);
- sa_release(id2id);
+ res = sstl_load(sstl, *stl_filename);
+ if (res != RES_OK) {
+ printf("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.vertex);
+ if (res != RES_OK) goto error;
+ res = read_triangles(&stl_desc, triangle2id, id2id, is_connection, desc_id,
+ &stardis->geometry.triangle);
+ if (res != RES_OK) goto error;
+
+end:
+ sa_release(id2id);
SSTL(ref_put(sstl));
return res;
+error:
+ goto end;
}
-
/*******************************************************************************
*
******************************************************************************/
static res_T
geometry_analyse
-(const char* medium_filename,
- const char* bc_filename,
- struct geometry* geometry,
- struct material** material,
- struct boundary** boundary)
+ (const char* medium_filename,
+ const char* bc_filename,
+ struct stardis* stardis)
{
res_T res = RES_OK;
FILE* input;
char* line = NULL;
+ unsigned sz = (unsigned)sa_size(stardis->descriptions);
struct htable_vertex vertex2id;
struct htable_triangle triangle2id;
+ struct htable_descriptions descriptions;
+ int descriptions_table_initialised = 0;
+ unsigned desc_id;
+ unsigned *p_desc;
+ ASSERT(sz == sa_size(stardis->descriptions));
- /*init geometry*/
- *geometry = NULL_GEOMETRY;
+ stardis->geometry = NULL_GEOMETRY;
- /*parse medium file*/
+ /* parse medium files */
input = fopen(medium_filename,"r");
if (!input){
- fprintf(stderr,"Can not open %s\n", medium_filename);
+ fprintf(stderr,"Cannot open %s\n", medium_filename);
res = RES_IO_ERR;
goto error;
}
- htable_vertex_init(NULL, &vertex2id);
- htable_triangle_init(NULL, &triangle2id);
+ htable_vertex_init(&stardis->allocator, &vertex2id);
+ htable_triangle_init(&stardis->allocator, &triangle2id);
+ htable_descriptions_init(&stardis->allocator, &descriptions);
+ descriptions_table_initialised = 1;
/* loop on media */
while (read_line(&line, input)){
char* stl_filename = NULL;
- struct material mat = NULL_MATERIAL;
+ struct description desc = NULL_DESCRIPTION__;
- res = parse_medium_line(line, &stl_filename, &mat);
+ res = parse_medium_line(line, &stl_filename, &desc);
if (res != RES_OK) goto error;
- sa_push(*material,mat);
- /*analyse medium stl*/
- res = read_stl(MEDIUM, geometry->medium_count, &stl_filename,
- &vertex2id, &triangle2id, geometry);
- if (res != RES_OK) goto error;
+ /* Deduplicate media: find if the same description already exists */
+ 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 */
+ htable_descriptions_set(&descriptions, &desc, &desc_id);
+ } else {
+ printf("Duplicate media description found: deduplicated.\n");
+ desc_id = *p_desc;
+ }
- geometry->medium_count += 1;
+ res = read_stl(sz-1, 0, &stl_filename, &vertex2id, &triangle2id, stardis);
if (stl_filename) free(stl_filename);
+ if (res != RES_OK) goto error;
}
fclose(input);
- /*parse boundary file*/
+ /* parse boundary files */
input = fopen(bc_filename,"r");
if (!input){
- fprintf(stderr,"Can not open %s\n", bc_filename);
+ fprintf(stderr,"Cannot open %s\n", bc_filename);
res = RES_IO_ERR;
goto error;
}
@@ -381,24 +643,33 @@ geometry_analyse
/* loop on boundaries */
while (read_line(&line, input)){
char* stl_filename = NULL;
- struct boundary bound = NULL_BOUNDARY;
-
- res = parse_boundary_line(line, &stl_filename, &bound);
- if (res != RES_OK) goto error;
- sa_push(*boundary,bound);
+ struct description desc = NULL_DESCRIPTION__;
- /*analyse medium stl*/
- res = read_stl(BOUNDARY, geometry->boundary_count, &stl_filename,
- &vertex2id, &triangle2id, geometry);
+ res = parse_boundary_line(line, &stl_filename, &desc);
if (res != RES_OK) goto error;
-
- geometry->boundary_count += 1;
+
+ /* Deduplicate media: find if the same description already exists */
+ 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 */
+ htable_descriptions_set(&descriptions, &desc, &desc_id);
+ } else {
+ printf("Duplicate boundary description found: deduplicated.\n");
+ desc_id = *p_desc;
+ }
+
+ res = read_stl(desc_id, 1, &stl_filename, &vertex2id, &triangle2id, stardis);
if (stl_filename) free(stl_filename);
+ if (res != RES_OK) goto error;
}
-
fclose(input);
exit:
+ if(descriptions_table_initialised)
+ htable_descriptions_release(&descriptions);
htable_vertex_release(&vertex2id);
htable_triangle_release(&triangle2id);
sa_release(line);
@@ -408,27 +679,6 @@ error:
}
static res_T
-check_consistency
-(struct boundary* boundary, double* probe)
-{
- res_T res = RES_OK;
- size_t i = 0;
-
- for (i=0; i<sa_size(boundary); ++i){
- if ( boundary[i].T && strchr(boundary[i].T,'t') && probe[3] == INF) {
- fprintf(stderr,"Transient boundary condition is inconsistent with INF time of the probe.\n");
- res = RES_BAD_ARG;
- goto error;
- }
- }
-
-exit:
- return res;
-error:
- goto exit;
-}
-
-static res_T
parse_camera
(char* cam_param, struct camera* cam)
{
@@ -487,11 +737,12 @@ stardis_init
{
res_T res = RES_OK;
+ res = mem_init_proxy_allocator(&stardis->allocator,
+ &mem_default_allocator);
+ if (res != RES_OK) goto error;
+
res = geometry_analyse(args->medium_filename,
- args->bc_filename,
- &stardis->geometry,
- &stardis->material,
- &stardis->boundary);
+ args->bc_filename, stardis);
if (res != RES_OK) goto error;
stardis->N = args->N;
@@ -504,59 +755,59 @@ stardis_init
stardis->radiative_temp[0] = args->radiative_temp[0];
stardis->radiative_temp[1] = args->radiative_temp[1];
- if (args->mode == PROBE_COMPUTE){
- res = check_consistency(stardis->boundary, stardis->probe);
- if (res != RES_OK) goto error;
- }
-
if (args->mode == IR_COMPUTE){
res = parse_camera(args->camera, &stardis->camera);
if (res != RES_OK) goto error;
}
+
exit:
return res;
error:
goto exit;
}
-
-res_T
+void
stardis_release(struct stardis* stardis)
{
size_t i = 0;
- res_T res = RES_OK;
-
- for (i=0; i<sa_size(stardis->material); ++i){
- free(stardis->material[i].Tinit);
- free(stardis->material[i].power);
- }
- sa_release(stardis->material);
- for (i=0; i<sa_size(stardis->boundary); ++i){
- free(stardis->boundary[i].T);
- free(stardis->boundary[i].flux);
+ for (i=0; i<sa_size(stardis->descriptions); ++i) {
+ struct description* desc = &stardis->descriptions[i];
+ switch (desc->type) {
+ case DESC_MAT_SOLID:
+ free(desc->d.solid.power);
+ free(desc->d.solid.Tinit);
+ break;
+ case DESC_MAT_FLUID:
+ free(desc->d.fluid.Tinit);
+ break;
+ case DESC_BOUND_H_FOR_SOLID:
+ case DESC_BOUND_H_FOR_FLUID:
+ free(desc->d.h_boundary.T);
+ break;
+ case DESC_BOUND_T_FOR_SOLID:
+ case DESC_BOUND_T_FOR_FLUID:
+ free(desc->d.t_boundary.T);
+ break;
+ case DESC_BOUND_F_FOR_SOLID:
+ free(desc->d.f_boundary.flux);
+ break;
+ default: FATAL("Invalid type.\n");
+ }
}
- sa_release(stardis->boundary);
-
- res = release_geometry(&stardis->geometry);
- if (res != RES_OK) goto error;
-
-exit:
- return res;
-error:
- goto exit;
+ sa_release(stardis->descriptions);
+ release_geometry(&stardis->geometry);
}
-
+/* TODO: rewrite dump! */
res_T
dump_vtk
-(FILE* output,
- const struct geometry* geometry,
- unsigned nbound)
+ (FILE* output,
+ const struct geometry* geometry)
{
res_T res = RES_OK;
- unsigned i,j = 0;
+ unsigned i;
struct vertex* vtx = geometry->vertex;
struct triangle* tri = geometry->triangle;
@@ -576,22 +827,15 @@ dump_vtk
fprintf(output,"3 %i %i %i\n",SPLIT3(tri[i].indices));
}
fprintf(output,"\nCELL_DATA %i \n", (int)sa_size(tri));
- fprintf(output,"SCALARS medium_front float 1\n");
+ fprintf(output,"SCALARS front_description float 1\n");
fprintf(output,"LOOKUP_TABLE default\n");
- for (i=0; i<sa_size(tri); ++i){
- fprintf(output,"%i\n",tri[i].medium_front);
+ for (i=0; i<sa_size(tri); ++i) {
+ fprintf(output,"%i\n",tri[i].front_description);
}
- fprintf(output,"SCALARS medium_back float 1\n");
+ fprintf(output,"SCALARS back_description float 1\n");
fprintf(output,"LOOKUP_TABLE default\n");
- for (i=0; i<sa_size(tri); ++i){
- fprintf(output,"%i\n",tri[i].medium_back);
- }
- for (j=0; j<nbound; ++j){
- fprintf(output,"SCALARS bound_id float 1\n");
- fprintf(output,"LOOKUP_TABLE default\n");
- for (i=0; i<sa_size(tri); ++i){
- fprintf(output,"%i\n",tri[i].bound_id);
- }
+ for (i=0; i<sa_size(tri); ++i) {
+ fprintf(output,"%i\n",tri[i].back_description);
}
exit:
diff --git a/src/stardis-app.h b/src/stardis-app.h
@@ -20,6 +20,20 @@
free(ARRAY);\
}
+/* Different types of descriptions */
+enum description_type {
+ DESC_MAT_SOLID,
+ DESC_MAT_FLUID,
+ DESC_BOUND_H_FOR_FLUID,
+ DESC_BOUND_H_FOR_SOLID,
+ DESC_BOUND_T_FOR_FLUID,
+ DESC_BOUND_T_FOR_SOLID,
+ DESC_BOUND_F_FOR_SOLID,
+ DESC_SOLID_FLUID_CONNECT,
+ DESCRIPTION_TYPE_COUNT__,
+ DESC_OUTSIDE
+};
+
struct vertex{
float xyz[3];
};
@@ -39,83 +53,709 @@ eq_vertex(const struct vertex* a, const struct vertex* b)
#define HTABLE_KEY_FUNCTOR_EQ eq_vertex
#include <rsys/hash_table.h>
-struct triangle{
- unsigned indices[3];
- unsigned sids[3];
- unsigned medium_front;
- unsigned medium_back;
- unsigned bound_id;
+typedef unsigned unsigned3[3];
+struct triangle {
+ unsigned3 indices;
+ unsigned front_description;
+ unsigned back_description;
+ unsigned connection_description;
};
-#define NULL_TRIANGLE__ {{0, 0, 0}, {0, 0, 0}, UINT_MAX, UINT_MAX, UINT_MAX}
+#define NULL_TRIANGLE__ {{0, 0, 0}, UINT_MAX, UINT_MAX, UINT_MAX }
static const struct triangle NULL_TRIANGLE = NULL_TRIANGLE__;
static char
-eq_triangle(const struct triangle* a, const struct triangle* b)
+eq_indices(const unsigned3* a, const unsigned3* b)
{
- return a->sids[0] == b->sids[0]
- && a->sids[1] == b->sids[1]
- && a->sids[2] == b->sids[2];
+ return (*a)[0] == (*b)[0] && (*a)[1] == (*b)[1] && (*a)[2] == (*b)[2];
}
-static size_t
-hash_triangle(const struct triangle* tri)
+static res_T
+cp_indices(unsigned3* dst, const unsigned3* src)
{
- return hash_fnv64(tri->sids, sizeof(unsigned[3]));
+ int i;
+ FOR_EACH(i, 0, 3) (*dst)[i] = (*src)[i];
+ return RES_OK;
}
-/* Declare the hash table that map a triangle to its index */
+/* Declare the hash table that maps a triangle id to its index */
#define HTABLE_NAME triangle
#define HTABLE_DATA unsigned
-#define HTABLE_KEY struct triangle
-#define HTABLE_KEY_FUNCTOR_EQ eq_triangle
-#define HTABLE_KEY_FUNCTOR_HASH hash_triangle
+#define HTABLE_KEY unsigned3
+#define HTABLE_KEY_FUNCTOR_EQ eq_indices
+#define HTABLE_KEY_FUNCTOR_COPY cp_indices
#include <rsys/hash_table.h>
-struct geometry{
+struct geometry {
struct vertex* vertex;
struct triangle* triangle;
- unsigned medium_count;
- unsigned boundary_count;
+ struct sdis_interface** interf_bytrg;
struct sdis_interface** interfaces;
};
-#define NULL_GEOMETRY__ {NULL, NULL, 0, 0, NULL}
+#define NULL_GEOMETRY__ {NULL, NULL, NULL, NULL }
static const struct geometry NULL_GEOMETRY = NULL_GEOMETRY__;
-static inline res_T
-release_geometry
-(struct geometry* geom)
+static inline void
+release_geometry(struct geometry* geom)
{
- res_T res = RES_OK;
+ 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 mat_fluid {
+ unsigned fluid_id;
+ double rho;
+ double cp;
+ char* Tinit;
+};
- sa_release(geom->vertex);
- sa_release(geom->triangle);
- sa_release(geom->interfaces);
+static char
+eq_fluid(const struct mat_fluid* a, const struct mat_fluid* b)
+{
+ if (a->fluid_id != b->fluid_id
+ || a->rho != b->rho
+ || a->cp != b->cp
+ || strcmp(a->Tinit, b->Tinit))
+ return 0;
+ return 1;
+}
- return res;
+static size_t
+hash_fluid(const struct mat_fluid* key)
+{
+#ifdef ARCH_32BITS
+ /* 32-bits Fowler/Noll/Vo hash function */
+ const uint32_t FNV32_PRIME =
+ (uint32_t) (((uint32_t) 1 << 24) + ((uint32_t) 1 << 8) + 0x93);
+ const uint32_t OFFSET32_BASIS = 2166136261u;
+ uint32_t hash = OFFSET32_BASIS;
+ size_t i;
+ ASSERT(key);
+ FOR_EACH(i, 0, sizeof(key->fluid_id)) {
+ hash = hash ^ (uint32_t)((unsigned char)((const char*)&key->fluid_id)[i]);
+ hash = hash * FNV32_PRIME;
+ }
+ FOR_EACH(i, 0, sizeof(key->rho)) {
+ hash = hash ^ (uint32_t)((unsigned char)((const char*)&key->rho)[i]);
+ hash = hash * FNV32_PRIME;
+ }
+ FOR_EACH(i, 0, sizeof(key->cp)) {
+ hash = hash ^ (uint32_t)((unsigned char)((const char*)&key->cp)[i]);
+ hash = hash * FNV32_PRIME;
+ }
+ FOR_EACH(i, 0, strlen(key->Tinit) * sizeof(*key->Tinit)) {
+ hash = hash ^ (uint32_t)((unsigned char)((const char*)key->Tinit)[i]);
+ hash = hash * FNV32_PRIME;
+ }
+ return hash;
+#elif defined(ARCH_64BITS)
+ /* 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->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;
+ }
+ FOR_EACH(i, 0, strlen(key->Tinit)*sizeof(*key->Tinit)) {
+ hash = hash ^ (uint64_t)((unsigned char) ((const char*)key->Tinit)[i]);
+ hash = hash * FNV64_PRIME;
+ }
+ return hash;
+#else
+#error "Unexpected architecture"
+#endif
}
-struct material{
+struct mat_solid {
+ unsigned solid_id;
double lambda;
double rho;
double cp;
double delta;
char* Tinit;
char* power;
+ int has_power;
};
-#define NULL_MATERIAL__ {0.0 ,0.0 ,0.0, 0.0, NULL, NULL}
-static const struct material NULL_MATERIAL = NULL_MATERIAL__;
-struct boundary{
+static char
+eq_solid(const struct mat_solid* a, const struct mat_solid* b)
+{
+ if (a->solid_id != b->solid_id
+ || a->lambda != b->lambda
+ || a->rho != b->rho
+ || a->cp != b->cp
+ || a->delta != b->delta
+ || strcmp(a->Tinit, b->Tinit)
+ || a->has_power != b->has_power)
+ return 0;
+ if (a->has_power && a->power != b->power) return 0;
+ return 1;
+}
+
+static size_t
+hash_solid(const struct mat_solid* key)
+{
+#ifdef ARCH_32BITS
+ /* 32-bits Fowler/Noll/Vo hash function */
+ const uint32_t FNV32_PRIME =
+ (uint32_t) (((uint32_t) 1 << 24) + ((uint32_t) 1 << 8) + 0x93);
+ const uint32_t OFFSET32_BASIS = 2166136261u;
+ uint32_t hash = OFFSET32_BASIS;
+ size_t i;
+ ASSERT(key);
+ FOR_EACH(i, 0, sizeof(key->solid_id)) {
+ hash = hash ^ (uint32_t)((unsigned char)((const char*)&key->solid_id)[i]);
+ hash = hash * FNV32_PRIME;
+ }
+ FOR_EACH(i, 0, sizeof(key->lambda)) {
+ hash = hash ^ (uint32_t)((unsigned char)((const char*)&key->lambda)[i]);
+ hash = hash * FNV32_PRIME;
+ }
+ FOR_EACH(i, 0, sizeof(key->rho)) {
+ hash = hash ^ (uint32_t)((unsigned char)((const char*)&key->rho)[i]);
+ hash = hash * FNV32_PRIME;
+ }
+ FOR_EACH(i, 0, sizeof(key->cp)) {
+ hash = hash ^ (uint32_t)((unsigned char)((const char*)&key->cp)[i]);
+ hash = hash * FNV32_PRIME;
+ }
+ FOR_EACH(i, 0, sizeof(key->delta)) {
+ hash = hash ^ (uint32_t)((unsigned char)((const char*)&key->delta)[i]);
+ hash = hash * FNV32_PRIME;
+}
+ FOR_EACH(i, 0, strlen(key->Tinit) * sizeof(*key->Tinit)) {
+ hash = hash ^ (uint32_t)((unsigned char)((const char*)key->Tinit)[i]);
+ hash = hash * FNV32_PRIME;
+ }
+ FOR_EACH(i, 0, sizeof(key->has_power)) {
+ hash = hash ^ (uint32_t)((unsigned char)((const char*)&key->has_power)[i]);
+ hash = hash * FNV32_PRIME;
+ }
+ if (!key->has_power) return hash;
+ FOR_EACH(i, 0, strlen(key->Tinit) * sizeof(*key->power)) {
+ hash = hash ^ (uint32_t) ((unsigned char) ((const char*) key->power)[i]);
+ hash = hash * FNV32_PRIME;
+ }
+ return hash;
+#elif defined(ARCH_64BITS)
+ /* 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->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;
+ }
+ FOR_EACH(i, 0, strlen(key->Tinit) * sizeof(*key->Tinit)) {
+ hash = hash ^ (uint64_t)((unsigned char)((const char*)key->Tinit)[i]);
+ hash = hash * FNV64_PRIME;
+ }
+ FOR_EACH(i, 0, sizeof(key->has_power)) {
+ hash = hash ^ (uint64_t) ((unsigned char) ((const char*) &key->has_power)[i]);
+ hash = hash * FNV64_PRIME;
+ }
+ if (!key->has_power) return hash;
+ FOR_EACH(i, 0, strlen(key->power) * sizeof(*key->power)) {
+ hash = hash ^ (uint64_t)((unsigned char)((const char*)key->power)[i]);
+ hash = hash * FNV64_PRIME;
+ }
+ return hash;
+#else
+#error "Unexpected architecture"
+#endif
+}
+
+struct h_boundary {
+ unsigned mat_id;
+ double emissivity;
+ double specular_fraction;
double hc;
+ double hc_max;
char* T;
+ int has_emissivity, has_hc;
+};
+
+static char
+eq_h_boundary
+ (const struct h_boundary* a,
+ const struct h_boundary* b)
+{
+ if (a->mat_id != b->mat_id
+ || a->specular_fraction != b->specular_fraction
+ || strcmp(a->T, b->T)
+ || 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 size_t
+hash_h_boundary(const struct h_boundary* key)
+{
+#ifdef ARCH_32BITS
+ /* 32-bits Fowler/Noll/Vo hash function */
+ const uint32_t FNV32_PRIME =
+ (uint32_t) (((uint32_t) 1 << 24) + ((uint32_t) 1 << 8) + 0x93);
+ const uint32_t OFFSET32_BASIS = 2166136261u;
+ uint32_t hash = OFFSET32_BASIS;
+ size_t i;
+ ASSERT(key);
+ FOR_EACH(i, 0, sizeof(key->mat_id)) {
+ hash = hash ^ (uint32_t)((unsigned char)((const char*)&key->mat_id)[i]);
+ hash = hash * FNV32_PRIME;
+ }
+ FOR_EACH(i, 0, sizeof(key->emissivity)) {
+ hash = hash ^ (uint32_t)((unsigned char)((const char*)&key->emissivity)[i]);
+ hash = hash * FNV32_PRIME;
+ }
+ FOR_EACH(i, 0, sizeof(key->specular_fraction)) {
+ hash = hash ^ (uint32_t)((unsigned char)((const char*)&key->specular_fraction)[i]);
+ hash = hash * FNV32_PRIME;
+ }
+ FOR_EACH(i, 0, strlen(key->T) * sizeof(*key->T)) {
+ hash = hash ^ (uint32_t)((unsigned char)((const char*)key->T)[i]);
+ hash = hash * FNV32_PRIME;
+ }
+ FOR_EACH(i, 0, sizeof(key->has_emissivity)) {
+ hash = hash ^ (uint32_t)((unsigned char)((const char*)&key->has_emissivity)[i]);
+ hash = hash * FNV32_PRIME;
+ }
+ FOR_EACH(i, 0, sizeof(key->has_hc)) {
+ hash = hash ^ (uint32_t) ((unsigned char) ((const char*) &key->has_hc)[i]);
+ hash = hash * FNV32_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 * FNV32_PRIME;
+ }
+ FOR_EACH(i, 0, sizeof(key->hc_max)) {
+ hash = hash ^ (uint32_t) ((unsigned char) ((const char*) &key->hc_max)[i]);
+ hash = hash * FNV32_PRIME;
+ }
+ return hash;
+#elif defined(ARCH_64BITS)
+ /* 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->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, strlen(key->T) * sizeof(*key->T)) {
+ hash = hash ^ (uint32_t)((unsigned char)((const char*)key->T)[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;
+#else
+#error "Unexpected architecture"
+#endif
+}
+
+struct t_boundary {
+ unsigned mat_id;
+ char* T;
+ struct te_expr* te_temperature;
+ double hc;
+ double hc_max;
+ int has_hc;
+};
+
+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 (a->mat_id != b->mat_id || strcmp(a->T, b->T))
+ 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;
+}
+
+static size_t
+hash_t_boundary
+ (const struct t_boundary* key,
+ const enum description_type type)
+{
+ ASSERT(type == DESC_BOUND_T_FOR_SOLID || type == DESC_BOUND_T_FOR_FLUID);
+#ifdef ARCH_32BITS
+ /* 32-bits Fowler/Noll/Vo hash function */
+ const uint32_t FNV32_PRIME =
+ (uint32_t) (((uint32_t) 1 << 24) + ((uint32_t) 1 << 8) + 0x93);
+ const uint32_t OFFSET32_BASIS = 2166136261u;
+ uint32_t hash = OFFSET32_BASIS;
+ size_t i;
+ ASSERT(key);
+ FOR_EACH(i, 0, sizeof(key->mat_id)) {
+ hash = hash ^ (uint32_t)((unsigned char)((const char*)&key->mat_id)[i]);
+ hash = hash * FNV32_PRIME;
+ }
+ FOR_EACH(i, 0, sizeof(key->has_hc)) {
+ hash = hash ^ (uint32_t) ((unsigned char) ((const char*) &key->has_hc)[i]);
+ hash = hash * FNV32_PRIME;
+ }
+ FOR_EACH(i, 0, strlen(key->T) * sizeof(*key->T)) {
+ hash = hash ^ (uint32_t) ((unsigned char) ((const char*) key->T)[i]);
+ hash = hash * FNV32_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 * FNV32_PRIME;
+ }
+ FOR_EACH(i, 0, sizeof(key->hc_max)) {
+ hash = hash ^ (uint32_t)((unsigned char)((const char*)&key->hc_max)[i]);
+ hash = hash * FNV32_PRIME;
+ }
+ return hash;
+#elif defined(ARCH_64BITS)
+ /* 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->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, strlen(key->T) * sizeof(*key->T)) {
+ hash = hash ^ (uint32_t) ((unsigned char) ((const char*) key->T)[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;
+#else
+#error "Unexpected architecture"
+#endif
+}
+
+struct f_boundary {
+ unsigned mat_id;
char* flux;
+ struct te_expr* te_flux;
+ double hc;
+ double hc_max;
+ int has_hc;
+};
+
+static char
+eq_f_boundary(const struct f_boundary* a, const struct f_boundary* b)
+{
+ if (a->mat_id != b->mat_id
+ || strcmp(a->flux, b->flux)
+ || 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;
+}
+
+static size_t
+hash_f_boundary(const struct f_boundary* key)
+{
+#ifdef ARCH_32BITS
+ /* 32-bits Fowler/Noll/Vo hash function */
+ const uint32_t FNV32_PRIME =
+ (uint32_t) (((uint32_t) 1 << 24) + ((uint32_t) 1 << 8) + 0x93);
+ const uint32_t OFFSET32_BASIS = 2166136261u;
+ uint32_t hash = OFFSET32_BASIS;
+ size_t i;
+ ASSERT(key);
+ FOR_EACH(i, 0, sizeof(key->mat_id)) {
+ hash = hash ^ (uint32_t)((unsigned char)((const char*)&key->mat_id)[i]);
+ hash = hash * FNV32_PRIME;
+ }
+ FOR_EACH(i, 0, sizeof(key->hc)) {
+ hash = hash ^ (uint32_t)((unsigned char)((const char*)&key->hc)[i]);
+ hash = hash * FNV32_PRIME;
+ }
+ FOR_EACH(i, 0, sizeof(key->hc_max)) {
+ hash = hash ^ (uint32_t)((unsigned char)((const char*)&key->hc_max)[i]);
+ hash = hash * FNV32_PRIME;
+ }
+ FOR_EACH(i, 0, strlen(key->flux) * sizeof(*key->flux)) {
+ hash = hash ^ (uint32_t)((unsigned char)((const char*)key->flux)[i]);
+ hash = hash * FNV32_PRIME;
+ }
+ FOR_EACH(i, 0, sizeof(key->has_hc)) {
+ hash = hash ^ (uint32_t)((unsigned char)((const char*)&key->has_hc)[i]);
+ hash = hash * FNV32_PRIME;
+ }
+ return hash;
+#elif defined(ARCH_64BITS)
+ /* 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->mat_id)) {
+ hash = hash ^ (uint32_t)((unsigned char)((const char*)&key->mat_id)[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, strlen(key->flux) * sizeof(*key->flux)) {
+ hash = hash ^ (uint32_t)((unsigned char)((const char*)key->flux)[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;
+#else
+#error "Unexpected architecture"
+#endif
+}
+
+struct solid_fluid_connect {
double emissivity;
double specular_fraction;
+ double hc;
+ double hc_max;
+ int has_emissivity, has_hc;
};
-#define NULL_BOUNDARY__ {-1, NULL, NULL, 0.0 , 0.0}
-static const struct boundary NULL_BOUNDARY = NULL_BOUNDARY__;
-struct camera{
+static char
+eq_sf_connect(const struct solid_fluid_connect* a, const struct solid_fluid_connect* b)
+{
+ return (char) (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 size_t
+hash_sf_connect(const struct solid_fluid_connect* key)
+{
+#ifdef ARCH_32BITS
+ /* 32-bits Fowler/Noll/Vo hash function */
+ const uint32_t FNV32_PRIME =
+ (uint32_t) (((uint32_t) 1 << 24) + ((uint32_t) 1 << 8) + 0x93);
+ const uint32_t OFFSET32_BASIS = 2166136261u;
+ uint32_t hash = OFFSET32_BASIS;
+ size_t i;
+ ASSERT(key);
+ FOR_EACH(i, 0, sizeof(key->emissivity)) {
+ hash = hash ^ (uint32_t)((unsigned char)((const char*)&key->emissivity)[i]);
+ hash = hash * FNV32_PRIME;
+ }
+ FOR_EACH(i, 0, sizeof(key->specular_fraction)) {
+ hash = hash ^ (uint32_t)((unsigned char)((const char*)&key->specular_fraction)[i]);
+ hash = hash * FNV32_PRIME;
+ }
+ FOR_EACH(i, 0, sizeof(key->hc)) {
+ hash = hash ^ (uint32_t)((unsigned char)((const char*)&key->hc)[i]);
+ hash = hash * FNV32_PRIME;
+ }
+ FOR_EACH(i, 0, sizeof(key->hc_max)) {
+ hash = hash ^ (uint32_t)((unsigned char)((const char*)&key->hc_max)[i]);
+ hash = hash * FNV32_PRIME;
+ }
+ FOR_EACH(i, 0, sizeof(key->has_emissivity)) {
+ hash = hash ^ (uint32_t)((unsigned char)((const char*)&key->has_emissivity)[i]);
+ hash = hash * FNV32_PRIME;
+ }
+ FOR_EACH(i, 0, sizeof(key->has_hc)) {
+ hash = hash ^ (uint32_t)((unsigned char)((const char*)&key->has_hc)[i]);
+ hash = hash * FNV32_PRIME;
+ }
+ return hash;
+#elif defined(ARCH_64BITS)
+ /* 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->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;
+#else
+#error "Unexpected architecture"
+#endif
+}
+
+struct description {
+ enum description_type type;
+ union d {
+ struct mat_fluid fluid;
+ struct mat_solid solid;
+ struct t_boundary t_boundary;
+ struct f_boundary f_boundary;
+ struct h_boundary h_boundary;
+ struct solid_fluid_connect sf_connect;
+ } d;
+};
+#define NULL_DESCRIPTION__ {DESCRIPTION_TYPE_COUNT__ };
+
+static char
+eq_description(const struct description* a, const struct description* b)
+{
+ if (a->type != b->type) return 0;
+ switch (a->type) {
+ case DESC_MAT_SOLID:
+ return eq_solid(&a->d.solid, &b->d.solid);
+ 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);
+ 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("Invalid type.\n");
+ }
+}
+
+static 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);
+ 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);
+ case DESC_BOUND_F_FOR_SOLID:
+ return hash_f_boundary(&key->d.f_boundary);
+ case DESC_SOLID_FLUID_CONNECT:
+ return hash_sf_connect(&key->d.sf_connect);
+ default: FATAL("Invalid type.\n");
+ }
+}
+
+struct camera {
double pos[3];
double tgt[3];
double up[3];
@@ -126,10 +766,9 @@ struct camera{
#define NULL_CAMERA__ {{1,1,1},{0,0,0},{0,0,1},30,4,{640,480}}
static const struct camera NULL_CAMERA = NULL_CAMERA__;
-struct stardis{
+struct stardis {
struct geometry geometry;
- struct material* material; /*array of materials*/
- struct boundary* boundary; /*array of boundaries*/
+ struct description* descriptions; /*array of materials and boundaries */
size_t N; /*number of MC realizations*/
unsigned nthreads;
@@ -137,8 +776,9 @@ struct stardis{
double scale_factor;
double radiative_temp[2];
struct camera camera;
+ struct mem_allocator allocator;
};
-#define NULL_STARDIS__ {NULL_GEOMETRY__, NULL, NULL, 0, 0, {0,0,0,0}, 0, {300,300}, NULL_CAMERA__}
+#define NULL_STARDIS__ {NULL_GEOMETRY__, NULL, 0, 0, {0,0,0,0}, 0, {300,300}, NULL_CAMERA__}
static const struct stardis NULL_STARDIS = NULL_STARDIS__;
extern res_T
@@ -149,10 +789,10 @@ stardis_init
extern res_T
stardis_compute(struct stardis* stardis, enum stardis_mode mode);
-extern res_T
+void
stardis_release(struct stardis* stardis);
extern res_T
-dump_vtk(FILE* output, const struct geometry* geometry, unsigned nbound);
+dump_vtk(FILE* output, const struct geometry* geometry);
#endif /*STARDIS-APP_H*/
diff --git a/src/stardis-compute.c b/src/stardis-compute.c
@@ -38,7 +38,32 @@ geometry_get_interface
void* context)
{
struct geometry* geom = context;
- *interf = geom->interfaces[itri];
+ *interf = geom->interf_bytrg[itri];
+}
+
+static res_T
+compile_expr_to_fn
+ (struct te_expr** f,
+ const char* math_expr,
+ const int t_allowed,
+ int* is_zero)
+{
+ te_variable vars[4] = {
+ TE_DEF_OFFSET("x", offsetof(struct sdis_rwalk_vertex, P[0])),
+ TE_DEF_OFFSET("y", offsetof(struct sdis_rwalk_vertex, P[1])),
+ TE_DEF_OFFSET("z", offsetof(struct sdis_rwalk_vertex, P[2])),
+ TE_DEF_OFFSET("t", offsetof(struct sdis_rwalk_vertex, time))
+ };
+
+ ASSERT(math_expr);
+ *f = te_compile(math_expr, vars, t_allowed ? 4 : 3, NULL);
+ if (!*f) {
+ if (!t_allowed && te_compile(math_expr, vars, 4, NULL))
+ printf("Expression use variable t and should not.\n");
+ return RES_BAD_ARG;
+ }
+ if (is_zero) *is_zero = !((*f)->type == TE_CONSTANT && (*f)->v.value == 0);
+ return RES_OK;
}
/*******************************************************************************
@@ -69,27 +94,22 @@ fluid_get_volumic_mass
return fluid_props->rho;
}
-res_T
-set_fluid_temperature_fn(struct fluid* fluid, const char* math_expr)
+static double
+fluid_get_temperature
+(const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
{
- te_variable vars[4] = {
- TE_OFFSET("x", offsetof(struct sdis_rwalk_vertex, P[0])),
- TE_OFFSET("y", offsetof(struct sdis_rwalk_vertex, P[1])),
- TE_OFFSET("z", offsetof(struct sdis_rwalk_vertex, P[2])),
- TE_OFFSET("t", offsetof(struct sdis_rwalk_vertex, time))
- };
-
- ASSERT(math_expr);
- fluid->temperature = te_compile(math_expr, vars, 4, NULL);
- if(!fluid->temperature) return RES_BAD_ARG;
- return RES_OK;
+ const struct fluid* fluid_props = sdis_data_cget(data);
+ return te_eval(fluid_props->temperature, vtx);
}
static double
-fluid_get_temperature
+fluid_get_tinit
(const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
{
const struct fluid* fluid_props = sdis_data_cget(data);
+ if(vtx->time > 0) {
+ return -1;
+ }
return te_eval(fluid_props->temperature, vtx);
}
@@ -161,24 +181,16 @@ solid_get_delta_boundary
}
#endif
-res_T
-set_solid_temperature_fn(struct solid* solid, const char* math_expr)
+static double
+solid_get_temperature
+(const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
{
- te_variable vars[4] = {
- TE_OFFSET("x", offsetof(struct sdis_rwalk_vertex, P[0])),
- TE_OFFSET("y", offsetof(struct sdis_rwalk_vertex, P[1])),
- TE_OFFSET("z", offsetof(struct sdis_rwalk_vertex, P[2])),
- TE_OFFSET("t", offsetof(struct sdis_rwalk_vertex, time))
- };
-
- ASSERT(math_expr);
- solid->t_init = te_compile(math_expr, vars, 4, NULL);
- if (!solid->t_init) return RES_BAD_ARG;
- return RES_OK;
+ const struct solid* solid_props = sdis_data_cget(data);
+ return te_eval(solid_props->t_init, vtx);
}
static double
-solid_get_temperature
+solid_get_tinit
(const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
{
const struct solid* solid_props = sdis_data_cget(data);
@@ -188,24 +200,6 @@ solid_get_temperature
return te_eval(solid_props->t_init, vtx);
}
-res_T
-set_solid_power_fn(struct solid* solid, const char* math_expr, int* has)
-{
- te_variable vars[4] = {
- TE_OFFSET("x", offsetof(struct sdis_rwalk_vertex, P[0])),
- TE_OFFSET("y", offsetof(struct sdis_rwalk_vertex, P[1])),
- TE_OFFSET("z", offsetof(struct sdis_rwalk_vertex, P[2])),
- TE_OFFSET("t", offsetof(struct sdis_rwalk_vertex, time))
- };
-
- ASSERT(math_expr);
- solid->power = te_compile(math_expr, vars, 4, NULL);
- if (!solid->power) return RES_BAD_ARG;
- if(has) /* has power? */
- *has = !(solid->power->type == TE_CONSTANT && solid->power->v.value == 0);
- return RES_OK;
-}
-
static double
solid_get_power
(const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
@@ -242,22 +236,6 @@ interface_get_convection_coef
return interface_props->hc;
}
-res_T
-set_interface_temperature_fn(struct intface* intface, const char* math_expr)
-{
- te_variable vars[4] = {
- TE_OFFSET("x", offsetof(struct sdis_interface_fragment, P[0])),
- TE_OFFSET("y", offsetof(struct sdis_interface_fragment, P[1])),
- TE_OFFSET("z", offsetof(struct sdis_interface_fragment, P[2])),
- TE_OFFSET("t", offsetof(struct sdis_interface_fragment, time))
- };
-
- ASSERT(math_expr);
- intface->temperature = te_compile(math_expr, vars, 4, NULL);
- if (!intface->temperature) return RES_BAD_ARG;
- return RES_OK;
-}
-
static double
interface_get_temperature
(const struct sdis_interface_fragment* frag, struct sdis_data* data)
@@ -270,22 +248,6 @@ interface_get_temperature
return te_eval(interface_props->temperature, frag);
}
-res_T
-set_interface_flux_fn(struct intface* intface, const char* math_expr)
-{
- te_variable vars[4] = {
- TE_OFFSET("x", offsetof(struct sdis_interface_fragment, P[0])),
- TE_OFFSET("y", offsetof(struct sdis_interface_fragment, P[1])),
- TE_OFFSET("z", offsetof(struct sdis_interface_fragment, P[2])),
- TE_OFFSET("t", offsetof(struct sdis_interface_fragment, time))
- };
-
- ASSERT(math_expr);
- intface->flux = te_compile(math_expr, vars, 4, NULL);
- if (!intface->flux) return RES_BAD_ARG;
- return RES_OK;
-}
-
static double
interface_get_flux
(const struct sdis_interface_fragment* frag, struct sdis_data* data)
@@ -314,7 +276,6 @@ interface_get_emissivity
return interface_props->emissivity;
}
-
static double
interface_get_alpha
(const struct sdis_interface_fragment* frag, struct sdis_data* data)
@@ -324,31 +285,37 @@ interface_get_alpha
return interface_props->alpha;
}
-
static res_T
select_probe_type
-(const struct sdis_scene* scn,
- const struct triangle* triangle,
- const struct material* mat,
- double pos[],
- size_t* iprim,
- double uv[])
+ (const struct sdis_scene* scn,
+ const struct triangle* triangle,
+ const struct description* descriptions,
+ double pos[3],
+ size_t* iprim,
+ double uv[2])
{
res_T res = RES_OK;
size_t i = 0;
- double projected_pos[3] = {0,0,0};
- double dp[3] = {0,0,0};
- double delta = 0;
- for (i=0; i<sa_size(triangle); ++i){
-
+ for (i = 0; i < sa_size(triangle); ++i) {
+ double delta = 0, dp[3], projected_pos[3];
res = sdis_scene_boundary_project_position(scn, i, pos, uv);
-
+ if (res != RES_OK) return res;
res = sdis_scene_get_boundary_position(scn, i, uv, projected_pos);
-
+ if (res != RES_OK) return res;
d3_sub(dp, pos, projected_pos);
- delta = mat[triangle[i].medium_front].delta;
+ if (triangle[i].front_description != UINT_MAX) {
+ const struct description *desc = descriptions + triangle[i].front_description;
+ if (desc->type == DESC_MAT_SOLID) delta = MMAX(desc->d.solid.delta, delta);
+
+ }
+ if (triangle[i].back_description != UINT_MAX) {
+ const struct description *desc = descriptions + triangle[i].back_description;
+ if (desc->type == DESC_MAT_SOLID) delta = MMAX(desc->d.solid.delta, delta);
+
+ }
+
if (d3_len(dp) < delta){
*iprim = i;
fprintf(stderr,"The probe is very close to the primitive %llu.\n",
@@ -400,183 +367,458 @@ dump_image(const struct sdis_accum_buffer* buf)
/*******************************************************************************
*
******************************************************************************/
-
-res_T
-stardis_compute(struct stardis* stardis, enum stardis_mode mode)
+
+struct int_descs {
+ unsigned front, back, connect; unsigned id;
+};
+#define INT_DESCS_NULL__ { UINT_MAX, UINT_MAX, UINT_MAX };
+static const struct int_descs INT_DESCS_NULL = INT_DESCS_NULL__;
+
+static inline char
+eq_desc(const struct int_descs* a, const struct int_descs* b)
{
- res_T res = RES_OK;
- struct sdis_device* dev = NULL;
- struct sdis_data* data = NULL;
- struct solid* solid_props = NULL;
- struct fluid* fluid_props = NULL;
- struct sdis_solid_shader solid_shader = SDIS_SOLID_SHADER_NULL;
- struct sdis_medium** solid_medium = NULL;
- struct sdis_fluid_shader fluid_shader = SDIS_FLUID_SHADER_NULL;
- struct sdis_medium** fluid_medium = NULL;
-
- struct sdis_interface_side_shader interface_shader_front =
- SDIS_INTERFACE_SIDE_SHADER_NULL;
- struct sdis_interface_side_shader interface_shader_back =
- SDIS_INTERFACE_SIDE_SHADER_NULL;
- struct sdis_interface_shader interface_shader_boundary =
- SDIS_INTERFACE_SHADER_NULL;
- struct sdis_interface_shader interface_shader_connection =
- SDIS_INTERFACE_SHADER_NULL;
- struct intface* interface_props = NULL;
- struct sdis_interface** interfaces = NULL;
+ return (char)(a->front == b->front && a->back == b->back
+ && a->connect == b->connect);
+}
- struct sdis_scene* scn = NULL;
- struct sdis_estimator* estimator = NULL;
- struct sdis_mc temperature;
- struct sdis_accum_buffer* buf = NULL;
- struct sdis_camera* cam = NULL;
+static INLINE size_t
+hash_desc(struct int_descs const* key)
+{
+#ifdef ARCH_32BITS
+ return (size_t)hash_fnv32(key, 3 * sizeof(unsigned));
+#elif defined(ARCH_64BITS)
+ return (size_t)hash_fnv64(key, 3 * sizeof(unsigned));
+#else
+#error "Unexpected architecture"
+#endif
+}
- int* bound2fluid = NULL;
- double pos[3] = {0,0,0};
- double time[2] = { 0, 0};
- size_t nfailures;
- unsigned i = 0;
+#include <rsys/hash_table.h>
+/* Declare the hash table that map a vertex to its index */
+#define HTABLE_NAME intface
+#define HTABLE_DATA struct sdis_interface*
+#define HTABLE_KEY struct int_descs
+#define HTABLE_KEY_FUNCTOR_EQ eq_desc
+#define HTABLE_KEY_FUNCTOR_HASH hash_desc
+#include <rsys/hash_table.h>
- SDIS(device_create(NULL, NULL, stardis->nthreads, 1, &dev));
+static res_T
+create_fluid
+ (struct sdis_device* dev,
+ const double rho,
+ const double cp,
+ const int t_allowed,
+ const char* t_expr,
+ struct sdis_medium*** media_ptr,
+ sdis_medium_getter_T get_temp,
+ unsigned* out_id)
+{
+ struct sdis_fluid_shader fluid_shader = SDIS_FLUID_SHADER_NULL;
+ struct sdis_data* data = NULL;
+ struct fluid* fluid_props;
+ size_t sz;
+ res_T res = RES_OK;
- /* Setup the fluid shader */
+ ASSERT(dev && rho >= 0 && cp >= 0 && t_expr && media_ptr && out_id);
fluid_shader.calorific_capacity = fluid_get_calorific_capacity;
fluid_shader.volumic_mass = fluid_get_volumic_mass;
- fluid_shader.temperature = fluid_get_temperature;
-
- /* create fluid for each hc boundary */
- for (i=0; i<stardis->geometry.boundary_count; ++i){
- SDIS(data_create(dev, sizeof(struct fluid), ALIGNOF(struct fluid),
- release_fluid_data, &data));
- fluid_props = sdis_data_get(data); /* Fetch the allocated memory space */
- fluid_props->cp = 0;
- fluid_props->rho = 0;
- res = set_fluid_temperature_fn(fluid_props, stardis->boundary[i].T);
- if (res != RES_OK) {
- fprintf(stderr, "Invalid temperature expression: %s\n", stardis->boundary[i].T);
- goto error;
- }
- SDIS(fluid_create(dev, &fluid_shader, data, sa_add(fluid_medium, 1)));
- SDIS(data_ref_put(data)); data = NULL;
- if (stardis->boundary[i].hc > -1) {
- size_t sz = sa_size(fluid_medium) - 1;
- ASSERT(sz < INT_MAX);
- sa_push(bound2fluid, (int)sz);
- } else {
- sa_push(bound2fluid, 0);
- }
+ fluid_shader.temperature = get_temp;
+ res = sdis_data_create(dev, sizeof(struct fluid), ALIGNOF(struct fluid),
+ release_fluid_data, &data);
+ if (res != RES_OK) goto error;
+ fluid_props = sdis_data_get(data); /* Fetch the allocated memory space */
+ fluid_props->cp = cp;
+ fluid_props->rho = rho;
+ res = compile_expr_to_fn(&fluid_props->temperature, t_expr, t_allowed, NULL);
+ if (res != RES_OK) {
+ fprintf(stderr, "Invalid initial temperature expression: %s\n",
+ t_expr);
+ goto error;
}
+ res = sdis_fluid_create(dev, &fluid_shader, data, sa_add(*media_ptr, 1));
+ if (res != RES_OK) goto error;
+ sz = sa_size(*media_ptr) - 1;
+ ASSERT(sz < INT_MAX);
+ *out_id = (unsigned)sz;
- /* Setup the solid shader */
+end:
+ if (data) SDIS(data_ref_put(data));
+ return res;
+error:
+ goto end;
+}
+
+static res_T
+create_solid
+ (struct sdis_device* dev,
+ const double lambda,
+ const double rho,
+ const double cp,
+ const double delta,
+ const int t_allowed,
+ const char* t_expr, /* Can be NULL if get_temp is NULL */
+ const char* power_expr, /* Can be NULL */
+ struct sdis_medium*** media_ptr,
+ sdis_medium_getter_T get_temp, /* Can be NULL */
+ int* out_has_power, /* Can be NULL */
+ unsigned* out_id)
+{
+ struct sdis_solid_shader solid_shader = SDIS_SOLID_SHADER_NULL;
+ struct sdis_data* data = NULL;
+ struct solid* solid_props;
+ size_t sz;
+ res_T res = RES_OK;
+
+ ASSERT(dev && lambda >= 0 && rho >= 0 && cp >= 0 && delta > 0
+ && media_ptr && out_id);
+ ASSERT(t_expr || !get_temp);
solid_shader.calorific_capacity = solid_get_calorific_capacity;
solid_shader.thermal_conductivity = solid_get_thermal_conductivity;
solid_shader.volumic_mass = solid_get_volumic_mass;
- solid_shader.delta_solid = solid_get_delta;
+ solid_shader.delta_solid = solid_get_delta;
#if Stardis_VERSION_MINOR == 3
solid_shader.delta_boundary = solid_get_delta_boundary;
#endif
- solid_shader.temperature = solid_get_temperature;
-
- /* Create the solid medium */
- for (i=0; i<stardis->geometry.medium_count; ++i){
- /* Create and setup the solid physical properties */
- int has;
- SDIS(data_create(dev, sizeof(struct solid), ALIGNOF(struct solid),
- release_solid_data, &data));
- solid_props = sdis_data_get(data); /* Fetch the allocated memory space */
- solid_props->cp = stardis->material[i].cp;
- solid_props->lambda = stardis->material[i].lambda;
- solid_props->rho = stardis->material[i].rho;
- solid_props->delta = stardis->material[i].delta;
- solid_props->power = NULL; /* just in case we go to error */
- res = set_solid_temperature_fn(solid_props, stardis->material[i].Tinit);
+ solid_shader.temperature = get_temp;
+ res = sdis_data_create(dev, sizeof(struct solid), ALIGNOF(struct solid),
+ release_solid_data, &data);
+ if (res != RES_OK) goto error;
+
+ solid_props = sdis_data_get(data); /* Fetch the allocated memory space */
+ solid_props->lambda = lambda;
+ solid_props->rho = rho;
+ solid_props->cp = cp;
+ solid_props->delta = delta;
+ solid_props->power = NULL; /* just in case we go to error */
+ if (t_expr) {
+ res = compile_expr_to_fn(&solid_props->t_init, t_expr, t_allowed, NULL);
if (res != RES_OK) {
fprintf(stderr, "Invalid initial temperature expression: %s\n",
- stardis->material[i].Tinit);
+ t_expr);
goto error;
}
- res = set_solid_power_fn(solid_props, stardis->material[i].power, &has);
+ }
+ if (power_expr) {
+ int has_power;
+ res = compile_expr_to_fn(&solid_props->power, power_expr, 1, &has_power);
if (res != RES_OK) {
fprintf(stderr, "Invalid volumic power expression: %s\n",
- stardis->material[i].power);
+ power_expr);
goto error;
}
- if (has) solid_shader.volumic_power = solid_get_power;
- SDIS(solid_create(dev, &solid_shader, data, sa_add(solid_medium, 1)));
- SDIS(data_ref_put(data)); data = NULL;
+ if (out_has_power) *out_has_power = has_power;
+ if(has_power) solid_shader.volumic_power = solid_get_power;
+ }
+ else {
+ if (out_has_power) *out_has_power = 0;
+ }
+ res = sdis_solid_create(dev, &solid_shader, data, sa_add(*media_ptr, 1));
+ if (res != RES_OK) goto error;
+ sz = sa_size(*media_ptr) - 1;
+ ASSERT(sz < INT_MAX);
+ *out_id = (unsigned)sz;
+
+end:
+ if (data) SDIS(data_ref_put(data));
+ return res;
+error:
+ goto end;
+}
+
+res_T
+stardis_compute(struct stardis* stardis, enum stardis_mode mode)
+{
+ res_T res = RES_OK;
+ struct sdis_device* dev = NULL;
+ struct sdis_medium** media = NULL;
+ struct sdis_interface** interfaces = NULL;
+ struct sdis_interface** interf_bytrg = NULL;
+
+ struct sdis_scene* scn = NULL;
+ struct sdis_estimator* estimator = NULL;
+ struct sdis_mc temperature;
+ struct sdis_accum_buffer* buf = NULL;
+ struct sdis_camera* cam = NULL;
+
+ struct htable_intface htable_interfaces; unsigned int_cpt = 0;
+ double pos[3] = {0,0,0};
+ double time[2] = { 0, 0};
+ size_t nfailures;
+ unsigned i = 0;
+ const int t_allowed = stardis->probe[3] < INF; /* Can condition use the variable t? */
+
+ res = sdis_device_create(NULL, &stardis->allocator, stardis->nthreads, 1, &dev);
+ if (res != RES_OK) goto error;
+
+ /* Create media and property holders found in descriptions */
+ for (i = 0; i < sa_size(stardis->descriptions); ++i) {
+ struct description* desc = stardis->descriptions + i;
+
+ switch (desc->type) {
+ case DESC_BOUND_H_FOR_SOLID:
+ /* Create an external fluid */
+ res = create_fluid(dev, 1, 1, t_allowed, desc->d.h_boundary.T, &media,
+ fluid_get_temperature, &desc->d.h_boundary.mat_id);
+ if (res != RES_OK) goto error;
+ break;
+ case DESC_BOUND_H_FOR_FLUID:
+ /* Create an external solid */
+ res = create_solid(dev, 1, 1, 1, 1, t_allowed, desc->d.h_boundary.T,
+ NULL, &media, solid_get_temperature, NULL, &desc->d.h_boundary.mat_id);
+ if (res != RES_OK) goto error;
+ break;
+ case DESC_BOUND_T_FOR_SOLID:
+ case DESC_BOUND_T_FOR_FLUID:
+ ASSERT(desc->d.t_boundary.T != NULL);
+ res = compile_expr_to_fn(&desc->d.t_boundary.te_temperature,
+ desc->d.t_boundary.T, 1, NULL);
+ if (res != RES_OK) {
+ fprintf(stderr, "Invalid boundary temperature expression: %s\n",
+ desc->d.h_boundary.T);
+ goto error;
+ }
+ /* Create an external solid */
+ res = create_solid(dev, 1, 1, 1, 1, t_allowed, NULL, NULL, &media, NULL,
+ NULL, &desc->d.t_boundary.mat_id);
+ if (res != RES_OK) goto error;
+ break;
+ case DESC_BOUND_F_FOR_SOLID:
+ ASSERT(desc->d.f_boundary.flux != NULL);
+ res = compile_expr_to_fn(&desc->d.f_boundary.te_flux,
+ desc->d.f_boundary.flux, 1, NULL);
+ if (res != RES_OK) {
+ fprintf(stderr, "Invalid boundary flux expression: %s\n",
+ desc->d.h_boundary.T);
+ goto error;
+ }
+ /* Create an external solid */
+ res = create_solid(dev, 1, 1, 1, 1, t_allowed, NULL, NULL, &media, NULL,
+ NULL, &desc->d.f_boundary.mat_id);
+ if (res != RES_OK) goto error;
+ break;
+ case DESC_SOLID_FLUID_CONNECT:
+ ASSERT(desc->d.sf_connect.specular_fraction >= 0
+ && desc->d.sf_connect.specular_fraction <= 1);
+ ASSERT(desc->d.sf_connect.emissivity >= 0);
+ ASSERT(desc->d.sf_connect.hc >= 0);
+ break;
+ case DESC_MAT_SOLID:
+ res = create_solid(dev,
+ stardis->descriptions[i].d.solid.lambda,
+ stardis->descriptions[i].d.solid.rho,
+ stardis->descriptions[i].d.solid.cp,
+ stardis->descriptions[i].d.solid.delta,
+ t_allowed,
+ stardis->descriptions[i].d.solid.Tinit,
+ stardis->descriptions[i].d.solid.power,
+ &media,
+ solid_get_tinit,
+ &stardis->descriptions[i].d.solid.has_power,
+ &desc->d.solid.solid_id);
+ if (res != RES_OK) goto error;
+ break;
+ case DESC_MAT_FLUID:
+ res = create_fluid(dev,
+ stardis->descriptions[i].d.fluid.rho,
+ stardis->descriptions[i].d.fluid.cp,
+ t_allowed,
+ stardis->descriptions[i].d.fluid.Tinit,
+ &media,
+ fluid_get_tinit,
+ &desc->d.fluid.fluid_id);
+ if (res != RES_OK) goto error;
+ break;
+ default: FATAL("Invalid type.\n");
+ }
}
- /* Setup the interface shader */
- interface_shader_front.temperature = interface_get_temperature;
- interface_shader_front.flux = interface_get_flux;
- interface_shader_front.emissivity = NULL;
- interface_shader_front.specular_fraction = NULL;
-
- interface_shader_back.temperature = interface_get_temperature;
- interface_shader_back.flux = NULL;
- interface_shader_back.emissivity = interface_get_emissivity;
- interface_shader_back.specular_fraction = interface_get_alpha;
-
- interface_shader_boundary.convection_coef = interface_get_convection_coef;
- interface_shader_boundary.front = interface_shader_front;
- interface_shader_boundary.back = interface_shader_back;
-
- for (i=0; i<sa_size(stardis->geometry.triangle); ++i){
- unsigned bound_id = stardis->geometry.triangle[i].bound_id;
-
- if (bound_id != UINT_MAX){ /* boundary interface */
- double hc = stardis->boundary[bound_id].hc;
- char* T = stardis->boundary[bound_id].T;
- char* flux = stardis->boundary[bound_id].flux;
- double emissivity =
- stardis->boundary[stardis->geometry.triangle[i].bound_id].emissivity;
- double alpha =
- stardis->boundary[stardis->geometry.triangle[i].bound_id].specular_fraction;
-
- SDIS(data_create(dev, sizeof(struct intface), ALIGNOF(struct intface),
+ /* Create interfaces */
+ htable_intface_init(&stardis->allocator, &htable_interfaces);
+ for (i=0; i<sa_size(stardis->geometry.triangle); ++i) {
+ struct triangle *trg = stardis->geometry.triangle + i;
+ struct int_descs int_descs = INT_DESCS_NULL;
+ struct sdis_interface** p_intface;
+ struct sdis_interface* intface;
+ struct sdis_medium* front_med = NULL;
+ struct sdis_medium* back_med = NULL;
+ struct sdis_interface_side_shader* fluid_side_shader = NULL;
+ unsigned fd = trg->front_description;
+ unsigned bd = trg->back_description;
+ unsigned cd = trg->connection_description;
+ int fluid_count = 0, solid_count = 0;
+ int solid_fluid_connection_count = 0, connection_count = 0, boundary_count = 0;
+
+ /* Create key */
+ int_descs.front = fd;
+ int_descs.back = bd;
+ int_descs.connect = trg->connection_description;
+
+ /* Search if interface already exists, or create it */
+ p_intface = htable_intface_find(&htable_interfaces, &int_descs);
+ if (p_intface) {
+ intface = *p_intface;
+ }
+ else {
+ /* create new interface and register */
+ struct sdis_data* data = NULL;
+ struct sdis_interface_shader interface_shader = SDIS_INTERFACE_SHADER_NULL;
+ struct intface* interface_props = NULL;
+ int front_defined = (fd != UINT_MAX);
+ int back_defined = (bd != UINT_MAX);
+ int connect_defined = (cd != UINT_MAX);
+ int incoherent = 0;
+ int_descs.id = int_cpt++;
+
+ SDIS(data_create(dev, sizeof(struct intface), ALIGNOF(struct intface),
release_interface_data, &data));
interface_props = sdis_data_get(data);
- if (hc > -1) {
- interface_props->hc = hc;
- interface_props->temperature = NULL;
- interface_props->flux = NULL;
- } else if (T != NULL){
- interface_props->flux = NULL;
- res = set_interface_temperature_fn(interface_props, T);
- if (res != RES_OK) {
- fprintf(stderr, "Invalid temperature expression: %s\n", T);
- goto error;
+ interface_props->temperature = NULL;
+ interface_props->flux = NULL;
+
+ if (front_defined) {
+ switch (stardis->descriptions[fd].type) {
+ case DESC_MAT_SOLID:
+ solid_count++;
+ front_med = media[fd];
+ break;
+ case DESC_MAT_FLUID:
+ fluid_count++;
+ front_med = media[fd];
+ fluid_side_shader = &interface_shader.front;
+ break;
+ default: FATAL("Invalid type.\n");
}
- } else if (flux != NULL){
- interface_props->temperature = NULL;
- res = set_interface_flux_fn(interface_props, flux);
- if (res != RES_OK) {
- fprintf(stderr, "Invalid flux expression: %s\n", flux);
- goto error;
+ }
+ if (back_defined) {
+ switch (stardis->descriptions[bd].type) {
+ case DESC_MAT_SOLID:
+ solid_count++;
+ back_med = media[bd];
+ break;
+ case DESC_MAT_FLUID:
+ fluid_count++;
+ back_med = media[bd];
+ fluid_side_shader = &interface_shader.back;
+ break;
+ default: FATAL("Invalid type.\n");
}
}
- interface_props->emissivity = emissivity;
- interface_props->alpha = alpha;
- SDIS(interface_create(dev,
- solid_medium[stardis->geometry.triangle[i].medium_front],
- fluid_medium[bound2fluid[bound_id]],
- &interface_shader_boundary, data, sa_add(interfaces,1)));
- SDIS(data_ref_put(data)); data = NULL;
+ if (connect_defined) {
+ const struct description* connect = stardis->descriptions + cd;
+ unsigned ext_id;
+ switch (connect->type) {
+ case DESC_BOUND_H_FOR_SOLID:
+ case DESC_BOUND_H_FOR_FLUID:
+ /* One among front and back should be defined, the other undefined;
+ * if not will raise an error below */
+ ext_id = connect->d.h_boundary.mat_id; /* External material id */
+ ASSERT(ext_id < sa_size(media));
+ ASSERT(sdis_medium_get_type(media[ext_id]) ==
+ (connect->type == DESC_BOUND_H_FOR_SOLID ? SDIS_FLUID : SDIS_SOLID));
+ if (front_defined && back_defined) incoherent = 1;
+ connection_count++;
+ boundary_count++;
+ 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_props->hc = connect->d.h_boundary.hc;
+ interface_props->emissivity = connect->d.h_boundary.emissivity;
+ interface_props->alpha = connect->d.h_boundary.specular_fraction;
+ break;
+ case DESC_BOUND_T_FOR_SOLID:
+ case DESC_BOUND_T_FOR_FLUID:
+ /* One among front and back should be defined, the other undefined;
+ * if not will raise an error below */
+ ext_id = connect->d.t_boundary.mat_id; /* External material id */
+ ASSERT(ext_id < sa_size(media));
+ ASSERT(sdis_medium_get_type(media[ext_id]) == SDIS_SOLID);
+ if (front_defined && back_defined) incoherent = 1;
+ connection_count++;
+ boundary_count++;
+ if (front_defined) back_med = media[ext_id];
+ else front_med = media[ext_id];
+ ASSERT(connect->d.t_boundary.te_temperature);
+ interface_props->temperature = connect->d.t_boundary.te_temperature;
+ if (connect->d.t_boundary.has_hc) {
+ 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_props->hc = connect->d.t_boundary.hc;
+ }
+ break;
+ case DESC_BOUND_F_FOR_SOLID:
+ /* One among front and back should be defined, the other undefined;
+ * if not will raise an error below */
+ connection_count++;
+ /* Reuse same solid for the external side */
+ if (front_defined) {
+ ASSERT(!back_defined);
+ back_med = front_med;
+ interface_shader.front.flux = interface_get_flux;
+ }
+ else {
+ ASSERT(back_defined);
+ front_med = back_med;
+ interface_shader.back.flux = interface_get_flux;
+ }
+ ASSERT(connect->d.f_boundary.te_flux);
+ interface_props->flux = connect->d.f_boundary.te_flux;
+ break;
+ case DESC_SOLID_FLUID_CONNECT:
+ /* Both front and back should be defined;
+ * if not will raise an error below */
+ connection_count++;
+ solid_fluid_connection_count++;
+ ASSERT(fluid_side_shader);
+ if (connect->d.sf_connect.has_hc) {
+ interface_shader.convection_coef = interface_get_convection_coef;
+ interface_shader.convection_coef_upper_bound = connect->d.sf_connect.hc_max;
+ interface_props->hc = connect->d.sf_connect.hc;
+ }
+ if (connect->d.sf_connect.has_emissivity) {
+ fluid_side_shader->emissivity = interface_get_emissivity;
+ interface_props->emissivity = connect->d.sf_connect.emissivity;
+ interface_props->alpha = connect->d.sf_connect.specular_fraction;
+ }
+ break;
+ default: FATAL("Invalid type.\n");
+ }
+ }
+
+ if(fluid_count == 2
+ || fluid_count + solid_count + connection_count < 2
+ || boundary_count ?
+ (fluid_count + solid_count != 1) : (fluid_count + solid_count != 2)
+ || (solid_fluid_connection_count && (fluid_count != 1 || solid_count != 1)))
+ incoherent = 1;
+ if (incoherent) {
+ /* Incoherent triangle description */
+ if(fluid_count == 2)
+ fprintf(stderr, "Incoherent triangle description\n");
+ res = RES_BAD_ARG;
+ goto error;
+ }
- } else { /* solid/solid interface */
- SDIS(interface_create(dev,
- solid_medium[stardis->geometry.triangle[i].medium_front],
- solid_medium[stardis->geometry.triangle[i].medium_back],
- &interface_shader_connection, NULL, sa_add(interfaces,1)));
+ res = sdis_interface_create(dev, front_med, back_med,
+ &interface_shader, data, &intface);
+ if (res != RES_OK) goto error;
+ SDIS(data_ref_put(data)); data = NULL;
+ res = htable_intface_set(&htable_interfaces, &int_descs, &intface);
+ if (res != RES_OK) goto error;
+ sa_push(interfaces, intface);
}
+ sa_push(interf_bytrg, intface);
}
+ stardis->geometry.interf_bytrg = interf_bytrg;
stardis->geometry.interfaces = interfaces;
- SDIS(scene_create(dev,
+ 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));
+ geometry_get_position, &stardis->geometry, &scn);
+ if (res != RES_OK) goto error;
if (mode == IR_COMPUTE){
@@ -627,7 +869,7 @@ stardis_compute(struct stardis* stardis, enum stardis_mode mode)
res = select_probe_type(scn,
stardis->geometry.triangle,
- stardis->material,
+ stardis->descriptions,
pos,
&iprim,
uv);
@@ -670,23 +912,9 @@ stardis_compute(struct stardis* stardis, enum stardis_mode mode)
}
end:
- for (i=0; i<sa_size(solid_medium); ++i){
- SDIS(medium_ref_put(solid_medium[i]));
- }
- sa_release(solid_medium);
-
- for (i=0; i<sa_size(fluid_medium); ++i){
- SDIS(medium_ref_put(fluid_medium[i]));
- }
- sa_release(fluid_medium);
-
- for (i=0; i<sa_size(interfaces); ++i){
- SDIS(interface_ref_put(interfaces[i]));
- }
- sa_release(interfaces);
- stardis->geometry.interfaces = NULL;
-
- sa_release(bound2fluid);
+ for (i=0; i<sa_size(media); ++i)
+ SDIS(medium_ref_put(media[i]));
+ sa_release(media);
if (cam) SDIS(camera_ref_put(cam));
if (buf) SDIS(accum_buffer_ref_put(buf));
@@ -696,6 +924,5 @@ end:
return res;
error:
- if(data) SDIS(data_ref_put(data));
goto end;
}
