htrdr

Solving radiative transfer in heterogeneous media
git clone git://git.meso-star.fr/htrdr.git
Log | Files | Refs | README | LICENSE

htrdr_combustion_laser.c (8840B)


      1 /* Copyright (C) 2018-2019, 2022-2025 Centre National de la Recherche Scientifique
      2  * Copyright (C) 2020-2022 Institut Mines Télécom Albi-Carmaux
      3  * Copyright (C) 2022-2025 Institut Pierre-Simon Laplace
      4  * Copyright (C) 2022-2025 Institut de Physique du Globe de Paris
      5  * Copyright (C) 2018-2025 |Méso|Star> (contact@meso-star.com)
      6  * Copyright (C) 2022-2025 Observatoire de Paris
      7  * Copyright (C) 2022-2025 Université de Reims Champagne-Ardenne
      8  * Copyright (C) 2022-2025 Université de Versaille Saint-Quentin
      9  * Copyright (C) 2018-2019, 2022-2025 Université Paul Sabatier
     10  *
     11  * This program is free software: you can redistribute it and/or modify
     12  * it under the terms of the GNU General Public License as published by
     13  * the Free Software Foundation, either version 3 of the License, or
     14  * (at your option) any later version.
     15  *
     16  * This program is distributed in the hope that it will be useful,
     17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
     18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
     19  * GNU General Public License for more details.
     20  *
     21  * You should have received a copy of the GNU General Public License
     22  * along with this program. If not, see <http://www.gnu.org/licenses/>. */
     23 
     24 #include "combustion/htrdr_combustion_laser.h"
     25 
     26 #include "core/htrdr.h"
     27 #include "core/htrdr_log.h"
     28 #include "core/htrdr_rectangle.h"
     29 
     30 #include <rsys/cstr.h>
     31 #include <rsys/double33.h>
     32 #include <rsys/mem_allocator.h>
     33 #include <rsys/ref_count.h>
     34 
     35 struct htrdr_combustion_laser {
     36   struct htrdr_rectangle* surface; /* Surface emission */
     37   double world2local[12];
     38   double low_ls[3], upp_ls[3]; /* Local space AABB */
     39   double flux_density; /* In W/m^2 */
     40   double wavelength; /* In nm */
     41   ref_T ref;
     42   struct htrdr* htrdr;
     43 };
     44 
     45 /*******************************************************************************
     46  * Helper functions
     47  ******************************************************************************/
     48 static void
     49 laser_release(ref_T* ref)
     50 {
     51   struct htrdr_combustion_laser* laser;
     52   struct htrdr* htrdr;
     53   ASSERT(ref);
     54   laser = CONTAINER_OF(ref, struct htrdr_combustion_laser, ref);
     55   if(laser->surface) htrdr_rectangle_ref_put(laser->surface);
     56   htrdr = laser->htrdr;
     57   MEM_RM(htrdr_get_allocator(htrdr), laser);
     58   htrdr_ref_put(htrdr);
     59 }
     60 
     61 /*******************************************************************************
     62  * Local functions
     63  ******************************************************************************/
     64 res_T
     65 htrdr_combustion_laser_create
     66   (struct htrdr* htrdr,
     67    const struct htrdr_combustion_laser_create_args* args,
     68    struct htrdr_combustion_laser** out_laser)
     69 {
     70   struct htrdr_combustion_laser* laser = NULL;
     71   res_T res = RES_OK;
     72   ASSERT(htrdr && args && out_laser);
     73 
     74   if(args->flux_density <= 0) {
     75     htrdr_log_err(htrdr, "Invalid flux density %g W.m^2\n", args->flux_density);
     76     res = RES_BAD_ARG;
     77     goto error;
     78   }
     79 
     80   if(args->wavelength <= 0) {
     81     htrdr_log_err(htrdr, "Invalid wavelength %g nm\n", args->wavelength);
     82     res = RES_BAD_ARG;
     83     goto error;
     84   }
     85 
     86   laser = MEM_CALLOC(htrdr_get_allocator(htrdr), 1, sizeof(*laser));
     87   if(!laser) {
     88     htrdr_log_err(htrdr, "Could not allocate the laser data structure.\n");
     89     res = RES_MEM_ERR;
     90     goto error;
     91   }
     92   ref_init(&laser->ref);
     93   htrdr_ref_get(htrdr);
     94   laser->htrdr = htrdr;
     95   laser->flux_density = args->flux_density;
     96   laser->wavelength = args->wavelength;
     97 
     98   res = htrdr_rectangle_create
     99     (laser->htrdr,
    100      args->surface.size,
    101      args->surface.position,
    102      args->surface.target,
    103      args->surface.up,
    104      &laser->surface);
    105   if(res != RES_OK) {
    106     htrdr_log_err(htrdr, "Could not create the laser surface emission -- %s.\n",
    107       res_to_cstr(res));
    108     goto error;
    109   }
    110 
    111   htrdr_rectangle_get_transform_inverse(laser->surface, laser->world2local);
    112 
    113   /* Define the laser sheet AABB in local space */
    114   laser->upp_ls[0] = args->surface.size[0] * 0.5;
    115   laser->upp_ls[1] = args->surface.size[1] * 0.5;
    116   laser->upp_ls[2] = 0;
    117   laser->low_ls[0] = -laser->upp_ls[0];
    118   laser->low_ls[1] = -laser->upp_ls[1];
    119   laser->upp_ls[2] = INF;
    120 
    121 exit:
    122   *out_laser = laser;
    123   return res;
    124 error:
    125   if(laser) { htrdr_combustion_laser_ref_put(laser); laser = NULL; }
    126   goto exit;
    127 }
    128 
    129 void
    130 htrdr_combustion_laser_ref_get(struct htrdr_combustion_laser* laser)
    131 {
    132   ASSERT(laser);
    133   ref_get(&laser->ref);
    134 }
    135 
    136 void
    137 htrdr_combustion_laser_ref_put(struct htrdr_combustion_laser* laser)
    138 {
    139   ASSERT(laser);
    140   ref_put(&laser->ref, laser_release);
    141 }
    142 
    143 void
    144 htrdr_combustion_laser_trace_ray
    145   (struct htrdr_combustion_laser* laser,
    146    const double pos[3],
    147    const double dir[3],
    148    const double range[2],
    149    double distance[2])
    150 {
    151   double pos_ls[3]; /* Position in local space */
    152   double dir_ls[3]; /* Direction in local space */
    153   double t_min;
    154   double t_max;
    155   int i;
    156   ASSERT(laser && pos && dir && range && distance);
    157 
    158   /* Transform the ray in local space */
    159   d33_muld3(pos_ls, laser->world2local, pos);
    160   d33_muld3(dir_ls, laser->world2local, dir);
    161   d3_add(pos_ls, pos_ls, laser->world2local+9);
    162 
    163   /* Reset the distance */
    164   distance[0] = INF;
    165   distance[1] = INF;
    166 
    167   /* Initialise the range */
    168   t_min = range[0];
    169   t_max = range[1];
    170 
    171   /* TODO one can optimise the Ray/AABB intersection test along the Z axis
    172    * whose AABB interval is in [0, inf) */
    173   FOR_EACH(i, 0, 3) {
    174     const double rcp_dir_ls = 1.0/dir_ls[i];
    175     double t0 = (laser->low_ls[i] - pos_ls[i]) * rcp_dir_ls;
    176     double t1 = (laser->upp_ls[i] - pos_ls[i]) * rcp_dir_ls;
    177     if(t0 > t1) SWAP(double, t0, t1);
    178     t_min = MMAX(t_min, t0);
    179     t_max = MMIN(t_max, t1);
    180     if(t_min > t_max) return; /* No intersection */
    181   }
    182 
    183   /* Save the hit distance */
    184   distance[0] = t_min;
    185   distance[1] = t_max;
    186 }
    187 
    188 void
    189 htrdr_combustion_laser_get_mesh
    190   (const struct htrdr_combustion_laser* laser,
    191    const double extent,
    192    struct htrdr_combustion_laser_mesh* mesh)
    193 {
    194   double transform[12];
    195   double low[3];
    196   double upp[3];
    197   int i;
    198   ASSERT(laser && extent > 0 && mesh);
    199 
    200   htrdr_rectangle_get_transform(laser->surface, transform);
    201 
    202   /* Define the laser sheet vertices in local space */
    203   low[0] = laser->low_ls[0];
    204   low[1] = laser->low_ls[1];
    205   low[2] = laser->low_ls[2];
    206   upp[0] = laser->upp_ls[0];
    207   upp[1] = laser->upp_ls[1];
    208   upp[2] = laser->low_ls[2] + extent;
    209 
    210   /* Define the mesh of the laser sheet in local space */
    211   d3(mesh->vertices + 0*3, low[0], low[1], low[2]);
    212   d3(mesh->vertices + 1*3, upp[0], low[1], low[2]);
    213   d3(mesh->vertices + 2*3, low[0], upp[1], low[2]);
    214   d3(mesh->vertices + 3*3, upp[0], upp[1], low[2]);
    215   d3(mesh->vertices + 4*3, low[0], low[1], upp[2]);
    216   d3(mesh->vertices + 5*3, upp[0], low[1], upp[2]);
    217   d3(mesh->vertices + 6*3, low[0], upp[1], upp[2]);
    218   d3(mesh->vertices + 7*3, upp[0], upp[1], upp[2]);
    219   mesh->nvertices = 8;
    220 
    221   /* Transform the laser sheet vertices in world space */
    222   FOR_EACH(i, 0, 8) {
    223     double* v = mesh->vertices + i*3;
    224     d33_muld3(v, transform, v);
    225     d3_add(v, v, transform+9);
    226   }
    227 
    228   /* Define the laser sheet triangles */
    229   mesh->triangles[0]  = 0; mesh->triangles[1]  = 2; mesh->triangles[2]  = 1;
    230   mesh->triangles[3]  = 1; mesh->triangles[4]  = 2; mesh->triangles[5]  = 3;
    231   mesh->triangles[6]  = 0; mesh->triangles[7]  = 4; mesh->triangles[8]  = 2;
    232   mesh->triangles[9]  = 2; mesh->triangles[10] = 4; mesh->triangles[11] = 6;
    233   mesh->triangles[12] = 1; mesh->triangles[13] = 3; mesh->triangles[14] = 5;
    234   mesh->triangles[15] = 5; mesh->triangles[16] = 3; mesh->triangles[17] = 7;
    235   mesh->triangles[18] = 1; mesh->triangles[19] = 5; mesh->triangles[20] = 0;
    236   mesh->triangles[21] = 0; mesh->triangles[22] = 5; mesh->triangles[23] = 4;
    237   mesh->triangles[24] = 3; mesh->triangles[25] = 2; mesh->triangles[26] = 7;
    238   mesh->triangles[27] = 7; mesh->triangles[28] = 2; mesh->triangles[29] = 6;
    239   mesh->ntriangles = 10;
    240 }
    241 
    242 void
    243 htrdr_combustion_laser_get_position
    244   (const struct htrdr_combustion_laser* laser,
    245    double pos[3])
    246 {
    247   ASSERT(laser && pos);
    248   htrdr_rectangle_get_center(laser->surface, pos);
    249 }
    250 
    251 void
    252 htrdr_combustion_laser_get_direction
    253   (const struct htrdr_combustion_laser* laser,
    254    double dir[3])
    255 {
    256   ASSERT(laser && dir);
    257   htrdr_rectangle_get_normal(laser->surface, dir);
    258 }
    259 
    260 double
    261 htrdr_combustion_laser_get_flux_density
    262   (const struct htrdr_combustion_laser* laser)
    263 {
    264   ASSERT(laser);
    265   return laser->flux_density;
    266 }
    267 
    268 double
    269 htrdr_combustion_laser_get_wavelength
    270   (const struct htrdr_combustion_laser* laser)
    271 {
    272   ASSERT(laser);
    273   return laser->wavelength;
    274 }
    275 
    276 double
    277 htrdr_combustion_laser_compute_surface_plane_distance
    278   (const struct htrdr_combustion_laser* laser,
    279    const double pos[3])
    280 {
    281   double center[3];
    282   double normal[3];
    283   double vec[3];
    284   ASSERT(laser && pos);
    285   htrdr_rectangle_get_normal(laser->surface, normal);
    286   htrdr_rectangle_get_center(laser->surface, center);
    287   d3_sub(vec, pos, center);
    288   return d3_dot(normal, vec);
    289 }