commit bd17c0e6867e6a84c0f3f40063d2d3bb52bf3a1c
parent 124f1ba33703524a312711f844859e4f012ac394
Author: Christophe Coustet <christophe.coustet@meso-star.com>
Date: Thu, 18 Jul 2019 10:25:29 +0200
Major refactoring: split huge files and functions
Diffstat:
| M | cmake/CMakeLists.txt | | | 12 | +++++++++++- |
| D | src/args.h | | | 373 | ------------------------------------------------------------------------------- |
| M | src/main.c | | | 24 | +++++++++++------------- |
| M | src/stardis-app.c | | | 603 | ++++++++----------------------------------------------------------------------- |
| M | src/stardis-app.h | | | 283 | +------------------------------------------------------------------------------ |
| M | src/stardis-compute.c | | | 2205 | ++++++++++++++++++------------------------------------------------------------- |
6 files changed, 572 insertions(+), 2928 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -57,11 +57,21 @@ set(VERSION ${VERSION_MAJOR}.${VERSION_MINOR}.${VERSION_PATCH})
set(SDIS_FILES_SRC
stardis-app.c
stardis-compute.c
+ stardis-fluid.c
+ stardis-intface.c
+ stardis-output.c
+ stardis-parsing.c
+ stardis-solid.c
main.c)
set(SDIS_FILES_INC
stardis-app.h
- args.h)
+ stardis-compute.h
+ stardis-fluid.h
+ stardis-intface.h
+ stardis-output.h
+ stardis-parsing.h
+ stardis-solid.h)
set(SDIS_FILES_DOC COPYING README.md)
diff --git a/src/args.h b/src/args.h
@@ -1,373 +0,0 @@
-/* Copyright (C) 2018 |Meso|Star> (contact@meso-star.com)*/
-
-#ifndef ARGS_H
-#define ARGS_H
-
-#include <getopt.h>
-#include <stdlib.h>
-
-#include <rsys/rsys.h>
-#include <rsys/cstr.h>
-#include <sdis_version.h>
-
-enum stardis_mode {
- UNDEF_MODE = 0,
- PROBE_COMPUTE = BIT(0),
- MEDIUM_COMPUTE = BIT(1),
- BOUNDARY_COMPUTE = BIT(2),
- GREEN_MODE = BIT(3),
- IR_COMPUTE = BIT(4),
- DUMP_VTK = BIT(5),
- PROBE_COMPUTE_ON_INTERFACE = BIT(6),
- COMPUTE_MODES = PROBE_COMPUTE | MEDIUM_COMPUTE | BOUNDARY_COMPUTE | IR_COMPUTE | DUMP_VTK | PROBE_COMPUTE_ON_INTERFACE
-};
-
-static char mode_option(const enum stardis_mode m) {
- int found = 0;
- char res = '?';
- if (m & PROBE_COMPUTE) { found++; res = 'p'; }
- if (m & MEDIUM_COMPUTE) { found++; res = 'm'; }
- if (m & BOUNDARY_COMPUTE) { found++; res = 'b'; }
- if (m & IR_COMPUTE) { found++; res = 'R'; }
- if (m & DUMP_VTK) { found++; res = 'd'; }
- if (m & PROBE_COMPUTE_ON_INTERFACE) { found++; res = 'P'; }
- ASSERT(found == 1);
- return res;
-}
-
-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__;
-
-static void
-print_help(char* prog)
-{
- printf("stardis-app built-on stardis library version %i.%i.%i\n",
- Stardis_VERSION_MAJOR,Stardis_VERSION_MINOR,Stardis_VERSION_PATCH);
- printf("usage: \n%s -M MEDIUM.txt -B BOUNDARY.txt\n", prog);
- printf("[ -p X,Y,Z,TIME | -m MEDIUM_NAME,TIME | -d | -b 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] [-s 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("\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 -d:\n");
- printf(" Dump the geometry in VTK format with medium front/back id and\n");
- printf(" boundary id.\n");
- printf("\n -b 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("\nOptionnal arguments\n");
- printf("-------------------\n");
- printf("\n -g:\n");
- printf(" Compute the green function.\n");
- printf(" Must be used in conjonction with one the -p, -m or -b options.\n");
- printf(" Provided TIME is silently ignored.\n");
- printf("\n -s 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");
-}
-
-
-static INLINE res_T
-parse_args(const int argc, char** argv, struct args* args)
-{
- int opt = 0;
- size_t len = 0;
- res_T res = RES_OK;
-
- if (argc == 1) {
- print_help(argv[0]);
- res = RES_BAD_ARG;
- goto error;
- }
-
- opterr = 0; /* No default error messages */
- while((opt = getopt(argc, argv, "hgn:t:b:B:M:m:p:P:dD:s: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 '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 & COMPUTE_MODES) {
- res = RES_BAD_ARG;
- fprintf(stderr, "Cannot specify multiple modes: -%c and -%c\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 & COMPUTE_MODES) {
- res = RES_BAD_ARG;
- fprintf(stderr, "Cannot specify multiple modes: -%c and -%c\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 'b': {
- int err = 0;
- char *ptr;
- if (args->mode & COMPUTE_MODES) {
- res = RES_BAD_ARG;
- fprintf(stderr, "Cannot specify multiple modes: -%c and -%c\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 & COMPUTE_MODES) {
- res = RES_BAD_ARG;
- fprintf(stderr, "Cannot specify multiple modes: -%c and -%c\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 & COMPUTE_MODES) {
- res = RES_BAD_ARG;
- fprintf(stderr, "Cannot specify multiple modes: -%c and -%c\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 's':
- 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 & COMPUTE_MODES) {
- res = RES_BAD_ARG;
- fprintf(stderr, "Cannot specify multiple modes: -%c and -%c\n",
- (char)opt, mode_option(args->mode));
- goto error;
- }
- args->mode |= IR_COMPUTE;
- args->camera = optarg;
- break;
- }
- }
-
-if (!args->medium_filename){
- fprintf(stderr, "Argument required -m\n");
- res = RES_BAD_ARG;
- goto error;
-}
-
-if (!args->bc_filename){
- fprintf(stderr, "Argument required -b\n");
- res = RES_BAD_ARG;
- goto error;
-}
-
-exit:
- return res;
-error:
- goto exit;
-}
-
-#endif /*ARGS_H*/
diff --git a/src/main.c b/src/main.c
@@ -1,12 +1,16 @@
/* Copyright (C) 2018 |Meso|Star> (contact@meso-star.com) */
#include "stardis-app.h"
+#include "stardis-output.h"
+
+#if !defined(ARCH_64BITS)
+#error "Please write the_version of hash functions needed by the current arch"
+#endif
int main(int argc, char** argv){
struct args args = ARGS_DEFAULT;
struct stardis stardis = NULL_STARDIS;
- size_t memsz = 0;
int err = 0;
res_T res = RES_OK;
@@ -25,21 +29,15 @@ int main(int argc, char** argv){
if (res != RES_OK) goto error;
goto exit;
}
-
- res = stardis_compute(&stardis, args.mode);
- if (res != RES_OK) goto error;
+ else if (args.mode & COMPUTE_MODES) {
+ res = stardis_compute(&stardis, args.mode);
+ if (res != RES_OK) goto error;
+ } else {
+ fprintf(stderr, "Unknown mode: %c.\n", mode_option(args.mode));
+ }
exit:
stardis_release(&stardis);
- if(stardis.allocator_initialized
- && (memsz = MEM_ALLOCATED_SIZE(&stardis.allocator)) != 0)
- {
- char dump[4096] = { '\0' };
- MEM_DUMP(&stardis.allocator, dump, sizeof(dump));
- fprintf(stderr, "%s\n", dump);
- fprintf(stderr, "Memory leaks: %lu Bytes\n", (unsigned long)memsz);
- }
- mem_shutdown_proxy_allocator(&stardis.allocator);
CHK(mem_allocated_size() == 0);
return err;
error:
diff --git a/src/stardis-app.c b/src/stardis-app.c
@@ -1,10 +1,7 @@
/* Copyright (C) 2018 |Meso|Star> (contact@meso-star.com) */
-#define _GNU_SOURCE 1
-#include <string.h>
-#include <ctype.h>
#include "stardis-app.h"
-
+#include "stardis-output.h"
#include <rsys/hash_table.h>
#define HTABLE_NAME descriptions
@@ -14,16 +11,25 @@
#define HTABLE_KEY_FUNCTOR_HASH hash_description
#include <rsys/hash_table.h>
+#define _GNU_SOURCE 1
+#include <string.h>
+#include <ctype.h>
+
/*******************************************************************************
- *
+ * Local Functions
******************************************************************************/
-static INLINE char*
-read_line(char** pb, FILE* stream)
+
+static char*
+read_line
+ (char** pb,
+ FILE* stream)
{
const int chunk_sz = 32;
char* buf = *pb;
size_t idx;
+ ASSERT(pb && stream);
+
if(!buf) buf = sa_add(buf, 128);
do {
@@ -48,444 +54,6 @@ read_line(char** pb, FILE* stream)
return buf;
}
-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;
-
- /* 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
-
-
-/* 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)"
-*/
-static res_T
-parse_medium_line(char* line, char** stl_filename, struct description* desc)
-{
- char* tk = NULL;
- res_T res = RES_OK;
-
-#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
- */
-static 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;
-
-#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;
-}
-
static res_T
read_vertices
(struct sstl_desc* desc,
@@ -496,6 +64,8 @@ read_vertices
res_T res = RES_OK;
unsigned vtx_index = 0;
+ ASSERT(desc && vertex2id && id2id && geometry);
+
for (vtx_index = 0; vtx_index < desc->vertices_count; ++vtx_index){
struct vertex vtx = NULL_VERTEX;
unsigned *found_id = NULL;
@@ -523,9 +93,14 @@ read_vertices
}
static int
-indices_to_key(struct unsigned3* key, const struct unsigned3* indices)
+indices_to_key
+ (struct unsigned3* key,
+ const struct unsigned3* indices)
{
int i, reversed = 0;
+
+ ASSERT(key && indices);
+
FOR_EACH(i, 0, 3) key->data[i] = indices->data[i];
if (key->data[0] > key->data[1]) {
reversed = !reversed;
@@ -557,6 +132,9 @@ read_triangles
unsigned* packed_indices = NULL;
float* packed_normals = NULL;
unsigned tri_index = 0, drop_cpt = 0;
+
+ ASSERT(stl_desc && triangle2id && stl_filename && id2id && geometry && boundary);
+
for (tri_index = 0; tri_index < stl_desc->triangles_count; ++tri_index) {
struct triangle tri = NULL_TRIANGLE;
const unsigned *found_id = NULL;
@@ -678,7 +256,7 @@ end:
goto end;
}
-res_T
+static res_T
read_stl
(const unsigned desc_id,
const enum content_type file_type,
@@ -692,6 +270,8 @@ read_stl
struct sstl_desc stl_desc;
unsigned* id2id = NULL;
+ ASSERT(stl_filename && vertex2id && triangle2id && stardis);
+
if (!stl_filename) return RES_BAD_ARG;
SSTL(create(NULL, &stardis->allocator, 0, &sstl));
@@ -717,10 +297,6 @@ error:
goto end;
}
-/*******************************************************************************
- *
- ******************************************************************************/
-
static res_T
geometry_analyse
(const char* medium_filename,
@@ -731,13 +307,17 @@ geometry_analyse
res_T res = RES_OK;
FILE* input = NULL;
char* line = NULL;
- unsigned sz = (unsigned)sa_size(stardis->descriptions);
+ unsigned sz;
struct htable_vertex vertex2id;
struct htable_triangle triangle2id;
struct htable_descriptions descriptions;
int tables_initialised = 0;
unsigned desc_id;
unsigned *p_desc;
+
+ ASSERT(medium_filename && bc_filename && stardis);
+
+ sz = (unsigned)sa_size(stardis->descriptions);
ASSERT(sz == sa_size(stardis->descriptions));
stardis->geometry = NULL_GEOMETRY;
@@ -846,59 +426,8 @@ error:
goto exit;
}
-static res_T
-parse_camera
- (char* cam_param, struct camera* cam)
-{
- char** line = NULL;
- char** opt = NULL;
- res_T res = RES_OK;
-
- 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;
-}
-
/*******************************************************************************
- *
+ * Public Functions
******************************************************************************/
res_T
@@ -908,9 +437,12 @@ stardis_init
{
res_T res = RES_OK;
+ ASSERT(args && stardis);
+
res = mem_init_proxy_allocator(&stardis->allocator,
&mem_default_allocator);
if (res != RES_OK) goto error;
+ stardis->allocator_initialized = 1;
res = geometry_analyse(args->medium_filename,
args->bc_filename, args->solve_boundary_filename, stardis);
@@ -951,8 +483,18 @@ error:
void
stardis_release(struct stardis* stardis)
{
- size_t i = 0;
+ size_t memsz, i = 0;
+
+ if(!stardis) return;
+ if (stardis->geometry.interfaces) {
+ for (i = 0; i < sa_size(stardis->geometry.interfaces); ++i)
+ SDIS(interface_ref_put(stardis->geometry.interfaces[i]));
+ sa_release(stardis->geometry.interfaces);
+ stardis->geometry.interfaces = NULL;
+ }
+ sa_release(stardis->geometry.interf_bytrg);
+ stardis->geometry.interf_bytrg = NULL;
for (i=0; i<sa_size(stardis->descriptions); ++i) {
struct description* desc = &stardis->descriptions[i];
switch (desc->type) {
@@ -985,48 +527,17 @@ stardis_release(struct stardis* stardis)
release_geometry(&stardis->geometry);
sa_release(stardis->boundary.primitives);
sa_release(stardis->boundary.sides);
-}
-
-/* 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);
+ if (stardis->allocator_initialized
+ && (memsz = MEM_ALLOCATED_SIZE(&stardis->allocator)) != 0)
+ {
+ char dump[4096] = { '\0' };
+ MEM_DUMP(&stardis->allocator, dump, sizeof(dump));
+ fprintf(stderr, "%s\n", dump);
+ fprintf(stderr, "Memory leaks: %lu Bytes\n", (unsigned long)memsz);
}
- 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);
+ if (stardis->allocator_initialized) {
+ mem_shutdown_proxy_allocator(&stardis->allocator);
+ stardis->allocator_initialized = 0;
}
-
-exit:
- return res;
-error:
- goto exit;
}
diff --git a/src/stardis-app.h b/src/stardis-app.h
@@ -3,6 +3,8 @@
#ifndef STARDIS_APP_H
#define STARDIS_APP_H
+#include "stardis-parsing.h"
+
#include <star/sstl.h>
#include <rsys/rsys.h>
#include <rsys/float3.h>
@@ -10,7 +12,6 @@
#include <rsys/hash_table.h>
#include <sdis.h>
#include <tinyexpr.h>
-#include "args.h"
#define FREE_AARRAY(ARRAY) if (ARRAY) {\
int i = 0; \
@@ -63,19 +64,6 @@ struct triangle {
#define NULL_TRIANGLE__ {{{0, 0, 0}}, UINT_MAX, UINT_MAX, UINT_MAX }
static const struct triangle NULL_TRIANGLE = NULL_TRIANGLE__;
-static INLINE 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));
-}
-
static INLINE char
eq_indices(const struct unsigned3* a, const struct unsigned3* b)
{
@@ -159,44 +147,6 @@ eq_fluid(const struct mat_fluid* a, const struct mat_fluid* b)
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->name)) {
- hash = hash ^ (uint32_t)((unsigned char)((const char*)&key->name)[i]);
- hash = hash * FNV32_PRIME;
- }
- 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;
- }
- hash = hash ^ (uint32_t)((unsigned char)(key->Tinit == 0));
- hash = hash * FNV32_PRIME;
- FOR_EACH(i, 0, key->Tinit ? strlen(key->Tinit) * sizeof(*key->Tinit) : 0) {
- hash = hash ^ (uint32_t)((unsigned char)((const char*)key->Tinit)[i]);
- hash = hash * FNV32_PRIME;
- }
- hash = hash ^ (uint32_t)((unsigned char)(key->Temp == 0));
- hash = hash * FNV32_PRIME;
- FOR_EACH(i, 0, key->Temp ? strlen(key->Temp) * sizeof(*key->Temp) : 0) {
- hash = hash ^ (uint32_t)((unsigned char)((const char*)key->Temp)[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);
@@ -233,9 +183,6 @@ hash_fluid(const struct mat_fluid* key)
hash = hash * FNV64_PRIME;
}
return hash;
-#else
-#error "Unexpected architecture"
-#endif
}
struct mat_solid {
@@ -286,61 +233,6 @@ eq_solid(const struct mat_solid* a, const struct mat_solid* b)
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->name)) {
- hash = hash ^ (uint32_t)((unsigned char)((const char*)&key->name)[i]);
- hash = hash * FNV32_PRIME;
- }
- 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;
- }
- hash = hash ^ (uint32_t)((unsigned char)(key->Tinit == 0));
- hash = hash * FNV32_PRIME;
- FOR_EACH(i, 0, key->Tinit ? strlen(key->Tinit) * sizeof(*key->Tinit) : 0) {
- hash = hash ^ (uint32_t)((unsigned char)((const char*)key->Tinit)[i]);
- hash = hash * FNV32_PRIME;
- }
- hash = hash ^ (uint32_t)((unsigned char)(key->Temp == 0));
- hash = hash * FNV32_PRIME;
- FOR_EACH(i, 0, key->Temp ? strlen(key->Temp) * sizeof(*key->Temp) : 0) {
- hash = hash ^ (uint32_t)((unsigned char)((const char*)key->Temp)[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);
@@ -394,9 +286,6 @@ hash_solid(const struct mat_solid* key)
hash = hash * FNV64_PRIME;
}
return hash;
-#else
-#error "Unexpected architecture"
-#endif
}
struct h_boundary {
@@ -451,53 +340,6 @@ eq_h_boundary
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->name)) {
- hash = hash ^ (uint32_t)((unsigned char)((const char*)&key->name)[i]);
- hash = hash * FNV32_PRIME;
- }
- 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);
@@ -543,9 +385,6 @@ hash_h_boundary(const struct h_boundary* key)
hash = hash * FNV64_PRIME;
}
return hash;
-#else
-#error "Unexpected architecture"
-#endif
}
struct t_boundary {
@@ -606,41 +445,6 @@ hash_t_boundary
{
(void)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->name)) {
- hash = hash ^ (uint32_t)((unsigned char)((const char*)&key->name)[i]);
- hash = hash * FNV32_PRIME;
- }
- 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);
@@ -674,9 +478,6 @@ hash_t_boundary
hash = hash * FNV64_PRIME;
}
return hash;
-#else
-#error "Unexpected architecture"
-#endif
}
struct f_boundary {
@@ -712,28 +513,6 @@ eq_f_boundary(const struct f_boundary* a, const struct f_boundary* b)
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->name)) {
- hash = hash ^ (uint32_t)((unsigned char)((const char*)&key->name)[i]);
- hash = hash * FNV32_PRIME;
- }
- 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, strlen(key->flux) * sizeof(*key->flux)) {
- hash = hash ^ (uint32_t)((unsigned char)((const char*)key->flux)[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);
@@ -754,9 +533,6 @@ hash_f_boundary(const struct f_boundary* key)
hash = hash * FNV64_PRIME;
}
return hash;
-#else
-#error "Unexpected architecture"
-#endif
}
struct solid_fluid_connect {
@@ -785,44 +561,6 @@ eq_sf_connect(const struct solid_fluid_connect* a, const struct solid_fluid_conn
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->name)) {
- hash = hash ^ (uint32_t)((unsigned char)((const char*)&key->name)[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, 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);
@@ -859,9 +597,6 @@ hash_sf_connect(const struct solid_fluid_connect* key)
hash = hash * FNV64_PRIME;
}
return hash;
-#else
-#error "Unexpected architecture"
-#endif
}
struct description {
@@ -885,7 +620,7 @@ init_description(struct description* desc) {
static INLINE void
print_description
(FILE* stream,
- struct description* desc)
+ const struct description* desc)
{
ASSERT(stream && desc);
switch (desc->type) {
@@ -1023,16 +758,4 @@ stardis_compute(struct stardis* stardis, enum stardis_mode mode);
extern void
stardis_release(struct stardis* stardis);
-extern res_T
-dump_vtk(FILE* output, const struct geometry* geometry);
-
-extern res_T
-read_stl
- (const unsigned desc_id,
- const enum content_type file_type,
- const char* stl_filename,
- struct htable_vertex* vertex2id,
- struct htable_triangle* triangle2id,
- struct stardis* stardis);
-
#endif /*STARDIS-APP_H*/
diff --git a/src/stardis-compute.c b/src/stardis-compute.c
@@ -1,38 +1,33 @@
/* Copyright (C) 2018 |Meso|Star> (contact@meso-star.com)*/
#include "stardis-app.h"
+#include "stardis-output.h"
+#include "stardis-compute.h"
+#include "stardis-fluid.h"
+#include "stardis-solid.h"
+#include "stardis-intface.h"
+
#include <sdis.h>
#include <sdis_version.h>
#include <rsys/double3.h>
#include <rsys/double2.h>
-#include <rsys/dynamic_array_uint.h>
#include <rsys/image.h>
-#include<star/senc.h>
-
-#include <rsys/hash_table.h>
-#define HTABLE_NAME vrtx_rank
-#define HTABLE_DATA unsigned
-#define HTABLE_KEY unsigned
-#include <rsys/hash_table.h>
-#define HTABLE_NAME weigth
-#define HTABLE_DATA double
-#define HTABLE_KEY unsigned
-#include <rsys/hash_table.h>
-
-struct w_ctx {
- struct description* desc;
- struct htable_weigth weigths;
+struct dummies {
+ struct sdis_medium* dummy_fluid;
+ unsigned dummy_fluid_id;
+ struct sdis_medium* dummy_solid;
+ unsigned dummy_solid_id;
};
+#define DUMMIES_NULL__ {\
+ NULL, UINT_MAX, NULL, UINT_MAX\
+}
+static const struct dummies DUMMIES_NULL = DUMMIES_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
-};
+/*******************************************************************************
+ * Local Functions
+ ******************************************************************************/
static void
geometry_get_position
@@ -41,6 +36,7 @@ geometry_get_position
void* context)
{
struct geometry* geom = context;
+ ASSERT(geom);
pos[0] = geom->vertex[ivert].xyz[0];
pos[1] = geom->vertex[ivert].xyz[1];
pos[2] = geom->vertex[ivert].xyz[2];
@@ -53,6 +49,7 @@ geometry_get_indices
void* context)
{
struct geometry* geom = context;
+ ASSERT(geom);
ids[0] = geom->triangle[itri].indices.data[0];
ids[1] = geom->triangle[itri].indices.data[1];
ids[2] = geom->triangle[itri].indices.data[2];
@@ -65,305 +62,16 @@ geometry_get_interface
void* context)
{
struct geometry* geom = context;
+ ASSERT(geom);
*interf = geom->interf_bytrg[itri];
}
-static res_T
-compile_expr_to_fn
- (struct te_expr** f,
- const char* math_expr,
- const int prohibited,
- 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))
- };
- int fst = 0, lst = 4;
- ASSERT(math_expr);
- ASSERT(prohibited == (prohibited & NO_VAR_ALLOWED)); /* Use only defined flags */
- if (prohibited & XYZ_PROHIBITED) fst = 3;
- if(prohibited & T_PROHIBITED) lst = 3;
- *f = te_compile(math_expr, vars+fst, lst-fst, NULL);
- if (!*f) {
- if ((prohibited != ALL_VARS_ALLOWED)
- && te_compile(math_expr, vars, 4, NULL))
- fprintf(stderr, "Expression %s use a probibited variable ", math_expr);
- if (prohibited == NO_VAR_ALLOWED)
- fprintf(stderr, "(no variable allowed)\n");
- else if (prohibited == XYZ_PROHIBITED)
- fprintf(stderr, "(xyz prohibited)\n");
- else {
- ASSERT(prohibited == T_PROHIBITED);
- fprintf(stderr, "(t prohibited)\n");
- }
- return RES_BAD_ARG;
- }
- if (is_zero) *is_zero = !((*f)->type == TE_CONSTANT && (*f)->v.value == 0);
- return RES_OK;
-}
-
-/*******************************************************************************
- * 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;
-};
-
-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);
-}
-
-/*******************************************************************************
- * 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;
-};
-
-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);
-}
-
-/*******************************************************************************
- * 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;
-};
-
-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 void
-release_interface_data(void* s)
-{
- struct intface* intface = (struct intface*)s;
- te_free(intface->temperature);
- te_free(intface->flux);
-}
-
-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 INLINE 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;
-}
-
+/* Check probe position and move it to the closest point on an interface
+ * if on_interface is set */
static res_T
select_probe_type
(const struct sdis_scene* scn,
- const struct triangle* triangle,
- const struct vertex* vertex,
- const struct description* descriptions,
+ const struct stardis* stardis,
const int on_interface,
double pos[3],
size_t* iprim,
@@ -374,7 +82,14 @@ select_probe_type
const struct description *min_desc = NULL;
size_t i = 0, found = SIZE_MAX;
enum description_type min_type = DESCRIPTION_TYPE_COUNT__;
-
+ const struct triangle* triangle;
+ const struct vertex* vertex;
+ const struct description* descriptions;
+
+ ASSERT(scn && stardis && pos && iprim && uv);
+ triangle = stardis->geometry.triangle;
+ vertex = stardis->geometry.vertex;
+ descriptions = stardis->descriptions;
for (i = 0; i < sa_size(triangle); ++i) {
const struct triangle* tr = triangle + i;
double dp[3], d, projected_pos[3], nn[3];
@@ -437,7 +152,7 @@ select_probe_type
}
if (!min_desc) {
- /* Don't manage pos is outside the model
+ /* Doesn't handle the case where pos is outside the model
* but very close and could be considered on the boundary */
fprintf(stderr, "The probe is outside the model.\n");
return RES_BAD_ARG;
@@ -469,1114 +184,551 @@ select_probe_type
return res;
}
-static void
-dump_image(const struct sdis_estimator_buffer* buf)
-{
- size_t definition[2];
- double* temps = NULL;
- size_t ix, iy;
-
- CHK(buf != NULL);
- CHK(sdis_estimator_buffer_get_definition(buf, definition) == RES_OK);
-
- temps = mem_alloc(definition[0] * definition[1] *sizeof(double));
- CHK(temps != NULL);
-
- /* 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;
- CHK(sdis_estimator_buffer_at(buf, ix, iy, &estimator) == RES_OK);
- CHK(sdis_estimator_get_temperature(estimator, &T) == RES_OK);
- 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;
- CHK(sdis_estimator_buffer_at(buf, ix, iy, &estimator) == RES_OK);
- CHK(sdis_estimator_get_temperature(estimator, &T) == RES_OK);
- 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;
- CHK(sdis_estimator_buffer_at(buf, ix, iy, &estimator) == RES_OK);
- CHK(sdis_estimator_get_realisation_time(estimator, &time) == RES_OK);
- 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;
- CHK(sdis_estimator_buffer_at(buf, ix, iy, &estimator) == RES_OK);
- CHK(sdis_estimator_get_realisation_time(estimator, &time) == RES_OK);
- 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;
- CHK(sdis_estimator_buffer_at(buf, ix, iy, &estimator) == RES_OK);
- CHK(sdis_estimator_get_failure_count(estimator, &nfails) == RES_OK);
- fprintf(stdout, "%zu\n", nfails);
- }
- }
- mem_rm(temps);
-}
-
-
-/*******************************************************************************
- *
- ******************************************************************************/
-
-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)
-{
-#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
-}
-
-#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>
-
static 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)
+create_holder
+ (struct counts* counts,
+ struct dummies* dummies,
+ struct description* description,
+ struct sdis_device* dev,
+ struct sdis_medium*** media,
+ const int is_green)
{
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);
+ switch (description->type) {
+ case DESC_BOUND_H_FOR_SOLID:
+ counts->hbound_count++;
+ /* Create an external fluid with bound temp as fluid temp */
+ res = create_fluid(dev, NULL, 1, 1, is_green, 1, NULL,
+ description->d.h_boundary.T, media, &description->d.h_boundary.mat_id);
+ if (res != RES_OK) goto error;
+ break;
+ case DESC_BOUND_H_FOR_FLUID:
+ /* Create an external solid with bound temp as solid temp */
+ counts->hbound_count++;
+ res = create_solid(dev, NULL, INF, 1, 1, 1, is_green, 1, NULL,
+ description->d.h_boundary.T, NULL, media, NULL,
+ &description->d.h_boundary.mat_id);
+ if (res != RES_OK) goto error;
+ break;
+ case DESC_BOUND_T_FOR_SOLID:
+ counts->tbound_count++;
+ ASSERT(description->d.t_boundary.T != NULL);
+ res = compile_expr_to_fn(&description->d.t_boundary.te_temperature,
+ description->d.t_boundary.T, prohibited, NULL);
if (res != RES_OK) {
- fprintf(stderr, "Invalid initial temperature expression: %s\n",
- tinit_expr);
+ fprintf(stderr, "Invalid boundary temperature expression: %s\n",
+ description->d.t_boundary.T);
goto error;
}
- }
- if (t_expr) {
- res = compile_expr_to_fn(&fluid_props->temp, t_expr, prohibited, NULL);
+ if (dummies->dummy_fluid) {
+ /* Reuse external dummy fluid */
+ description->d.t_boundary.mat_id = dummies->dummy_fluid_id;
+ }
+ else {
+ /* Create dummy fluid */
+ res = create_fluid(dev, NULL, 1, 1, is_green, 1, NULL,
+ NULL, media, &description->d.t_boundary.mat_id);
+ if (res != RES_OK) goto error;
+ dummies->dummy_fluid = sa_last(*media);
+ dummies->dummy_fluid_id = description->d.t_boundary.mat_id;
+ }
+ break;
+ case DESC_BOUND_T_FOR_FLUID:
+ counts->tbound_count++;
+ ASSERT(description->d.t_boundary.T != NULL);
+ res = compile_expr_to_fn(&description->d.t_boundary.te_temperature,
+ description->d.t_boundary.T, prohibited, NULL);
if (res != RES_OK) {
- fprintf(stderr, "Invalid temperature expression: %s\n",
- t_expr);
+ fprintf(stderr, "Invalid boundary temperature expression: %s\n",
+ description->d.t_boundary.T);
goto error;
}
+ if (dummies->dummy_solid) {
+ /* Reuse external dummy solid */
+ description->d.t_boundary.mat_id = dummies->dummy_solid_id;
+ }
+ else {
+ /* Create dummy solid */
+ res = create_solid(dev, NULL, 1, 1, 1, 1, is_green, 1, NULL,
+ NULL, NULL, media, NULL, &description->d.t_boundary.mat_id);
+ if (res != RES_OK) goto error;
+ dummies->dummy_solid = sa_last(*media);
+ dummies->dummy_solid_id = description->d.t_boundary.mat_id;
+ }
+ break;
+ case DESC_BOUND_F_FOR_SOLID:
+ counts->fbound_count++;
+ ASSERT(description->d.f_boundary.flux != NULL);
+ res = compile_expr_to_fn(&description->d.f_boundary.te_flux,
+ description->d.f_boundary.flux, 1, NULL);
+ if (res != RES_OK) {
+ fprintf(stderr, "Invalid boundary flux expression: %s\n",
+ description->d.f_boundary.flux);
+ goto error;
+ }
+ if (dummies->dummy_fluid) {
+ /* Reuse external dummy fluid */
+ description->d.f_boundary.mat_id = dummies->dummy_fluid_id;
+ }
+ else {
+ /* Create dummy fluid */
+ res = create_fluid(dev, NULL, 1, 1, is_green, 1, NULL, NULL,
+ media, &description->d.f_boundary.mat_id);
+ if (res != RES_OK) goto error;
+ dummies->dummy_fluid = sa_last(*media);
+ dummies->dummy_fluid_id = description->d.f_boundary.mat_id;
+ }
+ break;
+ case DESC_SOLID_FLUID_CONNECT:
+ counts->sfconnect_count++;
+ ASSERT(description->d.sf_connect.specular_fraction >= 0
+ && description->d.sf_connect.specular_fraction <= 1);
+ ASSERT(description->d.sf_connect.emissivity >= 0);
+ ASSERT(description->d.sf_connect.hc >= 0);
+ break;
+ case DESC_MAT_SOLID:
+ counts->smed_count++;
+ res = create_solid(dev,
+ description->d.solid.name,
+ description->d.solid.lambda,
+ description->d.solid.rho,
+ description->d.solid.cp,
+ description->d.solid.delta,
+ is_green,
+ 0,
+ description->d.solid.Tinit,
+ description->d.solid.Temp,
+ description->d.solid.power,
+ media,
+ &description->d.solid.has_power,
+ &description->d.solid.solid_id);
+ if (res != RES_OK) goto error;
+ break;
+ case DESC_MAT_FLUID:
+ counts->fmed_count++;
+ res = create_fluid(dev,
+ description->d.fluid.name,
+ description->d.fluid.rho,
+ description->d.fluid.cp,
+ is_green,
+ 0,
+ description->d.fluid.Tinit,
+ description->d.fluid.Temp,
+ media,
+ &description->d.fluid.fluid_id);
+ if (res != RES_OK) goto error;
+
+ break;
+ default: FATAL("Invalid type.\n");
}
- 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;
}
static 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)
+compute_probe
+ (struct sdis_scene* scn,
+ const struct counts* counts,
+ struct stardis* stardis,
+ enum stardis_mode mode)
{
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;
+ double time[2], pos[3], uv[2] = { 0,0 };
+ size_t iprim = SIZE_MAX;
+ struct sdis_green_function* green = NULL;
+ struct sdis_estimator* estimator = NULL;
+
+ ASSERT(scn && counts && stardis && (mode & PROBE_COMPUTE));
- 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);
+ d3_set(pos, stardis->probe);
+ res = select_probe_type(scn, stardis, 0, pos, &iprim, uv);
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) {
- 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;
+ if (mode & GREEN_MODE) {
+ ASSERT(iprim == SIZE_MAX);
+ res = sdis_solve_probe_green_function(scn,
+ stardis->N,
+ pos,
+ stardis->scale_factor,
+ stardis->radiative_temp[0],
+ stardis->radiative_temp[1],
+ &green);
+ if (res != RES_OK) goto error;
+ res = dump_green(green, counts, stardis->descriptions,
+ stardis->radiative_temp);
+ if (res != RES_OK) goto error;
+ } else {
+ ASSERT(iprim == SIZE_MAX);
+ d2_splat(time, stardis->probe[3]);
+ res = sdis_solve_probe(scn,
+ stardis->N,
+ pos,
+ time,
+ stardis->scale_factor,
+ stardis->radiative_temp[0],
+ stardis->radiative_temp[1],
+ stardis->dump_paths,
+ &estimator);
+ if (res != RES_OK) goto error;
+ res = print_single_MC_result(mode, estimator, stardis);
+ if (res != RES_OK) goto error;
}
- 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));
+ if (estimator) SDIS(estimator_ref_put(estimator));
+ if (green) SDIS(green_function_ref_put(green));
return res;
error:
goto end;
}
static res_T
-create_edge_file_if
+compute_probe_on_interface
(struct sdis_scene* scn,
- struct mem_allocator* allocator)
+ const struct counts* counts,
+ struct stardis* stardis,
+ enum stardis_mode mode)
{
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);
+ double time[2], pos[3], uv[2] = { 0,0 };
+ size_t iprim = SIZE_MAX;
+ struct sdis_estimator* estimator = NULL;
+ struct sdis_green_function* green = NULL;
- res = senc_descriptor_get_frontier_segments_count(analyze, &scount);
+ ASSERT(scn && counts && stardis && (mode & PROBE_COMPUTE_ON_INTERFACE));
+ d3_set(pos, stardis->probe);
+ res = select_probe_type(scn, stardis, 1, pos, &iprim, uv);
if (res != RES_OK) goto error;
- 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;
+ if (mode & GREEN_MODE) {
+ ASSERT(iprim == SIZE_MAX);
+ res = sdis_solve_probe_boundary_green_function(scn,
+ stardis->N,
+ iprim,
+ uv,
+ SDIS_FRONT,
+ stardis->scale_factor,
+ stardis->radiative_temp[0],
+ stardis->radiative_temp[1],
+ &green);
+ if (res != RES_OK) goto error;
+ res = dump_green(green, counts, stardis->descriptions,
+ stardis->radiative_temp);
+ if (res != RES_OK) goto error;
+ } else {
+ ASSERT(iprim == SIZE_MAX);
+ d2_splat(time, stardis->probe[3]);
+ res = sdis_solve_probe_boundary(scn,
+ stardis->N,
+ iprim,
+ uv,
+ time,
+ SDIS_FRONT, /* FIXME!!!!!!! */
+ stardis->scale_factor,
+ stardis->radiative_temp[0],
+ stardis->radiative_temp[1],
+ stardis->dump_paths,
+ &estimator);
+ if (res != RES_OK) goto error;
+ res = print_single_MC_result(mode, estimator, stardis);
+ if (res != RES_OK) goto error;
}
-
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);
- }
+ if (estimator) SDIS(estimator_ref_put(estimator));
+ if (green) SDIS(green_function_ref_put(green));
return res;
error:
goto end;
}
static res_T
-dump_path
- (const struct sdis_heat_path* path,
- void* context)
+compute_camera
+ (struct sdis_device* dev,
+ struct sdis_scene* scn,
+ const struct stardis* stardis,
+ const enum stardis_mode mode)
{
- FILE* stream = context;
- enum sdis_heat_path_flag status = SDIS_HEAT_PATH_NONE;
- size_t i, vcount;
+ res_T res = RES_OK;
+ double time[2];
+ size_t width, height;
+ struct sdis_estimator_buffer* buf = NULL;
+ struct sdis_camera* cam = NULL;
- /* 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;
- CHK(sdis_heat_path_get_vertex(path, i, &vtx) == RES_OK);
- 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;
- CHK(sdis_heat_path_get_vertex(path, i, &vtx) == RES_OK);
- 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;
- CHK(sdis_heat_path_get_vertex(path, i, &vtx) == RES_OK);
- 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;
- CHK(sdis_heat_path_get_vertex(path, i, &vtx) == RES_OK);
- 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");
- CHK(sdis_heat_path_get_status(path, &status) == RES_OK);
- 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 */
- return RES_OK;
-}
+ ASSERT(!(mode & GREEN_MODE) && (mode & IR_COMPUTE) && scn && stardis);
+ (void)mode;
-static res_T
-print_power_term
- (struct sdis_medium* mdm,
- const double power_term,
- void* ctx)
-{
- struct sdis_data* data = NULL;
- enum sdis_medium_type type;
- unsigned id;
- struct w_ctx* w_ctx = ctx;
- (void)w_ctx;
- CHK(mdm && ctx);
+ width = (size_t)stardis->camera.img[0];
+ height = (size_t)stardis->camera.img[1];
+ /* Setup the camera */
+ res = sdis_camera_create(dev, &cam);
+ if (res != RES_OK) goto error;
+ res = sdis_camera_set_proj_ratio(cam, (double)width / (double)height);
+ if (res != RES_OK) goto error;
+ res = sdis_camera_set_fov(cam, MDEG2RAD(stardis->camera.fov));
+ if (res != RES_OK) goto error;
+ res = sdis_camera_look_at(cam,
+ stardis->camera.pos,
+ stardis->camera.tgt,
+ stardis->camera.up);
+ if (res != RES_OK) goto error;
- data = sdis_medium_get_data(mdm);
- type = sdis_medium_get_type(mdm);
+ /* Launch the simulation */
+ time[0] = time[1] = stardis->camera.u.time;
+ res = sdis_solve_camera(scn,
+ cam,
+ time,
+ stardis->scale_factor,
+ stardis->radiative_temp[0],
+ stardis->radiative_temp[1],
+ width,
+ height,
+ (size_t)stardis->camera.spp,
+ stardis->dump_paths,
+ &buf);
+ if (res != RES_OK) goto error;
- switch (type) {
- case SDIS_FLUID: {
- FATAL("Unexpected power term in fluid");
-#if 0
- struct fluid* d__ = sdis_data_get(data);
- id = d__->id;
- desc = (struct description*)ctx + id;
- ASSERT(id == desc->d.fluid.fluid_id);
- printf("\tF\t%u\t%g", id, power_term);
- break;
-#endif
- }
- case SDIS_SOLID: {
- struct solid* d__ = sdis_data_get(data);
- id = d__->id;
- ASSERT(id == w_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_OK;
+ /* Write the image */
+ res = dump_image(buf);
+ if (res != RES_OK) goto error;
+
+end:
+ if (cam) SDIS(camera_ref_put(cam));
+ if (buf) SDIS(estimator_buffer_ref_put(buf));
+ return res;
+error:
+ goto end;
}
static res_T
-get_flux_terms
- (struct sdis_interface* interf,
- const enum sdis_side side,
- const double flux_term,
- void* ctx)
+compute_medium
+ (struct sdis_scene* scn,
+ struct sdis_medium** media,
+ struct stardis* stardis,
+ enum stardis_mode mode)
{
- struct sdis_data* data = NULL;
- struct intface* d__;
- unsigned id;
- struct w_ctx* w_ctx = ctx;
res_T res = RES_OK;
+ double time[2];
+ size_t sz, ii;
+ struct sdis_medium* medium = NULL;
+ struct sdis_estimator* estimator = NULL;
+ struct sdis_green_function* green = NULL;
- CHK(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;
+ ASSERT(scn && media && stardis && (mode & MEDIUM_COMPUTE));
- 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);
+ /* Find medium */
+ sz = sa_size(stardis->descriptions);
+ FOR_EACH(ii, 0, sz) {
+ struct description* desc = stardis->descriptions + ii;
+ if (desc->type == DESC_MAT_SOLID) {
+ if (0 == strcmp(stardis->solve_name, desc->d.solid.name)) {
+ ASSERT(sa_size(media) > ii);
+ medium = media[ii];
+ break;
+ }
}
- break;
+ else if (desc->type == DESC_MAT_FLUID) {
+ if (0 == strcmp(stardis->solve_name, desc->d.fluid.name)) {
+ ASSERT(sa_size(media) > ii);
+ medium = media[ii];
+ break;
+ }
+ }
+ }
+ if (medium == NULL) { /* Not found */
+ fprintf(stderr, "Cannot solve medium %s (unknown medium)\n",
+ stardis->solve_name);
+ res = RES_BAD_ARG;
+ goto error;
}
- default: FATAL("Unreachable code.\n"); break;
+
+ if (mode & GREEN_MODE) {
+ res = sdis_solve_medium_green_function(scn,
+ stardis->N,
+ medium,
+ stardis->scale_factor,
+ stardis->radiative_temp[0],
+ stardis->radiative_temp[1],
+ &green);
+ if (res != RES_OK) goto error;
+ } else {
+ d2_splat(time, stardis->probe[3]);
+ time[0] = time[1] = stardis->probe[3];
+ res = sdis_solve_medium(scn,
+ stardis->N,
+ medium,
+ time,
+ stardis->scale_factor,
+ stardis->radiative_temp[0],
+ stardis->radiative_temp[1],
+ stardis->dump_paths,
+ &estimator);
+ if (res != RES_OK) goto error;
+ res = print_single_MC_result(mode, estimator, stardis);
+ if (res != RES_OK) goto error;
}
+end:
+ if (estimator) SDIS(estimator_ref_put(estimator));
+ if (green) SDIS(green_function_ref_put(green));
return res;
+error:
+ goto end;
}
static res_T
-print_sample
- (struct sdis_green_path* path,
- void* ctx)
+compute_boundary
+ (struct sdis_scene* scn,
+ const struct counts* counts,
+ struct stardis* stardis,
+ enum stardis_mode mode)
{
- 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;
res_T res = RES_OK;
- CHK(path && ctx);
-
- res = sdis_green_path_get_limit_point(path, &pt);
- if (res != RES_OK) goto error;
-
+ double time[2], pos[3], uv[2] = { 0,0 };
+ size_t iprim = SIZE_MAX;
+ struct sdis_green_function* green = NULL;
+ struct sdis_estimator* estimator = NULL;
- /* 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
- */
+ ASSERT(scn && counts && stardis && (mode & BOUNDARY_COMPUTE));
- 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;
+ /* FIXME: MOVE TO BOUNDARY SETTING */
+ ASSERT(sa_size(stardis->boundary.primitives)
+ == sa_size(stardis->boundary.sides));
+ if (sa_size(stardis->boundary.sides) == 0) {
+ fprintf(stderr, "Cannot solve boundary %s (empty geometry)\n",
+ stardis->solve_name);
}
- 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);
+ d3_set(pos, stardis->probe);
+ res = select_probe_type(scn, stardis, 0, pos, &iprim, uv);
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);
+ if (mode & GREEN_MODE) {
+ ASSERT(iprim == SIZE_MAX);
+ res = sdis_solve_boundary_green_function(scn,
+ stardis->N,
+ stardis->boundary.primitives,
+ stardis->boundary.sides,
+ sa_size(stardis->boundary.sides),
+ stardis->scale_factor,
+ stardis->radiative_temp[0],
+ stardis->radiative_temp[1],
+ &green);
+ if (res != RES_OK) goto error;
+ res = dump_green(green, counts, stardis->descriptions,
+ stardis->radiative_temp);
+ if (res != RES_OK) goto error;
+ } else {
+ ASSERT(iprim == SIZE_MAX);
+ d2_splat(time, stardis->probe[3]);
+ res = sdis_solve_boundary(scn,
+ stardis->N,
+ stardis->boundary.primitives,
+ stardis->boundary.sides,
+ sa_size(stardis->boundary.sides),
+ time,
+ stardis->scale_factor,
+ stardis->radiative_temp[0],
+ stardis->radiative_temp[1],
+ stardis->dump_paths,
+ &estimator);
+ if (res != RES_OK) goto error;
+ res = print_single_MC_result(mode, estimator, stardis);
+ if (res != RES_OK) goto error;
+ }
+end:
+ if (estimator) SDIS(estimator_ref_put(estimator));
+ if (green) SDIS(green_function_ref_put(green));
+ return res;
+error:
+ goto end;
+}
- printf("\t%zu\t%zu", pcount, fcount);
+/*******************************************************************************
+ * Public Functions
+ ******************************************************************************/
- res = sdis_green_path_for_each_power_term(path, print_power_term, ctx);
- if (res != RES_OK) goto error;
+res_T
+compile_expr_to_fn
+ (struct te_expr** f,
+ const char* math_expr,
+ const int prohibited,
+ 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))
+ };
+ int fst = 0, lst = 4;
- 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");
+ ASSERT(f && math_expr);
+ ASSERT(prohibited == (prohibited & NO_VAR_ALLOWED)); /* Use only defined flags */
-error:
- return res;
+ if (prohibited & XYZ_PROHIBITED) fst = 3;
+ if (prohibited & T_PROHIBITED) lst = 3;
+ *f = te_compile(math_expr, vars + fst, lst - fst, NULL);
+ if (!*f) {
+ if (!te_compile(math_expr, vars, 4, NULL))
+ return RES_BAD_ARG; /* Expr is not valid */
+ /* Expr is valid */
+ if (prohibited != ALL_VARS_ALLOWED)
+ fprintf(stderr, "Expression \"%s\" use a probibited variable ", math_expr);
+ if (prohibited == NO_VAR_ALLOWED)
+ fprintf(stderr, "(no variable allowed)\n");
+ else if (prohibited == XYZ_PROHIBITED)
+ fprintf(stderr, "(xyz prohibited)\n");
+ else {
+ ASSERT(prohibited == T_PROHIBITED);
+ fprintf(stderr, "(t prohibited)\n");
+ }
+ return RES_BAD_OP; /* Expr is valid, but use a forbidden var */
+ }
+ if (is_zero) *is_zero = !((*f)->type == TE_CONSTANT && (*f)->v.value == 0);
+ return RES_OK;
}
res_T
-stardis_compute(struct stardis* stardis, enum stardis_mode mode)
+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_data* data = NULL;
-
struct sdis_scene* scn = NULL;
- struct sdis_estimator* estimator = NULL;
- struct sdis_mc temperature;
- struct sdis_estimator_buffer* buf = NULL;
- struct sdis_camera* cam = NULL;
-
- struct sdis_medium* dummy_fluid = NULL;
- unsigned dummy_fluid_id = UINT_MAX;
- struct sdis_medium* dummy_solid = NULL;
- unsigned dummy_solid_id = UINT_MAX;
-
+ struct dummies dummies = DUMMIES_NULL;
struct htable_intface htable_interfaces;
int htable_interfaces_initialised = 0;
- double pos[3] = { 0,0,0 };
- double time[2] = { 0, 0 };
- size_t nfailures;
- size_t smed_count = 0, fmed_count = 0,
- tbound_count = 0, hbound_count = 0, fbound_count = 0, sfconnect_count = 0;
+ struct counts counts = COUNTS_NULL;
unsigned i = 0;
+ ASSERT(stardis);
+
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;
- /* Could be less restrictive if green output included positions/dates */
- int prohibited = (mode & GREEN_MODE) ? NO_VAR_ALLOWED : ALL_VARS_ALLOWED;
-
- switch (desc->type) {
- case DESC_BOUND_H_FOR_SOLID:
- hbound_count++;
- /* Create an external fluid with bound temp as fluid temp */
- res = create_fluid(dev, NULL, 1, 1, (mode & GREEN_MODE), 1, NULL,
- desc->d.h_boundary.T, &media, &desc->d.h_boundary.mat_id);
- if (res != RES_OK) goto error;
- break;
- case DESC_BOUND_H_FOR_FLUID:
- /* Create an external solid with bound temp as solid temp */
- hbound_count++;
- res = create_solid(dev, NULL, INF, 1, 1, 1, (mode & GREEN_MODE), 1, NULL,
- desc->d.h_boundary.T, NULL, &media, NULL, &desc->d.h_boundary.mat_id);
- if (res != RES_OK) goto error;
- break;
- case DESC_BOUND_T_FOR_SOLID:
- tbound_count++;
- ASSERT(desc->d.t_boundary.T != NULL);
- res = compile_expr_to_fn(&desc->d.t_boundary.te_temperature,
- desc->d.t_boundary.T, prohibited, NULL);
- if (res != RES_OK) {
- fprintf(stderr, "Invalid boundary temperature expression: %s\n",
- desc->d.t_boundary.T);
- goto error;
- }
- if (dummy_fluid) {
- /* Reuse external dummy fluid */
- desc->d.t_boundary.mat_id = dummy_fluid_id;
- } else {
- /* Create dummy fluid */
- res = create_fluid(dev, NULL, 1, 1, (mode & GREEN_MODE), 1, NULL,
- NULL, &media, &desc->d.t_boundary.mat_id);
- if (res != RES_OK) goto error;
- dummy_fluid = sa_last(media);
- dummy_fluid_id = desc->d.t_boundary.mat_id;
- }
- break;
- case DESC_BOUND_T_FOR_FLUID:
- tbound_count++;
- ASSERT(desc->d.t_boundary.T != NULL);
- res = compile_expr_to_fn(&desc->d.t_boundary.te_temperature,
- desc->d.t_boundary.T, prohibited, NULL);
- if (res != RES_OK) {
- fprintf(stderr, "Invalid boundary temperature expression: %s\n",
- desc->d.t_boundary.T);
- goto error;
- }
- if (dummy_solid) {
- /* Reuse external dummy solid */
- desc->d.t_boundary.mat_id = dummy_solid_id;
- } else {
- /* Create dummy solid */
- res = create_solid(dev, NULL, 1, 1, 1, 1, (mode & GREEN_MODE), 1, NULL,
- NULL, NULL, &media, NULL, &desc->d.t_boundary.mat_id);
- if (res != RES_OK) goto error;
- dummy_solid = sa_last(media);
- dummy_solid_id = desc->d.t_boundary.mat_id;
- }
- break;
- case DESC_BOUND_F_FOR_SOLID:
- fbound_count++;
- 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.f_boundary.flux);
- goto error;
- }
- if (dummy_fluid) {
- /* Reuse external dummy fluid */
- desc->d.f_boundary.mat_id = dummy_fluid_id;
- } else {
- /* Create dummy fluid */
- res = create_fluid(dev, NULL, 1, 1, (mode & GREEN_MODE), 1, NULL, NULL,
- &media, &desc->d.f_boundary.mat_id);
- if (res != RES_OK) goto error;
- dummy_fluid = sa_last(media);
- dummy_fluid_id = desc->d.f_boundary.mat_id;
- }
- break;
- case DESC_SOLID_FLUID_CONNECT:
- sfconnect_count++;
- 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:
- smed_count++;
- res = create_solid(dev,
- stardis->descriptions[i].d.solid.name,
- stardis->descriptions[i].d.solid.lambda,
- stardis->descriptions[i].d.solid.rho,
- stardis->descriptions[i].d.solid.cp,
- stardis->descriptions[i].d.solid.delta,
- (mode & GREEN_MODE),
- 0,
- stardis->descriptions[i].d.solid.Tinit,
- stardis->descriptions[i].d.solid.Temp,
- stardis->descriptions[i].d.solid.power,
- &media,
- &stardis->descriptions[i].d.solid.has_power,
- &desc->d.solid.solid_id);
- if (res != RES_OK) goto error;
- break;
- case DESC_MAT_FLUID:
- fmed_count++;
- res = create_fluid(dev,
- stardis->descriptions[i].d.fluid.name,
- stardis->descriptions[i].d.fluid.rho,
- stardis->descriptions[i].d.fluid.cp,
- (mode & GREEN_MODE),
- 0,
- stardis->descriptions[i].d.fluid.Tinit,
- stardis->descriptions[i].d.fluid.Temp,
- &media,
- &desc->d.fluid.fluid_id);
- if (res != RES_OK) goto error;
-
- break;
- default: FATAL("Invalid type.\n");
- }
+ res = create_holder(&counts, &dummies, stardis->descriptions + i,
+ dev, &media, mode & GREEN_MODE);
+ if (res != RES_OK) goto error;
}
/* Create interfaces */
htable_intface_init(&stardis->allocator, &htable_interfaces);
htable_interfaces_initialised = 1;
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_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 (stardis->descriptions[fd].type) {
- case DESC_MAT_SOLID:
- id = stardis->descriptions[fd].d.solid.solid_id;
- solid_count++;
- front_med = media[id];
- break;
- case DESC_MAT_FLUID:
- fluid_count++;
- id = stardis->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 (stardis->descriptions[bd].type) {
- case DESC_MAT_SOLID:
- id = stardis->descriptions[bd].d.solid.solid_id;
- solid_count++;
- back_med = media[id];
- break;
- case DESC_MAT_FLUID:
- fluid_count++;
- id = stardis->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 = stardis->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", i);
- print_trg_as_obj(stderr, stardis->geometry.vertex, trg->indices.data);
- fprintf(stderr, "Front: ");
- if (!front_defined) fprintf(stderr, "undefined\n");
- else print_description(stderr, &stardis->descriptions[fd]);
- fprintf(stderr, "Back: ");
- if (!back_defined) fprintf(stderr, "undefined\n");
- else print_description(stderr, &stardis->descriptions[bd]);
- fprintf(stderr, "Connection: ");
- if (!connect_defined) fprintf(stderr, "undefined\n");
- else print_description(stderr, &stardis->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", i);
- goto error;
- }
- sa_push(interfaces, intface);
- res = htable_intface_set(&htable_interfaces, &int_descs, &intface);
- if (res != RES_OK) goto error;
- }
- sa_push(interf_bytrg, intface);
+ res = create_intface(dev, i, &htable_interfaces, &stardis->geometry,
+ stardis->descriptions, media);
+ if (res != RES_OK) goto error;
}
- stardis->geometry.interf_bytrg = interf_bytrg;
- stardis->geometry.interfaces = interfaces;
+ /* Create scene */
res = sdis_scene_create(dev,
sa_size(stardis->geometry.triangle),
geometry_get_indices, geometry_get_interface,
@@ -1584,408 +736,31 @@ stardis_compute(struct stardis* stardis, enum stardis_mode mode)
geometry_get_position, &stardis->geometry, &scn);
if (res != RES_OK) goto error;
+ /* If the scene has holes dump them */
res = create_edge_file_if(scn, &stardis->allocator);
if (res != RES_OK) goto error;
- if (mode & GREEN_MODE) {
- struct sdis_green_function* green = NULL;
- size_t ok_count, failed_count;
- struct w_ctx w_ctx;
- ASSERT((mode & PROBE_COMPUTE) || (mode & PROBE_COMPUTE_ON_INTERFACE)
- || (mode & BOUNDARY_COMPUTE) || (mode & MEDIUM_COMPUTE));
-
- if (mode & PROBE_COMPUTE || mode & PROBE_COMPUTE_ON_INTERFACE) {
- double uv[2] = { 0,0 };
- size_t iprim = SIZE_MAX;
- /* Launch the probe simulation */
- pos[0] = stardis->probe[0];
- pos[1] = stardis->probe[1];
- pos[2] = stardis->probe[2];
-
- res = select_probe_type(scn,
- stardis->geometry.triangle,
- stardis->geometry.vertex,
- stardis->descriptions,
- (mode & PROBE_COMPUTE_ON_INTERFACE),
- pos,
- &iprim,
- uv);
- if (res != RES_OK) goto error;
-
- if (iprim == SIZE_MAX) {
- res = sdis_solve_probe_green_function(scn,
- stardis->N,
- pos,
- stardis->scale_factor,
- stardis->radiative_temp[0],
- stardis->radiative_temp[1],
- &green);
- } else {
- res = sdis_solve_probe_boundary_green_function(scn,
- stardis->N,
- iprim,
- uv,
- SDIS_FRONT,
- stardis->scale_factor,
- stardis->radiative_temp[0],
- stardis->radiative_temp[1],
- &green);
- }
- }
- else if (mode & BOUNDARY_COMPUTE) {
- ASSERT(sa_size(stardis->boundary.primitives)
- == sa_size(stardis->boundary.sides));
- if (sa_size(stardis->boundary.sides) == 0) {
- fprintf(stderr, "Cannot solve boundary %s (empty geometry)\n",
- stardis->solve_name);
- }
-
- res = sdis_solve_boundary_green_function(scn,
- stardis->N,
- stardis->boundary.primitives,
- stardis->boundary.sides,
- sa_size(stardis->boundary.sides),
- stardis->scale_factor,
- stardis->radiative_temp[0],
- stardis->radiative_temp[1],
- &green);
- }
- else if (mode & MEDIUM_COMPUTE) {
- size_t sz, ii;
- struct sdis_medium* medium = NULL;
- /* Find medium */
- sz = sa_size(stardis->descriptions);
- FOR_EACH(ii, 0, sz) {
- struct description* desc = stardis->descriptions + ii;
- if (desc->type == DESC_MAT_SOLID) {
- if (0 == strcmp(stardis->solve_name, desc->d.solid.name)) {
- ASSERT(sa_size(media) > ii);
- medium = media[ii];
- break;
- }
- }
- else if (desc->type == DESC_MAT_FLUID) {
- if (0 == strcmp(stardis->solve_name, desc->d.fluid.name)) {
- ASSERT(sa_size(media) > ii);
- medium = media[ii];
- break;
- }
- }
- }
- if (medium == NULL) {
- /* Not found */
- fprintf(stderr, "Cannot solve medium %s (unknown medium)\n",
- stardis->solve_name);
- res = RES_BAD_ARG;
- goto error;
- }
- res = sdis_solve_medium_green_function(scn,
- stardis->N,
- medium,
- stardis->scale_factor,
- stardis->radiative_temp[0],
- stardis->radiative_temp[1],
- &green);
- }
- if (res != RES_OK) goto error;
-
- /* Output the green function */
- CHK(sdis_green_function_get_invalid_paths_count(green, &failed_count) == RES_OK);
- CHK(sdis_green_function_get_paths_count(green, &ok_count) == RES_OK);
-
- /* 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",
- smed_count, fmed_count,
- tbound_count, hbound_count, fbound_count,
- ok_count, failed_count);
-
- /* List Media */
- if (smed_count) {
- printf("# Solids\n");
- printf("# ID Name lambda rho cp power\n");
- FOR_EACH(i, 0, (unsigned)sa_size(stardis->descriptions)) {
- struct description* desc = stardis->descriptions + i;
- 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 (fmed_count) {
- printf("# Fluids\n");
- printf("# ID Name rho cp\n");
- FOR_EACH(i, 0, sa_size(stardis->descriptions)) {
- struct description* desc = stardis->descriptions + i;
- 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 (tbound_count) {
- printf("# T Boundaries\n");
- printf("# ID Name temperature\n");
- FOR_EACH(i, 0, sa_size(stardis->descriptions)) {
- struct description* desc = stardis->descriptions + i;
- 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 (hbound_count) {
- printf("# H Boundaries\n");
- printf("# ID Name emissivity specular_fraction hc hc_max T_env\n");
- FOR_EACH(i, 0, sa_size(stardis->descriptions)) {
- struct description* desc = stardis->descriptions + i;
- 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 (fbound_count) {
- printf("# F Boundaries\n");
- printf("# ID Name flux\n");
- FOR_EACH(i, 0, sa_size(stardis->descriptions)) {
- struct description* desc = stardis->descriptions + i;
- 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(stardis->descriptions),
- stardis->radiative_temp[0], stardis->radiative_temp[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 = stardis->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);
- } else {
- if (mode & IR_COMPUTE) {
- size_t width = (size_t)stardis->camera.img[0];
- size_t height = (size_t)stardis->camera.img[1];
-
- /* Setup the camera */
- SDIS(camera_create(dev, &cam));
- SDIS(camera_set_proj_ratio(cam, (double)width / (double)height));
- SDIS(camera_set_fov(cam, MDEG2RAD(stardis->camera.fov)));
- SDIS(camera_look_at(cam,
- stardis->camera.pos,
- stardis->camera.tgt,
- stardis->camera.up));
-
- /* Launch the simulation */
- time[0] = time[1] = stardis->camera.u.time;
- res = sdis_solve_camera(scn,
- cam,
- time,
- stardis->scale_factor,
- stardis->radiative_temp[0],
- stardis->radiative_temp[1],
- width,
- height,
- (size_t)stardis->camera.spp,
- stardis->dump_paths,
- &buf);
- if (res != RES_OK) goto error;
-
- /* Write the image */
- dump_image(buf);
- }
- else if (mode & PROBE_COMPUTE || mode & PROBE_COMPUTE_ON_INTERFACE) {
- double uv[2] = { 0,0 };
- size_t iprim = SIZE_MAX;
- /* Launch the probe simulation */
- pos[0] = stardis->probe[0];
- pos[1] = stardis->probe[1];
- pos[2] = stardis->probe[2];
- time[0] = time[1] = stardis->probe[3];
-
- res = select_probe_type(scn,
- stardis->geometry.triangle,
- stardis->geometry.vertex,
- stardis->descriptions,
- (mode & PROBE_COMPUTE_ON_INTERFACE),
- pos,
- &iprim,
- uv);
- if (res != RES_OK) goto error;
-
- if (iprim == SIZE_MAX) {
- res = sdis_solve_probe(scn,
- stardis->N,
- pos,
- time,
- stardis->scale_factor,
- stardis->radiative_temp[0],
- stardis->radiative_temp[1],
- stardis->dump_paths,
- &estimator);
- } else {
- res = sdis_solve_probe_boundary(scn,
- stardis->N,
- iprim,
- uv,
- time,
- SDIS_FRONT,
- stardis->scale_factor,
- stardis->radiative_temp[0],
- stardis->radiative_temp[1],
- stardis->dump_paths,
- &estimator);
- }
- }
- else if (mode & MEDIUM_COMPUTE) {
- size_t sz, ii;
- struct sdis_medium* medium = NULL;
- /* Find medium */
- sz = sa_size(stardis->descriptions);
- FOR_EACH(ii, 0, sz) {
- struct description* desc = stardis->descriptions + ii;
- if (desc->type == DESC_MAT_SOLID) {
- if (0 == strcmp(stardis->solve_name, desc->d.solid.name)) {
- ASSERT(sa_size(media) > ii);
- medium = media[ii];
- break;
- }
- }
- else if (desc->type == DESC_MAT_FLUID) {
- if (0 == strcmp(stardis->solve_name, desc->d.fluid.name)) {
- ASSERT(sa_size(media) > ii);
- medium = media[ii];
- break;
- }
- }
- }
- if (medium == NULL) {
- /* Not found */
- fprintf(stderr, "Cannot solve medium %s (unknown medium)\n",
- stardis->solve_name);
- res = RES_BAD_ARG;
- goto error;
- }
- time[0] = time[1] = stardis->probe[3];
- res = sdis_solve_medium(scn,
- stardis->N,
- medium,
- time,
- stardis->scale_factor,
- stardis->radiative_temp[0],
- stardis->radiative_temp[1],
- stardis->dump_paths,
- &estimator);
- }
- else if (mode & BOUNDARY_COMPUTE) {
- ASSERT(sa_size(stardis->boundary.primitives)
- == sa_size(stardis->boundary.sides));
- if (sa_size(stardis->boundary.sides) == 0) {
- fprintf(stderr, "Cannot solve boundary %s (empty geometry)\n",
- stardis->solve_name);
- }
-
- time[0] = time[1] = stardis->probe[3];
- res = sdis_solve_boundary(scn,
- stardis->N,
- stardis->boundary.primitives,
- stardis->boundary.sides,
- sa_size(stardis->boundary.sides),
- time,
- stardis->scale_factor,
- stardis->radiative_temp[0],
- stardis->radiative_temp[1],
- stardis->dump_paths,
- &estimator);
- } else {
- CHK(0);
- }
- if (res != RES_OK) goto error;
- if (mode & PROBE_COMPUTE
- || mode & PROBE_COMPUTE_ON_INTERFACE
- || mode & MEDIUM_COMPUTE
- || mode & BOUNDARY_COMPUTE) {
- /* Fetch the estimation data */
- SDIS(estimator_get_temperature(estimator, &temperature));
- SDIS(estimator_get_failure_count(estimator, &nfailures));
-
- /* Print the results */
- switch (mode & COMPUTE_MODES) {
- case PROBE_COMPUTE:
- case PROBE_COMPUTE_ON_INTERFACE:
- printf("Temperature at t=[%g, %g, %g, %g] = %g +/- %g\n",
- pos[0], pos[1], pos[2], time[0],
- 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, time[0],
- 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, time[0],
- 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);
-
- /* Dump paths according to user settings */
- sdis_estimator_for_each_path(estimator, dump_path, stdout);
- }
- }
+ /* Compute */
+ if (mode & PROBE_COMPUTE)
+ res = compute_probe(scn, &counts, stardis, mode);
+ else if (mode & PROBE_COMPUTE_ON_INTERFACE)
+ res = compute_probe_on_interface(scn, &counts, stardis, mode);
+ else if (mode & IR_COMPUTE)
+ res = compute_camera(dev, scn, stardis, mode);
+ else if (mode & MEDIUM_COMPUTE)
+ res = compute_medium(scn, media, stardis, mode);
+ else if (mode & BOUNDARY_COMPUTE)
+ res = compute_boundary(scn, &counts, stardis, mode);
+ else FATAL("Unknown mode.\n");
+ if (res != RES_OK) goto error;
end:
- if (data) SDIS(data_ref_put(data));
- if (interfaces) {
- for (i = 0; i < sa_size(interfaces); ++i)
- SDIS(interface_ref_put(interfaces[i]));
- sa_release(interfaces);
- stardis->geometry.interfaces = NULL;
- }
- sa_release(interf_bytrg);
- stardis->geometry.interf_bytrg = NULL;
- if (htable_interfaces_initialised) htable_intface_release(&htable_interfaces);
+ if (htable_interfaces_initialised)
+ htable_intface_release(&htable_interfaces);
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(estimator_buffer_ref_put(buf));
- if (estimator) SDIS(estimator_ref_put(estimator));
if (scn) SDIS(scene_ref_put(scn));
if (dev) SDIS(device_ref_put(dev));
