rnatm_mesh.c (9587B)
1 /* Copyright (C) 2022, 2023, 2025 Centre National de la Recherche Scientifique 2 * Copyright (C) 2022, 2023, 2025 Institut Pierre-Simon Laplace 3 * Copyright (C) 2022, 2023, 2025 Institut de Physique du Globe de Paris 4 * Copyright (C) 2022, 2023, 2025 |Méso|Star> (contact@meso-star.com) 5 * Copyright (C) 2022, 2023, 2025 Observatoire de Paris 6 * Copyright (C) 2022, 2023, 2025 Université de Reims Champagne-Ardenne 7 * Copyright (C) 2022, 2023, 2025 Université de Versaille Saint-Quentin 8 * Copyright (C) 2022, 2023, 2025 Université Paul Sabatier 9 * 10 * This program is free software: you can redistribute it and/or modify 11 * it under the terms of the GNU General Public License as published by 12 * the Free Software Foundation, either version 3 of the License, or 13 * (at your option) any later version. 14 * 15 * This program is distributed in the hope that it will be useful, 16 * but WITHOUT ANY WARRANTY; without even the implied warranty of 17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 18 * GNU General Public License for more details. 19 * 20 * You should have received a copy of the GNU General Public License 21 * along with this program. If not, see <http://www.gnu.org/licenses/>. */ 22 23 #include "rnatm.h" 24 #include "rnatm_c.h" 25 #include "rnatm_log.h" 26 27 #include <star/smsh.h> 28 #include <star/suvm.h> 29 30 #include <rsys/clock_time.h> 31 32 /******************************************************************************* 33 * Helper functions 34 ******************************************************************************/ 35 static res_T 36 check_smsh_desc(struct rnatm* atm, const struct smsh_desc* desc) 37 { 38 res_T res = RES_OK; 39 ASSERT(atm && desc); 40 41 if(desc->dnode != 3 || desc->dcell != 4) { 42 log_err(atm, 43 "An atmosphere mesh must be a 3D tetrahedral mesh " 44 "(dimension of the mesh: %u; dimension of the vertices: %u)\n", 45 desc->dnode, desc->dcell); 46 res = RES_BAD_ARG; 47 goto error; 48 } 49 50 exit: 51 return res; 52 error: 53 goto exit; 54 } 55 56 static INLINE void 57 tetrahedron_get_indices(const size_t itetra, size_t ids[4], void* ctx) 58 { 59 const struct smsh_desc* desc = ctx; 60 const uint64_t* indices = NULL; 61 ASSERT(ctx && ids && itetra < desc->ncells && desc->dcell == 4); 62 indices = smsh_desc_get_cell(desc, itetra); 63 ids[0] = (size_t)indices[0]; 64 ids[1] = (size_t)indices[1]; 65 ids[2] = (size_t)indices[2]; 66 ids[3] = (size_t)indices[3]; 67 } 68 69 static INLINE void 70 vertex_get_position(const size_t ivert, double pos[3], void* ctx) 71 { 72 struct smsh_desc* desc = ctx; 73 const double* position = NULL; 74 ASSERT(ctx && pos && ivert < desc->nnodes && desc->dnode == 3); 75 position = smsh_desc_get_node(desc, ivert); 76 pos[0] = position[0]; 77 pos[1] = position[1]; 78 pos[2] = position[2]; 79 } 80 81 static res_T 82 setup_uvm 83 (struct rnatm* atm, 84 const struct rnatm_create_args* args, 85 const char* filename, 86 struct suvm_device* suvm, 87 struct smsh* smsh, 88 struct suvm_volume** out_volume, 89 size_t* ntetrahedra, 90 size_t* nvertices) 91 { 92 struct suvm_tetrahedral_mesh_args mesh_args = SUVM_TETRAHEDRAL_MESH_ARGS_NULL; 93 struct smsh_load_args smsh_load_args = SMSH_LOAD_ARGS_NULL; 94 struct smsh_desc smsh_desc = SMSH_DESC_NULL; 95 struct suvm_volume* volume = NULL; 96 res_T res = RES_OK; 97 ASSERT(atm && args && filename && suvm && smsh && out_volume); 98 99 /* Load and retrieve the Star-Mesh data */ 100 smsh_load_args.path = filename; 101 smsh_load_args.memory_mapping = 0; 102 res = smsh_load(smsh, &smsh_load_args); 103 if(res != RES_OK) goto error; 104 res = smsh_get_desc(smsh, &smsh_desc); 105 if(res != RES_OK) goto error; 106 res = check_smsh_desc(atm, &smsh_desc); 107 if(res != RES_OK) goto error; 108 109 /* Partition the unstructured volumetric mesh */ 110 mesh_args.ntetrahedra = smsh_desc.ncells; 111 mesh_args.nvertices = smsh_desc.nnodes; 112 mesh_args.get_indices = tetrahedron_get_indices; 113 mesh_args.get_position = vertex_get_position; 114 mesh_args.tetrahedron_data = SUVM_DATA_NULL; /* Tetra data are not in SUVM */ 115 mesh_args.vertex_data = SUVM_DATA_NULL; /* Vertex data are not in SUVM */ 116 mesh_args.precompute_normals = args->precompute_normals; 117 mesh_args.context = &smsh_desc; 118 res = suvm_tetrahedral_mesh_create(suvm, &mesh_args, &volume); 119 if(res != RES_OK) goto error; 120 121 exit: 122 *out_volume = volume; 123 *ntetrahedra = smsh_desc.ncells; 124 *nvertices = smsh_desc.nnodes; 125 return res; 126 error: 127 if(volume) { SUVM(volume_ref_put(volume)); volume = NULL; } 128 goto exit; 129 } 130 131 /******************************************************************************* 132 * Exported functions 133 ******************************************************************************/ 134 res_T 135 rnatm_fetch_cell 136 (const struct rnatm* atm, 137 const double pos[3], 138 const size_t cpnt, 139 struct rnatm_cell_pos* cell) 140 { 141 res_T res = RES_OK; 142 143 if(!atm || !pos || !cell) { 144 res = RES_BAD_ARG; 145 goto error; 146 } 147 148 cell->component = cpnt; 149 150 /* Gas */ 151 if(cpnt == RNATM_GAS) { 152 res = suvm_volume_at 153 (atm->gas.volume, pos, &cell->prim, cell->barycentric_coords); 154 if(res != RES_OK) { 155 log_err(atm, "Error retrieving gas cell at %g, %g, %g\n", 156 SPLIT3(pos)); 157 goto error; 158 } 159 160 /* Aerosol */ 161 } else if(cpnt < rnatm_get_aerosols_count(atm)) { 162 const struct aerosol* aerosol = darray_aerosol_cdata_get(&atm->aerosols) + cpnt; 163 res = suvm_volume_at 164 (aerosol->volume, pos, &cell->prim, cell->barycentric_coords); 165 if(res != RES_OK) { 166 log_err(atm, "Error retrieving %s cell at %g, %g, %g\n", 167 str_cget(&aerosol->name), SPLIT3(pos)); 168 goto error; 169 } 170 171 /* Invalid component */ 172 } else { 173 res = RES_BAD_ARG; 174 goto error; 175 } 176 177 exit: 178 return res; 179 error: 180 goto exit; 181 } 182 183 res_T 184 rnatm_fetch_cell_list 185 (const struct rnatm* atm, 186 const double pos[3], 187 struct rnatm_cell_pos* cells, 188 size_t* ncells) 189 { 190 size_t cpnt = RNATM_GAS; 191 size_t naerosols = 0; 192 size_t i = 0; 193 res_T res = RES_OK; 194 195 if(!atm || !pos || !cells) { 196 res = RES_BAD_ARG; 197 goto error; 198 } 199 200 naerosols = darray_aerosol_size_get(&atm->aerosols); 201 ASSERT(naerosols+1/*gas*/ < RNATM_MAX_COMPONENTS_COUNT); 202 203 do { 204 res = rnatm_fetch_cell(atm, pos, cpnt, cells+i); 205 if(res != RES_OK) goto error; 206 207 ++i; 208 209 } while(++cpnt < naerosols); 210 ASSERT(i == naerosols+1); 211 212 if(ncells) *ncells = naerosols + 1; 213 214 exit: 215 return res; 216 error: 217 goto exit; 218 } 219 220 /******************************************************************************* 221 * Local functions 222 ******************************************************************************/ 223 res_T 224 setup_meshes(struct rnatm* atm, const struct rnatm_create_args* args) 225 { 226 char buf[128]; 227 struct time t0, t1; 228 double gas_low[3]; 229 double gas_upp[3]; 230 struct suvm_device* suvm = NULL; 231 struct smsh_create_args smsh_args = SMSH_CREATE_ARGS_DEFAULT; 232 struct smsh* smsh = NULL; 233 size_t i; 234 res_T res = RES_OK; 235 ASSERT(atm && args); 236 237 log_info(atm, "load and structure the atmosphere meshes\n"); 238 time_current(&t0); 239 240 /* Create the Star-UVM device */ 241 res = suvm_device_create(atm->logger, atm->allocator, atm->verbose, &suvm); 242 if(res != RES_OK) goto error; 243 244 /* Create the Star-Mesh loader */ 245 smsh_args.logger = atm->logger; 246 smsh_args.allocator = atm->allocator; 247 smsh_args.verbose = atm->verbose; 248 res = smsh_create(&smsh_args, &smsh); 249 if(res != RES_OK) goto error; 250 251 /* Load and structure gas volumetric mesh */ 252 log_info(atm, "gas mesh: %s\n", args->gas.smsh_filename); 253 res = setup_uvm(atm, args, args->gas.smsh_filename, suvm, smsh, 254 &atm->gas.volume, &atm->gas.ntetrahedra, &atm->gas.nvertices); 255 if(res != RES_OK) goto error; 256 res = suvm_volume_get_aabb(atm->gas.volume, gas_low, gas_upp); 257 if(res != RES_OK) goto error; 258 259 /* Load and structure aerosol volumetric meshes */ 260 FOR_EACH(i, 0, args->naerosols) { 261 double aerosol_low[3]; 262 double aerosol_upp[3]; 263 struct aerosol* aerosol = darray_aerosol_data_get(&atm->aerosols)+i; 264 const char* aerosol_name = str_cget(&aerosol->name); 265 const char* filename = args->aerosols[i].smsh_filename; 266 267 /* Load and structure the aerosol mesh */ 268 log_info(atm, "%s mesh: %s\n", aerosol_name, filename); 269 res = setup_uvm(atm, args, filename, suvm, smsh, &aerosol->volume, 270 &aerosol->ntetrahedra, &aerosol->nvertices); 271 if(res != RES_OK) goto error; 272 res = suvm_volume_get_aabb(aerosol->volume, aerosol_low, aerosol_upp); 273 if(res != RES_OK) goto error; 274 275 /* Check that the aerosol is included in the gas */ 276 if(gas_low[0] > aerosol_low[0] 277 || gas_low[1] > aerosol_low[1] 278 || gas_low[2] > aerosol_low[2] 279 || gas_upp[0] < aerosol_upp[0] 280 || gas_upp[1] < aerosol_upp[1] 281 || gas_upp[2] < aerosol_upp[2]) { 282 log_err(atm, 283 "The %s may not be included in the gas " 284 "(gas AABB: {%g, %g, %g} - {%g, %g, %g}; " 285 "%s AABB: {%g, %g, %g} - {%g, %g, %g})\n", 286 aerosol_name, 287 SPLIT3(gas_low), 288 SPLIT3(gas_upp), 289 aerosol_name, 290 SPLIT3(aerosol_low), 291 SPLIT3(aerosol_upp)); 292 res = RES_BAD_ARG; 293 goto error; 294 } 295 } 296 297 /* Dump elapsed time */ 298 time_sub(&t0, time_current(&t1), &t0); 299 time_dump(&t0, TIME_ALL, NULL, buf, sizeof(buf)); 300 log_info(atm, "atmopshere meshes setuped in %s\n", buf); 301 302 exit: 303 if(smsh) SMSH(ref_put(smsh)); 304 if(suvm) SUVM(device_ref_put(suvm)); 305 return res; 306 error: 307 if(atm->gas.volume) { 308 SUVM(volume_ref_put(atm->gas.volume)); 309 atm->gas.volume = NULL; 310 } 311 FOR_EACH(i, 0, args->naerosols) { 312 struct aerosol* aerosol = darray_aerosol_data_get(&atm->aerosols)+i; 313 if(aerosol->volume) { 314 SUVM(volume_ref_put(aerosol->volume)); 315 aerosol->volume = NULL; 316 } 317 } 318 goto exit; 319 }