commit 6df6aa046f4e153d16df6d0fb83056d5710c9723
parent fd8791d4293cb259b045f101581e41e2d6969231
Author: Christophe Coustet <christophe.coustet@meso-star.com>
Date: Thu, 11 Apr 2019 17:34:11 +0200
Add solve_boundary capability and a time arg to the rendering mode
Diffstat:
5 files changed, 210 insertions(+), 68 deletions(-)
diff --git a/src/args.h b/src/args.h
@@ -10,16 +10,17 @@
#include <rsys/cstr.h>
#include <sdis_version.h>
-enum stardis_mode{
+enum stardis_mode {
PROBE_COMPUTE,
MEDIUM_COMPUTE,
+ BOUNDARY_COMPUTE,
IR_COMPUTE,
DUMP_VTK,
UNDEF_MODE
};
static const char MODE_OPTION[UNDEF_MODE] = {
- 'p', 'm', 'R', 'd'
+ 'p', 'm', 'b', 'R', 'd'
};
enum dump_path_type {
@@ -33,9 +34,10 @@ struct args{
char* medium_filename;
char* bc_filename;
char* medium_name;
+ char* solve_boundary_filename;
size_t N;
unsigned nthreads;
- union u {
+ union u { /* Trick to allow static initialization with INF */
struct pt { double p[3]; uint64_t t; } pt;
double probe[4];
} u;
@@ -52,7 +54,7 @@ struct args{
#define DEFAULT_NTHREADS 1
#endif
#define ARGS_DEFAULT__ {\
- NULL, NULL, NULL, 10000, DEFAULT_NTHREADS,\
+ 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__;
@@ -62,56 +64,64 @@ 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 \n");
- printf(" | -R spp = SPP:fov = FOV : up = XUP, YUP, ZUP : pos = X, Y, Z : tgt = XT, YT, ZT : img = WxH]\n");
- printf("[-s SCALE_FACTOR] [-D {all | error | success}]\n");
+ 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("[-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("\nArguments\n");
+ printf("\n -h: print this help.\n");
+ printf("\nMandatory arguments\n");
printf("---------\n");
- printf("\n -M MEDIUM.txt :\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("\n -B BOUNDARY.txt:\n");
printf(" BOUNDARY.txt is a text file which contains the description of the\n");
printf(" boundary conditions.\n");
- printf("\nOptionnal arguments\n");
+ printf("\nExclusive arguments\n");
printf("-------------------\n");
- printf("\n -p X:Y:Z:TIME :\n");
- printf(" this is the position and the time of the probe.\n");
- printf(" default value: 0.0,0.0,0.0,INF.\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. If some boundary conditions are\n");
printf(" time-dependant, TIME can not be INF.\n");
- printf("\n -m medium_name,TIME :\n");
- printf(" this is the medium name and the time of the by-medium integration.\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 can not be INF.\n");
- printf("\n -s SCALE_FACTOR :\n");
- printf(" default value: 1.0.\n");
- printf("\n -d :\n");
- printf(" dump geometry in VTK format with medium front/back id and\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 -D { all | error | success } :\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 can not 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 -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("\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("\n -t NUM_OF_THREADS:\n");
printf(" default value: all threads available.\n");
- printf("\n -r AMBIENT_RAD_TEMP,REFERENCE_TEMP :\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 theradiative transfer.\n");
- printf("\n -R 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(" 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");
}
@@ -128,7 +138,7 @@ parse_args(const int argc, char** argv, struct args* args)
goto error;
}
- while((opt = getopt(argc, argv, "hn:t:B:M:m:p:dD:s:r:R:")) != -1) {
+ while((opt = getopt(argc, argv, "hn:t:b:B:M:m:p:dD:s:r:R:")) != -1) {
switch(opt) {
case 'h':
print_help(argv[0]);
@@ -180,6 +190,33 @@ parse_args(const int argc, char** argv, struct args* args)
}
break;
+ case 'b': {
+ int err = 0;
+ char *ptr;
+ if (args->mode != UNDEF_MODE) {
+ 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 != strchr(optarg, ','))
+ err = 1;
+ else {
+ *ptr = '\0';
+ ptr++;
+ args->solve_boundary_filename = optarg;
+ err = RES_OK != cstr_to_double(ptr, args->u.probe + 3);
+ }
+ if (err) {
+ res = RES_BAD_ARG;
+ fprintf(stderr, "Invalid argument -b %s\n", optarg);
+ goto error;
+ }
+ break;
+ }
case 'm': {
char* ptr;
if (args->mode != UNDEF_MODE) {
diff --git a/src/main.c b/src/main.c
@@ -20,7 +20,7 @@ int main(int argc, char** argv){
res = stardis_init(&args, &stardis);
if (res != RES_OK) goto error;
- if (args.mode == DUMP_VTK){
+ if (args.mode == DUMP_VTK) {
res = dump_vtk(stdout, &stardis.geometry);
if (res != RES_OK) goto error;
goto exit;
diff --git a/src/stardis-app.c b/src/stardis-app.c
@@ -526,9 +526,10 @@ read_triangles
struct htable_triangle* triangle2id,
const char* stl_filename,
const unsigned* id2id,
- const int is_connection,
+ const enum content_type file_type,
const unsigned desc_id,
- struct geometry* geometry)
+ struct geometry* geometry,
+ struct solve_boundary* boundary)
{
res_T res = RES_OK;
unsigned* packed_indices = NULL;
@@ -577,6 +578,14 @@ read_triangles
ASSERT(id < sa_size(geometry->triangle));
} else {
size_t sz;
+ if (file_type == CONTENT_BOUNDARY) {
+ /* Cannot create new geometry! */
+ fprintf(stderr,
+ "Triangle %u from file %s cannot compute temperature on this surface\n"
+ "Not part of the model\n",
+ tri_index, stl_filename);
+ return RES_BAD_ARG;
+ }
sz = htable_triangle_size_get(triangle2id);
ASSERT(sz == sa_size(geometry->triangle));
ASSERT(sz < UINT_MAX);
@@ -586,20 +595,29 @@ read_triangles
sa_push(geometry->triangle, tri);
}
- if (is_connection) {
+ switch (file_type) {
+ case CONTENT_CONNECTION:
side_desc_ptr = &geometry->triangle[id].connection_description;
- }
- else {
+ break;
+ case CONTENT_GEOMETRY:
/* Is the triangle in triangle list reversed wrt stl_desc? */
side_desc_ptr = (current_side == SDIS_FRONT) ?
&geometry->triangle[id].front_description
: &geometry->triangle[id].back_description;
+ break;
+ case CONTENT_BOUNDARY:
+ sa_push(boundary->primitives, id);
+ sa_push(boundary->sides, current_side);
+ continue; /* Next triangle */
+ break;
+ default: FATAL("Invalid type.\n");
}
if (*side_desc_ptr != UINT_MAX && *side_desc_ptr != desc_id) {
/* Already described with a different description! */
fprintf(stderr, "Triangle %u %s with 2 different descriptions\n",
tri_index,
- is_connection ? "connection" : (current_side == SDIS_FRONT ? "front" : "back"));
+ (file_type == CONTENT_CONNECTION
+ ? "connection" : (current_side == SDIS_FRONT ? "front" : "back")));
print_trg_as_obj(stderr, (const struct vertex*)stl_desc->vertices, &tri);
return RES_BAD_ARG;
}
@@ -628,11 +646,11 @@ read_triangles
return res;
}
-static res_T
+res_T
read_stl
(const unsigned desc_id,
- const int is_connection,
- char* const* stl_filename,
+ const enum content_type file_type,
+ const char* stl_filename,
struct htable_vertex* vertex2id,
struct htable_triangle* triangle2id,
struct stardis* stardis)
@@ -646,17 +664,17 @@ read_stl
SSTL(create(NULL, &stardis->allocator, 0, &sstl));
- res = sstl_load(sstl, *stl_filename);
+ res = sstl_load(sstl, stl_filename);
if (res != RES_OK) {
- fprintf(stderr, "Cannot read STL file: %s\n", *stl_filename);
+ fprintf(stderr, "Cannot read STL file: %s\n", stl_filename);
goto error;
}
SSTL(get_desc(sstl, &stl_desc));
res = read_vertices(&stl_desc, vertex2id, &id2id, &stardis->geometry);
if (res != RES_OK) goto error;
- res = read_triangles(&stl_desc, triangle2id, *stl_filename, id2id, is_connection,
- desc_id, &stardis->geometry);
+ res = read_triangles(&stl_desc, triangle2id, stl_filename, id2id, file_type,
+ desc_id, &stardis->geometry, &stardis->boundary);
if (res != RES_OK) goto error;
end:
@@ -675,6 +693,7 @@ static res_T
geometry_analyse
(const char* medium_filename,
const char* bc_filename,
+ const char* bound_filename, /* Can be NULL */
struct stardis* stardis)
{
res_T res = RES_OK;
@@ -726,7 +745,8 @@ geometry_analyse
desc_id = *p_desc;
}
- res = read_stl(sz-1, 0, &stl_filename, &vertex2id, &triangle2id, stardis);
+ res = read_stl(sz-1, CONTENT_GEOMETRY, stl_filename, &vertex2id,
+ &triangle2id, stardis);
if (stl_filename) free(stl_filename);
if (res != RES_OK) goto error;
}
@@ -764,13 +784,21 @@ geometry_analyse
desc_id = *p_desc;
}
- res = read_stl(desc_id, 1, &stl_filename, &vertex2id, &triangle2id, stardis);
+ res = read_stl(desc_id, CONTENT_CONNECTION, stl_filename, &vertex2id,
+ &triangle2id, stardis);
if (stl_filename) free(stl_filename);
if (res != RES_OK) goto error;
}
fclose(input);
input = NULL;
+ /* Process vtk file describing boundary compute */
+ if (bound_filename) {
+ res = read_stl(0, CONTENT_BOUNDARY, bound_filename, &vertex2id,
+ &triangle2id, stardis);
+ if (res != RES_OK) goto error;
+ }
+
exit:
if(input) fclose(input);
if (tables_initialised) {
@@ -786,7 +814,7 @@ error:
static res_T
parse_camera
-(char* cam_param, struct camera* cam)
+ (char* cam_param, struct camera* cam)
{
char** line = NULL;
char** opt = NULL;
@@ -801,7 +829,10 @@ parse_camera
{
size_t len = 0;
opt = split_line(line[i],'=');
- if (strcmp(opt[0],"fov")==0){
+ 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) {
@@ -848,7 +879,7 @@ stardis_init
if (res != RES_OK) goto error;
res = geometry_analyse(args->medium_filename,
- args->bc_filename, stardis);
+ args->bc_filename, args->solve_boundary_filename, stardis);
if (res != RES_OK) goto error;
stardis->N = args->N;
@@ -871,7 +902,10 @@ stardis_init
if (res != RES_OK) goto error;
}
else if (args->mode == MEDIUM_COMPUTE) {
- strcpy(stardis->medium_name, args->medium_name);
+ strcpy(stardis->solve_name, args->medium_name);
+ }
+ else if (args->mode == BOUNDARY_COMPUTE) {
+ strcpy(stardis->solve_name, args->solve_boundary_filename);
}
exit:
@@ -911,6 +945,8 @@ stardis_release(struct stardis* stardis)
}
sa_release(stardis->descriptions);
release_geometry(&stardis->geometry);
+ sa_release(stardis->boundary.primitives);
+ sa_release(stardis->boundary.sides);
}
/* TODO: rewrite dump! */
diff --git a/src/stardis-app.h b/src/stardis-app.h
@@ -900,17 +900,32 @@ struct camera {
double fov;
unsigned spp;
unsigned img[2];
+ union { /* Trick to allow static initialization with INF */
+ uint64_t t;
+ double time;
+ } u;
};
-#define NULL_CAMERA__ {{1,1,1},{0,0,0},{0,0,1},30,4,{640,480}}
+#define NULL_CAMERA__ {\
+ {1,1,1},{0,0,0},{0,0,1},30,4,{640,480},\
+ { 0x7FF0000000000000 } /* time=INF */\
+}
static const struct camera NULL_CAMERA = NULL_CAMERA__;
+struct solve_boundary {
+ size_t* primitives;
+ enum sdis_side* sides;
+};
+
struct stardis {
struct geometry geometry;
struct description* descriptions; /*array of materials and boundaries */
double probe[4]; /* x,y,z,t of probe when mode is PROBE_COMPUTE */
struct camera camera; /* camera when mode is IR_COMPUTE */
- char medium_name[64]; /* medium name when mode is MEDIUM_COMPUTE */
+ char solve_name[64]; /* medium name when mode is MEDIUM_COMPUTE,
+ boundary file name when BOUNDARY_COMPUTE*/
+ struct solve_boundary boundary; /* Boundary to solve when mode
+ is BOUNDARY_COMPUTE */
size_t N; /*number of MC realizations*/
unsigned nthreads;
@@ -921,11 +936,20 @@ struct stardis {
enum sdis_heat_path_flag dump_paths;
};
#define NULL_ALLOCATOR__ {NULL}
-#define NULL_STARDIS__ {NULL_GEOMETRY__, NULL,\
- {0,0,0,0}, NULL_CAMERA__, "",\
+#define NULL_STARDIS__ {\
+ NULL_GEOMETRY__, NULL,\
+ {0,0,0,0}, NULL_CAMERA__, "", {NULL,NULL},\
0, 0, 0, {300,300}, NULL_ALLOCATOR__, 0, SDIS_HEAT_PATH_NONE}
static const struct stardis NULL_STARDIS = NULL_STARDIS__;
+enum content_type {
+ CONTENT_GEOMETRY,
+ CONTENT_CONNECTION,
+ CONTENT_BOUNDARY,
+ CONTENT_TYPES_COUNT
+};
+
+
extern res_T
stardis_init
(const struct args* args,
@@ -934,10 +958,19 @@ stardis_init
extern res_T
stardis_compute(struct stardis* stardis, enum stardis_mode mode);
-void
+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
@@ -1145,7 +1145,7 @@ stardis_compute(struct stardis* stardis, enum stardis_mode mode)
/* Launch the simulation */
res = sdis_solve_camera(scn,
cam,
- INF,
+ stardis->camera.u.time,
stardis->scale_factor,
stardis->radiative_temp[0],
stardis->radiative_temp[1],
@@ -1208,14 +1208,14 @@ stardis_compute(struct stardis* stardis, enum stardis_mode mode)
FOR_EACH(ii, 0, sz) {
struct description* desc = stardis->descriptions + ii;
if (desc->type == DESC_MAT_SOLID) {
- if (0 == strcmp(stardis->medium_name, desc->d.solid.name)) {
+ 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->medium_name, desc->d.fluid.name)) {
+ if (0 == strcmp(stardis->solve_name, desc->d.fluid.name)) {
ASSERT(sa_size(media) > ii);
medium = media[ii];
break;
@@ -1225,7 +1225,7 @@ stardis_compute(struct stardis* stardis, enum stardis_mode mode)
if (medium == NULL) {
/* Not found */
fprintf(stderr, "Cannot solve medium %s (unknown medium)\n",
- stardis->medium_name);
+ stardis->solve_name);
res = RES_BAD_ARG;
goto error;
}
@@ -1239,24 +1239,60 @@ stardis_compute(struct stardis* stardis, enum stardis_mode mode)
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 == MEDIUM_COMPUTE) {
+ if (mode == PROBE_COMPUTE
+ || 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 */
- if (mode == PROBE_COMPUTE)
+ switch (mode) {
+ case PROBE_COMPUTE:
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 */
- else
+ break;
+ case MEDIUM_COMPUTE:
printf("Temperature in medium %s at t=%g = %g +/- %g\n",
- stardis->medium_name, time[0],
+ 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);
