stardis-test

sadist_lib_solid_sphere.c (11049B)

---

 /* Copyright (C) 2024 |Méso|Star> (contact@meso-star.com)
  *
  * This program is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
  * the Free Software Foundation, either version 3 of the License, or
  * (at your option) any later version.
  *
  * This program is distributed in the hope that it will be useful,
  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  * GNU General Public License for more details.
  *
  * You should have received a copy of the GNU General Public License
  * along with this program. If not, see <http://www.gnu.org/licenses/>. */
 
 #include <stardis/stardis-prog-properties.h>
 
 #include <star/s3d.h>
 #include <star/ssp.h>
 #include <star/sstl.h>
 
 #include <rsys/cstr.h>
 #include <rsys/double3.h>
 #include <rsys/mem_allocator.h>
 
 #include <getopt.h>
 
 struct args {
   const char* stl;
   double lambda; /* [W/m/K] */
   double rho; /* [kg/m^3] */
   double cp; /* [J/K/kg] */
   double radius; /* [m/fp_to_meter] */
 };
 #define ARGS_DEFAULT__ {NULL,1,1,1,1}
 static const struct args ARGS_DEFAULT = ARGS_DEFAULT__;
 
 /* Define a solid sphere */
 struct solid_sphere {
   struct s3d_scene_view* view;
   double lambda; /* [W/m/K] */
   double rho; /* [kg/m^3] */
   double cp; /* [J/K/kg] */
   double radius; /* [m/fp_to_meter] */
 };
 #define SOLID_SPHERE_NULL__ {NULL,0,0,0,0}
 static const struct solid_sphere SOLID_SPHERE_NULL = SOLID_SPHERE_NULL__;
 
 /*******************************************************************************
  * Helper functions
  ******************************************************************************/
 static void
 get_indices(const unsigned itri, unsigned ids[3], void* ctx)
 {
   const struct sstl_desc* desc = ctx;
   ASSERT(desc && itri < desc->triangles_count);
   ids[0] = desc->indices[itri*3+0];
   ids[1] = desc->indices[itri*3+1];
   ids[2] = desc->indices[itri*3+2];
 }
 
 static void
 get_position(const unsigned ivert, float pos[3], void* ctx)
 {
   const struct sstl_desc* desc = ctx;
   ASSERT(desc && ivert < desc->vertices_count);
   pos[0] = desc->vertices[ivert*3+0];
   pos[1] = desc->vertices[ivert*3+1];
   pos[2] = desc->vertices[ivert*3+2];
 }
 
 static struct s3d_scene_view*
 create_scene_view(const char* filename)
 {
   /* Star-STL */
   struct sstl_desc desc;
   struct sstl* sstl = NULL;
 
   /* Star3D */
   struct s3d_vertex_data vdata = S3D_VERTEX_DATA_NULL;
   struct s3d_device* dev = NULL;
   struct s3d_scene* scn = NULL;
   struct s3d_shape* shape = NULL;
   struct s3d_scene_view* view = NULL;
 
   res_T res = RES_OK;
   ASSERT(filename);
 
   if((res = sstl_create(NULL, NULL, 1, &sstl)) != RES_OK) goto error;
   if((res = sstl_load(sstl, filename)) != RES_OK) goto error;
   if((res = sstl_get_desc(sstl, &desc)) != RES_OK) goto error;
 
   if((res = s3d_device_create(NULL, NULL, 0, &dev)) != RES_OK) goto error;
   if((res = s3d_scene_create(dev, &scn)) != RES_OK) goto error;
   if((res = s3d_shape_create_mesh(dev, &shape)) != RES_OK) goto error;
   if((res = s3d_scene_attach_shape(scn, shape)) != RES_OK) goto error;
 
   vdata.usage = S3D_POSITION;
   vdata.type = S3D_FLOAT3;
   vdata.get = get_position;
   res = s3d_mesh_setup_indexed_vertices(shape,
     (unsigned int)desc.triangles_count, get_indices,
     (unsigned int)desc.vertices_count, &vdata, 1, &desc);
   if(res != RES_OK) goto error;
 
   res = s3d_scene_view_create(scn, S3D_TRACE|S3D_GET_PRIMITIVE, &view);
   if(res != RES_OK) goto error;
 
 exit:
   if(sstl) SSTL(ref_put(sstl));
   if(dev) S3D(device_ref_put(dev));
   if(scn) S3D(scene_ref_put(scn));
   if(shape) S3D(shape_ref_put(shape));
   return view;
 error:
   if(view) { S3D(scene_view_ref_put(view)); view = NULL; }
   goto exit;
 }
 
 static void
 print_usage(FILE* stream, const char* name)
 {
   ASSERT(name);
   fprintf(stream,
     "usage: %s -m mesh [-c capa] [-h] [-l lambda] [-R radius] [-r rho]\n",
     name);
 }
 
 static res_T
 parse_args
   (const struct stardis_description_create_context* ctx,
    struct args* args,
    const int argc,
    char* argv[])
 {
   int opt = 0;
   res_T res = RES_OK;
 
   /* Check pre-conditions */
   ASSERT(ctx && args);
   *args = ARGS_DEFAULT;
 
   optind = 1;
   while((opt = getopt(argc, argv, "c:hl:m:R:r:")) != -1) {
     switch(opt) {
       case 'c':
         res = cstr_to_double(optarg, &args->cp);
         if(res == RES_OK && args->cp <= 0) res = RES_BAD_ARG;
         break;
       case 'h':
         print_usage(stdout, argv[0]);
         break;
       case 'l':
         res = cstr_to_double(optarg, &args->lambda);
         if(res == RES_OK && args->lambda <= 0) res = RES_BAD_ARG;
         break;
       case 'm':
         args->stl = optarg;
         break;
       case 'R':
         res = cstr_to_double(optarg, &args->radius);
         if(res == RES_OK && args->radius <= 0) res = RES_BAD_ARG;
         break;
       case 'r':
         res = cstr_to_double(optarg, &args->rho);
         if(res == RES_OK && args->rho <= 0) res = RES_BAD_ARG;
         break;
       default: res = RES_BAD_ARG; break;
     }
     if(res != RES_OK) {
       if(optarg) {
         fprintf(stderr, "%s: invalid option argument '%s' -- '%c'\n",
           ctx->name, optarg, opt);
       }
       goto error;
     }
   }
 
   if(!args->stl) {
     fprintf(stderr, "%s: missing mesh -- option '-m'\n", argv[0]);
     res = RES_BAD_ARG;
     goto error;
   }
 
 exit:
   return res;
 error:
   goto exit;
 }
 
 /*******************************************************************************
  * Create data
  ******************************************************************************/
 void*
 stardis_create_data
   (const struct stardis_description_create_context *ctx,
    void* libdata,
    size_t argc,
    char* argv[])
 {
   struct args args = ARGS_DEFAULT__;
   struct solid_sphere* solid = NULL;
   res_T res = RES_OK;
   (void)libdata;
 
   solid = mem_alloc(sizeof(*solid));
   if(!solid) {
     fprintf(stderr, "%s:%d: error allocating the solid sphere.\n",
       __FILE__, __LINE__);
     goto error;
   }
   *solid = SOLID_SPHERE_NULL;
 
   res = parse_args(ctx, &args, (int)argc, argv);
   if(res != RES_OK) goto error;
 
   solid->lambda = args.lambda;
   solid->rho = args.rho;
   solid->cp = args.cp;
   solid->radius = args.radius;
 
   if(!(solid->view = create_scene_view(args.stl))) goto error;
 
 exit:
   return solid;
 error:
   if(solid) {
     stardis_release_data(solid);
     solid = NULL;
   }
   goto exit;
 }
 
 void
 stardis_release_data(void* data)
 {
   struct solid_sphere* solid = data;
   ASSERT(solid);
   if(solid->view) S3D(scene_view_ref_put(solid->view));
   mem_rm(solid);
 }
 
 /*******************************************************************************
  * Solid properties
  ******************************************************************************/
 double
 stardis_calorific_capacity(const struct stardis_vertex* vtx, void* data)
 {
   const struct solid_sphere* sphere = data;
   (void)vtx;
   return sphere->cp;
 }
 
 double
 stardis_volumic_mass(const struct stardis_vertex* vtx, void* data)
 {
   const struct solid_sphere* sphere = data;
   (void)vtx;
   return sphere->rho;
 }
 
 double
 stardis_conductivity(const struct stardis_vertex* vtx, void* data)
 {
   const struct solid_sphere* sphere = data;
   (void)vtx;
   return sphere->lambda;
 }
 
 double
 stardis_delta_solid(const struct stardis_vertex* vtx, void* data)
 {
   const struct solid_sphere* sphere = data;
   (void)vtx;
   return sphere->radius/20.0;
 }
 
 double
 stardis_volumic_power(const struct stardis_vertex* vtx, void* data)
 {
   (void)vtx, (void)data;
   return STARDIS_VOLUMIC_POWER_NONE;
 }
 
 double
 stardis_medium_temperature(const struct stardis_vertex* vtx, void* data)
 {
   (void)vtx, (void)data;
   return STARDIS_TEMPERATURE_NONE;
 }
 
 int
 stardis_sample_conductive_path
   (struct sdis_scene* scn,
    struct ssp_rng* rng,
    struct stardis_path* path,
    void* data)
 {
   const struct solid_sphere* sphere = data;
   struct s3d_hit hit = S3D_HIT_NULL;
   struct s3d_attrib attr0, attr1, attr2;
 
   double pos[3];
   double epsilon = 0;
   ASSERT(scn && rng && path && data);
   (void)scn;
 
   epsilon = sphere->radius * 1.e-6; /* Epsilon shell */
 
   d3_set(pos, path->vtx.P);
 
   do {
     /* Distance from the geometry center to the current position */
     const double dst = d3_len(pos);
 
     double dir[3] = {0,0,0};
     double r = 0; /* Radius */
 
     r = sphere->radius - dst;
     CHK(dst > 0);
 
     if(r > epsilon) {
       /* Uniformly sample a new position on the surrounding sphere */
       ssp_ran_sphere_uniform(rng, dir, NULL);
 
       /* Move to the new position */
       d3_muld(dir, dir, r);
       d3_add(pos, pos, dir);
 
     /* The current position is in the epsilon shell:
      * move it to the nearest interface position */
     } else {
       float posf[3];
 
       d3_set(dir, pos);
       d3_normalize(dir, dir);
       d3_muld(pos, dir, sphere->radius);
 
       /* Map the position to the sphere geometry */
       f3_set_d3(posf, pos);
       S3D(scene_view_closest_point(sphere->view, posf, (float)INF, NULL, &hit));
     }
 
   /* The calculation is performed in steady state, so the path necessarily stops
    * at a boundary */
   } while(S3D_HIT_NONE(&hit));
 
   /* Setup the path state */
   d3_set(path->vtx.P, pos);
   path->weight = 0;
   path->at_limit = 0;
 
   /* Configure the intersecting primitive. The vertices must be exactly the same
    * (i.e. bit-compatible) as those used internally by Stardis, as they are used
    * as the key to retrieve the corresponding Stardis triangle. Note that
    * Stardis also uses Star-STL internally to load its geometry. As a result,
    * the in-memory representation of the loaded data is the same between Stardis
    * and the current code, i.e. single-precision floating point. And Star-3D
    * also uses the same representation. We can therefore use Star-3D data
    * directly to define the coordinates of the vertices of the intersecting
    * primitive */
   S3D(triangle_get_vertex_attrib(&hit.prim, 0, S3D_POSITION, &attr0));
   S3D(triangle_get_vertex_attrib(&hit.prim, 1, S3D_POSITION, &attr1));
   S3D(triangle_get_vertex_attrib(&hit.prim, 2, S3D_POSITION, &attr2));
   d3_set_f3(path->tri.vtx0, attr0.value);
   d3_set_f3(path->tri.vtx1, attr1.value);
   d3_set_f3(path->tri.vtx2, attr2.value);
 
   return RES_OK;
 }
 
 /*******************************************************************************
  * Legal notices
  ******************************************************************************/
 const char*
 get_copyright_notice(void* data)
 {
   (void)data; /* Avoid "unused variable" warnings */
   return "Copyright (C) 2024 |Méso|Star> (contact@meso-star.com)";
 }
 
 const char*
 get_license_short(void* data)
 {
   (void)data; /* Avoid "unused variable" warnings */
   return "GNU GPL version 3 or later <http://www.gnu.org/licenses/>";
 }
 
 const char*
 get_license_text(void* data)
 {
   (void)data; /* Avoid "unused variable" warnings */
   return
     "This is free software released under the GPL v3+ license: GNU GPL\n"
     "version 3 or later. You are welcome to redistribute it under certain\n"
     "conditions; refer to <http://www.gnu.org/licenses/> for details.";
 }