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 
55 static
56 inline 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 
93 static
94 inline double MMG5_lenedg33_ani(MMG5_pMesh mesh ,MMG5_pSol met, int ia,
95  MMG5_pTetra pt)
96 {
97  int ip1,ip2;
98  char 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 
124 static
125 inline double MMG5_lenedgspl33_ani(MMG5_pMesh mesh ,MMG5_pSol met, int ia,
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 
157 static
158 inline 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 
199 static
200 inline double MMG5_lenedg_ani(MMG5_pMesh mesh ,MMG5_pSol met, int ia,
201  MMG5_pTetra pt)
202 {
203  int ip1,ip2;
204  char 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 
229 static
230 inline 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 
239 static
240 inline 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 
261 static
262 inline 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 
282 static
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)/6.;
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 
327 static
328 inline 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 
386 static
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 
416 static
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_lenedgCoor_ani
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
MG_BDY
#define MG_BDY
Definition: mmgcommon.h:145
MMG5_Point::c
double c[3]
Definition: libmmgtypes.h:215
MMG5_moymet
int MMG5_moymet(MMG5_pMesh mesh, MMG5_pSol met, MMG5_pTetra pt, double *m1)
Definition: anisosiz_3d.c:69
MMG5_orvol
double MMG5_orvol(MMG5_pPoint point, int *v)
Definition: tools.c:836
MMG5_lenSurfEdg_iso
static double MMG5_lenSurfEdg_iso(MMG5_pMesh mesh, MMG5_pSol met, int ip1, int ip2, char isedg)
Definition: inlined_functions.h:291
MMG5_lenSurfEdg_ani
static double MMG5_lenSurfEdg_ani(MMG5_pMesh mesh, MMG5_pSol met, int np0, int np1, char isedg)
Definition: inlined_functions.h:198
MMG5_Tetra
Definition: libmmgtypes.h:339
MMG5_Sol
Definition: libmmgtypes.h:563
MG_SIN
#define MG_SIN(tag)
Definition: mmgcommon.h:163
MMG5_lenedg_ani
static double MMG5_lenedg_ani(MMG5_pMesh mesh, MMG5_pSol met, int ia, MMG5_pTetra pt)
Definition: inlined_functions_3d.h:200
MMG5_caltet_iso_4pt
static double MMG5_caltet_iso_4pt(double *a, double *b, double *c, double *d)
Definition: inlined_functions_3d.h:328
MG_NOM
#define MG_NOM
Definition: mmgcommon.h:144
MMG5_lenedg_iso
static double MMG5_lenedg_iso(MMG5_pMesh mesh, MMG5_pSol met, int ia, MMG5_pTetra pt)
Definition: inlined_functions_3d.h:230
MMG5_lenedgspl_iso
static double MMG5_lenedgspl_iso(MMG5_pMesh mesh, MMG5_pSol met, int ia, MMG5_pTetra pt)
Definition: inlined_functions_3d.h:240
MMG5_cenrad_iso
int MMG5_cenrad_iso(MMG5_pMesh mesh, double *ct, double *c, double *rad)
Definition: cenrad_3d.c:45
MMG5_xTetra::tag
int16_t tag[6]
Definition: libmmgtypes.h:363
MMG5_Mesh::point
MMG5_pPoint point
Definition: libmmgtypes.h:542
inlined_functions.h
inlined Functions
MMG3D_ALPHAD
#define MMG3D_ALPHAD
Definition: mmg3d.h:118
MMG5_Mesh::xtetra
MMG5_pxTetra xtetra
Definition: libmmgtypes.h:545
MMG5_orcal
static double MMG5_orcal(MMG5_pMesh mesh, MMG5_pSol met, int iel)
Definition: inlined_functions_3d.h:262
MMG5_lenedgspl_ani
static double MMG5_lenedgspl_ani(MMG5_pMesh mesh, MMG5_pSol met, int ia, MMG5_pTetra pt)
Definition: inlined_functions_3d.h:158
mesh
MMG5_pMesh * mesh
Definition: API_functionsf_3d.c:65
MG_GEO
#define MG_GEO
Definition: mmgcommon.h:142
MMG5_Point
Structure to store points of a MMG mesh.
Definition: libmmgtypes.h:214
MMG5_Sol::m
double * m
Definition: libmmgtypes.h:571
MMG5_caltet_ani
static double MMG5_caltet_ani(MMG5_pMesh mesh, MMG5_pSol met, MMG5_pTetra pt)
Definition: inlined_functions_3d.h:417
MG_MIN
#define MG_MIN(a, b)
Definition: mmgcommon.h:137
MMG5_Point::tag
int16_t tag
Definition: libmmgtypes.h:223
MMG5_lenedg33_ani
static double MMG5_lenedg33_ani(MMG5_pMesh mesh, MMG5_pSol met, int ia, MMG5_pTetra pt)
Definition: inlined_functions_3d.h:94
MMG5_Mesh::tetra
MMG5_pTetra tetra
Definition: libmmgtypes.h:544
MMG5_EPSD2
#define MMG5_EPSD2
Definition: mmgcommon.h:96
MMG5_lenedgspl33_ani
static double MMG5_lenedgspl33_ani(MMG5_pMesh mesh, MMG5_pSol met, int ia, MMG5_pTetra pt)
Definition: inlined_functions_3d.h:125
MMG5_lenSurfEdg33_ani
static double MMG5_lenSurfEdg33_ani(MMG5_pMesh mesh, MMG5_pSol met, int np0, int np1, char isedg)
Definition: inlined_functions.h:266
MMG5_Mesh
MMG mesh structure.
Definition: libmmgtypes.h:509
MMG3D_caltetLES_iso
static double MMG3D_caltetLES_iso(MMG5_pMesh mesh, MMG5_pSol met, MMG5_pTetra pt)
Definition: inlined_functions_3d.h:283
MMG5_iare
static const unsigned char MMG5_iare[6][2]
vertices of extremities of the edges of the tetra
Definition: mmg3d.h:153
mmg3d.h
MMG5_caltet_iso
static double MMG5_caltet_iso(MMG5_pMesh mesh, MMG5_pSol met, MMG5_pTetra pt)
Definition: inlined_functions_3d.h:387
MMG5_caltet
double(* MMG5_caltet)(MMG5_pMesh mesh, MMG5_pSol met, MMG5_pTetra pt)
Definition: mmg3dexterns.c:6
MMG5_Tetra::xt
int xt
Definition: libmmgtypes.h:345
MMG5_Tetra::v
int v[4]
Definition: libmmgtypes.h:341