mmg2d
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
53static inline
54double 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
197static inline
198double 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
265static 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
290static
291inline 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_2d.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