Mmg
Simplicial remeshers (mesh adaptation, isovalue discretization, lagrangian movement)
anisosiz_s.c
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
36#include "libmmgs_private.h"
37#include "libmmgs.h"
38#include "inlined_functions_private.h"
39#include "mmgsexterns_private.h"
40#include "mmgexterns_private.h"
41
54static int MMG5_defmetsin(MMG5_pMesh mesh,MMG5_pSol met,MMG5_int it,int ip) {
55 MMG5_pTria pt;
56 MMG5_pPoint p0;
57 MMG5_pPar par;
58 double *m,n[3],isqhmin,isqhmax,b0[3],b1[3],ps1,tau[3];
59 double ntau2,gammasec[3];
60 double c[3],kappa,maxkappa,alpha,hausd,hausd_v;
61 MMG5_int list[MMGS_LMAX+2],k,iel,idp,init_s;
62 int ilist,i;
63 uint8_t i0,i1,i2;
64
65 pt = &mesh->tria[it];
66 idp = pt->v[ip];
67 p0 = &mesh->point[idp];
68
69 /* local parameters at vertex: useless for now because new points are created
70 * without reference (inside the domain) */
71 hausd_v = mesh->info.hausd;
72 isqhmin = mesh->info.hmin;
73 isqhmax = mesh->info.hmax;
74
75 int8_t dummy;
76 ilist = MMG5_boulet(mesh,it,ip,list,1,&dummy);
77 if ( ilist < 1 ) return 0;
78
79 maxkappa = 0.0;
80 for (k=0; k<ilist; k++) {
81 iel = list[k] / 3;
82 i0 = list[k] % 3;
83 i1 = MMG5_inxt2[i0];
84 i2 = MMG5_iprv2[i0];
85 pt = &mesh->tria[iel];
86
87 /* Computation of the two control points associated to edge p0p1 with
88 * p1=mesh->point[pt->v[i1]]: p0 is singular */
89 MMG5_nortri(mesh,pt,n);
90 if ( MG_EDG(pt->tag[i2]) )
91 MMG5_bezierEdge(mesh,idp,pt->v[i1],b0,b1,1,n);
92 else
93 MMG5_bezierEdge(mesh,idp,pt->v[i1],b0,b1,0,n);
94
95 /* tangent vector */
96 tau[0] = 3.0*(b0[0] - p0->c[0]);
97 tau[1] = 3.0*(b0[1] - p0->c[1]);
98 tau[2] = 3.0*(b0[2] - p0->c[2]);
99 ntau2 = tau[0]*tau[0] + tau[1]*tau[1] + tau[2]*tau[2];
100
101 /* 2nd order derivative */
102 gammasec[0] = 6.0*p0->c[0] - 12.0*b0[0] + 6.0*b1[0];
103 gammasec[1] = 6.0*p0->c[1] - 12.0*b0[1] + 6.0*b1[1];
104 gammasec[2] = 6.0*p0->c[2] - 12.0*b0[2] + 6.0*b1[2];
105 if ( ntau2 < MMG5_EPSD ) continue;
106 ntau2 = 1.0 / ntau2;
107
108 /* derivative via the normal parametrization */
109 ps1 = gammasec[0]*tau[0] + gammasec[1]*tau[1] + gammasec[2]*tau[2];
110 c[0] = gammasec[0] - ps1*tau[0]*ntau2;
111 c[1] = gammasec[1] - ps1*tau[1]*ntau2;
112 c[2] = gammasec[2] - ps1*tau[2]*ntau2;
113
114 kappa = ntau2 * sqrt(c[0]*c[0] + c[1]*c[1] + c[2]*c[2]);
115
116 /* local parameters for triangle */
117 hausd = hausd_v;
118 init_s = 0;
119 for (i=0; i<mesh->info.npar; i++) {
120 par = &mesh->info.par[i];
121 if ( (par->elt == MMG5_Triangle) && (pt->ref == par->ref ) ) {
122 hausd = par->hausd;
123 if ( !init_s ) {
124 isqhmin = par->hmin;
125 isqhmax = par->hmax;
126 init_s = 1;
127 }
128 else {
129 isqhmin = MG_MAX(par->hmin,isqhmin);
130 isqhmax = MG_MIN(par->hmax,isqhmax);
131 }
132 }
133 }
134 maxkappa = MG_MAX(kappa/hausd,maxkappa);
135 }
136
137 isqhmin = 1.0 / (isqhmin*isqhmin);
138 isqhmax = 1.0 / (isqhmax*isqhmax);
139
140 alpha = 1.0 / 8.0 * maxkappa;
141 alpha = MG_MIN(alpha,isqhmin);
142 alpha = MG_MAX(alpha,isqhmax);
143
144 m = &met->m[6*idp];
145 memset(m,0,6*sizeof(double));
146 m[0] = m[3] = m[5] = alpha;
147
148 return 1;
149}
150
173static int MMG5_defmetrid(MMG5_pMesh mesh,MMG5_pSol met,MMG5_int it,int ip) {
174 MMG5_pTria pt;
175 MMG5_pPoint p0,p1,p2;
176 MMG5_Bezier b;
177 MMG5_pPar par;
178 MMG5_int k,iel,idp,*list,list1[MMGS_LMAX+2];
179 int ilist1,ilist2,ilist;
180 MMG5_int list2[MMGS_LMAX+2],iprid[2],isloc;
181 int ier;
182 double *m,isqhmin,isqhmax,*n1,*n2,*n,*t,trot[2],u[2];
183 double r[3][3],lispoi[3*MMGS_LMAX+1],ux,uy,uz,det,bcu[3];
184 double detg,detd;
185 uint8_t i,i0,i1,i2;
186 static int8_t mmgWarn0=0;
187
188 pt = &mesh->tria[it];
189 idp = pt->v[ip];
190 p0 = &mesh->point[idp];
191
192 /* local parameters */
193 isqhmin = mesh->info.hmin;
194 isqhmax = mesh->info.hmax;
195 isloc = 0;
196 for (k=0; k<mesh->info.npar; k++) {
197 par = &mesh->info.par[k];
198 if ( (par->elt == MMG5_Triangle) && (pt->ref == par->ref ) ) {
199 if ( !isloc ) {
200 isqhmin = par->hmin;
201 isqhmax = par->hmax;
202 isloc = 1;
203 }
204 else {
205 isqhmin = MG_MAX(isqhmin,par->hmin);
206 isqhmax = MG_MIN(isqhmax,par->hmax);
207 }
208 }
209 }
210
211 isqhmin = 1.0 / (isqhmin*isqhmin);
212 isqhmax = 1.0 / (isqhmax*isqhmax);
213
214 n1 = &mesh->xpoint[p0->xp].n1[0];
215 n2 = &mesh->xpoint[p0->xp].n2[0];
216 t = p0->n;
217
218 m = &met->m[6*idp];
219 memset(m,0,6*sizeof(double));
220 m[0] = isqhmax;
221 m[1] = isqhmax;
222 m[2] = isqhmax;
223 m[3] = isqhmax;
224 m[4] = isqhmax;
225
226 ier = bouletrid(mesh,it,ip,&ilist1,list1,&ilist2,list2,&iprid[0],&iprid[1]);
227 if ( !ier ) {
228 if ( !mmgWarn0 ) {
229 mmgWarn0 = 1;
230 fprintf(stderr,"\n ## Error: %s: unable to compute the two balls at at"
231 " least 1 ridge point.\n",__func__);
232 }
233 return 0;
234 }
235
236 /* Specific size in direction of t */
237 m[0] = MG_MAX(m[0],MMG5_ridSizeInTangentDir(mesh,p0,idp,iprid,isqhmin,isqhmax));
238
239 /* Characteristic sizes in directions u1 and u2 */
240 for (i=0; i<2; i++) {
241 if ( i==0 ) {
242 n = n1;
243 ilist = ilist1;
244 list = &list1[0];
245 }
246 else {
247 n = n2;
248 ilist = ilist2;
249 list = &(list2[0]);
250 }
251 MMG5_rotmatrix(n,r);
252
253 /* Apply rotation to the half-ball under consideration */
254 i1 = 0;
255 for (k=0; k<ilist; k++) {
256 iel = list[k] / 3;
257 i0 = list[k] % 3;
258 i1 = MMG5_inxt2[i0];
259 pt = &mesh->tria[iel];
260 p1 = &mesh->point[pt->v[i1]];
261
262 ux = p1->c[0] - p0->c[0];
263 uy = p1->c[1] - p0->c[1];
264 uz = p1->c[2] - p0->c[2];
265
266 lispoi[3*k+1] = r[0][0]*ux + r[0][1]*uy + r[0][2]*uz;
267 lispoi[3*k+2] = r[1][0]*ux + r[1][1]*uy + r[1][2]*uz;
268 lispoi[3*k+3] = r[2][0]*ux + r[2][1]*uy + r[2][2]*uz;
269 }
270
271 /* last point : the half-ball is open : ilist tria, and ilist +1 points ;
272 lists are enumerated in direct order */
273 i2 = MMG5_inxt2[i1];
274 p2 = &mesh->point[pt->v[i2]];
275
276 ux = p2->c[0] - p0->c[0];
277 uy = p2->c[1] - p0->c[1];
278 uz = p2->c[2] - p0->c[2];
279
280 lispoi[3*ilist+1] = r[0][0]*ux + r[0][1]*uy + r[0][2]*uz;
281 lispoi[3*ilist+2] = r[1][0]*ux + r[1][1]*uy + r[1][2]*uz;
282 lispoi[3*ilist+3] = r[2][0]*ux + r[2][1]*uy + r[2][2]*uz;
283
284 /* At this point, lispoi contains all the points of the half-ball of p0, rotated
285 so that t_{p_0}S = [z = 0] */
286
287 /* Rotated tangent vector (trot[2] = 0), and third direction */
288 trot[0] = r[0][0]*t[0] + r[0][1]*t[1] + r[0][2]*t[2];
289 trot[1] = r[1][0]*t[0] + r[1][1]*t[1] + r[1][2]*t[2];
290
291 u[0] = -trot[1];
292 u[1] = trot[0];
293
294 /* Travel half-ball at p0 and stop at first triangle containing u */
295 for (k=0; k<ilist; k++) {
296 detg = lispoi[3*k+1]*u[1] - lispoi[3*k+2]*u[0];
297 detd = u[0]*lispoi[3*(k+1)+2] - u[1]*lispoi[3*(k+1)+1];
298 if ( detg > 0.0 && detd > 0.0 ) break;
299 }
300
301 /* If triangle not found, try with -u */
302 if ( k == ilist ) {
303 u[0] *= -1.0;
304 u[1] *= -1.0;
305
306 for (k=0; k<ilist; k++) {
307 detg = lispoi[3*k+1]*u[1] - lispoi[3*k+2]*u[0];
308 detd = u[0]*lispoi[3*(k+1)+2] - u[1]*lispoi[3*(k+1)+1];
309 if ( detg > 0.0 && detd > 0.0 ) break;
310 }
311 }
312 if ( k == ilist ) continue;
313
314 iel = list[k] / 3;
315 i0 = list[k] % 3;
316 pt = &mesh->tria[iel];
317 if ( !MMG5_bezierCP(mesh,pt,&b,1) ) continue;
318
319 /* Barycentric coordinates of vector u in tria iel */
320 detg = lispoi[3*k+1]*u[1] - lispoi[3*k+2]*u[0];
321 detd = u[0]*lispoi[3*(k+1)+2] - u[1]*lispoi[3*(k+1)+1];
322 det = detg + detd;
323 if ( det < MMG5_EPSD ) continue;
324
325 det = 1.0 / det;
326 bcu[0] = 0.0;
327 bcu[1] = u[0]*lispoi[3*(k+1)+2] - u[1]*lispoi[3*(k+1)+1];
328 bcu[1] *= det;
329 assert(bcu[1] <= 1.0);
330 bcu[2] = 1.0 - bcu[1];
331
332 /* Computation of tangent vector and second derivative of curve t \mapsto
333 b(tbcu) (not in rotated frame) */
334 m[i+1] = MG_MAX(m[i+1],
335 MMG5_ridSizeInNormalDir(mesh,i0,bcu,&b,isqhmin,isqhmax));
336 }
337
338 return 1;
339}
340
352static int MMG5_defmetref(MMG5_pMesh mesh,MMG5_pSol met,MMG5_int it,int ip) {
353 MMG5_pTria pt;
354 MMG5_pPoint p0,p1;
355 MMG5_Bezier b;
356 MMG5_pPar par;
357 MMG5_int list[MMGS_LMAX+2],k,iel,ipref[2],idp,isloc;
358 int i,ilist;
359 double *m,isqhmin,isqhmax,*n,r[3][3],lispoi[3*MMGS_LMAX+1];
360 double ux,uy,uz,det2d,intm[3],c[3];
361 double tAA[6],tAb[3],hausd;
362 uint8_t i0,i1,i2;
363 static int8_t mmgWarn0=0;
364
365 ipref[0] = ipref[1] = 0;
366 pt = &mesh->tria[it];
367 idp = pt->v[ip];
368 p0 = &mesh->point[idp];
369
370 int8_t dummy;
371 ilist = MMG5_boulet(mesh,it,ip,list,1,&dummy);
372 if ( ilist < 1 )
373 return 0;
374
375 /* Computation of the rotation matrix T_p0 S -> [z = 0] */
376 n = &mesh->xpoint[p0->xp].n1[0];
377 assert ( n[0]*n[0] + n[1]*n[1] + n[2]*n[2] > MMG5_EPSD2 );
378
379 MMG5_rotmatrix(n,r);
380 m = &met->m[6*idp];
381
382 /* Apply rotation \circ translation to the whole ball */
383 assert ( ilist );
384 for (k=0; k<ilist; k++) {
385 iel = list[k] / 3;
386 i0 = list[k] % 3;
387 i1 = MMG5_inxt2[i0];
388 i2 = MMG5_iprv2[i0];
389 pt = &mesh->tria[iel];
390 p1 = &mesh->point[pt->v[i1]];
391
392 /* Store the two ending points of ref curves */
393 if ( MG_REF & pt->tag[i1] ) {
394 if ( !ipref[0] ) {
395 ipref[0] = pt->v[i2];
396 }
397 else if ( !ipref[1] && (pt->v[i2] != ipref[0]) ) {
398 ipref[1] = pt->v[i2];
399 }
400 else if ( (pt->v[i2] != ipref[0]) && (pt->v[i2] != ipref[1]) ) {
401 if ( !mmgWarn0 ) {
402 mmgWarn0 = 1;
403 fprintf(stderr,"\n ## Warning: %s: at least 1 metric not computed:"
404 " non singular point at intersection of 3 ref edges.\n",
405 __func__);
406 }
407 return 0;
408 }
409 }
410
411 if ( MG_REF & pt->tag[i2] ) {
412 if ( !ipref[0] ) {
413 ipref[0] = pt->v[i1];
414 }
415 else if ( !ipref[1] && (pt->v[i1] != ipref[0]) ) {
416 ipref[1] = pt->v[i1];
417 }
418 else if ( (pt->v[i1] != ipref[0]) && (pt->v[i1] != ipref[1]) ) {
419 if ( !mmgWarn0 ) {
420 mmgWarn0 = 1;
421 fprintf(stderr,"\n ## Warning: %s: at least 1 metric not computed:"
422 " non singular point at intersection of 3 ref edges.\n",
423 __func__);
424 }
425 return 0;
426 }
427 }
428
429 ux = p1->c[0] - p0->c[0];
430 uy = p1->c[1] - p0->c[1];
431 uz = p1->c[2] - p0->c[2];
432
433 lispoi[3*k+1] = r[0][0]*ux + r[0][1]*uy + r[0][2]*uz;
434 lispoi[3*k+2] = r[1][0]*ux + r[1][1]*uy + r[1][2]*uz;
435 lispoi[3*k+3] = r[2][0]*ux + r[2][1]*uy + r[2][2]*uz;
436 }
437
438 /* list goes modulo ilist */
439 lispoi[3*ilist+1] = lispoi[1];
440 lispoi[3*ilist+2] = lispoi[2];
441 lispoi[3*ilist+3] = lispoi[3];
442
443 /* Check all projections over tangent plane. */
444 for (k=0; k<ilist-1; k++) {
445 det2d = lispoi[3*k+1]*lispoi[3*(k+1)+2] - lispoi[3*k+2]*lispoi[3*(k+1)+1];
446 if ( det2d < 0.0 ) {
447 return 0;
448 }
449 }
450 det2d = lispoi[3*(ilist-1)+1]*lispoi[3*0+2] - lispoi[3*(ilist-1)+2]*lispoi[3*0+1];
451 if ( det2d < 0.0 ) {
452 return 0;
453 }
454 assert(ipref[0] && ipref[1]);
455
456 /* At this point, lispoi contains all the points of the ball of p0, rotated
457 so that t_{p_0}S = [z = 0], ipref1 and ipref2 are the indices of other ref points. */
458
459 /* Second step : reconstitution of the curvature tensor at p0 in the tangent plane,
460 with a quadric fitting approach */
461 memset(intm,0x00,3*sizeof(double));
462 memset(tAA,0x00,6*sizeof(double));
463 memset(tAb,0x00,3*sizeof(double));
464
465 hausd = mesh->info.hausd;
466 isqhmin = mesh->info.hmin;
467 isqhmax = mesh->info.hmax;
468 isloc = 0;
469 for (k=0; k<ilist; k++) {
470 /* Approximation of the curvature in the normal section associated to tau : by assumption,
471 p1 is either regular, either on a ridge (or a singularity), but p0p1 is not ridge*/
472 iel = list[k] / 3;
473 i0 = list[k] % 3;
474 pt = &mesh->tria[iel];
475 MMG5_bezierCP(mesh,pt,&b,1);
476
477
478 /* 1. Fill matrice tAA and second member tAb with A=(\sum X_{P_i}^2 \sum
479 * Y_{P_i}^2 \sum X_{P_i}Y_{P_i}) and b=\sum Z_{P_i} with P_i the physical
480 * points at edge [i0;i1] extremities and middle.
481 * 2. Compute the physical coor \a c of the curve edge's
482 * mid-point.
483 */
484 MMG5_fillDefmetregSys(k,p0,i0,b,r,c,lispoi,tAA,tAb);
485
486 /* local parameters */
487 for (i=0; i<mesh->info.npar; i++) {
488 par = &mesh->info.par[i];
489 if ( (par->elt == MMG5_Triangle) && (pt->ref == par->ref ) ) {
490 if ( !isloc ) {
491 hausd = par->hausd;
492 isqhmin = par->hmin;
493 isqhmax = par->hmax;
494 isloc = 1;
495 }
496 else {
497 // take the wanted value between the two local parameters asked
498 // by the user.
499 hausd = MG_MIN(hausd,par->hausd);
500 isqhmin = MG_MAX(isqhmin,par->hmin);
501 isqhmax = MG_MIN(isqhmax,par->hmax);
502 }
503 }
504 }
505 }
506
507 isqhmin = 1.0 / (isqhmin*isqhmin);
508 isqhmax = 1.0 / (isqhmax*isqhmax);
509
510 return MMG5_solveDefmetrefSys(mesh,p0,ipref,r,c,tAA,tAb,m,
511 isqhmin,isqhmax,hausd);
512}
513
529int MMGS_surfballRotation(MMG5_pMesh mesh,MMG5_pPoint p0,MMG5_int *list,int ilist,
530 double r[3][3],double *lispoi,double n[3]) {
531 MMG5_pTria pt;
532 MMG5_pPoint p1;
533 double ux,uy,uz,area;
534 MMG5_int iel;
535 int i0,i1,k;
536
537 /* Computation of the rotation matrix T_p0 S -> [z = 0] */
538 assert ( n[0]*n[0] + n[1]*n[1] + n[2]*n[2] > MMG5_EPSD2 );
539
540 if ( !MMG5_rotmatrix(n,r) ) {
541 return 0;
542 }
543
544 /* Apply rotation \circ translation to the whole ball */
545 assert ( ilist );
546 for (k=0; k<ilist; k++) {
547 iel = list[k] / 3;
548 i0 = list[k] % 3;
549 i1 = MMG5_inxt2[i0];
550 pt = &mesh->tria[iel];
551 p1 = &mesh->point[pt->v[i1]];
552
553 ux = p1->c[0] - p0->c[0];
554 uy = p1->c[1] - p0->c[1];
555 uz = p1->c[2] - p0->c[2];
556
557 lispoi[3*k+1] = r[0][0]*ux + r[0][1]*uy + r[0][2]*uz;
558 lispoi[3*k+2] = r[1][0]*ux + r[1][1]*uy + r[1][2]*uz;
559 lispoi[3*k+3] = r[2][0]*ux + r[2][1]*uy + r[2][2]*uz;
560 }
561
562 /* list goes modulo ilist */
563 lispoi[3*ilist+1] = lispoi[1];
564 lispoi[3*ilist+2] = lispoi[2];
565 lispoi[3*ilist+3] = lispoi[3];
566
567 /* Check all projections over tangent plane. */
568 for (k=0; k<ilist-1; k++) {
569 area = lispoi[3*k+1]*lispoi[3*(k+1)+2] - lispoi[3*k+2]*lispoi[3*(k+1)+1];
570 if ( area <= 0.0 ) {
571 return 0;
572 }
573 }
574 area = lispoi[3*(ilist-1)+1]*lispoi[3*0+2] - lispoi[3*(ilist-1)+2]*lispoi[3*0+1];
575 if ( area <= 0.0 ) {
576 return 0;
577 }
578
579 return 1;
580}
581
593static int MMG5_defmetreg(MMG5_pMesh mesh,MMG5_pSol met,MMG5_int it,int ip) {
594 MMG5_pTria pt;
595 MMG5_pPoint p0;
596 MMG5_Bezier b;
597 MMG5_pPar par;
598 MMG5_int list[MMGS_LMAX+2],iel,idp,isloc;
599 int ilist,k,i;
600 double *m,r[3][3],lispoi[3*MMGS_LMAX+1];
601 double c[3],isqhmin,isqhmax;
602 double tAA[6],tAb[3],hausd;
603 uint8_t i0;
604
605 pt = &mesh->tria[it];
606 idp = pt->v[ip];
607 p0 = &mesh->point[idp];
608 m = &met->m[6*idp];
609
610 int8_t dummy;
611 ilist = MMG5_boulet(mesh,it,ip,list,1,&dummy);
612 if ( ilist < 1 )
613 return 0;
614
615 /* Rotation of the ball of p0 */
616 if ( !MMGS_surfballRotation(mesh,p0,list,ilist,r,lispoi,p0->n) ) {
617 return 0;
618 }
619
620 /* At this point, lispoi contains all the points of the ball of p0, rotated
621 so that t_{p_0}S = [z = 0] */
622
623 /* Second step : reconstitution of the curvature tensor at p0 in the tangent
624 plane, with a quadric fitting approach */
625 memset(tAA,0x00,6*sizeof(double));
626 memset(tAb,0x00,3*sizeof(double));
627
628 hausd = mesh->info.hausd;
629 isqhmin = mesh->info.hmin;
630 isqhmax = mesh->info.hmax;
631 isloc = 0;
632 for (k=0; k<ilist; k++) {
633 /* Approximation of the curvature in the normal section associated to tau :
634 by assumption, p1 is either regular, either on a ridge (or a
635 singularity), but p0p1 is not ridge*/
636 iel = list[k] / 3;
637 i0 = list[k] % 3;
638 pt = &mesh->tria[iel];
639 MMG5_bezierCP(mesh,pt,&b,1);
640
641 /* 1. Fill matrice tAA and second member tAb with A=(\sum X_{P_i}^2 \sum
642 * Y_{P_i}^2 \sum X_{P_i}Y_{P_i}) and b=\sum Z_{P_i} with P_i the physical
643 * points at edge [i0;i1] extremities and middle.
644 * 2. Compute the physical coor \a c of the curve edge's
645 * mid-point.
646 */
647 MMG5_fillDefmetregSys(k,p0,i0,b,r,c,lispoi,tAA,tAb);
648
649 /* local parameters */
650 for (i=0; i<mesh->info.npar; i++) {
651 par = &mesh->info.par[i];
652 if ( (par->elt == MMG5_Triangle) && (pt->ref == par->ref ) ) {
653 if ( !isloc ) {
654 hausd = par->hausd;
655 isqhmin = par->hmin;
656 isqhmax = par->hmax;
657 isloc = 1;
658 }
659 else {
660 // take the minimum value between the two local parameters asked by
661 // the user.
662 hausd = MG_MIN(hausd,par->hausd);
663 isqhmin = MG_MAX(isqhmin,par->hmin);
664 isqhmax = MG_MIN(isqhmax,par->hmax);
665 }
666 }
667 }
668 }
669
670 isqhmin = 1.0 / (isqhmin*isqhmin);
671 isqhmax = 1.0 / (isqhmax*isqhmax);
672
673 /* 2. Solve tAA * tmp_m = tAb and fill m with tmp_m (after rotation) */
674 return(MMG5_solveDefmetregSys( mesh,r, c, tAA, tAb, m, isqhmin, isqhmax,
675 hausd));
676}
677
691static inline
692int MMGS_intextmet(MMG5_pMesh mesh,MMG5_pSol met,MMG5_int np,double me[6]) {
693 MMG5_pPoint p0;
694 double *n;
695 double dummy_n[3];
696
697 p0 = &mesh->point[np];
698
699 dummy_n[0] = dummy_n[1] = dummy_n[2] = 0.;
700
701 if ( MG_SIN(p0->tag) || (p0->tag & MG_NOM) ) {
702 n = &dummy_n[0];
703 }
704 else {
705 n = &p0->n[0];
706 }
707
708 return MMG5_mmgIntextmet(mesh,met,np,me,n);
709
710}
711
735int MMGS_defsiz_ani(MMG5_pMesh mesh,MMG5_pSol met) {
736 MMG5_pTria pt;
737 MMG5_pPoint ppt;
738 double mm[6];
739 MMG5_int k;
740 int8_t ismet;
741 int8_t i;
742 static int8_t mmgErr=0;
743
744 if ( !MMG5_defsiz_startingMessage (mesh,met,__func__) ) {
745 return 0;
746 }
747
748 for (k=1; k<=mesh->np; k++) {
749 ppt = &mesh->point[k];
750 ppt->flag = 0;
751 ppt->s = 0;
752 }
753
754 if ( met->m ) {
755 assert ( met->np );
756 ismet = 1;
757 }
758 else {
759 ismet = 0;
760
761 MMG5_calelt = MMG5_caltri_ani;
762 MMG5_lenSurfEdg = MMG5_lenSurfEdg_ani;
763
764 if ( !MMGS_Set_solSize(mesh,met,MMG5_Vertex,mesh->np,3) )
765 return 0;
766 }
767
771 if ( !mesh->info.nosizreq ) {
772 if ( !MMGS_set_metricAtPointsOnReqEdges ( mesh,met,ismet ) ) {
773 return 0;
774 }
775 }
776
778 for (k=1; k<=mesh->nt; k++) {
779 pt = &mesh->tria[k];
780 if ( !MG_EOK(pt) || pt->ref < 0 ) continue;
781
782 for (i=0; i<3; i++) {
783 ppt = &mesh->point[pt->v[i]];
784 if ( ppt->flag || !MG_VOK(ppt) ) continue;
785 if ( ismet ) memcpy(mm,&met->m[6*(pt->v[i])],6*sizeof(double));
786
787 if ( MS_SIN(ppt->tag) ) {
788 if ( !MMG5_defmetsin(mesh,met,k,i) ) continue;
789 }
790 else if ( ppt->tag & MG_GEO ) {
791 if ( !MMG5_defmetrid(mesh,met,k,i)) continue;
792 }
793 else if ( ppt->tag & MG_REF ) {
794 if ( !MMG5_defmetref(mesh,met,k,i) ) continue;
795 }
796 else if ( ppt->tag ) continue;
797 else {
798 if ( !MMG5_defmetreg(mesh,met,k,i) ) continue;
799 }
800 if ( ismet ) {
801 if ( !MMGS_intextmet(mesh,met,pt->v[i],mm) ) {
802 if ( !mmgErr ) {
803 fprintf(stderr,"\n ## Error: %s: unable to intersect metrics"
804 " at point %" MMG5_PRId ".\n",__func__,
805 MMGS_indPt(mesh,pt->v[i]));
806 mmgErr = 1;
807 }
808 return 0;
809 }
810 }
811 ppt->flag = 1;
812 }
813 }
814 /* Now the metric storage at ridges is the "mmg" one. */
815 mesh->info.metRidTyp = 1;
816
821 MMG5_defUninitSize ( mesh,met,ismet );
822
823 return 1;
824}
825
835int MMGS_gradsiz_ani(MMG5_pMesh mesh,MMG5_pSol met) {
836 MMG5_pPoint p1;
837 double *m,mv;
838 MMG5_int k;
839 int it;
840
841 if ( abs(mesh->info.imprim) > 5 || mesh->info.ddebug )
842 fprintf(stdout," ** Anisotropic mesh gradation\n");
843
844 /* First step : make ridges iso in each apairing direction */
845 // remark ALGIANE: a mettre à plat : on veut vraiment faire ça ? Pour une demi
846 // sphère, on n'a pas envie de garder une métrique "infinie du coté du plan et
847 // petite dans la direction de la sphère?? "
848 for (k=1; k<= mesh->np; k++) {
849 p1 = &mesh->point[k];
850 if ( !MG_VOK(p1) ) continue;
851 if ( MS_SIN(p1->tag) ) continue;
852 if ( !(p1->tag & MG_GEO) ) continue;
853
854 m = &met->m[6*k];
855 mv = MG_MAX(m[1],m[2]);
856 m[1] = mv;
857 m[2] = mv;
858 mv = MG_MAX(m[3],m[4]);
859 m[3] = mv;
860 m[4] = mv;
861 }
862
863 /* Second step : standard gradation procedure */
864 MMG5_gradsiz_ani(mesh,met,&it);
865
866 return 1;
867}
MMGS_Set_solSize
int MMGS_Set_solSize(MMG5_pMesh mesh, MMG5_pSol sol, int typEntity, MMG5_int np, int typSol)
Initialize an array of solution fields: set dimension, types and number of fields.
Definition: API_functions_s.c:93
ier
int ier
Definition: API_functionsf_2d.c:761
mesh
MMG5_pMesh * mesh
Definition: API_functionsf_2d.c:66
MMG5_fillDefmetregSys
void MMG5_fillDefmetregSys(MMG5_int k, MMG5_pPoint p0, int i0, MMG5_Bezier b, double r[3][3], double c[3], double *lispoi, double tAA[6], double tAb[3])
Definition: anisosiz.c:290
MMG5_solveDefmetregSys
int MMG5_solveDefmetregSys(MMG5_pMesh mesh, double r[3][3], double c[3], double tAA[6], double tAb[3], double *m, double isqhmin, double isqhmax, double hausd)
Definition: anisosiz.c:436
MMG5_defUninitSize
void MMG5_defUninitSize(MMG5_pMesh mesh, MMG5_pSol met, int8_t ismet)
Definition: anisosiz.c:228
MMG5_ridSizeInNormalDir
double MMG5_ridSizeInNormalDir(MMG5_pMesh mesh, int i0, double *bcu, MMG5_Bezier *b, double isqhmin, double isqhmax)
Definition: anisosiz.c:836
MMG5_solveDefmetrefSys
int MMG5_solveDefmetrefSys(MMG5_pMesh mesh, MMG5_pPoint p0, MMG5_int ipref[2], double r[3][3], double c[3], double tAA[6], double tAb[3], double *m, double isqhmin, double isqhmax, double hausd)
Definition: anisosiz.c:538
MMG5_gradsiz_ani
MMG5_int MMG5_gradsiz_ani(MMG5_pMesh mesh, MMG5_pSol met, int *it)
Definition: anisosiz.c:2259
MMG5_ridSizeInTangentDir
double MMG5_ridSizeInTangentDir(MMG5_pMesh mesh, MMG5_pPoint p0, MMG5_int idp, MMG5_int *iprid, double isqhmin, double isqhmax)
Definition: anisosiz.c:764
MMG5_defmetreg
static int MMG5_defmetreg(MMG5_pMesh mesh, MMG5_pSol met, MMG5_int it, int ip)
Definition: anisosiz_s.c:593
MMGS_gradsiz_ani
int MMGS_gradsiz_ani(MMG5_pMesh mesh, MMG5_pSol met)
Definition: anisosiz_s.c:835
MMGS_intextmet
static int MMGS_intextmet(MMG5_pMesh mesh, MMG5_pSol met, MMG5_int np, double me[6])
Definition: anisosiz_s.c:692
MMGS_defsiz_ani
int MMGS_defsiz_ani(MMG5_pMesh mesh, MMG5_pSol met)
Definition: anisosiz_s.c:735
MMG5_defmetrid
static int MMG5_defmetrid(MMG5_pMesh mesh, MMG5_pSol met, MMG5_int it, int ip)
Definition: anisosiz_s.c:173
MMGS_surfballRotation
int MMGS_surfballRotation(MMG5_pMesh mesh, MMG5_pPoint p0, MMG5_int *list, int ilist, double r[3][3], double *lispoi, double n[3])
Definition: anisosiz_s.c:529
MMG5_defmetsin
static int MMG5_defmetsin(MMG5_pMesh mesh, MMG5_pSol met, MMG5_int it, int ip)
Definition: anisosiz_s.c:54
MMG5_defmetref
static int MMG5_defmetref(MMG5_pMesh mesh, MMG5_pSol met, MMG5_int it, int ip)
Definition: anisosiz_s.c:352
MMG5_bezierEdge
void MMG5_bezierEdge(MMG5_pMesh mesh, MMG5_int i0, MMG5_int i1, double b0[3], double b1[3], int8_t isrid, double v[3])
Definition: bezier.c:51
MMG5_boulet
int MMG5_boulet(MMG5_pMesh mesh, MMG5_int start, int ip, MMG5_int *list, int8_t s, int8_t *opn)
Definition: boulep.c:363
bouletrid
int bouletrid(MMG5_pMesh mesh, MMG5_int start, MMG5_int ip, int *il1, MMG5_int *l1, int *il2, MMG5_int *l2, MMG5_int *ip0, MMG5_int *ip1)
Definition: boulep_s.c:201
MMG5_EPSD
#define MMG5_EPSD
Definition: eigenv_private.h:31
MMGS_indPt
MMG5_int MMGS_indPt(MMG5_pMesh mesh, MMG5_int kp)
Definition: gentools_s.c:139
inlined_functions_private.h
MMG5_lenSurfEdg_ani
static double MMG5_lenSurfEdg_ani(MMG5_pMesh mesh, MMG5_pSol met, MMG5_int np0, MMG5_int np1, int8_t isedg)
Definition: inlined_functions_private.h:199
MMG5_defsiz_startingMessage
int MMG5_defsiz_startingMessage(MMG5_pMesh mesh, MMG5_pSol met, const char *funcname)
Definition: isosiz.c:77
MMGS_set_metricAtPointsOnReqEdges
int MMGS_set_metricAtPointsOnReqEdges(MMG5_pMesh mesh, MMG5_pSol met, int8_t ismet)
Definition: isosiz_s.c:90
libmmgs.h
API headers and documentation for the mmgs library.
MMGS_LMAX
#define MMGS_LMAX
Definition: libmmgs.h:98
libmmgs_private.h
MS_SIN
#define MS_SIN(tag)
Definition: libmmgs_private.h:49
MMG5_Vertex
@ MMG5_Vertex
Definition: libmmgtypes.h:230
MMG5_Triangle
@ MMG5_Triangle
Definition: libmmgtypes.h:232
MMG5_mmgIntextmet
int MMG5_mmgIntextmet(MMG5_pMesh mesh, MMG5_pSol met, MMG5_int np, double me[6], double n[3])
Definition: mettools.c:1139
MG_GEO
#define MG_GEO
Definition: mmgcommon_private.h:146
MMG5_caltri_ani
double MMG5_caltri_ani(MMG5_pMesh mesh, MMG5_pSol met, MMG5_pTria ptt)
Definition: quality.c:115
MG_EOK
#define MG_EOK(pt)
Definition: mmgcommon_private.h:167
MG_MIN
#define MG_MIN(a, b)
Definition: mmgcommon_private.h:141
MG_MAX
#define MG_MAX(a, b)
Definition: mmgcommon_private.h:140
MG_EDG
#define MG_EDG(tag)
Definition: mmgcommon_private.h:173
MMG5_iprv2
static const uint8_t MMG5_iprv2[3]
Definition: mmgcommon_private.h:586
MG_SIN
#define MG_SIN(tag)
Definition: mmgcommon_private.h:169
MMG5_inxt2
static const uint8_t MMG5_inxt2[6]
Definition: mmgcommon_private.h:585
MG_VOK
#define MG_VOK(ppt)
Definition: mmgcommon_private.h:166
MMG5_nortri
int MMG5_nortri(MMG5_pMesh mesh, MMG5_pTria pt, double *n)
Definition: tools.c:209
MMG5_rotmatrix
int MMG5_rotmatrix(double n[3], double r[3][3])
Definition: tools.c:467
MG_NOM
#define MG_NOM
Definition: mmgcommon_private.h:148
MMG5_EPSD2
#define MMG5_EPSD2
Definition: mmgcommon_private.h:95
MG_REF
#define MG_REF
Definition: mmgcommon_private.h:145
mmgexterns_private.h
mmgsexterns_private.h
MMG5_Bezier
Definition: mmgcommon_private.h:611
MMG5_Info::ddebug
int8_t ddebug
Definition: libmmgtypes.h:539
MMG5_Info::hmin
double hmin
Definition: libmmgtypes.h:525
MMG5_Info::npar
int npar
Definition: libmmgtypes.h:530
MMG5_Info::metRidTyp
uint8_t metRidTyp
Definition: libmmgtypes.h:554
MMG5_Info::hmax
double hmax
Definition: libmmgtypes.h:525
MMG5_Info::nosizreq
uint8_t nosizreq
Definition: libmmgtypes.h:553
MMG5_Info::par
MMG5_pPar par
Definition: libmmgtypes.h:524
MMG5_Info::hausd
double hausd
Definition: libmmgtypes.h:525
MMG5_Info::imprim
int imprim
Definition: libmmgtypes.h:536
MMG5_Mesh
MMG mesh structure.
Definition: libmmgtypes.h:613
MMG5_Mesh::info
MMG5_Info info
Definition: libmmgtypes.h:659
MMG5_Mesh::point
MMG5_pPoint point
Definition: libmmgtypes.h:649
MMG5_Mesh::xpoint
MMG5_pxPoint xpoint
Definition: libmmgtypes.h:650
MMG5_Mesh::nt
MMG5_int nt
Definition: libmmgtypes.h:620
MMG5_Mesh::tria
MMG5_pTria tria
Definition: libmmgtypes.h:655
MMG5_Mesh::np
MMG5_int np
Definition: libmmgtypes.h:620
MMG5_Par
Local parameters for a specific entity and reference.
Definition: libmmgtypes.h:263
MMG5_Par::hmin
double hmin
Definition: libmmgtypes.h:264
MMG5_Par::hmax
double hmax
Definition: libmmgtypes.h:265
MMG5_Par::hausd
double hausd
Definition: libmmgtypes.h:266
MMG5_Par::ref
MMG5_int ref
Definition: libmmgtypes.h:267
MMG5_Par::elt
int8_t elt
Definition: libmmgtypes.h:268
MMG5_Point
Structure to store vertices of an MMG mesh.
Definition: libmmgtypes.h:276
MMG5_Point::n
double n[3]
Definition: libmmgtypes.h:278
MMG5_Point::c
double c[3]
Definition: libmmgtypes.h:277
MMG5_Point::tag
uint16_t tag
Definition: libmmgtypes.h:290
MMG5_Point::xp
MMG5_int xp
Definition: libmmgtypes.h:285
MMG5_Point::s
MMG5_int s
Definition: libmmgtypes.h:289
MMG5_Point::flag
MMG5_int flag
Definition: libmmgtypes.h:288
MMG5_Sol
Definition: libmmgtypes.h:671
MMG5_Sol::m
double * m
Definition: libmmgtypes.h:680
MMG5_Sol::np
MMG5_int np
Definition: libmmgtypes.h:674
MMG5_Tria
Structure to store triangles of a MMG mesh.
Definition: libmmgtypes.h:338
MMG5_Tria::ref
MMG5_int ref
Definition: libmmgtypes.h:341
MMG5_Tria::tag
uint16_t tag[3]
Definition: libmmgtypes.h:348
MMG5_Tria::v
MMG5_int v[3]
Definition: libmmgtypes.h:340
MMG5_xPoint::n2
double n2[3]
Definition: libmmgtypes.h:301
MMG5_xPoint::n1
double n1[3]
Definition: libmmgtypes.h:301
Generated on Tue Oct 29 2024 11:49:58 for Mmg by doxygen 1.9.6