commit 386113e9f5f05a0bedac52823651d311352306f4
parent bd17c0e6867e6a84c0f3f40063d2d3bb52bf3a1c
Author: Christophe Coustet <christophe.coustet@meso-star.com>
Date: Thu, 18 Jul 2019 12:11:11 +0200
Add missing files
Diffstat:
12 files changed, 2785 insertions(+), 0 deletions(-)
diff --git a/src/stardis-compute.h b/src/stardis-compute.h
@@ -0,0 +1,47 @@
+/* Copyright (C) 2018 |Meso|Star> (contact@meso-star.com)*/
+
+#ifndef SDIS_COMPUTE_H
+#define SDIS_COMPUTE_H
+
+#include <rsys/rsys.h>
+#include <rsys/hash_table.h>
+
+#define HTABLE_NAME weigth
+#define HTABLE_DATA double
+#define HTABLE_KEY unsigned
+#include <rsys/hash_table.h>
+
+struct te_expr;
+struct stardis;
+enum stardis_mode;
+
+struct counts {
+ size_t smed_count, fmed_count, tbound_count, hbound_count,
+ fbound_count, sfconnect_count;
+};
+#define COUNTS_NULL__ {\
+ 0, 0, 0, 0, 0, 0\
+ }
+static const struct counts COUNTS_NULL = COUNTS_NULL__;
+
+/* A type to limit variable use in tinyexpr expressions */
+enum var_prohibited_t {
+ ALL_VARS_ALLOWED = 0,
+ T_PROHIBITED = BIT(4),
+ XYZ_PROHIBITED = BIT(5),
+ NO_VAR_ALLOWED = T_PROHIBITED | XYZ_PROHIBITED
+};
+
+extern res_T
+compile_expr_to_fn
+ (struct te_expr** f,
+ const char* math_expr,
+ const int prohibited,
+ int* is_zero);
+
+extern res_T
+stardis_compute
+ (struct stardis* stardis,
+ enum stardis_mode mode);
+
+#endif
+\ No newline at end of file
diff --git a/src/stardis-fluid.c b/src/stardis-fluid.c
@@ -0,0 +1,142 @@
+#include "stardis-fluid.h"
+#include "stardis-compute.h"
+
+#include <sdis.h>
+#include <tinyexpr.h>
+#include<rsys/stretchy_array.h>
+
+/*******************************************************************************
+ * Local Functions
+ ******************************************************************************/
+
+static double
+fluid_get_calorific_capacity
+ (const struct sdis_rwalk_vertex* vtx,
+ struct sdis_data* data)
+{
+ const struct fluid* fluid_props = sdis_data_cget(data);
+ (void)vtx;
+ return fluid_props->cp;
+}
+
+static double
+fluid_get_volumic_mass
+ (const struct sdis_rwalk_vertex* vtx,
+ struct sdis_data* data)
+{
+ const struct fluid* fluid_props = sdis_data_cget(data);
+ (void)vtx;
+ return fluid_props->rho;
+}
+
+static double
+fluid_get_temperature
+ (const struct sdis_rwalk_vertex* vtx,
+ struct sdis_data* data)
+{
+ const struct fluid* fluid_props = sdis_data_cget(data);
+ char msg[128];
+ ASSERT(fluid_props->t_init || fluid_props->temp);
+ if (fluid_props->is_green || vtx->time > fluid_props->t0) {
+ /* Always use temp for Green mode, regardless of time */
+ if (fluid_props->temp)
+ return te_eval(fluid_props->temp, vtx);
+ else return -1;
+ }
+ /* Time is t0: use t_init */
+ if (fluid_props->t_init)
+ return te_eval(fluid_props->t_init, vtx);
+ /* Must have had t_init defined: error! */
+ if (fluid_props->name[0])
+ sprintf(msg,
+ "fluid_get_temperature: getting undefined Tinit (fluid '%s')\n",
+ fluid_props->name);
+ else
+ sprintf(msg, "fluid_get_temperature: getting undefined Tinit\n");
+ FATAL(msg);
+}
+
+static void
+release_fluid_data
+ (void* f)
+{
+ struct fluid* fluid = (struct fluid*)f;
+ te_free(fluid->t_init);
+ te_free(fluid->temp);
+}
+
+/*******************************************************************************
+ * Public Functions
+ ******************************************************************************/
+
+res_T
+create_fluid
+ (struct sdis_device* dev,
+ const char* name,
+ const double rho,
+ const double cp,
+ const int is_green,
+ const int is_outside,
+ const char* tinit_expr,
+ const char* t_expr,
+ struct sdis_medium*** media_ptr,
+ unsigned* out_id)
+{
+ res_T res = RES_OK;
+ struct sdis_fluid_shader fluid_shader = SDIS_FLUID_SHADER_NULL;
+ struct sdis_data* data = NULL;
+ struct fluid* fluid_props;
+ size_t sz;
+ /* Could be less restrictive if green output included positions/dates */
+ int prohibited = is_green ? NO_VAR_ALLOWED : ALL_VARS_ALLOWED;
+
+ ASSERT(dev && rho >= 0 && cp >= 0 && 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;
+ res = sdis_data_create(dev, sizeof(struct fluid), ALIGNOF(struct fluid),
+ release_fluid_data, &data);
+ if (res != RES_OK) goto error;
+
+ sz = sa_size(*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';
+ fluid_props->cp = cp;
+ fluid_props->rho = rho;
+ fluid_props->t_init = NULL;
+ fluid_props->temp = NULL;
+ fluid_props->is_green = is_green;
+ fluid_props->is_outside = is_outside;
+ fluid_props->t0 = 0;
+ fluid_props->id = (unsigned)sz;
+
+ if (tinit_expr) {
+ res = compile_expr_to_fn(&fluid_props->t_init, tinit_expr,
+ prohibited | T_PROHIBITED, NULL);
+ if (res != RES_OK) {
+ if (res == RES_BAD_ARG)
+ fprintf(stderr, "Invalid initial temperature expression: %s\n",
+ tinit_expr);
+ goto error;
+ }
+ }
+ if (t_expr) {
+ res = compile_expr_to_fn(&fluid_props->temp, t_expr, prohibited, NULL);
+ if (res != RES_OK) {
+ fprintf(stderr, "Invalid 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;
+ *out_id = fluid_props->id;
+
+end:
+ if (data) SDIS(data_ref_put(data));
+ return res;
+error:
+ goto end;
+}
diff --git a/src/stardis-fluid.h b/src/stardis-fluid.h
@@ -0,0 +1,43 @@
+/* Copyright (C) 2018 |Meso|Star> (contact@meso-star.com)*/
+
+#ifndef SDIS_FLUID_H
+#define SDIS_FLUID_H
+
+#include <sdis.h>
+#include <rsys/rsys.h>
+
+struct te_expr;
+struct sdis_device;
+struct sdis_medium;
+
+/*******************************************************************************
+ * Fluid data
+ ******************************************************************************/
+struct fluid {
+ char name[32];
+ double cp; /* Calorific capacity */
+ double rho; /* Volumic mass */
+ /* Compute mode */
+ int is_green, is_outside;
+ double t0;
+ /* TinyExpr stuff to compute temperature */
+ struct te_expr *temp;
+ struct te_expr *t_init;
+ /* ID */
+ unsigned id;
+};
+
+extern res_T
+create_fluid
+ (struct sdis_device* dev,
+ const char* name,
+ const double rho,
+ const double cp,
+ const int is_green,
+ const int is_outside,
+ const char* tinit_expr,
+ const char* t_expr,
+ struct sdis_medium*** media_ptr,
+ unsigned* out_id);
+
+#endif
diff --git a/src/stardis-intface.c b/src/stardis-intface.c
@@ -0,0 +1,360 @@
+#include "stardis-intface.h"
+#include "stardis-app.h"
+#include "stardis-output.h"
+
+#include <tinyexpr.h>
+#include <sdis.h>
+#include <rsys/stretchy_array.h>
+
+/*******************************************************************************
+ * Local Functions
+ ******************************************************************************/
+
+static double
+interface_get_convection_coef
+ (const struct sdis_interface_fragment* frag,
+ struct sdis_data* data)
+{
+ const struct intface* interface_props = sdis_data_cget(data);
+ (void)frag;
+ return interface_props->hc;
+}
+
+static double
+interface_get_temperature
+ (const struct sdis_interface_fragment* frag,
+ struct sdis_data* data)
+{
+ const struct intface* interface_props = sdis_data_cget(data);
+
+ if (interface_props->temperature == NULL)
+ return -1;
+
+ return te_eval(interface_props->temperature, frag);
+}
+
+static double
+interface_get_flux
+ (const struct sdis_interface_fragment* frag,
+ struct sdis_data* data)
+{
+ const struct intface* interface_props = sdis_data_cget(data);
+
+ if (interface_props->flux == NULL)
+ return SDIS_FLUX_NONE;
+
+ return te_eval(interface_props->flux, frag);
+}
+
+static double
+interface_get_emissivity
+ (const struct sdis_interface_fragment* frag,
+ struct sdis_data* data)
+{
+ const struct intface* interface_props = sdis_data_cget(data);
+ (void)frag;
+ return interface_props->emissivity;
+}
+
+static double
+interface_get_alpha
+ (const struct sdis_interface_fragment* frag,
+ struct sdis_data* data)
+{
+ const struct intface* interface_props = sdis_data_cget(data);
+ (void)frag;
+ return interface_props->alpha;
+}
+
+static void
+release_interface_data
+ (void* s)
+{
+ struct intface* intface = (struct intface*)s;
+ te_free(intface->temperature);
+ te_free(intface->flux);
+}
+
+/*******************************************************************************
+ * Public Functions
+ ******************************************************************************/
+
+res_T
+create_intface
+ (struct sdis_device* dev,
+ unsigned tr_idx,
+ struct htable_intface* htable_interfaces,
+ struct geometry* geometry,
+ const struct description* descriptions,
+ struct sdis_medium* const* media)
+{
+ struct triangle *trg;
+ 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;
+ struct sdis_data* data = NULL;
+ unsigned fd, bd, cd;
+ int fluid_count = 0, solid_count = 0;
+ int solid_fluid_connection_count = 0, connection_count = 0, boundary_count = 0;
+ res_T res = RES_OK;
+
+ ASSERT(dev && htable_interfaces && geometry && descriptions && media);
+
+ trg = geometry->triangle + tr_idx;
+ fd = trg->front_description;
+ bd = trg->back_description;
+ cd = trg->connection_description;
+
+ /* 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_interface_shader interface_shader = SDIS_INTERFACE_SHADER_NULL;
+ struct intface* interface_props = NULL;
+ unsigned id;
+ int front_defined = (fd != UINT_MAX);
+ int back_defined = (bd != UINT_MAX);
+ int connect_defined = (cd != UINT_MAX);
+
+ SDIS(data_create(dev, sizeof(struct intface), ALIGNOF(struct intface),
+ release_interface_data, &data));
+ interface_props = sdis_data_get(data);
+ interface_props->temperature = NULL;
+ interface_props->flux = NULL;
+ interface_props->front_boundary_id = UINT_MAX;
+ interface_props->back_boundary_id = UINT_MAX;
+
+ if (front_defined) {
+ switch (descriptions[fd].type) {
+ case DESC_MAT_SOLID:
+ id = descriptions[fd].d.solid.solid_id;
+ solid_count++;
+ front_med = media[id];
+ break;
+ case DESC_MAT_FLUID:
+ fluid_count++;
+ id = descriptions[fd].d.fluid.fluid_id;
+ front_med = media[id];
+ fluid_side_shader = &interface_shader.front;
+ break;
+ default: FATAL("Invalid type.\n");
+ }
+ }
+ if (back_defined) {
+ switch (descriptions[bd].type) {
+ case DESC_MAT_SOLID:
+ id = descriptions[bd].d.solid.solid_id;
+ solid_count++;
+ back_med = media[id];
+ break;
+ case DESC_MAT_FLUID:
+ fluid_count++;
+ id = descriptions[bd].d.fluid.fluid_id;
+ back_med = media[id];
+ fluid_side_shader = &interface_shader.back;
+ break;
+ default: FATAL("Invalid type.\n");
+ }
+ }
+ 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) {
+ ASSERT(front_defined && back_defined);
+ fprintf(stderr,
+ "Can only define a DESC_SOLID_FLUID_CONNECT between a fluid and a solid\n");
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ }
+ else {
+ if (front_defined == back_defined) {
+ fprintf(stderr,
+ "Cannot define a boundary between 2 %s regions\n",
+ front_defined ? "defined" : "undefined");
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ def_medium = front_defined ? front_med : back_med;
+ 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) {
+ fprintf(stderr,
+ "Can only define a DESC_BOUND_H_FOR_FLUID boundary for a fluid region\n");
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ type_checked = 1;
+ /* fall through */
+ case DESC_BOUND_H_FOR_SOLID:
+ 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");
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ 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));
+ 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;
+ ASSERT(!fluid_side_shader); /* Cause defined medium is solid */
+ if (connect->d.h_boundary.has_emissivity) {
+ fluid_side_shader =
+ front_defined ? &interface_shader.back : &interface_shader.front;
+ fluid_side_shader->emissivity = interface_get_emissivity;
+ interface_props->emissivity = connect->d.h_boundary.emissivity;
+ interface_props->alpha = connect->d.h_boundary.specular_fraction;
+ }
+ break;
+ case DESC_BOUND_T_FOR_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;
+ goto error;
+ }
+ type_checked = 1;
+ /* fall through */
+ case DESC_BOUND_T_FOR_SOLID:
+ 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");
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ 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_FLUID);
+ connection_count++;
+ boundary_count++;
+ if (front_defined) {
+ back_med = media[ext_id];
+ /* We set the known T inside
+ * TODO: should be outside to allow contact resistances */
+ interface_shader.front.temperature = interface_get_temperature;
+ }
+ else {
+ front_med = media[ext_id];
+ /* We set the known T inside
+ * TODO: should be outside to allow contact resistances */
+ interface_shader.back.temperature = interface_get_temperature;
+ }
+ ASSERT(connect->d.t_boundary.te_temperature);
+ interface_props->temperature
+ = te_duplicate(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:
+ 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;
+ goto error;
+ }
+ connection_count++;
+ boundary_count++;
+ if (front_defined) {
+ back_med = media[connect->d.f_boundary.mat_id];
+ interface_shader.front.flux = interface_get_flux;
+ }
+ else {
+ front_med = media[connect->d.f_boundary.mat_id];
+ interface_shader.back.flux = interface_get_flux;
+ }
+ ASSERT(connect->d.f_boundary.te_flux);
+ interface_props->flux = te_duplicate(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 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");
+ else print_description(stderr, descriptions + fd);
+ fprintf(stderr, "Back: ");
+ if (!back_defined) fprintf(stderr, "undefined\n");
+ else print_description(stderr, descriptions + bd);
+ fprintf(stderr, "Connection: ");
+ if (!connect_defined) fprintf(stderr, "undefined\n");
+ else print_description(stderr, descriptions + cd);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ res = sdis_interface_create(dev, front_med, back_med,
+ &interface_shader, data, &intface);
+ SDIS(data_ref_put(data)); data = NULL;
+ 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;
+ }
+ sa_push(geometry->interf_bytrg, intface);
+
+end:
+ return res;
+error:
+ goto end;
+}
+\ No newline at end of file
diff --git a/src/stardis-intface.h b/src/stardis-intface.h
@@ -0,0 +1,67 @@
+/* Copyright (C) 2018 |Meso|Star> (contact@meso-star.com)*/
+
+#ifndef SDIS_INTFACE_H
+#define SDIS_INTFACE_H
+
+#include <rsys/hash_table.h>
+
+#include <limits.h>
+
+struct te_expr;
+struct sdis_data;
+struct sdis_device;
+struct geometry;
+struct description;
+struct sdis_medium;
+
+/*******************************************************************************
+ * Interface data
+ ******************************************************************************/
+struct intface {
+ double hc; /* Convection coefficient */
+ double emissivity;
+ double alpha;
+ /* TinyExpr stuff to compute temperature & flux */
+ struct te_expr *temperature;
+ struct te_expr *flux;
+ /* IDs */
+ unsigned front_boundary_id, back_boundary_id;
+};
+
+/* Declare the hash table that map an interface to its descriptor */
+struct int_descs {
+ unsigned front, back, connect;
+};
+#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)
+{
+ return (char)(a->front == b->front && a->back == b->back
+ && a->connect == b->connect);
+}
+
+static INLINE size_t
+hash_desc(struct int_descs const* key)
+{
+ return (size_t)hash_fnv64(key, 3 * sizeof(unsigned));
+}
+
+#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>
+
+extern res_T
+create_intface
+ (struct sdis_device* dev,
+ unsigned tr_idx,
+ struct htable_intface* htable_interfaces,
+ struct geometry* geometry,
+ const struct description* descriptions,
+ struct sdis_medium* const* media);
+
+#endif
diff --git a/src/stardis-output.c b/src/stardis-output.c
@@ -0,0 +1,775 @@
+#include "stardis-output.h"
+#include "stardis-compute.h"
+#include "stardis-fluid.h"
+#include "stardis-solid.h"
+#include "stardis-intface.h"
+#include "stardis-app.h"
+
+#include <sdis.h>
+
+#include<star/senc.h>
+
+#include <rsys/math.h>
+#include <rsys/mem_allocator.h>
+#include <rsys/dynamic_array_uint.h>
+
+#include <stdio.h>
+
+struct w_ctx {
+ const struct description* desc;
+ struct htable_weigth weigths;
+};
+
+#include <rsys/hash_table.h>
+#define HTABLE_NAME vrtx_rank
+#define HTABLE_DATA unsigned
+#define HTABLE_KEY unsigned
+#include <rsys/hash_table.h>
+
+/*******************************************************************************
+ * Local Functions
+ ******************************************************************************/
+
+static res_T
+print_power_term
+ (struct sdis_medium* mdm,
+ const double power_term,
+ void* ctx)
+{
+ res_T res = RES_OK;
+ struct sdis_data* data = NULL;
+ enum sdis_medium_type type;
+ unsigned id;
+ ASSERT(mdm && ctx); (void)ctx;
+
+ data = sdis_medium_get_data(mdm);
+ type = sdis_medium_get_type(mdm);
+
+ switch (type) {
+ case SDIS_FLUID: {
+ FATAL("Unexpected power term in fluid");
+ }
+ case SDIS_SOLID: {
+ struct solid* d__ = sdis_data_get(data);
+ id = d__->id;
+ ASSERT(id == ((struct w_ctx*)ctx)->desc[id].d.solid.solid_id);
+ printf("\tS\t%u\t%g", id, power_term);
+ break;
+ }
+ default: FATAL("Unreachable code.\n"); break;
+ }
+ return res;
+}
+
+static res_T
+get_flux_terms
+ (struct sdis_interface* interf,
+ const enum sdis_side side,
+ const double flux_term,
+ void* ctx)
+{
+ res_T res = RES_OK;
+ struct sdis_data* data = NULL;
+ struct intface* d__;
+ unsigned id;
+ struct w_ctx* w_ctx = ctx;
+
+ 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;
+
+ switch (w_ctx->desc[id].type) {
+ case DESC_BOUND_T_FOR_SOLID:
+ case DESC_BOUND_T_FOR_FLUID:
+ case DESC_BOUND_H_FOR_SOLID:
+ case DESC_BOUND_H_FOR_FLUID:
+ 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);
+ w = htable_weigth_find(&w_ctx->weigths, &id);
+ if (w)
+ *w += flux_term;
+ else {
+ res = htable_weigth_set(&w_ctx->weigths, &id, &flux_term);
+ }
+ break;
+ }
+ default: FATAL("Unreachable code.\n"); break;
+ }
+ return res;
+}
+
+res_T
+dump_path
+ (const struct sdis_heat_path* path,
+ void* context)
+{
+ res_T res = RES_OK;
+ FILE* stream = context;
+ enum sdis_heat_path_flag status = SDIS_HEAT_PATH_NONE;
+ size_t i, vcount;
+
+ /* Header */
+ fprintf(stream, "---\n");
+ fprintf(stream, "# vtk DataFile Version 2.0\n");
+ fprintf(stream, "Heat paths\n");
+ fprintf(stream, "ASCII\n");
+ fprintf(stream, "DATASET POLYDATA\n");
+ /* Write path positions */
+ sdis_heat_path_get_vertices_count(path, &vcount);
+ fprintf(stream, "POINTS %lu double\n", (unsigned long)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;
+ fprintf(stream, "%g %g %g\n", SPLIT3(vtx.P));
+ }
+ /* Write the segment of the path */
+ fprintf(stream, "LINES %lu %lu\n", 1lu, (unsigned long)(1 + vcount));
+ fprintf(stream, "%lu", (unsigned long)vcount);
+ FOR_EACH(i, 0, vcount) fprintf(stream, " %lu", (unsigned long)i);
+ fprintf(stream, "\n");
+ fprintf(stream, "POINT_DATA %lu\n", (unsigned long)vcount);
+ /* Write the type of the random walk vertices */
+ fprintf(stream, "SCALARS Vertex_Type float 1\n");
+ 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;
+ 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;
+ case SDIS_HEAT_VERTEX_RADIATIVE: fprintf(stream, "1.0\n"); break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+ }
+ fprintf(stream, "LOOKUP_TABLE vertex_type 3\n");
+ fprintf(stream, "0.0 1.0 1.0 1.0\n"); /* 0.0 = Magenta: conduction */
+ fprintf(stream, "1.0 1.0 0.0 1.0\n"); /* 0.5 = Yellow: convection */
+ fprintf(stream, "1.0 0.0 1.0 1.0\n"); /* 1.0 = Purple: radiative */
+ /* Write the weights of the random walk vertices */
+ fprintf(stream, "SCALARS Weight double 1\n");
+ 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;
+ fprintf(stream, "%g\n", vtx.weight);
+ }
+ /* Write the time of the random walk vertices */
+ fprintf(stream, "SCALARS Time double 1\n");
+ 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;
+ fprintf(stream, "%g\n", IS_INF(vtx.time) ? FLT_MAX : vtx.time);
+ }
+ /* Write path type */
+ fprintf(stream, "CELL_DATA %lu\n", 1lu);
+ 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;
+ switch (status) {
+ case SDIS_HEAT_PATH_SUCCEED: fprintf(stream, "0.0\n"); break;
+ case SDIS_HEAT_PATH_FAILED: fprintf(stream, "1.0\n"); break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+ fprintf(stream, "LOOKUP_TABLE path_type 2\n");
+ fprintf(stream, "0.0 0.0 1.0 1.0\n"); /* 0.0 = Blue: success */
+ fprintf(stream, "1.0 0.0 0.0 1.0\n"); /* 1.0 = Red: failure */
+
+end:
+ return res;
+error:
+ goto end;
+}
+
+/*******************************************************************************
+ * Public Functions
+ ******************************************************************************/
+
+res_T
+print_sample
+ (struct sdis_green_path* path,
+ void* ctx)
+{
+ res_T res = RES_OK;
+ struct sdis_point pt = SDIS_POINT_NULL;
+ struct sdis_data* data = NULL;
+ enum sdis_medium_type type;
+ struct htable_weigth_iterator it, end;
+ unsigned id;
+ size_t pcount, fcount;
+ struct w_ctx* w_ctx = ctx;
+ CHK(path && ctx);
+
+ res = sdis_green_path_get_limit_point(path, &pt);
+ if (res != RES_OK) goto error;
+
+ /* For each path, prints:
+ * # end #power_terms #flux_terms power_term_1 ... power_term_n flux_term_1 ... flux_term_n
+ * with:
+ * - end = end_type end_id; end_type = T | H | X | R | F | S
+ * - power_term_i = power_type_i power_id_i factor_i
+ * - flux_term_i = flux_id_i factor_i
+ */
+
+ switch (pt.type) {
+ case SDIS_FRAGMENT: {
+ struct intface* d__;
+ data = sdis_interface_get_data(pt.data.itfrag.intface);
+ d__ = sdis_data_get(data);
+ 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) {
+ case DESC_BOUND_T_FOR_SOLID:
+ case DESC_BOUND_T_FOR_FLUID:
+ ASSERT(id == w_ctx->desc[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);
+ printf("H\t%u", id);
+ break;
+ case DESC_BOUND_F_FOR_SOLID:
+ ASSERT(id == w_ctx->desc[id].d.f_boundary.mat_id);
+ printf("X\t%u", id);
+ break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+ break;
+ }
+ 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) {
+ /* Radiative output */
+ printf("R\t%u", (unsigned)sa_size(w_ctx->desc));
+ }
+ 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)
+ /* If outside the model and in a fluid with known temperature,
+ * its a fluid attached to a H boundary */
+ printf("H\t%u", id);
+ /* In a standard fluid with known temperature */
+ else printf("F\t%u", id);
+
+ }
+ else {
+ 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(!d__->is_outside); /* FIXME: what if in external solid? */
+ printf("S\t%u", id);
+ }
+ break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+
+ res = sdis_green_function_get_power_terms_count(path, &pcount);
+ if (res != RES_OK) goto error;
+
+ 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;
+ 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;
+
+ htable_weigth_begin(&w_ctx->weigths, &it);
+ htable_weigth_end(&w_ctx->weigths, &end);
+ while (!htable_weigth_iterator_eq(&it, &end)) {
+ double* w = htable_weigth_iterator_data_get(&it);
+ unsigned* k = htable_weigth_iterator_key_get(&it);
+ printf("\t%u\t%g", *k, *w);
+ htable_weigth_iterator_next(&it);
+ }
+ printf("\n");
+
+end:
+ return res;
+error:
+ goto end;
+}
+
+res_T
+dump_image
+ (const struct sdis_estimator_buffer* buf)
+{
+ res_T res = RES_OK;
+ size_t definition[2];
+ double* temps = NULL;
+ size_t ix, iy;
+
+ CHK(buf != NULL);
+ res = sdis_estimator_buffer_get_definition(buf, definition);
+ if (res != RES_OK) goto error;
+
+ temps = mem_alloc(definition[0] * definition[1] * sizeof(double));
+ if (temps == NULL) {
+ res = RES_MEM_ERR;
+ goto error;
+ }
+
+ /* Compute the per pixel temperature */
+ fprintf(stdout, "# vtk DataFile Version 2.0\nvtk output\nASCII\nDATASET STRUCTURED_POINTS\n");
+ fprintf(stdout, "DIMENSIONS %zu %zu 1\n", definition[0], definition[1]);
+ fprintf(stdout, "ORIGIN 0 0 0\n");
+ fprintf(stdout, "SPACING 1 1 1\n");
+ fprintf(stdout, "POINT_DATA %zu\n", definition[0] * definition[1]);
+ fprintf(stdout, "SCALARS temperature_estimate float 1\n");
+ fprintf(stdout, "LOOKUP_TABLE default\n");
+ FOR_EACH(iy, 0, definition[1]) {
+ 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;
+ fprintf(stdout, "%f\n", T.E);
+ }
+ }
+ fprintf(stdout, "SCALARS temperature_std_dev float 1\n");
+ fprintf(stdout, "LOOKUP_TABLE default\n");
+ FOR_EACH(iy, 0, definition[1]) {
+ 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;
+ fprintf(stdout, "%f\n", T.SE);
+ }
+ }
+ fprintf(stdout, "SCALARS computation_time float 1\n");
+ fprintf(stdout, "LOOKUP_TABLE default\n");
+ FOR_EACH(iy, 0, definition[1]) {
+ 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;
+ fprintf(stdout, "%f\n", time.E);
+ }
+ }
+ fprintf(stdout, "SCALARS computation_time_std_dev float 1\n");
+ fprintf(stdout, "LOOKUP_TABLE default\n");
+ FOR_EACH(iy, 0, definition[1]) {
+ 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;
+ fprintf(stdout, "%f\n", time.SE);
+ }
+ }
+ fprintf(stdout, "SCALARS temperature_failures_count int 1\n");
+ fprintf(stdout, "LOOKUP_TABLE default\n");
+ FOR_EACH(iy, 0, definition[1]) {
+ 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;
+ fprintf(stdout, "%zu\n", nfails);
+ }
+ }
+ mem_rm(temps);
+
+end:
+ return res;
+error:
+ goto end;
+}
+
+res_T
+dump_green
+ (struct sdis_green_function* green,
+ const struct counts* counts,
+ const struct description* descriptions,
+ const double radiative_temps[2])
+{
+ res_T res = RES_OK;
+ size_t ok_count, failed_count;
+ struct w_ctx w_ctx;
+ unsigned i;
+
+ 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;
+
+ /* Output counts */
+ printf("---BEGIN GREEN---\n");
+ printf("# #solids #fluids #t_boundaries #h_boundaries #f_boundaries #ok #failures\n");
+ printf("%zu %zu %zu %zu %zu %zu %zu\n",
+ counts->smed_count, counts->fmed_count,
+ counts->tbound_count, counts->hbound_count, counts->fbound_count,
+ ok_count, failed_count);
+
+ /* List Media */
+ 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;
+ const struct mat_solid* sl;
+ if (desc->type != DESC_MAT_SOLID) continue;
+ sl = &desc->d.solid;
+ printf("%u\t%s\t%g\t%g\t%g\t%s\n",
+ i, sl->name, sl->lambda, sl->rho, sl->cp, (sl->has_power ? sl->power : "0"));
+ }
+ }
+ 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;
+ 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);
+ }
+ }
+
+ /* List Boundaries */
+ 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;
+ const struct t_boundary* bd;
+ 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%s\n", i, bd->name, bd->T);
+ }
+ }
+ 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;
+ const struct h_boundary* bd;
+ 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%s\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->T);
+ }
+ }
+ 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;
+ const struct f_boundary* bd;
+ if (desc->type != DESC_BOUND_F_FOR_SOLID) continue;
+ bd = &desc->d.f_boundary;
+ printf("%u\t%s\t%s\n",
+ i, bd->name, bd->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("# Samples\n");
+ printf("# end #power_terms #flux_terms power_term_1 ... power_term_n flux_term_1 ... flux_term_n\n");
+ printf("# end = end_type end_id; end_type = T | H | X | R | F | S\n");
+ printf("# power_term_i = power_type_i power_id_i factor_i\n");
+ printf("# flux_term_i = flux_id_i factor_i\n");
+
+ w_ctx.desc = descriptions;
+ htable_weigth_init(NULL, &w_ctx.weigths);
+
+ res = sdis_green_function_for_each_path(green, print_sample, &w_ctx);
+ htable_weigth_release(&w_ctx.weigths);
+ if (res != RES_OK) goto error;
+
+
+ printf("---END GREEN---\n");
+
+ CHK(sdis_green_function_ref_put(green) == RES_OK);
+
+end:
+ return res;
+error:
+ goto end;
+}
+
+res_T
+create_edge_file_if
+ (struct sdis_scene* scn,
+ struct mem_allocator* allocator)
+{
+ 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;
+ struct senc_descriptor* analyze = NULL;
+ struct darray_uint comps;
+ struct darray_uint edges;
+ unsigned cc_count = 0;
+ unsigned *cc, *ee, uimax = UINT_MAX;
+
+ res = sdis_scene_get_analysis(scn, &analyze);
+ if (res != RES_OK) goto error;
+ res = sdis_scene_release_analysis(scn);
+ if (res != RES_OK) goto error;
+
+ ASSERT(scn && allocator);
+
+ res = senc_descriptor_get_frontier_segments_count(analyze, &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");
+
+ htable_vrtx_rank_init(allocator, &ranks);
+ darray_uint_init(allocator, &comps);
+ darray_uint_init(allocator, &edges);
+ structs_initialized = 1;
+
+ /* Exit vertices, with a new numbering scheme */
+ vcount = 0;
+ FOR_EACH(i, 0, scount) {
+ unsigned v, edge[2], rk[2];
+ res = senc_descriptor_get_frontier_segment(analyze, 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++;
+ }
+ }
+ }
+
+ /* 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 = senc_descriptor_get_global_vertex(analyze, 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 = senc_descriptor_get_frontier_segment(analyze, 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;
+ }
+
+end:
+ if (obj) fclose(obj);
+ if (analyze) senc_descriptor_ref_put(analyze);
+ if (structs_initialized) {
+ htable_vrtx_rank_release(&ranks);
+ darray_uint_release(&comps);
+ darray_uint_release(&edges);
+ }
+ return res;
+error:
+ goto end;
+}
+
+res_T
+print_single_MC_result
+ (const enum stardis_mode mode,
+ struct sdis_estimator* estimator,
+ struct stardis* stardis)
+{
+ res_T res = RES_OK;
+ struct sdis_mc temperature;
+ size_t nfailures;
+
+ ASSERT(estimator && stardis);
+
+ /* Fetch the estimation data */
+ res = sdis_estimator_get_temperature(estimator, &temperature);
+ if (res != RES_OK) goto error;
+ res = sdis_estimator_get_failure_count(estimator, &nfailures);
+ if (res != RES_OK) goto error;
+
+ /* Print the results */
+ switch (mode & COMPUTE_MODES) {
+ case PROBE_COMPUTE:
+ case PROBE_COMPUTE_ON_INTERFACE:
+ printf("Temperature at [%g, %g, %g], t=%g = %g +/- %g\n",
+ SPLIT4(stardis->probe),
+ temperature.E, /* Expected value */
+ temperature.SE); /* Standard error */
+ break;
+ case MEDIUM_COMPUTE:
+ printf("Temperature in medium %s at t=%g = %g +/- %g\n",
+ stardis->solve_name, stardis->probe[3],
+ temperature.E, /* Expected value */
+ temperature.SE); /* Standard error */
+ break;
+ case BOUNDARY_COMPUTE:
+ printf("Temperature at boundary %s at t=%g = %g +/- %g\n",
+ stardis->solve_name, stardis->probe[3],
+ temperature.E, /* Expected value */
+ temperature.SE); /* Standard error */
+ break;
+ default: FATAL("Invalid mode.");
+ }
+ printf("#failures: %lu/%lu\n",
+ (unsigned long)nfailures,
+ (unsigned long)stardis->N);
+ if (nfailures)
+ fprintf(stderr, "#failures: %lu/%lu\n",
+ (unsigned long)nfailures,
+ (unsigned long)stardis->N);
+
+ /* Dump paths according to user settings */
+ res = sdis_estimator_for_each_path(estimator, dump_path, stdout);
+ if (res != RES_OK) goto error;
+end:
+ return res;
+error:
+ goto end;
+}
+
+/* TODO: rewrite dump! */
+res_T
+dump_vtk
+ (FILE* output,
+ const struct geometry* geometry)
+{
+ res_T res = RES_OK;
+ unsigned i;
+ struct vertex* vtx = geometry->vertex;
+ struct triangle* tri = geometry->triangle;
+
+ if (!output) {
+ fprintf(stderr, "Invalid output file to dump vtk.\n");
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ fprintf(output, "# vtk DataFile Version 2.0\nvtk output\nASCII\nDATASET POLYDATA\n");
+ fprintf(output, "POINTS %i float\n\n", (int)sa_size(vtx));
+ for (i = 0; i < sa_size(vtx); ++i) {
+ fprintf(output, "%f %f %f\n", SPLIT3(vtx[i].xyz));
+ }
+ fprintf(output, "\nPOLYGONS %i %i\n", (int)sa_size(tri), (int)(4 * sa_size(tri)));
+ for (i = 0; i < sa_size(tri); ++i) {
+ fprintf(output, "3 %i %i %i\n", SPLIT3(tri[i].indices.data));
+ }
+ fprintf(output, "\nCELL_DATA %i \n", (int)sa_size(tri));
+ 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].front_description);
+ }
+ 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].back_description);
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+void
+print_trg_as_obj
+ (FILE* stream,
+ const struct vertex* vertices,
+ const unsigned* indices)
+{
+ 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));
+}
diff --git a/src/stardis-output.h b/src/stardis-output.h
@@ -0,0 +1,62 @@
+/* Copyright (C) 2018 |Meso|Star> (contact@meso-star.com)*/
+
+#ifndef SDIS_OUTPUT_H
+#define SDIS_OUTPUT_H
+
+#include <rsys/rsys.h>
+#include <stdio.h>
+
+struct sdis_green_path;
+struct sdis_heat_path;
+struct sdis_estimator_buffer;
+struct sdis_green_function;
+struct counts;
+struct description;
+struct sdis_scene;
+struct mem_allocator;
+enum stardis_mode;
+struct sdis_estimator;
+struct stardis;
+struct geometry;
+struct vertex;
+
+extern res_T
+print_sample
+ (struct sdis_green_path* path,
+ void* ctx);
+
+extern res_T
+dump_image
+ (const struct sdis_estimator_buffer* buf);
+
+extern res_T
+dump_green
+ (struct sdis_green_function* green,
+ const struct counts* counts,
+ const struct description* descriptions,
+ const double radiative_temps[2]);
+
+extern res_T
+create_edge_file_if
+ (struct sdis_scene* scn,
+ struct mem_allocator* allocator);
+
+extern res_T
+print_single_MC_result
+ (const enum stardis_mode mode,
+ struct sdis_estimator* estimator,
+ struct stardis* stardis);
+
+/* TODO: rewrite dump! */
+extern res_T
+dump_vtk
+ (FILE* output,
+ const struct geometry* geometry);
+
+extern void
+print_trg_as_obj
+ (FILE* stream,
+ const struct vertex* vertices,
+ const unsigned* indices);
+
+#endif
diff --git a/src/stardis-parsing.c b/src/stardis-parsing.c
@@ -0,0 +1,899 @@
+#include "stardis-parsing.h"
+#include "stardis-app.h"
+#include "stardis-version.h"
+
+#include <rsys/cstr.h>
+#include <sdis_version.h>
+
+#include <getopt.h>
+#include <stdlib.h>
+#include <stdio.h>
+
+/*******************************************************************************
+ * Local Functions
+ ******************************************************************************/
+
+static char**
+split_line
+ (char* a_str,
+ const char a_delim)
+{
+ char** result = 0;
+ size_t count = 0;
+ char* tmp = a_str;
+ char* last_comma = 0;
+ char delim[2];
+ delim[0] = a_delim;
+ delim[1] = 0;
+
+ ASSERT(a_str);
+
+ /* Count how many elements will be extracted. */
+ while (*tmp)
+ {
+ if (a_delim == *tmp)
+ {
+ count++;
+ last_comma = tmp;
+ }
+ tmp++;
+ }
+
+ /* Add space for trailing token. */
+ count += last_comma < (a_str + strlen(a_str) - 1);
+
+ /* Add space for terminating null string so caller
+ knows where the list of returned strings ends. */
+ count++;
+
+ result = malloc(sizeof(char*) * count);
+
+ if (result)
+ {
+ size_t idx = 0;
+ char* token = strtok(a_str, delim);
+
+ while (token)
+ {
+ ASSERT(idx < count);
+#ifdef COMPILER_CL
+ *(result + idx++) = _strdup(token);
+#else
+ *(result + idx++) = strdup(token);
+#endif
+ token = strtok(0, delim);
+ }
+ ASSERT(idx == count - 1);
+ *(result + idx) = 0;
+ }
+
+ return result;
+}
+
+#ifdef COMPILER_GCC
+static char*
+_strupr
+ (char* s)
+{
+ char* ptr;
+ for (ptr = s; *ptr; ++ptr) {
+ int tmp = toupper(*ptr);
+ ASSERT(tmp == (char)tmp);
+ *ptr = (char)tmp;
+ }
+ return s;
+}
+#endif
+
+/*******************************************************************************
+ * Public Functions
+ ******************************************************************************/
+
+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 & DUMP_VTK) { found++; res = 'd'; }
+ ASSERT(found == 1);
+ return res;
+}
+
+void print_multiple_modes
+ (FILE* stream,
+ const int modes)
+{
+ int b = 0, fst = 1;
+ enum stardis_mode m = UNDEF_MODE;
+ ASSERT(stream);
+ do {
+ m = BIT(b++);
+ if (m & modes) {
+ fprintf(stream, (fst ? " -%c" : ", -%c"), mode_option(m));
+ fst = 0;
+ }
+ } while (m <= modes);
+}
+
+void
+print_version
+ ()
+{
+ printf("stardis-app version %i.%i.%i built on stardis library version %i.%i.%i\n",
+ StardisApp_VERSION_MAJOR, StardisApp_VERSION_MINOR, StardisApp_VERSION_PATCH,
+ Stardis_VERSION_MAJOR, Stardis_VERSION_MINOR, Stardis_VERSION_PATCH);
+}
+
+void
+print_help
+ (char* prog)
+{
+ 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\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("\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("\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.vtk,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.vtk. 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 -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("\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("\n -g:\n");
+ printf(" Compute the green function.\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 -n NUM_OF_REALIZATIONS:\n");
+ printf(" Number of Monte-Carlo realizations.\n");
+ printf(" default value: 10000.\n");
+ printf("\n -t NUM_OF_THREADS:\n");
+ printf(" Number of threads; default value is all threads available.\n");
+ printf("\n -r AMBIENT_RAD_TEMP,REFERENCE_TEMP:\n");
+ printf(" Parameters for the radiative transfer.\n");
+ printf(" AMBIENT_RAD_TEMP is the ambient radiative temperature .\n");
+ printf(" REFERENCE_TEMP is the temperature used for the linearization\n");
+ printf(" of the radiative transfer.\n");
+}
+
+res_T
+parse_args
+ (const int argc,
+ char** argv,
+ struct args* args)
+{
+ int opt = 0;
+ size_t len = 0;
+ res_T res = RES_OK;
+
+ ASSERT(argv && args);
+
+ 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:B:M:m:p:P:dD:s:S:f:r:R:")) != -1) {
+ switch (opt) {
+ case '?': /* Unreconised option */
+ res = RES_BAD_ARG;
+ fprintf(stderr, "Invalid option -%c\n", optopt);
+ goto error;
+
+ case 'h':
+ print_help(argv[0]);
+ args->just_help = 1;
+ goto exit;
+
+ case 'v':
+ print_version();
+ args->just_help = 1;
+ goto exit;
+
+ case 'g':
+ args->mode |= GREEN_MODE;
+ break;
+
+ case 'n': {
+ unsigned long n;
+ res = cstr_to_ulong(optarg, &n);
+ if (res != RES_OK || n == 0) {
+ res = RES_BAD_ARG;
+ fprintf(stderr, "Invalid argument -%c %s\n", opt, optarg);
+ goto error;
+ }
+ args->N = n;
+ break;
+ }
+
+ case 'B':
+ args->bc_filename = optarg;
+ break;
+
+ case 'M':
+ args->medium_filename = optarg;
+ break;
+
+ case 't':
+ res = cstr_to_uint(optarg, &args->nthreads);
+ if (res != RES_OK
+ || args->nthreads <= 0) {
+ res = RES_BAD_ARG;
+ fprintf(stderr, "Invalid argument -%c %s\n", opt, optarg);
+ goto error;
+ }
+ break;
+
+ case 'p':
+ 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
+ || res != RES_OK
+ || args->u.probe[3] < 0)
+ {
+ res = RES_BAD_ARG;
+ fprintf(stderr, "Invalid argument -%c %s\n", opt, optarg);
+ goto error;
+ }
+ break;
+
+ case 'P':
+ 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
+ || res != RES_OK
+ || args->u.probe[3] < 0)
+ {
+ res = RES_BAD_ARG;
+ fprintf(stderr, "Invalid argument -%c %s\n", opt, optarg);
+ goto error;
+ }
+ break;
+
+ case 's':{
+ int err = 0;
+ char *ptr;
+ 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 |= BOUNDARY_COMPUTE;
+ ptr = strrchr(optarg, ',');
+ /* Single ',' allowed */
+ 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);
+ }
+ if (err) {
+ res = RES_BAD_ARG;
+ fprintf(stderr, "Invalid argument -%c %s\n", opt, optarg);
+ goto error;
+ }
+ break;
+ }
+
+ case 'm': {
+ char* ptr;
+ 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 |= MEDIUM_COMPUTE;
+ ptr = strpbrk(optarg, ",");
+ if (!ptr || ptr == optarg)
+ res = RES_BAD_ARG;
+ else res = cstr_to_double(ptr + 1, args->u.probe + 3);
+ if (res != RES_OK) {
+ fprintf(stderr, "Invalid argument -%c %s\n", opt, optarg);
+ goto error;
+ }
+ *ptr = '\0';
+ args->medium_name = optarg;
+ break;
+ }
+
+ 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")) {
+ args->dump_paths = DUMP_ALL;
+ }
+ else if (0 == strcmp(optarg, "error")) {
+ args->dump_paths = DUMP_ERROR;
+ }
+ else if (0 == strcmp(optarg, "success")) {
+ args->dump_paths = DUMP_SUCCESS;
+ }
+ else {
+ res = RES_BAD_ARG;
+ fprintf(stderr, "Invalid argument -%c %s\n", opt, optarg);
+ goto error;
+ }
+ break;
+
+ case 'f':
+ cstr_to_double(optarg, &args->scale_factor);
+ if (res != RES_OK
+ || args->scale_factor <= 0) {
+ res = RES_BAD_ARG;
+ fprintf(stderr, "Invalid argument -%c %s\n", opt, optarg);
+ goto error;
+ }
+ break;
+
+ case 'r':
+ cstr_to_list_double(optarg, ',', args->radiative_temp, &len, 2);
+ if (res != RES_OK
+ || len != 2
+ || args->radiative_temp[0] <= 0
+ || args->radiative_temp[1] <= 0) {
+ res = RES_BAD_ARG;
+ fprintf(stderr, "Invalid argument -%c %s\n", opt, optarg);
+ goto error;
+ }
+ break;
+
+ case 'R':
+ 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 |= IR_COMPUTE;
+ args->camera = optarg;
+ break;
+ }
+ }
+
+ if (args->mode & GREEN_MODE
+ && (args->mode != (args->mode & GREEN_COMPATIBLE_MODES)))
+ {
+ fprintf(stderr, "Option -g can only be used in conjunction with:");
+ print_multiple_modes(stderr, GREEN_COMPATIBLE_MODES);
+ fprintf(stderr, "\n");
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ if (!args->medium_filename) {
+ fprintf(stderr, "Missing mandatory argument: -M MEDIUM.txt\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;
+ }
+
+exit:
+ return res;
+error:
+ fprintf(stderr, "Use the -h option to get help.\n");
+ goto exit;
+}
+
+res_T
+parse_camera
+ (char* cam_param,
+ struct camera* cam)
+{
+ char** line = NULL;
+ char** opt = NULL;
+ res_T res = RES_OK;
+
+ ASSERT(cam_param && cam);
+
+ line = split_line(cam_param, ':');
+
+ if (line)
+ {
+ int i = 0;
+ for (i = 0; *(line + i); i++)
+ {
+ size_t len = 0;
+ opt = split_line(line[i], '=');
+ if (strcmp(opt[0], "t") == 0) {
+ res = cstr_to_double(opt[1], &cam->u.time);
+ if (res != RES_OK) goto error;
+ }
+ 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], "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], "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], "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], "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], "spp") == 0) {
+ res = cstr_to_uint(opt[1], &cam->spp);
+ if (res != RES_OK) goto error;
+ }
+ }
+ }
+
+exit:
+ FREE_AARRAY(line)
+ FREE_AARRAY(opt)
+
+ return res;
+error:
+ goto exit;
+}
+
+/* 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)
+{
+ char* tk = NULL;
+ 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));
+
+ CHK_TOK(strtok(NULL, " "), "file name");
+ *stl_filename = malloc(strlen(tk) + 1);
+ strcpy(*stl_filename, tk);
+
+ 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;
+ }
+ /* Depending of the number of spaces before '"' strtok should
+ * be called once or twice */
+ CHK_TOK(strtok(NULL, "\""), "Tinit");
+ if (*(tk + strspn(tk, " \t")) == '\0') {
+ CHK_TOK(strtok(NULL, "\""), "Tinit");
+ }
+ desc->d.solid.Tinit = malloc(strlen(tk) + 1);
+ strcpy(desc->d.solid.Tinit, tk);
+ /* Closing " fot Tinit; can return NULL if line ends just after "Tinit" */
+ tk = strtok(NULL, "\"");
+
+ /* Volumic Power is optional */
+ if (tk && *(tk + strspn(tk, " \t")) == '\0') {
+ /* Depending of the number of spaces before '"' strtok should
+ * be called once or twice */
+ tk = strtok(NULL, "\"");
+ }
+ if (tk) {
+ desc->d.solid.power = malloc(strlen(tk) + 1);
+ strcpy(desc->d.solid.power, tk);
+ /* Can be changed aftwerwards if compiled to constant 0 */
+ desc->d.solid.has_power = 1;
+ }
+ }
+ 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));
+
+ CHK_TOK(strtok(NULL, " "), "file name");
+ *stl_filename = malloc(strlen(tk) + 1);
+ strcpy(*stl_filename, tk);
+
+ 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");
+ if (*(tk + strspn(tk, " \t")) == '\0') {
+ CHK_TOK(strtok(NULL, "\""), "Tinit");
+ }
+ desc->d.fluid.Tinit = malloc(strlen(tk) + 1);
+ strcpy(desc->d.fluid.Tinit, tk);
+ }
+ else {
+ res = RES_BAD_ARG;
+ fprintf(stderr, "Invalid medium type: %s\n", tk);
+ }
+
+#undef CHK_TOK
+ exit :
+ return res;
+}
+
+/* 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)
+{
+ 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;
+
+ 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"));
+ }
+ }
+ }
+ }
+ }
+
+ 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;
+
+ 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;
+ }
+ /* Depending of the number of spaces before '"' strtok should
+ * be called once or twice */
+ CHK_TOK(strtok(NULL, "\""), "temperature");
+ if (*(tk + strspn(tk, " \t")) == '\0') {
+ CHK_TOK(strtok(NULL, "\""), "temperature");
+ }
+ desc->d.h_boundary.T = malloc(strlen(tk) + 1);
+ strcpy(desc->d.h_boundary.T, tk);
+ }
+ 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));
+
+ /* 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);
+
+ /* Depending of the number of spaces before '"' strtok should
+ * be called once or twice */
+ CHK_TOK(strtok(NULL, "\""), "temperature");
+ if (*(tk + strspn(tk, " \t")) == '\0') {
+ CHK_TOK(strtok(NULL, "\""), "temperature");
+ }
+ 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(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;
+ }
+ 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");
+ if (*(tk + strspn(tk, " \t")) == '\0') {
+ CHK_TOK(strtok(NULL, "\""), "flux");
+ }
+ desc->d.f_boundary.flux = malloc(strlen(tk) + 1);
+ strcpy(desc->d.f_boundary.flux, tk);
+ }
+ 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));
+
+ 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.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.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;
+ }
+ }
+ else {
+ fprintf(stderr, "Invalid boundary type: %s", tk);
+ res = RES_BAD_ARG;
+ goto exit;
+ }
+
+#undef CHK_TOK
+ exit :
+ return res;
+}
diff --git a/src/stardis-parsing.h b/src/stardis-parsing.h
@@ -0,0 +1,114 @@
+/* Copyright (C) 2018 |Meso|Star> (contact@meso-star.com)*/
+
+#ifndef STARDIS_PARSING_H
+#define STARDIS_PARSING_H
+
+#include <sdis.h>
+
+#include <rsys/rsys.h>
+
+struct camera;
+struct description;
+
+enum stardis_mode {
+ UNDEF_MODE = 0,
+ PROBE_COMPUTE = BIT(0),
+ PROBE_COMPUTE_ON_INTERFACE = BIT(1),
+ MEDIUM_COMPUTE = BIT(2),
+ BOUNDARY_COMPUTE = BIT(3),
+ IR_COMPUTE = BIT(4),
+ DUMP_VTK = BIT(6),
+ GREEN_MODE = BIT(7),
+
+ GREEN_COMPATIBLE_MODES = PROBE_COMPUTE | PROBE_COMPUTE_ON_INTERFACE | MEDIUM_COMPUTE
+ | BOUNDARY_COMPUTE,
+
+ COMPUTE_MODES = GREEN_COMPATIBLE_MODES | IR_COMPUTE,
+
+ NON_COMPUTE_MODES = UNDEF_MODE | DUMP_VTK | GREEN_MODE,
+
+ EXCLUSIVE_MODES = COMPUTE_MODES | DUMP_VTK
+};
+
+STATIC_ASSERT(GREEN_COMPATIBLE_MODES == (COMPUTE_MODES & GREEN_COMPATIBLE_MODES),
+ Cannot_have_a_GREEN_COMPATIBLE_MODE_that_is_not_a_COMPUTE_MODE);
+
+enum dump_path_type {
+ DUMP_NONE = 0,
+ DUMP_SUCCESS = BIT(0),
+ DUMP_ERROR = BIT(1),
+ DUMP_ALL = DUMP_SUCCESS | DUMP_ERROR,
+};
+
+struct args {
+ char* medium_filename;
+ char* bc_filename;
+ char* medium_name;
+ char* solve_boundary_filename;
+ size_t N;
+ unsigned nthreads;
+ union u { /* Trick to allow static initialization with INF */
+ struct pt { double p[3]; uint64_t t; } pt;
+ double probe[4];
+ } u;
+ enum stardis_mode mode;
+ double scale_factor;
+ double radiative_temp[2];
+ char* camera;
+ enum dump_path_type dump_paths;
+ int just_help;
+};
+
+#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 char
+mode_option
+ (const enum stardis_mode m);
+
+extern void
+print_multiple_modes
+ (FILE* stream,
+ const int modes);
+
+extern void
+print_version
+ ();
+
+extern void
+print_help
+ (char* prog);
+
+extern res_T
+parse_args
+ (const int argc,
+ char** argv,
+ struct args* args);
+
+extern res_T
+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
+ (char* line,
+ char** stl_filename,
+ struct description* desc);
+
+#endif /*ARGS_H*/
diff --git a/src/stardis-solid.c b/src/stardis-solid.c
@@ -0,0 +1,211 @@
+#include "stardis-solid.h"
+#include "stardis-compute.h"
+
+#include <sdis.h>
+#include <tinyexpr.h>
+#include<rsys/stretchy_array.h>
+
+/*******************************************************************************
+ * Local Functions
+ ******************************************************************************/
+
+static double
+solid_get_calorific_capacity
+ (const struct sdis_rwalk_vertex* vtx,
+ struct sdis_data* data)
+{
+ const struct solid* solid_props = sdis_data_cget(data);
+ (void)vtx;
+ return solid_props->cp;
+}
+
+static double
+solid_get_thermal_conductivity
+ (const struct sdis_rwalk_vertex* vtx,
+ struct sdis_data* data)
+{
+ const struct solid* solid_props = sdis_data_cget(data);
+ (void)vtx;
+ return solid_props->lambda;
+}
+
+static double
+solid_get_volumic_mass
+ (const struct sdis_rwalk_vertex* vtx,
+ struct sdis_data* data)
+{
+ const struct solid* solid_props = sdis_data_cget(data);
+ (void)vtx;
+ return solid_props->rho;
+}
+
+static double
+solid_get_delta
+ (const struct sdis_rwalk_vertex* vtx,
+ struct sdis_data* data)
+{
+ const struct solid* solid_props = sdis_data_cget(data);
+ (void)vtx;
+ return solid_props->delta;
+}
+
+#if Stardis_VERSION_MINOR == 3
+static double
+solid_get_delta_boundary
+ (const struct sdis_rwalk_vertex* vtx,
+ struct sdis_data* data)
+{
+ const struct solid* solid_props = sdis_data_cget(data);
+ return solid_props->delta
+ /* Emperical scale factor that ensures that delta_boundary > delta withouht
+ * being an exact multiple of delta. */
+ /** 2.1 */;
+}
+#endif
+
+static double
+solid_get_temperature
+ (const struct sdis_rwalk_vertex* vtx,
+ struct sdis_data* data)
+{
+ const struct solid* solid_props = sdis_data_cget(data);
+ char msg[128];
+ ASSERT(solid_props->t_init || solid_props->temp);
+ if (solid_props->is_green || vtx->time > solid_props->t0) {
+ /* Always use temp for Green mode, regardless of time */
+ if (solid_props->temp)
+ return te_eval(solid_props->temp, vtx);
+ else return -1;
+ }
+ /* Time is t0: use t_init */
+ if (solid_props->t_init)
+ return te_eval(solid_props->t_init, vtx);
+ /* Must have had t_init defined: error! */
+ if (solid_props->name[0])
+ sprintf(msg,
+ "solid_get_temperature: getting undefined Tinit (solid '%s')\n",
+ solid_props->name);
+ else
+ sprintf(msg, "solid_get_temperature: getting undefined Tinit\n");
+ FATAL(msg);
+}
+
+static double
+solid_get_power
+ (const struct sdis_rwalk_vertex* vtx,
+ struct sdis_data* data)
+{
+ const struct solid* solid_props = sdis_data_cget(data);
+ return te_eval(solid_props->power, vtx);
+}
+
+static void
+release_solid_data
+ (void* s)
+{
+ struct solid* solid = (struct solid*)s;
+ te_free(solid->t_init);
+ te_free(solid->power);
+}
+
+/*******************************************************************************
+ * Public Functions
+ ******************************************************************************/
+
+res_T
+create_solid
+ (struct sdis_device* dev,
+ const char* name,
+ const double lambda,
+ const double rho,
+ const double cp,
+ const double delta,
+ const int is_green,
+ const int is_outside,
+ const char* tinit_expr, /* Can be NULL if not used by getter */
+ const char* t_expr, /* Can be NULL if not used by getter */
+ const char* power_expr, /* Can be NULL */
+ struct sdis_medium*** media_ptr,
+ int* out_has_power, /* Can be NULL */
+ unsigned* out_id)
+{
+ res_T res = RES_OK;
+ struct sdis_solid_shader solid_shader = SDIS_SOLID_SHADER_NULL;
+ struct sdis_data* data = NULL;
+ struct solid* solid_props;
+ size_t sz;
+ /* Could be less restrictive if green output included positions/dates */
+ int prohibited = is_green ? NO_VAR_ALLOWED : ALL_VARS_ALLOWED;
+
+ ASSERT(dev && lambda >= 0 && rho >= 0 && cp >= 0 && delta > 0
+ && media_ptr && out_id);
+ 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;
+#if Stardis_VERSION_MINOR == 3
+ 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),
+ release_solid_data, &data);
+ if (res != RES_OK) goto error;
+
+ sz = sa_size(*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';
+ solid_props->lambda = lambda;
+ solid_props->rho = rho;
+ solid_props->cp = cp;
+ solid_props->delta = delta;
+ solid_props->t_init = NULL;
+ solid_props->temp = NULL;
+ solid_props->power = NULL;
+ solid_props->is_green = is_green;
+ solid_props->is_outside = is_outside;
+ solid_props->t0 = 0;
+ solid_props->id = (unsigned)sz;
+ if (tinit_expr) {
+ res = compile_expr_to_fn(&solid_props->t_init, tinit_expr,
+ prohibited | T_PROHIBITED, NULL);
+ if (res != RES_OK) {
+ if (res == RES_BAD_ARG)
+ fprintf(stderr, "Invalid initial temperature expression: %s\n",
+ tinit_expr);
+ goto error;
+ }
+ }
+ if (t_expr) {
+ res = compile_expr_to_fn(&solid_props->temp, t_expr, prohibited, NULL);
+ if (res != RES_OK) {
+ fprintf(stderr, "Invalid temperature expression: %s\n",
+ t_expr);
+ goto error;
+ }
+ }
+ if (power_expr) {
+ int has_power;
+ res = compile_expr_to_fn(&solid_props->power, power_expr, prohibited, &has_power);
+ if (res != RES_OK) {
+ fprintf(stderr, "Invalid volumic power expression: %s\n",
+ power_expr);
+ goto error;
+ }
+ 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;
+ *out_id = solid_props->id;
+
+end:
+ if (data) SDIS(data_ref_put(data));
+ return res;
+error:
+ goto end;
+}
diff --git a/src/stardis-solid.h b/src/stardis-solid.h
@@ -0,0 +1,49 @@
+/* Copyright (C) 2018 |Meso|Star> (contact@meso-star.com)*/
+
+#ifndef SDIS_SOLID_H
+#define SDIS_SOLID_H
+
+#include <rsys/rsys.h>
+
+struct te_expr;
+struct sdis_device;
+struct sdis_medium;
+
+/*******************************************************************************
+ * Solid data
+ ******************************************************************************/
+struct solid {
+ char name[32];
+ double cp; /* Calorific capacity */
+ double lambda; /* Conductivity */
+ double rho; /* Volumic mass */
+ double delta; /* Numerical parameter */
+ /* Compute mode */
+ int is_green, is_outside;
+ double t0;
+ /* TinyExpr stuff to compute temperature & power */
+ struct te_expr *temp;
+ struct te_expr *t_init;
+ struct te_expr *power;
+ /* ID */
+ unsigned id;
+};
+
+extern res_T
+create_solid
+ (struct sdis_device* dev,
+ const char* name,
+ const double lambda,
+ const double rho,
+ const double cp,
+ const double delta,
+ const int is_green,
+ const int is_outside,
+ const char* tinit_expr, /* Can be NULL if not used by getter */
+ const char* t_expr, /* Can be NULL if not used by getter */
+ const char* power_expr, /* Can be NULL */
+ struct sdis_medium*** media_ptr,
+ int* out_has_power, /* Can be NULL */
+ unsigned* out_id);
+
+#endif
diff --git a/src/stardis-version.h b/src/stardis-version.h
@@ -0,0 +1,14 @@
+/* Copying and distribution of this file, with or without modification,
+ * are permitted in any medium without royalty provided the copyright
+ * notice and this notice are preserved. This file is offered as-is,
+ * without any warranty. */
+
+#ifndef StardisApp_VERSION_H
+#define StardisApp_VERSION_H
+
+#define StardisApp_VERSION_MAJOR 0
+#define StardisApp_VERSION_MINOR 1
+#define StardisApp_VERSION_PATCH 0
+
+#endif /* StardisApp_VERSION_H */
+
