mmgs
inlined_functions.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 "mmgcommon.h"
36 
37 #ifndef _INLINED_FUNC_H
38 #define _INLINED_FUNC_H
39 
53 static inline
54 double MMG5_lenEdg(MMG5_pMesh mesh,int np0,int np1,
55  double *m0,double *m1,int8_t isedg) {
56  MMG5_pPoint p0,p1;
57  double gammaprim0[3],gammaprim1[3],t[3],*n1,*n2,ux,uy,uz,ps1,ps2,l0,l1;
58  static int8_t mmgWarn=0;
59 
60  p0 = &mesh->point[np0];
61  p1 = &mesh->point[np1];
62 
63  ux = p1->c[0] - p0->c[0];
64  uy = p1->c[1] - p0->c[1];
65  uz = p1->c[2] - p0->c[2];
66 
67  /* computation of the two tangent vectors to the underlying curve of [i0i1] */
68  if ( MG_SIN(p0->tag) || (MG_NOM & p0->tag) ) {
69  gammaprim0[0] = ux;
70  gammaprim0[1] = uy;
71  gammaprim0[2] = uz;
72  }
73  else if ( isedg ) {
74  memcpy(t,p0->n,3*sizeof(double));
75  ps1 = ux*t[0] + uy*t[1] + uz*t[2];
76  gammaprim0[0] = ps1*t[0];
77  gammaprim0[1] = ps1*t[1];
78  gammaprim0[2] = ps1*t[2];
79  }
80  else {
81  if ( MG_GEO & p0->tag ) {
82  //assert(p0->xp);
83  n1 = &mesh->xpoint[p0->xp].n1[0];
84  n2 = &mesh->xpoint[p0->xp].n2[0];
85  ps1 = ux*n1[0] + uy*n1[1] + uz*n1[2];
86  ps2 = ux*n2[0] + uy*n2[1] + uz*n2[2];
87 
88  if ( fabs(ps2) < fabs(ps1) ) {
89  n1 = &mesh->xpoint[p0->xp].n2[0];
90  ps1 = ps2;
91  }
92  }
93  else if ( MG_REF & p0->tag || MG_BDY & p0->tag ) {
94  // ( MG_BDY & p0->tag ) => mmg3d
95  n1 = &mesh->xpoint[p0->xp].n1[0];
96  ps1 = ux*n1[0] + uy*n1[1] + uz*n1[2];
97  }
98  else {
99  // we come from mmgs because in mmg3d the boundary points are tagged
100  // MG_BDY.
101  n1 = &(p0->n[0]);
102  ps1 = ux*n1[0] + uy*n1[1] + uz*n1[2];
103  }
104  gammaprim0[0] = ux - ps1*n1[0];
105  gammaprim0[1] = uy - ps1*n1[1];
106  gammaprim0[2] = uz - ps1*n1[2];
107  }
108 
109  if ( MG_SIN(p1->tag) || (MG_NOM & p1->tag) ) {
110  gammaprim1[0] = -ux;
111  gammaprim1[1] = -uy;
112  gammaprim1[2] = -uz;
113  }
114  else if ( isedg ) {
115  memcpy(t,p1->n,3*sizeof(double));
116  ps1 = -ux*t[0] - uy*t[1] - uz*t[2];
117  gammaprim1[0] = ps1*t[0];
118  gammaprim1[1] = ps1*t[1];
119  gammaprim1[2] = ps1*t[2];
120  }
121  else {
122  if ( MG_GEO & p1->tag ) {
123  n1 = &mesh->xpoint[p1->xp].n1[0];
124  n2 = &mesh->xpoint[p1->xp].n2[0];
125  ps1 = -ux*n1[0] - uy*n1[1] - uz*n1[2];
126  ps2 = -ux*n2[0] - uy*n2[1] - uz*n2[2];
127 
128  if ( fabs(ps2) < fabs(ps1) ) {
129  n1 = &mesh->xpoint[p1->xp].n2[0];
130  ps1 = ps2;
131  }
132  }
133  else if ( MG_REF & p1->tag || MG_BDY & p1->tag ) {
134  // ( MG_BDY & p1->tag ) => mmg3d )
135  n1 = &mesh->xpoint[p1->xp].n1[0];
136  ps1 = - ux*n1[0] - uy*n1[1] - uz*n1[2];
137  }
138  else {
139  // we come from mmgs because in mmg3d the boundary points are tagged
140  // MG_BDY.
141  n1 = &(p1->n[0]);
142  ps1 = -ux*n1[0] - uy*n1[1] - uz*n1[2];
143  }
144  gammaprim1[0] = - ux - ps1*n1[0];
145  gammaprim1[1] = - uy - ps1*n1[1];
146  gammaprim1[2] = - uz - ps1*n1[2];
147  }
148 
149  /* computation of the length of the two tangent vectors in their respective
150  * tangent plane */
151  /* l_ab = int_a^b sqrt(m_ij d_t x_i(t) d_t x_j(t) ) : evaluated by a 2-point
152  * quadrature method. */
153  l0 = m0[0]*gammaprim0[0]*gammaprim0[0] + m0[3]*gammaprim0[1]*gammaprim0[1] \
154  + m0[5]*gammaprim0[2]*gammaprim0[2] \
155  + 2.0*m0[1]*gammaprim0[0]*gammaprim0[1] + 2.0*m0[2]*gammaprim0[0]*gammaprim0[2] \
156  + 2.0*m0[4]*gammaprim0[1]*gammaprim0[2];
157 
158  l1 = m1[0]*gammaprim1[0]*gammaprim1[0] + m1[3]*gammaprim1[1]*gammaprim1[1] \
159  + m1[5]*gammaprim1[2]*gammaprim1[2] \
160  +2.0*m1[1]*gammaprim1[0]*gammaprim1[1] + 2.0*m1[2]*gammaprim1[0]*gammaprim1[2] \
161  + 2.0*m1[4]*gammaprim1[1]*gammaprim1[2];
162 
163  if( l0 < 0.) {
164  if ( !mmgWarn ) {
165  fprintf(stderr," ## Warning: %s: at least 1 negative edge length "
166  "(%e)\n",__func__,l0);
167  mmgWarn = 1;
168  }
169  return 0.;
170  }
171  if(l1 < 0.) {
172  if ( !mmgWarn ) {
173  fprintf(stderr," ## Warning: %s: at least 1 negative edge length "
174  "(%e)\n",__func__,l1);
175  mmgWarn = 1;
176  }
177  return 0.;
178  }
179  l0 = 0.5*(sqrt(l0) + sqrt(l1));
180 
181  return l0;
182 }
183 
197 static inline
198 double MMG5_lenSurfEdg_ani(MMG5_pMesh mesh,MMG5_pSol met,int np0,int np1,int8_t isedg) {
199  MMG5_pPoint p0,p1;
200  double *m0,*m1,met0[6],met1[6],ux,uy,uz,rbasis[3][3];
201  static int8_t mmgWarn = 0;
202 
203  p0 = &mesh->point[np0];
204  p1 = &mesh->point[np1];
205 
206  ux = p1->c[0] - p0->c[0];
207  uy = p1->c[1] - p0->c[1];
208  uz = p1->c[2] - p0->c[2];
209 
210  /* Set metrics */
211  if ( MG_SIN(p0->tag) || (MG_NOM & p0->tag)) {
212  m0 = &met->m[6*np0];
213  }
214  else if ( MG_GEO & p0->tag ) {
215  /* Note that rbasis isn't used here */
216  if ( !MMG5_buildridmet(mesh,met,np0,ux,uy,uz,met0,rbasis) ) {
217  if ( !mmgWarn ) {
218  fprintf(stderr," ## Warning: %s: a- unable to compute at least 1 ridge"
219  " metric.\n",__func__);
220  mmgWarn = 1;
221  }
222  return 0.;
223  }
224  m0 = met0;
225  }
226  else {
227  m0 = &met->m[6*np0];
228  }
229 
230  if ( MG_SIN(p1->tag) || (MG_NOM & p1->tag)) {
231  m1 = &met->m[6*np1];
232  }
233  else if ( MG_GEO & p1->tag ) {
234  /* Note that rbasis isn't used here */
235  if ( !MMG5_buildridmet(mesh,met,np1,ux,uy,uz,met1,rbasis) ) {
236  if ( !mmgWarn ) {
237  fprintf(stderr," ## Warning: %s: b- unable to compute at least 1 ridge"
238  " metric.\n",__func__);
239  mmgWarn = 1;
240  }
241  return 0.;
242  }
243  m1 = met1;
244  }
245  else {
246  m1 = &met->m[6*np1];
247  }
248 
249  return MMG5_lenEdg(mesh,np0,np1,m0,m1,isedg);
250 }
251 
252 
265 static inline
267  int np0,int np1,int8_t isedg) {
268  double *m0,*m1;
269 
270  /* Set metrics */
271  m0 = &met->m[6*np0];
272  m1 = &met->m[6*np1];
273 
274  return MMG5_lenEdg(mesh,np0,np1,m0,m1,isedg);
275 }
276 
290 static
291 inline double MMG5_lenSurfEdg_iso(MMG5_pMesh mesh,MMG5_pSol met,int ip1,int ip2, int8_t isedg) {
292  MMG5_pPoint p1,p2;
293  double h1,h2,l,r,len;
294 
295  p1 = &mesh->point[ip1];
296  p2 = &mesh->point[ip2];
297  h1 = met->m[ip1];
298  h2 = met->m[ip2];
299  l = (p2->c[0]-p1->c[0])*(p2->c[0]-p1->c[0]) + (p2->c[1]-p1->c[1])*(p2->c[1]-p1->c[1]) \
300  + (p2->c[2]-p1->c[2])*(p2->c[2]-p1->c[2]);
301  l = sqrt(l);
302  r = h2 / h1 - 1.0;
303  len = fabs(r) < MMG5_EPS ? l / h1 : l / (h2-h1) * log1p(r);
304 
305  return len;
306 }
307 
308 #endif
MMG5_pMesh * mesh
Definition: API_functionsf_s.c:63
#define MMG5_EPS
Definition: eigenv.h:32
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_lenEdg(MMG5_pMesh mesh, int np0, int np1, double *m0, double *m1, int8_t isedg)
Definition: inlined_functions.h:54
int MMG5_buildridmet(MMG5_pMesh mesh, MMG5_pSol met, int np0, double ux, double uy, double uz, double mr[6], double r[3][3])
Definition: mettools.c:127
#define MG_GEO
Definition: mmgcommon.h:138
#define MG_SIN(tag)
Definition: mmgcommon.h:161
#define MG_BDY
Definition: mmgcommon.h:141
#define MG_NOM
Definition: mmgcommon.h:140
#define MG_REF
Definition: mmgcommon.h:137
MMG mesh structure.
Definition: libmmgtypes.h:575
MMG5_pPoint point
Definition: libmmgtypes.h:612
MMG5_pxPoint xpoint
Definition: libmmgtypes.h:613
Structure to store points of a MMG mesh.
Definition: libmmgtypes.h:252
double n[3]
Definition: libmmgtypes.h:254
int16_t tag
Definition: libmmgtypes.h:264
int xp
Definition: libmmgtypes.h:259
double c[3]
Definition: libmmgtypes.h:253
Definition: libmmgtypes.h:633
double * m
Definition: libmmgtypes.h:642
double n2[3]
Definition: libmmgtypes.h:275
double n1[3]
Definition: libmmgtypes.h:275