stardis-solver

Solve coupled heat transfers
git clone git://git.meso-star.fr/stardis-solver.git
Log | Files | Refs | README | LICENSE

commit ed81eccde198e58150af3471986200424a86ca42
parent 30310f0f9aaac5bdf556138d2ee8e32c57639153
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Thu,  2 May 2019 16:11:11 +0200

Implement a more efficient 'scene_get_medium' procedure

Diffstat:
Msrc/sdis_heat_path_boundary_Xd.h | 11+++++++----
Msrc/sdis_heat_path_conductive_Xd.h | 9+++++----
Msrc/sdis_heat_path_convective_Xd.h | 5++++-
Msrc/sdis_heat_path_radiative_Xd.h | 7+++++--
Msrc/sdis_scene.c | 11+++++++++++
Msrc/sdis_scene_Xd.h | 103+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++--------------
Msrc/sdis_scene_c.h | 17+++++++++++++++++
7 files changed, 134 insertions(+), 29 deletions(-)

diff --git a/src/sdis_heat_path_boundary_Xd.h b/src/sdis_heat_path_boundary_Xd.h @@ -90,6 +90,7 @@ XD(select_reinjection_dir) struct sdis_interface* interf; struct sdis_medium* mdm0; struct sdis_medium* mdm1; + struct hit_filter_data filter_data; struct sXd(hit) hit; struct sXd(hit) hit0; struct sXd(hit) hit1; @@ -109,13 +110,15 @@ XD(select_reinjection_dir) f2(range, 0, FLT_MAX); fX_set_dX(org, rwalk->vtx.P); - SXD(scene_view_trace_ray(scn->sXd(view), org, dir0, range, &rwalk->hit, &hit0)); - SXD(scene_view_trace_ray(scn->sXd(view), org, dir1, range, &rwalk->hit, &hit1)); + filter_data.XD(hit) = rwalk->hit; + filter_data.epsilon = delta * 0.01; + SXD(scene_view_trace_ray(scn->sXd(view), org, dir0, range, &filter_data, &hit0)); + SXD(scene_view_trace_ray(scn->sXd(view), org, dir1, range, &filter_data, &hit1)); /* Retrieve the medium at the reinjection pos along dir0 */ if(SXD_HIT_NONE(&hit0)) { XD(move_pos)(dX(set)(tmp, rwalk->vtx.P), dir0, (float)delta); - res = scene_get_medium(scn, tmp, NULL, &mdm0); + res = scene_get_medium_in_closed_boundaries(scn, tmp, &mdm0); if(res != RES_OK) goto error; } else { interf = scene_get_interface(scn, hit0.prim.prim_id); @@ -126,7 +129,7 @@ XD(select_reinjection_dir) /* Retrieve the medium at the reinjection pos along dir1 */ if(SXD_HIT_NONE(&hit1)) { XD(move_pos)(dX(set)(tmp, rwalk->vtx.P), dir1, (float)delta); - res = scene_get_medium(scn, tmp, NULL, &mdm1); + res = scene_get_medium_in_closed_boundaries(scn, tmp, &mdm1); if(res != RES_OK) goto error; } else { interf = scene_get_interface(scn, hit1.prim.prim_id); diff --git a/src/sdis_heat_path_conductive_Xd.h b/src/sdis_heat_path_conductive_Xd.h @@ -134,7 +134,6 @@ XD(sample_next_step_robust) fX_set_dX(org, pos); do { - struct get_medium_info info; double pos_next[DIM]; /* Compute the next step */ @@ -144,7 +143,7 @@ XD(sample_next_step_robust) /* Retrieve the medium of the next step */ if(hit0->distance > delta) { XD(move_pos)(dX(set)(pos_next, pos), dir0, delta); - res = scene_get_medium(scn, pos_next, &info, &mdm); + res = scene_get_medium_in_closed_boundaries(scn, pos_next, &mdm); if(res != RES_OK) goto error; } else { struct sdis_interface* interf; @@ -156,6 +155,7 @@ XD(sample_next_step_robust) /* Check medium consistency */ if(current_mdm != mdm) { +#if 0 #if DIM == 2 log_err(scn->dev, "%s: inconsistent medium during the solid random walk at {%g, %g}.\n", @@ -165,6 +165,7 @@ XD(sample_next_step_robust) "%s: inconsistent medium during the solid random walk at {%g, %g, %g}.\n", FUNC_NAME, SPLIT3(pos)); #endif +#endif } } while(current_mdm != mdm && ++iattempt < MAX_ATTEMPTS); @@ -214,8 +215,8 @@ XD(conductive_path) (void)ctx, (void)istep; /* Check the random walk consistency */ - CHK(scene_get_medium(scn, rwalk->vtx.P, NULL, &mdm) == RES_OK); - if(mdm != rwalk->mdm) { + res = scene_get_medium_in_closed_boundaries(scn, rwalk->vtx.P, &mdm); + if(res != RES_OK || mdm != rwalk->mdm) { log_err(scn->dev, "%s: invalid solid random walk. " "Unexpected medium at {%g, %g, %g}.\n", FUNC_NAME, SPLIT3(rwalk->vtx.P)); res = RES_BAD_OP_IRRECOVERABLE; diff --git a/src/sdis_heat_path_convective_Xd.h b/src/sdis_heat_path_convective_Xd.h @@ -34,6 +34,7 @@ XD(register_heat_vertex_in_fluid) const double weight) { struct sdis_rwalk_vertex vtx = SDIS_RWALK_VERTEX_NULL; + struct hit_filter_data filter_data; const float empirical_dst = 0.1f; const float range[2] = {0, FLT_MAX}; float org[DIM]; @@ -50,7 +51,9 @@ XD(register_heat_vertex_in_fluid) fX(set)(dir, rwalk->hit.normal); if(rwalk->hit_side == SDIS_BACK) fX(minus)(dir, dir); - SXD(scene_view_trace_ray(scn->sXd(view), org, dir, range, &rwalk->hit, &hit)); + filter_data.XD(hit) = rwalk->hit; + filter_data.epsilon = 1.e-6; + SXD(scene_view_trace_ray(scn->sXd(view), org, dir, range, &filter_data, &hit)); dst = SXD_HIT_NONE(&hit) ? empirical_dst : hit.distance * 0.5f; vtx = rwalk->vtx; diff --git a/src/sdis_heat_path_radiative_Xd.h b/src/sdis_heat_path_radiative_Xd.h @@ -52,6 +52,7 @@ XD(trace_radiative_path) /* Launch the radiative random walk */ for(;;) { const struct sdis_interface* interf = NULL; + struct hit_filter_data filter_data; struct sdis_interface_fragment frag = SDIS_INTERFACE_FRAGMENT_NULL; struct sdis_medium* chk_mdm = NULL; double alpha; @@ -63,12 +64,14 @@ XD(trace_radiative_path) fX_set_dX(pos, rwalk->vtx.P); /* Trace the radiative ray */ + filter_data.XD(hit) = rwalk->hit; + filter_data.epsilon = 1.e-6; #if (SDIS_XD_DIMENSION == 2) SXD(scene_view_trace_ray_3d - (scn->sXd(view), pos, dir, range, &rwalk->hit, &rwalk->hit)); + (scn->sXd(view), pos, dir, range, &filter_data, &rwalk->hit)); #else SXD(scene_view_trace_ray - (scn->sXd(view), pos, dir, range, &rwalk->hit, &rwalk->hit)); + (scn->sXd(view), pos, dir, range, &filter_data, &rwalk->hit)); #endif if(SXD_HIT_NONE(&rwalk->hit)) { /* Fetch the ambient radiative temperature */ rwalk->hit_side = SDIS_SIDE_NULL__; diff --git a/src/sdis_scene.c b/src/sdis_scene.c @@ -373,3 +373,14 @@ scene_get_medium : scene_get_medium_3d(scn, pos, info, out_medium); } +res_T +scene_get_medium_in_closed_boundaries + (const struct sdis_scene* scn, + const double pos[], + struct sdis_medium** out_medium) +{ + return scene_is_2d(scn) + ? scene_get_medium_in_closed_boundaries_2d(scn, pos, out_medium) + : scene_get_medium_in_closed_boundaries_3d(scn, pos, out_medium); +} + diff --git a/src/sdis_scene_Xd.h b/src/sdis_scene_Xd.h @@ -166,9 +166,8 @@ clear_properties(struct sdis_scene* scn) * approximative way is to test that its position have at least one barycentric * coordinate roughly equal to 0 or 1. */ static FINLINE int -hit_on_vertex(const struct s2d_hit* hit) +hit_on_vertex(const struct s2d_hit* hit, const float on_vertex_eps) { - const float on_vertex_eps = 1.e-4f; float v; ASSERT(hit && !S2D_HIT_NONE(hit)); v = 1.f - hit->u; @@ -183,9 +182,8 @@ hit_on_vertex(const struct s2d_hit* hit) * simple but approximative way is to test that its position have at least one * barycentric coordinate roughly equal to 0 or 1. */ static FINLINE int -hit_on_edge(const struct s3d_hit* hit) +hit_on_edge(const struct s3d_hit* hit, const float on_edge_eps) { - const float on_edge_eps = 1.e-4f; float w; ASSERT(hit && !S3D_HIT_NONE(hit)); w = 1.f - hit->uv[0] - hit->uv[1]; @@ -206,23 +204,21 @@ XD(hit_filter_function) const float org[DIM], const float dir[DIM], void* ray_data, - void* filter_data) + void* global_data) { - const struct sXd(hit)* hit_from = ray_data; - (void)org, (void)dir, (void)filter_data; + const struct hit_filter_data* filter_data = ray_data; + const struct sXd(hit)* hit_from = &filter_data->XD(hit); + (void)org, (void)dir, (void)global_data; - if(!hit_from || SXD_HIT_NONE(hit_from)) return 0; /* No filtering */ + if(!ray_data || SXD_HIT_NONE(hit_from)) return 0; /* No filtering */ if(SXD_PRIMITIVE_EQ(&hit_from->prim, &hit->prim)) return 1; - if(eq_epsf(hit->distance, 0, 1.e-6f)) { + if(eq_epsf(hit->distance, 0, (float)filter_data->epsilon)) { /* If the targeted point is near of the origin, check that it lies on an * edge/vertex shared by the 2 primitives. */ -#if DIM == 2 - return hit_on_vertex(hit_from) && hit_on_vertex(hit); -#else - return hit_on_edge(hit_from) && hit_on_edge(hit); -#endif + return HIT_ON_BOUNDARY(hit_from, 1.e-4f) + && HIT_ON_BOUNDARY(hit, 1e-4f); } return 0; } @@ -801,6 +797,7 @@ XD(scene_get_medium) size_t iprim, nprims; size_t nfailures = 0; const size_t max_failures = 10; + float P[DIM]; /* Range of the parametric coordinate into which positions are challenged */ #if DIM == 2 float st[3]; @@ -821,15 +818,24 @@ XD(scene_get_medium) f2(st[2], 5.f/12.f, 5.f/12.f); #endif + fX_set_dX(P, pos); + SXD(scene_view_primitives_count(scn->sXd(view), &nprims)); FOR_EACH(iprim, 0, nprims) { struct sXd(hit) hit; struct sXd(attrib) attr; struct sXd(primitive) prim; const float range[2] = {0.f, FLT_MAX}; - float N[DIM], P[DIM], dir[DIM], cos_N_dir; + float N[DIM], dir[DIM], cos_N_dir; size_t istep = 0; + if(iprim && (iprim % 100) == 0) { + log_warn(scn->dev, + "%s: performance issue. Up to %lu primitives were tested to define the " + "current medium at {%g, %g, %g}.\n", + FUNC_NAME, (unsigned long)iprim, SPLIT3(P)); + } + do { /* Retrieve a position onto the primitive */ SXD(scene_view_get_primitive(scn->sXd(view), (unsigned)iprim, &prim)); @@ -837,7 +843,7 @@ XD(scene_get_medium) /* Trace a ray from the random walk vertex toward the retrieved primitive * position */ - fX(normalize)(dir, fX(sub)(dir, attr.value, fX_set_dX(P, pos))); + fX(normalize)(dir, fX(sub)(dir, attr.value, P)); SXD(scene_view_trace_ray(scn->sXd(view), P, dir, range, NULL, &hit)); /* Unforeseen error. One has to intersect a primitive ! */ @@ -852,7 +858,7 @@ XD(scene_get_medium) } /* Discard the hit if it is on a vertex/edge, and target a new position * onto the current primitive */ - } while((SXD_HIT_NONE(&hit) || HIT_ON_BOUNDARY(&hit)) + } while((SXD_HIT_NONE(&hit) || HIT_ON_BOUNDARY(&hit, 1.e-4f)) && ++istep < nsteps); /* The hits of all targeted positions on the current primitive are on @@ -863,7 +869,7 @@ XD(scene_get_medium) cos_N_dir = fX(dot)(N, dir); /* Not too close and not roughly orthognonal */ - if(hit.distance > 1.e-6 && absf(cos_N_dir) > 1.e-1f) { + if(hit.distance > 1.e-6 && absf(cos_N_dir) > 1.e-2f) { const struct sdis_interface* interf; interf = scene_get_interface(scn, hit.prim.prim_id); medium = interface_get_medium @@ -893,6 +899,67 @@ error: #endif goto exit; } + +static INLINE res_T +XD(scene_get_medium_in_closed_boundaries) + (const struct sdis_scene* scn, + const double pos[DIM], + struct sdis_medium** out_medium) +{ + struct sdis_medium* medium = NULL; + const float dirs[6][3] = {{1,0,0},{-1,0,0},{0,1,0},{0,-1,0},{0,0,1},{0,0,-1}}; + float P[DIM]; + int idir; + res_T res = RES_OK; + ASSERT(scn && pos); + + fX_set_dX(P, pos); + FOR_EACH(idir, 0, 2*DIM) { + struct sXd(hit) hit; + const float range[2] = {0.f, FLT_MAX}; + float N[DIM], cos_N_dir; + + /* Trace a ray from the random walk vertex toward the retrieved primitive + * position */ + SXD(scene_view_trace_ray(scn->sXd(view), P, dirs[idir], range, NULL, &hit)); + + /* Unforeseen error. One has to intersect a primitive ! */ + if(SXD_HIT_NONE(&hit)) continue; + + /* Discard a hits if it lies on an edge/point */ + if(HIT_ON_BOUNDARY(&hit, 1.e-4f)) continue; + + fX(normalize)(N, hit.normal); + cos_N_dir = fX(dot)(N, dirs[idir]); + + /* Not too close and not roughly orthognonal */ + if(hit.distance > 1.e-6 && absf(cos_N_dir) > 1.e-2f) { + const struct sdis_interface* interf; + interf = scene_get_interface(scn, hit.prim.prim_id); + medium = interface_get_medium + (interf, cos_N_dir < 0 ? SDIS_FRONT : SDIS_BACK); + break; + } + } + if(idir >= 2*DIM) { + res = XD(scene_get_medium)(scn, pos, NULL, &medium); + if(res != RES_OK) goto error; + } + +exit: + *out_medium = medium; + return res; +error: +#if DIM == 2 + log_err(scn->dev, "%s: could not retrieve the medium at {%g, %g}.\n", + FUNC_NAME, SPLIT2(pos)); +#else + log_err(scn->dev, "%s: could not retrieve the medium at {%g, %g, %g}.\n", + FUNC_NAME, SPLIT3(pos)); +#endif + goto exit; +} + #undef SDIS_SCENE_DIMENSION #undef DIM #undef sencXd diff --git a/src/sdis_scene_c.h b/src/sdis_scene_c.h @@ -31,6 +31,12 @@ struct prim_prop { unsigned back_enclosure; /* Id of the back facing enclosure */ }; +struct hit_filter_data { + struct s2d_hit hit_2d; + struct s3d_hit hit_3d; + double epsilon; /* Threshold defining roughly equal intersections */ +}; + struct get_medium_info { /* Targeted position */ float pos_tgt[3]; @@ -223,6 +229,17 @@ scene_get_medium struct get_medium_info* info, /* May be NULL */ struct sdis_medium** medium); +/* This function assumes that the tested position lies into finite enclosure. + * The medium into which it lies is thus retrieved by tracing rays along the + * main axis. For possible infinite enclosure, one has to use the + * `scene_get_medium' function instead that, in counterpart, can be more time + * consuming. */ +extern LOCAL_SYM res_T +scene_get_medium_in_closed_boundaries + (const struct sdis_scene* scn, + const double position[], + struct sdis_medium** medium); + static INLINE void scene_get_enclosure_ids (const struct sdis_scene* scn,