commit 735156957b5f8fab315eafcca2d87e9816f1a274
parent 2858b60e4500d6060d9fdb03117c872176d80ab7
Author: Christophe Coustet <christophe.coustet@meso-star.com>
Date: Tue, 4 Jun 2019 12:13:08 +0200
Add probe on interface option
Probe point is no longer moved to interface if close, but is moved or
not according of the option used (-p VS -P).
Diffstat:
| M | src/args.h | | | 69 | ++++++++++++++++++++++++++++++++++++++++++++++++++++----------------- |
| M | src/stardis-compute.c | | | 116 | ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++--------------- |
2 files changed, 147 insertions(+), 38 deletions(-)
diff --git a/src/args.h b/src/args.h
@@ -18,7 +18,8 @@ enum stardis_mode {
GREEN_MODE = BIT(3),
IR_COMPUTE = BIT(4),
DUMP_VTK = BIT(5),
- COMPUTE_MODES = PROBE_COMPUTE | MEDIUM_COMPUTE | BOUNDARY_COMPUTE | IR_COMPUTE | DUMP_VTK
+ 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) {
@@ -29,6 +30,7 @@ static char mode_option(const enum stardis_mode 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;
}
@@ -93,8 +95,12 @@ print_help(char* prog)
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. If some boundary conditions are\n");
- printf(" time-dependant, TIME cannot be INF.\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");
@@ -152,8 +158,14 @@ parse_args(const int argc, char** argv, struct args* args)
goto error;
}
- while((opt = getopt(argc, argv, "hgn:t:b:B:M:m:p:dD:s:r:R:")) != -1) {
+ 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;
@@ -164,12 +176,13 @@ parse_args(const int argc, char** argv, struct args* args)
break;
case 'n': {
- args->N = strtoull(optarg, NULL, 10);
- if (args->N <= 0) {
+ long long n = strtoll(optarg, NULL, 10);
+ if (n <= 0 || n > SIZE_MAX) {
res = RES_BAD_ARG;
- fprintf(stderr, "Invalid argument -n %s\n", optarg);
+ fprintf(stderr, "Invalid argument -%c %s\n", opt, optarg);
goto error;
}
+ args->N = n;
break;
}
@@ -186,7 +199,7 @@ parse_args(const int argc, char** argv, struct args* args)
if (res != RES_OK
|| args->nthreads <= 0){
res = RES_BAD_ARG;
- fprintf(stderr, "Invalid argument -t %s\n", optarg);
+ fprintf(stderr, "Invalid argument -%c %s\n", opt, optarg);
goto error;
}
break;
@@ -200,10 +213,31 @@ parse_args(const int argc, char** argv, struct args* args)
}
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) {
+ || res != RES_OK
+ || args->u.probe[3] < 0)
+ {
res = RES_BAD_ARG;
- fprintf(stderr, "Invalid argument -p %s\n", optarg);
+ fprintf(stderr, "Invalid argument -%c %s\n", opt, optarg);
goto error;
}
break;
@@ -220,21 +254,22 @@ parse_args(const int argc, char** argv, struct args* args)
args->mode |= BOUNDARY_COMPUTE;
ptr = strrchr(optarg, ',');
/* Single ',' allowed */
- if (ptr != strchr(optarg, ','))
+ if (!ptr || ptr != strchr(optarg, ','))
err = 1;
else {
*ptr = '\0';
ptr++;
args->solve_boundary_filename = optarg;
- err = RES_OK != cstr_to_double(ptr, args->u.probe + 3);
+ 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 -b %s\n", optarg);
+ fprintf(stderr, "Invalid argument -%c %s\n", opt, optarg);
goto error;
}
break;
}
+
case 'm': {
char* ptr;
if (args->mode & COMPUTE_MODES) {
@@ -249,7 +284,7 @@ parse_args(const int argc, char** argv, struct args* args)
res = RES_BAD_ARG;
else res = cstr_to_double(ptr+1, args->u.probe+3);
if (res != RES_OK) {
- fprintf(stderr, "Invalid argument -m %s\n", optarg);
+ fprintf(stderr, "Invalid argument -%c %s\n", opt, optarg);
goto error;
}
*ptr = '\0';
@@ -276,7 +311,7 @@ parse_args(const int argc, char** argv, struct args* args)
args->dump_paths = DUMP_SUCCESS;
} else {
res = RES_BAD_ARG;
- fprintf(stderr, "Invalid argument -D %s\n", optarg);
+ fprintf(stderr, "Invalid argument -%c %s\n", opt, optarg);
goto error;
}
break;
@@ -286,7 +321,7 @@ parse_args(const int argc, char** argv, struct args* args)
if (res != RES_OK
|| args->scale_factor <=0){
res = RES_BAD_ARG;
- fprintf(stderr, "Invalid argument -s %s\n", optarg);
+ fprintf(stderr, "Invalid argument -%c %s\n", opt, optarg);
goto error;
}
break;
@@ -298,7 +333,7 @@ parse_args(const int argc, char** argv, struct args* args)
|| args->radiative_temp[0] <=0
|| args->radiative_temp[1] <=0){
res = RES_BAD_ARG;
- fprintf(stderr, "Invalid argument -r %s\n", optarg);
+ fprintf(stderr, "Invalid argument -%c %s\n", opt, optarg);
goto error;
}
break;
diff --git a/src/stardis-compute.c b/src/stardis-compute.c
@@ -4,6 +4,7 @@
#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>
@@ -332,41 +333,110 @@ static res_T
select_probe_type
(const struct sdis_scene* scn,
const struct triangle* triangle,
+ const struct vertex* vertex,
const struct description* descriptions,
+ const int on_interface,
double pos[3],
size_t* iprim,
double uv[2])
{
res_T res = RES_OK;
- size_t i = 0;
+ double min_d = DBL_MAX, min_delta = DBL_MAX, new_pos[3], tmp_uv[2], min_uv[2];
+ const struct description *min_desc = NULL;
+ size_t i = 0, found = SIZE_MAX;
+ enum description_type min_type = DESCRIPTION_TYPE_COUNT__;
for (i = 0; i < sa_size(triangle); ++i) {
- double delta = 0, dp[3], projected_pos[3];
- res = sdis_scene_boundary_project_position(scn, i, pos, uv);
+ const struct triangle* tr = triangle + i;
+ double dp[3], d, projected_pos[3], nn[3];
+ float v0[3], v1[3], v2[3], e1[3], e2[3], n[3];
+ const struct description *desc = NULL;
+ enum sdis_side s;
+
+ res = sdis_scene_boundary_project_position(scn, i, pos, tmp_uv);
if (res != RES_OK) return res;
- res = sdis_scene_get_boundary_position(scn, i, uv, projected_pos);
+ res = sdis_scene_get_boundary_position(scn, i, tmp_uv, projected_pos);
if (res != RES_OK) return res;
- d3_sub(dp, pos, projected_pos);
-
- if (triangle[i].front_description != UINT_MAX) {
- const struct description *desc = descriptions + triangle[i].front_description;
- if (desc->type == DESC_MAT_SOLID) delta = MMAX(desc->d.solid.delta, delta);
+ d3_sub(dp, projected_pos, pos);
+ d = d3_len(dp);
+
+ if (d == 0) {
+ /* Best possible match */
+ found = i;
+ fprintf(stderr,
+ "The probe is on the primitive %llu.\n", (long long int)*iprim);
+ break;
}
- if (triangle[i].back_description != UINT_MAX) {
- const struct description *desc = descriptions + triangle[i].back_description;
- if (desc->type == DESC_MAT_SOLID) delta = MMAX(desc->d.solid.delta, delta);
+ if (d >= min_d) {
+ /* No improvement */
+ continue;
+ }
+ min_d = d;
+ d2_set(min_uv, tmp_uv);
+
+ /* Find side */
+ f3_set(v0, vertex[tr->indices.data[0]].xyz);
+ f3_set(v1, vertex[tr->indices.data[1]].xyz);
+ f3_set(v2, vertex[tr->indices.data[2]].xyz);
+ f3_sub(e1, v1, v0);
+ f3_sub(e2, v2, v0);
+ f3_cross(n, e1, e2);
+
+ /* Stardis solver convention is triangle normal is back side */
+ d3_set_f3(nn, n);
+ s = (d3_dot(nn, dp) < 0) ? SDIS_BACK : SDIS_FRONT;
+
+ if (s == SDIS_FRONT) {
+ if (tr->front_description != UINT_MAX) {
+ desc = descriptions + tr->front_description;
+ }
+ } else {
+ if (tr->back_description != UINT_MAX) {
+ desc = descriptions + tr->back_description;
+ }
}
- if (d3_len(dp) < delta){
- *iprim = i;
- fprintf(stderr,"The probe is very close to the primitive %llu.\n",
- (long long int)*iprim );
- fprintf(stderr,"So the probe is moved to this primitive at the position: %g %g %g\n",
- SPLIT3(projected_pos));
- d3_set(pos, projected_pos);
- break;
+ min_desc = desc;
+ if (!desc) continue;
+
+ min_type = desc->type;
+ if (desc->type == DESC_MAT_SOLID) {
+ found = i;
+ min_delta = desc->d.solid.delta;
+ d3_set(new_pos, projected_pos);
}
}
+
+ if (!min_desc) {
+ /* Don't manage pos is outside the model
+ * but very close and could be considered on the boundary */
+ fprintf(stderr, "The probe is outside the model.\n");
+ return RES_BAD_ARG;
+ }
+
+ if (on_interface) {
+ *iprim = found;
+ if (min_type == DESC_MAT_FLUID)
+ fprintf(stderr, "Probe was on a fluid.\n");
+ fprintf(stderr, "The probe is moved to %g,%g,%g (primitive %llu)\n",
+ SPLIT3(new_pos), (long long int)*iprim);
+ if (min_type == DESC_MAT_SOLID && min_d > min_delta)
+ fprintf(stderr, "WARNING: The probe move is %.1f delta long.\n", min_d / min_delta);
+ d3_set(pos, new_pos);
+ d2_set(uv, min_uv);
+ } else {
+ if (min_type == DESC_MAT_FLUID)
+ fprintf(stderr,
+ "The probe is in a fluid: computing fluid temperature, not using a specific position.\n");
+ if (min_type == DESC_MAT_SOLID && min_d < min_delta) {
+ fprintf(stderr,
+ "WARNING: The probe is very close to the primitive %llu (%.1f delta).\n",
+ (long long int)found, min_d / min_delta);
+ fprintf(stderr,
+ "To have this probe moved on the closest interface use -P instead of -p\n");
+ }
+ }
+
return res;
}
@@ -1427,7 +1497,9 @@ stardis_compute(struct stardis* stardis, enum stardis_mode mode)
res = select_probe_type(scn,
stardis->geometry.triangle,
+ stardis->geometry.vertex,
stardis->descriptions,
+ (mode & PROBE_COMPUTE_ON_INTERFACE),
pos,
&iprim,
uv);
@@ -1646,7 +1718,9 @@ stardis_compute(struct stardis* stardis, enum stardis_mode mode)
res = select_probe_type(scn,
stardis->geometry.triangle,
+ stardis->geometry.vertex,
stardis->descriptions,
+ (mode & PROBE_COMPUTE_ON_INTERFACE),
pos,
&iprim,
uv);
@@ -1748,7 +1822,7 @@ stardis_compute(struct stardis* stardis, enum stardis_mode mode)
SDIS(estimator_get_failure_count(estimator, &nfailures));
/* Print the results */
- switch (mode) {
+ switch (mode & COMPUTE_MODES) {
case PROBE_COMPUTE:
printf("Temperature at t=[%g, %g, %g, %g] = %g +/- %g\n",
pos[0], pos[1], pos[2], time[0],
