Mmg
Simplicial remeshers (mesh adaptation, isovalue discretization, lagrangian movement)
movpt_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 <math.h>
37
38#include "libmmgs_private.h"
39#include "mmgexterns_private.h"
40
52int movintpt_iso(MMG5_pMesh mesh,MMG5_pSol met,MMG5_int *list,int ilist) {
53 MMG5_pPoint p0,p1,ppt0;
54 MMG5_pTria pt,pt0;
55 MMG5_Bezier b;
56 double aa,bb,ab,ll,l,mlon,devmean,GV[3],gv[2],cosalpha,sinalpha,r[3][3],*n,lispoi[3*MMG5_TRIA_LMAX+3];
57 double ux,uy,uz,det2d,detloc,step,lambda[3],uv[2],o[3],no[3],to[3],Vold,Vnew,calold,calnew,caltmp;
58 MMG5_int k,iel,ipp,ibeg,iend;
59 int ier,kel,npt;
60 int8_t i0,i1,i2;
61 static int8_t mmgErr0=0,mmgErr1=0;
62
63 step = 0.1;
64 Vold = 0.0;
65 Vnew = 0.0;
66
67 k = list[0] / 3;
68 i0 = list[0] % 3;
69 i1 = MMG5_inxt2[i0];
70 pt = &mesh->tria[k];
71 ibeg = pt->v[i1];
72 ipp = pt->v[i0];
73 p0 = &mesh->point[ipp]; /* point to move */
74
75 k = list[ilist-1] / 3;
76 i0 = list[ilist-1] % 3;
77 i1 = MMG5_inxt2[i0];
78 i2 = MMG5_inxt2[i1];
79 pt = &mesh->tria[k];
80 iend = pt->v[i2];
81
82 /* check for open ball */
83 if ( iend != ibeg ) return 0;
84
85 npt = ilist; // number of POINTS in the ball = number of triangles. Each point is counted as the
86 // i1 of its associated triangle
87
88 /* Step 1 : computation of the gradient of variance as a function from R^3 to R
89 compute mlon := 1/n \sum_{i=1,...,n}{||p-p_i||^2} */
90 mlon = 0.0;
91 for (k=0; k<ilist; k++) {
92 iel = list[k] / 3;
93 i0 = list[k] % 3;
94 i1 = MMG5_inxt2[i0];
95 pt = &mesh->tria[iel];
96 p1 = &mesh->point[pt->v[i1]];
97 mlon += (p1->c[0]-p0->c[0])*(p1->c[0]-p0->c[0]) + (p1->c[1]-p0->c[1])*(p1->c[1]-p0->c[1]) \
98 + (p1->c[2]-p0->c[2])*(p1->c[2]-p0->c[2]);
99 }
100 mlon = mlon / npt;
101
102 /* compute GV := 4/n * sum_{i=1,...,n}{(||p-p_i||^2 -m)(p-pi)} and
103 Vold = sum_{i=1,...,n}{(||p-p_i||^2 -m)^2} */
104 GV[0] = GV[1] = GV[2] = 0.0;
105 for (k=0; k<ilist; k++) {
106 iel = list[k] / 3;
107 i0 = list[k] % 3;
108 i1 = MMG5_inxt2[i0];
109 pt = &mesh->tria[iel];
110 p1 = &mesh->point[pt->v[i1]];
111
112 devmean = (p1->c[0]-p0->c[0])*(p1->c[0]-p0->c[0]) + (p1->c[1]-p0->c[1])*(p1->c[1]-p0->c[1]) \
113 + (p1->c[2]-p0->c[2])*(p1->c[2]-p0->c[2]) - mlon;
114 GV[0] += devmean*(p1->c[0]-p0->c[0]);
115 GV[1] += devmean*(p1->c[1]-p0->c[1]);
116 GV[2] += devmean*(p1->c[2]-p0->c[2]);
117 Vold += devmean*devmean;
118 }
119
120 GV[0] *= (4.0 / npt);
121 GV[1] *= (4.0 / npt);
122 GV[2] *= (4.0 / npt);
123 /* Vold *= (1.0 / npt); */
124
125 /* Step 2 : computation of the rotation matrix T_p0 S -> [z = 0] */
126 n = &p0->n[0]; //once again, that depends on what kind of point is considered.
127 aa = n[0]*n[0];
128 bb = n[1]*n[1];
129 ab = n[0]*n[1];
130 ll = aa+bb;
131 cosalpha = n[2];
132 sinalpha = sqrt(1.0- MG_MIN(1.0,cosalpha*cosalpha));
133
134 /* No rotation needed in this case */
135 if ( ll < MMG5_EPS ) {
136 if ( n[2] > 0.0 ) {
137 r[0][0] = 1.0 ; r[0][1] = 0.0 ; r[0][2] = 0.0;
138 r[1][0] = 0.0 ; r[1][1] = 1.0 ; r[1][2] = 0.0;
139 r[2][0] = 0.0 ; r[2][1] = 0.0 ; r[2][2] = 1.0;
140 }
141 else {
142 r[0][0] = -1.0 ; r[0][1] = 0.0 ; r[0][2] = 0.0;
143 r[1][0] = 0.0 ; r[1][1] = 1.0 ; r[1][2] = 0.0;
144 r[2][0] = 0.0 ; r[2][1] = 0.0 ; r[2][2] = -1.0;
145 }
146 }
147 else {
148 l = sqrt(ll);
149
150 r[0][0] = (aa*cosalpha + bb)/ll;
151 r[0][1] = ab*(cosalpha-1)/ll;
152 r[0][2] = -n[0]*sinalpha/l;
153 r[1][0] = r[0][1];
154 r[1][1] = (bb*cosalpha + aa)/ll;
155 r[1][2] = -n[1]*sinalpha/l;
156 r[2][0] = n[0]*sinalpha/l;
157 r[2][1] = n[1]*sinalpha/l;
158 r[2][2] = cosalpha;
159 }
160
161 /* apply rotation to vector GV */
162 gv[0] = r[0][0]*GV[0] + r[0][1]*GV[1] + r[0][2]*GV[2];
163 gv[1] = r[1][0]*GV[0] + r[1][1]*GV[1] + r[1][2]*GV[2];
164
165 /* Apply rotation \circ translation to the whole ball */
166 assert ( ilist > 0 && ilist < MMG5_TRIA_LMAX );
167 for (k=0; k<ilist; k++) {
168 iel = list[k] / 3;
169 i0 = list[k] % 3;
170 i1 = MMG5_inxt2[i0];
171 pt = &mesh->tria[iel];
172 p1 = &mesh->point[pt->v[i1]];
173
174 ux = p1->c[0] - p0->c[0];
175 uy = p1->c[1] - p0->c[1];
176 uz = p1->c[2] - p0->c[2];
177
178 lispoi[3*k+1] = r[0][0]*ux + r[0][1]*uy + r[0][2]*uz;
179 lispoi[3*k+2] = r[1][0]*ux + r[1][1]*uy + r[1][2]*uz;
180 lispoi[3*k+3] = r[2][0]*ux + r[2][1]*uy + r[2][2]*uz;
181 }
182
183 /* list goes modulo ilist */
184 lispoi[3*ilist+1] = lispoi[1];
185 lispoi[3*ilist+2] = lispoi[2];
186 lispoi[3*ilist+3] = lispoi[3];
187
188 gv[0] = 0.0;
189 gv[1] = 0.0;
190 for (k=0; k<ilist; k++) {
191 gv[0] += lispoi[3*k+1];
192 gv[1] += lispoi[3*k+2];
193 }
194 gv[0] *= (1.0 / npt);
195 gv[1] *= (1.0 / npt);
196
197 /* Check all projections over tangent plane. */
198 for (k=0; k<ilist-1; k++) {
199 det2d = lispoi[3*k+1]*lispoi[3*(k+1)+2] - lispoi[3*k+2]*lispoi[3*(k+1)+1];
200 if ( det2d < 0.0 ) return 0;
201 }
202 det2d = lispoi[3*(ilist-1)+1]*lispoi[3*0+2] - lispoi[3*(ilist-1)+2]*lispoi[3*0+1];
203 if ( det2d < 0.0 ) return 0;
204
205 /* Step 3 : locate new point in the ball, and compute its barycentric coordinates */
206 det2d = lispoi[1]*gv[1] - lispoi[2]*gv[0];
207 kel = 0;
208 if ( det2d >= 0.0 ) {
209 for (k=0; k<ilist; k++) {
210 detloc = gv[0]*lispoi[3*(k+1)+2] - gv[1]*lispoi[3*(k+1)+1]; //orientation with respect to p2
211 if ( detloc >= 0.0 ) {
212 kel = k;
213 break;
214 }
215 }
216 if ( k == ilist ) return 0;
217 }
218 else {
219 for (k=ilist-1; k>=0; k--) {
220 detloc = gv[1]*lispoi[3*k+1] - gv[0]*lispoi[3*k+2]; //orientation with respect to p2
221 if ( detloc >= 0.0 ) {
222 kel = k;
223 break;
224 }
225 }
226 if ( k == -1 ) return 0;
227 }
228
229 /* Sizing of time step : make sure point does not go out corresponding triangle. */
230 det2d = -gv[1]*(lispoi[3*(kel+1)+1] - lispoi[3*(kel)+1]) \
231 + gv[0]*(lispoi[3*(kel+1)+2] - lispoi[3*(kel)+2]);
232 if ( fabs(det2d) < MMG5_EPSD ) return 0;
233
234 det2d = 1/det2d;
235 step *= det2d;
236 det2d = lispoi[3*(kel)+1]*(lispoi[3*(kel+1)+2] - lispoi[3*(kel)+2]) - \
237 lispoi[3*(kel)+2 ]*(lispoi[3*(kel+1)+1] - lispoi[3*(kel)+1]);
238 step *= det2d;
239 step = fabs(step);
240 gv[0] *= step;
241 gv[1] *= step;
242
243 /* Computation of the barycentric coordinates of the new point in the corresponding triangle. */
244 det2d = lispoi[3*kel+1]*lispoi[3*(kel+1)+2] - lispoi[3*kel+2]*lispoi[3*(kel+1)+1];
245 if ( det2d < MMG5_EPSD ) return 0;
246 det2d = 1.0 / det2d;
247 lambda[1] = lispoi[3*(kel+1)+2]*gv[0] - lispoi[3*(kel+1)+1]*gv[1];
248 lambda[2] = -lispoi[3*(kel)+2]*gv[0] + lispoi[3*(kel)+1]*gv[1];
249 lambda[1]*= (det2d);
250 lambda[2]*= (det2d);
251 lambda[0] = 1.0 - lambda[1] - lambda[2];
252
253 /* Step 4 : come back to original problem, and compute patch in triangle iel */
254 iel = list[kel]/3;
255 i0 = list[kel]%3;
256 pt = &mesh->tria[iel];
257
258 ier = MMG5_bezierCP(mesh,pt,&b,1);
259 if ( !ier ) {
260 if( !mmgErr0 ) {
261 mmgErr0 = 1;
262 fprintf(stderr,"\n ## Warning: %s: function MMG5_bezierCP return 0.\n",
263 __func__);
264 }
265 return 0;
266 }
267
268 /* Now, for Bezier interpolation, one should identify which of i,i1,i2 is 0,1,2
269 recall uv[0] = barycentric coord associated to pt->v[1], uv[1] associated to pt->v[2] */
270 if ( i0 == 0 ) {
271 uv[0] = lambda[1];
272 uv[1] = lambda[2];
273 }
274 else if ( i0 == 1 ) {
275 uv[0] = lambda[0];
276 uv[1] = lambda[1];
277 }
278 else {
279 uv[0] = lambda[2];
280 uv[1] = lambda[0];
281 }
282
283 ier = MMGS_bezierInt(&b,uv,o,no,to);
284 if ( !ier ) {
285 if( !mmgErr1 ) {
286 mmgErr1 = 1;
287 fprintf(stderr," ## Warning: %s: function MMGS_bezierInt return 0.\n",
288 __func__);
289 }
290 return 0;
291 }
292
293 /* First test : check whether variance has been decreased */
294 mlon = 0.0;
295 for (k=0; k<ilist; k++) {
296 iel = list[k] / 3;
297 i0 = list[k] % 3;
298 i1 = MMG5_inxt2[i0];
299 pt = &mesh->tria[iel];
300 p1 = &mesh->point[pt->v[i1]];
301 mlon += (p1->c[0]-o[0])*(p1->c[0]-o[0]) + (p1->c[1]-o[1])*(p1->c[1]-o[1]) \
302 + (p1->c[2]-o[2])*(p1->c[2]-o[2]);
303 }
304 mlon = mlon / npt;
305
306 for (k=0; k<ilist; k++) {
307 iel = list[k] / 3;
308 i0 = list[k] % 3;
309 i1 = MMG5_inxt2[i0];
310 pt = &mesh->tria[iel];
311 p1 = &mesh->point[pt->v[i1]];
312 devmean = (p1->c[0]-o[0])*(p1->c[0]-o[0]) + (p1->c[1]-o[1])*(p1->c[1]-o[1]) \
313 + (p1->c[2]-o[2])*(p1->c[2]-o[2]) - mlon;
314 Vnew += devmean*devmean;
315 }
316 /* Vnew *= (1.0 / npt); */
317 /* if ( Vold < Vnew ) return 0; */
318
319 /* Second test : check whether geometric approximation has not been too much degraded */
320 ppt0 = &mesh->point[0];
321 ppt0->tag = p0->tag;
322 ppt0->c[0] = o[0];
323 ppt0->c[1] = o[1];
324 ppt0->c[2] = o[2];
325
326 ppt0->n[0] = no[0];
327 ppt0->n[1] = no[1];
328 ppt0->n[2] = no[2];
329
330 calold = calnew = DBL_MAX;
331 for (k= 0; k<ilist; k++) {
332 iel = list[k] / 3;
333 i0 = list[k] % 3;
334 pt = &mesh->tria[iel];
335 pt0 = &mesh->tria[0];
336 memcpy(pt0,pt,sizeof(MMG5_Tria));
337 pt0->v[i0] = 0;
338 caltmp = caleltsig_iso(mesh,NULL,iel);
339 calold = MG_MIN(calold,caltmp);
340 caltmp = caleltsig_iso(mesh,NULL,0);
341 if ( caltmp < MMG5_NULKAL ) return 0;
342 calnew = MG_MIN(calnew,caltmp);
343 }
344 if ( calold < MMG5_EPSOK && calnew <= calold ) return 0;
345 else if (calnew < MMG5_EPSOK) return 0;
346 else if ( calnew < 0.3*calold ) return 0;
347
348 /* Finally, update coordinates and normals of point, if new position is accepted : */
349 p0->c[0] = o[0];
350 p0->c[1] = o[1];
351 p0->c[2] = o[2];
352
353 p0->n[0] = no[0];
354 p0->n[1] = no[1];
355 p0->n[2] = no[2];
356
357 return 1;
358}
359
375int MMGS_paramDisp(MMG5_pMesh mesh,MMG5_int it,int8_t isrid,MMG5_int ip0,MMG5_int ip,
376 double step,double o[3]) {
377 MMG5_pTria pt;
378 MMG5_Bezier b;
379 double uv[2],nn1[3],to[3];
380 int ier;
381
382 /* move towards ip */
383 pt = &mesh->tria[it];
384
385 ier = MMG5_bezierCP(mesh,pt,&b,1);
386 assert(ier);
387
388 /* fill table uv */
389 if ( pt->v[0] == ip0 ) {
390 if ( pt->v[1] == ip ) {
391 uv[0] = step;
392 uv[1] = 0.0;
393 }
394 else if ( pt->v[2] == ip ) {
395 uv[0] = 0.0;
396 uv[1] = step;
397 }
398 }
399 else if ( pt->v[0] == ip ) {
400 if ( pt->v[1] == ip0 ) {
401 uv[0] = 1.0 - step;
402 uv[1] = 0.0;
403 }
404 else if ( pt->v[2] == ip0 ) {
405 uv[0] = 0.0;
406 uv[1] = 1.0-step;
407 }
408 }
409 else {
410 if ( pt->v[1] == ip0 ) {
411 uv[0] = 1.0 - step;
412 uv[1] = step;
413 }
414 else if ( pt->v[2] == ip0 ) {
415 uv[0] = step;
416 uv[1] = 1.0-step;
417 }
418 }
419 ier = MMGS_bezierInt(&b,uv,o,nn1,to);
420 assert(ier);
421
422 return ier;
423}
424
448static
449int MMGS_update_normalAndTangent(MMG5_pMesh mesh,MMG5_pPoint p0,MMG5_pPoint p,
450 double llold,double lam0,double lam1,double lam2,
451 double no1[3],double no2[3],
452 double np1[3],double np2[3],
453 double nn1[3],double nn2[3],double to[3] ) {
454
455 double dd1,dd2,ddt,ps2;
456
457 nn1[0] = no1[0]+np1[0];
458 nn1[1] = no1[1]+np1[1];
459 nn1[2] = no1[2]+np1[2];
460
461 nn2[0] = no2[0]+np2[0];
462 nn2[1] = no2[1]+np2[1];
463 nn2[2] = no2[2]+np2[2];
464
465 ps2 = (p->c[0]-p0->c[0])*nn1[0]+(p->c[1]-p0->c[1])*nn1[1]+(p->c[2]-p0->c[2])*nn1[2];
466 if ( llold < MMG5_EPSD2 ) {
467 return 0;
468 }
469
470 ps2 *= (2.0 / llold);
471 nn1[0] -= ps2*(p->c[0]-p0->c[0]);
472 nn1[1] -= ps2*(p->c[1]-p0->c[1]);
473 nn1[2] -= ps2*(p->c[2]-p0->c[2]);
474
475 ps2 = (p->c[0]-p0->c[0])*nn2[0]+(p->c[1]-p0->c[1])*nn2[1]+(p->c[2]-p0->c[2])*nn2[2];
476 ps2 *= (2.0/llold);
477 nn2[0] -= ps2*(p->c[0]-p0->c[0]);
478 nn2[1] -= ps2*(p->c[1]-p0->c[1]);
479 nn2[2] -= ps2*(p->c[2]-p0->c[2]);
480
481 dd1 = nn1[0]*nn1[0] + nn1[1]*nn1[1] + nn1[2]*nn1[2];
482 dd2 = nn2[0]*nn2[0] + nn2[1]*nn2[1] + nn2[2]*nn2[2];
483 if ( (dd1 < MMG5_EPSD2) || (dd2<MMG5_EPSD2) ) {
484 return 0;
485 }
486
487 dd1 = 1.0 / sqrt(dd1);
488 nn1[0] = dd1*nn1[0];
489 nn1[1] = dd1*nn1[1];
490 nn1[2] = dd1*nn1[2];
491
492 dd2 = 1.0 / sqrt(dd2);
493 nn2[0] = dd2*nn2[0];
494 nn2[1] = dd2*nn2[1];
495 nn2[2] = dd2*nn2[2];
496
497 /* At this point, the control points associated to the interpolation formula for normals
498 have been computed .*/
499 nn1[0] = lam0*no1[0] + lam1*nn1[0] + lam2*np1[0];
500 nn1[1] = lam0*no1[1] + lam1*nn1[1] + lam2*np1[1];
501 nn1[2] = lam0*no1[2] + lam1*nn1[2] + lam2*np1[2];
502
503 nn2[0] = lam0*no2[0] + lam1*nn2[0] + lam2*np2[0];
504 nn2[1] = lam0*no2[1] + lam1*nn2[1] + lam2*np2[1];
505 nn2[2] = lam0*no2[2] + lam1*nn2[2] + lam2*np2[2];
506
507 to[0] = nn1[1]*nn2[2]-nn1[2]*nn2[1];
508 to[1] = nn1[2]*nn2[0]-nn1[0]*nn2[2];
509 to[2] = nn1[0]*nn2[1]-nn1[1]*nn2[0];
510
511 ddt = to[0]*to[0] + to[1]*to[1] + to[2]*to[2];
512 dd1 = nn1[0]*nn1[0] + nn1[1]*nn1[1] + nn1[2]*nn1[2];
513 dd2 = nn2[0]*nn2[0] + nn2[1]*nn2[1] + nn2[2]*nn2[2];
514
515 if ( (dd1 < MMG5_EPSD2) || (dd2<MMG5_EPSD2) || (ddt < MMG5_EPSD2) ) {
516 return 0;
517 }
518
519 dd1 = 1.0 / sqrt(dd1);
520 nn1[0] = dd1*nn1[0];
521 nn1[1] = dd1*nn1[1];
522 nn1[2] = dd1*nn1[2];
523
524 dd2 = 1.0 / sqrt(dd2);
525 nn2[0] = dd2*nn2[0];
526 nn2[1] = dd2*nn2[1];
527 nn2[2] = dd2*nn2[2];
528
529 ddt = 1.0 / sqrt(ddt);
530 to[0] = ddt*to[0];
531 to[1] = ddt*to[1];
532 to[2] = ddt*to[2];
533
534 return 1;
535}
536
555int MMGS_moveTowardPoint(MMG5_pMesh mesh,MMG5_pPoint p0,MMG5_pPoint p,
556 double llold,double lam0,double lam1,double lam2,
557 double nn1[3],double nn2[3],double to[3]) {
558
559 double *no1,*no2,*np1,*np2,psn11,psn12;
560
561 no1 = &mesh->xpoint[p0->xp].n1[0];
562 no2 = &mesh->xpoint[p0->xp].n2[0];
563
564 if ( MS_SIN(p->tag) ) {
565 np1 = &mesh->xpoint[p0->xp].n1[0];
566 np2 = &mesh->xpoint[p0->xp].n2[0];
567 }
568 else {
569 np1 = &mesh->xpoint[p->xp].n1[0];
570 np2 = &mesh->xpoint[p->xp].n2[0];
571 }
572 psn11 = no1[0]*np1[0] + no1[1]*np1[1] + no1[2]*np1[2];
573 psn12 = no1[0]*np2[0] + no1[1]*np2[1] + no1[2]*np2[2];
574
575 /* no1 goes with np1, no2 with np2 */
576 if ( fabs(1.0-fabs(psn11)) < fabs(1.0-fabs(psn12)) ){
577 if ( !MMGS_update_normalAndTangent( mesh,p0,p,llold,lam0,lam1,lam2,
578 no1,no2,np1,np2,nn1,nn2,to ) ) {
579 return 0;
580 }
581 }
582
583 /* no1 goes with np2 and no2 with np1 */
584 else {
585 if ( !MMGS_update_normalAndTangent( mesh,p0,p,llold,lam0,lam1,lam2,
586 no1,no2,np2,np1,nn1,nn2,to ) ) {
587 return 0;
588 }
589 }
590
591 return 1;
592}
593
594/* compute movement of a ridge point whose ball (consisting of triangles) is passed */
595int movridpt_iso(MMG5_pMesh mesh,MMG5_pSol met,MMG5_int *list,int ilist) {
596 MMG5_pTria pt,pt0;
597 MMG5_pxPoint go;
598 MMG5_pPoint p0,p1,p2,ppt0;
599 double step,ll1old,ll1new,ll2old,ll2new,o[3];
600 double nn1[3],nn2[3],to[3],calold,calnew,lam0,lam1,lam2;
601 MMG5_int k,iel,ip,ip0,ip1,ip2,it,it1,it2;
602 int8_t i0,i1,i2,isrid1,isrid2,isrid;
603
604 step = 0.1;
605 isrid1 = 0 ; isrid2 = 0;
606 it1 = 0 ; it2 = 0;
607 ip1 = 0 ; ip2 = 0;
608
609 /* First step : make sure 2 ridge edges pass through point, and recoved them, with neighbouring triangles */
610 for (k=0; k<ilist; k++) {
611 iel = list[k] / 3;
612 i0 = list[k] % 3;
613 i1 = MMG5_inxt2[i0];
614 i2 = MMG5_inxt2[i1];
615 pt = &mesh->tria[iel];
616
617 if ( MG_EDG(pt->tag[i1]) ) {
618 if ( !it1 ) {
619 it1 = iel;
620 ip1 = pt->v[i2]; // edge(i1) = (p0p2)
621 if ( pt->tag[i1] & MG_GEO ) isrid1 = 1;
622 }
623 else if ( it1 && !it2 ) {
624 if ( ip1 != pt->v[i2] ) {
625 it2 = iel;
626 ip2 = pt->v[i2]; // edge(i1) = (p0p2)
627 if ( pt->tag[i1] & MG_GEO ) isrid2 = 1;
628 }
629 }
630 else if ( it1 && it2 && (pt->v[i2] != ip1) && (pt->v[i2] != ip2) ) {
631 return 0;
632 }
633 }
634
635 if ( MG_EDG(pt->tag[i2]) ) {
636 if ( !it1 ) {
637 it1 = iel;
638 ip1 = pt->v[i1]; // edge(i2) = (p0p1)
639 if ( pt->tag[i2] & MG_GEO ) isrid1 = 1;
640 }
641 else if ( it1 && !it2 ) {
642 if ( ip1 != pt->v[i1] ) {
643 it2 = iel;
644 ip2 = pt->v[i1]; // edge(i1) = (p0p2)
645 if ( pt->tag[i2] & MG_GEO ) isrid2 = 1;
646 }
647 }
648 else if ( it1 && it2 && (pt->v[i1] != ip1) && (pt->v[i1] != ip2) ) {
649 return 0;
650 }
651 }
652 }
653
654 /* Second step : compute variance of lengths of edges (p0p1),(p0p2) and its gradient */
655 iel = list[0] / 3;
656 i0 = list[0] % 3;
657 pt = &mesh->tria[iel];
658 ip0 = pt->v[i0];
659 p0 = &mesh->point[ip0];
660 p1 = &mesh->point[ip1];
661 p2 = &mesh->point[ip2];
662
663 ll1old = (p1->c[0]-p0->c[0])*(p1->c[0]-p0->c[0])
664 + (p1->c[1]-p0->c[1])*(p1->c[1]-p0->c[1])
665 + (p1->c[2]-p0->c[2])*(p1->c[2]-p0->c[2]);
666
667 ll2old = (p2->c[0] -p0->c[0])*(p2->c[0]-p0->c[0])
668 + (p2->c[1]-p0->c[1])*(p2->c[1]-p0->c[1])
669 + (p2->c[2]-p0->c[2])*(p2->c[2]-p0->c[2]);
670
671 if ( (!ll1old) || (!ll2old) ) return 0;
672
673 if ( ll1old < ll2old ) { //move towards p2
674 ip = ip2;
675 isrid = isrid2;
676 it = it2;
677 }
678 else {
679 ip = ip1;
680 isrid = isrid1;
681 it = it1;
682 }
683
684
685 /* Third step : infer arc length of displacement, parameterized over edges */
686 if ( !MMGS_paramDisp ( mesh,it,isrid,ip0,ip,step,o) ) {
687 return 0;
688 }
689
690 /* check whether proposed move is admissible */
691 ll1new = (p1->c[0]-o[0])*(p1->c[0]-o[0])
692 + (p1->c[1]-o[1])* (p1->c[1]-o[1])
693 + (p1->c[2]-o[2])* (p1->c[2]-o[2]);
694
695 ll2new = (p2->c[0]-o[0])*(p2->c[0]-o[0])
696 + (p2->c[1]-o[1])*(p2->c[1]-o[1])
697 + (p2->c[2]-o[2])*(p2->c[2]-o[2]);
698
699 if( fabs(ll2new -ll1new) >= fabs(ll2old -ll1old) ) return 0;
700
701 /* normal and tangent updates */
702 // Bezier basis function of order 2
703 lam0 = (1.0-step)*(1.0-step);
704 lam1 = 2.0*step*(1.0-step);
705 lam2 = step*step;
706
707 /* move toward p2 */
708 if ( ll2old > ll1old ) {
709 if ( !MMGS_moveTowardPoint(mesh,p0,p2,ll2old,lam0,lam1,lam2,nn1,nn2,to) ) {
710 return 0;
711 }
712 }
713 else {
714 /* move toward p1 */
715 if ( !MMGS_moveTowardPoint( mesh,p0,p1,ll1old,lam0,lam1,lam2,nn1,nn2,to) ) {
716 return 0;
717 }
718 }
719
720 /* check proposed motion */
721 ppt0 = &mesh->point[0];
722 ppt0->xp = 0;
723 ppt0->tag = p0->tag;
724
725 go = &mesh->xpoint[0];
726 go->n1[0] = nn1[0];
727 go->n1[1] = nn1[1];
728 go->n1[2] = nn1[2];
729
730 if ( isrid ) {
731 go->n2[0] = nn2[0];
732 go->n2[1] = nn2[1];
733 go->n2[2] = nn2[2];
734 }
735 ppt0->c[0] = o[0];
736 ppt0->c[1] = o[1];
737 ppt0->c[2] = o[2];
738
739 ppt0->n[0] = to[0];
740 ppt0->n[1] = to[1];
741 ppt0->n[2] = to[2];
742
743 for (k=0; k<ilist; k++) {
744 iel = list[k] / 3;
745 i0 = list[k] % 3;
746 pt = &mesh->tria[iel];
747 pt0 = &mesh->tria[0];
748 memcpy(pt0,pt,sizeof(MMG5_Tria));
749 pt0->v[i0] = 0;
750 calold = caleltsig_iso(mesh,NULL,iel);
751 calnew = caleltsig_iso(mesh,NULL,0);
752
753 if ( (calnew < 0.001) && (calnew<calold) ) return 0;
754 //if ( chkedg(mesh,0) ) return 0;
755 }
756
757 /* coordinates, normals, tangents update */
758 p0->c[0] = o[0];
759 p0->c[1] = o[1];
760 p0->c[2] = o[2];
761
762 mesh->xpoint[p0->xp].n1[0] = nn1[0];
763 mesh->xpoint[p0->xp].n1[1] = nn1[1];
764 mesh->xpoint[p0->xp].n1[2] = nn1[2];
765
766 if ( isrid ) {
767 mesh->xpoint[p0->xp].n2[0] = nn2[0];
768 mesh->xpoint[p0->xp].n2[1] = nn2[1];
769 mesh->xpoint[p0->xp].n2[2] = nn2[2];
770 }
771 p0->n[0] = to[0];
772 p0->n[1] = to[1];
773 p0->n[2] = to[2];
774
775 return 1;
776}
ier
int ier
Definition: API_functionsf_2d.c:743
mesh
MMG5_pMesh * mesh
Definition: API_functionsf_2d.c:66
MMGS_bezierInt
int MMGS_bezierInt(MMG5_pBezier pb, double uv[2], double o[3], double no[3], double to[3])
Definition: bezier_s.c:207
MMG5_EPSD
#define MMG5_EPSD
Definition: eigenv_private.h:31
MMG5_EPS
#define MMG5_EPS
Definition: eigenv_private.h:32
libmmgs_private.h
caleltsig_iso
double caleltsig_iso(MMG5_pMesh mesh, MMG5_pSol met, MMG5_int iel)
Definition: quality_s.c:141
MS_SIN
#define MS_SIN(tag)
Definition: libmmgs_private.h:49
MG_GEO
#define MG_GEO
Definition: mmgcommon_private.h:145
MMG5_EPSOK
#define MMG5_EPSOK
Definition: mmgcommon_private.h:97
MG_MIN
#define MG_MIN(a, b)
Definition: mmgcommon_private.h:140
MG_EDG
#define MG_EDG(tag)
Definition: mmgcommon_private.h:171
MMG5_NULKAL
#define MMG5_NULKAL
Definition: mmgcommon_private.h:98
MMG5_TRIA_LMAX
#define MMG5_TRIA_LMAX
Definition: mmgcommon_private.h:59
MMG5_inxt2
static const uint8_t MMG5_inxt2[6]
Definition: mmgcommon_private.h:583
MMG5_EPSD2
#define MMG5_EPSD2
Definition: mmgcommon_private.h:95
mmgexterns_private.h
MMGS_moveTowardPoint
int MMGS_moveTowardPoint(MMG5_pMesh mesh, MMG5_pPoint p0, MMG5_pPoint p, double llold, double lam0, double lam1, double lam2, double nn1[3], double nn2[3], double to[3])
Definition: movpt_s.c:555
movridpt_iso
int movridpt_iso(MMG5_pMesh mesh, MMG5_pSol met, MMG5_int *list, int ilist)
Definition: movpt_s.c:595
movintpt_iso
int movintpt_iso(MMG5_pMesh mesh, MMG5_pSol met, MMG5_int *list, int ilist)
Definition: movpt_s.c:52
MMGS_update_normalAndTangent
static int MMGS_update_normalAndTangent(MMG5_pMesh mesh, MMG5_pPoint p0, MMG5_pPoint p, double llold, double lam0, double lam1, double lam2, double no1[3], double no2[3], double np1[3], double np2[3], double nn1[3], double nn2[3], double to[3])
Definition: movpt_s.c:449
MMGS_paramDisp
int MMGS_paramDisp(MMG5_pMesh mesh, MMG5_int it, int8_t isrid, MMG5_int ip0, MMG5_int ip, double step, double o[3])
Definition: movpt_s.c:375
MMG5_Bezier
Definition: mmgcommon_private.h:609
MMG5_Mesh
MMG mesh structure.
Definition: libmmgtypes.h:605
MMG5_Mesh::point
MMG5_pPoint point
Definition: libmmgtypes.h:641
MMG5_Mesh::xpoint
MMG5_pxPoint xpoint
Definition: libmmgtypes.h:642
MMG5_Mesh::tria
MMG5_pTria tria
Definition: libmmgtypes.h:647
MMG5_Point
Structure to store points of a MMG mesh.
Definition: libmmgtypes.h:270
MMG5_Point::n
double n[3]
Definition: libmmgtypes.h:272
MMG5_Point::tag
int16_t tag
Definition: libmmgtypes.h:284
MMG5_Point::c
double c[3]
Definition: libmmgtypes.h:271
MMG5_Point::xp
MMG5_int xp
Definition: libmmgtypes.h:279
MMG5_Sol
Definition: libmmgtypes.h:662
MMG5_Tria
Definition: libmmgtypes.h:332
MMG5_Tria::tag
int16_t tag[3]
Definition: libmmgtypes.h:342
MMG5_Tria::v
MMG5_int v[3]
Definition: libmmgtypes.h:334
MMG5_xPoint
Structure to store surface points of a MMG mesh.
Definition: libmmgtypes.h:294
MMG5_xPoint::n2
double n2[3]
Definition: libmmgtypes.h:295
MMG5_xPoint::n1
double n1[3]
Definition: libmmgtypes.h:295
Generated on Tue Dec 13 2022 00:19:54 for Mmg by doxygen 1.9.5