mmg3d
inlined_functions_3d.h
Go to the documentation of this file.
1/* =============================================================================
2** This file is part of the mmg software package for the tetrahedral
3** mesh modification.
4** Copyright (c) Bx INP/CNRS/Inria/UBordeaux/UPMC, 2004-
5**
6** mmg is free software: you can redistribute it and/or modify it
7** under the terms of the GNU Lesser General Public License as published
8** by the Free Software Foundation, either version 3 of the License, or
9** (at your option) any later version.
10**
11** mmg is distributed in the hope that it will be useful, but WITHOUT
12** ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13** FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
14** License for more details.
15**
16** You should have received a copy of the GNU Lesser General Public
17** License and of the GNU General Public License along with mmg (in
18** files COPYING.LESSER and COPYING). If not, see
19** <http://www.gnu.org/licenses/>. Please read their terms carefully and
20** use this copy of the mmg distribution only if you accept them.
21** =============================================================================
22*/
23
35#include "mmg3d.h"
36#include "inlined_functions.h"
37
38#ifndef _INLINED_FUNCT_3D_H
39#define _INLINED_FUNCT_3D_H
40
41
55static
56inline double MMG5_lenedgCoor_ani(double *ca,double *cb,double *sa,double *sb) {
57 double ux,uy,uz,dd1,dd2,len;
58
59 ux = cb[0] - ca[0];
60 uy = cb[1] - ca[1];
61 uz = cb[2] - ca[2];
62
63 dd1 = sa[0]*ux*ux + sa[3]*uy*uy + sa[5]*uz*uz \
64 + 2.0*(sa[1]*ux*uy + sa[2]*ux*uz + sa[4]*uy*uz);
65 if ( dd1 <= 0.0 ) dd1 = 0.0;
66
67 dd2 = sb[0]*ux*ux + sb[3]*uy*uy + sb[5]*uz*uz \
68 + 2.0*(sb[1]*ux*uy + sb[2]*ux*uz + sb[4]*uy*uz);
69 if ( dd2 <= 0.0 ) dd2 = 0.0;
70
71 /*longueur approchee*/
72 /*precision a 3.5 10e-3 pres*/
73 if(fabs(dd1-dd2) < 0.05 ) {
74 len = sqrt(0.5*(dd1+dd2));
75 return len;
76 }
77 len = (sqrt(dd1)+sqrt(dd2)+4.0*sqrt(0.5*(dd1+dd2))) / 6.0;
78
79 return len;
80}
81
93static
94inline double MMG5_lenedg33_ani(MMG5_pMesh mesh ,MMG5_pSol met, int ia,
95 MMG5_pTetra pt)
96{
97 int ip1,ip2;
98 int8_t isedg;
99
100 ip1 = pt->v[MMG5_iare[ia][0]];
101 ip2 = pt->v[MMG5_iare[ia][1]];
102
103 if ( pt->xt && (mesh->xtetra[pt->xt].tag[ia] & MG_BDY)) {
104 isedg = ( mesh->xtetra[pt->xt].tag[ia] & MG_GEO);
105 return MMG5_lenSurfEdg33_ani(mesh, met, ip1, ip2, isedg);
106 } else {
107 return MMG5_lenedgCoor_ani(mesh->point[ip1].c,mesh->point[ip2].c,
108 &met->m[6*ip1],&met->m[6*ip2]);
109 }
110 return 0.0;
111}
112
124static
126 MMG5_pTetra pt)
127{
128 MMG5_pPoint pp1,pp2;
129 double *m1,*m2;
130 int ip1,ip2;
131
132 ip1 = pt->v[MMG5_iare[ia][0]];
133 ip2 = pt->v[MMG5_iare[ia][1]];
134
135 pp1 = &mesh->point[ip1];
136 pp2 = &mesh->point[ip2];
137
138 m1 = &met->m[6*ip1];
139 m2 = &met->m[6*ip2];
140
141 return MMG5_lenedgCoor_ani(pp1->c,pp2->c,m1,m2);
142}
143
144
145
157static
158inline double MMG5_lenedgspl_ani(MMG5_pMesh mesh ,MMG5_pSol met, int ia,
159 MMG5_pTetra pt)
160{
161 MMG5_pPoint pp1,pp2;
162 double m1[6],m2[6];
163 int ip1,ip2,i;
164
165 ip1 = pt->v[MMG5_iare[ia][0]];
166 ip2 = pt->v[MMG5_iare[ia][1]];
167
168 pp1 = &mesh->point[ip1];
169 pp2 = &mesh->point[ip2];
170
171 if ( !(MG_SIN(pp1->tag) || (MG_NOM & pp1->tag)) && (pp1->tag & MG_GEO) ) {
172 if ( !MMG5_moymet(mesh,met,pt,m1) ) return 0;
173 } else {
174 for ( i=0; i<6; ++i )
175 m1[i] = met->m[6*ip1+i];
176 }
177
178 if ( !(MG_SIN(pp2->tag)|| (MG_NOM & pp2->tag)) && (pp2->tag & MG_GEO) ) {
179 if ( !MMG5_moymet(mesh,met,pt,m2) ) return 0;
180 } else {
181 for ( i=0; i<6; ++i )
182 m2[i] = met->m[6*ip2+i];
183 }
184
185 return MMG5_lenedgCoor_ani(pp1->c,pp2->c,m1,m2);
186}
187
199static
200inline double MMG5_lenedg_ani(MMG5_pMesh mesh ,MMG5_pSol met, int ia,
201 MMG5_pTetra pt)
202{
203 int ip1,ip2;
204 int8_t isedg;
205
206 ip1 = pt->v[MMG5_iare[ia][0]];
207 ip2 = pt->v[MMG5_iare[ia][1]];
208
209 if ( pt->xt && (mesh->xtetra[pt->xt].tag[ia] & MG_BDY)) {
210 isedg = ( mesh->xtetra[pt->xt].tag[ia] & MG_GEO);
211 return MMG5_lenSurfEdg_ani(mesh, met, ip1, ip2, isedg);
212 } else {
213 return MMG5_lenedgspl_ani(mesh ,met, ia, pt);
214 }
215 return 0.0;
216}
217
229static
230inline double MMG5_lenedg_iso(MMG5_pMesh mesh,MMG5_pSol met,int ia,
231 MMG5_pTetra pt) {
232 int ip1,ip2;
233
234 ip1 = pt->v[MMG5_iare[ia][0]];
235 ip2 = pt->v[MMG5_iare[ia][1]];
236 return MMG5_lenSurfEdg_iso(mesh,met,ip1,ip2,0);
237}
238
239static
240inline double MMG5_lenedgspl_iso(MMG5_pMesh mesh ,MMG5_pSol met, int ia,
241 MMG5_pTetra pt) {
242 int ip1,ip2;
243
244 ip1 = pt->v[MMG5_iare[ia][0]];
245 ip2 = pt->v[MMG5_iare[ia][1]];
246
247 return MMG5_lenSurfEdg_iso(mesh,met,ip1,ip2,0);
248
249}
250
251
261static
262inline double MMG5_orcal(MMG5_pMesh mesh,MMG5_pSol met,int iel) {
263 MMG5_pTetra pt;
264
265 pt = &mesh->tetra[iel];
266
267 return MMG5_caltet(mesh,met,pt);
268}
269
270
282static
284 double ct[12],cs[3],rad,Vref,V,cal;
285 int j,l;
286
287 for (j=0,l=0; j<4; j++,l+=3) {
288 memcpy(&ct[l],mesh->point[pt->v[j]].c,3*sizeof(double));
289 }
290
291 if(!MMG5_cenrad_iso(mesh,ct,cs,&rad)) {
292 return 0.0;
293 }
294
295 assert(rad>0.);
296
297 /* Vref volume */
298 Vref = 8.*sqrt(3)/27.*rad*sqrt(rad);
299
300 V = MMG5_orvol(mesh->point,pt->v) * MMG3D_DET2VOL;
301
302 if ( V<0. ) {
303 return 0.0;
304 }
305
306 cal = V/Vref; //1-Qles in order to have the best quality equal to 1
307
308 /* Prevent calities > 1 due to the machin epsilon */
309 cal = MG_MIN (1., cal);
310
311 // the normalization by ALPHAD
312 //in order to be coherent with the other quality measure
313 return cal/MMG3D_ALPHAD;
314}
315
327static
328inline double MMG5_caltet_iso_4pt(double *a, double *b, double *c, double *d) {
329 double abx,aby,abz,acx,acy,acz,adx,ady,adz,bcx,bcy,bcz,bdx,bdy,bdz;
330 double cdx,cdy,cdz;
331 double vol,v1,v2,v3,rap;
332
333 /* volume: (ac^ad).ab/6. Note that here we compute 6*vol. */
334 abx = b[0] - a[0];
335 aby = b[1] - a[1];
336 abz = b[2] - a[2];
337 rap = abx*abx + aby*aby + abz*abz;
338
339 acx = c[0] - a[0];
340 acy = c[1] - a[1];
341 acz = c[2] - a[2];
342 rap += acx*acx + acy*acy + acz*acz;
343
344 adx = d[0] - a[0];
345 ady = d[1] - a[1];
346 adz = d[2] - a[2];
347 rap += adx*adx + ady*ady + adz*adz;
348
349 v1 = acy*adz - acz*ady;
350 v2 = acz*adx - acx*adz;
351 v3 = acx*ady - acy*adx;
352 vol = abx * v1 + aby * v2 + abz * v3;
353 if ( vol < MMG5_EPSD2 ) return 0.0;
354
355 bcx = c[0] - b[0];
356 bcy = c[1] - b[1];
357 bcz = c[2] - b[2];
358 rap += bcx*bcx + bcy*bcy + bcz*bcz;
359
360 bdx = d[0] - b[0];
361 bdy = d[1] - b[1];
362 bdz = d[2] - b[2];
363 rap += bdx*bdx + bdy*bdy + bdz*bdz;
364
365 cdx = d[0] - c[0];
366 cdy = d[1] - c[1];
367 cdz = d[2] - c[2];
368 rap += cdx*cdx + cdy*cdy + cdz*cdz;
369 if ( rap < MMG5_EPSD2 ) return 0.0;
370
371 /* quality = 6*vol / len^3/2. Q = 1/(12 sqrt(3)) for the regular tetra of length 1. */
372 rap = rap * sqrt(rap);
373 return vol / rap;
374}
375
386static
388 double *a,*b,*c,*d;
389 int ia, ib, ic, id;
390
391 ia = pt->v[0];
392 ib = pt->v[1];
393 ic = pt->v[2];
394 id = pt->v[3];
395
396 a = mesh->point[ia].c;
397 b = mesh->point[ib].c;
398 c = mesh->point[ic].c;
399 d = mesh->point[id].c;
400
401 return MMG5_caltet_iso_4pt(a,b,c,d);
402
403}
404
416static
418 double cal,abx,aby,abz,acx,acy,acz,adx,ady,adz;
419 double h1,h2,h3,h4,h5,h6,det,vol,rap,v1,v2,v3,num;
420 double bcx,bcy,bcz,bdx,bdy,bdz,cdx,cdy,cdz;
421 double *a,*b,*c,*d;
422 double mm[6];
423 int ip[4];
424
425 ip[0] = pt->v[0];
426 ip[1] = pt->v[1];
427 ip[2] = pt->v[2];
428 ip[3] = pt->v[3];
429
430 /* average metric */
431 if ( !MMG5_moymet(mesh,met,pt,&mm[0]) )
432 return (0.0);
433
434 a = mesh->point[ip[0]].c;
435 b = mesh->point[ip[1]].c;
436 c = mesh->point[ip[2]].c;
437 d = mesh->point[ip[3]].c;
438
439 /* volume */
440 abx = b[0] - a[0];
441 aby = b[1] - a[1];
442 abz = b[2] - a[2];
443
444 acx = c[0] - a[0];
445 acy = c[1] - a[1];
446 acz = c[2] - a[2];
447
448 adx = d[0] - a[0];
449 ady = d[1] - a[1];
450 adz = d[2] - a[2];
451
452 bcx = c[0] - b[0];
453 bcy = c[1] - b[1];
454 bcz = c[2] - b[2];
455
456 bdx = d[0] - b[0];
457 bdy = d[1] - b[1];
458 bdz = d[2] - b[2];
459
460 cdx = d[0] - c[0];
461 cdy = d[1] - c[1];
462 cdz = d[2] - c[2];
463
464 v1 = acy*adz - acz*ady;
465 v2 = acz*adx - acx*adz;
466 v3 = acx*ady - acy*adx;
467 vol = abx * v1 + aby * v2 + abz * v3;
468 if ( vol <= 0. ) return 0.0;
469
470 det = mm[0] * ( mm[3]*mm[5] - mm[4]*mm[4]) \
471 - mm[1] * ( mm[1]*mm[5] - mm[2]*mm[4]) \
472 + mm[2] * ( mm[1]*mm[4] - mm[2]*mm[3]);
473 if ( det < MMG5_EPSD2 ) {
474 return 0.0;
475 }
476 det = sqrt(det) * vol;
477
478 /* edge lengths */
479 h1 = mm[0]*abx*abx + mm[3]*aby*aby + mm[5]*abz*abz
480 + 2.0*(mm[1]*abx*aby + mm[2]*abx*abz + mm[4]*aby*abz);
481 h2 = mm[0]*acx*acx + mm[3]*acy*acy + mm[5]*acz*acz
482 + 2.0*(mm[1]*acx*acy + mm[2]*acx*acz + mm[4]*acy*acz);
483 h3 = mm[0]*adx*adx + mm[3]*ady*ady + mm[5]*adz*adz
484 + 2.0*(mm[1]*adx*ady + mm[2]*adx*adz + mm[4]*ady*adz);
485 h4 = mm[0]*bcx*bcx + mm[3]*bcy*bcy + mm[5]*bcz*bcz
486 + 2.0*(mm[1]*bcx*bcy + mm[2]*bcx*bcz + mm[4]*bcy*bcz);
487 h5 = mm[0]*bdx*bdx + mm[3]*bdy*bdy + mm[5]*bdz*bdz
488 + 2.0*(mm[1]*bdx*bdy + mm[2]*bdx*bdz + mm[4]*bdy*bdz);
489 h6 = mm[0]*cdx*cdx + mm[3]*cdy*cdy + mm[5]*cdz*cdz
490 + 2.0*(mm[1]*cdx*cdy + mm[2]*cdx*cdz + mm[4]*cdy*cdz);
491
492 /* quality */
493 rap = h1 + h2 + h3 + h4 + h5 + h6;
494
495 num = sqrt(rap) * rap;
496
497 cal = det / num;
498
499 return cal;
500}
501
502#endif
MMG5_pMesh * mesh
Definition: API_functionsf_3d.c:65
int MMG5_moymet(MMG5_pMesh mesh, MMG5_pSol met, MMG5_pTetra pt, double *m1)
Definition: anisosiz_3d.c:70
int MMG5_cenrad_iso(MMG5_pMesh mesh, double *ct, double *c, double *rad)
Definition: cenrad_3d.c:45
inlined Functions
static double MMG5_lenSurfEdg_ani(MMG5_pMesh mesh, MMG5_pSol met, int np0, int np1, int8_t isedg)
Definition: inlined_functions.h:198
static double MMG5_lenSurfEdg33_ani(MMG5_pMesh mesh, MMG5_pSol met, int np0, int np1, int8_t isedg)
Definition: inlined_functions.h:266
static double MMG5_lenSurfEdg_iso(MMG5_pMesh mesh, MMG5_pSol met, int ip1, int ip2, int8_t isedg)
Definition: inlined_functions.h:291
static double MMG5_caltet_iso(MMG5_pMesh mesh, MMG5_pSol met, MMG5_pTetra pt)
Definition: inlined_functions_3d.h:387
static double MMG5_caltet_ani(MMG5_pMesh mesh, MMG5_pSol met, MMG5_pTetra pt)
Definition: inlined_functions_3d.h:417
static double MMG5_caltet_iso_4pt(double *a, double *b, double *c, double *d)
Definition: inlined_functions_3d.h:328
static double MMG5_lenedg_iso(MMG5_pMesh mesh, MMG5_pSol met, int ia, MMG5_pTetra pt)
Definition: inlined_functions_3d.h:230
static double MMG5_lenedg33_ani(MMG5_pMesh mesh, MMG5_pSol met, int ia, MMG5_pTetra pt)
Definition: inlined_functions_3d.h:94
static double MMG3D_caltetLES_iso(MMG5_pMesh mesh, MMG5_pSol met, MMG5_pTetra pt)
Definition: inlined_functions_3d.h:283
static double MMG5_orcal(MMG5_pMesh mesh, MMG5_pSol met, int iel)
Definition: inlined_functions_3d.h:262
static double MMG5_lenedgspl_iso(MMG5_pMesh mesh, MMG5_pSol met, int ia, MMG5_pTetra pt)
Definition: inlined_functions_3d.h:240
static double MMG5_lenedg_ani(MMG5_pMesh mesh, MMG5_pSol met, int ia, MMG5_pTetra pt)
Definition: inlined_functions_3d.h:200
static double MMG5_lenedgspl_ani(MMG5_pMesh mesh, MMG5_pSol met, int ia, MMG5_pTetra pt)
Definition: inlined_functions_3d.h:158
static double MMG5_lenedgCoor_ani(double *ca, double *cb, double *sa, double *sb)
Compute edge length from edge's coordinates.
Definition: inlined_functions_3d.h:56
static double MMG5_lenedgspl33_ani(MMG5_pMesh mesh, MMG5_pSol met, int ia, MMG5_pTetra pt)
Definition: inlined_functions_3d.h:125
#define MMG3D_ALPHAD
Definition: mmg3d.h:119
#define MMG3D_DET2VOL
Definition: mmg3d.h:128
static const uint8_t MMG5_iare[6][2]
vertices of extremities of the edges of the tetra
Definition: mmg3d.h:156
double(* MMG5_caltet)(MMG5_pMesh mesh, MMG5_pSol met, MMG5_pTetra pt)
Definition: mmg3dexterns.c:7
#define MG_GEO
Definition: mmgcommon.h:138
#define MG_MIN(a, b)
Definition: mmgcommon.h:133
#define MG_SIN(tag)
Definition: mmgcommon.h:161
#define MG_BDY
Definition: mmgcommon.h:141
#define MG_NOM
Definition: mmgcommon.h:140
double MMG5_orvol(MMG5_pPoint point, int *v)
Definition: tools.c:838
#define MMG5_EPSD2
Definition: mmgcommon.h:86
MMG mesh structure.
Definition: libmmgtypes.h:575
MMG5_pPoint point
Definition: libmmgtypes.h:612
MMG5_pTetra tetra
Definition: libmmgtypes.h:614
MMG5_pxTetra xtetra
Definition: libmmgtypes.h:615
Structure to store points of a MMG mesh.
Definition: libmmgtypes.h:252
int16_t tag
Definition: libmmgtypes.h:264
double c[3]
Definition: libmmgtypes.h:253
Definition: libmmgtypes.h:633
double * m
Definition: libmmgtypes.h:642
Definition: libmmgtypes.h:381
int v[4]
Definition: libmmgtypes.h:383
int xt
Definition: libmmgtypes.h:387
int16_t tag[6]
Definition: libmmgtypes.h:405