schiff_mesh.c (18687B)
1 /* Copyright (C) 2015, 2016, 2026 Centre National de la Recherche Scientifique 2 * Copyright (C) 2026 Clermont Auvergne INP 3 * Copyright (C) 2026 Institut Mines Télécom Albi-Carmaux 4 * Copyright (C) 2017, 2019-2021, 2026 |Méso|Star> (contact@meso-star.com) 5 * Copyright (C) 2026 Université de Lorraine 6 * Copyright (C) 2026 Université de Toulouse 7 * 8 * This program is free software: you can redistribute it and/or modify 9 * it under the terms of the GNU General Public License as published by 10 * the Free Software Foundation, either version 3 of the License, or 11 * (at your option) any later version. 12 * 13 * This program is distributed in the hope that it will be useful, 14 * but WITHOUT ANY WARRANTY; without even the implied warranty of 15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 16 * GNU General Public License for more details. 17 * 18 * You should have received a copy of the GNU General Public License 19 * along with this program. If not, see <http://www.gnu.org/licenses/>. */ 20 21 #include "schiff_mesh.h" 22 23 #include <rsys/float2.h> 24 #include <rsys/float3.h> 25 #include <rsys/stretchy_array.h> 26 27 /******************************************************************************* 28 * Helper functions 29 ******************************************************************************/ 30 /* Assume that the vertices are arranged in "phi major" order and that the 31 * lower and upper polar vertices are the two last ones. */ 32 static void 33 setup_sphere_indices 34 (struct darray_uint* indices, 35 const unsigned nthetas, 36 const unsigned nphis) 37 { 38 unsigned i, itheta, iphi; 39 ASSERT(indices && nthetas && nphis); 40 41 /* Define the indices of the contour primitives */ 42 i = 0; 43 FOR_EACH(itheta, 0, nthetas) { 44 const unsigned itheta0 = itheta * (nphis - 1); 45 const unsigned itheta1 = ((itheta + 1) % nthetas) * (nphis - 1); 46 FOR_EACH(iphi, 0, nphis-2) { 47 const unsigned iphi0 = iphi + 0; 48 const unsigned iphi1 = iphi + 1; 49 50 darray_uint_data_get(indices)[i++] = itheta0 + iphi0; 51 darray_uint_data_get(indices)[i++] = itheta0 + iphi1; 52 darray_uint_data_get(indices)[i++] = itheta1 + iphi0; 53 54 darray_uint_data_get(indices)[i++] = itheta1 + iphi0; 55 darray_uint_data_get(indices)[i++] = itheta0 + iphi1; 56 darray_uint_data_get(indices)[i++] = itheta1 + iphi1; 57 } 58 } 59 60 /* Define the indices of the polar primitives */ 61 FOR_EACH(itheta, 0, nthetas) { 62 const unsigned itheta0 = itheta * (nphis - 1); 63 const unsigned itheta1 = ((itheta + 1) % nthetas) * (nphis - 1); 64 65 darray_uint_data_get(indices)[i++] = nthetas * (nphis - 1); 66 darray_uint_data_get(indices)[i++] = itheta0; 67 darray_uint_data_get(indices)[i++] = itheta1; 68 69 darray_uint_data_get(indices)[i++] = nthetas * (nphis - 1) + 1; 70 darray_uint_data_get(indices)[i++] = itheta1 + (nphis - 2); 71 darray_uint_data_get(indices)[i++] = itheta0 + (nphis - 2); 72 } 73 } 74 75 /******************************************************************************* 76 * Local functions 77 ******************************************************************************/ 78 res_T 79 schiff_mesh_init_sphere 80 (struct mem_allocator* allocator, 81 struct schiff_mesh* sphere, 82 const unsigned nthetas) 83 { 84 /* Theta in [0, 2PI[ and phi in [0, PI[ */ 85 const unsigned nphis = (unsigned)(((double)nthetas + 0.5) / 2.0); 86 const double step_theta = 2*PI / (double)nthetas; 87 const double step_phi = PI / (double)nphis; 88 size_t nverts; 89 size_t ntris; 90 size_t i; 91 unsigned itheta, iphi; 92 res_T res = RES_OK; 93 ASSERT(allocator && sphere && nthetas); 94 95 sphere->coordinates = SCHIFF_CARTESIAN; 96 darray_float_init(allocator, &sphere->vertices.cartesian); 97 darray_uint_init(allocator, &sphere->indices); 98 darray_sincos_init(allocator, &sphere->thetas); 99 darray_sincos_init(allocator, &sphere->phis); 100 101 nverts = nthetas*(nphis-1)/* #contour verts */ + 2/* polar verts */; 102 ntris = 2*nthetas*(nphis-2)/* #contour tris */ + 2*nthetas/* #polar tris */; 103 104 res = darray_float_resize(&sphere->vertices.cartesian, nverts * 3/*#coords*/); 105 if(res != RES_OK) goto error; 106 res = darray_uint_resize(&sphere->indices, ntris * 3/*#indices per tri*/); 107 if(res != RES_OK) goto error; 108 res = darray_sincos_resize(&sphere->thetas, nthetas); 109 if(res != RES_OK) goto error; 110 res = darray_sincos_resize(&sphere->phis, nphis); 111 if(res != RES_OK) goto error; 112 113 /* Precompute the cosine/sinus of the theta/phi angles */ 114 FOR_EACH(itheta, 0, nthetas) { 115 const double theta = -PI + (double)itheta * step_theta; 116 darray_sincos_data_get(&sphere->thetas)[itheta].angle = theta; 117 darray_sincos_data_get(&sphere->thetas)[itheta].sinus = sin(theta); 118 darray_sincos_data_get(&sphere->thetas)[itheta].cosine = cos(theta); 119 } 120 FOR_EACH(iphi, 0, nphis-1) { 121 const double phi = -PI/2 + (double)(iphi + 1) * step_phi; 122 darray_sincos_data_get(&sphere->phis)[iphi].angle = phi; 123 darray_sincos_data_get(&sphere->phis)[iphi].sinus = sin(phi); 124 darray_sincos_data_get(&sphere->phis)[iphi].cosine = cos(phi); 125 } 126 127 /* Build the contour vertices */ 128 i = 0; 129 FOR_EACH(itheta, 0, nthetas) { 130 const struct sin_cos* theta = darray_sincos_data_get(&sphere->thetas) + itheta; 131 FOR_EACH(iphi, 0, nphis-1) { 132 const struct sin_cos* phi = darray_sincos_data_get(&sphere->phis) + iphi; 133 darray_float_data_get(&sphere->vertices.cartesian)[i++] = 134 (float)(theta->cosine * phi->cosine); 135 darray_float_data_get(&sphere->vertices.cartesian)[i++] = 136 (float)(theta->sinus * phi->cosine); 137 darray_float_data_get(&sphere->vertices.cartesian)[i++] = 138 (float)phi->sinus; 139 } 140 } 141 /* Setup polar vertices */ 142 f3(darray_float_data_get(&sphere->vertices.cartesian) + i+0, 0.f, 0.f,-1.f); 143 f3(darray_float_data_get(&sphere->vertices.cartesian) + i+3, 0.f, 0.f, 1.f); 144 145 /* Define the indices of the sphere */ 146 setup_sphere_indices(&sphere->indices, nthetas, nphis); 147 148 sphere->is_init = 1; 149 exit: 150 darray_sincos_release(&sphere->thetas); 151 darray_sincos_release(&sphere->phis); 152 return res; 153 error: 154 schiff_mesh_release(sphere); 155 goto exit; 156 } 157 158 res_T 159 schiff_mesh_init_sphere_polar 160 (struct mem_allocator* allocator, 161 struct schiff_mesh* sphere, 162 const unsigned nthetas) 163 { 164 /* Theta in [0, 2PI[ and phi in [0, PI[ */ 165 const unsigned nphis = (unsigned)(((double)nthetas + 0.5) / 2.0); 166 const double step_theta = 2*PI / (double)nthetas; 167 const double step_phi = PI / (double)nphis; 168 size_t nverts; 169 size_t ntris; 170 size_t i; 171 unsigned itheta, iphi; 172 res_T res = RES_OK; 173 ASSERT(allocator && sphere && nthetas); 174 175 sphere->coordinates = SCHIFF_POLAR; 176 darray_uint_init(allocator, &sphere->vertices.polar); 177 darray_uint_init(allocator, &sphere->indices); 178 darray_sincos_init(allocator, &sphere->thetas); 179 darray_sincos_init(allocator, &sphere->phis); 180 181 nverts = nthetas*(nphis-1)/* #contour verts */ + 2/* polar verts */; 182 ntris = 2*nthetas*(nphis-2)/* #contour tris */ + 2*nthetas/* #polar tris */; 183 184 res = darray_uint_resize(&sphere->vertices.polar, nverts * 2/* theta phi*/); 185 if(res != RES_OK) goto error; 186 res = darray_uint_resize(&sphere->indices, ntris * 3/*#indices per tri*/); 187 if(res != RES_OK) goto error; 188 res = darray_sincos_resize(&sphere->thetas, nthetas); 189 if(res != RES_OK) goto error; 190 res = darray_sincos_resize(&sphere->phis, nphis + 1/*Polar*/); 191 if(res != RES_OK) goto error; 192 193 /* Cosine/Sinus of the theta/phi contour angles */ 194 FOR_EACH(itheta, 0, nthetas) { 195 const double theta = -PI + (double)itheta * step_theta; 196 darray_sincos_data_get(&sphere->thetas)[itheta].angle = theta; 197 darray_sincos_data_get(&sphere->thetas)[itheta].sinus = (float)sin(theta); 198 darray_sincos_data_get(&sphere->thetas)[itheta].cosine = (float)cos(theta); 199 } 200 FOR_EACH(iphi, 0, nphis-1) { 201 const double phi = -PI/2 + (double)(iphi + 1) * step_phi; 202 darray_sincos_data_get(&sphere->phis)[iphi].angle = (float)phi; 203 darray_sincos_data_get(&sphere->phis)[iphi].sinus = (float)sin(phi); 204 darray_sincos_data_get(&sphere->phis)[iphi].cosine = (float)cos(phi); 205 } 206 207 /* Cosine/Sinus of the theta/phi polar angles */ 208 darray_sincos_data_get(&sphere->phis)[nphis-1].angle =-(float)PI/2; 209 darray_sincos_data_get(&sphere->phis)[nphis-1].sinus =-1.f; 210 darray_sincos_data_get(&sphere->phis)[nphis-1].cosine = 0.f; 211 darray_sincos_data_get(&sphere->phis)[nphis-0].angle = (float)PI/2; 212 darray_sincos_data_get(&sphere->phis)[nphis-0].sinus = 1.f; 213 darray_sincos_data_get(&sphere->phis)[nphis-0].cosine = 0.f; 214 215 /* Build the contour vertices */ 216 i = 0; 217 FOR_EACH(itheta, 0, nthetas) { 218 FOR_EACH(iphi, 0, nphis-1) { 219 darray_uint_data_get(&sphere->vertices.polar)[i++] = itheta; 220 darray_uint_data_get(&sphere->vertices.polar)[i++] = iphi; 221 } 222 } 223 /* Setup polar vertices */ 224 darray_uint_data_get(&sphere->vertices.polar)[i++] = 0; 225 darray_uint_data_get(&sphere->vertices.polar)[i++] = nphis-1; 226 darray_uint_data_get(&sphere->vertices.polar)[i++] = 0; 227 darray_uint_data_get(&sphere->vertices.polar)[i++] = nphis-0; 228 229 /* Define the indices of the sphere */ 230 setup_sphere_indices(&sphere->indices, nthetas, nphis); 231 232 sphere->is_init = 1; 233 exit: 234 return res; 235 error: 236 schiff_mesh_release(sphere); 237 goto exit; 238 } 239 240 res_T 241 schiff_mesh_init_cylinder 242 (struct mem_allocator* allocator, 243 struct schiff_mesh* cylinder, 244 const unsigned nsteps) 245 { 246 const double step = 2*PI / (double)nsteps; 247 size_t nverts; 248 size_t ntris; 249 unsigned i; 250 res_T res = RES_OK; 251 ASSERT(allocator && cylinder && nsteps); 252 253 cylinder->coordinates = SCHIFF_CARTESIAN; 254 darray_float_init(allocator, &cylinder->vertices.cartesian); 255 darray_uint_init(allocator, &cylinder->indices); 256 257 nverts = nsteps*2/* #contour verts */ + 2/* #polar verts */; 258 ntris = nsteps*2/* #contour tris */ + 2*nsteps/* #cap tris */; 259 260 res = darray_float_resize(&cylinder->vertices.cartesian, nverts*3/*#coords*/); 261 if(res != RES_OK) goto error; 262 res = darray_uint_resize(&cylinder->indices, ntris*3/*#indices per tri*/); 263 if(res != RES_OK) goto error; 264 265 /* Generate the vertex coordinates */ 266 FOR_EACH(i, 0, nsteps) { 267 const float theta = (float)(i* step); 268 const float x = (float)cos(theta); 269 const float y = (float)sin(theta); 270 f3(darray_float_data_get(&cylinder->vertices.cartesian)+(i*2+0)*3, x,y,0); 271 f3(darray_float_data_get(&cylinder->vertices.cartesian)+(i*2+1)*3, x,y,1); 272 } 273 274 /* "Polar" vertices */ 275 f3(darray_float_data_get(&cylinder->vertices.cartesian)+(i*2+0)*3, 0,0,0); 276 f3(darray_float_data_get(&cylinder->vertices.cartesian)+(i*2+1)*3, 0,0,1); 277 278 /* Contour primitives */ 279 FOR_EACH(i, 0, nsteps) { 280 const unsigned id = i * 2; 281 unsigned* iprim = darray_uint_data_get(&cylinder->indices) + id*3; 282 283 iprim[0] = (id + 0); 284 iprim[1] = (id + 1); 285 iprim[2] = (id + 2) % (nsteps*2); 286 287 iprim += 3; 288 289 iprim[0] = (id + 2) % (nsteps*2); 290 iprim[1] = (id + 1); 291 iprim[2] = (id + 3) % (nsteps*2); 292 } 293 294 /* Cap primitives */ 295 FOR_EACH(i, 0, nsteps) { 296 const unsigned id = i* 2; 297 unsigned* iprim = darray_uint_data_get(&cylinder->indices) + (nsteps + i)*6; 298 299 iprim[0] = (nsteps * 2); 300 iprim[1] = (id + 0); 301 iprim[2] = (id + 2) % (nsteps*2); 302 303 iprim += 3; 304 305 iprim[0] = (nsteps * 2) + 1; 306 iprim[1] = (id + 3) % (nsteps*2); 307 iprim[2] = (id + 1); 308 } 309 310 cylinder->is_init = 1; 311 exit: 312 return res; 313 error: 314 schiff_mesh_release(cylinder); 315 goto exit; 316 } 317 318 res_T 319 schiff_mesh_init_helical_pipe 320 (struct mem_allocator* allocator, 321 struct schiff_mesh* helical_pipe, 322 const unsigned nsteps_helicoid, 323 const unsigned nsteps_circle) 324 325 { 326 unsigned* indices; 327 unsigned* bottom_cap; 328 unsigned* top_cap; 329 size_t ihelicoid, icircle; 330 size_t nverts, ntris; 331 size_t ibottom_cap_index, itop_cap_index; 332 size_t ibottom_cap_vertex, itop_cap_vertex; 333 size_t ilast_meridian_vertex; 334 res_T res; 335 ASSERT(allocator && helical_pipe && nsteps_helicoid>=2 && nsteps_circle>=4); 336 337 helical_pipe->coordinates = SCHIFF_NO_COORDINATE; 338 darray_uint_init(allocator, &helical_pipe->indices); 339 340 nverts = (nsteps_helicoid+1)*nsteps_circle/*#contour*/ + 2/*#cap*/; 341 ntris = (nsteps_helicoid)*nsteps_circle*2/*#contour*/ + 2*nsteps_circle/*#cap*/; 342 343 res = darray_uint_resize(&helical_pipe->indices, ntris*3/*#indices per tri*/); 344 if(res != RES_OK) goto error; 345 346 indices = darray_uint_data_get(&helical_pipe->indices); 347 348 /* Setup the indices toward the contour primitives */ 349 FOR_EACH(ihelicoid, 0, nsteps_helicoid) { 350 /* Offset toward the first index of the current meridian */ 351 const size_t offset = ihelicoid * nsteps_circle; 352 353 FOR_EACH(icircle, 0, nsteps_circle) { 354 size_t prim_ids[4]; 355 const size_t id = icircle + offset; 356 unsigned* iprim = indices + id*3/*#ids per prim*/*2/*#prims per step*/; 357 358 /* Indices of the quad */ 359 prim_ids[0] = icircle + offset; 360 prim_ids[1] = icircle + nsteps_circle + offset; 361 prim_ids[2] = (icircle + 1) % nsteps_circle + nsteps_circle + offset; 362 prim_ids[3] = (icircle + 1) % nsteps_circle + offset; 363 /* First triangle of the quad */ 364 iprim[0] = (unsigned)prim_ids[0]; 365 iprim[1] = (unsigned)prim_ids[3]; 366 iprim[2] = (unsigned)prim_ids[1]; 367 /* Second triangule of the quad */ 368 iprim[3] = (unsigned)prim_ids[1]; 369 iprim[4] = (unsigned)prim_ids[3]; 370 iprim[5] = (unsigned)prim_ids[2]; 371 } 372 } 373 374 /* Define the index of the cap vertices */ 375 ibottom_cap_vertex = nverts - 2; /* Index toward the bottom cap vertex */ 376 itop_cap_vertex = nverts - 1; /* Index toward the top cap vertex */ 377 378 /* Index of the first vertex of the last meridian circle */ 379 ilast_meridian_vertex = nsteps_helicoid * nsteps_circle; 380 381 ibottom_cap_index = /* Index of the 1st bottom cap index */ 382 nsteps_helicoid 383 * nsteps_circle 384 * 2 /*#prims per step*/ 385 * 3 /*#ids per prim*/; 386 itop_cap_index = /* Index of the 1st top cap index */ 387 ibottom_cap_index 388 + nsteps_circle * 3/*# ids per primitive */; 389 390 bottom_cap = indices + ibottom_cap_index ; 391 top_cap = indices + itop_cap_index; 392 393 FOR_EACH(icircle, 0, nsteps_circle) { 394 const size_t id = icircle * 3/*#ids per prim*/; 395 396 bottom_cap[id + 0] = (unsigned)icircle; 397 bottom_cap[id + 1] = (unsigned)ibottom_cap_vertex; 398 bottom_cap[id + 2] = (unsigned)((icircle + 1) % nsteps_circle); 399 400 top_cap[id + 0] = (unsigned)((icircle+1)%nsteps_circle + ilast_meridian_vertex); 401 top_cap[id + 1] = (unsigned)itop_cap_vertex; 402 top_cap[id + 2] = (unsigned)(icircle + ilast_meridian_vertex); 403 } 404 405 helical_pipe->is_init = 1; 406 exit: 407 return res; 408 error: 409 schiff_mesh_release(helical_pipe); 410 goto exit; 411 } 412 413 res_T 414 schiff_mesh_helical_pipe_create_vertices 415 (const struct schiff_mesh* mesh, 416 const double pitch, 417 const double height, 418 const double hradius, /* Helicoid radius */ 419 const double cradius, /* Circle radius */ 420 const unsigned nsteps_helicoid, 421 const unsigned nsteps_circle, 422 size_t* nvertices, 423 float** out_vertices) 424 { 425 double* meridian = NULL; 426 float* vertices = NULL; 427 double c; 428 double phi_max; 429 double step_helicoid, step_circle; 430 double rcp_sqrt_hradius2_add_c2; /* Precomputed value */ 431 size_t ihelicoid, icircle; 432 size_t icoord, icoord_top, icoord_bottom; 433 size_t nverts; 434 res_T res = RES_OK; 435 436 ASSERT(mesh && pitch > 0 && height > 0 && hradius > 0 && cradius > 0); 437 ASSERT(nvertices && out_vertices); 438 ASSERT(mesh->coordinates == SCHIFF_NO_COORDINATE); 439 ASSERT(nsteps_helicoid * nsteps_circle * 2 * 3 + nsteps_circle*3*2/*cap*/ 440 == darray_uint_size_get(&mesh->indices)); 441 (void)mesh; 442 443 nverts = (nsteps_helicoid+1) * nsteps_circle/*#contour*/ + 2/*#cap*/; 444 if(!sa_add(vertices, nverts * 3/*#coords per vertex*/)) { 445 res = RES_MEM_ERR; 446 goto error; 447 } 448 if(!sa_add(meridian, nsteps_circle * 3/*#coords per vertex*/)) { 449 res = RES_MEM_ERR; 450 goto error; 451 } 452 453 c = pitch / (2*PI); 454 phi_max = height * 2*PI / pitch; 455 step_helicoid = phi_max / (double)nsteps_helicoid; 456 step_circle = 2*PI / (double)nsteps_circle; 457 rcp_sqrt_hradius2_add_c2 = 1.0 / sqrt(hradius*hradius + c*c); 458 459 /* Compute the meridian positions */ 460 FOR_EACH(icircle, 0, nsteps_circle) { 461 const double theta = (double)icircle * step_circle; 462 const double cos_theta = cos(theta); 463 const double sin_theta = sin(theta); 464 const size_t id = icircle *3; 465 466 meridian[id + 0] = cradius * cos_theta + hradius; 467 meridian[id + 1] = -cradius * c * rcp_sqrt_hradius2_add_c2 * sin_theta; 468 meridian[id + 2] = cradius*hradius * rcp_sqrt_hradius2_add_c2 * sin_theta; 469 } 470 471 icoord = 0; 472 icoord_bottom = 473 (nsteps_helicoid + 1) 474 * nsteps_circle 475 * 3/*#coords per vertex*/; 476 icoord_top = icoord_bottom + 3/*#coords per vertex */; 477 478 FOR_EACH(ihelicoid, 0, nsteps_helicoid + 1) { 479 const double phi = (double)ihelicoid * step_helicoid; 480 const double cos_phi = cos(phi); 481 const double sin_phi = sin(phi); 482 483 FOR_EACH(icircle, 0, nsteps_circle) { 484 const double* pos = meridian + icircle*3; 485 vertices[icoord + 0] = (float)(pos[0]*cos_phi - pos[1]*sin_phi); 486 vertices[icoord + 1] = (float)(pos[0]*sin_phi + pos[1]*cos_phi); 487 vertices[icoord + 2] = (float)(pos[2] + c * phi); 488 icoord += 3; 489 } 490 491 /* Top cap vertex */ 492 if(ihelicoid == 0) { 493 vertices[icoord_bottom + 0] = (float)(hradius * cos_phi); 494 vertices[icoord_bottom + 1] = (float)(hradius * sin_phi); 495 vertices[icoord_bottom + 2] = (float)(c*phi); 496 497 /* Bottom cap vertex */ 498 } else if(ihelicoid == nsteps_helicoid) { 499 vertices[icoord_top + 0] = (float)(hradius * cos_phi); 500 vertices[icoord_top + 1] = (float)(hradius * sin_phi); 501 vertices[icoord_top + 2] = (float)(c * phi); 502 } 503 } 504 505 exit: 506 if(meridian) sa_release(meridian); 507 *nvertices = nverts; 508 *out_vertices = vertices; 509 return res; 510 error: 511 if(vertices) { 512 sa_release(vertices); 513 vertices = NULL; 514 nverts = 0; 515 } 516 goto exit; 517 } 518 519 void 520 schiff_mesh_helical_pipe_destroy_vertices 521 (const struct schiff_mesh* mesh, float* out_vertices) 522 { 523 ASSERT(mesh && out_vertices); 524 (void)mesh; 525 sa_release(out_vertices); 526 } 527 528 void 529 schiff_mesh_release(struct schiff_mesh* mesh) 530 { 531 ASSERT(mesh); 532 if(!mesh->is_init) 533 return; /* The mesh is not initialised */ 534 535 switch(mesh->coordinates) { 536 case SCHIFF_CARTESIAN: 537 darray_float_release(&mesh->vertices.cartesian); 538 break; 539 case SCHIFF_POLAR: 540 darray_uint_release(&mesh->vertices.polar); 541 darray_sincos_release(&mesh->thetas); 542 darray_sincos_release(&mesh->phis); 543 break; 544 case SCHIFF_NO_COORDINATE: /* Do nothing */ break; 545 default: FATAL("Unreachable code\n"); break; 546 } 547 darray_uint_release(&mesh->indices); 548 mesh->coordinates = SCHIFF_NO_COORDINATE; 549 mesh->is_init = 0; 550 } 551