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: gemio.c,v 1.24 2010/08/12 05:18:57 fetk Exp $"
21 * ***************************************************************************
22 */
23
24 /*
25 * ***************************************************************************
26 * File: Gemio.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: gemio.c,v 1.24 2010/08/12 05:18:57 fetk Exp $")
37
38 /* comment and whitespace characters */
39 VPRIVATE char *MCwhiteChars = " =,;\t\n";
40 VPRIVATE char *MCcommChars = "#%";
41
42 /* the MC copyright */
43 /* (NOTE: Only 509 char-length constant strings supported in ISO C) */
44 VPRIVATE char *MCcopyRightPCT1 =
45 "%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n\
46 %%% MC = < Manifold Code >\n\
47 %%% Copyright (C) 1994-- Michael Holst\n\
48 %%%\n";
49
50 VPRIVATE char *MCcopyRightPCT2 =
51 "%%% This program is free software; you can redistribute it and/or modify it\n\
52 %%% under the terms of the GNU General Public License as published by the\n\
53 %%% Free Software Foundation; either version 2 of the License, or (at your\n\
54 %%% option) any later version.\n\
55 %%%\n";
56
57 VPRIVATE char *MCcopyRightPCT3 =
58 "%%% This program is distributed in the hope that it will be useful,\n\
59 %%% but WITHOUT ANY WARRANTY; without even the implied warranty of\n\
60 %%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.\n\
61 %%% See the GNU General Public License for more details.\n\
62 %%%\n";
63
64 VPRIVATE char *MCcopyRightPCT4 =
65 "%%% You should have received a copy of the GNU General Public License along\n\
66 %%% with this program; if not, write to the Free Software Foundation, Inc.,\n\
67 %%% 675 Mass Ave, Cambridge, MA 02139, USA.\n\
68 %%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n";
69
70 /*
71 * ***************************************************************************
72 * Class Gem: Inlineable methods
73 * ***************************************************************************
74 */
75 #if !defined(VINLINE_GEM)
76
77 #endif /* if !defined(VINLINE_GEM) */
78 /*
79 * ***************************************************************************
80 * Class Gem: Non-inlineable methods
81 * ***************************************************************************
82 */
83
84 /*
85 * ***************************************************************************
86 * Routine: Gem_read
87 *
88 * Purpose: Read in the user-specified initial vertex-simplex mesh.
89 * provided to us in MC-Simplex-Format (MCSF), and transform
90 * into our internal datastructures.
91 *
92 * Parms: key ==> input format type
93 * key=0 ==> simplex format
94 * key=1 ==> edge format
95 * key=2 ==> simplex-nabor format
96 *
97 * Notes: The user provides the following information about a domain and an
98 * initial simplex-triangulation:
99 *
100 * D = DIMENSION = spatial dimension of problem (1, 2, or 3)
101 * N = NVERTICES = total number of vertices in mesh
102 * L = NSIMPLICES = total number of simplices in the mesh
103 *
104 * double vertex* = List of all vertex coordinates in the form:
105 *
106 * vertex[ N * (3+D) ] = {
107 * { id1, chart1, v1_x1, ..., v1_xD },
108 * { id2, chart2, v2_x1, ..., v2_xD },
109 * ...
110 * { idN, chartN, vN_x1, ..., vN_xD }
111 * }
112 *
113 * Here, vj_xk is the kth component (float or double) of the vertex
114 * with global vertex number j. The "idj" flag is a 32-bit integer
115 * "name" for vertex j, "chartj" is a 32-bit number representing the
116 * "chart" with which to interpret the coordinates vj_xk, in the sense
117 * of the charts of an atlas of manifold domain.
118 *
119 * int simplex* = List of all simplices by vertex number in the form:
120 *
121 * simplex[ L * (2 + 2*(D+1)) ] = {
122 * { id1, g1, m1, f1_1, ..., f1_{D+1}, s1_v1, ..., s1_v{D+1} },
123 * { id2, g2, m2, f2_1, ..., f2_{D+1}, s2_v1, ..., s2_v{D+1} },
124 * ...
125 * { idL, gL, mL, fL_1, ..., f2_{D+1}, sL_v1, ..., sL_v{D+1} }
126 * }
127 *
128 * Here, sj_vk is a 32-bit integer giving the global vertex number
129 * making up the kth vertex of the simplex with global number j. The
130 * "idj" flag is a 32-bit integer "name" for simplex j, "gj" is a
131 * 32-bit group number associated with simplex j (for grouping subsets
132 * of simplices together for various reasons), "mj" is a 32-bit integer
133 * containing the information about the simplex j such as its material
134 * type, and "fj_k" is a 32-bit integer containing the information
135 * about each face k of simplex j such as their boundary types (each
136 * face opposite vertex k in simplex j).
137 *
138 * Thus, in 2D, a simplex (triangle) is specified by 3 consecutive
139 * vertex numbers, and in 3D a simplex (tetrahedra) is specified by 4
140 * consecutive vertex numbers. The physical coordinates of any vertex
141 * k in the simplex array are given in the vertex array in the
142 * appropriate row of the array.
143 *
144 * NOTE: The ordering of the vertices in a simplex is *extremely*
145 * important here; see the note below.
146 *
147 * Ordering of the vertices in a simplex:
148 *
149 * 1D: Well, this is pretty straight-forward; we will order the vertices
150 * from left-to-right in each simplex; this will produce the correct
151 * sign in integration by parts.
152 *
153 * 2D: All closed triangles must be specified by three consecutive vertices
154 * in simplex and must be counter-clockwise-ordered by their vertices,
155 * as seen from the "up" side of the 2D body/shell/surface. This
156 * produces the correct surface-normals from the right-hand-rule for
157 * surface (line in 2D) integrals.
158 *
159 * 3D: All closed tetrahedra must be specified by four consecutive vertices
160 * in "simplex" in the following way: The first three vertices must
161 * represent any one of the four faces as a counter-clockwise-ordered
162 * triangle, as seen from INSIDE the tetrahedra. I.e., you can think
163 * about this first triangle as lying in the plane, and you are
164 * standing on the plane looking down at it. The fourth vertex
165 * specified following the first three in the simplex is then some
166 * height above the plane containing the first counter-clockwise
167 * ordered triangle. This specification allows the remaining three
168 * (of the four) triangles making up any tetrahedra to be correctly
169 * specified in a counter-clockwise, inward-facing manner, so that the
170 * correct surface-normals can be calculated consistently.
171 *
172 * Author: Michael Holst
173 * ***************************************************************************
174 */
Gem_read(Gem * thee,int key,Vio * sock)175 VPUBLIC int Gem_read(Gem *thee, int key, Vio *sock)
176 {
177 int i, j, theDim, theDimII, inum;
178 int numVV, indexV, chartV;
179 int numSS, indexS, groupS, materialS, faceTypeS[4], sv[4];
180 int vnum, vtp, vtpI;
181 int fnum[3], ftp[3], ftpB[3];
182 double x[3];
183 char tok[VMAX_BUFSIZE];
184 VV *vx;
185 SS *sm;
186
187 /* some error checking */
188 VJMPERR1( VNULL != sock );
189 VJMPERR2( (0<=key) && (key<=2) );
190
191 if (key == 1) {
192 return Gem_readEdge(thee,sock);
193 } else if (key == 2) {
194 Vnm_print(2,"Gem_read: Cannot handle simplex-nabor format\n");
195 VJMPERR2( 0 );
196 }
197
198 /* deal with socket or file setup */
199 Vio_setWhiteChars(sock, MCwhiteChars);
200 Vio_setCommChars(sock, MCcommChars);
201
202 /* get the "mcsf_begin=1" tag */
203 VJMPERR2( 1==Vio_scanf(sock,"%s",tok) );
204 VJMPERR2( !strcmp(tok,"mcsf_begin") );
205 VJMPERR2( 1==Vio_scanf(sock,"%d",&inum) );
206 VJMPERR2( inum == 1 );
207
208 /* read the intrinsic dimension */
209 VJMPERR2( 1==Vio_scanf(sock,"%s",tok) );
210 VJMPERR2( !strcmp(tok,"dim") );
211 VJMPERR2( 1==Vio_scanf(sock,"%d",&theDim) );
212 VJMPERR2( (theDim == 2) || (theDim == 3) );
213
214 /* read the imbedding dimension */
215 VJMPERR2( 1==Vio_scanf(sock,"%s",tok) );
216 VJMPERR2( !strcmp(tok,"dimii") );
217 VJMPERR2( 1==Vio_scanf(sock,"%d",&theDimII) );
218 VJMPERR2( (theDimII == 2) || (theDimII == 3) );
219
220 /* read and set the number of vertices information */
221 VJMPERR2( 1==Vio_scanf(sock,"%s",tok) );
222 VJMPERR2( !strcmp(tok,"vertices") );
223 VJMPERR2( 1==Vio_scanf(sock,"%d",&numVV) );
224 thee->numVV0 = numVV;
225
226 /* read and set the number of simplices information */
227 VJMPERR2( 1==Vio_scanf(sock,"%s",tok) );
228 VJMPERR2( !strcmp(tok,"simplices") );
229 VJMPERR2( 1==Vio_scanf(sock,"%d",&numSS) );
230
231 /* reset ourself to handle the new manifold dimensions */
232 Vnm_print(0,"Gem_read: theDim=%d, theDimII=%d, numVV=%d, numSS=%d\n",
233 theDim, theDimII, numVV, numSS);
234 Vnm_print(0,"Gem_read: Reseting manifold structures.\n");
235 Gem_reset( thee, theDim, theDimII );
236
237 /* read in the vertices */
238 Vnm_print(0,"Gem_read: Reading ..vertices..");
239 VJMPERR2( 1==Vio_scanf(sock,"%s",tok) );
240 VJMPERR2( !strcmp(tok,"vert") );
241 VJMPERR2( 1==Vio_scanf(sock,"%s",tok) );
242 VJMPERR2( !strcmp(tok,"[") );
243 for (i=0; i<numVV; i++) {
244
245 /* grab the three integer info fields */
246 VJMPERR2( 1==Vio_scanf(sock,"%d",&indexV) );
247 VJMPERR2( 1==Vio_scanf(sock,"%d",&chartV) );
248
249 /* grab the "dimii" coordinates */
250 for (j=0; j<theDimII; j++)
251 VJMPERR2( 1==Vio_scanf(sock,"%le",&x[j]) );
252
253 /* trivially embed any 2D coordinates into 3D */
254 if (theDimII == 2) x[2] = 0.;
255
256 /* create the new vertex */
257 vx = Gem_createAndInitVV(thee);
258 VV_setReality(vx, 0);
259 VV_setDim(vx, theDim);
260 VV_setClass(vx, 0);
261 VV_setType(vx, 0);
262 VV_setId(vx, i);
263 VV_setChart(vx, chartV);
264
265 /* set the vertex coordinates */
266 for (j=0; j<3; j++)
267 VV_setCoord(vx, j, x[j]);
268 }
269 VJMPERR2( 1==Vio_scanf(sock,"%s",tok) );
270 VJMPERR2( !strcmp(tok,"]") );
271
272 /* read in the simplices */
273 Vnm_print(0,"..simplices..");
274 VJMPERR2( 1==Vio_scanf(sock,"%s",tok) );
275 VJMPERR2( !strcmp(tok,"simp") );
276 VJMPERR2( 1==Vio_scanf(sock,"%s",tok) );
277 VJMPERR2( !strcmp(tok,"[") );
278 for (i=0; i<numSS; i++) {
279
280 /* grab the three integer info fields */
281 VJMPERR2( 1==Vio_scanf(sock,"%d",&indexS) );
282 VJMPERR2( 1==Vio_scanf(sock,"%d",&groupS) );
283 VJMPERR2( 1==Vio_scanf(sock,"%d",&materialS) );
284
285 /* grab the "dim+1" face types */
286 for (j=0; j<theDim+1; j++)
287 VJMPERR2( 1==Vio_scanf(sock,"%d",&faceTypeS[j]) );
288
289 /* grab the "dim+1" vertex numbers */
290 for (j=0; j<theDim+1; j++)
291 VJMPERR2( 1==Vio_scanf(sock,"%d",&sv[j]) );
292
293 /* create the new simplex */
294 sm = Gem_createAndInitSS(thee);
295 SS_setReality(sm, 0);
296 SS_setDim(sm, theDim);
297 SS_setClass(sm, theDim);
298 SS_setType(sm, materialS);
299 SS_setId(sm, i);
300 SS_setChart(sm, groupS);
301
302 /* set the simplex face types and vertex labels */
303 for (j=0; j<theDim+1; j++) {
304 SS_setFaceType( sm, j, faceTypeS[j] );
305 if (VBOUNDARY( faceTypeS[j] )) {
306 (thee->numBF)++;
307 }
308 vx = Gem_VV(thee, sv[j] );
309 SS_setVertex( sm, j, vx );
310 }
311
312 /* calculate (our contribution to) vertex types from our face types */
313 for (j=0; j<theDim+1; j++) {
314
315 /* get the vertex in question */
316 vnum = j;
317 vx = SS_vertex( sm, vnum );
318
319 /* get face numbers of two/three faces which touch vertex vnum */
320 fnum[0] = vmapOV3[vnum][0];
321 fnum[1] = vmapOV3[vnum][1];
322 fnum[2] = vmapOV3[vnum][2]; /* 2D: third face always interior */
323
324 /* some shorthand notation... */
325 vtp = VV_type(vx);
326 vtpI = VINTERIOR( VV_type(vx) );
327 ftp[0] = SS_faceType(sm,fnum[0]);
328 ftp[1] = SS_faceType(sm,fnum[1]);
329 ftp[2] = SS_faceType(sm,fnum[2]);
330 ftpB[0] = VBOUNDARY( SS_faceType(sm,fnum[0]) );
331 ftpB[1] = VBOUNDARY( SS_faceType(sm,fnum[1]) );
332 ftpB[2] = VBOUNDARY( SS_faceType(sm,fnum[2]) );
333
334 /* if any of the faces are Boundary, then mark vertex Boundary */
335 if ( ftpB[0] || ftpB[1] || ftpB[2] ) {
336
337 /* deal with existing vertex type */
338 if (vtpI) {
339 (thee->numBV)++;
340 } else {
341 /* this is okay; may already be boundary */
342 }
343
344 /* okay, determine max boundary flag (including vtp) */
345 if (ftpB[0]) vtp = VMAX2(vtp,ftp[0]);
346 if (ftpB[1]) vtp = VMAX2(vtp,ftp[1]);
347 if (ftpB[2]) vtp = VMAX2(vtp,ftp[2]);
348
349 /* set the type */
350 VV_setType(vx, vtp);
351
352 /* otherwise, if all faces Interior, just leave vertex alone */
353 }
354 }
355
356 /* build the ringed vertex datastructure */
357 SS_buildRing(sm);
358 }
359 VJMPERR2( 1==Vio_scanf(sock,"%s",tok) );
360 VJMPERR2( !strcmp(tok,"]") );
361
362 /* get the "mcsf_end" tag */
363 VJMPERR2( 1==Vio_scanf(sock,"%s",tok) );
364 VJMPERR2( !strcmp(tok,"mcsf_end") );
365 VJMPERR2( 1==Vio_scanf(sock,"%d",&inum) );
366 VJMPERR2( inum == 1 );
367
368 /* finish the i/o */
369 Vnm_print(0,"..done.\n");
370
371 /* create initial edge markings in the simplices */
372 Gem_markEdges(thee);
373
374 /* count v/e/f/s and check the mesh */
375 Gem_countChk(thee);
376
377 /* return with no errors */
378 return 1;
379
380 VERROR1:
381 Vnm_print(2,"Gem_read: Problem with I/O socket.\n");
382 Vnm_print(2,"Gem_read: Reseting manifold structures (bailing out).\n");
383 Gem_reset( thee, 0, 0 );
384 return 0;
385
386 VERROR2:
387 Vnm_print(2,"Gem_read: Format problem with input\n");
388 Vnm_print(2,"Gem_read: Reseting manifold structures (bailing out).\n");
389 Gem_reset( thee, 0, 0 );
390 return 0;
391 }
392
393 /*
394 * ***************************************************************************
395 * Routine: Gem_write
396 *
397 * Purpose: Write out a 2-simplex or 3-simplex mesh in "MCSF" format.
398 *
399 * Parms: key ==> output format type
400 * key=0 ==> simplex format
401 * key=1 ==> edge format
402 * key=2 ==> simplex-nabor format
403 * fkey ==> simplex write option
404 * fkey=0 ==> write simplices
405 * fkey=1 ==> write only simplex boundary faces
406 * fkey=2 ==> write only simplices that have
407 * at least one boundary face
408 * (NOT IMPLIMENTED HERE)
409 *
410 * Author: Michael Holst
411 * ***************************************************************************
412 */
Gem_write(Gem * thee,int key,Vio * sock,int fkey)413 VPUBLIC void Gem_write(Gem *thee, int key, Vio *sock, int fkey)
414 {
415 int i, j, k, nid[4];
416 VV *vx;
417 SS *sm;
418
419 /* some error checking */
420 VJMPERR1( VNULL != sock );
421 VJMPERR2( Gem_numVV(thee) > 0 );
422 VJMPERR2( (0<=key) && (key<=2) );
423 VJMPERR2( (0<=fkey) && (fkey<=2) );
424
425 if (key == 1) {
426 Gem_writeEdge(thee, sock);
427 } else if ((fkey==1) && (Gem_dim(thee) == 3)) {
428 Gem_writeFace3d(thee, sock);
429 } else if (fkey==2) {
430 Vnm_print(2,"Gem_write: fkey=<%d> not implemented (not useful).\n",
431 fkey);
432 } else { /* (fkey==0) */
433
434 /* deal with socket or file setup */
435 Vio_setWhiteChars(sock, MCwhiteChars);
436 Vio_setCommChars(sock, MCcommChars);
437
438 /* write out the header */
439 Vio_printf(sock,"%s",MCcopyRightPCT1);
440 Vio_printf(sock,"%s",MCcopyRightPCT2);
441 Vio_printf(sock,"%s",MCcopyRightPCT3);
442 Vio_printf(sock,"%s",MCcopyRightPCT4);
443 Vio_printf(sock,"\n");
444 Vio_printf(sock,"%s=%d;\n\n","mcsf_begin",1);
445
446 Vio_printf(sock," %s=%d; %% intrinsic manifold dimension\n",
447 "dim",Gem_dim(thee));
448
449 Vio_printf(sock," %s=%d; %% imbedding manifold dimension\n",
450 "dimii",Gem_dimII(thee));
451
452 Vio_printf(sock," %s=%d; %% number of vertices\n",
453 "vertices",Gem_numVV(thee));
454
455 Vio_printf(sock,"%s=%d; %% number of simplices\n",
456 "simplices",Gem_numSS(thee));
457
458 /* write out the vertex data */
459 Vio_printf(sock,"\n%s=%s\n","vert","[");
460 Vio_printf(sock,"%%-------- ---- -----------------"
461 " ----------------- -----------------\n");
462 Vio_printf(sock,
463 "%% Vert-ID Chrt X0-Coordinate X1-Coordinate X2-Coordinate\n");
464 Vio_printf(sock,"%%-------- ---- -----------------"
465 " ----------------- -----------------\n");
466 for (i=0; i<Gem_numVV(thee); i++) {
467 vx = Gem_VV(thee,i);
468 Vio_printf(sock,"%-9d %-4d", VV_id(vx), VV_chart(vx) );
469 for (j=0; j<Gem_dimII(thee); j++)
470 Vio_printf(sock," %17.10le", VV_coord(vx,j) );
471 Vio_printf(sock,"\n");
472 }
473 Vio_printf(sock,"%s;\n","]");
474
475 /* write out the simplex data */
476 Vio_printf(sock,"\n%s=%s\n","simp","[");
477 Vio_printf(sock,"%%-------- ---- ---- -------------------"
478 " ---------------------------------------\n");
479 Vio_printf(sock,
480 "%% Simp-ID Grp Mat Face-Types Vertex-Numbers\n");
481 Vio_printf(sock,"%%-------- ---- ---- -------------------"
482 " ---------------------------------------\n");
483 for (i=0; i<Gem_numSS(thee); i++) {
484 sm = Gem_SS(thee,i);
485 Vio_printf(sock,"%-9d %-4d %-4d ", i, SS_chart(sm), SS_type(sm));
486 if (Gem_dimVV(thee) == 3) {
487 Vio_printf(sock,"%-4d %-4d %-4d %d %d %d",
488 SS_faceType(sm,0),
489 SS_faceType(sm,1),
490 SS_faceType(sm,2),
491 VV_id( SS_vertex(sm,0) ),
492 VV_id( SS_vertex(sm,1) ),
493 VV_id( SS_vertex(sm,2) )
494 );
495 if (key == 0) {
496 Vio_printf(sock,"\n");
497 } else if (key == 2) {
498 for (k=0; k<3; k++) {
499 if (SS_nabor(sm,k) != VNULL) {
500 nid[k] = SS_id( SS_nabor(sm,k) );
501 } else {
502 nid[k] = -1;
503 }
504 }
505 Vio_printf(sock," %d %d %d\n",
506 nid[0], nid[1], nid[2]);
507 } else { VASSERT(0); }
508 } else if (Gem_dimVV(thee) == 4) {
509 Vio_printf(sock,"%-4d %-4d %-4d %-4d %d %d %d %d",
510 SS_faceType(sm,0),
511 SS_faceType(sm,1),
512 SS_faceType(sm,2),
513 SS_faceType(sm,3),
514 VV_id( SS_vertex(sm,0) ),
515 VV_id( SS_vertex(sm,1) ),
516 VV_id( SS_vertex(sm,2) ),
517 VV_id( SS_vertex(sm,3) )
518 );
519 if (key == 0) {
520 Vio_printf(sock,"\n");
521 } else if (key == 2) {
522 for (k=0; k<4; k++) {
523 if (SS_nabor(sm,k) != VNULL) {
524 nid[k] = SS_id( SS_nabor(sm,k) );
525 } else {
526 nid[k] = -1;
527 }
528 }
529 Vio_printf(sock," %d %d %d %d\n",
530 nid[0], nid[1], nid[2], nid[3]);
531 } else { VASSERT(0); }
532 } else { VASSERT(0); }
533 }
534 Vio_printf(sock,"%s;\n","]");
535
536 /* finish the file */
537 Vio_printf(sock,"\n%s=%d;\n\n","mcsf_end",1);
538 }
539
540 /* return with no errors */
541 return;
542
543 VERROR1:
544 Vnm_print(2,"Gem_write: Problem with I/O socket.\n");
545
546 VERROR2:
547 Vnm_print(2,"Gem_write: Detected a problem (bailing out).\n");
548 return;
549 }
550
551 /*
552 * ***************************************************************************
553 * Routine: Gem_writeFace3d
554 *
555 * Purpose: Write out the faces of a 3-simplex mesh as a complete and legal
556 * 2-simplex mesh in "MCSF" format (described above for Gem_read).
557 *
558 * Parms: iodev ==> output device type (file/buff/unix/inet)
559 * iofmt ==> output device format (ascii/xdr)
560 * thost ==> output hostname (for sockets)
561 * fname ==> output file/buff/unix/inet name
562 *
563 * Author: Michael Holst
564 * ***************************************************************************
565 */
Gem_writeFace3d(Gem * thee,Vio * sock)566 VPUBLIC void Gem_writeFace3d(Gem *thee, Vio *sock)
567 {
568 int i, j, ii, numVV, numSS, *vTmp;
569 VV *vx;
570 SS *sm;
571
572 /* some error checking */
573 VJMPERR1( VNULL != sock );
574 VJMPERR2( Gem_numVV(thee) > 0 );
575 VJMPERR2( Gem_dim(thee) == 3 );
576 VJMPERR2( Gem_dimII(thee) == 3 );
577
578 /* number of boundary vertices and faces */
579 numVV = Gem_numBV(thee);
580 numSS = Gem_numBF(thee);
581
582 /* deal with socket or file setup */
583 Vio_setWhiteChars(sock, MCwhiteChars);
584 Vio_setCommChars(sock, MCcommChars);
585
586 /* write out the header */
587 Vio_printf(sock,"%s",MCcopyRightPCT1);
588 Vio_printf(sock,"%s",MCcopyRightPCT2);
589 Vio_printf(sock,"%s",MCcopyRightPCT3);
590 Vio_printf(sock,"%s",MCcopyRightPCT4);
591 Vio_printf(sock,"\n");
592 Vio_printf(sock,"%s=%d;\n\n","mcsf_begin",1);
593
594 Vio_printf(sock," %s=%d; %% intrinsic manifold dimension\n",
595 "dim",2);
596
597 Vio_printf(sock," %s=%d; %% imbedding manifold dimension\n",
598 "dimii",3);
599
600 Vio_printf(sock," %s=%d; %% number of vertices\n",
601 "vertices",numVV);
602
603 Vio_printf(sock,"%s=%d; %% number of simplices\n",
604 "simplices",numSS);
605
606 /* write out the vertex data */
607 Vio_printf(sock,"\n%s=%s\n","vert","[");
608 Vio_printf(sock,"%%-------- ---- -----------------"
609 " ----------------- -----------------\n");
610 Vio_printf(sock,
611 "%% Vert-ID Chrt X0-Coordinate X1-Coordinate X2-Coordinate\n");
612 Vio_printf(sock,"%%-------- ---- -----------------"
613 " ----------------- -----------------\n");
614 vTmp = Vmem_malloc( thee->vmem, Gem_numVV(thee), sizeof(int) );
615 ii = 0;
616 for (i=0; i<Gem_numVV(thee); i++) {
617 vx = Gem_VV(thee,i);
618 if (VBOUNDARY( VV_type(vx) )) {
619 vTmp[VV_id(vx)] = ii;
620 Vio_printf(sock,"%-9d %-4d", ii, VV_chart(vx) );
621 for (j=0; j<Gem_dimII(thee); j++)
622 Vio_printf(sock," %17.10le", VV_coord(vx,j) );
623 Vio_printf(sock,"\n");
624 ii++;
625 }
626 }
627 VASSERT( ii == numVV );
628 Vio_printf(sock,"%s;\n","]");
629
630 /* write out the simplex data */
631 Vio_printf(sock,"\n%s=%s\n","simp","[");
632 Vio_printf(sock,"%%-------- ---- ---- -------------------"
633 " ---------------------------------------\n");
634 Vio_printf(sock,
635 "%% Simp-ID Grp Mat Face-Types Vertex-Numbers\n");
636 Vio_printf(sock,"%%-------- ---- ---- -------------------"
637 " ---------------------------------------\n");
638 ii = 0;
639 for (i=0; i<Gem_numSS(thee); i++) {
640 sm = Gem_SS(thee,i);
641 for (j=0; j<4; j++) {
642 if (VBOUNDARY( SS_faceType(sm,j) )) {
643 Vio_printf(sock,"%-9d %-4d %-4d ",
644 ii, SS_chart(sm), SS_type(sm));
645 Vio_printf(sock,"%-4d %-4d %-4d %d %d %d\n",
646 0, 0, 0,
647 vTmp[ VV_id( SS_vertex(sm,vmapF[j][0]) ) ],
648 vTmp[ VV_id( SS_vertex(sm,vmapF[j][1]) ) ],
649 vTmp[ VV_id( SS_vertex(sm,vmapF[j][2]) ) ] );
650 ii++;
651 }
652 }
653 }
654 VASSERT( ii == numSS );
655 Vio_printf(sock,"%s;\n","]");
656
657 /* finish the file */
658 Vio_printf(sock,"\n%s=%d;\n\n","mcsf_end",1);
659
660 /* release the temporary storage */
661 Vmem_free( thee->vmem, Gem_numVV(thee), sizeof(int), (void**)&vTmp );
662
663 /* return with no errors */
664 Vnm_print(0,"Gem_writeFace3d: Finished writing output\n");
665 return;
666
667 VERROR1:
668 Vnm_print(2,"Gem_writeFace3d: Problem with I/O socket.\n");
669
670 VERROR2:
671 Vnm_print(2,"Gem_writeFace3d: Detected a problem (bailing out).\n");
672 return;
673 }
674
675 /*
676 * ***************************************************************************
677 * Routine: Gem_readEdge
678 *
679 * Purpose: Read in the user-specified initial vertex-edge mesh
680 * provided to us in MC-Edge-Format (MCEF), and transform
681 * into our internal datastructures.
682 *
683 * Notes: The user provides the following information about a domain and an
684 * initial edge-based triangulation:
685 *
686 * D = DIMENSION = spatial dimension of problem (1, 2, or 3)
687 * N = NVERTICES = total number of vertices in mesh
688 * L = NEDGES = total number of edges in the mesh
689 *
690 * double vertex* = List of all vertex coordinates in the form:
691 *
692 * vertex[ N * (3+D) ] = {
693 * { id1, proc1, info1, v1_x1, ..., v1_xD },
694 * { id2, proc2, info2, v2_x1, ..., v2_xD },
695 * ...
696 * { idN, procN, infoN, vN_x1, ..., vN_xD }
697 * }
698 *
699 * Here, vj_xk is the kth component (float or double) of the vertex
700 * with global vertex number j. The "idj" flag is a 32-bit integer
701 * "name" for vertex j, "procj" is a 32-bit color or processor number
702 * associated with vertex j, and "infoj" is a 32-bit integer containing
703 * information about vertex j.
704 *
705 * int edge* = List of all edges by vertex number in the form:
706 *
707 * edge[ L * (2 + 2*(D+1)) ] = {
708 * { id1, proc1, info1, e1_v1, e1_v2 },
709 * { id2, proc2, info2, e2_v1, e2_v2 },
710 * ...
711 * { idL, procL, infoL, eL_v1, eL_v2 }
712 * }
713 *
714 * Here, ej_vk is a 32-bit integer giving the global vertex number
715 * making up the kth vertex of the edge with global number j. The
716 * "idj" flag is a 32-bit integer "name" for edge j, "procj" is a
717 * 32-bit color or processor number associated with edge j,
718 * "infoj" is a 32-bit integer containing the information about the
719 * edge j.
720 *
721 * Notes: To recover a simplex (triangle) mesh from the vertices and edges,
722 * we must traverse them and rebuild the simplex structure, and do it
723 * all in linear time.
724 *
725 * The three-step simplex recovery algorithm is as follows:
726 *
727 * (1) Definite vertices and edges from the input data.
728 *
729 * (2) Create all possible simplices as follows:
730 *
731 * For (vx=firstVV; vx!=lastVV; vx=nextVV) {
732 * | Build all simplices which use vx as a vertex as follows:
733 * | 2D CASE: For each distinct pair of vertices connected
734 * | by an edge with vx, if this pair are also
735 * | connected by an edge, then the three make a
736 * | triangle. If the triangle has not already
737 * | been created, then do so and add to the
738 * | simplex rings for each of the three vertices.
739 * | 3D CASE: For each distinct trio of vertices connected
740 * | by an edge with vx, if this trio forms a triangle,
741 * | then the foursome make a tetrahedron. If the
742 * | tetrahedron has not already been created, then do
743 * | so and add to the simplex rings for each of the
744 * | four vertices.
745 * EndFor
746 *
747 * (3) Remove a few "bad" simplices which were created incorrectly
748 * in Step (2) above as follows:
749 *
750 * For (sm=firstSS; sm!=lastSS; sm=nextSS) {
751 * | Remove all "bad" simplices, defined to be those satisfying:
752 * | 2D CASE: All interior edges have multiple nabors, and
753 * | all boundary edges have at least one nabor.
754 * | 3D CASE: All interior faces have multiple nabors, and
755 * | all boundary faces have at least one nabor.
756 * EndFor
757 *
758 * NOTE: The crucial Step (3) of the above algorithm only works
759 * correctly if we correctly identify the boundary faces of the
760 * simplices, which is only possible if the input edge-based mesh
761 * was given with boundary edge information. (We can construct
762 * correct face types from the types of the one or three edges
763 * forming the face in 2D or 3D respectively.)
764 *
765 * We attempt to recover simplex face types from the given
766 * edge types. However, (at least) one degenerate cases is not
767 * recoverable: an isolated Neumann face surounded by Dirichlet faces
768 * will appear simply as three Dirichlet edges in the edge-based mesh,
769 * and we turn this into a Dirichlet face. Note that if a Neumann
770 * face consists of more than one simplex face, then it will be
771 * recovered correctly from the edge types.
772 *
773 * Note also that if a non-simplex mesh is provided as input in
774 * as an edge-based mesh, we will build the edges as specified,
775 * but our attempt to build simplices will only recover those
776 * simplices which actually exist in the edge mesh. Any other
777 * non-simplex polyedra will appear as holes in the final
778 * simplex mesh.
779 *
780 * Author: Michael Holst
781 * ***************************************************************************
782 */
Gem_readEdge(Gem * thee,Vio * sock)783 VPUBLIC int Gem_readEdge(Gem *thee, Vio *sock)
784 {
785 int i, j, theDim, theDimII, inum;
786 int numVV, indexV, chartV;
787 int numEE, indexE, groupE, typeE, se[2];
788 double x[3];
789 char tok[VMAX_BUFSIZE];
790 VV *vx;
791 EE *eg;
792
793 /* some error checking */
794 VJMPERR1( VNULL != sock );
795
796 /* deal with socket or file setup */
797 Vio_setWhiteChars(sock, MCwhiteChars);
798 Vio_setCommChars(sock, MCcommChars);
799
800 /* get the "mcef_begin=1" tag */
801 VJMPERR2( 1==Vio_scanf(sock,"%s",tok) );
802 VJMPERR2( !strcmp(tok,"mcef_begin") );
803 VJMPERR2( 1==Vio_scanf(sock,"%d",&inum) );
804 VJMPERR2( inum == 1 );
805
806 /* read the intrinsic dimension */
807 VJMPERR2( 1==Vio_scanf(sock,"%s",tok) );
808 VJMPERR2( !strcmp(tok,"dim") );
809 VJMPERR2( 1==Vio_scanf(sock,"%d", &theDim) );
810 VJMPERR2( (theDim == 2) || (theDim == 3) );
811
812 /* read the imbedding dimension */
813 VJMPERR2( 1==Vio_scanf(sock,"%s",tok) );
814 VJMPERR2( !strcmp(tok,"dimii") );
815 VJMPERR2( 1==Vio_scanf(sock,"%d", &theDimII) );
816 VJMPERR2( (theDimII == 2) || (theDimII == 3) );
817
818 /* read and set the number of vertices information */
819 VJMPERR2( 1==Vio_scanf(sock,"%s",tok) );
820 VJMPERR2( !strcmp(tok,"vertices") );
821 VJMPERR2( 1==Vio_scanf(sock,"%d", &numVV) );
822 thee->numVV0 = numVV;
823
824 /* read and set the number of edges information */
825 VJMPERR2( 1==Vio_scanf(sock,"%s",tok) );
826 VJMPERR2( !strcmp(tok,"edges") );
827 VJMPERR2( 1==Vio_scanf(sock,"%d", &numEE) );
828
829 /* reset ourself to handle the new manifold dimensions */
830 Vnm_print(0,"Gem_readEdge: theDim=%d, theDimII=%d, numVV=%d, numEE=%d\n",
831 theDim, theDimII, numVV, numEE);
832 Vnm_print(0,"Gem_readEdge: Reseting manifold structures.\n");
833 Gem_reset( thee, theDim, theDimII );
834
835 /* read in the vertices */
836 Vnm_print(0,"Gem_readEdge: Reading ..vertices..");
837 VJMPERR2( 1==Vio_scanf(sock,"%s",tok) );
838 VJMPERR2( !strcmp(tok,"vert") );
839 VJMPERR2( 1==Vio_scanf(sock,"%s",tok) );
840 VJMPERR2( !strcmp(tok,"[") );
841 for (i=0; i<numVV; i++) {
842
843 /* grab the three integer info fields */
844 VJMPERR2( 1==Vio_scanf(sock,"%d",&indexV) );
845 VJMPERR2( 1==Vio_scanf(sock,"%d",&chartV) );
846
847 /* grab the "dimii" coordinates */
848 for (j=0; j<theDimII; j++)
849 VJMPERR2( 1==Vio_scanf(sock,"%le",&x[j]) );
850
851 /* trivially embed any 2D coordinates into 3D */
852 if (theDimII == 2) x[2] = 0.;
853
854 /* create the new vertex */
855 vx = Gem_createAndInitVV(thee);
856 VV_setReality(vx, 0);
857 VV_setDim(vx, theDim);
858 VV_setClass(vx, 0);
859 VV_setType(vx, 0);
860 VV_setId(vx, i);
861 VV_setChart(vx, chartV);
862
863 /* set the vertex coordinates */
864 for (j=0; j<3; j++)
865 VV_setCoord(vx, j, x[j]);
866 }
867 VJMPERR2( 1==Vio_scanf(sock,"%s",tok) );
868 VJMPERR2( !strcmp(tok,"]") );
869
870 /* read in the edges */
871 Vnm_print(0,"..edges..");
872 VJMPERR2( 1==Vio_scanf(sock,"%s",tok) );
873 VJMPERR2( !strcmp(tok,"edge") );
874 VJMPERR2( 1==Vio_scanf(sock,"%s",tok) );
875 VJMPERR2( !strcmp(tok,"[") );
876 for (i=0; i<numEE; i++) {
877
878 /* grab the three integer info fields */
879 VJMPERR2( 1==Vio_scanf(sock,"%d",&indexE) );
880 VJMPERR2( 1==Vio_scanf(sock,"%d",&groupE) );
881 VJMPERR2( 1==Vio_scanf(sock,"%d",&typeE) );
882
883 /* grab the "2" vertex numbers */
884 for (j=0; j<2; j++)
885 VJMPERR2( 1==Vio_scanf(sock,"%d",&se[j]) );
886
887 /* create the new edge */
888 eg = Gem_createAndInitEE(thee);
889 EE_setReality(eg, 0);
890 EE_setDim(eg, theDim);
891 EE_setClass(eg, theDim);
892 EE_setType(eg, typeE);
893 EE_setId(eg, i);
894 EE_setChart(eg, groupE);
895
896 /* set the vertex labels */
897 for (j=0; j<2; j++) {
898 vx = Gem_VV(thee, se[j] );
899 EE_setVertex( eg, j, vx );
900 }
901
902 /* build the ringed vertex datastructure */
903 EE_buildRing(eg);
904 }
905 VJMPERR2( 1==Vio_scanf(sock,"%s",tok) );
906 VJMPERR2( !strcmp(tok,"]") );
907
908 /* get the "mcef_end" tag */
909 VJMPERR2( 1==Vio_scanf(sock,"%s",tok) );
910 VJMPERR2( !strcmp(tok,"mcef_end") );
911 VJMPERR2( 1==Vio_scanf(sock,"%d",&inum) );
912 VJMPERR2( inum == 1 );
913
914 /* finish the i/o */
915 Vnm_print(0,"..done.\n");
916
917 /* close the input file */
918 Vnm_print(0,"Gem_readEdge: Finished reading input\n");
919
920 /* build the simplices from the edges */
921 Gem_buildSfromE(thee);
922
923 /* destroy the edges */
924 Gem_destroyEdges(thee);
925
926 /* create initial edge markings in the simplices */
927 Gem_markEdges(thee);
928
929 /* count v/e/f/s and check the mesh */
930 Gem_countChk(thee);
931
932 /* return with no errors */
933 return 1;
934
935 VERROR1:
936 Vnm_print(2,"Gem_readEdge: Problem with I/O socket.\n");
937 Vnm_print(2,"Gem_readEdge: Reseting manifold structures (bailing out).\n");
938 Gem_reset( thee, 0, 0 );
939 return 0;
940
941 VERROR2:
942 Vnm_print(2,"Gem_readEdge: Format problem with input\n");
943 Vnm_print(2,"Gem_readEdge: Reseting manifold structures (bailing out).\n");
944 Gem_reset( thee, 0, 0 );
945 return 0;
946 }
947
948 /*
949 * ***************************************************************************
950 * Routine: Gem_writeEdge
951 *
952 * Purpose: Write out the edges of a 2- or 3-simplex mesh
953 * in "MCEF" (MC-Edge-Format).
954 *
955 * Author: Michael Holst
956 * ***************************************************************************
957 */
Gem_writeEdge(Gem * thee,Vio * sock)958 VPUBLIC void Gem_writeEdge(Gem *thee, Vio *sock)
959 {
960 int i, j;
961 VV *vx;
962 EE *eg;
963
964 /* some error checking */
965 VJMPERR1( VNULL != sock );
966 VJMPERR2( Gem_numVV(thee) > 0 );
967
968 /* create edges if we don't have them yet */
969 Gem_createEdges(thee);
970
971 /* deal with socket or file setup */
972 Vio_setWhiteChars(sock, MCwhiteChars);
973 Vio_setCommChars(sock, MCcommChars);
974
975 /* write out the header */
976 Vio_printf(sock,"%s",MCcopyRightPCT1);
977 Vio_printf(sock,"%s",MCcopyRightPCT2);
978 Vio_printf(sock,"%s",MCcopyRightPCT3);
979 Vio_printf(sock,"%s",MCcopyRightPCT4);
980 Vio_printf(sock,"\n");
981 Vio_printf(sock,"%s=%d;\n\n","mcef_begin",1);
982
983 Vio_printf(sock," %s=%d; %% intrinsic manifold dimension\n",
984 "dim",Gem_dim(thee));
985
986 Vio_printf(sock," %s=%d; %% imbedding manifold dimension\n",
987 "dimii",Gem_dimII(thee));
988
989 Vio_printf(sock," %s=%d; %% number of vertices\n",
990 "vertices",Gem_numVV(thee));
991
992 Vio_printf(sock," %s=%d; %% number of edges\n",
993 "edges",Gem_numEE(thee));
994
995 /* write out the vertex data */
996 Vio_printf(sock,"\n%s=%s\n","vert","[");
997 Vio_printf(sock,"%%-------- ---- -----------------"
998 " ----------------- -----------------\n");
999 Vio_printf(sock,
1000 "%% Vert-ID Chrt X0-Coordinate X1-Coordinate X2-Coordinate\n");
1001 Vio_printf(sock,"%%-------- ---- -----------------"
1002 " ----------------- -----------------\n");
1003 for (i=0; i<Gem_numVV(thee); i++) {
1004 vx = Gem_VV(thee,i);
1005 Vio_printf(sock,"%-9d %-4d", VV_id(vx), VV_chart(vx) );
1006 for (j=0; j<Gem_dimII(thee); j++)
1007 Vio_printf(sock," %17.10le", VV_coord(vx,j) );
1008 Vio_printf(sock,"\n");
1009 }
1010 Vio_printf(sock,"%s;\n","]");
1011
1012 /* write out the edge data */
1013 Vio_printf(sock,"\n%s=%s\n","edge","[");
1014 Vio_printf(sock,"%%-------- ---- ----"
1015 " ---------------------------------------\n");
1016 Vio_printf(sock,
1017 "%% Edge-ID Grp Mat Vertex-Numbers\n");
1018 Vio_printf(sock,"%%-------- ---- ----"
1019 " ---------------------------------------\n");
1020 for (i=0; i<Gem_numEE(thee); i++) {
1021 eg = Gem_EE(thee,i);
1022 Vio_printf(sock,"%-9d %-4d %-4d %d %d\n",
1023 i, EE_chart(eg), EE_type(eg),
1024 VV_id( EE_vertex(eg,0) ),
1025 VV_id( EE_vertex(eg,1) )
1026 );
1027 }
1028 Vio_printf(sock,"%s;\n","]");
1029
1030 /* finish the file */
1031 Vio_printf(sock,"\n%s=%d;\n\n","mcef_end",1);
1032
1033 /* destroy edges */
1034 Gem_destroyEdges(thee);
1035
1036 /* return with no errors */
1037 Vnm_print(0,"Gem_writeEdge: Finished writing output\n");
1038 return;
1039
1040 VERROR1:
1041 Vnm_print(2,"Gem_writeEdge: Problem with I/O socket.\n");
1042
1043 VERROR2:
1044 Vnm_print(2,"Gem_writeEdge: Detected a problem (bailing out).\n");
1045 return;
1046 }
1047
1048 /*
1049 * ***************************************************************************
1050 * Routine: Gem_buildSfromE
1051 *
1052 * Purpose: Build a vertex-simplex mesh from a vertex-edge mesh.
1053 *
1054 * Notes: Below is an email I sent to R. Bank outlining the idea of the
1055 * algorithm. I wanted to do this completely topologically, using
1056 * no geometry information, so that it would work for abstract
1057 * simplex manifolds.
1058 *
1059 * ---------------------------------------------------------------------------
1060 * Randy,
1061 *
1062 * I realized on the drive that your 2-manifold example is not actually
1063 * a problem after all.
1064 *
1065 * Consider first the planar situation we were worrying about, e.g., take a tet
1066 * and push the fourth vertex down into the plane of the other three vertices,
1067 * and you have the three (good) triangles inside one (bad) triangle.
1068 * As we agreed, any situation like this in the plane is managable (e.g. we can
1069 * detect the "bad" triangle) as long as we known what the boundary edges are
1070 * for the entire mesh.
1071 *
1072 * I conjecture that the "bad" simplices in this case and other 2D cases, as
1073 * well as the 3D case, are precisely those whose faces (edges in 2D) satisfy
1074 * both of the following conditions. (If this is a theorem, then everything
1075 * in this note is rigorous.)
1076 * (a) All interior faces are shared with >1 other simplex
1077 * (b) All boundary faces are shared with >0 other simplex
1078 *
1079 * Consider now your example, e.g. the four surface triangles of a tetrahedron
1080 * as a 2-manifold without boundary. We would like to be able to recover all
1081 * four surface triangles from the six edges, and we don't want to toss out one
1082 * of the triangles as we did above. The key difference is that above, the
1083 * problem triangle either has one or more boundary edges, OR it is imbedded in
1084 * the interior of a mesh, so its interior edges have neighoring triangles.
1085 *
1086 * That doesn't happen for the tet surface example. We would first build all
1087 * possible triangles from triples of edges. According to the above
1088 * "bad simplex" rule, all four of the triangles are "good", since they all
1089 * have exactly one naboring triangle sharing each edge. So we get the correct
1090 * triangulation of the surface of the tet.
1091 *
1092 * I conjecture that the following algorithm will rebuild d-simplices
1093 * (d=2 or d=3) from edges, using only topological information, with no
1094 * geometry information (and no floating point arithmetic at all, for that
1095 * matter). The edge-based input mesh has to satisfy three properties:
1096 *
1097 * 1. The dimension "d" of the mesh is specified (either d=2 or d=3)
1098 * 2. The edge-based mesh was built from an underlying d-simplex mesh
1099 * 3. The boundary edges are marked as such
1100 *
1101 * The two-step simplex reconstruction algorithm is then as follows:
1102 *
1103 * 1. For each vertex "v" do:
1104 * For all vertices connected to "v" by a single edge do:
1105 * Build every possible d-simplex
1106 * EndFor
1107 * EndFor
1108 * 2. For each simplex "s" that was built in step 1 do:
1109 * 3. If d=3, calculate all "face" types of "s" based on edge types
1110 * 4. Remove simplex "s" from list of simplices if (a) AND (b) hold:
1111 * (a) ALL inter faces (edges if d=2) shared by >1 other simplices
1112 * (b) ALL bndry faces (edges if d=2) shared by >0 other simplices
1113 * EndFor
1114 *
1115 * If you can come up with a triangulation of any 2-manifold, with or without
1116 * boundary, with or without holes, etc, which can't be turned into an edge
1117 * mesh and then recovered correctly using the above algorithm, then I'll
1118 * give you a dollar... -mike
1119 *
1120 * Author: Michael Holst
1121 * ***************************************************************************
1122 */
Gem_buildSfromE(Gem * thee)1123 VPUBLIC void Gem_buildSfromE(Gem *thee)
1124 {
1125 int i, j, k, l, m, p, q, r, neg, theDim, goodFace, conformCount;
1126 int groupS, materialS, faceTypeS[4], sv[4];
1127 VV *vx, *v[4], *vtx[100];
1128 EE *eg;
1129 SS *sm, *sm0, *sm1, *sm2;
1130
1131 theDim = Gem_dim(thee);
1132
1133 /* build simplices from edges */
1134 for (i=0; i<Gem_numVV(thee); i++) {
1135 vx = Gem_VV(thee,i);
1136
1137 /* make a list of all vertices which are edge nabors of vx */
1138 neg = 0;
1139 for (eg=VV_firstEE(vx); eg!=VNULL; eg=EE_link(eg,vx)) {
1140 VASSERT( neg < 100 );
1141 vtx[neg] = EE_otherVertex(eg,vx);
1142 neg++;
1143 }
1144
1145 /* 2D: consider all duos of these nabors; do they make an edge? */
1146 if (theDim == 2) {
1147 for (p=0; p<neg; p++) {
1148 for (q=0; q<p; q++) {
1149 if (VNULL!=VV_commonEdge(vtx[p],vtx[q])) {
1150 if (VNULL==VV_commonSimplex3(vx,vtx[p],vtx[q],VNULL)) {
1151
1152 sv[0] = VV_id(vx);
1153 sv[1] = VV_id(vtx[p]);
1154 sv[2] = VV_id(vtx[q]);
1155
1156 /* group, material, and facetypes */
1157 groupS = 0;
1158 materialS = 0;
1159 for (j=0; j<theDim+1; j++) faceTypeS[j]=0;
1160
1161 /* create the new simplex */
1162 sm = Gem_createAndInitSS(thee);
1163 SS_setReality(sm, 0);
1164 SS_setDim(sm, theDim);
1165 SS_setClass(sm, theDim);
1166 SS_setType(sm, materialS);
1167 SS_setId( sm, Gem_numSS(thee)-1 );
1168 SS_setChart(sm, groupS);
1169
1170 /* set simplex face types and vertex labels */
1171 for (j=0; j<theDim+1; j++) {
1172 SS_setFaceType( sm, j, faceTypeS[j] );
1173 if (VBOUNDARY( faceTypeS[j] )) {
1174 (thee->numBF)++;
1175 }
1176 SS_setVertex(sm,j,Gem_VV(thee,sv[j]));
1177 }
1178
1179 /* orientation check */
1180 if (Gem_simplexVolume(thee,sm) < 0) {
1181 SS_setVertex(sm,0,Gem_VV(thee,sv[1]));
1182 SS_setVertex(sm,1,Gem_VV(thee,sv[0]));
1183 }
1184
1185 /* build the ringed vertex datastructure */
1186 SS_buildRing(sm);
1187 }
1188 }
1189 }
1190 }
1191
1192 /* 3D: consider all trios of these nabors; do they make a triangle? */
1193 } else if (theDim == 3) {
1194 for (p=0; p<neg; p++) {
1195 for (q=0; q<p; q++) {
1196 for (r=0; r<q; r++) {
1197 if ( (VNULL != VV_commonEdge(vtx[p],vtx[q]))
1198 && (VNULL != VV_commonEdge(vtx[p],vtx[r]))
1199 && (VNULL != VV_commonEdge(vtx[q],vtx[r])) ) {
1200 if (VNULL == VV_commonSimplex4(vx,
1201 vtx[p],vtx[q],vtx[r],VNULL)) {
1202
1203 sv[0] = VV_id(vx);
1204 sv[1] = VV_id(vtx[p]);
1205 sv[2] = VV_id(vtx[q]);
1206 sv[3] = VV_id(vtx[r]);
1207
1208 /* group, material, and facetypes */
1209 groupS = 0;
1210 materialS = 0;
1211 for (j=0; j<theDim+1; j++) faceTypeS[j]=0;
1212
1213 /* create the new simplex */
1214 sm = Gem_createAndInitSS(thee);
1215 SS_setReality(sm, 0);
1216 SS_setDim(sm, theDim);
1217 SS_setClass(sm, theDim);
1218 SS_setType(sm, materialS);
1219 SS_setId( sm, Gem_numSS(thee)-1 );
1220 SS_setChart(sm, groupS);
1221
1222 /* set simplex face types/vertex labels */
1223 for (j=0; j<theDim+1; j++) {
1224 SS_setFaceType( sm, j, faceTypeS[j] );
1225 if (VBOUNDARY( faceTypeS[j] )) {
1226 (thee->numBF)++;
1227 }
1228 SS_setVertex(sm,j,Gem_VV(thee,sv[j]));
1229 }
1230
1231 /* orientation check */
1232 if (Gem_simplexVolume(thee,sm) < 0) {
1233 SS_setVertex(sm,0,Gem_VV(thee,sv[1]));
1234 SS_setVertex(sm,1,Gem_VV(thee,sv[0]));
1235 }
1236
1237 /* build the ringed vertex datastructure */
1238 SS_buildRing(sm);
1239 }
1240 }
1241 }
1242 }
1243 }
1244 } else { VASSERT(0); }
1245 }
1246
1247 /*
1248 * All that remains is to handle the one degenerate situation.
1249 * This occurs when a group of simplices is contained entirely inside
1250 * another simplex. When this happens, the above simplex construction
1251 * algorithm will incorrectly add the enclosing simplex, which really
1252 * was not part of the original mesh. This of course creates some
1253 * overlapping/nonconforming simplices. It can be detected and corrected
1254 * however by noting that these simplices will be precisely those
1255 * for which ALL of their interior faces have more than one face nabor,
1256 * and all of their boundary faces also have an unexpected nabor.
1257 * The naboring simplices also have some faces with multiple nabors,
1258 * but only the bad simplices have ALL of the faces with bad nabor
1259 * structure. Therefore, we just need to traverse the simplices once
1260 * and remove those which have multiple face nabors on all faces.
1261 *
1262 * NOTE: This requires that we have actually reconstructed the boundary
1263 * face information correctly from the input boundary edge information.
1264 * If the input edge mesh did not contain boundary edge information,
1265 * we have no way to remove these degenerate "extra" simplices.
1266 */
1267 i=0;
1268 while (i < Gem_numSS(thee)) {
1269 sm = Gem_SS(thee,i);
1270
1271 /* get local vertices */
1272 for (j=0; j<Gem_dimVV(thee); j++)
1273 v[j] = SS_vertex(sm,j);
1274
1275 /* check for a bad simplex (they don't have a single good face) */
1276 goodFace=0;
1277 for (j=0; j<Gem_dimVV(thee); j++) {
1278 k=(j+1) % Gem_dimVV(thee);
1279 l=(k+1) % Gem_dimVV(thee);
1280 m=(l+1) % Gem_dimVV(thee);
1281 conformCount = 0;
1282 for (sm0=VV_firstSS(v[k]); sm0!=VNULL;sm0=SS_link(sm0,v[k])) {
1283 for (sm1=VV_firstSS(v[l]);sm1!=VNULL;sm1=SS_link(sm1,v[l])){
1284 if (Gem_dim(thee) == 2) {
1285 if ((sm0!=sm) && (sm0==sm1)) conformCount++;
1286 } else {
1287 for (sm2=VV_firstSS(v[m]); sm2!=VNULL;
1288 sm2=SS_link(sm2,v[m])) {
1289 if ((sm0!=sm) && (sm0==sm1) && (sm0==sm2)) {
1290 conformCount++;
1291 }
1292 }
1293 }
1294 }
1295 }
1296 /* if this face has less that two nabors, it is "good" */
1297 if ( conformCount < 2 ) goodFace++;
1298 }
1299
1300 /* bad simplices have either 0 good faces if they are interior, */
1301 /* or 1 good face (and that one good face is a boundary face). */
1302 if ( (goodFace == 0) || (goodFace == 1) ) {
1303 Vnm_print(0,"Gem_buildSfromE: simplex <%d> is bad..", i);
1304 Vnm_print(0,"..removing (not an error)");
1305
1306 /* get last simplex in list */
1307 sm0 = Gem_lastSS(thee);
1308
1309 /* if sm is not last in the list, exchange it with the last one */
1310 SS_meltRing(sm);
1311 if ( sm != sm0 ) {
1312 SS_meltRing(sm0);
1313
1314 SS_setType(sm, SS_type(sm0));
1315 SS_setChart(sm, SS_chart(sm0));
1316
1317 /* set simplex face types/vertex labels */
1318 for (j=0; j<theDim+1; j++) {
1319 SS_setVertex( sm, j, SS_vertex(sm0,j) );
1320 SS_setFaceType( sm, j, SS_faceType(sm0,j) );
1321 }
1322
1323 /* build the ringed vertex datastructure */
1324 SS_buildRing(sm);
1325 }
1326
1327 /* now remove the last guy */
1328 Gem_destroySS(thee);
1329
1330 Vnm_print(0,"..done.\n");
1331 }
1332
1333 /* iterate */
1334 i++;
1335 }
1336
1337 /* finally, make a boundary */
1338 Gem_makeBnd(thee, 1);
1339
1340 /* return */
1341 return;
1342 }
1343
1344 /*
1345 * ***************************************************************************
1346 * Routine: Gem_writeBrep
1347 *
1348 * Purpose: Write out the boundary edges or boundary faces of a
1349 * 2-simplex or 3-simplex mesh in a "BREP" format for input
1350 * into Barry Joe's 2D/3D Geompak.
1351 *
1352 * Author: Michael Holst
1353 * ***************************************************************************
1354 */
Gem_writeBrep(Gem * thee,Vio * sock)1355 VPUBLIC void Gem_writeBrep(Gem *thee, Vio *sock)
1356 {
1357 if (Gem_dim(thee) == 2) {
1358 Gem_writeBrep2to3(thee, sock);
1359 } else if (Gem_dim(thee) == 3) {
1360 Gem_writeBrep3(thee, sock);
1361 } else { VASSERT(0); }
1362 }
1363
1364 /*
1365 * ***************************************************************************
1366 * Routine: Gem_writeBrep2
1367 *
1368 * Purpose: Write out boundary edges of a 2-simplex mesh
1369 * in a "BREP" format for input into Barry Joe's 2D Geompak.
1370 *
1371 * Author: Michael Holst
1372 * ***************************************************************************
1373 */
Gem_writeBrep2(Gem * thee,Vio * sock)1374 VPUBLIC void Gem_writeBrep2(Gem *thee, Vio *sock)
1375 {
1376 /* some error checking */
1377 VJMPERR1( VNULL != sock );
1378 VJMPERR2( (Gem_numVV(thee) > 0) && (Gem_numSS(thee) > 0) );
1379 VJMPERR2( Gem_dim(thee) == 2 );
1380
1381 Vnm_print(2,"Gem_writeBrep2: this routine is broken.\n");
1382
1383 /* return with no errors */
1384 return;
1385
1386 VERROR1:
1387 Vnm_print(2,"Gem_writeBrep2: Problem with I/O socket.\n");
1388
1389 VERROR2:
1390 Vnm_print(2,"Gem_writeBrep2: Detected a problem (bailing out).\n");
1391 return;
1392 }
1393
1394 /*
1395 * ***************************************************************************
1396 * Routine: Gem_writeBrep3
1397 *
1398 * Purpose: Write out boundary faces of a 3-simplex mesh
1399 * in a "BREP" format for input into Barry Joe's 3D Geompak.
1400 *
1401 * Author: Michael Holst
1402 * ***************************************************************************
1403 */
Gem_writeBrep3(Gem * thee,Vio * sock)1404 VPUBLIC void Gem_writeBrep3(Gem *thee, Vio *sock)
1405 {
1406 int i, j, ii, numVV, numSS, *vTmp;
1407 VV *vx;
1408 SS *sm;
1409
1410 /* some error checking */
1411 VJMPERR1( VNULL != sock );
1412 VJMPERR2( (Gem_numVV(thee) > 0) && (Gem_numSS(thee) > 0) );
1413 VJMPERR2( Gem_dim(thee) == 3 );
1414
1415 /* number of boundary vertices and faces */
1416 numVV = Gem_numBV(thee);
1417 numSS = Gem_numBF(thee);
1418
1419 /* deal with socket or file setup */
1420 Vio_setWhiteChars(sock, MCwhiteChars);
1421 Vio_setCommChars(sock, MCcommChars);
1422
1423 /* write out the header */
1424 Vio_printf(sock,"%s\n","BREP");
1425 Vio_printf(sock,
1426 "1.0e-10 30.0 20.0 40.0 0.05 30.0 0.25 0.5 5 5 0 0.8\n");
1427 Vio_printf(sock,"%d %d 0 1 0 0 0\n\n", numVV, numSS);
1428
1429 /* write out the vertices */
1430 vTmp = Vmem_malloc( thee->vmem, Gem_numVV(thee), sizeof(int) );
1431 ii=0;
1432 for (i=0; i<Gem_numVV(thee); i++) {
1433 vx = Gem_VV(thee,i);
1434 if (VBOUNDARY( VV_type(vx) )) {
1435 vTmp[VV_id(vx)] = ii;
1436 Vio_printf(sock,"%13.6e %13.6e %13.6e\n",
1437 VV_coord(vx,0), VV_coord(vx,1), VV_coord(vx,2) );
1438 ii++;
1439 }
1440 }
1441 Vio_printf(sock,"\n");
1442
1443 /* write out the face starts */
1444 ii=0;
1445 for (i=0; i<Gem_numSS(thee); i++) {
1446 sm = Gem_SS(thee,i);
1447 for (j=0; j<=Gem_dimVV(thee); j++) {
1448 if (VBOUNDARY( SS_faceType(sm,j) )) {
1449 Vio_printf(sock," %4d",3*ii+1);
1450 ii++;
1451 }
1452 }
1453 }
1454 Vio_printf(sock," %4d\n",3*ii+1);
1455
1456 /* write out the face numbers */
1457 ii=0;
1458 for (i=0; i<Gem_numSS(thee); i++) {
1459 sm = Gem_SS(thee,i);
1460 for (j=0; j<=Gem_dimVV(thee); j++) {
1461 if (VBOUNDARY( SS_faceType(sm,j) )) {
1462 Vio_printf(sock," %4d",ii+1);
1463 ii++;
1464 }
1465 }
1466 }
1467 Vio_printf(sock,"\n\n");
1468
1469 /* write out the faces */
1470 for (i=0; i<Gem_numSS(thee); i++) {
1471 sm = Gem_SS(thee,i);
1472 for (j=0; j<=Gem_dimVV(thee); j++) {
1473 if (VBOUNDARY( SS_faceType(sm,j) )) {
1474 Vio_printf(sock,"%4d %4d %4d\n",
1475 vTmp[ VV_id( SS_vertex(sm,vmapF[j][1]) ) ] + 1,
1476 vTmp[ VV_id( SS_vertex(sm,vmapF[j][0]) ) ] + 1,
1477 vTmp[ VV_id( SS_vertex(sm,vmapF[j][2]) ) ] + 1 );
1478 }
1479 }
1480 }
1481 Vio_printf(sock,"\n");
1482
1483 /* write out the face numbers */
1484 Vio_printf(sock," %4d %4d\n",1,numSS+1);
1485 ii=0;
1486 for (i=0; i<Gem_numSS(thee); i++) {
1487 sm = Gem_SS(thee,i);
1488 for (j=0; j<=Gem_dimVV(thee); j++) {
1489 if (VBOUNDARY( SS_faceType(sm,j) )) {
1490 Vio_printf(sock," %4d",ii+1);
1491 ii++;
1492 }
1493 }
1494 }
1495 Vio_printf(sock,"\n");
1496
1497 /* release the temporary storage */
1498 Vmem_free( thee->vmem, Gem_numVV(thee), sizeof(int), (void**)&vTmp );
1499
1500 /* return with no errors */
1501 Vnm_print(0,"Gem_writeBrep3: Finished writing output\n");
1502 return;
1503
1504 VERROR1:
1505 Vnm_print(2,"Gem_writeBrep3: Problem with I/O socket.\n");
1506
1507 VERROR2:
1508 Vnm_print(2,"Gem_writeBrep3: Detected a problem (bailing out).\n");
1509 return;
1510 }
1511
1512 /*
1513 * ***************************************************************************
1514 * Routine: Gem_writeBrep2to3
1515 *
1516 * Purpose: Write out triangles of a 2-manifold as a 3-simplex boundary mesh
1517 * in a "BREP" format for input into Barry Joe's 3D Geompak.
1518 *
1519 * Author: Michael Holst
1520 * ***************************************************************************
1521 */
Gem_writeBrep2to3(Gem * thee,Vio * sock)1522 VPUBLIC void Gem_writeBrep2to3(Gem *thee, Vio *sock)
1523 {
1524 int i, ii, numVV, numSS;
1525 VV *vx;
1526 SS *sm;
1527
1528 /* some error checking */
1529 VJMPERR1( VNULL != sock );
1530 VJMPERR2( (Gem_numVV(thee) > 0) && (Gem_numSS(thee) > 0) );
1531 VJMPERR2( Gem_dim(thee) == 2 );
1532
1533 /* number of boundary vertices and faces */
1534 numVV = Gem_numVV(thee);
1535 numSS = Gem_numSS(thee);
1536
1537 /* deal with socket or file setup */
1538 Vio_setWhiteChars(sock, MCwhiteChars);
1539 Vio_setCommChars(sock, MCcommChars);
1540
1541 /* write out the header */
1542 Vio_printf(sock,"%s\n","BREP");
1543 Vio_printf(sock,
1544 "1.0e-10 30.0 20.0 40.0 0.05 30.0 0.25 0.5 5 5 0 0.8\n");
1545 Vio_printf(sock,"%d %d 0 1 0 0 0\n\n", numVV, numSS);
1546
1547 /* write out the vertices */
1548 for (i=0; i<Gem_numVV(thee); i++) {
1549 vx = Gem_VV(thee,i);
1550 Vio_printf(sock,"%13.6e %13.6e %13.6e\n",
1551 VV_coord(vx,0), VV_coord(vx,1), VV_coord(vx,2) );
1552 }
1553 Vio_printf(sock,"\n");
1554
1555 /* write out the face starts */
1556 ii=0;
1557 for (i=0; i<Gem_numSS(thee); i++) {
1558 sm = Gem_SS(thee,i);
1559 Vio_printf(sock," %4d",3*ii+1);
1560 ii++;
1561 }
1562 Vio_printf(sock," %4d\n",3*ii+1);
1563
1564 /* write out the face numbers */
1565 ii=0;
1566 for (i=0; i<Gem_numSS(thee); i++) {
1567 sm = Gem_SS(thee,i);
1568 Vio_printf(sock," %4d",ii+1);
1569 ii++;
1570 }
1571 Vio_printf(sock,"\n\n");
1572
1573 /* write out the faces */
1574 for (i=0; i<Gem_numSS(thee); i++) {
1575 sm = Gem_SS(thee,i);
1576 Vio_printf(sock,"%4d %4d %4d\n",
1577 VV_id( SS_vertex(sm,0) ) + 1,
1578 VV_id( SS_vertex(sm,1) ) + 1,
1579 VV_id( SS_vertex(sm,2) ) + 1 );
1580 }
1581 Vio_printf(sock,"\n");
1582
1583 /* write out the face numbers */
1584 Vio_printf(sock," %4d %4d\n",1,numSS+1);
1585 ii=0;
1586 for (i=0; i<Gem_numSS(thee); i++) {
1587 sm = Gem_SS(thee,i);
1588 Vio_printf(sock," %4d",ii+1);
1589 ii++;
1590 }
1591 Vio_printf(sock,"\n");
1592
1593 /* return with no errors */
1594 Vnm_print(0,"Gem_writeBrep2to3: Finished writing output\n");
1595 return;
1596
1597 VERROR1:
1598 Vnm_print(2,"Gem_writeBrep2to3: Problem with I/O socket.\n");
1599
1600 VERROR2:
1601 Vnm_print(2,"Gem_writeBrep2to3: Detected a problem (bailing out).\n");
1602 return;
1603 }
1604
1605