1 /*
2 * ***************************************************************************
3 * MC = < Manifold Code >
4 * Copyright (C) 1994-- Michael Holst
5 *
6 * This library is free software; you can redistribute it and/or
7 * modify it under the terms of the GNU Lesser General Public
8 * License as published by the Free Software Foundation; either
9 * version 2.1 of the License, or (at your option) any later version.
10 *
11 * This library is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 * Lesser General Public License for more details.
15 *
16 * You should have received a copy of the GNU Lesser General Public
17 * License along with this library; if not, write to the Free Software
18 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 *
20 * rcsid="$Id: gemg.c,v 1.17 2010/08/12 05:18:56 fetk Exp $"
21 * ***************************************************************************
22 */
23
24 /*
25 * ***************************************************************************
26 * File: Gemg.c
27 *
28 * Purpose: Class Gem: methods.
29 *
30 * Author: Michael Holst
31 * ***************************************************************************
32 */
33
34 #include "gem_p.h"
35
36 VEMBED(rcsid="$Id: gemg.c,v 1.17 2010/08/12 05:18:56 fetk Exp $")
37
38 /*
39 * ***************************************************************************
40 * Class Gem: Inlineable methods
41 * ***************************************************************************
42 */
43 #if !defined(VINLINE_GEM)
44
45 #endif /* if !defined(VINLINE_GEM) */
46 /*
47 * ***************************************************************************
48 * Class Gem: Non-inlineable methods
49 * ***************************************************************************
50 */
51
52 /*
53 * ***************************************************************************
54 * Routine: Gem_simplexInfo
55 *
56 * Purpose: Build the basic master-to-element transformation information.
57 *
58 * Notes:
59 *
60 * gchart ==> unified common chart for vertex coordinates
61 * chart[4] ==> individual charts for vertex coordinates
62 *
63 * D ==> jacobian determinant of the transformation
64 * Dcook ==> jacobian determinant of cooked trans
65 * faceD[4] ==> face jacobian determinants
66 *
67 * ff[3][3], bb[3] ==> affine trans from master to arbitrary el
68 * gg[3][3], cc[3] ==> affine trans from arbitrary el to master
69 *
70 * loc[4][3] ==> local ordering of vertices for each face
71 * vx[4][3] ==> vertex coordinate labels
72 * nvec[4][3] ==> normal vectors to the faces
73 * evec[6][3] ==> edge vectors
74 * elen[6] ==> edge vector lengths
75 *
76 * dimV ==> number of vertices in the d-simplex
77 * dimE ==> number of edges in the d-simplex
78 * dimF ==> number of faces in the d-simplex
79 * dimS ==> number of simplices in the d-simplex (=1)
80 *
81 * sid ==> global simplex ID
82 * vid[4] ==> global vertex IDs
83 * fid[4] ==> global face IDs
84 * eid[6] ==> global edge IDs
85 *
86 * stype ==> global simplex type
87 * vtype[4] ==> global vertex types
88 * ftype[4] ==> global face types
89 * etype[6] ==> LOCAL edge types
90 *
91 * *s ==> pointer to the simplex
92 * *v[4] ==> pointers to vertices of the simplex
93 *
94 * Author: Michael Holst
95 * ***************************************************************************
96 */
Gem_simplexInfo(Gem * thee,SS * sm,TT * t)97 VPUBLIC void Gem_simplexInfo(Gem *thee, SS *sm, TT *t)
98 {
99 unsigned int i, j;
100 double area, masterArea, vx[4][3], v0[3], v1[3], v2[3];
101
102 /* this simplex */
103 t->s = sm;
104
105 /* basic simplex info */
106 t->dimV = SS_dimVV(t->s);
107 t->dimE = SS_dimEE(t->s);
108 t->dimF = SS_dimVV(t->s);
109 t->dimS = 1;
110 t->D = 1.;
111 t->sid = SS_id(t->s);
112 t->stype = SS_type(t->s);
113
114 /* vertex and face info */
115 for (j=0; j<3; j++) {
116 t->bc[j] = 0.;
117 }
118 for (i=0; i<SS_dimVV(t->s); i++) {
119
120 t->v[i] = SS_vertex(t->s,i);
121
122 t->chart[i] = VV_chart(t->v[i]);
123 t->faceD[i] = 1.;
124
125 t->vid[i] = VV_id(t->v[i]);
126 t->vtype[i] = VV_type(t->v[i]);
127
128 # if defined(VG_ELEMENT)
129 t->fid[i] = SS_faceNumber(t->s,i);
130 # endif
131
132 t->ftype[i] = SS_faceType(t->s,i);
133
134 for (j=0; j<3; j++) {
135 vx[i][j] = VV_coord(t->v[i],j);
136 t->bc[j] += vx[i][j];
137 t->nvec[i][j] = 0.;
138 }
139 for (j=0; j<SS_dim(t->s); j++)
140 t->loc[i][j] = vmapF[i][j];
141 }
142 for (j=0; j<3; j++) {
143 t->bc[j] /= SS_dimVV(t->s);
144 }
145
146 /* edge info */
147 for (i=0; i<SS_dimEE(t->s); i++) {
148 # if defined(VG_ELEMENT)
149 t->eid[i] = SS_edgeNumber(t->s,i);
150 # endif
151 t->etype[i] = 0;
152 /* t->etype[i] = SS_edgeType(t->s,i); */
153 /*
154 * *** If we have the edge list rather than explicit storage of edge
155 * *** info in the simplex, then get the edge info as follows:
156 * j = vmapEI[i][0];
157 * k = vmapEI[i][1];
158 * eg = SS_edge(t->s,j,k);
159 * t->eid[i] = EE_id( eg );
160 * t->etype[i] = EE_type( eg );
161 */
162 }
163
164 /* call oneChart to get the vertex coordinates w.r.t. a single chart */
165 /* (IF the user is using multiple coordinate systems...) */
166 thee->pde->oneChart(Gem_dim(thee), Gem_dimII(thee), t->stype,
167 t->chart, vx, SS_dimVV(t->s));
168
169 /* chart choice for all vertices of this simplex */
170 t->gchart = t->chart[0];
171
172 /* set vertex coordinates for simplex now w.r.t. a single chart */
173 for (i=0; i<SS_dimVV(t->s); i++)
174 for (j=0; j<3; j++)
175 t->vx[i][j] = vx[i][j];
176
177 /* 2-simplex: get triangle area plus (virtual) 4th vertex along the way */
178 if (SS_dim(t->s) == 2) {
179
180 /* produce a unit normal */
181 for (i=0; i<3; i++) {
182 v0[i] = t->vx[1][i] - t->vx[0][i];
183 v1[i] = t->vx[2][i] - t->vx[0][i];
184 }
185 Vec3_xcry(v2, v0, v1);
186 area = Vec3_nrm2(v2);
187 Vec3_scal(v2, 1./area);
188
189 /* triangle area is 1/2 parallelogram area given by cross-product */
190 area *= 0.5;
191
192 /* transformation jacobian is an area ratio: area(this)/area(master) */
193 masterArea = 0.5;
194 t->D = area / masterArea;
195 /* VASSERT( t->D != 0. ); */
196 /* VASSERT( t->D > 0. ); */
197
198 /* cookup a "good" (virtual) fourth vertex from baricenter up normal */
199 /* first get the baricenter */
200 for (j=0; j<3; j++) {
201 t->vx[3][j] = 0.;
202 for (i=0; i<3; i++)
203 t->vx[3][j] += t->vx[i][j];
204 t->vx[3][j] /= 3.;
205 }
206 /* now move up the normal */
207 for (j=0; j<3; j++)
208 t->vx[3][j] += v2[j];
209 }
210
211 /* we have 4 vertex coords of a (some) tetrahedron (all in one chart) */
212 /* SO, now get all 6 edge vectors and 6 edge lengths (squared) */
213 for (i=0; i<6; i++) {
214 t->elen[i] = 0.0;
215 for (j=0; j<3; j++) {
216 t->evec[i][j] = t->vx[vmapEI[i][1]][j] - t->vx[vmapEI[i][0]][j];
217 t->elen[i] += ( t->evec[i][j]*t->evec[i][j] );
218 }
219 }
220
221 /* initialize transformations to/from master simplex and this simplex */
222 for (i=0; i<3; i++) {
223 t->bb[i] = t->vx[0][i];
224 t->cc[i] = 0.0;
225 for (j=0; j<3; j++) t->gg[i][j] = 0.;
226 }
227 t->ff[0][0] = t->evec[ vmapE[1][0] ][0];
228 t->ff[1][0] = t->evec[ vmapE[1][0] ][1];
229 t->ff[2][0] = t->evec[ vmapE[1][0] ][2];
230 t->ff[0][1] = t->evec[ vmapE[2][0] ][0];
231 t->ff[1][1] = t->evec[ vmapE[2][0] ][1];
232 t->ff[2][1] = t->evec[ vmapE[2][0] ][2];
233 t->ff[0][2] = t->evec[ vmapE[3][0] ][0];
234 t->ff[1][2] = t->evec[ vmapE[3][0] ][1];
235 t->ff[2][2] = t->evec[ vmapE[3][0] ][2];
236
237 /* jacobian of the transformation */
238 t->Dcook = (t->ff[0][0]*t->ff[1][1]*t->ff[2][2])
239 + (t->ff[0][1]*t->ff[1][2]*t->ff[2][0])
240 + (t->ff[0][2]*t->ff[1][0]*t->ff[2][1])
241 - (t->ff[2][0]*t->ff[1][1]*t->ff[0][2])
242 - (t->ff[1][2]*t->ff[2][1]*t->ff[0][0])
243 - (t->ff[2][2]*t->ff[1][0]*t->ff[0][1]);
244
245 /* 3-simplex: tetrahedron volume is jacobian of the transformation */
246 if (SS_dim(t->s) == 3) t->D = t->Dcook;
247
248 /* 2D PLANER HACK ONLY: make sign reflect orientation of triangle */
249 #if 1
250 if (SS_dim(t->s) == 2)
251 if ( (t->ff[0][0]*t->ff[1][1] - t->ff[1][0]*t->ff[0][1]) < 0. )
252 t->D = - t->D;
253 #endif
254 }
255
256 /*
257 * ***************************************************************************
258 * Routine: Gem_buildVolumeTrans
259 *
260 * Purpose: Build the complete master-to-element transformation.
261 *
262 * Notes: We just call Gem_simplexInfo to build the basic information,
263 * and then we compute the transformation and some additional
264 * things like the inverse transformations and various determinants.
265 *
266 * Author: Michael Holst
267 * ***************************************************************************
268 */
Gem_buildVolumeTrans(Gem * thee,SS * sm,TT * t)269 VPUBLIC void Gem_buildVolumeTrans(Gem *thee, SS *sm, TT *t)
270 {
271 unsigned int i, j;
272 double Di;
273
274 /* get jacobian with vertex coordinates w.r.t. a single chart */
275 Gem_simplexInfo(thee, sm, t);
276 /* VASSERT( t->Dcook > 0. ); */
277 Di = 1. / t->Dcook;
278
279 /* linear part of inverse transformation from this simplex to master */
280 t->gg[0][0] = Di * (t->ff[1][1]*t->ff[2][2] - t->ff[2][1]*t->ff[1][2]);
281 t->gg[0][1] = Di * (t->ff[2][1]*t->ff[0][2] - t->ff[0][1]*t->ff[2][2]);
282 t->gg[0][2] = Di * (t->ff[0][1]*t->ff[1][2] - t->ff[1][1]*t->ff[0][2]);
283 t->gg[1][0] = Di * (t->ff[2][0]*t->ff[1][2] - t->ff[1][0]*t->ff[2][2]);
284 t->gg[1][1] = Di * (t->ff[0][0]*t->ff[2][2] - t->ff[2][0]*t->ff[0][2]);
285 t->gg[1][2] = Di * (t->ff[1][0]*t->ff[0][2] - t->ff[0][0]*t->ff[1][2]);
286 t->gg[2][0] = Di * (t->ff[1][0]*t->ff[2][1] - t->ff[2][0]*t->ff[1][1]);
287 t->gg[2][1] = Di * (t->ff[2][0]*t->ff[0][1] - t->ff[0][0]*t->ff[2][1]);
288 t->gg[2][2] = Di * (t->ff[0][0]*t->ff[1][1] - t->ff[1][0]*t->ff[0][1]);
289
290 /* affine part of inverse transformation from this simplex to master */
291 for (i=0; i<3; i++) {
292 t->cc[i] = 0.0;
293 for (j=0; j<3; j++)
294 t->cc[i] -= t->gg[i][j] * t->bb[j];
295 }
296 }
297
298 /*
299 * ***************************************************************************
300 * Routine: Gem_buildSurfaceTrans
301 *
302 * Purpose: Build the complete masterFace-to-elementFace transformation.
303 *
304 * Notes: We ASSUME that Gem_simplexInfo has already been called to build
305 * the basic information in the TT structure we are given.
306 * We then compute some additional things here for a SINGLE face
307 * specified by "iface", such as the inverse transformations and
308 * various determinants.
309 *
310 * We do not actually have to assume that Gem_buildVolumeTrans
311 * has been previously called, only that Gem_simplexInfo has
312 * been called. However, it seems that all situations that
313 * occur result in Gem_buildVolumeTrans being called for an
314 * element before Gem_buildSurfTrans is called on any face.
315 *
316 * Author: Michael Holst
317 * ***************************************************************************
318 */
Gem_buildSurfaceTrans(Gem * thee,int iface,TT * t)319 VPUBLIC void Gem_buildSurfaceTrans(Gem *thee, int iface, TT *t)
320 {
321 double length, masterLength, area, masterArea, v0[3], v1[3], v2[3];
322
323 /* consistency check */
324 VASSERT( (iface>=0) && (iface<(int)SS_dimVV(t->s)) );
325
326 /* Recover el num and local node nums in simplex forming edge; */
327 /* recall that iface is the local node num OPPOSITE the edge */
328 if (SS_dim(t->s) == 2) {
329
330 /* the unit normal vector */
331 v0[0] = t->vx[ t->loc[iface][0] ][0];
332 v0[1] = t->vx[ t->loc[iface][0] ][1];
333 v0[2] = t->vx[ t->loc[iface][0] ][2];
334
335 v1[0] = t->vx[ t->loc[iface][1] ][0];
336 v1[1] = t->vx[ t->loc[iface][1] ][1];
337 v1[2] = t->vx[ t->loc[iface][1] ][2];
338
339 /*
340 * ********************************************************************
341 * Approximate normal for a 3-vector in the plane of the face
342 * NOTE: requires that "loc[iface][2]" be defined; CURRENTLY NOT!
343 * v2[0] = t->vx[ t->loc[iface][2] ][0];
344 * v2[1] = t->vx[ t->loc[iface][2] ][1];
345 * v2[2] = t->vx[ t->loc[iface][2] ][2];
346 * Vec3_copy(normal, v0);
347 * Vec3_axpy(normal, v1, (1.0));
348 * Vec3_scal(normal, 0.5);
349 * Vec3_axpy(normal, v2, -1.);
350 * ********************************************************************
351 */
352 /* Exact normal for a 2-vector in the 2-d plane */
353 Vec3_copy(v2, v0);
354 Vec3_axpy(v2, v1, -1.);
355 length = Vec3_nrm2(v2);
356 t->nvec[iface][0] = -v2[1];
357 t->nvec[iface][1] = v2[0];
358 t->nvec[iface][2] = 0.0;
359 area = Vec3_nrm2(t->nvec[iface]);
360 /* VASSERT( area != 0.0 ); */
361 Vec3_scal(t->nvec[iface], 1./area);
362
363 /* master edge length depends on the edge */
364 if (iface == 0) {
365 masterLength = VSQRT(2.0);
366 } else {
367 masterLength = 1.0;
368 }
369
370 /* transformation jacobian is a length ratio: len(this)/len(master) */
371 t->faceD[iface] = length / masterLength;
372
373 } else { /* (SS_dim(t->s) == 3) */
374
375 /* the unit normal vector */
376 v0[0] = t->vx[ t->loc[iface][1] ][0] - t->vx[ t->loc[iface][0] ][0];
377 v0[1] = t->vx[ t->loc[iface][1] ][1] - t->vx[ t->loc[iface][0] ][1];
378 v0[2] = t->vx[ t->loc[iface][1] ][2] - t->vx[ t->loc[iface][0] ][2];
379 v1[0] = t->vx[ t->loc[iface][2] ][0] - t->vx[ t->loc[iface][0] ][0];
380 v1[1] = t->vx[ t->loc[iface][2] ][1] - t->vx[ t->loc[iface][0] ][1];
381 v1[2] = t->vx[ t->loc[iface][2] ][2] - t->vx[ t->loc[iface][0] ][2];
382 Vec3_xcry(v2, v0, v1);
383 area = Vec3_nrm2(v2);
384 Vec3_copy(t->nvec[iface], v2);
385 /* VASSERT( area != 0.0 ); */
386 Vec3_scal(t->nvec[iface], 1./area);
387
388 /* tri-face area is 1/2 parallelogram area given by cross-product */
389 area *= 0.5;
390
391 /* master tri-face area depends on the face */
392 if (iface == 0) {
393 masterArea = 0.5 * VSQRT(3.0);
394 } else {
395 masterArea = 0.5;
396 }
397
398 /* transformation jacobian is a volume ratio: vol(this)/vol(master) */
399 t->faceD[iface] = area / masterArea;
400 }
401 /* VASSERT ( (t->faceD[iface]) > 0. ); */
402 }
403
404 /*
405 * ***************************************************************************
406 * Routine: Gem_edgeLength
407 *
408 * Purpose: Calculate the edge lengths of a simplex.
409 *
410 * Author: Michael Holst
411 * ***************************************************************************
412 */
Gem_edgeLength(Gem * thee,VV * v0,VV * v1)413 VPUBLIC double Gem_edgeLength(Gem *thee, VV *v0, VV *v1)
414 {
415 int j, edgeType, chart[4];
416 double value, vx[4][3];
417
418 /* setup the chart and vertex info */
419 chart[0] = VV_chart(v0);
420 chart[1] = VV_chart(v1);
421 for (j=0; j<Gem_dimII(thee); j++) {
422 vx[0][j] = VV_coord(v0,j);
423 vx[1][j] = VV_coord(v1,j);
424 }
425
426 /* have no idea what the edge type really is */
427 edgeType = 0;
428
429 /* call oneChart to get vertex coordinates w.r.t. a single chart */
430 thee->pde->oneChart(Gem_dim(thee),Gem_dimII(thee), edgeType,chart,vx,2);
431
432 /* okay, now return the edge length */
433 value = 0.;
434 for (j=0; j<Gem_dimII(thee); j++)
435 value += VSQR( vx[0][j] - vx[1][j] );
436 return VSQRT( value );
437 }
438
439 /*
440 * ***************************************************************************
441 * Routine: Gem_newestVertex
442 *
443 * Purpose: Determine the edge of a simplex opposite the newest vertex.
444 *
445 * Author: Michael Holst
446 * ***************************************************************************
447 */
Gem_newestVertex(Gem * thee,SS * sm,int face,double * len)448 VPUBLIC int Gem_newestVertex(Gem *thee, SS *sm, int face, double *len)
449 {
450 unsigned int p;
451 int v0,v1,v2,localNV;
452 long v0_id,v1_id,v2_id,nv;
453
454 /* determine the newest vertex in this simplex */
455 nv = -1;
456 localNV = -1;
457 for (p=0; p<SS_dimVV(sm); p++) {
458 v0_id = VV_id( SS_vertex(sm, p) );
459 if (v0_id > nv) {
460 nv = v0_id;
461 localNV = p;
462 }
463 }
464 VASSERT( (localNV >= 0) );
465
466 /* assign one of the opposite edges (there is only one in 2D) */
467 v0 = vmapOV3[localNV][0];
468 v1 = vmapOV3[localNV][1];
469 if (Gem_dim(thee) == 3) {
470
471 /* get the third vertex of opposite face */
472 v2 = vmapOV3[localNV][2];
473
474 /* get the global vertex numbers */
475 v0_id = VV_id ( SS_vertex(sm, v0) );
476 v1_id = VV_id ( SS_vertex(sm, v1) );
477 v2_id = VV_id ( SS_vertex(sm, v2) );
478
479 /* if v0 is oldest of (v0,v1), check if v2 is older than v1 */
480 if (v0_id < v1_id) {
481 if (v2_id < v1_id) {
482 v1 = v2;
483 }
484
485 /* if v1 is oldest of (v0,v1), check if v2 is older than v0 */
486 } else {
487 if (v2_id < v0_id) {
488 v0 = v2;
489 }
490 }
491 }
492
493 /* return the result */
494 return vmapE[v0][v1];
495 }
496
497 /*
498 * ***************************************************************************
499 * Routine: Gem_longestEdge
500 *
501 * Purpose: Determine the longest edge of a simplex or a simplex face.
502 *
503 * Note: It is critical to have a consistent tie-breaking rule in order
504 * to guarantee that recursive refinement procedures:
505 * (1) produce conforming meshes (two simplices will refine the
506 * same edge of a shared face)
507 * (2) terminate in finite steps (due to (1)).
508 *
509 * If the following statement is true then A is longer than B:
510 * ( (len(A).ge.len(B)) AND NOT(len(B).eq.len(A)) )
511 *
512 * Here len(A) is the length of A, .eq. is a "fuzzy" equals,
513 * and .ge. is a "fuzzy" greater than or equal to. To define the
514 * fuzzy comparisons we start by defining "greater than or equal to"
515 * ( x.ge.y ) as the logical statement ( x*(1+TOL) >= y ).
516 * Using x.ge.y we then define (x.eq.y) as (x.ge.y AND y.ge.x).
517 *
518 * In the case of a tie, the longest is defined as the edge with
519 * the smallest global edge number (as determined by the global
520 * vertex numbers).
521 *
522 * Note that by expanding the .eq. definition and simplifying we
523 * can rewrite the test "A is longer than B" as follows:
524 * ( (len(A).ge.len(B)) AND NOT(len(B).ge.len(A)) )
525 *
526 * Author: Michael Holst and Stephen Bond
527 * ***************************************************************************
528 */
Gem_longestEdge(Gem * thee,SS * sm,int face,double * len)529 VPUBLIC int Gem_longestEdge(Gem *thee, SS *sm, int face, double *len)
530 {
531 unsigned int p,q,pp,qq,ii,jj,q_lim;
532 int localENumI,localENumJ,vxp_id,vxq_id;
533 long kk, vertexMax, globalENum;
534 double tmpLen,edgeLen,TOL;
535 VV *vxp, *vxq;
536
537 VASSERT( (face==-1) || ((Gem_dim(thee)==3) && (face >=0) && (face <= 3)) );
538
539 localENumI = -1;
540 localENumJ = -1;
541 vertexMax = 100000000; /* conservative max of vertex numbers */
542 globalENum = -1; /* the current longest edge global number */
543 edgeLen = -VVLARGE; /* the current longest edge length */
544 TOL = VVSMALL; /* tolerance for our floating point tests */
545
546 if (face == -1) {
547 q_lim = SS_dimVV(sm);
548 } else {
549 q_lim = SS_dimVV(sm) - 1;
550 }
551 for (q=1; q<q_lim; q++) {
552 for (p=0; p<q; p++) {
553
554 if (face == -1) {
555 pp = p;
556 qq = q;
557 } else {
558 pp = vmapF[face][p];
559 qq = vmapF[face][q];
560 }
561 vxp = SS_vertex(sm, pp);
562 vxq = SS_vertex(sm, qq);
563 vxp_id = VV_id( vxp );
564 vxq_id = VV_id( vxq );
565
566 tmpLen = Gem_edgeLength(thee, vxp, vxq);
567 VASSERT( tmpLen > 0. );
568
569 /* if ( len(this_edge).ge.len(curr_max_edge) ) */
570 if ( tmpLen*(1.+TOL) >= edgeLen ) {
571
572 ii = VMAX2( vxp_id, vxq_id );
573 jj = VMIN2( vxp_id, vxq_id );
574 kk = (vertexMax * ii) + jj;
575
576 /* if ( NOT(len(curr_max_edge).ge.len(this_edge)) */
577 /* OR num(this_edge) < num(curr_max_edge) ) */
578 if ( ( !(edgeLen*(1.+TOL) >= tmpLen) )
579 || (kk < globalENum) ) {
580
581 localENumI = pp;
582 localENumJ = qq;
583 globalENum = kk;
584 edgeLen = tmpLen;
585
586 }
587 }
588 }
589 }
590
591 *len = edgeLen;
592 VASSERT( (localENumI>=0) && (localENumJ>=0) );
593 return vmapE[localENumI][localENumJ];
594 }
595
596 /*
597 * ***************************************************************************
598 * Routine: Gem_shortestEdge
599 *
600 * Purpose: Determine the shortest edge of a simplex or a simplex face.
601 *
602 * Note: It is critical to have a consistent tie-breaking rule in order
603 * to guarantee that recursive refinement procedures:
604 * (1) produce conforming meshes (two simplices will refine the
605 * same edge of a shared face)
606 * (2) terminate in finite steps (due to (1)).
607 *
608 * If the following statement is true then A is shorter than B:
609 * ( (len(B).ge.len(A)) AND NOT(len(A).eq.len(B)) )
610 *
611 * Here len(A) is the length of A, .eq. is a "fuzzy" equals,
612 * and .ge. is a "fuzzy" greater than or equal to. To define the
613 * fuzzy comparisons we start by defining "greater than or equal to"
614 * ( x.ge.y ) as the logical statement ( x*(1+TOL) >= y ).
615 * Using x.ge.y we can define (x.eq.y) as (x.ge.y AND y.ge.x).
616 *
617 * In the case of a tie, the shortest is defined as the edge with
618 * the largest global edge number (as determined by the global
619 * vertex numbers).
620 *
621 * Note that by expanding the .eq. definition and simplifying we
622 * can rewrite the test "A is shorter than B" as follows:
623 * ( (len(B).ge.len(A)) AND NOT(len(A).ge.len(B)) )
624 *
625 * Author: Michael Holst and Stephen Bond
626 * ***************************************************************************
627 */
Gem_shortestEdge(Gem * thee,SS * sm,int face,double * len)628 VPUBLIC int Gem_shortestEdge(Gem *thee, SS *sm, int face, double *len)
629 {
630 unsigned int p,q,pp,qq,ii,jj,q_lim;
631 int localENumI,localENumJ,vxp_id,vxq_id;
632 long kk, vertexMax, globalENum;
633 double tmpLen,edgeLen,TOL;
634 VV *vxp, *vxq;
635
636 localENumI = -1;
637 localENumJ = -1;
638 vertexMax = 100000000; /* conservative max of vertex numbers */
639 globalENum = -1; /* the current longest edge global number */
640 edgeLen = VVLARGE; /* the current longest edge length */
641 TOL = VVSMALL; /* tolerance for our floating point tests */
642
643 if (face == -1) {
644 q_lim = SS_dimVV(sm);
645 } else {
646 q_lim = SS_dimVV(sm) - 1;
647 }
648 for (q=1; q<q_lim; q++) {
649 for (p=0; p<q; p++) {
650
651 if (face == -1) {
652 pp = p;
653 qq = q;
654 } else {
655 pp = vmapF[face][p];
656 qq = vmapF[face][q];
657 }
658 vxp = SS_vertex(sm, pp);
659 vxq = SS_vertex(sm, qq);
660 vxp_id = VV_id( vxp );
661 vxq_id = VV_id( vxq );
662
663 tmpLen = Gem_edgeLength(thee, vxp, vxq);
664 VASSERT( tmpLen > 0. );
665
666 /* if ( len(curr_min_edge).ge.len(this_edge) ) */
667 if ( edgeLen*(1.+TOL) >= tmpLen ) {
668
669 ii = VMAX2( vxp_id, vxq_id );
670 jj = VMIN2( vxp_id, vxq_id );
671 kk = (vertexMax * ii) + jj;
672
673 /* if ( NOT(len(this_edge).ge.len(curr_min_edge)) */
674 /* OR num(this_edge) > num(curr_min_edge) ) */
675 if ( ( !( tmpLen*(1.+TOL) >= edgeLen ) )
676 || (kk > globalENum) ) {
677 localENumI = pp;
678 localENumJ = qq;
679 globalENum = kk;
680 edgeLen = tmpLen;
681 }
682 }
683 }
684 }
685 *len = edgeLen;
686 VASSERT( (localENumI>=0) && (localENumJ>=0) );
687 return vmapE[localENumI][localENumJ];
688 }
689
690 /*
691 * ***************************************************************************
692 * Routine: Gem_edgeRatio
693 *
694 * Purpose: Calculate ratio of longest-to-shortest edge of a simplex.
695 *
696 * Author: Michael Holst
697 * ***************************************************************************
698 */
Gem_edgeRatio(Gem * thee,SS * sm)699 VPUBLIC double Gem_edgeRatio(Gem *thee, SS *sm)
700 {
701 double lLong, lShort;
702 Gem_longestEdge(thee, sm, -1, &lLong);
703 Gem_shortestEdge(thee, sm, -1, &lShort);
704 return lLong / lShort;
705 }
706
707 /*
708 * ***************************************************************************
709 * Routine: Gem_simplexVolume
710 *
711 * Purpose: Calculate the determinant of the transformation from the
712 * master element to this element.
713 *
714 * Notes: We just call Gem_simplexInfo to build the basic information,
715 * and then we build the transformation and compute the determinant
716 * here.
717 *
718 * For a manifold, we need to know the orientation of the manifold
719 * in order to decide what is counter-clock-wise, and what is
720 * clockwise, in terms of vertex orderings. One will lead to a
721 * positive volume, and the other to a negative volume, when we
722 * compute volume using the determinant of the jacobian of the
723 * affine transformation to the master element. We assume here
724 * that the uniform chart computed by Gem_simplexInfo is such that
725 * the orientation in the chart reflects the correct orientation
726 * of the manifold (locally).
727 *
728 * Author: Michael Holst
729 * ***************************************************************************
730 */
Gem_simplexVolume(Gem * thee,SS * sm)731 VPUBLIC double Gem_simplexVolume(Gem *thee, SS *sm)
732 {
733 TT t;
734 Gem_simplexInfo(thee, sm, &t);
735 return t.D;
736 }
737
738 /*
739 * ***************************************************************************
740 * Routine: Gem_shapeChk
741 *
742 * Purpose: Traverse the simplices and check their shapes.
743 *
744 * Author: Michael Holst
745 * ***************************************************************************
746 */
Gem_shapeChk(Gem * thee)747 VPUBLIC void Gem_shapeChk(Gem *thee)
748 {
749 int i;
750 double f,g;
751 double tmp0,tmp1,tmp2,tmp3,tmp4;
752 double min0,min1,min2,min3,min4;
753 double max0,max1,max2,max3,max4;
754 SS *sm;
755
756 /* go through all simplices and setup the initial edge markings */
757 min0 = VVLARGE;
758 max0 = VVSMALL;
759 min1 = VVLARGE;
760 max1 = VVSMALL;
761 min2 = VVLARGE;
762 max2 = VVSMALL;
763 min3 = VVLARGE;
764 max3 = VVSMALL;
765 min4 = VVLARGE;
766 max4 = VVSMALL;
767 for (i=0; i<Gem_numSS(thee); i++) {
768 sm = Gem_SS(thee,i);
769 tmp0 = Gem_shapeMeasure( thee, sm, &f, &g );
770 tmp1 = Gem_edgeRatio( thee, sm );
771 tmp2 = Gem_simplexVolume( thee, sm );
772 Gem_shortestEdge(thee, sm, -1, &tmp3);
773 Gem_longestEdge(thee, sm, -1, &tmp4);
774 if (tmp0 < min0) min0 = tmp0;
775 if (tmp0 > max0) max0 = tmp0;
776 if (tmp1 < min1) min1 = tmp1;
777 if (tmp1 > max1) max1 = tmp1;
778 if (tmp2 < min2) min2 = tmp2;
779 if (tmp2 > max2) max2 = tmp2;
780 if (tmp3 < min3) min3 = tmp3;
781 if (tmp3 > max3) max3 = tmp3;
782 if (tmp4 < min4) min4 = tmp4;
783 if (tmp4 > max4) max4 = tmp4;
784 }
785 Vnm_print(0,"Gem_shapeChk: Shape statistics for <%d> simplices:\n",
786 Gem_numSS(thee));
787 Vnm_print(0,"Gem_shapeChk: Joe-Liu: [best,worst] = [%8.2e,%8.2e]\n",
788 max0,min0);
789 Vnm_print(0,"Gem_shapeChk: Edge-Ratio: [best,worst] = [%8.2e,%8.2e]\n",
790 min1,max1);
791 Vnm_print(0,"Gem_shapeChk: Volume: [small,large] = [%8.2e,%8.2e]\n",
792 min2,max2);
793 Vnm_print(0,"Gem_shapeChk: Short Edge: [small,large] = [%8.2e,%8.2e]\n",
794 min3,max3);
795 Vnm_print(0,"Gem_shapeChk: Long Edge: [small,large] = [%8.2e,%8.2e]\n",
796 min4,max4);
797 }
798
799 /*
800 * ***************************************************************************
801 * Routine: Gem_markEdges
802 *
803 * Purpose: Produce the initial edge markings in all of the simplices.
804 *
805 * Author: Michael Holst
806 * ***************************************************************************
807 */
Gem_markEdges(Gem * thee)808 VPUBLIC void Gem_markEdges(Gem *thee)
809 {
810 int i, refEg, ref_0, ref_1, ref_2, ref_3;
811 double len, len_0, len_1, len_2, len_3;
812 SS *sm;
813 VV *v0, *v1, *v2, *v3;
814
815 /* go through all simplices and setup the initial edge markings */
816 Vnm_print(0,"Gem_markEdges: marking edges..");
817
818 for (i=0; i<Gem_numSS(thee); i++) {
819 sm = Gem_SS(thee,i);
820
821 if (Gem_dim(thee) == 2) {
822
823 refEg = Gem_longestEdge(thee,sm,-1,&len);
824 SS_setRefinementEdge(sm,refEg);
825
826 } else if (Gem_dim(thee) == 3) {
827
828 /* (3D ONLY -- for Arnold/Mukherjee refinment) */
829 /* determine longest edge of each face; these are marked edges */
830
831 /* FACE 0 (face opposite vertex 0) */
832 v1 = SS_vertex(sm, vmapF[0][0]);
833 v2 = SS_vertex(sm, vmapF[0][1]);
834 v3 = SS_vertex(sm, vmapF[0][2]);
835 ref_0 = Gem_longestEdge(thee,sm,0,&len_0);
836
837 /* FACE 1 (face opposite vertex 1) */
838 v2 = SS_vertex(sm, vmapF[1][0]);
839 v0 = SS_vertex(sm, vmapF[1][1]);
840 v3 = SS_vertex(sm, vmapF[1][2]);
841 ref_1 = Gem_longestEdge(thee,sm,1,&len_1);
842
843 /* FACE 2 (face opposite vertex 2) */
844 v0 = SS_vertex(sm, vmapF[2][0]);
845 v1 = SS_vertex(sm, vmapF[2][1]);
846 v3 = SS_vertex(sm, vmapF[2][2]);
847 ref_2 = Gem_longestEdge(thee,sm,2,&len_2);
848
849 /* FACE 3 (face opposite vertex 3) */
850 v0 = SS_vertex(sm, vmapF[3][0]);
851 v2 = SS_vertex(sm, vmapF[3][1]);
852 v1 = SS_vertex(sm, vmapF[3][2]);
853 ref_3 = Gem_longestEdge(thee,sm,3,&len_3);
854
855 /* now mark edges appropriately */
856 refEg = Gem_longestEdge(thee,sm,-1,&len);
857 if (refEg == ref_0) {
858 SS_setRefinementEdge(sm,ref_0);
859 SS_setMarkedEdge1(sm,ref_1);
860 SS_setMarkedEdge2(sm,ref_2);
861 SS_setMarkedEdge3(sm,ref_3);
862 } else if (refEg == ref_1) {
863 SS_setRefinementEdge(sm,ref_1);
864 SS_setMarkedEdge1(sm,ref_0);
865 SS_setMarkedEdge2(sm,ref_2);
866 SS_setMarkedEdge3(sm,ref_3);
867 } else if (refEg == ref_2) {
868 SS_setRefinementEdge(sm,ref_2);
869 SS_setMarkedEdge1(sm,ref_0);
870 SS_setMarkedEdge2(sm,ref_1);
871 SS_setMarkedEdge3(sm,ref_3);
872 } else if (refEg == ref_3) {
873 SS_setRefinementEdge(sm,ref_3);
874 SS_setMarkedEdge1(sm,ref_0);
875 SS_setMarkedEdge2(sm,ref_1);
876 SS_setMarkedEdge3(sm,ref_2);
877 } else {
878 VASSERT( 0 );
879 }
880 #if 0
881 Vnm_print(2,"Gem_markEdges: RefEg = <%d> <%le>\n",refEg,len);
882 Vnm_print(2,"Gem_markEdges: ref_0 = <%d> <%le>\n",ref_0,len_0);
883 Vnm_print(2,"Gem_markEdges: ref_1 = <%d> <%le>\n",ref_1,len_1);
884 Vnm_print(2,"Gem_markEdges: ref_2 = <%d> <%le>\n",ref_2,len_2);
885 Vnm_print(2,"Gem_markEdges: ref_3 = <%d> <%le>\n",ref_3,len_3);
886 #endif
887 }
888 }
889
890 Vnm_print(0,"..done.\n");
891 }
892
893 /*
894 * ***************************************************************************
895 * Routine: Gem_reorderSV
896 *
897 * Purpose: Go through simplices and enforce a vertex ordering that
898 * will produce a positive determinant.
899 *
900 * Notes: This relies on an embedding of R2 into R3 (or R3 into R4)
901 * and breaks e.g. if this is a non-orientable 2-manifold, etc.
902 *
903 * Author: Michael Holst
904 * ***************************************************************************
905 */
Gem_reorderSV(Gem * thee)906 VPUBLIC void Gem_reorderSV(Gem *thee)
907 {
908 SS *sm;
909 int i;
910
911 /* go through all simplices and reorder vertices if necessary */
912 Vnm_print(0,"Gem_simpVertReorder: reordering vertices in each simplex..");
913
914 for (i=0; i<Gem_numSS(thee); i++) {
915 sm = Gem_SS(thee,i);
916 if (Gem_simplexVolume(thee,sm) < 0.) {
917 SS_reverse(sm);
918 }
919 VASSERT( Gem_simplexVolume(thee,sm) > 0.);
920 }
921 Vnm_print(0,"..done.\n");
922 }
923
924 /*
925 * ***************************************************************************
926 * Routine: Gem_smoothMeshLaplace
927 *
928 * Purpose: Smooth the mesh using simple Laplace smoothing.
929 *
930 * Notes: We don't need the simplex rings here, but we need the edge rings.
931 *
932 * Author: Michael Holst
933 * ***************************************************************************
934 */
Gem_smoothMeshLaplace(Gem * thee)935 VPUBLIC void Gem_smoothMeshLaplace(Gem *thee)
936 {
937 double v[50][3], w[50], x[3];
938 int i, numNabors;
939 VV *vx, *vx2;
940 EE *eg;
941
942 Gem_createEdges(thee);
943
944 /* go through all vertices (can go forwards now!) and relax hole nodes */
945 for (vx=Gem_firstVV(thee); vx!=VNULL; vx=Gem_nextVV(thee)) {
946
947 /* get coordinates of this vertex */
948 #if 0
949 for (i=0; i<Gem_dimII(thee); i++)
950 c[i] = VV_coord(vx,i);
951 #endif
952
953 /* initialize */
954 numNabors=0;
955 for (i=0; i<50; i++) w[i] = 0.0;
956
957 /* get coordinates of all vertices which share edges */
958 for (eg=VV_firstEE(vx); eg!=VNULL; eg=EE_link(eg,vx)) {
959 vx2 = EE_otherVertex(eg,vx);
960 w[numNabors] = 1.0;
961 for (i=0; i<Gem_dimII(thee); i++)
962 v[numNabors][i] = VV_coord(vx2,i);
963 numNabors++;
964 VASSERT( numNabors <= 50 );
965 }
966
967 /* if interior node, do a simple smoothing */
968 if ( VINTERIOR( VV_type(vx) ) ) {
969
970 /* do a simple averaging of the vertex positions */
971 #if 0
972 Vec3_copy(x, c);
973 for (i=0; i<numNabors; i++)
974 Vec3_axpy(x, &(v[i][0]), w[i]);
975 Vec3_scal(x, (1.0/(numNabors+1)));
976 #endif
977 Vec3_init(x, 0.0);
978 for (i=0; i<numNabors; i++)
979 Vec3_axpy(x, &(v[i][0]), w[i]);
980 Vec3_scal( x, (1.0/numNabors) );
981
982 /* set the new vertex positions */
983 for (i=0; i<Gem_dimII(thee); i++)
984 VV_setCoord(vx, i, x[i]);
985
986 } else {
987
988 } /* end if this is a vertex we need to smooth */
989 } /* end loop over vertices */
990
991 Gem_destroyEdges(thee);
992 }
993
994 /*
995 * ***************************************************************************
996 * Routine: Gem_smoothMesh
997 *
998 * Purpose: Smooth the mesh.
999 *
1000 * Author: Michael Holst
1001 * ***************************************************************************
1002 */
1003 #define VRINGMAX 1000
Gem_smoothMesh(Gem * thee)1004 VPUBLIC void Gem_smoothMesh(Gem *thee)
1005 {
1006 int i, j, k, numSring, maxSring, maxVring, maxVBring, numVring, numVBring;
1007 int found;
1008 double c[3], p[3], bc[3], mass, totalMass;
1009 VV *vx, *vxt, *vring[VRINGMAX];
1010 /* VV *vbring[VRINGMAX]; */
1011 SS *sm, *sring[VRINGMAX];
1012
1013 /* go through all vertices and adjust them to improve the mesh */
1014 maxSring = 0;
1015 maxVring = 0;
1016 maxVBring = 0;
1017 for (vx=Gem_firstVV(thee); vx!=VNULL; vx=Gem_nextVV(thee)) {
1018
1019 /* get coordinates of this vertex */
1020 for (i=0; i<Gem_dimII(thee); i++)
1021 c[i] = VV_coord(vx,i);
1022
1023 /* get the ring of simplices around this vertex */
1024 numSring=0;
1025 sring[numSring] = VV_firstSS(vx);
1026 while (sring[numSring] != VNULL) {
1027 numSring++;
1028 sring[numSring] = SS_link(sring[numSring-1],vx);
1029 }
1030 VASSERT( numSring > 0 );
1031 VASSERT( numSring <= VRINGMAX );
1032 if (numSring > maxSring) maxSring = numSring;
1033
1034 /* get the ring of vertices around this vertex */
1035 numVring=0;
1036 numVBring=0;
1037 vring[numVring] = VNULL;
1038 /* vbring[numVBring] = VNULL; */
1039 for (i=0; i<numSring; i++) {
1040 sm = sring[i];
1041 for (j=0; j<Gem_dimVV(thee); j++) {
1042 vxt = SS_vertex(sm,j);
1043 found = 0;
1044 for (k=0; k<numVring; k++)
1045 if ((vxt == vring[k]) || (vxt == vx))
1046 found = 1;
1047 if (!found) {
1048 vring[numVring] = vxt;
1049 numVring++;
1050 if ( VBOUNDARY( VV_type(vxt) )
1051 && ( VV_type(vxt) == VV_type(vx) ) ) {
1052 /* vbring[numVBring] = vxt; */
1053 numVBring++;
1054 }
1055 }
1056 }
1057 }
1058 VASSERT( numVring > 0 );
1059 VASSERT( numVring <= VRINGMAX );
1060 VASSERT( numVBring <= VRINGMAX );
1061 if (numVring > maxVring) maxVring = numVring;
1062 if (numVBring > maxVBring) maxVBring = numVBring;
1063
1064 /* if interior node, do a simple smoothing */
1065 if ( VINTERIOR( VV_type(vx) ) ) {
1066
1067 /* find centroid of entire vertex cloud */
1068 for (j=0; j<Gem_dimII(thee); j++) {
1069 bc[j] = 0.;
1070 totalMass = 0.;
1071 for (i=0; i<numVring; i++) {
1072 if ( VBOUNDARY( VV_type(vring[i]) )
1073 && (VV_type(vring[i]) > 1) ) {
1074 mass = 0.66;
1075 } else {
1076 mass = 1.0;
1077 }
1078 bc[j] += mass * VV_coord(vring[i],j);
1079 totalMass += mass;
1080 }
1081 /* bc[j] /= totalMass; */
1082 p[j] = ( totalMass * c[j] + bc[j] ) / ( 2.0 * totalMass );
1083 VV_setCoord( vx, j, p[j] );
1084 }
1085
1086 } /* end if this is a vertex we need to smooth */
1087
1088 } /* end loop over vertices */
1089
1090 Vnm_print(0, "Gem_smoothMesh: maxSring = <%d>\n", maxSring);
1091 Vnm_print(0, "Gem_smoothMesh: maxVring = <%d>\n", maxVring);
1092 Vnm_print(0, "Gem_smoothMesh: maxVBring = <%d>\n", maxVBring);
1093 }
1094
1095 /*
1096 * ***************************************************************************
1097 * Routine: Gem_smoothMeshBnd
1098 *
1099 * Purpose: Smooth the boundary mesh.
1100 *
1101 * Author: Michael Holst
1102 * ***************************************************************************
1103 */
Gem_smoothMeshBnd(Gem * thee)1104 VPUBLIC void Gem_smoothMeshBnd(Gem *thee)
1105 {
1106 int i, j, k, numSring, maxSring, maxVring, maxVBring, numVring, numVBring;
1107 int found;
1108 double c[3], p[3], bc[3], mass, totalMass;
1109 VV *vx, *vxt, *vring[VRINGMAX];
1110 /* VV *vbring[VRINGMAX]; */
1111 SS *sm, *sring[VRINGMAX];
1112
1113 /* go through all vertices and adjust them to improve the mesh */
1114 maxSring = 0;
1115 maxVring = 0;
1116 maxVBring = 0;
1117 for (vx=Gem_firstVV(thee); vx!=VNULL; vx=Gem_nextVV(thee)) {
1118
1119 /* get coordinates of this vertex */
1120 for (i=0; i<Gem_dimII(thee); i++)
1121 c[i] = VV_coord(vx,i);
1122
1123 /* get the ring of simplices around this vertex */
1124 numSring=0;
1125 sring[numSring] = VV_firstSS(vx);
1126 while (sring[numSring] != VNULL) {
1127 numSring++;
1128 sring[numSring] = SS_link(sring[numSring-1],vx);
1129 }
1130 VASSERT( numSring > 0 );
1131 VASSERT( numSring <= VRINGMAX );
1132 if (numSring > maxSring) maxSring = numSring;
1133
1134 /* get the ring of vertices around this vertex */
1135 numVring=0;
1136 numVBring=0;
1137 vring[numVring] = VNULL;
1138 /* vbring[numVBring] = VNULL; */
1139 for (i=0; i<numSring; i++) {
1140 sm = sring[i];
1141 for (j=0; j<Gem_dimVV(thee); j++) {
1142 vxt = SS_vertex(sm,j);
1143 found = 0;
1144 for (k=0; k<numVring; k++)
1145 if ((vxt == vring[k]) || (vxt == vx))
1146 found = 1;
1147 if (!found) {
1148 vring[numVring] = vxt;
1149 numVring++;
1150 if ( VBOUNDARY( VV_type(vxt) )
1151 && ( VV_type(vxt) == VV_type(vx) ) ) {
1152 /* vbring[numVBring] = vxt; */
1153 numVBring++;
1154 }
1155 }
1156 }
1157 }
1158 VASSERT( numVring > 0 );
1159 VASSERT( numVring <= VRINGMAX );
1160 VASSERT( numVBring <= VRINGMAX );
1161 if (numVring > maxVring) maxVring = numVring;
1162 if (numVBring > maxVBring) maxVBring = numVBring;
1163
1164 /* if boundary node, do a simple smoothing */
1165 if ( VBOUNDARY( VV_type(vx) ) ) {
1166
1167 /* find centroid of entire vertex cloud */
1168 for (j=0; j<Gem_dimII(thee); j++) {
1169 bc[j] = 0.;
1170 totalMass = 0.;
1171 for (i=0; i<numVring; i++) {
1172 if ( VBOUNDARY( VV_type(vring[i]) ) ) {
1173 mass = 1.0;
1174 } else {
1175 mass = 0.66;
1176 }
1177 bc[j] += mass * VV_coord(vring[i],j);
1178 totalMass += mass;
1179 }
1180 /* bc[j] /= totalMass; */
1181 p[j] = ( totalMass * c[j] + bc[j] ) / ( 2.0 * totalMass );
1182 }
1183 thee->pde->mapBoundary( Gem_dim(thee), Gem_dimII(thee),
1184 VV_type(vx), VV_chart(vx), p );
1185 for (j=0; j<Gem_dimII(thee); j++)
1186 VV_setCoord( vx, j, p[j] );
1187
1188 } /* end if this is a vertex we need to smooth */
1189
1190 } /* end loop over vertices */
1191
1192 Vnm_print(0, "Gem_smoothMeshBnd: maxSring = <%d>\n", maxSring);
1193 Vnm_print(0, "Gem_smoothMeshBnd: maxVring = <%d>\n", maxVring);
1194 Vnm_print(0, "Gem_smoothMeshBnd: maxVBring = <%d>\n", maxVBring);
1195 }
1196
1197 /*
1198 * ***************************************************************************
1199 * Routine: Gem_shapeMeasure
1200 *
1201 * Purpose: Calculate the shape quality measure for this simplex.
1202 *
1203 * Notes: Let |s| denote the volume (which may be negative) of a given
1204 * d-simplex "s", let v_i (i=0,...,d) denote the vertices of s,
1205 * and let e_{ij} denote the d-vectors representing the 3 or 6 edges
1206 * of s that connect v_i to v_j. We compute the following shape
1207 * quality measure for the simplex s:
1208 *
1209 * f(s,d) 2^{2(1-1/d)} * 3^{(d-1)/2} * |s|^{2/d}
1210 * meas(s,d) = ------ = --------------------------------------
1211 * g(s,d) \sum_{0<=i<j<=d} |e_{ij}|^2
1212 *
1213 * 2D Notes: The shape function meas(s,2) is (nearly) the same one used by
1214 * Randy Bank and Kent Smith in their joint paper on mesh smoothing:
1215 *
1216 * f(s,2) 2 * 3^{1/2} * |s|
1217 * meas(s,2) = ------ = ---------------------------
1218 * g(s,2) \sum_{0<=i<j<=2} |e_{ij}|^2
1219 *
1220 * It has the property that its maximal value of 1 is obtained
1221 * for an equilateral triangle, and it is scaling invariant, i.e.,
1222 * we don't have to worry about the size of the triangle.
1223 * Their original function was normalized (their numerator was
1224 * 2*f(s,2) above) to yield a value of 1 for a equalateral triangle
1225 * (with volume 1); this is modified to yield a maximal value
1226 * of 1 for the unit triangle (with volume 1/2) to work better with
1227 * the unit triangle code. To effect this slightly different
1228 * normalization, the numerator of the quality function was changed
1229 * to 2(3)^{1/2}|s|.
1230 *
1231 * 3D Notes: The shape function meas(s,3) is (nearly) the same one used by
1232 * Joe and Liu in their paper on quality measures for tetrahedra:
1233 *
1234 * f(s,3) 2^{4/3} * 3 * |s|^{2/3}
1235 * meas(s,3) = ------ = ---------------------------
1236 * g(s,3) \sum_{0<=i<j<=3} |e_{ij}|^2
1237 *
1238 * It is also scaling invariant, so we don't have to worry about
1239 * the size of the tetrahedron. Their original function was
1240 * normalized (their numerator was 12(3|s|)^{2/3}) to yield a
1241 * maximal value of 1 for a regular tetrahedron (with volume 1);
1242 * this was modified to yield a maximal value of 1 for the unit
1243 * tetrahedron (with volume 1/6) to work better with the unit
1244 * tetrahedron code. To effect this slightly different
1245 * normalization, the numerator of the quality function was changed
1246 * to f(s,3) = 12(3|s|/6)^{2/3} = 2^{4/3}*3*|s|^{2/3}, as above.
1247 *
1248 * Author: Michael Holst
1249 * ***************************************************************************
1250 */
Gem_shapeMeasure(Gem * thee,SS * sm,double * f,double * g)1251 VPUBLIC double Gem_shapeMeasure(Gem *thee, SS *sm, double *f, double *g)
1252 {
1253 int i;
1254 double d, vol, siggn, value;
1255 TT t;
1256
1257 /* get simplex volume and its sign */
1258 Gem_simplexInfo(thee, sm, &t);
1259 siggn = VSIGN( 1., t.D );
1260 vol = VABS( t.D );
1261
1262 /* compute volume scaling */
1263 d = (double)Gem_dim(thee);
1264 (*f) = siggn * VPOW( 2., 2.*(1.-1./d) )
1265 * VPOW( 3., (d-1.)*.5 )
1266 * VPOW( vol, 2./d );
1267
1268 /* accumlate (squared) edge lengths */
1269 (*g) = 0.;
1270 for (i=0; i<Gem_dimEE(thee); i++)
1271 (*g) += t.elen[i];
1272
1273 if ((*g) > 0.) {
1274 value = ((*f) / (*g));
1275 } else {
1276 value = 0.;
1277 }
1278
1279 return value;
1280 }
1281
1282 /*
1283 * ***************************************************************************
1284 * Routine: Gem_shapeMeasureGrad
1285 *
1286 * Purpose: Calculate gradient of the shape quality measure for this simplex,
1287 * where the last vertex (vertex d) is treated as the set of
1288 * independent variables.
1289 *
1290 * Author: Michael Holst
1291 * ***************************************************************************
1292 */
Gem_shapeMeasureGrad(Gem * thee,SS * sm,int vmap[],double dm[])1293 VPUBLIC void Gem_shapeMeasureGrad(Gem *thee, SS *sm, int vmap[], double dm[])
1294 {
1295 int i, d;
1296 double f, g, vol, siggn, fac, df[3], dg[3], x[4], y[4], z[4];
1297 TT t;
1298
1299 d = Gem_dim(thee);
1300 (void)Gem_shapeMeasure(thee, sm, &f, &g);
1301
1302 if (d==2) {
1303 for (i=0; i<d+1; i++) {
1304 x[i] = VV_coord( SS_vertex( sm, vmap[i] ), 0 );
1305 y[i] = VV_coord( SS_vertex( sm, vmap[i] ), 1 );
1306 }
1307 fac = 2. * VSQRT(3.);
1308 df[0] = fac * ( y[0] - y[1] );
1309 df[1] = fac * ( x[1] - x[0] );
1310 df[2] = 0.;
1311 dg[0] = 2. * (x[2] - x[0]) + 2. * (x[2] - x[1]);
1312 dg[1] = 2. * (y[2] - y[0]) + 2. * (y[2] - y[1]);
1313 dg[2] = 0.;
1314 } else if (d==3) {
1315 for (i=0; i<d+1; i++) {
1316 x[i] = VV_coord( SS_vertex( sm, vmap[i] ), 0 );
1317 y[i] = VV_coord( SS_vertex( sm, vmap[i] ), 1 );
1318 z[i] = VV_coord( SS_vertex( sm, vmap[i] ), 2 );
1319 }
1320 Gem_simplexInfo(thee, sm, &t);
1321 siggn = VSIGN( 1., t.D );
1322 vol = VABS( t.D );
1323 fac = siggn * 2. * VPOW(2.,4./3.) * VPOW(vol,-5./3.);
1324 df[0] = fac * ( (y[1]-y[0])*(z[2]-z[0]) - (y[2]-y[0])*(z[1]-z[0]) );
1325 df[1] = fac * ( (x[2]-x[0])*(z[1]-z[0]) - (x[1]-x[0])*(z[2]-z[0]) );
1326 df[2] = fac * ( (x[1]-x[0])*(y[2]-y[0]) - (x[2]-x[0])*(y[1]-y[0]) );
1327 dg[0] = 2. * (x[3] - x[0]) + 2. * (x[3] - x[1]) + 2. * (x[3] - x[2]);
1328 dg[1] = 2. * (y[3] - y[0]) + 2. * (y[3] - y[1]) + 2. * (y[3] - y[2]);
1329 dg[2] = 2. * (z[3] - z[0]) + 2. * (z[3] - z[1]) + 2. * (z[3] - z[2]);
1330 } else { VASSERT(0); }
1331
1332 fac = 1. / (g*g);
1333 for (i=0; i<3; i++)
1334 dm[i] = fac * ( g*df[i] - f*dg[i] );
1335 }
1336
1337 /*
1338 * ***************************************************************************
1339 * Routine: Gem_smoothMeshOpt
1340 *
1341 * Purpose: Smooth the mesh using a volume optimization approach.
1342 *
1343 * Author: Michael Holst
1344 * ***************************************************************************
1345 */
Gem_smoothMeshOpt(Gem * thee)1346 VPUBLIC void Gem_smoothMeshOpt(Gem *thee)
1347 {
1348 int i, j, k, done, imap, numSring, minIdx, maxSring, maxVring;
1349 int iters, itmax, numVring, numVBring, maxVBring, found, map[4];
1350 double step, f, g, minShape, longEdgeLen;
1351 double dm[3], a[3], b[3], c[3], shape[VRINGMAX];
1352 VV *vx, *vxt, *vring[VRINGMAX];
1353 /* VV *vbring[VRINGMAX]; */
1354 SS *sm, *sring[VRINGMAX];
1355
1356 /* go through all vertices and adjust them to improve the mesh */
1357 maxSring = 0;
1358 maxVring = 0;
1359 maxVBring = 0;
1360 for (vx=Gem_firstVV(thee); vx!=VNULL; vx=Gem_nextVV(thee)) {
1361
1362 /* get coordinates of this vertex */
1363 for (i=0; i<Gem_dimII(thee); i++)
1364 c[i] = VV_coord(vx,i);
1365
1366 /* get the ring of simplices around this vertex */
1367 numSring=0;
1368 sring[numSring] = VV_firstSS(vx);
1369 while (sring[numSring] != VNULL) {
1370 numSring++;
1371 sring[numSring] = SS_link(sring[numSring-1],vx);
1372 }
1373 VASSERT( numSring > 0 );
1374 VASSERT( numSring <= VRINGMAX );
1375 if (numSring > maxSring) maxSring = numSring;
1376
1377 /* get the ring of vertices around this vertex */
1378 numVring=0;
1379 numVBring=0;
1380 vring[numVring] = VNULL;
1381 /* vbring[numVBring] = VNULL; */
1382 for (i=0; i<numSring; i++) {
1383 sm = sring[i];
1384 for (j=0; j<Gem_dimVV(thee); j++) {
1385 vxt = SS_vertex(sm,j);
1386 found = 0;
1387 for (k=0; k<numVring; k++)
1388 if ((vxt == vring[k]) || (vxt == vx))
1389 found = 1;
1390 if (!found) {
1391 vring[numVring] = vxt;
1392 numVring++;
1393 if ( VBOUNDARY( VV_type(vxt) )
1394 && ( VV_type(vxt) == VV_type(vx) ) ) {
1395 /* vbring[numVBring] = vxt; */
1396 numVBring++;
1397 }
1398 }
1399 }
1400 }
1401 VASSERT( numVring > 0 );
1402 VASSERT( numVring <= VRINGMAX );
1403 VASSERT( numVBring <= VRINGMAX );
1404 if (numVring > maxVring) maxVring = numVring;
1405 if (numVBring > maxVBring) maxVBring = numVBring;
1406
1407 /* calculate shape measures; collect vertices */
1408 minShape = 1.1;
1409 minIdx = -1;
1410 for (i=0; i<numSring; i++) {
1411 shape[i] = Gem_shapeMeasure(thee, sring[i], &f, &g);
1412 if (shape[i] < minShape) {
1413 minShape = shape[i];
1414 minIdx = i;
1415 }
1416 }
1417 VASSERT( minShape <= 1.0 );
1418 VASSERT( minIdx >= 0.0 );
1419
1420 /* collect info about the worst simplex */
1421 imap = -1;
1422 sm = sring[minIdx];
1423 for (i=0; i<Gem_dimVV(thee); i++) {
1424 vxt = SS_vertex(sm,i);
1425 if (vxt == vx) imap = i;
1426 }
1427 VASSERT( imap >= 0 );
1428
1429 /* make last vertex the movable one for simplicity */
1430 for (i=0; i<Gem_dimVV(thee); i++)
1431 map[i] = (imap+i+1) % Gem_dimVV(thee);
1432
1433 /* if interior node, do a simple smoothing */
1434 if ( VINTERIOR( VV_type(vx) ) ) {
1435
1436 /* get the current measure and gradient */
1437 Gem_shapeMeasureGrad(thee, sm, map, dm);
1438
1439 /* set the initial end-points of the bisection interval */
1440 (void)Gem_longestEdge(thee, sm, -1, &longEdgeLen);
1441 step = longEdgeLen;
1442 for (j=0; j<Gem_dimII(thee); j++) {
1443 a[j] = c[j];
1444 b[j] = a[j] + step * dm[j];
1445 VV_setCoord( vx, j, b[j] );
1446 }
1447
1448 /* bisect along gradient to find optimal shape for simplex */
1449 /* (find point the improves worst shape in the entire ring) */
1450 itmax = 8;
1451 iters = 0;
1452 done = 0;
1453 while (!done) {
1454 VJMPERR1( !Vnm_sigInt() );
1455 iters++;
1456 done = 1;
1457 for (j=0; j<numSring; j++) {
1458 shape[j] = Gem_shapeMeasure(thee, sring[j], &f, &g);
1459 if (shape[j] <= minShape) done = 0;
1460 }
1461 if (iters > itmax) {
1462 done = 1;
1463 for (j=0; j<Gem_dimII(thee); j++)
1464 VV_setCoord( vx, j, c[j] );
1465 } else {
1466 step *= 0.5;
1467 for (j=0; j<Gem_dimII(thee); j++) {
1468 b[j] = a[j] + step * dm[j];
1469 VV_setCoord( vx, j, b[j] );
1470 }
1471 }
1472 }
1473
1474 } else {
1475
1476 } /* end if this is a vertex we need to smooth */
1477
1478 } /* end loop over vertices */
1479
1480 VERROR1:
1481 Vnm_print(0, "Gem_smoothMeshOpt: maxSring = <%d>\n", maxSring);
1482 Vnm_print(0, "Gem_smoothMeshOpt: maxVring = <%d>\n", maxVring);
1483 Vnm_print(0, "Gem_smoothMeshOpt: maxVBring = <%d>\n", maxVBring);
1484 }
1485
1486 /*
1487 * ***************************************************************************
1488 * Routine: Gem_buildCharts
1489 *
1490 * Purpose: Unify the charts of vertices.
1491 *
1492 * Author: Michael Holst
1493 * ***************************************************************************
1494 */
Gem_buildCharts(Gem * thee)1495 VPUBLIC void Gem_buildCharts(Gem *thee)
1496 {
1497 int j, chartV;
1498 VV *vx;
1499 double v[4][3];
1500
1501 for (vx=Gem_firstVV(thee); vx!=VNULL; vx=Gem_nextVV(thee)) {
1502 chartV = VV_chart(vx);
1503 for (j=0; j<Gem_dimII(thee); j++)
1504 v[0][j] = VV_coord(vx,j);
1505 thee->pde->oneChart(Gem_dim(thee), Gem_dimII(thee), 1, &chartV, v, 1);
1506 VV_setChart(vx, chartV);
1507 }
1508 }
1509
1510 /*
1511 * ***************************************************************************
1512 * Routine: Gem_clearCharts
1513 *
1514 * Purpose: Unify the charts of vertices.
1515 *
1516 * Author: Michael Holst
1517 * ***************************************************************************
1518 */
Gem_clearCharts(Gem * thee)1519 VPUBLIC void Gem_clearCharts(Gem *thee)
1520 {
1521 VV *vx;
1522
1523 for (vx=Gem_firstVV(thee); vx!=VNULL; vx=Gem_nextVV(thee)) {
1524 VV_setChart(vx, 0);
1525 }
1526 }
1527
1528