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: gemchk.c,v 1.17 2010/08/12 05:18:53 fetk Exp $"
21 * ***************************************************************************
22 */
23
24 /*
25 * ***************************************************************************
26 * File: Gemchk.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: gemchk.c,v 1.17 2010/08/12 05:18:53 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_memChk
55 *
56 * Purpose: Print the exact current malloc usage.
57 *
58 * Author: Michael Holst
59 * ***************************************************************************
60 */
Gem_memChk(Gem * thee)61 VPUBLIC void Gem_memChk(Gem *thee)
62 {
63 int i;
64
65 if (thee->iMadeVmem) Vmem_print(thee->vmem);
66 if (thee->vertices != VNULL) Vset_memChk(thee->vertices);
67 if (thee->edges != VNULL) Vset_memChk(thee->edges);
68 if (thee->simplices != VNULL) Vset_memChk(thee->simplices);
69 for (i=0; i<VMAXSQ; i++) {
70 if (thee->sQueM[i] != VNULL) Vset_memChk(thee->sQueM[i]);
71 }
72 }
73
74 /*
75 * ***************************************************************************
76 * Routine: Gem_memChkVV
77 *
78 * Purpose: Print the exact current malloc usage: vertices.
79 *
80 * Author: Michael Holst
81 * ***************************************************************************
82 */
Gem_memChkVV(Gem * thee,int * num,int * size,int * vecUse,int * vecMal,int * vecOhd)83 VPUBLIC void Gem_memChkVV(Gem *thee, int *num,
84 int *size, int *vecUse, int *vecMal, int *vecOhd)
85 {
86 Vset_check(thee->vertices, num, size, vecUse, vecMal, vecOhd);
87 }
88
89 /*
90 * ***************************************************************************
91 * Routine: Gem_memChkEE
92 *
93 * Purpose: Print the exact current malloc usage: edges.
94 *
95 * Author: Michael Holst
96 * ***************************************************************************
97 */
Gem_memChkEE(Gem * thee,int * num,int * size,int * vecUse,int * vecMal,int * vecOhd)98 VPUBLIC void Gem_memChkEE(Gem *thee, int *num,
99 int *size, int *vecUse, int *vecMal, int *vecOhd)
100 {
101 Vset_check(thee->edges, num, size, vecUse, vecMal, vecOhd);
102 }
103
104 /*
105 * ***************************************************************************
106 * Routine: Gem_memChkSS
107 *
108 * Purpose: Print the exact current malloc usage: simplices.
109 *
110 * Author: Michael Holst
111 * ***************************************************************************
112 */
Gem_memChkSS(Gem * thee,int * num,int * size,int * vecUse,int * vecMal,int * vecOhd)113 VPUBLIC void Gem_memChkSS(Gem *thee, int *num,
114 int *size, int *vecUse, int *vecMal, int *vecOhd)
115 {
116 Vset_check(thee->simplices, num, size, vecUse, vecMal, vecOhd);
117 }
118
119 /*
120 * ***************************************************************************
121 * Routine: Gem_memChkSQ
122 *
123 * Purpose: Print the exact current malloc usage: simplex queues.
124 *
125 * Author: Michael Holst
126 * ***************************************************************************
127 */
Gem_memChkSQ(Gem * thee,int currentQ,int * tnum,int * tsize,int * tVecUse,int * tVecMal,int * tVecOhd)128 VPUBLIC void Gem_memChkSQ(Gem *thee, int currentQ,
129 int *tnum, int *tsize, int *tVecUse, int *tVecMal, int *tVecOhd)
130 {
131 Vset_check(thee->sQueM[currentQ], tnum, tsize, tVecUse, tVecMal, tVecOhd);
132 }
133
134 /*
135 * ***************************************************************************
136 * Routine: Gem_memChkMore
137 *
138 * Purpose: Estimate the current RAM usage.
139 *
140 * Author: Michael Holst
141 * ***************************************************************************
142 */
Gem_memChkMore(Gem * thee)143 VPUBLIC void Gem_memChkMore(Gem *thee)
144 {
145 int i, geom_ramUse, geom_ramMal, geom_ramOhd;
146
147 int v_num, v_size, v_ramUse, v_ramMal, v_ramOhd;
148 int e_num, e_size, e_ramUse, e_ramMal, e_ramOhd;
149 int s_num, s_size, s_ramUse, s_ramMal, s_ramOhd;
150 int sq_num, sq_size, sq_ramUse, sq_ramMal, sq_ramOhd;
151
152 Vnm_print(0, "Gem_memChkMore: Geometry memory usage\n");
153
154 /* get sizes of vertices, edges, simplices */
155 Gem_memChkVV(thee,
156 &v_num, &v_size, &v_ramUse, &v_ramMal, &v_ramOhd);
157 Gem_memChkEE(thee,
158 &e_num, &e_size, &e_ramUse, &e_ramMal, &e_ramOhd);
159 Gem_memChkSS(thee,
160 &s_num, &s_size, &s_ramUse, &s_ramMal, &s_ramOhd);
161 Vnm_print(0, "%9d VV ;%3d bytes/VV =%9d bytes[%9d:t=%9d:b+%6d:o]\n",
162 v_num, v_size, v_ramUse, (v_ramMal+v_ramOhd), v_ramMal, v_ramOhd);
163 Vnm_print(0, "%9d EE ;%3d bytes/EE =%9d bytes[%9d:t=%9d:b+%6d:o]\n",
164 e_num, e_size, e_ramUse, (e_ramMal+e_ramOhd), e_ramMal, e_ramOhd);
165 Vnm_print(0, "%9d SS ;%3d bytes/SS =%9d bytes[%9d:t=%9d:b+%6d:o]\n",
166 s_num, s_size, s_ramUse, (s_ramMal+s_ramOhd), s_ramMal, s_ramOhd);
167
168 /* get simplex Q sizes and also total all of the sizes */
169 geom_ramUse = v_ramUse + s_ramUse + e_ramUse;
170 geom_ramMal = v_ramMal + s_ramMal + e_ramMal;
171 geom_ramOhd = v_ramOhd + s_ramOhd + e_ramOhd;
172 for (i=0; i<VMAXSQ; i++) {
173 Gem_memChkSQ(thee,
174 i,&sq_num,&sq_size,&sq_ramUse,&sq_ramMal,&sq_ramOhd);
175 Vnm_print(0, "%9d Q%d ;%3d bytes/Q0 =%9d bytes[%9d:t=%9d:b+%6d:o]\n",
176 sq_num, i, sq_size, sq_ramUse, (sq_ramMal+sq_ramOhd),
177 sq_ramMal, sq_ramOhd);
178 geom_ramUse += sq_ramUse;
179 geom_ramMal += sq_ramMal;
180 geom_ramOhd += sq_ramOhd;
181 }
182 Vnm_print(0, " TOTALS =%9d bytes[%9d:t=%9d:b+%6d:o]\n",
183 geom_ramUse, (geom_ramMal+geom_ramOhd), geom_ramMal, geom_ramOhd);
184 }
185
186 /*
187 * ***************************************************************************
188 * Routine: Gem_speedChk
189 *
190 * Purpose: Calculate the cost to traverse the various structures.
191 *
192 * Author: Michael Holst
193 * ***************************************************************************
194 */
Gem_speedChk(Gem * thee)195 VPUBLIC void Gem_speedChk(Gem *thee)
196 {
197 const int N = 100000;
198 const int loop = 1000;
199 int i, j, itotal;
200 double rtotal;
201 double *d0;
202 int *i0;
203 VV *v0, *vx;
204
205 Vnm_print(0,"Gem_speedChk: Creating normal double vec...\n");
206 Vnm_tstart(90, "double vec creation");
207 d0 = Vmem_malloc( thee->vmem, N, sizeof(double) );
208 Vnm_tstop(90, "double vec creation");
209
210 Vnm_print(0,"Gem_speedChk: Creating normal int vec...\n");
211 Vnm_tstart(90, "int vec creation");
212 i0 = Vmem_malloc( thee->vmem, N, sizeof(int) );
213 Vnm_tstop(90, "int vec creation");
214
215 Vnm_print(0,"Gem_speedChk: Creating normal VV vec...\n");
216 Vnm_tstart(90, "VV vec creation");
217 v0 = Vmem_malloc( thee->vmem, N, sizeof(VV) );
218 Vnm_tstop(90, "VV vec creation");
219
220 Vnm_print(0,"Gem_speedChk: Creating Vset VV vec...\n");
221 Vnm_tstart(90, "VV Vset creation");
222 for (i=0; i<N; i++)
223 vx = Gem_createAndInitVV(thee);
224 Vnm_tstop(90, "VV Vset creation");
225
226 Vnm_print(0,"Gem_speedChk: Looping over double vec...\n");
227 Vnm_print(0,"Gem_speedChk: N, loop = %d, %d\n", N, loop);
228 Vnm_tstart(90, "double vec loop");
229 rtotal = 0.0;
230 for (j=0; j<loop; j++)
231 for (i=0; i<N; i++)
232 rtotal += ( j * d0[i] * (j+d0[i]) );
233 Vnm_tstop(90, "double vec loop");
234
235 Vnm_print(0,"Gem_speedChk: Looping over int vec...\n");
236 Vnm_print(0,"Gem_speedChk: N, loop = %d, %d\n", N, loop);
237 Vnm_tstart(90, "int vec loop");
238 itotal = 0;
239 for (j=0; j<loop; j++)
240 for (i=0; i<N; i++)
241 itotal += ( j * i0[i] * (j+i0[i]) );
242 Vnm_tstop(90, "int vec loop");
243
244 Vnm_print(0,"Gem_speedChk: Looping over VV vec...\n");
245 Vnm_print(0,"Gem_speedChk: N, loop = %d, %d\n", N, loop);
246 Vnm_tstart(90, "VV vec loop");
247 itotal = 0;
248 for (j=0; j<loop; j++)
249 for (i=0; i<N; i++)
250 itotal += ( j * VV_id(&v0[i]) * (j+VV_id(&v0[i])) );
251 Vnm_tstop(90, "VV vec loop");
252
253 Vnm_print(0,"Gem_speedChk: Looping (index) over Vset VV vec...\n");
254 Vnm_print(0,"Gem_speedChk: N, loop = %d, %d\n", N, loop);
255 Vnm_tstart(90, "VV Vset loop");
256 itotal = 0;
257 for (j=0; j<loop; j++)
258 for (i=0; i<N; i++)
259 itotal += ( j * VV_id(Gem_VV(thee,i))
260 * (j + VV_id(Gem_VV(thee,i)) ) );
261 Vnm_tstop(90, "VV Vset loop");
262
263 Vnm_print(0,"Gem_speedChk: Looping (iterator) over Vset VV vec...\n");
264 Vnm_print(0,"Gem_speedChk: N, loop = %d, %d\n", N, loop);
265 Vnm_tstart(90, "VV Vset loop");
266 itotal = 0;
267 for (j=0; j<loop; j++)
268 for (vx=Gem_firstVV(thee);vx!=VNULL;vx=Gem_nextVV(thee))
269 itotal += ( j * VV_id(vx) * (j + VV_id(vx)) );
270 Vnm_tstop(90, "VV Vset loop");
271
272 Vnm_print(0,"Gem_speedChk: Freeing normal double vec...\n");
273 Vnm_tstart(90, "Free double vec");
274 Vmem_free( thee->vmem, N, sizeof(double), (void**)&d0 );
275 Vnm_tstop(90, "Free double vec");
276
277 Vnm_print(0,"Gem_speedChk: Freeing normal int vec...\n");
278 Vnm_tstart(90, "Free int vec");
279 Vmem_free( thee->vmem, N, sizeof(int), (void**)&i0 );
280 Vnm_tstop(90, "Free int vec");
281
282 Vnm_print(0,"Gem_speedChk: Freeing normal VV vec...\n");
283 Vnm_tstart(90, "Free VV vec");
284 Vmem_free( thee->vmem, N, sizeof(VV), (void**)&v0 );
285 Vnm_tstop(90, "Free VV vec");
286
287 Vnm_print(0,"Gem_speedChk: Freeing Vset VV vec...\n");
288 Vnm_tstart(90, "Free VV Vset");
289 for (i=0; i<N; i++)
290 Gem_destroyVV(thee);
291 Vnm_tstop(90, "Free VV Vset");
292 }
293
294 /*
295 * ***************************************************************************
296 * Routine: Gem_formChk
297 *
298 * Purpose: Check the self-consistency of the geometry datastructures.
299 *
300 * Notes: key==0 --> check: min (just vertices and simplices)
301 * key==1 --> check: min + simplex ring
302 * key==2 --> check: min + simplex ring + edge ring
303 * key==3 --> check: min + simplex ring + edge ring + conform
304 *
305 * Author: Michael Holst
306 * ***************************************************************************
307 */
Gem_formChk(Gem * thee,int key)308 VPUBLIC void Gem_formChk(Gem *thee, int key)
309 {
310 int i, j, k;
311 int l, m, bndVert, bndFace;
312 int edgeOrderProb, conformCount, fType[4];
313 double svol;
314 VV *vx, *v[4];
315 SS *sm, *sm0, *sm1, *sm2;
316 EE *eg;
317
318 /* input check and some i/o */
319 VASSERT( (key >= 0) && (key <= 3) );
320
321 /* go through all vertices and check for consistency */
322 Vnm_tstart(80, "form check");
323 Vnm_print(0,"Gem_formChk: Topology check on: "
324 "<%d> SS, <%d> EE, <%d> VV:\n",
325 Gem_numSS(thee), Gem_numEE(thee), Gem_numVV(thee));
326 Vnm_print(0,"Gem_formChk: ..VV..\n");
327 bndVert = 0;
328 for (i=0; i<Gem_numVV(thee); i++) {
329 vx = Gem_VV(thee,i);
330 if ( (i>0) && (i % VPRTKEY) == 0 ) Vnm_print(0,"[CV:%d]",i);
331 VJMPERR1( !Vnm_sigInt() );
332
333 /* check basic data */
334 if ((int)VV_id(vx)!=i)
335 Vnm_print(2,"Gem_formChk: VV_id BAD vtx <%d>\n", i);
336 if (key == 0) {
337 if (VV_firstSS(vx) != VNULL)
338 Vnm_print(2,"Gem_formChk: firstS BAD vtx <%d>\n",i);
339 if (VV_firstEE(vx) != VNULL)
340 Vnm_print(2,"Gem_formChk: firstE BAD vtx <%d>\n",i);
341 } else if (key == 1) {
342 if (VV_firstSS(vx) == VNULL)
343 Vnm_print(2,"Gem_formChk: firstS BAD vtx <%d>\n",i);
344 if (VV_firstEE(vx) != VNULL)
345 Vnm_print(2,"Gem_formChk: firstE BAD vtx <%d>\n",i);
346 } else if ((key == 2) || (key == 3)) {
347 if (VV_firstSS(vx) == VNULL)
348 Vnm_print(2,"Gem_formChk: firstS BAD vtx <%d>\n",i);
349 if ((VV_firstEE(vx) == VNULL) && (Gem_numEE(thee) > 0))
350 Vnm_print(2,"Gem_formChk: firstE BAD vtx <%d>\n",i);
351 } else { VASSERT(0); }
352
353 /* boundary data */
354 if ( VBOUNDARY( VV_type(vx) ) ) bndVert++;
355 }
356 if (bndVert != Gem_numBV(thee))
357 Vnm_print(2,"Gem_formChk: bndVert BAD\n");
358
359 /* go through all edges and check for consistency */
360 if (key >= 2) {
361 Vnm_print(0,"Gem_formChk: ..EE..\n");
362 edgeOrderProb = 0;
363 for (i=0; i<Gem_numEE(thee); i++) {
364 eg = Gem_EE(thee,i);
365 if ( (i>0) && (i % VPRTKEY) == 0 ) Vnm_print(0,"[CE:%d]",i);
366 VJMPERR1( !Vnm_sigInt() );
367
368 /* check basic data */
369 /* (edge re-orderings ok; may not be error) */
370 if ((int)EE_id(eg) != i) edgeOrderProb = 1;
371 if ( EE_vertex(eg,0) == VNULL )
372 Vnm_print(2,"Gem_formChk: vertex0 BAD edge <%d>\n",i);
373 if ( EE_vertex(eg,1) == VNULL )
374 Vnm_print(2,"Gem_formChk: vertex1 BAD edge <%d>\n",i);
375
376 /* check edge ring consistencies */
377 if ( ! VV_edgeInRing(EE_vertex(eg,0), eg) )
378 Vnm_print(2,"Gem_formChk: edgeInRing0 BAD edge <%d>\n",i);
379 if ( ! VV_edgeInRing(EE_vertex(eg,1), eg) )
380 Vnm_print(2,"Gem_formChk: edgeInRing1 BAD edge <%d>\n",i);
381
382 /* check midpoint data; should be VNULL */
383 if ( EE_midPoint(eg) != VNULL )
384 Vnm_print(2,"Gem_formChk: midPoint BAD edge <%d>\n",i);
385 }
386 if (edgeOrderProb)
387 Vnm_print(0,"Gem_formChk: WARNING noncons edge order..)\n");
388 }
389
390 /* go through all simplices and check for consistency */
391 Vnm_print(0,"Gem_formChk: ..SS..\n");
392 bndFace = 0;
393 for (i=0; i<Gem_numSS(thee); i++) {
394 sm = Gem_SS(thee,i);
395 if ( (i>0) && (i % VPRTKEY) == 0 ) Vnm_print(0,"[CS:%d]",i);
396 VJMPERR1( !Vnm_sigInt() );
397
398 if ( (int)SS_id(sm) != i )
399 Vnm_print(2,"Gem_formChk: SS_id BAD sim <%d>\n",i);
400
401 /* check for well-ordering of vertices through positive volume */
402 svol = Gem_simplexVolume(thee,sm);
403 if (svol <= VSMALL)
404 Vnm_print(2,"Gem_formChk: Degenerate sim <%d> has volume <%g>\n",
405 i, svol);
406
407 /* another check for well-ordering of vertices via positive volume */
408 if ( !Gem_orient(thee,sm) )
409 Vnm_print(2,"Gem_formChk: sim <%d> is badly ordered\n", i);
410
411 /* no simplices should be marked for refinement/unrefinement */
412 if ( SS_refineKey(sm,0) != 0 )
413 Vnm_print(2,"Gem_formChk: SS_refineKey BAD sim <%d>\n",i);
414 if ( SS_refineKey(sm,1) != 0 )
415 Vnm_print(2,"Gem_formChk: SS_refineKey BAD sim <%d>\n",i);
416
417 /* get simplex face info */
418 for (j=0; j<4; j++) fType[j] = SS_faceType(sm,j);
419
420 /* check vertex data and simplex ring structure */
421 for (j=0; j<Gem_dimVV(thee); j++) {
422 v[j] = SS_vertex(sm,j);
423 if ( v[j] == VNULL )
424 Vnm_print(2,"Gem_formChk: v[%d] BAD sim <%d>\n",j,i);
425 if (key >= 1) {
426 if ( ! VV_simplexInRing( v[j], sm ) )
427 Vnm_print(2,"Gem_formChk: sInR BAD sim <%d>\n",i);
428 }
429 }
430
431 /* check simplex ring, boundary faces count, vertex type */
432 for (j=0; j<Gem_dimVV(thee); j++) {
433
434 /* boundary vertex check */
435 if ( VBOUNDARY(fType[j]) ) {
436 bndFace++;
437 for (k=1; k<Gem_dimVV(thee); k++) {
438 l=(j+k) % Gem_dimVV(thee);
439 if ( VINTERIOR( VV_type(SS_vertex(sm,l)) ) )
440 Vnm_print(2,"Gem_formChk: Dv[%d] BAD sim <%d>\n",
441 l,i);
442 }
443 }
444
445 /* face check with all nabors; verify we are conforming */
446 if (key >= 3) {
447 k=(j+1) % Gem_dimVV(thee);
448 l=(k+1) % Gem_dimVV(thee);
449 m=(l+1) % Gem_dimVV(thee);
450 conformCount = 0;
451 for (sm0=VV_firstSS(v[k]); sm0!=VNULL;sm0=SS_link(sm0,v[k])) {
452 for (sm1=VV_firstSS(v[l]);sm1!=VNULL;sm1=SS_link(sm1,v[l])){
453 if (Gem_dim(thee) == 2) {
454 if ((sm0!=sm) && (sm0==sm1)) conformCount++;
455 } else {
456 for (sm2=VV_firstSS(v[m]); sm2!=VNULL;
457 sm2=SS_link(sm2,v[m])) {
458 if ((sm0!=sm) && (sm0==sm1) && (sm0==sm2)) {
459 conformCount++;
460 }
461 }
462 }
463 }
464 }
465 if ( VBOUNDARY(fType[j]) ) {
466 if ( conformCount != 0 ) {
467 Vnm_print(2,"Gem_formChk: conform BAD sim <%d>",i);
468 Vnm_print(2,"..face <%d> should NOT have a nabor\n",j);
469 }
470 } else {
471 if ( conformCount < 1 ) {
472 Vnm_print(2,"Gem_formChk: conform BAD sim <%d>",i);
473 Vnm_print(2,"..face <%d> should have a nabor\n",j);
474 }
475 if ( conformCount > 1 ) {
476 Vnm_print(2,"Gem_formChk: conform BAD sim <%d>",i);
477 Vnm_print(2,"..face <%d> has more than one nabor\n",j);
478 }
479 }
480 }
481
482 }
483 }
484 if ( bndFace != Gem_numBF(thee) )
485 Vnm_print(2,"Gem_formChk: numBF count BAD\n");
486 Vnm_print(0,"Gem_formChk: ..done.\n");
487
488 /* return with no errors */
489 Vnm_tstop(80, "form check");
490 return;
491
492 VERROR1:
493 Vnm_print(0,"Gem_formChk: ..done. [Early termination was forced]\n");
494 Vnm_tstop(80, "form check");
495 return;
496 }
497
498 /*
499 * ***************************************************************************
500 * Routine: Gem_contentChk
501 *
502 * Purpose: Print out contents of all geometry structures.
503 *
504 * Author: Michael Holst
505 * ***************************************************************************
506 */
Gem_contentChk(Gem * thee)507 VPUBLIC void Gem_contentChk(Gem *thee)
508 {
509 int i, j, k, ii;
510 SS *sm;
511 EE *eg;
512 VV *vx;
513
514 Vnm_print(0,"Gem_contentChk: printing data.\n");
515 for (i=0; i<Gem_numSS(thee); i++) {
516 sm = Gem_SS(thee,i);
517 Vnm_print(0,"SIMPLEX [%d]\n", SS_id(sm));
518 Vnm_print(0," g.bits =<%d>\n",sm->g.bits);
519 Vnm_print(0," g.uid =<%d>\n",sm->g.uid);
520 Vnm_print(0," d.tags =<%d>\n",sm->d.tags);
521 Vnm_print(0," d.faces=<%d>\n",sm->d.faces);
522
523 Vnm_print(0," d.vPtr =");
524 for (ii=0; ii<4; ii++) {
525 j = -1;
526 if (sm->d.vPtr[ii] != VNULL) j = VV_id(sm->d.vPtr[ii]);
527 Vnm_print(0," <%d>",j);
528 }
529 Vnm_print(0,"\n");
530
531 Vnm_print(0," d.sPtr =");
532 for (ii=0; ii<4; ii++) {
533 j = -1;
534 if (sm->d.sPtr[ii] != VNULL) j = SS_id(sm->d.sPtr[ii]);
535 Vnm_print(0," <%d>",j);
536 }
537 Vnm_print(0,"\n");
538
539 # if defined(VG_ELEMENT)
540 Vnm_print(0," d.fNum =");
541 for (ii=0; ii<4; ii++) {
542 j = sm->d.fNum[ii];
543 Vnm_print(0," <%d>",j);
544 }
545 Vnm_print(0,"\n");
546
547 Vnm_print(0," d.eNum =");
548 for (ii=0; ii<6; ii++) {
549 j = sm->d.eNum[ii];
550 Vnm_print(0," <%d>",j);
551 }
552 Vnm_print(0,"\n");
553 # endif
554
555 }
556 for (i=0; i<Gem_numEE(thee); i++) {
557 eg = Gem_EE(thee,i);
558 Vnm_print(0,"EDGE [%d]\n", EE_id(eg));
559 Vnm_print(0," g.bits=<%d>\n",eg->g.bits);
560 Vnm_print(0," g.uid =<%d>\n",eg->g.uid);
561 k = 0;
562 if (eg->d.midPtr != VNULL) k = VV_id(eg->d.midPtr);
563 Vnm_print(0," d.midPtr=<%d>\n", k);
564
565 Vnm_print(0," d.vPtr =");
566 for (ii=0; ii<2; ii++) {
567 j = 0;
568 if (eg->d.vPtr[ii] != 0) j = VV_id(eg->d.vPtr[ii]);
569 Vnm_print(0," <%d>",j);
570 }
571 Vnm_print(0,"\n");
572
573 Vnm_print(0," d.ePtr =");
574 for (ii=0; ii<2; ii++) {
575 j = 0;
576 if (eg->d.ePtr[ii] != 0) j = EE_id(eg->d.ePtr[ii]);
577 Vnm_print(0," <%d>",j);
578 }
579 Vnm_print(0,"\n");
580 }
581 for (i=0; i<Gem_numVV(thee); i++) {
582 vx = Gem_VV(thee,i);
583 Vnm_print(0,"VERTEX [%d]\n", VV_id(vx));
584 Vnm_print(0," g.bits=<%d>\n",vx->g.bits);
585 Vnm_print(0," g.uid =<%d>\n",vx->g.uid);
586 j = 0;
587 k = 0;
588 if (vx->d.ePtr != VNULL) j = EE_id(vx->d.ePtr);
589 if (vx->d.sPtr != VNULL) k = SS_id(vx->d.sPtr);
590 Vnm_print(0," d.ePtr=<%d>\n", j);
591 Vnm_print(0," d.sPtr=<%d>\n", k);
592 Vnm_print(0," d.x =");
593 for (ii=0; ii<3; ii++)
594 Vnm_print(0," <%e>", vx->d.x[ii]);
595 Vnm_print(0,"\n");
596 }
597 }
598
599 /*
600 * ***************************************************************************
601 * Routine: Gem_ramClear
602 *
603 * Purpose: Clear some structures.
604 *
605 * Notes: key==0 --> clear down to: min (vertices and simplices)
606 * key==1 --> clear down to: min + simplex ring
607 * key==2 --> clear down to: min + simplex ring + edges + edge ring
608 *
609 * Author: Michael Holst
610 * ***************************************************************************
611 */
Gem_ramClear(Gem * thee,int key)612 VPUBLIC void Gem_ramClear(Gem *thee, int key)
613 {
614 /* input check and some i/o */
615 VASSERT( (-1 <= key) && (key <= 2) );
616
617 /* clear the geometry structures */
618 if (key == 2) {
619 Vnm_print(0, "Gem_ramClear: clearing queues..");
620 Gem_resetSQ(thee,0);
621 Gem_resetSQ(thee,1);
622 Vnm_print(0, "..done.\n");
623 } else if (key == 1) {
624 Vnm_print(0, "Gem_ramClear: clearing edges..");
625 Gem_resetEE(thee);
626 Vnm_print(0, "..done.\n");
627 } else if (key == 0) {
628 } else if (key == -1) {
629 Vnm_print(0, "Gem_ramClear: clearing vertices..");
630 Gem_resetVV(thee);
631 Vnm_print(0, "..done.\n");
632 Vnm_print(0, "Gem_ramClear: clearing simplices..");
633 Gem_resetSS(thee);
634 Vnm_print(0, "..done.\n");
635 }
636 }
637
638 /*
639 * ***************************************************************************
640 * Routine: Gem_makeBnd
641 *
642 * Purpose: Force naborless faces to become boundary faces.
643 *
644 * Notes: btype==0 --> create interior boundary faces
645 * btype==1 --> create boundary faces of type "1"
646 * btype==2 --> create boundary faces of type "2"
647 *
648 * Author: Michael Holst
649 * ***************************************************************************
650 */
Gem_makeBnd(Gem * thee,int btype)651 VPUBLIC void Gem_makeBnd(Gem *thee, int btype)
652 {
653 int i, j, k, l, m, nabors;
654 VV *v[4];
655 SS *sm, *sm0, *sm1, *sm2;
656
657 /* input check and some i/o */
658 VASSERT( (0 <= btype) && (btype <= 2) );
659
660 /* go through all simplices and zero all boundary faces */
661 Vnm_print(0,"Gem_makeBnd: zeroing boundary faces/vertices..");
662 Gem_setNumBF(thee, 0);
663 Gem_setNumBV(thee, 0);
664 for (i=0; i<Gem_numSS(thee); i++) {
665 sm = Gem_SS(thee,i);
666 if ( (i>0) && (i % VPRTKEY) == 0 ) Vnm_print(0,"[BS:%d]",i);
667
668 /* get local vertices */
669 for (j=0; j<Gem_dimVV(thee); j++)
670 v[j] = SS_vertex(sm,j);
671
672 /* reset all vertices and faces to interior type */
673 for (j=0; j<Gem_dimVV(thee); j++) {
674 /* the other three local vertex/face numbers besides "j" */
675 k=(j+1) % Gem_dimVV(thee);
676 l=(k+1) % Gem_dimVV(thee);
677 m=(l+1) % Gem_dimVV(thee);
678 SS_setFaceType(sm, j, 0);
679 VV_setType(v[k], 0);
680 VV_setType(v[l], 0);
681 if (Gem_dim(thee) == 3) VV_setType(v[m], 0);
682 }
683 }
684 Vnm_print(0,"..done.\n");
685
686 /* are we done */
687 if (btype == 0) return;
688
689 /* okay now make a boundary */
690 Vnm_print(0,"Gem_makeBnd: rebuilding boundary faces/vertices..");
691 for (i=0; i<Gem_numSS(thee); i++) {
692 sm = Gem_SS(thee,i);
693 if ( (i>0) && (i % VPRTKEY) == 0 ) Vnm_print(0,"[BS:%d]",i);
694
695 /* get local vertices */
696 for (j=0; j<Gem_dimVV(thee); j++)
697 v[j] = SS_vertex(sm,j);
698
699 /* rebuild everything */
700 for (j=0; j<Gem_dimVV(thee); j++) {
701
702 /* the other three local vertex/face numbers besides "j" */
703 k=(j+1) % Gem_dimVV(thee);
704 l=(k+1) % Gem_dimVV(thee);
705 m=(l+1) % Gem_dimVV(thee);
706
707 /* look for a face nabor sharing face "j" (opposite vertex "j") */
708 nabors = 0;
709 for (sm0=VV_firstSS(v[k]); sm0!=VNULL;sm0=SS_link(sm0,v[k])) {
710 for (sm1=VV_firstSS(v[l]); sm1!=VNULL; sm1=SS_link(sm1,v[l])) {
711 if (Gem_dim(thee) == 2) {
712 if ((sm0!=sm) && (sm0==sm1)) nabors++;
713 } else {
714 for (sm2=VV_firstSS(v[m]); sm2!=VNULL;
715 sm2=SS_link(sm2,v[m])) {
716 if ((sm0!=sm) && (sm0==sm1) && (sm0==sm2)) {
717 nabors++;
718 }
719 }
720 }
721 }
722 }
723
724 /* if no one there, then face "j" is actually a boundary face */
725 if (nabors == 0) {
726 SS_setFaceType(sm, j, btype);
727 Gem_numBFpp(thee);
728 if (VINTERIOR( VV_type(v[k])) ) {
729 VV_setType(v[k], btype);
730 Gem_numBVpp(thee);
731 }
732 if (VINTERIOR( VV_type(v[l])) ) {
733 VV_setType(v[l], btype);
734 Gem_numBVpp(thee);
735 }
736 if (Gem_dim(thee) == 3) {
737 if (VINTERIOR( VV_type(v[m])) ) {
738 VV_setType(v[m], btype);
739 Gem_numBVpp(thee);
740 }
741 }
742 }
743 }
744 }
745 Vnm_print(0,"..done.\n");
746 }
747
748 /*
749 * ***************************************************************************
750 * Routine: Gem_makeBndExt
751 *
752 * Purpose: Mark selected boundary faces in a special way.
753 *
754 * Author: Michael Holst
755 * ***************************************************************************
756 */
Gem_makeBndExt(Gem * thee,int key)757 VPUBLIC void Gem_makeBndExt(Gem *thee, int key)
758 {
759 int i, j, k, l, m, p, q, nabors, btype, done, btypeGeneric;
760 VV *v[4];
761 SS *sm, *sm0, *sm1, *sm2;
762 double x[4][3], xchk;
763
764 /* go through all simplices and zero all boundary faces */
765 Vnm_print(0,"Gem_makeBnd: zeroing boundary faces/vertices..");
766 Gem_setNumBF(thee, 0);
767 Gem_setNumBV(thee, 0);
768 for (i=0; i<Gem_numSS(thee); i++) {
769 sm = Gem_SS(thee,i);
770 if ( (i>0) && (i % VPRTKEY) == 0 ) Vnm_print(0,"[BS:%d]",i);
771
772 /* get local vertices */
773 for (j=0; j<Gem_dimVV(thee); j++)
774 v[j] = SS_vertex(sm,j);
775
776 /* reset all vertices and faces to interior type */
777 for (j=0; j<Gem_dimVV(thee); j++) {
778 /* the other three local vertex/face numbers besides "j" */
779 k=(j+1) % Gem_dimVV(thee);
780 l=(k+1) % Gem_dimVV(thee);
781 m=(l+1) % Gem_dimVV(thee);
782 SS_setFaceType(sm, j, 0);
783 VV_setType(v[k], 0);
784 VV_setType(v[l], 0);
785 if (Gem_dim(thee) == 3) VV_setType(v[m], 0);
786 }
787 }
788 Vnm_print(0,"..done.\n");
789
790 /* okay now make a boundary */
791 Vnm_print(0,"Gem_makeBnd: rebuilding boundary faces/vertices..");
792 for (i=0; i<Gem_numSS(thee); i++) {
793 sm = Gem_SS(thee,i);
794 if ( (i>0) && (i % VPRTKEY) == 0 ) Vnm_print(0,"[BS:%d]",i);
795
796 /* get local vertices */
797 for (j=0; j<Gem_dimVV(thee); j++)
798 v[j] = SS_vertex(sm,j);
799
800 /* rebuild everything */
801 for (j=0; j<Gem_dimVV(thee); j++) {
802
803 /* the other three local vertex/face numbers besides "j" */
804 k=(j+1) % Gem_dimVV(thee);
805 l=(k+1) % Gem_dimVV(thee);
806 m=(l+1) % Gem_dimVV(thee);
807
808 /* look for a face nabor sharing face "j" (opposite vertex "j") */
809 nabors = 0;
810 for (sm0=VV_firstSS(v[k]); sm0!=VNULL;sm0=SS_link(sm0,v[k])) {
811 for (sm1=VV_firstSS(v[l]); sm1!=VNULL; sm1=SS_link(sm1,v[l])) {
812 if (Gem_dim(thee) == 2) {
813 if ((sm0!=sm) && (sm0==sm1)) nabors++;
814 } else {
815 for (sm2=VV_firstSS(v[m]); sm2!=VNULL;
816 sm2=SS_link(sm2,v[m])) {
817 if ((sm0!=sm) && (sm0==sm1) && (sm0==sm2)) {
818 nabors++;
819 }
820 }
821 }
822 }
823 }
824
825 /* if no one there, then face "j" is actually a boundary face */
826 if (nabors == 0) {
827
828 /* grab coordinates of the vertices of this face */
829 for (q=0; q<Gem_dim(thee); q++) {
830 x[0][q] = VV_coord(v[k],q);
831 }
832 for (q=0; q<Gem_dim(thee); q++) {
833 x[1][q] = VV_coord(v[l],q);
834 }
835 if (Gem_dim(thee) == 3) {
836 for (q=0; q<Gem_dim(thee); q++) {
837 x[2][q] = VV_coord(v[m],q);
838 }
839 }
840
841 /* default is interior; should not occur! */
842 btypeGeneric = 18;
843 done = 0;
844 btype = btypeGeneric;
845
846 /* ---------- check for base marking ---------- */
847 xchk = 0.0;
848 for (p=0; p<Gem_dim(thee); p++) {
849 xchk += VABS( x[p][1] - (-1.0) );
850 }
851 if (xchk < VSMALL) {
852 done = 1;
853 btype = 1;
854 }
855
856 /* ---------- check for base marking again ---------- */
857 xchk = 0.0;
858 for (p=0; p<Gem_dim(thee); p++) {
859 xchk += VABS( x[p][1] - ( 0.0) );
860 }
861 if (xchk < VSMALL) {
862 done = 1;
863 btype = 18;
864 }
865
866 /* ---------- check for first section ---------- */
867 if (!done) {
868 done = 1;
869 btype = 2;
870 for (p=0; p<Gem_dim(thee); p++) {
871 if (! ( ( 1.9 <= x[p][0])
872 && ( 6.1 >= x[p][0])
873 && (-VSMALL <= x[p][1])
874 && (-1.1 <= x[p][2])
875 && ( 1.1 >= x[p][2]) )) {
876 done = 0;
877 btype = btypeGeneric;
878 }
879 }
880
881 if (done) {
882 xchk = 0.0;
883 for (p=0; p<Gem_dim(thee); p++) {
884 xchk += VABS( x[p][1] - 10.0 );
885 }
886 if (xchk < VSMALL) {
887 btype = 10;
888 }
889 }
890 }
891
892 /* ---------- check for second section ---------- */
893 if (!done) {
894 done = 1;
895 btype = 4;
896 for (p=0; p<Gem_dim(thee); p++) {
897 if (! ( ( 7.9 <= x[p][0])
898 && (12.1 >= x[p][0])
899 && (-VSMALL <= x[p][1])
900 && (-1.1 <= x[p][2])
901 && ( 1.1 >= x[p][2]) )) {
902 done = 0;
903 btype = btypeGeneric;
904 }
905 }
906
907 if (done) {
908 xchk = 0.0;
909 for (p=0; p<Gem_dim(thee); p++) {
910 xchk += VABS( x[p][1] - 10.0 );
911 }
912 if (xchk < VSMALL) {
913 btype = 12;
914 }
915 }
916 }
917
918
919
920 /* ---------- check for third section ---------- */
921 if (!done) {
922 done = 1;
923 btype = 6;
924 for (p=0; p<Gem_dim(thee); p++) {
925 if (! ( (13.9 <= x[p][0])
926 && (18.1 >= x[p][0])
927 && (-VSMALL <= x[p][1])
928 && (-1.1 <= x[p][2])
929 && ( 1.1 >= x[p][2]) )) {
930 done = 0;
931 btype = btypeGeneric;
932 }
933 }
934
935 if (done) {
936 xchk = 0.0;
937 for (p=0; p<Gem_dim(thee); p++) {
938 xchk += VABS( x[p][1] - 10.0 );
939 }
940 if (xchk < VSMALL) {
941 btype = 14;
942 }
943 }
944 }
945
946
947 /* ---------- check for fourth section ---------- */
948 if (!done) {
949 done = 1;
950 btype = 8;
951 for (p=0; p<Gem_dim(thee); p++) {
952 if (! ( (19.9 <= x[p][0])
953 && (24.1 >= x[p][0])
954 && (-VSMALL <= x[p][1])
955 && (-1.1 <= x[p][2])
956 && ( 1.1 >= x[p][2]) )) {
957 done = 0;
958 btype = btypeGeneric;
959 }
960 }
961
962 if (done) {
963 xchk = 0.0;
964 for (p=0; p<Gem_dim(thee); p++) {
965 xchk += VABS( x[p][1] - 10.0 );
966 }
967 if (xchk < VSMALL) {
968 btype = 16;
969 }
970 }
971 }
972
973 /* should have been marked with SOME boundary type */
974 VASSERT( 0 != btype );
975
976 /* set the facetype */
977 SS_setFaceType(sm, j, btype);
978 Gem_numBFpp(thee);
979
980 /* set the vertex types (dirichlet overrides robin) */
981 if (!VDIRICHLET( VV_type(v[k])) ) {
982 if (VINTERIOR( VV_type(v[k])) ) {
983 Gem_numBVpp(thee);
984 }
985 VV_setType(v[k], btype);
986 }
987 if (!VDIRICHLET( VV_type(v[l])) ) {
988 if (VINTERIOR( VV_type(v[l])) ) {
989 Gem_numBVpp(thee);
990 }
991 VV_setType(v[l], btype);
992 }
993 if (Gem_dim(thee) == 3) {
994 if (!VDIRICHLET( VV_type(v[m])) ) {
995 if (VINTERIOR( VV_type(v[m])) ) {
996 Gem_numBVpp(thee);
997 }
998 VV_setType(v[m], btype);
999 }
1000 }
1001 }
1002 }
1003 }
1004 Vnm_print(0,"..done.\n");
1005 }
1006
1007