1 /* *****************************************************************
2 MESQUITE -- The Mesh Quality Improvement Toolkit
3
4 Copyright 2004 Lawrence Livermore National Laboratory. Under
5 the terms of Contract B545069 with the University of Wisconsin --
6 Madison, Lawrence Livermore National Laboratory retains certain
7 rights in this software.
8
9 This library is free software; you can redistribute it and/or
10 modify it under the terms of the GNU Lesser General Public
11 License as published by the Free Software Foundation; either
12 version 2.1 of the License, or (at your option) any later version.
13
14 This library is distributed in the hope that it will be useful,
15 but WITHOUT ANY WARRANTY; without even the implied warranty of
16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 Lesser General Public License for more details.
18
19 You should have received a copy of the GNU Lesser General Public License
20 (lgpl.txt) along with this library; if not, write to the Free Software
21 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
22
23 kraftche@cae.wisc.edu
24
25 ***************************************************************** */
26
27 #include "MsqError.hpp"
28 #include "TopologyInfo.hpp"
29
30 #include <string.h>
31 #include <assert.h>
32
33 namespace MBMesquite {
34
35 TopologyInfo TopologyInfo::instance;
36
37 const char long_polygon_name[] = "Polygon";
38 const char long_triangle_name[] = "Triangle";
39 const char long_quadrilateral_name[] = "Quadrilateral";
40 const char long_polyhedron_name[] = "Polyhedron";
41 const char long_tetrahedron_name[] = "Tetrahedron";
42 const char long_hexahedron_name[] = "Hexahedron";
43 const char long_prism_name[] = "Prism";
44 const char long_pyramid_name[] = "Pyramd";
45 const char long_septahedron_name[] = "Septahedron";
46 const char short_polygon_name[] = "Polygon";
47 const char short_triangle_name[] = "Tri";
48 const char short_quadrilateral_name[] = "Quad";
49 const char short_polyhedron_name[] = "Polyhedron";
50 const char short_tetrahedron_name[] = "Tet";
51 const char short_hexahedron_name[] = "Hex";
52 const char short_prism_name[] = "Pri";
53 const char short_pyramid_name[] = "Pyr";
54 const char short_septahedron_name[] = "Sept";
55
TopologyInfo()56 TopologyInfo::TopologyInfo()
57 {
58 memset( dimMap, 0, sizeof(dimMap) );
59 memset( adjMap, 0, sizeof(adjMap) );
60 memset( edgeMap, 0, sizeof(edgeMap) );
61 memset( faceMap, 0, sizeof(faceMap) );
62 memset( vertAdjMap, 0, sizeof(vertAdjMap) );
63 memset( shortNames, 0, sizeof(shortNames) );
64 memset( longNames, 0, sizeof(longNames) );
65
66 longNames[POLYGON] = long_polygon_name;
67 longNames[TRIANGLE] = long_triangle_name;
68 longNames[QUADRILATERAL] = long_quadrilateral_name;
69 longNames[POLYHEDRON] = long_polyhedron_name;
70 longNames[TETRAHEDRON] = long_tetrahedron_name;
71 longNames[HEXAHEDRON] = long_hexahedron_name;
72 longNames[PRISM] = long_prism_name;
73 longNames[PYRAMID] = long_pyramid_name;
74 longNames[SEPTAHEDRON] = long_septahedron_name;
75
76 shortNames[POLYGON] = short_polygon_name;
77 shortNames[TRIANGLE] = short_triangle_name;
78 shortNames[QUADRILATERAL] = short_quadrilateral_name;
79 shortNames[POLYHEDRON] = short_polyhedron_name;
80 shortNames[TETRAHEDRON] = short_tetrahedron_name;
81 shortNames[HEXAHEDRON] = short_hexahedron_name;
82 shortNames[PRISM] = short_prism_name;
83 shortNames[PYRAMID] = short_pyramid_name;
84 shortNames[SEPTAHEDRON] = short_septahedron_name;
85
86 dimMap[POLYGON ] = 2;
87 dimMap[TRIANGLE] = 2;
88 dimMap[QUADRILATERAL] = 2;
89 dimMap[POLYHEDRON] = 3;
90 dimMap[TETRAHEDRON] = 3;
91 dimMap[HEXAHEDRON] = 3;
92 dimMap[PRISM] = 3;
93 dimMap[PYRAMID] = 3;
94 dimMap[SEPTAHEDRON] = 3;
95
96 adjMap[TRIANGLE][0] = 3;
97 adjMap[TRIANGLE][1] = 3;
98 adjMap[TRIANGLE][2] = 1;
99 adjMap[TRIANGLE][3] = 0;
100
101 adjMap[QUADRILATERAL][0] = 4;
102 adjMap[QUADRILATERAL][1] = 4;
103 adjMap[QUADRILATERAL][2] = 1;
104 adjMap[QUADRILATERAL][3] = 0;
105
106 adjMap[TETRAHEDRON][0] = 4;
107 adjMap[TETRAHEDRON][1] = 6;
108 adjMap[TETRAHEDRON][2] = 4;
109 adjMap[TETRAHEDRON][3] = 1;
110
111 adjMap[HEXAHEDRON][0] = 8;
112 adjMap[HEXAHEDRON][1] = 12;
113 adjMap[HEXAHEDRON][2] = 6;
114 adjMap[HEXAHEDRON][3] = 1;
115
116 adjMap[PRISM][0] = 6;
117 adjMap[PRISM][1] = 9;
118 adjMap[PRISM][2] = 5;
119 adjMap[PRISM][3] = 1;
120
121 adjMap[PYRAMID][0] = 5;
122 adjMap[PYRAMID][1] = 8;
123 adjMap[PYRAMID][2] = 5;
124 adjMap[PYRAMID][3] = 1;
125
126 adjMap[SEPTAHEDRON][0] = 7;
127 adjMap[SEPTAHEDRON][1] = 11;
128 adjMap[SEPTAHEDRON][2] = 6; /* See description in TSTT mesh interface doc */
129 adjMap[SEPTAHEDRON][3] = 1;
130
131 int side;
132 for (side = 0; side < 3; ++side)
133 {
134 edgeMap[TRIANGLE-FIRST_FACE][side][0] = side;
135 edgeMap[TRIANGLE-FIRST_FACE][side][1] = (side+1)%3;
136 }
137 for (side = 0; side < 4; ++side)
138 {
139 edgeMap[QUADRILATERAL-FIRST_FACE][side][0] = side;
140 edgeMap[QUADRILATERAL-FIRST_FACE][side][1] = (side+1)%4;
141 }
142 for (side = 0; side < 3; ++side)
143 {
144 edgeMap[TETRAHEDRON-FIRST_FACE][side][0] = side;
145 edgeMap[TETRAHEDRON-FIRST_FACE][side][1] = (side+1)%3;
146 }
147 for (side = 3; side < 6; ++side)
148 {
149 edgeMap[TETRAHEDRON-FIRST_FACE][side][0] = side -3 ;
150 edgeMap[TETRAHEDRON-FIRST_FACE][side][1] = 3;
151 }
152 for (side = 0; side < 4; ++side)
153 {
154 edgeMap[HEXAHEDRON-FIRST_FACE][side][0] = side;
155 edgeMap[HEXAHEDRON-FIRST_FACE][side][1] = (side+1)%4;
156 }
157 for (side = 4; side < 8; ++side)
158 {
159 edgeMap[HEXAHEDRON-FIRST_FACE][side][0] = side - 4;
160 edgeMap[HEXAHEDRON-FIRST_FACE][side][1] = side;
161 }
162 for (side = 8; side < 12; ++side)
163 {
164 edgeMap[HEXAHEDRON-FIRST_FACE][side][0] = side - 4;
165 edgeMap[HEXAHEDRON-FIRST_FACE][side][1] = 4+(side+1)%4;
166 }
167 for (side = 0; side < 3; ++side)
168 {
169 edgeMap[PRISM-FIRST_FACE][side][0] = side;
170 edgeMap[PRISM-FIRST_FACE][side][1] = (side+1)%3;
171 }
172 for (side = 3; side < 6; ++side)
173 {
174 edgeMap[PRISM-FIRST_FACE][side][0] = side - 3;
175 edgeMap[PRISM-FIRST_FACE][side][1] = side;
176 }
177 for (side = 6; side < 9; ++side)
178 {
179 edgeMap[PRISM-FIRST_FACE][side][0] = side-3;
180 edgeMap[PRISM-FIRST_FACE][side][1] = 3+(side+1)%3;
181 }
182 for (side = 0; side < 4; ++side)
183 {
184 edgeMap[PYRAMID-FIRST_FACE][side][0] = side;
185 edgeMap[PYRAMID-FIRST_FACE][side][1] = (side+1)%4;
186 }
187 for (side = 4; side < 8; ++side)
188 {
189 edgeMap[PYRAMID-FIRST_FACE][side][0] = side - 4;
190 edgeMap[PYRAMID-FIRST_FACE][side][1] = 4;
191 }
192
193 for (side = 0; side < 3; ++side)
194 {
195 faceMap[TETRAHEDRON-FIRST_VOL][side][0] = 3;
196 faceMap[TETRAHEDRON-FIRST_VOL][side][1] = side;
197 faceMap[TETRAHEDRON-FIRST_VOL][side][2] = (side+1)%3;
198 faceMap[TETRAHEDRON-FIRST_VOL][side][3] = 3;
199 }
200 faceMap[TETRAHEDRON-FIRST_VOL][3][0] = 3;
201 faceMap[TETRAHEDRON-FIRST_VOL][3][1] = 2;
202 faceMap[TETRAHEDRON-FIRST_VOL][3][2] = 1;
203 faceMap[TETRAHEDRON-FIRST_VOL][3][3] = 0;
204
205 for (side = 0; side < 4; ++side)
206 {
207 faceMap[HEXAHEDRON-FIRST_VOL][side][0] = 4;
208 faceMap[HEXAHEDRON-FIRST_VOL][side][1] = side;
209 faceMap[HEXAHEDRON-FIRST_VOL][side][2] = (side+1)%4;
210 faceMap[HEXAHEDRON-FIRST_VOL][side][3] = 4+(side+1)%4;
211 faceMap[HEXAHEDRON-FIRST_VOL][side][4] = side + 4;
212 }
213 faceMap[HEXAHEDRON-FIRST_VOL][4][0] = 4;
214 faceMap[HEXAHEDRON-FIRST_VOL][4][1] = 3;
215 faceMap[HEXAHEDRON-FIRST_VOL][4][2] = 2;
216 faceMap[HEXAHEDRON-FIRST_VOL][4][3] = 1;
217 faceMap[HEXAHEDRON-FIRST_VOL][4][4] = 0;
218 faceMap[HEXAHEDRON-FIRST_VOL][5][0] = 4;
219 faceMap[HEXAHEDRON-FIRST_VOL][5][1] = 4;
220 faceMap[HEXAHEDRON-FIRST_VOL][5][2] = 5;
221 faceMap[HEXAHEDRON-FIRST_VOL][5][3] = 6;
222 faceMap[HEXAHEDRON-FIRST_VOL][5][4] = 7;
223
224 for (side = 0; side < 4; ++side)
225 {
226 faceMap[PYRAMID-FIRST_VOL][side][0] = 3;
227 faceMap[PYRAMID-FIRST_VOL][side][1] = side;
228 faceMap[PYRAMID-FIRST_VOL][side][2] = (side+1)%4;
229 faceMap[PYRAMID-FIRST_VOL][side][3] = 4;
230 }
231 faceMap[PYRAMID-FIRST_VOL][4][0] = 4;
232 faceMap[PYRAMID-FIRST_VOL][4][1] = 3;
233 faceMap[PYRAMID-FIRST_VOL][4][2] = 2;
234 faceMap[PYRAMID-FIRST_VOL][4][3] = 1;
235 faceMap[PYRAMID-FIRST_VOL][4][4] = 0;
236
237 for (side = 0; side < 3; ++side)
238 {
239 faceMap[PRISM-FIRST_VOL][side][0] = 4;
240 faceMap[PRISM-FIRST_VOL][side][1] = side;
241 faceMap[PRISM-FIRST_VOL][side][2] = (side+1)%3;
242 faceMap[PRISM-FIRST_VOL][side][3] = 3+(side+1)%3;
243 faceMap[PRISM-FIRST_VOL][side][4] = side + 3;
244 }
245 faceMap[PRISM-FIRST_VOL][3][0] = 3;
246 faceMap[PRISM-FIRST_VOL][3][1] = 2;
247 faceMap[PRISM-FIRST_VOL][3][2] = 1;
248 faceMap[PRISM-FIRST_VOL][3][3] = 0;
249 faceMap[PRISM-FIRST_VOL][4][0] = 3;
250 faceMap[PRISM-FIRST_VOL][4][1] = 3;
251 faceMap[PRISM-FIRST_VOL][4][2] = 4;
252 faceMap[PRISM-FIRST_VOL][4][3] = 5;
253
254 int i;
255 for (i = 0; i < 3; ++i)
256 {
257 vertAdjMap[TRIANGLE-FIRST_FACE][i][0] = 2;
258 vertAdjMap[TRIANGLE-FIRST_FACE][i][1] = (i+1)%3;
259 vertAdjMap[TRIANGLE-FIRST_FACE][i][2] = (i+2)%3;
260 }
261
262 for (i = 0; i < 4; ++i)
263 {
264 vertAdjMap[QUADRILATERAL-FIRST_FACE][i][0] = 2;
265 vertAdjMap[QUADRILATERAL-FIRST_FACE][i][1] = (i+1)%4;
266 vertAdjMap[QUADRILATERAL-FIRST_FACE][i][2] = (i+3)%4;
267 }
268
269
270 unsigned tet_corner_data[] = { 1, 2, 3,
271 0, 3, 2,
272 3, 0, 1,
273 2, 1, 0 };
274 for (i = 0; i < 4; ++i)
275 {
276 vertAdjMap[TETRAHEDRON-FIRST_FACE][i][0] = 3;
277 for (unsigned j = 0; j < 3; ++j)
278 vertAdjMap[TETRAHEDRON-FIRST_FACE][i][j+1] = tet_corner_data[3*i+j];
279 }
280
281 for (i = 0; i < 4; ++i)
282 {
283 vertAdjMap[PYRAMID-FIRST_FACE][i][0] = 3;
284 vertAdjMap[PYRAMID-FIRST_FACE][i][1] = (i+1)%4;
285 vertAdjMap[PYRAMID-FIRST_FACE][i][2] = (i+3)%4;
286 vertAdjMap[PYRAMID-FIRST_FACE][i][3] = 4;
287 }
288 vertAdjMap[PYRAMID-FIRST_FACE][4][0] = 4;
289 for (i = 0; i < 4; i++)
290 vertAdjMap[PYRAMID-FIRST_FACE][4][i+1] = 3 - i;
291
292 for (i = 0; i < 4; ++i)
293 {
294 vertAdjMap[HEXAHEDRON-FIRST_FACE][i][0] = 3;
295 vertAdjMap[HEXAHEDRON-FIRST_FACE][i][1] = (i+1)%4;
296 vertAdjMap[HEXAHEDRON-FIRST_FACE][i][2] = (i+3)%4;
297 vertAdjMap[HEXAHEDRON-FIRST_FACE][i][3] = i+4;
298 }
299 for (i = 4; i < 8; ++i)
300 {
301 vertAdjMap[HEXAHEDRON-FIRST_FACE][i][0] = 3;
302 vertAdjMap[HEXAHEDRON-FIRST_FACE][i][1] = (i+3)%4+4;
303 vertAdjMap[HEXAHEDRON-FIRST_FACE][i][2] = (i+1)%4+4;
304 vertAdjMap[HEXAHEDRON-FIRST_FACE][i][3] = i-4;
305 }
306
307 for (i = 0; i < 3; ++i)
308 {
309 vertAdjMap[PRISM-FIRST_FACE][i][0] = 3;
310 vertAdjMap[PRISM-FIRST_FACE][i][1] = (i+1)%3;
311 vertAdjMap[PRISM-FIRST_FACE][i][2] = (i+2)%3;
312 vertAdjMap[PRISM-FIRST_FACE][i][3] = i+3;
313 }
314 for (i = 3; i < 6; ++i)
315 {
316 vertAdjMap[PRISM-FIRST_FACE][i][0] = 3;
317 vertAdjMap[PRISM-FIRST_FACE][i][1] = (i+2)%3+3;
318 vertAdjMap[PRISM-FIRST_FACE][i][2] = (i+1)%3+3;
319 vertAdjMap[PRISM-FIRST_FACE][i][3] = i-3;
320 }
321
322 // Build reverse vertex-vertex adjacency index map
323 const EntityTopology types[] = { TRIANGLE,
324 QUADRILATERAL,
325 TETRAHEDRON,
326 PYRAMID,
327 PRISM,
328 HEXAHEDRON };
329 const int num_types = sizeof(types)/sizeof(types[0]);
330 for (i = 0; i < num_types; ++i)
331 {
332 const unsigned num_vert = corners( types[i] );
333 for (unsigned v = 0; v < num_vert; ++v)
334 {
335 unsigned num_v_adj;
336 const unsigned* v_adj = adjacent_vertices( types[i], v, num_v_adj );
337 unsigned* reverse = revVertAdjIdx[types[i]-FIRST_FACE][v];
338 reverse[0] = num_v_adj;
339
340 for (unsigned j = 0; j < num_v_adj; ++j)
341 {
342 unsigned num_j_adj, k;
343 const unsigned* j_adj = adjacent_vertices( types[i], v_adj[j], num_j_adj );
344 for (k = 0; k < num_j_adj && j_adj[k] != v; ++k);
345 assert( k < num_j_adj ); // If this fails, vertAdjMap is corrupt!
346 reverse[j+1] = k;
347 }
348 }
349 }
350 }
351
higher_order(EntityTopology topo,unsigned num_nodes,bool & midedge,bool & midface,bool & midvol,MsqError & err)352 void TopologyInfo::higher_order( EntityTopology topo,
353 unsigned num_nodes,
354 bool& midedge,
355 bool& midface,
356 bool& midvol,
357 MsqError& err )
358 {
359 int ho = higher_order( topo, num_nodes, err );
360 midedge = (bool)( (ho & (1<<1)) >> 1);
361 midface = (bool)( (ho & (1<<2)) >> 2);
362 midvol = (bool)( (ho & (1<<3)) >> 3);
363 }
364
365
higher_order(EntityTopology topo,unsigned num_nodes,MsqError & err)366 int TopologyInfo::higher_order( EntityTopology topo,
367 unsigned num_nodes,
368 MsqError& err )
369 {
370 int result = 0;
371 if (topo == POLYGON) // polygons currently do not have higher order elements
372 return 0;
373
374 if (topo >= MIXED || num_nodes < instance.adjMap[topo][0])
375 {
376 MSQ_SETERR(err)("Invalid element topology", MsqError::INVALID_ARG);
377 return 0;
378 }
379
380 unsigned dim = instance.dimMap[topo];
381 assert( num_nodes >= instance.adjMap[topo][0] );
382 unsigned nodes = num_nodes - instance.adjMap[topo][0];
383 unsigned edges = instance.adjMap[topo][1];
384 unsigned faces = instance.adjMap[topo][2];
385 if (edges && nodes >= edges)
386 {
387 nodes -= edges;
388 result |= 1<<1;
389 }
390 if (faces && nodes >= faces)
391 {
392 nodes -= faces;
393 result |= 1<<2;
394 }
395 if (1 == nodes)
396 {
397 if (2 == dim)
398 {
399 nodes -= 1;
400 result |= 1<<2;
401 }
402 else if(3 == dim)
403 {
404 nodes -= 1;
405 result |= 1<<3;
406 }
407 }
408
409 if (nodes)
410 {
411 MSQ_SETERR(err)("Invalid element topology", MsqError::INVALID_STATE);
412 }
413
414 return result;
415 }
416
higher_order_from_side(EntityTopology topo,unsigned num_nodes,unsigned side_dimension,unsigned side_number,MsqError & err)417 int TopologyInfo::higher_order_from_side( EntityTopology topo,
418 unsigned num_nodes,
419 unsigned side_dimension,
420 unsigned side_number,
421 MsqError& err )
422 {
423 bool mids[4] = { true };
424 higher_order( topo, num_nodes, mids[1], mids[2], mids[3], err );
425 MSQ_ERRZERO(err);
426
427 if (side_dimension > dimension(topo) ||
428 side_number > adjacent(topo, side_dimension)) {
429 MSQ_SETERR(err)(MsqError::INVALID_ARG,"Invalid side number: %u\n", side_number );
430 return 0;
431 }
432
433 if (!mids[side_dimension])
434 return -1;
435
436 int result = side_number;
437 switch (side_dimension) {
438 case 3: if (mids[2]) result += faces(topo);
439 case 2: if (mids[1]) result += edges(topo);
440 case 1: result += corners(topo);
441 case 0: break;
442 default:
443 MSQ_SETERR(err)(MsqError::INVALID_ARG,"Invalid dimension: %u\n", side_dimension );
444 return 0;
445 }
446 return result;
447 }
448
side_from_higher_order(EntityTopology topo,unsigned num_nodes,unsigned node_number,unsigned & side_dim_out,unsigned & side_num_out,MsqError & err)449 void TopologyInfo::side_from_higher_order( EntityTopology topo,
450 unsigned num_nodes,
451 unsigned node_number,
452 unsigned& side_dim_out,
453 unsigned& side_num_out,
454 MsqError& err )
455 {
456 bool midedge, midface, midvol;
457 higher_order( topo, num_nodes, midedge, midface, midvol, err );
458 MSQ_ERRRTN(err);
459 side_num_out = node_number;
460
461 if (side_num_out < corners(topo)) {
462 side_dim_out = 0;
463 return;
464 }
465 side_num_out -= corners(topo);
466
467 if (midedge) {
468 if (side_num_out < edges(topo)) {
469 side_dim_out = 1;
470 return;
471 }
472 side_num_out -= edges(topo);
473 }
474
475 if (midface) {
476 if (side_num_out < faces(topo)) {
477 side_dim_out = 2;
478 return;
479 }
480 side_num_out -= faces(topo);
481 }
482
483 if (midvol && side_num_out == 0) {
484 side_dim_out = 3;
485 return;
486 }
487
488 MSQ_SETERR(err)(MsqError::INVALID_ARG,"Invalid node index\n");
489 }
490
edge_vertices(EntityTopology topo,unsigned edge,MsqError & err)491 const unsigned* TopologyInfo::edge_vertices( EntityTopology topo,
492 unsigned edge,
493 MsqError& err)
494 {
495 if (topo < (EntityTopology)FIRST_FACE ||
496 topo > (EntityTopology)LAST_VOL ||
497 edge >= edges( topo ) )
498 {
499 MSQ_SETERR(err)(MsqError::INVALID_ARG);
500 topo = (EntityTopology)FIRST_FACE;
501 edge = 0;
502 }
503
504 return instance.edgeMap[topo-FIRST_FACE][edge];
505 }
506
edge_vertices(EntityTopology topo,unsigned edge)507 const unsigned* TopologyInfo::edge_vertices( EntityTopology topo,
508 unsigned edge )
509 {
510 if (topo < (EntityTopology)FIRST_FACE ||
511 topo > (EntityTopology)LAST_VOL ||
512 edge >= edges( topo ) )
513 {
514 return 0;
515 }
516 return instance.edgeMap[topo-FIRST_FACE][edge];
517 }
518
face_vertices(EntityTopology topo,unsigned face,unsigned & length,MsqError & err)519 const unsigned* TopologyInfo::face_vertices( EntityTopology topo,
520 unsigned face,
521 unsigned& length,
522 MsqError& err )
523 {
524 if (topo < (EntityTopology)FIRST_VOL ||
525 topo > (EntityTopology)LAST_VOL ||
526 face >= faces( topo ) )
527 {
528 MSQ_SETERR(err)(MsqError::INVALID_ARG);
529 topo = (EntityTopology)FIRST_VOL;
530 face = 0;
531 }
532
533 length = instance.faceMap[topo-FIRST_VOL][face][0];
534 return instance.faceMap[topo-FIRST_VOL][face] + 1;
535 }
face_vertices(EntityTopology topo,unsigned face,unsigned & length)536 const unsigned* TopologyInfo::face_vertices( EntityTopology topo,
537 unsigned face,
538 unsigned& length )
539 {
540 if (topo < (EntityTopology)FIRST_VOL ||
541 topo > (EntityTopology)LAST_VOL ||
542 face >= faces( topo ) )
543 {
544 return 0;
545 }
546
547 length = instance.faceMap[topo-FIRST_VOL][face][0];
548 return instance.faceMap[topo-FIRST_VOL][face] + 1;
549 }
550
551
552
553
side_vertices(EntityTopology topo,unsigned dim,unsigned side,unsigned & count_out,MsqError & err)554 const unsigned* TopologyInfo::side_vertices( EntityTopology topo,
555 unsigned dim,
556 unsigned side,
557 unsigned& count_out,
558 MsqError& err )
559 {
560 static const unsigned all[] = { 0, 1, 2, 3, 4, 5, 6, 7 };
561 const unsigned* result;
562
563 if (dim != 0 && dim == dimension(topo))
564 {
565 count_out = corners( topo );
566 result = all;
567 }
568 else if (dim == 1)
569 {
570 count_out = 2;
571 result = edge_vertices( topo, side, err );
572 }
573 else if( dim == 2)
574 {
575 result = face_vertices( topo, side, count_out, err );
576 }
577 else
578 {
579 MSQ_SETERR(err)(MsqError::INVALID_ARG);
580 count_out = 0;
581 result = 0;
582 }
583 return result;
584 }
side_vertices(EntityTopology topo,unsigned dim,unsigned side,unsigned & count_out)585 const unsigned* TopologyInfo::side_vertices( EntityTopology topo,
586 unsigned dim,
587 unsigned side,
588 unsigned& count_out )
589 {
590 static const unsigned all[] = { 0, 1, 2, 3, 4, 5, 6, 7 };
591 const unsigned* result;
592
593 if (dim != 0 && dim == dimension(topo))
594 {
595 count_out = corners( topo );
596 result = all;
597 }
598 else if (dim == 1)
599 {
600 count_out = 2;
601 result = edge_vertices( topo, side );
602 }
603 else if( dim == 2)
604 {
605 result = face_vertices( topo, side, count_out );
606 }
607 else
608 {
609 result = 0;
610 }
611 return result;
612 }
613
614
615
side_number(EntityTopology topo,unsigned num_nodes,unsigned node_index,unsigned & side_dim_out,unsigned & side_num_out,MsqError & err)616 void TopologyInfo::side_number( EntityTopology topo,
617 unsigned num_nodes,
618 unsigned node_index,
619 unsigned& side_dim_out,
620 unsigned& side_num_out,
621 MsqError& err )
622 {
623 if (topo >= (EntityTopology)MIXED || num_nodes < instance.adjMap[topo][0])
624 {
625 MSQ_SETERR(err)("Invalid element topology", MsqError::INVALID_ARG);
626 return;
627 }
628
629 unsigned nodes = instance.adjMap[topo][0];
630 unsigned edges = instance.adjMap[topo][1];
631 unsigned faces = instance.adjMap[topo][2];
632 side_num_out = node_index;
633
634 if (side_num_out < nodes)
635 {
636 side_dim_out = 0;
637 return;
638 }
639 num_nodes -= nodes;
640 side_num_out -= nodes;
641
642 if (edges && num_nodes >= edges)
643 {
644 if (side_num_out < edges)
645 {
646 side_dim_out = 1;
647 return;
648 }
649 num_nodes -= edges;
650 side_num_out -= edges;
651 }
652 if (faces && num_nodes >= faces)
653 {
654 if (side_num_out < faces)
655 {
656 side_dim_out = 2;
657 return;
658 }
659 num_nodes -= faces;
660 side_num_out -= faces;
661 }
662 if (side_num_out == 0)
663 {
664 side_dim_out = instance.dimMap[topo];
665 side_num_out = 0;
666 return;
667 }
668
669 MSQ_SETERR(err)(MsqError::INVALID_ARG);
670 }
671
672
673
adjacent_vertices(EntityTopology topo,unsigned index,unsigned & num_adj_out)674 const unsigned* TopologyInfo::adjacent_vertices( EntityTopology topo,
675 unsigned index,
676 unsigned& num_adj_out )
677 {
678 const unsigned count = corners( topo );
679 if (!count || index >= count)
680 {
681 num_adj_out = 0;
682 return 0;
683 }
684
685 const unsigned* vect = instance.vertAdjMap[topo-FIRST_FACE][index];
686 num_adj_out = vect[0];
687 return vect + 1;
688 }
689
reverse_vertex_adjacency_offsets(EntityTopology topo,unsigned index,unsigned & num_adj_out)690 const unsigned* TopologyInfo::reverse_vertex_adjacency_offsets(
691 EntityTopology topo,
692 unsigned index,
693 unsigned& num_adj_out )
694 {
695 const unsigned count = corners( topo );
696 if (!count || index >= count)
697 {
698 num_adj_out = 0;
699 return 0;
700 }
701
702 const unsigned* vect = instance.revVertAdjIdx[topo-FIRST_FACE][index];
703 num_adj_out = vect[0];
704 return vect + 1;
705 }
706
compare_sides(const size_t * verts1,EntityTopology type1,unsigned side1,const size_t * verts2,EntityTopology type2,unsigned side2,unsigned side_dim,MsqError & err)707 bool TopologyInfo::compare_sides( const size_t* verts1,
708 EntityTopology type1,
709 unsigned side1,
710 const size_t* verts2,
711 EntityTopology type2,
712 unsigned side2,
713 unsigned side_dim,
714 MsqError& err )
715 {
716 const unsigned *conn1, *conn2;
717 unsigned len1, len2;
718
719 conn1 = side_vertices( type1, side_dim, side1, len1, err );
720 MSQ_ERRZERO(err);
721 conn2 = side_vertices( type2, side_dim, side2, len2, err );
722 MSQ_ERRZERO(err);
723
724 // obviously not the same if different number of vertices
725 // (triangular face cannot match quadrilateral face)
726 if (len1 != len2)
727 return false;
728
729 // Find location (i) in vertices of side of second element
730 // that matches the first vertex in the side of the first
731 // element.
732 unsigned i, j;
733 for (i = 0; i < len2; ++i)
734 if (verts1[conn1[0]] == verts2[conn2[i]])
735 break;
736 // If not found, then no match
737 if (i == len2)
738 return false;
739
740 // Try comparing side connectivity in forward order
741 for (j = 1; j < len1; ++j)
742 if (verts1[conn1[j]] != verts2[conn2[(i+j)%len2]])
743 break;
744 // If they match, we're done
745 if (j == len1)
746 return true;
747
748 // Try comparing in reverse order
749 for (j = 1; j < len1; ++j)
750 if (verts1[conn1[j]] != verts2[conn2[(i+len2-j)%len2]])
751 return false;
752 // If here, matched in reverse order
753 return true;
754 }
755
find_edge(EntityTopology topo,const unsigned * side_vertices,bool & reversed_out,MsqError & err)756 unsigned TopologyInfo::find_edge( EntityTopology topo,
757 const unsigned* side_vertices,
758 bool& reversed_out,
759 MsqError& err )
760 {
761 if (dimension(topo) <= 1) {
762 MSQ_SETERR(err)(MsqError::INVALID_ARG,"Invalid element dimension");
763 return (unsigned)-1;
764 }
765
766 for (unsigned i = 0; i < edges(topo); ++i) {
767 const unsigned* edge = edge_vertices( topo, i, err );
768 MSQ_ERRZERO(err);
769
770 if (edge[0] == side_vertices[0] &&
771 edge[1] == side_vertices[1]) {
772 reversed_out = false;
773 return i;
774 }
775
776 if (edge[0] == side_vertices[1] &&
777 edge[1] == side_vertices[0]) {
778 reversed_out = true;
779 return i;
780 }
781 }
782
783 MSQ_SETERR(err)(MsqError::INVALID_ARG,"No such edge");
784 return (unsigned)-1;
785 }
786
find_face(EntityTopology topo,const unsigned * side_vertices,unsigned num_vertices,bool & reversed_out,MsqError & err)787 unsigned TopologyInfo::find_face( EntityTopology topo,
788 const unsigned* side_vertices,
789 unsigned num_vertices,
790 bool& reversed_out,
791 MsqError& err )
792 {
793 if (dimension(topo) <= 2) {
794 MSQ_SETERR(err)(MsqError::INVALID_ARG,"Invalid element dimension");
795 return (unsigned)-1;
796 }
797
798 for (unsigned i = 0; i < faces(topo); ++i) {
799 unsigned j, n, offset;
800 const unsigned* face = face_vertices( topo, i, n, err );
801 MSQ_ERRZERO(err);
802 if (n != num_vertices)
803 continue;
804
805 for (offset = 0; offset < num_vertices; ++offset)
806 if (face[offset] == side_vertices[0])
807 break;
808 if (offset == num_vertices)
809 continue;
810
811 for (j = 1; j < num_vertices; ++j)
812 if (side_vertices[j] != face[(offset + j)%num_vertices])
813 break;
814 if (j == num_vertices) {
815 reversed_out = false;
816 return i;
817 }
818
819 for (j = 1; j < num_vertices; ++j)
820 if (side_vertices[j] != face[(offset + num_vertices - j)%num_vertices])
821 break;
822 if (j == num_vertices) {
823 reversed_out = true;
824 return i;
825 }
826 }
827
828 MSQ_SETERR(err)(MsqError::INVALID_ARG,"No such face");
829 return (unsigned)-1;
830 }
831
832
find_side(EntityTopology topo,const unsigned * side_vertices,unsigned num_vertices,unsigned & dimension_out,unsigned & number_out,bool & reversed_out,MsqError & err)833 void TopologyInfo::find_side( EntityTopology topo,
834 const unsigned* side_vertices,
835 unsigned num_vertices,
836 unsigned& dimension_out,
837 unsigned& number_out,
838 bool& reversed_out,
839 MsqError& err )
840 {
841 switch (num_vertices) {
842 case 1:
843 dimension_out = 0;
844 number_out = *side_vertices;
845 reversed_out = false;
846 if (*side_vertices >= corners(topo))
847 MSQ_SETERR(err)(MsqError::INVALID_ARG,"Invalid corner number: %u\n", *side_vertices);
848 break;
849 case 2:
850 dimension_out = 1;
851 number_out = find_edge( topo, side_vertices, reversed_out, err );
852 MSQ_CHKERR(err);
853 break;
854 case 3:
855 case 4:
856 dimension_out = 2;
857 number_out = find_face( topo, side_vertices, num_vertices, reversed_out, err );
858 MSQ_CHKERR(err);
859 break;
860 default:
861 MSQ_SETERR(err)(MsqError::UNSUPPORTED_ELEMENT, "Invalid number of side vertices: %u\n", num_vertices );
862 break;
863 }
864 }
865
866
867
868
869 } //namepsace Mesquite
870