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: gemgen.c,v 1.18 2010/08/12 05:18:56 fetk Exp $"
21 * ***************************************************************************
22 */
23
24 /*
25 * ***************************************************************************
26 * File: Gemgen.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: gemgen.c,v 1.18 2010/08/12 05:18:56 fetk Exp $")
37
38 VPRIVATE int predInit = 0;
39
40 /*
41 * ***************************************************************************
42 * Class Gem: Inlineable methods
43 * ***************************************************************************
44 */
45 #if !defined(VINLINE_GEM)
46
47 #endif /* if !defined(VINLINE_GEM) */
48 /*
49 * ***************************************************************************
50 * Class Gem: Non-inlineable methods
51 * ***************************************************************************
52 */
53
54 /*
55 * ***************************************************************************
56 * Routine: Gem_pointInSimplex
57 *
58 * Purpose: Determine whether or a not a point is in a simplex.
59 *
60 * Author: Michael Holst
61 * ***************************************************************************
62 */
Gem_pointInSimplex(Gem * thee,SS * sm,double x[])63 VPUBLIC int Gem_pointInSimplex(Gem *thee, SS *sm, double x[])
64 {
65 int i, j;
66 double a[3], b[3], c[3], d[3], bc[3], len, z;
67
68 /* first do a quick check */
69 for (j=0; j<3; j++) {
70 bc[j] = 0.;
71 }
72 for (i=0; i<SS_dimVV(sm); i++) {
73 for (j=0; j<3; j++) {
74 bc[j] += VV_coord(SS_vertex(sm,i),j);
75 }
76 }
77 for (j=0; j<3; j++) {
78 bc[j] /= SS_dimVV(sm);
79 }
80 (void)Gem_longestEdge(thee, sm, -1, &len);
81 z = len + VSMALL;
82 if ( (VABS(x[0]-bc[0])>z)||(VABS(x[1]-bc[1])>z)||(VABS(x[2]-bc[2])>z) ) {
83 return 0;
84 }
85
86 /* OK, that didn't work; now do a more careful check */
87 Gem_predinit(thee);
88
89 if (Gem_dim(thee)==2) {
90 for (j=0; j<2; j++) {
91 a[j] = VV_coord(SS_vertex(sm,0),j);
92 b[j] = VV_coord(SS_vertex(sm,1),j);
93 c[j] = VV_coord(SS_vertex(sm,2),j);
94 }
95 if ( (Vpred_orient2d(a,b,x) >= -VSMALL)
96 && (Vpred_orient2d(b,c,x) >= -VSMALL)
97 && (Vpred_orient2d(c,a,x) >= -VSMALL) ) {
98 return 1;
99 } else {
100 return 0;
101 }
102 } else if (Gem_dim(thee)==3) {
103 for (j=0; j<3; j++) {
104 a[j] = VV_coord(SS_vertex(sm,1),j);
105 b[j] = VV_coord(SS_vertex(sm,0),j);
106 c[j] = VV_coord(SS_vertex(sm,2),j);
107 d[j] = VV_coord(SS_vertex(sm,3),j);
108 }
109 if ( (Vpred_orient3d(a,b,c,x) >= -VSMALL)
110 && (Vpred_orient3d(a,c,d,x) >= -VSMALL)
111 && (Vpred_orient3d(a,d,b,x) >= -VSMALL)
112 && (Vpred_orient3d(b,d,c,x) >= -VSMALL) ) {
113 return 1;
114 } else {
115 return 0;
116 }
117 } else {
118 VASSERT(0);
119 return 0;
120 }
121 }
122
123 /*
124 * ***************************************************************************
125 * Routine: Gem_pointInSimplexVal
126 *
127 * Purpose: Evaluate basis functions at a point in a simplex.
128 *
129 * Author: Michael Holst
130 * ***************************************************************************
131 */
Gem_pointInSimplexVal(Gem * thee,SS * sm,double x[],double phi[],double phix[][3])132 VPUBLIC int Gem_pointInSimplexVal(Gem *thee, SS *sm, double x[],
133 double phi[], double phix[][3])
134 {
135 int i, j, found;
136 double xm[3];
137 TT t;
138
139 /* Affine transformation from master element to this element */
140 Gem_buildVolumeTrans(thee,sm,&t);
141
142 /* Map current point to master element */
143 for (i=0; i<Gem_dimII(thee); i++) {
144 xm[i] = t.cc[i];
145 for (j=0; j<Gem_dimII(thee); j++)
146 xm[i] += t.gg[i][j]*x[j];
147 }
148
149 /* Evaluate master element basis functions at the point */
150 phi[0] = 1. - xm[0] - xm[1] - xm[2];
151 phi[1] = xm[0];
152 phi[2] = xm[1];
153 phi[3] = xm[2];
154
155 /* Evaluate master element basis functions gradients at the point */
156 phix[0][0] = -1.;
157 phix[0][1] = -1.;
158 phix[0][2] = -1.;
159 phix[1][0] = 1.;
160 phix[1][1] = 0.;
161 phix[1][2] = 0.;
162 phix[2][0] = 0.;
163 phix[2][1] = 1.;
164 phix[2][2] = 0.;
165 phix[3][0] = 0.;
166 phix[3][1] = 0.;
167 phix[3][2] = 1.;
168
169 /* Check if point lies in master element */
170 found = 1;
171 for (j=0; j<Gem_dimVV(thee); j++) {
172 if ( phi[j] <= -VVSMALL ) found = 0;
173 }
174
175 /* return */
176 return found;
177 }
178
179 /*
180 * ***************************************************************************
181 * Routine: Gem_findSimplex
182 *
183 * Purpose: Find a simplex containing a given vertex.
184 *
185 * Author: Michael Holst
186 * ***************************************************************************
187 */
Gem_findSimplex(Gem * thee,VV * vx)188 VPUBLIC SS* Gem_findSimplex(Gem *thee, VV *vx)
189 {
190 int i, j;
191 double x[3];
192 SS *sm;
193
194 /* grab coordinates of the point */
195 for (j=0; j<Gem_dimII(thee); j++)
196 x[j] = VV_coord( vx, j );
197
198 /* find the simplex in which this vertex lies */
199 for (i=0; i<Gem_numSS(thee); i++) {
200 sm = Gem_SS(thee,i);
201 if ( Gem_pointInSimplex(thee, sm, x) )
202 return sm;
203 }
204 return VNULL;
205 }
206
207 /*
208 * ***************************************************************************
209 * Routine: Gem_predinit
210 *
211 * Purpose: Initialize the geometric predicates.
212 *
213 * Author: Michael Holst
214 * ***************************************************************************
215 */
Gem_predinit(Gem * thee)216 VPUBLIC void Gem_predinit(Gem *thee)
217 {
218 if ( !predInit ) {
219 predInit = 1;
220 Vpred_exactinit();
221 }
222 }
223
224 /*
225 * ***************************************************************************
226 * Routine: Gem_orient
227 *
228 * Purpose: Determine the orientation of the vertices in a simplex.
229 *
230 * Author: Michael Holst
231 * ***************************************************************************
232 */
Gem_orient(Gem * thee,SS * sm)233 VPUBLIC int Gem_orient(Gem *thee, SS *sm)
234 {
235 int j;
236 double a[3], b[3], c[3], d[3];
237
238 Gem_predinit(thee);
239
240 if (Gem_dim(thee)==2) {
241 for (j=0; j<2; j++) {
242 a[j] = VV_coord(SS_vertex(sm,0),j);
243 b[j] = VV_coord(SS_vertex(sm,1),j);
244 c[j] = VV_coord(SS_vertex(sm,2),j);
245 }
246 if ( Vpred_orient2d(a,b,c) >= VVSMALL ) {
247 return 1;
248 } else if ( Vpred_orient2d(a,b,c) <= -VVSMALL ) {
249 return -1;
250 } else {
251 return 0;
252 }
253 } else if (Gem_dim(thee)==3) {
254 for (j=0; j<3; j++) {
255 a[j] = VV_coord(SS_vertex(sm,1),j);
256 b[j] = VV_coord(SS_vertex(sm,0),j);
257 c[j] = VV_coord(SS_vertex(sm,2),j);
258 d[j] = VV_coord(SS_vertex(sm,3),j);
259 }
260 if ( Vpred_orient3d(a,b,c,d) >= VVSMALL ) {
261 return 1;
262 } else if ( Vpred_orient3d(a,b,c,d) <= -VVSMALL ) {
263 return -1;
264 } else {
265 return 0;
266 }
267 } else {
268 VASSERT(0);
269 return 0;
270 }
271 }
272
273 /*
274 * ***************************************************************************
275 * Routine: Gem_orientVal
276 *
277 * Purpose: Determine the orientation value of the vertices in a simplex.
278 *
279 * Author: Michael Holst
280 * ***************************************************************************
281 */
Gem_orientVal(Gem * thee,SS * sm)282 VPUBLIC double Gem_orientVal(Gem *thee, SS *sm)
283 {
284 int j;
285 double a[3], b[3], c[3], d[3];
286
287 Gem_predinit(thee);
288
289 if (Gem_dim(thee)==2) {
290 for (j=0; j<2; j++) {
291 a[j] = VV_coord(SS_vertex(sm,0),j);
292 b[j] = VV_coord(SS_vertex(sm,1),j);
293 c[j] = VV_coord(SS_vertex(sm,2),j);
294 }
295 return Vpred_orient2d(a,b,c);
296 } else if (Gem_dim(thee)==3) {
297 for (j=0; j<3; j++) {
298 a[j] = VV_coord(SS_vertex(sm,1),j);
299 b[j] = VV_coord(SS_vertex(sm,0),j);
300 c[j] = VV_coord(SS_vertex(sm,2),j);
301 d[j] = VV_coord(SS_vertex(sm,3),j);
302 }
303 return Vpred_orient3d(a,b,c,d);
304 } else {
305 VASSERT(0);
306 return 0;
307 }
308 }
309
310 /*
311 * ***************************************************************************
312 * Routine: Gem_inSphere
313 *
314 * Purpose: Determine if a vertex lies in a sphere of other vertices.
315 *
316 * Author: Michael Holst
317 * ***************************************************************************
318 */
Gem_inSphere(Gem * thee,SS * sm,int sm_facet,VV * vx,VV * vxnb)319 VPUBLIC int Gem_inSphere(Gem *thee, SS *sm, int sm_facet, VV *vx, VV *vxnb)
320 {
321 int j;
322 double a[3], b[3], c[3], d[3], e[3];
323
324 Gem_predinit(thee);
325
326 if (Gem_dim(thee)==2) {
327 for (j=0; j<2; j++) {
328 a[j] = VV_coord(vx,j);
329 b[j] = VV_coord(SS_vertex(sm,SS_faceVertexNumber(sm,sm_facet,0)),j);
330 c[j] = VV_coord(SS_vertex(sm,SS_faceVertexNumber(sm,sm_facet,1)),j);
331 d[j] = VV_coord(vxnb,j);
332 }
333 if ( Vpred_incircle(a, b, c, d) >= VSMALL ) {
334 return 1;
335 } else if ( Vpred_incircle(a, b, c, d) <= -VSMALL ) {
336 return -1;
337 } else {
338 return 0;
339 }
340 } else if (Gem_dim(thee)==3) {
341 for (j=0; j<3; j++) {
342 a[j] = VV_coord(vx,j);
343 b[j] = VV_coord(SS_vertex(sm,SS_faceVertexNumber(sm,sm_facet,1)),j);
344 c[j] = VV_coord(SS_vertex(sm,SS_faceVertexNumber(sm,sm_facet,0)),j);
345 d[j] = VV_coord(SS_vertex(sm,SS_faceVertexNumber(sm,sm_facet,2)),j);
346 e[j] = VV_coord(vxnb,j);
347 }
348 if ( Vpred_insphere(a, b, c, d, e) >= VSMALL ) {
349 return 1;
350 } else if ( Vpred_insphere(a, b, c, d, e) <= -VSMALL ) {
351 return -1;
352 } else {
353 return 0;
354 }
355 } else {
356 VASSERT(0);
357 return 0;
358 }
359 }
360
361 /*
362 * ***************************************************************************
363 * Routine: Gem_delaunay
364 *
365 * Purpose: Incremental flip Delaunay generator.
366 *
367 * Notes: We use an incremental flip Delaunay mesh generator.
368 * In 2D, this is based on the standard edge-flipping.
369 * In 3D, this is based on 2-to-3 flips (shared face-to-edge flip),
370 * and 3-to-2 flips (a restricted edge-to-face flip).
371 *
372 * The 2D version of the algorithm here is similar to several
373 * edge-flip-based Delaunay algorithms in the literature.
374 *
375 * The 3D version of the algorithm here is similar to Barry Joe's
376 * and Ernst Mucke's. However, the ringed-vertex datastructure used
377 * here leads to a more concise implementation, as compared to
378 * implementations using e.g. Mucke's edge-facet datastructure.
379 *
380 * Both this algorithm and similar incremental flip algorithms in
381 * the 3D case are based on several recent theoretical results
382 * due to Barry Joe, which guarantee the following:
383 *
384 * (1) The flipping algorithm to re-establish Delaunay-ness
385 * always terminates in a finite number of steps.
386 *
387 * (2) The algorithm works regardless of the flipping order.
388 *
389 * (3) Only the exterior faces of the star region of the new
390 * vertex (what Mucke calls "link facets", because these
391 * "facets" link two triangles or tets together), and in the
392 * 3D case the three edges which make up those faces, need
393 * be tested and then possibly flipped.
394 *
395 * The algorithm is as follows:
396 *
397 * (1) Given N inputs points, a single enclosing simplex is
398 * formed by adding d+1 additional points, and then forming
399 * that single simplex.
400 *
401 * (2) The N points are then added to the mesh one at a time
402 * by locating the simplex containing each point, adding
403 * the point, and splitting the containing simplex into
404 * d+1 children (note that unless the points are in
405 * "general" position, this may give rise to degenerate
406 * simplices).
407 *
408 * (3) The link-facets of the newly added vertex are checked
409 * for Delaunayness. The link-facets are the faces
410 * opposite the new vertex in each simplex that uses
411 * the new vertex. (This is the boundary of the support
412 * region for a finite element basis function, for example.)
413 * If a non-Delaunay face is located, we attempt to flip
414 * one of the three edges of the face using a 3-to-2
415 * (edge-to-face) flip. We only flip such as edge if the
416 * simplex ring about the edge has length three. If no
417 * such edge flips are possible, we do a 2-to-3 flip
418 * (face-to-edge) if this is possible (it is only possible
419 * if the two tets sharing the face form a convex region).
420 * If no flipping is possible, we temporarily ignore this
421 * non-Delaunay face and move on to the next face.
422 *
423 * (4) Step (3) above terminates in a finite number of steps
424 * thanks to the theoretical results of Barry Joe.
425 *
426 * Author: Michael Holst
427 * ***************************************************************************
428 */
429 #define VXNUM 20
Gem_delaunay(Gem * thee)430 VPUBLIC void Gem_delaunay(Gem *thee)
431 {
432 int i, j, theDim, vxid, reality, class, type, color;
433 double max, *defX[3] = { VNULL, VNULL, VNULL };
434 VV *vx, *v[4], *vb[4];
435 SS *sm, *lastS, *fsm[4];
436 Vio *sock;
437
438 /* dimension */
439 theDim = Gem_dim(thee);
440 Vnm_print(2,"Gem_delaunay: dimension=<%d>\n", theDim);
441
442 /* reset the boundary vertex/face counters */
443 thee->numBV = 0;
444 thee->numBF = 0;
445
446 /* reset boundary information and first simplex ptrs, also get bbox */
447 max = -VLARGE;
448 for (vxid=0; vxid<Gem_numVV(thee); vxid++) {
449 vx = Gem_VV(thee,vxid);
450 /* VV_setType(vx, 0); */
451 VV_setFirstSS(vx, VNULL);
452 for (j=0; j<Gem_dimII(thee); j++)
453 if (VABS(VV_coord(vx,j)) > max) max = VABS(VV_coord(vx,j));
454 }
455
456 /* make bbox-vertices and then the bbox-simplex from the bounding box */
457 for (i=0; i<Gem_dimVV(thee); i++) {
458
459 /* create the new vertex */
460 vx = Gem_createAndInitVV(thee);
461 vb[i] = vx;
462 VV_setType(vx, 1);
463
464 /* set the vertex coordinates */
465 if (Gem_dim(thee) == 2) {
466 if (i==0) {
467 VV_setCoord(vx, 0, -10.0*max);
468 VV_setCoord(vx, 1, -10.0*max);
469 VV_setCoord(vx, 2, 0.0*max);
470 } else if (i==1) {
471 VV_setCoord(vx, 0, 10.0*max);
472 VV_setCoord(vx, 1, -10.0*max);
473 VV_setCoord(vx, 2, 0.0*max);
474 } else if (i==2) {
475 VV_setCoord(vx, 0, 0.0*max);
476 VV_setCoord(vx, 1, 10.0*max);
477 VV_setCoord(vx, 2, 0.0*max);
478 } else { VASSERT(0); }
479 } else if (Gem_dim(thee) == 3) {
480 if (i==0) {
481 VV_setCoord(vx, 0, 20.0*max);
482 VV_setCoord(vx, 1, 0.0*max);
483 VV_setCoord(vx, 2, 0.0*max);
484 } else if (i==1) {
485 VV_setCoord(vx, 0, 0.0*max);
486 VV_setCoord(vx, 1, 20.0*max);
487 VV_setCoord(vx, 2, 0.0*max);
488 } else if (i==2) {
489 VV_setCoord(vx, 0, 0.0*max);
490 VV_setCoord(vx, 1, 0.0*max);
491 VV_setCoord(vx, 2, 20.0*max);
492 } else if (i==3) {
493 VV_setCoord(vx, 0, -20.0*max);
494 VV_setCoord(vx, 1, -20.0*max);
495 VV_setCoord(vx, 2, -20.0*max);
496 } else { VASSERT(0); }
497 } else { VASSERT(0); }
498 }
499
500 /* delete any existing simplices */
501 if ( thee->simplices != VNULL ) Vset_dtor( &(thee->simplices) );
502 thee->simplices = Vset_ctor( VNULL, "SS", sizeof( SS ), VMAX_OBJECTS, 0 );
503
504 /* make the bbox simplex */
505 sm = Gem_createAndInitSS(thee);
506 SS_setDim(sm, theDim);
507 SS_setClass(sm, theDim);
508 /* set the simplex face types and vertex labels */
509 for (j=0; j<Gem_dimVV(thee); j++) {
510 SS_setFaceType( sm, j, 0 );
511 SS_setVertex( sm, j, vb[j] );
512 }
513 /* build the ringed vertex datastructure */
514 SS_buildRing(sm);
515
516 /* make sure we have the right orientation on initial guy */
517 if ( Gem_orient(thee,sm) < 0 ) SS_reverse(sm);
518 VASSERT( Gem_orient(thee,sm) > 0 );
519
520 /* make a picture */
521 sock = Vio_socketOpen("w","INET", "ASC", "localhost", "1");
522 Gem_writeGV(thee,sock,0,0,0,1.0,0,defX);
523 Vio_socketClose(&sock);
524 #if 0
525 Vnm_print(2,"HIT RETURN TO CONTINUE...");
526 c=getchar();
527 #endif
528
529 /* insert verts (except last dimVV guys in bbox) one at a time into mesh */
530 for (vxid=0; vxid<Gem_numVV(thee)-Gem_dimVV(thee); vxid++) {
531 vx = Gem_VV(thee,vxid);
532 if (vxid < VXNUM) {
533 Vnm_print(2,"Gem_delaunay: inserting V<%d>\n", VV_id(vx) );
534 }
535
536 /* find encasing simplex, remove it, replace with d+1 children */
537 if ( (sm=Gem_findSimplex(thee,vx)) != VNULL ) {
538
539 /* initialize and get simplices' vertices and its face types */
540 for (i=0; i<4; i++)
541 fsm[i] = VNULL;
542 for (i=0; i<4; i++) {
543 v[i] = VNULL;
544 }
545 for (i=0; i<Gem_dimVV(thee); i++)
546 v[i] = SS_vertex( sm, i );
547
548 /* get various attibutes of parent that children should inherit */
549 reality = SS_reality(sm); /* we must force inheritance here */
550 class = SS_class(sm); /* we must force inheritance here */
551 type = SS_type(sm); /* we must force inheritance here */
552 color = SS_chart(sm); /* we must force inheritance here */
553
554 /* remove sm from simplex rings, reinit for reuse as first child */
555 SS_meltRing(sm);
556 SS_reinit(sm);
557
558 /* first simplex -- new vx plays the role of the old v0 */
559 fsm[0] = sm; /* first child re-uses the parent simplex */
560 SS_setReality(fsm[0], reality);
561 SS_setClass(fsm[0], class);
562 SS_setType(fsm[0], type);
563 SS_setChart(fsm[0], color);
564 SS_setVertex(fsm[0], 0, vx);
565 SS_setVertex(fsm[0], 1, v[0]);
566 SS_setVertex(fsm[0], 2, v[1]);
567 SS_setFaceType(fsm[0], 0, 0);
568 SS_setFaceType(fsm[0], 1, 0);
569 SS_setFaceType(fsm[0], 2, 0);
570 if (Gem_dim(thee)==3) {
571 SS_setVertex(fsm[0], 3, v[3]);
572 SS_setFaceType(fsm[0], 3, 0);
573 } else {
574 SS_setVertex(fsm[0], 3, VNULL);
575 SS_setFaceType(fsm[0], 3, 0);
576 }
577 SS_buildRing(fsm[0]);
578
579 /* second simplex -- new vx plays the role of the old v1 */
580 fsm[1] = Gem_createAndInitSS(thee);
581 SS_setReality(fsm[1], reality);
582 SS_setClass(fsm[1], class);
583 SS_setType(fsm[1], type);
584 SS_setChart(fsm[1], color);
585 SS_setVertex(fsm[1], 0, vx);
586 SS_setVertex(fsm[1], 1, v[1]);
587 SS_setVertex(fsm[1], 2, v[2]);
588 SS_setFaceType(fsm[1], 0, 0);
589 SS_setFaceType(fsm[1], 1, 0);
590 SS_setFaceType(fsm[1], 2, 0);
591 if (Gem_dim(thee)==3) {
592 SS_setVertex(fsm[1], 3, v[3]);
593 SS_setFaceType(fsm[1], 3, 0);
594 } else {
595 SS_setVertex(fsm[1], 3, VNULL);
596 SS_setFaceType(fsm[1], 3, 0);
597 }
598 SS_buildRing(fsm[1]);
599
600 /* third simplex -- new vx plays the role of the old v2 */
601 fsm[2] = Gem_createAndInitSS(thee);
602 SS_setReality(fsm[2], reality);
603 SS_setClass(fsm[2], class);
604 SS_setType(fsm[2], type);
605 SS_setChart(fsm[2], color);
606 SS_setVertex(fsm[2], 0, vx);
607 SS_setVertex(fsm[2], 1, v[2]);
608 SS_setVertex(fsm[2], 2, v[0]);
609 SS_setFaceType(fsm[2], 0, 0);
610 SS_setFaceType(fsm[2], 1, 0);
611 SS_setFaceType(fsm[2], 2, 0);
612 if (Gem_dim(thee)==3) {
613 SS_setVertex(fsm[2], 3, v[3]);
614 SS_setFaceType(fsm[2], 3, 0);
615 } else {
616 SS_setVertex(fsm[2], 3, VNULL);
617 SS_setFaceType(fsm[2], 3, 0);
618 }
619 SS_buildRing(fsm[2]);
620
621 /* fourth simplex -- new vx plays the role of the old v3 */
622 if (Gem_dim(thee)==3) {
623 fsm[3] = Gem_createAndInitSS(thee);
624 SS_setReality(fsm[3], reality);
625 SS_setClass(fsm[3], class);
626 SS_setType(fsm[3], type);
627 SS_setChart(fsm[3], color);
628 SS_setVertex(fsm[3], 0, v[0]);
629 SS_setVertex(fsm[3], 1, v[1]);
630 SS_setVertex(fsm[3], 2, v[2]);
631 SS_setVertex(fsm[3], 3, vx);
632 SS_setFaceType(fsm[3], 0, 0);
633 SS_setFaceType(fsm[3], 1, 0);
634 SS_setFaceType(fsm[3], 2, 0);
635 SS_setFaceType(fsm[3], 3, 0);
636 SS_buildRing(fsm[3]);
637 }
638
639 } else { VASSERT(0); }
640
641 /* make a picture */
642 if (vxid < VXNUM) {
643 #if 0
644 sock = Vio_socketOpen("w","INET", "ASC", "localhost", "1");
645 Gem_writeGV(thee,sock,0,0,0,1.0,0,defX);
646 Vio_socketClose(&sock);
647 Vnm_print(2,"HIT RETURN TO CONTINUE...");
648 c=getchar();
649 #endif
650 }
651
652 /* flip link facets of this point until we are delaunay */
653 Gem_flip(thee, vx);
654
655 /* make a picture */
656 /* if ( (vxid < VXNUM) && ((vxid % 10) == 0) ){ */
657 if (vxid < VXNUM) {
658 sock = Vio_socketOpen("w","INET", "ASC", "localhost", "1");
659 Gem_writeGV(thee,sock,0,0,0,1.0,0,defX);
660 Vio_socketClose(&sock);
661 #if 0
662 Vnm_print(2,"HIT RETURN TO CONTINUE...");
663 c=getchar();
664 #endif
665 }
666 }
667
668 /* delete the entire rings of simplices which use the bbox vertices */
669 for (i=0; i<Gem_dimVV(thee); i++) {
670
671 /* delete the current ring */
672 while ( (sm=VV_firstSS(vb[i])) != VNULL ) {
673
674 /* pull sm out of its simplex rings */
675 SS_meltRing(sm);
676 SS_reinit(sm);
677
678 /* if sm is not last in the list, exchange it with the last one */
679 if ( sm != (lastS=Gem_lastSS(thee)) ) {
680 SS_meltRing(lastS);
681 SS_setType(sm, SS_type(lastS));
682 SS_setChart(sm, SS_chart(lastS));
683 for (j=0; j<Gem_dimVV(thee); j++) {
684 SS_setVertex( sm, j, SS_vertex(lastS,j) );
685 SS_setFaceType( sm, j, SS_faceType(lastS,j) );
686 }
687 SS_reinit(lastS);
688 SS_buildRing(sm);
689 }
690
691 /* now remove the last guy */
692 Gem_destroySS(thee);
693 }
694 }
695
696 /* delete the bbox vertices */
697 for (i=0; i<Gem_dimVV(thee); i++)
698 Gem_destroyVV(thee);
699
700 /* one pass through simplices to destroy non-convex simplices */
701 #if 0
702 smid = 0;
703 while (smid < Gem_numSS(thee)) {
704 sm = Gem_SS(thee,smid);
705 found = 0;
706 for (j=0; j<Gem_dimVV(thee); j++)
707 if (VV_type(SS_vertex(sm,j)) == 0)
708 found = 1;
709 if (!found) {
710 /* pull sm out of its simplex rings */
711 SS_meltRing(sm);
712 SS_reinit(sm);
713
714 /* if sm is not last in the list, exchange it with the last one */
715 if ( sm != (lastS=Gem_lastSS(thee)) ) {
716 SS_meltRing(lastS);
717 SS_setType(sm, SS_type(lastS));
718 SS_setChart(sm, SS_chart(lastS));
719 for (j=0; j<Gem_dimVV(thee); j++) {
720 SS_setVertex( sm, j, SS_vertex(lastS,j) );
721 SS_setFaceType( sm, j, SS_faceType(lastS,j) );
722 }
723 SS_reinit(lastS);
724 SS_buildRing(sm);
725 } else {
726 smid++;
727 }
728
729 /* now remove the last guy */
730 Gem_destroySS(thee);
731 } else {
732 smid++;
733 }
734 }
735 #endif
736
737 /* make a boundary */
738 Gem_makeBnd(thee,1);
739
740 /* make a picture */
741 Vnm_print(2,"DRAWING FINAL MESH\n");
742 sock = Vio_socketOpen("w","INET", "ASC", "localhost", "1");
743 Gem_writeGV(thee,sock,0,0,0,1.0,0,defX);
744 Vio_socketClose(&sock);
745 }
746
747 /*
748 * ***************************************************************************
749 * Routine: Gem_flip
750 *
751 * Purpose: Edge or face flip for the incremental flip algorithm.
752 *
753 * Author: Michael Holst
754 * ***************************************************************************
755 */
Gem_flip(Gem * thee,VV * vx)756 VPUBLIC void Gem_flip(Gem *thee, VV *vx)
757 {
758 int i, j, k, currentQ, qid, sm0_facet, sm1_facet, sm2_facet, smnb_facet;
759 int reality, class, type, color;
760 int scnt, found, jmax, imin, imax, swapped, smid;
761 SS *sm0, *sm1, *sm2, *sm3, *smnb, *ss[100], *ss2[100], *lastS, *sm;
762 SS sx0, sx1, sx2;
763 VV *vx0, *vx1, *vx2, *vxnb, *vv[100];
764
765 /* check the simplex Qs for emptyness */
766 currentQ = 0;
767 if (Gem_numSQ(thee,currentQ) > 0) {
768 Vnm_print(2,"Gem_flip: non-empty flip Q%d....clearing..",currentQ);
769 Gem_resetSQ(thee,currentQ);
770 Vnm_print(2,"..done.\n");
771 }
772 VASSERT( Gem_numSQ(thee,currentQ) == 0 );
773 if (Gem_numSQ(thee,!currentQ) > 0) {
774 Vnm_print(2,"Gem_flip: non-empty flip Q%d....clearing..",!currentQ);
775 Gem_resetSQ(thee,!currentQ);
776 Vnm_print(2,"..done.\n");
777 }
778 VASSERT( Gem_numSQ(thee,!currentQ) == 0 );
779
780 /* place all non-Delaunay link facets into the initial Q */
781 for (sm0=VV_firstSS(vx); sm0!=VNULL; sm0=SS_link(sm0,vx)) {
782 sm0_facet = SS_vptr2localVnum(sm0,vx);
783 sm1 = SS_nabor(sm0,sm0_facet);
784 if (sm1 != VNULL) {
785 sm1_facet = SS_sharedFaceLocNum(sm1,sm0);
786 vxnb = SS_vertex(sm1,sm1_facet);
787 if ( Gem_inSphere(thee,sm0,sm0_facet,vx,vxnb) > 0 )
788 Gem_appendSQ(thee, currentQ, sm0);
789 }
790 }
791
792 /* traverse current Q and flip non-Delaunay link-facets until Q is empty */
793 while ( Gem_numSQ(thee,currentQ) > 0 ) {
794
795 /* traverse Q, flipping non-Delaunay link-facets */
796 for (qid=0; qid<Gem_numSQ(thee,currentQ); qid++) {
797
798 /* retrieve the Q link-facet pair */
799 sm0 = Gem_SQ(thee,currentQ, qid);
800 sm0_facet = SS_vptr2localVnum(sm0,vx);
801 /* VASSERT( Gem_orient(thee,sm0) >= 0 ); */
802
803 /* if non-null nabor */
804 sm1 = SS_nabor(sm0,sm0_facet);
805 if (sm1 != VNULL) {
806 sm1_facet = SS_sharedFaceLocNum(sm1,sm0);
807 vxnb = SS_vertex(sm1,sm1_facet);
808
809 /* if still non-delaunay */
810 if ( Gem_inSphere(thee,sm0,sm0_facet,vx,vxnb) > 0 ) {
811
812 /* 2D: just flip the edge (if region is convex) */
813 if (Gem_dim(thee) == 2) {
814
815 /* 2-2 swap */
816 for (i=0; i<3; i++) {
817 SS_setVertex(&sx0, i, SS_vertex(sm0,i));
818 SS_setVertex(&sx1, i, SS_vertex(sm1,i));
819 }
820 SS_setVertex(&sx0, vmapF[sm0_facet][0], vxnb);
821 SS_setVertex(&sx1, vmapF[sm1_facet][0], vx);
822
823 /* if this region was non-convex, swap back and re-Q sm0 */
824 if ( (Gem_orient(thee,&sx0) < 0)
825 || (Gem_orient(thee,&sx1) < 0) ) {
826
827 smnb = SS_nabor(sm0,sm0_facet);
828 if ( smnb != VNULL ) {
829 smnb_facet = SS_sharedFaceLocNum(smnb,sm0);
830 vxnb = SS_vertex(smnb,smnb_facet);
831 /* if facet non-delaunay, place sm0 on second Q */
832 if ( Gem_inSphere(thee,sm0,sm0_facet,vx,vxnb) > 0 )
833 Gem_appendSQ(thee, !currentQ, sm0);
834 }
835
836 /* region convex; check two new faces for delaunayness */
837 } else {
838 SS_meltRing(sm0);
839 SS_meltRing(sm1);
840 SS_setVertex(sm0, vmapF[sm0_facet][0], vxnb);
841 SS_setVertex(sm1, vmapF[sm1_facet][0], vx);
842 SS_buildRing(sm0);
843 SS_buildRing(sm1);
844 VASSERT( Gem_orient(thee,sm0) >= 0 );
845 VASSERT( Gem_orient(thee,sm1) >= 0 );
846
847 /* check for delaunayness of sm0 */
848 sm0_facet = SS_vptr2localVnum(sm0,vx);
849 smnb = SS_nabor(sm0,sm0_facet);
850 if ( smnb != VNULL ) {
851 smnb_facet = SS_sharedFaceLocNum(smnb,sm0);
852 vxnb = SS_vertex(smnb,smnb_facet);
853 /* if facet non-delaunay, place sm0 on second Q */
854 if ( Gem_inSphere(thee,sm0,sm0_facet,vx,vxnb) > 0 )
855 Gem_appendSQ(thee, !currentQ, sm0);
856 }
857
858 /* check for delaunayness of sm1 */
859 sm1_facet = SS_vptr2localVnum(sm1,vx);
860 smnb = SS_nabor(sm1,sm1_facet);
861 if ( smnb != VNULL ) {
862 smnb_facet = SS_sharedFaceLocNum(smnb,sm1);
863 vxnb = SS_vertex(smnb,smnb_facet);
864 /* if facet non-delaunay, place sm1 on second Q */
865 if ( Gem_inSphere(thee,sm1,sm1_facet,vx,vxnb) > 0 )
866 Gem_appendSQ(thee, !currentQ, sm1);
867 }
868 }
869
870 /* 3D: try 3-ring edge flip first, then faces (if convex) */
871 } else if (Gem_dim(thee) == 3) {
872
873 vx0 = SS_vertex(sm0, vmapF[sm0_facet][0]);
874 vx1 = SS_vertex(sm0, vmapF[sm0_facet][1]);
875 vx2 = SS_vertex(sm0, vmapF[sm0_facet][2]);
876
877 /* 3-2 (edge-to-face) swap */
878 swapped = 0;
879 scnt = 0;
880
881 /* edge (0-1) attempt */
882 if (scnt != 3) {
883 scnt = 0;
884 vv[0] = vx0;
885 vv[1] = vx1;
886 for (sm2=VV_firstSS(vv[0]); sm2!=VNULL;
887 sm2=SS_link(sm2,vv[0])) {
888 for (sm3=VV_firstSS(vv[1]); sm3!=VNULL;
889 sm3=SS_link(sm3,vv[1])) {
890 if (sm2 == sm3) {
891 ss[scnt] = sm2;
892 scnt++;
893 }
894 }
895 }
896 VASSERT( scnt <= 100 );
897 }
898
899 /* edge (1-2) attempt */
900 if (scnt != 3) {
901 scnt = 0;
902 vv[0] = vx1;
903 vv[1] = vx2;
904 for (sm2=VV_firstSS(vv[0]); sm2!=VNULL;
905 sm2=SS_link(sm2,vv[0])) {
906 for (sm3=VV_firstSS(vv[1]); sm3!=VNULL;
907 sm3=SS_link(sm3,vv[1])) {
908 if (sm2 == sm3) {
909 ss[scnt] = sm2;
910 scnt++;
911 }
912 }
913 }
914 VASSERT( scnt <= 100 );
915 }
916
917 /* edge (2-0) attempt */
918 if (scnt != 3) {
919 scnt = 0;
920 vv[0] = vx2;
921 vv[1] = vx0;
922 for (sm2=VV_firstSS(vv[0]); sm2!=VNULL;
923 sm2=SS_link(sm2,vv[0])) {
924 for (sm3=VV_firstSS(vv[1]); sm3!=VNULL;
925 sm3=SS_link(sm3,vv[1])) {
926 if (sm2 == sm3) {
927 ss[scnt] = sm2;
928 scnt++;
929 }
930 }
931 }
932 VASSERT( scnt <= 100 );
933 }
934
935 /* can we do one of the 3-2 (edge-to-face) swaps */
936 if (scnt == 3) {
937
938 Vnm_print(2,"STARTING: 3-2 (edge-to-face) flip\n");
939
940 /* grab/count the unique vertices in this situation */
941 vv[2] = vx;
942 vv[3] = vxnb;
943 jmax = 4;
944 for (i=0; i<3; i++) {
945 sm2 = ss[i];
946 for (k=0; k<4; k++) {
947 found = 0;
948 for (j=0; j<jmax; j++)
949 if (SS_vertex(sm2, k) == vv[j]) found = 1;
950 if (!found) {
951 vv[jmax] = SS_vertex(sm2, k);
952 jmax++;
953 }
954 }
955 }
956
957 /* if jmax==5, then 3-simplex ring closed; can flip */
958 if (jmax == 5) {
959
960 /* order simps touching vx first in simplex list */
961 imin = 0;
962 imax = 2;
963 for (i=0; i<3; i++) {
964 found = 0;
965 for (j=0; j<4; j++)
966 if (SS_vertex(ss[i],j) == vx) found = 1;
967 if (found) {
968 ss2[imin] = ss[i];
969 imin++;
970 } else {
971 ss2[imax] = ss[i];
972 imax--;
973 }
974 }
975 VASSERT( (1 <= imin) && (imin <= 2) );
976 VASSERT( imin == (imax+1) );
977 sm0 = ss2[0];
978 sm1 = ss2[1];
979 sm2 = ss2[2];
980
981 /* sm2 not touch vx; pull out of simplex rings */
982 SS_meltRing(sm2);
983 SS_reinit(sm2);
984
985 Vnm_print(2,"FINISHING: 3-2 (edge-to-face) flip\n");
986 SS_meltRing(sm0);
987 SS_meltRing(sm1);
988 SS_setVertex(sm0, 0, vv[0]);
989 SS_setVertex(sm0, 1, vv[2]);
990 SS_setVertex(sm0, 2, vv[3]);
991 SS_setVertex(sm0, 3, vv[4]);
992 SS_setVertex(sm1, 0, vv[1]);
993 SS_setVertex(sm1, 1, vv[2]);
994 SS_setVertex(sm1, 2, vv[3]);
995 SS_setVertex(sm1, 3, vv[4]);
996 SS_buildRing(sm0);
997 SS_buildRing(sm1);
998 if ( Gem_orient(thee,sm0) < 0 ) SS_reverse(sm0);
999 if ( Gem_orient(thee,sm1) < 0 ) SS_reverse(sm1);
1000 VASSERT( Gem_orient(thee,sm0) > 0 );
1001 VASSERT( Gem_orient(thee,sm1) > 0 );
1002
1003 swapped = 1;
1004 } else {
1005 Vnm_print(2,"BACKING-OUT: 3-to-2 flip\n");
1006 }
1007 }
1008
1009 /* 2-3 (face-to-edge) swap */
1010 if (!swapped) {
1011
1012 Vnm_print(2,"STARTING: 2-3 (face-to-edge) flip\n");
1013
1014 /* attempt the face-to-edge swap */
1015 SS_setVertex(&sx0, 0, vx);
1016 SS_setVertex(&sx0, 1, vx0);
1017 SS_setVertex(&sx0, 2, vx1);
1018 SS_setVertex(&sx0, 3, vxnb);
1019 SS_setVertex(&sx1, 0, vx);
1020 SS_setVertex(&sx1, 1, vx1);
1021 SS_setVertex(&sx1, 2, vx2);
1022 SS_setVertex(&sx1, 3, vxnb);
1023 SS_setVertex(&sx2, 0, vx);
1024 SS_setVertex(&sx2, 1, vx2);
1025 SS_setVertex(&sx2, 2, vx0);
1026 SS_setVertex(&sx2, 3, vxnb);
1027 if ( Gem_orient(thee,&sx0) < 0 ) {
1028 Vnm_print(2,"-- sx0 BAD: <%g> --\n",
1029 Gem_orientVal(thee,&sx0));
1030 }
1031 if ( Gem_orient(thee,&sx1) < 0 ) {
1032 Vnm_print(2,"-- sx1 BAD: <%g> --\n",
1033 Gem_orientVal(thee,&sx1));
1034 }
1035 if ( Gem_orient(thee,&sx2) < 0 ) {
1036 Vnm_print(2,"-- sx2 BAD: <%g> --\n",
1037 Gem_orientVal(thee,&sx2));
1038 }
1039
1040 /* if region was non-convex, swap back and re-Q sm0 */
1041 if ( (Gem_orient(thee,&sx0) < 0)
1042 || (Gem_orient(thee,&sx1) < 0)
1043 || (Gem_orient(thee,&sx2) < 0) ) {
1044 Vnm_print(2,"BACKING-OUT: 2-to-3 flip\n");
1045 smnb = SS_nabor(sm0,sm0_facet);
1046 if ( smnb != VNULL ) {
1047 smnb_facet = SS_sharedFaceLocNum(smnb,sm0);
1048 vxnb = SS_vertex(smnb,smnb_facet);
1049 /* if facet non-delaunay, place sm0 on Q2 */
1050 if (Gem_inSphere(thee,sm0,sm0_facet,vx,vxnb)>0)
1051 Gem_appendSQ(thee, !currentQ, sm0);
1052 }
1053
1054 /* region convex; check 3 new faces for delaunayness */
1055 } else {
1056 Vnm_print(2,"FINISHING: 2-3 (face-to-edge) flip\n");
1057
1058 /* get various attibutes of parent */
1059 reality = SS_reality(sm0);
1060 class = SS_class(sm0);
1061 type = SS_type(sm0);
1062 color = SS_chart(sm0);
1063
1064 /* grab a third simplex */
1065 sm2 = Gem_createAndInitSS(thee);
1066 SS_setReality(sm2, reality);
1067 SS_setClass(sm2, class);
1068 SS_setType(sm2, type);
1069 SS_setChart(sm2, color);
1070 SS_setFaceType(sm2, 0, 0);
1071 SS_setFaceType(sm2, 1, 0);
1072 SS_setFaceType(sm2, 2, 0);
1073 SS_setFaceType(sm2, 3, 0);
1074
1075 /* set things up */
1076 SS_meltRing(sm0);
1077 SS_meltRing(sm1);
1078 Vnm_print(2,"vx=<%d> v0=<%d> v1=<%d> v2=<%d> vxnb=<%d>\n",
1079 VV_id(vx), VV_id(vx0), VV_id(vx1), VV_id(vx2), VV_id(vxnb));
1080 SS_setVertex(sm0, 0, vx);
1081 SS_setVertex(sm0, 1, vx0);
1082 SS_setVertex(sm0, 2, vx1);
1083 SS_setVertex(sm0, 3, vxnb);
1084 SS_setVertex(sm1, 0, vx);
1085 SS_setVertex(sm1, 1, vx1);
1086 SS_setVertex(sm1, 2, vx2);
1087 SS_setVertex(sm1, 3, vxnb);
1088 SS_setVertex(sm2, 0, vx);
1089 SS_setVertex(sm2, 1, vx2);
1090 SS_setVertex(sm2, 2, vx0);
1091 SS_setVertex(sm2, 3, vxnb);
1092 SS_buildRing(sm0);
1093 SS_buildRing(sm1);
1094 SS_buildRing(sm2);
1095 VASSERT( Gem_orient(thee,sm0) >= 0 );
1096 VASSERT( Gem_orient(thee,sm1) >= 0 );
1097 VASSERT( Gem_orient(thee,sm2) >= 0 );
1098
1099 /* check for delaunayness of sm0 */
1100 sm0_facet = SS_vptr2localVnum(sm0,vx);
1101 smnb = SS_nabor(sm0,sm0_facet);
1102 if ( smnb != VNULL ) {
1103 smnb_facet = SS_sharedFaceLocNum(smnb,sm0);
1104 vxnb = SS_vertex(smnb,smnb_facet);
1105 /* if facet non-delaunay, place sm0 on Q2 */
1106 if (Gem_inSphere(thee,sm0,sm0_facet,vx,vxnb)>0)
1107 Gem_appendSQ(thee, !currentQ, sm0);
1108 }
1109
1110 /* check for delaunayness of sm1 */
1111 sm1_facet = SS_vptr2localVnum(sm1,vx);
1112 smnb = SS_nabor(sm1,sm1_facet);
1113 if ( smnb != VNULL ) {
1114 smnb_facet = SS_sharedFaceLocNum(smnb,sm1);
1115 vxnb = SS_vertex(smnb,smnb_facet);
1116 /* if facet non-delaunay, place sm1 on Q2 */
1117 if (Gem_inSphere(thee,sm1,sm1_facet,vx,vxnb)>0)
1118 Gem_appendSQ(thee, !currentQ, sm1);
1119 }
1120
1121 /* check for delaunayness of sm2 */
1122 sm2_facet = SS_vptr2localVnum(sm2,vx);
1123 smnb = SS_nabor(sm2,sm2_facet);
1124 if ( smnb != VNULL ) {
1125 smnb_facet = SS_sharedFaceLocNum(smnb,sm2);
1126 vxnb = SS_vertex(smnb,smnb_facet);
1127 /* if facet non-delaunay, place sm2 on Q2 */
1128 if (Gem_inSphere(thee,sm2,sm2_facet,vx,vxnb)>0)
1129 Gem_appendSQ(thee, !currentQ, sm2);
1130 }
1131 }
1132 }
1133
1134 } else { VASSERT(0); } /* if (Gem_dim(thee) == 2) { */
1135 } /* if ( Gem_inSphere(thee,sm0,sm0_facet,vx,vxnb) > 0 ) { */
1136 } /* if (sm1 != VNULL) { */
1137 } /* for (qid=0; qid<Gem_numSQ(thee,currentQ); qid++) { */
1138
1139 /* reset the current Q and then swap the Qs */
1140 Gem_resetSQ(thee,currentQ);
1141 currentQ = !currentQ;
1142
1143 } /* while ( Gem_numSQ(thee,currentQ) > 0 ) { */
1144
1145 /* all Qs should be empty at this point */
1146 for (i=0; i<VMAXSQ; i++)
1147 VASSERT( Gem_numSQ(thee,i) == 0 );
1148
1149 /* clean up deleted simplices */
1150 smid = 0;
1151 while (smid < Gem_numSS(thee)) {
1152 sm = Gem_SS(thee,smid);
1153 if (SS_vertex(sm,0) == VNULL) {
1154
1155 /* if sm is not last in the list, exchange it with the last one */
1156 lastS=Gem_lastSS(thee);
1157 if ( (sm != lastS) && (SS_vertex(lastS,0) != VNULL) ) {
1158 SS_meltRing(lastS);
1159 SS_setType(sm, SS_type(lastS));
1160 SS_setChart(sm, SS_chart(lastS));
1161 for (j=0; j<Gem_dimVV(thee); j++) {
1162 SS_setVertex( sm, j, SS_vertex(lastS,j) );
1163 SS_setFaceType( sm, j, SS_faceType(lastS,j) );
1164 }
1165 SS_buildRing(sm);
1166 SS_reinit(lastS);
1167 } else if (SS_vertex(lastS,0) == VNULL) {
1168 /* don't increment */
1169 } else {
1170 smid++;
1171 }
1172
1173 /* now remove the last guy */
1174 Gem_destroySS(thee);
1175 } else {
1176 smid++;
1177 }
1178 }
1179
1180 for (i=0; i<Gem_numSS(thee); i++) {
1181 sm = Gem_SS(thee,i);
1182 VASSERT( SS_vertex(sm,0) != VNULL );
1183 }
1184 }
1185
1186 /*
1187 * ***************************************************************************
1188 * Routine: Gem_deform
1189 *
1190 * Purpose: Edge or face flip for the incremental flip algorithm.
1191 *
1192 * Author: Michael Holst
1193 * ***************************************************************************
1194 */
Gem_deform(Gem * thee,int dimX,double * defX[MAXV])1195 VPUBLIC int Gem_deform(Gem *thee, int dimX, double *defX[MAXV])
1196 {
1197 int i, j;
1198 double val;
1199 VV *vx;
1200
1201 /* sanity check */
1202 VWARN( dimX == Gem_dimII(thee) );
1203
1204 /* deform (actually displace) the vertices */
1205 for (i=0; i<Gem_numVV(thee); i++) {
1206 vx = Gem_VV(thee,i);
1207 for (j=0; j<dimX; j++) {
1208 val = VV_coord(vx,j) + defX[j][i];
1209 VV_setCoord(vx, j, val);
1210 }
1211 }
1212
1213 return 0;
1214 }
1215
1216