1 /**
2 * MOAB, a Mesh-Oriented datABase, is a software component for creating,
3 * storing and accessing finite element mesh data.
4 *
5 * Copyright 2004 Sandia Corporation. Under the terms of Contract
6 * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
7 * retains certain 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 */
15
16 #include "moab/CN.hpp"
17 #include "MBCNArrays.hpp"
18 #include "MBCN.h"
19 #include <assert.h>
20 #include <string.h>
21 #include <iterator>
22
23 namespace moab {
24
25 const char *CN::entityTypeNames[] = {
26 "Vertex",
27 "Edge",
28 "Tri",
29 "Quad",
30 "Polygon",
31 "Tet",
32 "Pyramid",
33 "Prism",
34 "Knife",
35 "Hex",
36 "Polyhedron",
37 "EntitySet",
38 "MaxType"
39 };
40
41 short int CN::numberBasis = 0;
42
43 short int CN::permuteVec[MBMAXTYPE][3][MAX_SUB_ENTITIES + 1];
44 short int CN::revPermuteVec[MBMAXTYPE][3][MAX_SUB_ENTITIES + 1];
45
46 const DimensionPair CN::TypeDimensionMap[] =
47 {
48 DimensionPair(MBVERTEX, MBVERTEX),
49 DimensionPair(MBEDGE, MBEDGE),
50 DimensionPair(MBTRI, MBPOLYGON),
51 DimensionPair(MBTET, MBPOLYHEDRON),
52 DimensionPair(MBENTITYSET, MBENTITYSET),
53 DimensionPair(MBMAXTYPE, MBMAXTYPE)
54 };
55
56 short CN::increasingInts[] = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9,
57 10,11,12,13,14,15,16,17,18,19,
58 20,21,22,23,24,25,26,27,28,29,
59 30,31,32,33,34,35,36,37,38,39 };
60
61 //! set the basis of the numbering system; may or may not do things besides setting the
62 //! member variable
SetBasis(const int in_basis)63 void CN::SetBasis(const int in_basis)
64 {
65 numberBasis = in_basis;
66 }
67
68 //! return a type for the given name
EntityTypeFromName(const char * name)69 EntityType CN::EntityTypeFromName(const char *name)
70 {
71 for (EntityType i = MBVERTEX; i < MBMAXTYPE; i++) {
72 if (0 == strcmp(name, entityTypeNames[i]))
73 return i;
74 }
75
76 return MBMAXTYPE;
77 }
78
SubEntityNodeIndices(const EntityType this_topo,const int num_nodes,const int sub_dimension,const int sub_index,EntityType & subentity_topo,int & num_sub_entity_nodes,int sub_entity_conn[])79 void CN::SubEntityNodeIndices( const EntityType this_topo,
80 const int num_nodes,
81 const int sub_dimension,
82 const int sub_index,
83 EntityType& subentity_topo,
84 int& num_sub_entity_nodes,
85 int sub_entity_conn[] )
86 {
87 // If asked for a node, the special case...
88 if (sub_dimension == 0) {
89 assert( sub_index < num_nodes );
90 subentity_topo = MBVERTEX;
91 num_sub_entity_nodes = 1;
92 sub_entity_conn[0] = sub_index;
93 return;
94 }
95
96 const int ho_bits = HasMidNodes( this_topo, num_nodes );
97 subentity_topo = SubEntityType(this_topo, sub_dimension, sub_index);
98 num_sub_entity_nodes = VerticesPerEntity(subentity_topo);
99 const short* corners = mConnectivityMap[this_topo][sub_dimension-1].conn[sub_index];
100 std::copy( corners, corners+num_sub_entity_nodes, sub_entity_conn );
101
102 int sub_sub_corners[MAX_SUB_ENTITY_VERTICES];
103 int side, sense, offset;
104 for (int dim = 1; dim <= sub_dimension; ++dim) {
105 if (!(ho_bits & (1<<dim)))
106 continue;
107
108 const short num_mid = NumSubEntities( subentity_topo, dim );
109 for (int i = 0; i < num_mid; ++i) {
110 const EntityType sub_sub_topo = SubEntityType( subentity_topo, dim, i );
111 const int sub_sub_num_vert = VerticesPerEntity( sub_sub_topo );
112 SubEntityVertexIndices( subentity_topo, dim, i, sub_sub_corners );
113
114 for (int j = 0; j < sub_sub_num_vert; ++j)
115 sub_sub_corners[j] = corners[sub_sub_corners[j]];
116 SideNumber( this_topo, sub_sub_corners, sub_sub_num_vert, dim, side, sense, offset );
117 sub_entity_conn[num_sub_entity_nodes++] = HONodeIndex( this_topo, num_nodes, dim, side );
118 }
119 }
120 }
121
122
123 //! return the vertices of the specified sub entity
124 //! \param parent_conn Connectivity of parent entity
125 //! \param parent_type Entity type of parent entity
126 //! \param sub_dimension Dimension of sub-entity being queried
127 //! \param sub_index Index of sub-entity being queried
128 //! \param sub_entity_conn Connectivity of sub-entity, based on parent_conn and canonical
129 //! ordering for parent_type
130 //! \param num_sub_vertices Number of vertices in sub-entity
SubEntityConn(const void * parent_conn,const EntityType parent_type,const int sub_dimension,const int sub_index,void * sub_entity_conn,int & num_sub_vertices)131 void CN::SubEntityConn(const void *parent_conn, const EntityType parent_type,
132 const int sub_dimension,
133 const int sub_index,
134 void *sub_entity_conn, int &num_sub_vertices)
135 {
136 static int sub_indices[MAX_SUB_ENTITY_VERTICES];
137
138 SubEntityVertexIndices(parent_type, sub_dimension, sub_index, sub_indices);
139
140 num_sub_vertices = VerticesPerEntity(SubEntityType(parent_type, sub_dimension, sub_index));
141 void **parent_conn_ptr = static_cast<void **>(const_cast<void *>(parent_conn));
142 void **sub_conn_ptr = static_cast<void **>(sub_entity_conn);
143 for (int i = 0; i < num_sub_vertices; i++)
144 sub_conn_ptr[i] = parent_conn_ptr[sub_indices[i]];
145 }
146
147 //! given an entity and a target dimension & side number, get that entity
AdjacentSubEntities(const EntityType this_type,const int * source_indices,const int num_source_indices,const int source_dim,const int target_dim,std::vector<int> & index_list,const int operation_type)148 short int CN::AdjacentSubEntities(const EntityType this_type,
149 const int *source_indices,
150 const int num_source_indices,
151 const int source_dim,
152 const int target_dim,
153 std::vector<int> &index_list,
154 const int operation_type)
155 {
156 // first get all the vertex indices
157 std::vector<int> tmp_indices;
158 const int* it1 = source_indices;
159
160 assert(source_dim <= 3 &&
161 target_dim >= 0 && target_dim <= 3 &&
162 // make sure we're not stepping off the end of the array;
163 ((source_dim > 0 &&
164 *it1 < mConnectivityMap[this_type][source_dim-1].num_sub_elements) ||
165 (source_dim == 0 &&
166 *it1 < mConnectivityMap[this_type][Dimension(this_type)-1].num_corners_per_sub_element[0])) &&
167 *it1 >= 0);
168
169
170 #define MUC CN::mUpConnMap[this_type][source_dim][target_dim]
171
172 // if we're looking for the vertices of a single side, return them in
173 // the canonical ordering; otherwise, return them in sorted order
174 if (num_source_indices == 1 && 0 == target_dim && source_dim != target_dim) {
175
176 // element of mConnectivityMap should be for this type and for one
177 // less than source_dim, which should give the connectivity of that sub element
178 const ConnMap &cm = mConnectivityMap[this_type][source_dim-1];
179 std::copy(cm.conn[source_indices[0]],
180 cm.conn[source_indices[0]]+cm.num_corners_per_sub_element[source_indices[0]],
181 std::back_inserter(index_list));
182 return 0;
183 }
184
185 // now go through source indices, folding adjacencies into target list
186 for (it1 = source_indices; it1 != source_indices+num_source_indices; it1++) {
187 // *it1 is the side index
188 // at start of iteration, index_list has the target list
189
190 // if a union, or first iteration and index list was empty, copy the list
191 if (operation_type == CN::UNION ||
192 (it1 == source_indices && index_list.empty())) {
193 std::copy(MUC.targets_per_source_element[*it1],
194 MUC.targets_per_source_element[*it1]+
195 MUC.num_targets_per_source_element[*it1],
196 std::back_inserter(index_list));
197 }
198 else {
199 // else we're intersecting, and have a non-empty list; intersect with this target list
200 tmp_indices.clear();
201 for (int i = MUC.num_targets_per_source_element[*it1]-1; i>= 0; i--)
202 if (std::find(index_list.begin(), index_list.end(), MUC.targets_per_source_element[*it1][i]) !=
203 index_list.end())
204 tmp_indices.push_back(MUC.targets_per_source_element[*it1][i]);
205 // std::set_intersection(MUC.targets_per_source_element[*it1],
206 // MUC.targets_per_source_element[*it1]+
207 // MUC.num_targets_per_source_element[*it1],
208 // index_list.begin(), index_list.end(),
209 // std::back_inserter(tmp_indices));
210 index_list.swap(tmp_indices);
211
212 // if we're at this point and the list is empty, the intersection will be NULL;
213 // return if so
214 if (index_list.empty()) return 0;
215 }
216 }
217
218 if (operation_type == CN::UNION && num_source_indices != 1) {
219 // need to sort then unique the list
220 std::sort(index_list.begin(), index_list.end());
221 index_list.erase(std::unique(index_list.begin(), index_list.end()),
222 index_list.end());
223 }
224
225 return 0;
226 }
227
228 template <typename T> static
side_number(const T * parent_conn,const EntityType parent_type,const T * child_conn,const int child_num_verts,const int child_dim,int & side_no,int & sense,int & offset)229 short int side_number(const T *parent_conn,
230 const EntityType parent_type,
231 const T *child_conn,
232 const int child_num_verts,
233 const int child_dim,
234 int &side_no,
235 int &sense,
236 int &offset)
237 {
238 int parent_num_verts = CN::VerticesPerEntity(parent_type);
239 int side_indices[8];
240 assert(sizeof(side_indices)/sizeof(side_indices[0]) >= (size_t)child_num_verts);
241
242 for (int i = 0; i < child_num_verts; i++) {
243 side_indices[i] = std::find(parent_conn, parent_conn+parent_num_verts, child_conn[i]) - parent_conn;
244 if (side_indices[i] == parent_num_verts)
245 return -1;
246 }
247
248 return CN::SideNumber(parent_type, &side_indices[0], child_num_verts,
249 child_dim, side_no, sense, offset);
250 }
251
SideNumber(const EntityType parent_type,const int * parent_conn,const int * child_conn,const int child_num_verts,const int child_dim,int & side_no,int & sense,int & offset)252 short int CN::SideNumber(const EntityType parent_type, const int *parent_conn,
253 const int *child_conn, const int child_num_verts,
254 const int child_dim,
255 int &side_no, int &sense, int &offset)
256 {
257 return side_number(parent_conn, parent_type, child_conn, child_num_verts,
258 child_dim, side_no, sense, offset);
259 }
260
SideNumber(const EntityType parent_type,const unsigned int * parent_conn,const unsigned int * child_conn,const int child_num_verts,const int child_dim,int & side_no,int & sense,int & offset)261 short int CN::SideNumber(const EntityType parent_type, const unsigned int *parent_conn,
262 const unsigned int *child_conn, const int child_num_verts,
263 const int child_dim,
264 int &side_no, int &sense, int &offset)
265 {
266 return side_number(parent_conn, parent_type, child_conn, child_num_verts,
267 child_dim, side_no, sense, offset);
268 }
SideNumber(const EntityType parent_type,const long * parent_conn,const long * child_conn,const int child_num_verts,const int child_dim,int & side_no,int & sense,int & offset)269 short int CN::SideNumber(const EntityType parent_type, const long *parent_conn,
270 const long *child_conn, const int child_num_verts,
271 const int child_dim,
272 int &side_no, int &sense, int &offset)
273 {
274 return side_number(parent_conn, parent_type, child_conn, child_num_verts,
275 child_dim, side_no, sense, offset);
276 }
SideNumber(const EntityType parent_type,const unsigned long * parent_conn,const unsigned long * child_conn,const int child_num_verts,const int child_dim,int & side_no,int & sense,int & offset)277 short int CN::SideNumber(const EntityType parent_type, const unsigned long *parent_conn,
278 const unsigned long *child_conn, const int child_num_verts,
279 const int child_dim,
280 int &side_no, int &sense, int &offset)
281 {
282 return side_number(parent_conn, parent_type, child_conn, child_num_verts,
283 child_dim, side_no, sense, offset);
284 }
285
SideNumber(const EntityType parent_type,const unsigned long long * parent_conn,const unsigned long long * child_conn,const int child_num_verts,const int child_dim,int & side_no,int & sense,int & offset)286 short int CN::SideNumber(const EntityType parent_type, const unsigned long long *parent_conn,
287 const unsigned long long *child_conn, const int child_num_verts,
288 const int child_dim,
289 int &side_no, int &sense, int &offset)
290
291 {
292 return side_number(parent_conn, parent_type, child_conn, child_num_verts,
293 child_dim, side_no, sense, offset);
294 }
295
SideNumber(const EntityType parent_type,void * const * parent_conn,void * const * child_conn,const int child_num_verts,const int child_dim,int & side_no,int & sense,int & offset)296 short int CN::SideNumber(const EntityType parent_type, void * const *parent_conn,
297 void * const *child_conn, const int child_num_verts,
298 const int child_dim,
299 int &side_no, int &sense, int &offset)
300 {
301 return side_number(parent_conn, parent_type, child_conn, child_num_verts,
302 child_dim, side_no, sense, offset);
303 }
304
SideNumber(const EntityType parent_type,const int * child_conn_indices,const int child_num_verts,const int child_dim,int & side_no,int & sense,int & offset)305 short int CN::SideNumber( const EntityType parent_type,
306 const int *child_conn_indices,
307 const int child_num_verts,
308 const int child_dim,
309 int &side_no,
310 int &sense,
311 int &offset )
312 {
313 int parent_dim = Dimension(parent_type);
314 int parent_num_verts = VerticesPerEntity(parent_type);
315
316 // degenerate case (vertex), output == input
317 if (child_dim == 0) {
318 if (child_num_verts != 1)
319 return -1;
320 side_no = *child_conn_indices;
321 sense = offset = 0;
322 }
323
324 // given a parent and child element, find the corresponding side number
325
326 // dim_diff should be -1, 0 or 1 (same dimension, one less dimension, two less, resp.)
327 if (child_dim > parent_dim || child_dim < 0)
328 return -1;
329
330 // different types of same dimension won't be the same
331 if (parent_dim == child_dim &&
332 parent_num_verts != child_num_verts) {
333 side_no = -1;
334 sense = 0;
335 return 0;
336 }
337
338 // loop over the sub-elements, comparing to child connectivity
339 int sub_conn_indices[10];
340 for (int i = 0; i < NumSubEntities(parent_type, child_dim); i++) {
341 int sub_size = VerticesPerEntity(SubEntityType(parent_type, child_dim, i));
342 if (sub_size != child_num_verts)
343 continue;
344
345 SubEntityVertexIndices(parent_type, child_dim, i, sub_conn_indices);
346 bool they_match = ConnectivityMatch(child_conn_indices,
347 sub_conn_indices, sub_size,
348 sense, offset);
349 if (they_match) {
350 side_no = i;
351 return 0;
352 }
353 }
354
355 // if we've gotten here, we don't match
356 side_no = -1;
357
358 // return value is no success
359 return 1;
360 }
361
362 //! return the dimension and index of the opposite side, given parent entity type and child
363 //! dimension and index. This function is only defined for certain types of parent/child types:
364 //! (Parent, Child dim->Opposite dim):
365 //! (Tri, 1->0), (Tri, 0->1), (Quad, 1->1), (Quad, 0->0),
366 //! (Tet, 2->0), (Tet, 1->1), (Tet, 0->2),
367 //! (Hex, 2->2), (Hex, 1->1)(diagonally across element), (Hex, 0->0) (diagonally across element)
368 //! All other parent types and child dimensions return an error.
369 //!
370 //! \param parent_type The type of parent element
371 //! \param child_type The type of child element
372 //! \param child_index The index of the child element
373 //! \param opposite_index The index of the opposite element
374 //! \return status Returns 0 if successful, -1 if not
OppositeSide(const EntityType parent_type,const int child_index,const int child_dim,int & opposite_index,int & opposite_dim)375 short int CN::OppositeSide(const EntityType parent_type,
376 const int child_index,
377 const int child_dim,
378 int &opposite_index,
379 int &opposite_dim)
380 {
381 switch (parent_type) {
382 case MBEDGE:
383 if (0 != child_dim) return -1;
384 else opposite_index = 1-child_index;
385 opposite_dim = 0;
386 break;
387
388 case MBTRI:
389 switch (child_dim) {
390 case 0:
391 opposite_dim = 1;
392 opposite_index = (child_index+1)%3;
393 break;
394 case 1:
395 opposite_dim = 0;
396 opposite_index = (child_index+2)%3;
397 break;
398 default:
399 return -1;
400 }
401 break;
402
403 case MBQUAD:
404 switch (child_dim) {
405 case 0:
406 case 1:
407 opposite_dim = child_dim;
408 opposite_index = (child_index+2)%4;
409 break;
410 default:
411 return -1;
412 }
413 break;
414
415 case MBTET:
416 switch (child_dim) {
417 case 0:
418 opposite_dim = 2;
419 opposite_index = (child_index+1)%3 + 2*(child_index/3);
420 break;
421 case 1:
422 opposite_dim = 1;
423 opposite_index = child_index < 3
424 ? 3 + (child_index + 2)%3
425 : (child_index + 1)%3;
426 break;
427 case 2:
428 opposite_dim = 0;
429 opposite_index = (child_index+2)%3 + child_index/3;
430 break;
431 default:
432 return -1;
433 }
434 break;
435 case MBHEX:
436 opposite_dim = child_dim;
437 switch (child_dim) {
438 case 0:
439 opposite_index = child_index < 4
440 ? 4 + (child_index + 2) % 4
441 : (child_index - 2) % 4;
442 break;
443 case 1:
444 opposite_index = 4*(2-child_index/4) + (child_index+2)%4;
445 break;
446 case 2:
447 opposite_index = child_index < 4
448 ? (child_index + 2) % 4
449 : 9 - child_index;
450 break;
451 default:
452 return -1;
453 }
454 break;
455
456
457 default:
458 return -1;
459 }
460
461 return 0;
462 }
463
464 template <typename T>
connectivity_match(const T * conn1_i,const T * conn2_i,const int num_vertices,int & direct,int & offset)465 inline bool connectivity_match( const T* conn1_i,
466 const T* conn2_i,
467 const int num_vertices,
468 int& direct, int& offset )
469 {
470
471 bool they_match;
472
473 // special test for 2 handles, since we don't want to wrap the list in this
474 // case
475 if (num_vertices == 2) {
476 they_match = false;
477 if (conn1_i[0] == conn2_i[0] && conn1_i[1] == conn2_i[1]) {
478 direct = 1;
479 they_match = true;
480 offset = 0;
481 }
482 else if (conn1_i[0] == conn2_i[1] && conn1_i[1] == conn2_i[0]) {
483 they_match = true;
484 direct = -1;
485 offset = 1;
486 }
487 }
488
489 else {
490 const T *iter;
491 iter = std::find(&conn2_i[0], &conn2_i[num_vertices], conn1_i[0]);
492 if(iter == &conn2_i[num_vertices])
493 return false;
494
495 they_match = true;
496
497 offset = iter - conn2_i;
498 int i;
499
500 // first compare forward
501 for(i = 1; i<num_vertices; ++i)
502 {
503 if(conn1_i[i] != conn2_i[(offset+i)%num_vertices])
504 {
505 they_match = false;
506 break;
507 }
508 }
509
510 if(they_match == true)
511 {
512 direct = 1;
513 return they_match;
514 }
515
516 they_match = true;
517
518 // then compare reverse
519 for(i = 1; i<num_vertices; i++)
520 {
521 if(conn1_i[i] != conn2_i[(offset+num_vertices-i)%num_vertices])
522 {
523 they_match = false;
524 break;
525 }
526 }
527 if (they_match)
528 {
529 direct = -1;
530 }
531 }
532
533 return they_match;
534 }
535
536
ConnectivityMatch(const int * conn1_i,const int * conn2_i,const int num_vertices,int & direct,int & offset)537 bool CN::ConnectivityMatch( const int *conn1_i,
538 const int *conn2_i,
539 const int num_vertices,
540 int &direct, int &offset )
541 {
542 return connectivity_match<int>(conn1_i, conn2_i, num_vertices, direct, offset );
543 }
544
ConnectivityMatch(const unsigned int * conn1_i,const unsigned int * conn2_i,const int num_vertices,int & direct,int & offset)545 bool CN::ConnectivityMatch( const unsigned int *conn1_i,
546 const unsigned int *conn2_i,
547 const int num_vertices,
548 int &direct, int &offset )
549 {
550 return connectivity_match<unsigned int>(conn1_i, conn2_i, num_vertices, direct, offset );
551 }
552
ConnectivityMatch(const long * conn1_i,const long * conn2_i,const int num_vertices,int & direct,int & offset)553 bool CN::ConnectivityMatch( const long *conn1_i,
554 const long *conn2_i,
555 const int num_vertices,
556 int &direct, int &offset )
557 {
558 return connectivity_match<long>(conn1_i, conn2_i, num_vertices, direct, offset );
559 }
560
ConnectivityMatch(const unsigned long * conn1_i,const unsigned long * conn2_i,const int num_vertices,int & direct,int & offset)561 bool CN::ConnectivityMatch( const unsigned long *conn1_i,
562 const unsigned long *conn2_i,
563 const int num_vertices,
564 int &direct, int &offset )
565 {
566 return connectivity_match<unsigned long>(conn1_i, conn2_i, num_vertices, direct, offset );
567 }
568
ConnectivityMatch(const unsigned long long * conn1_i,const unsigned long long * conn2_i,const int num_vertices,int & direct,int & offset)569 bool CN::ConnectivityMatch( const unsigned long long *conn1_i,
570 const unsigned long long *conn2_i,
571 const int num_vertices,
572 int &direct, int &offset )
573 {
574 return connectivity_match<unsigned long long>(conn1_i, conn2_i, num_vertices, direct, offset );
575 }
576
ConnectivityMatch(void * const * conn1_i,void * const * conn2_i,const int num_vertices,int & direct,int & offset)577 bool CN::ConnectivityMatch( void* const *conn1_i,
578 void* const *conn2_i,
579 const int num_vertices,
580 int &direct, int &offset )
581 {
582 return connectivity_match<void*>(conn1_i, conn2_i, num_vertices, direct, offset );
583 }
584
585
586
587 //! for an entity of this type and a specified subfacet (dimension and index), return
588 //! the index of the higher order node for that entity in this entity's connectivity array
HONodeIndex(const EntityType this_type,const int num_verts,const int subfacet_dim,const int subfacet_index)589 short int CN::HONodeIndex(const EntityType this_type, const int num_verts,
590 const int subfacet_dim, const int subfacet_index)
591 {
592 int i;
593 int has_mids[4];
594 HasMidNodes(this_type, num_verts, has_mids);
595
596 // if we have no mid nodes on the subfacet_dim, we have no index
597 if (subfacet_index != -1 && !has_mids[subfacet_dim]) return -1;
598
599 // put start index at last index (one less than the number of vertices
600 // plus the index basis)
601 int index = VerticesPerEntity(this_type) - 1 + numberBasis;
602
603 // for each subfacet dimension less than the target subfacet dim which has mid nodes,
604 // add the number of subfacets of that dimension to the index
605 for (i = 1; i < subfacet_dim; i++)
606 if (has_mids[i]) index += NumSubEntities(this_type, i);
607
608
609 // now add the index of this subfacet, or one if we're asking about the entity as a whole
610 if (subfacet_index == -1 && has_mids[subfacet_dim])
611 // want the index of the last ho node on this subfacet
612 index += NumSubEntities(this_type, subfacet_dim);
613
614 else if (subfacet_index != -1 && has_mids[subfacet_dim])
615 index += subfacet_index + 1 - numberBasis;
616
617 // that's it
618 return index;
619 }
620
621 //! given data about an element and a vertex in that element, return the dimension
622 //! and index of the sub-entity that the vertex resolves. If it does not resolve a
623 //! sub-entity, either because it's a corner node or it's not in the element, -1 is
624 //! returned in both return values
HONodeParent(EntityType elem_type,int num_verts,int ho_index,int & parent_dim,int & parent_index)625 void CN::HONodeParent( EntityType elem_type,
626 int num_verts,
627 int ho_index,
628 int& parent_dim,
629 int& parent_index )
630 {
631 // begin with error values
632 parent_dim = parent_index = -1;
633
634 // given the number of verts and the element type, get the hasmidnodes solution
635 int has_mids[4];
636 HasMidNodes(elem_type, num_verts, has_mids);
637
638 int index = VerticesPerEntity(elem_type)-1;
639 const int dim = Dimension(elem_type);
640
641 // keep a running sum of the ho node indices for this type of element, and stop
642 // when you get to the dimension which has the ho node
643 for (int i = 1; i < dim; i++) {
644 if (has_mids[i]) {
645 if (ho_index <= index + NumSubEntities(elem_type, i)) {
646 // the ho_index resolves an entity of dimension i, so set the return values
647 // and break out of the loop
648 parent_dim = i;
649 parent_index = ho_index - index - 1;
650 return;
651 }
652 else {
653 index += NumSubEntities(elem_type, i);
654 }
655 }
656 }
657
658 // mid region node case
659 if( has_mids[dim] && ho_index == index+1 ) {
660 parent_dim = dim;
661 parent_index = 0;
662 }
663 }
664
665 } // namespace moab
666
667
668 using moab::CN;
669 using moab::EntityType;
670
671 //! get the basis of the numbering system
MBCN_GetBasis(int * rval)672 void MBCN_GetBasis(int *rval) {*rval = CN::GetBasis();}
673
674 //! set the basis of the numbering system
MBCN_SetBasis(const int in_basis)675 void MBCN_SetBasis(const int in_basis) {CN::SetBasis(in_basis);}
676
677 //! return the string type name for this type
MBCN_EntityTypeName(const int this_type,char * rval,int rval_len)678 void MBCN_EntityTypeName(const int this_type, char *rval, int rval_len)
679 {
680 const char *rval_tmp = CN::EntityTypeName((EntityType)this_type);
681 int rval_len_tmp = strlen(rval_tmp);
682 rval_len_tmp = (rval_len_tmp < rval_len ? rval_len_tmp : rval_len);
683 strncpy(rval, rval_tmp, rval_len_tmp);
684 }
685
686 //! given a name, find the corresponding entity type
MBCN_EntityTypeFromName(const char * name,int * rval)687 void MBCN_EntityTypeFromName(const char *name, int *rval)
688 {
689 *rval = CN::EntityTypeFromName(name);
690 }
691
692 //! return the topological entity dimension
MBCN_Dimension(const int t,int * rval)693 void MBCN_Dimension(const int t, int *rval)
694 {
695 *rval = CN::Dimension((EntityType)t);
696 }
697
698 //! return the number of (corner) vertices contained in the specified type.
MBCN_VerticesPerEntity(const int t,int * rval)699 void MBCN_VerticesPerEntity(const int t, int *rval)
700 {
701 *rval = CN::VerticesPerEntity((EntityType)t);
702 }
703
704 //! return the number of sub-entities bounding the entity.
MBCN_NumSubEntities(const int t,const int d,int * rval)705 void MBCN_NumSubEntities(const int t, const int d, int *rval)
706 {
707 *rval = CN::NumSubEntities((EntityType)t, d);
708 }
709
710 //! return the type of a particular sub-entity.
711 //! \param this_type Type of entity for which sub-entity type is being queried
712 //! \param sub_dimension Topological dimension of sub-entity whose type is being queried
713 //! \param index Index of sub-entity whose type is being queried
714 //! \return type Entity type of sub-entity with specified dimension and index
MBCN_SubEntityType(const int this_type,const int sub_dimension,const int index,int * rval)715 void MBCN_SubEntityType(const int this_type,
716 const int sub_dimension,
717 const int index, int *rval)
718
719 {
720
721 *rval = CN::SubEntityType((EntityType)this_type, sub_dimension, index);
722
723 }
724
725
726 //! return the vertex indices of the specified sub-entity.
727 //! \param this_type Type of entity for which sub-entity connectivity is being queried
728 //! \param sub_dimension Dimension of sub-entity
729 //! \param sub_index Index of sub-entity
730 //! \param sub_entity_conn Connectivity of sub-entity (returned to calling function)
MBCN_SubEntityVertexIndices(const int this_type,const int sub_dimension,const int sub_index,int sub_entity_conn[])731 void MBCN_SubEntityVertexIndices(const int this_type,
732 const int sub_dimension,
733 const int sub_index,
734 int sub_entity_conn[])
735 {
736 CN::SubEntityVertexIndices((EntityType)this_type, sub_dimension,
737 sub_index, sub_entity_conn);
738 }
739
740 //! return the vertices of the specified sub entity
741 //! \param parent_conn Connectivity of parent entity
742 //! \param parent_type Entity type of parent entity
743 //! \param sub_dimension Dimension of sub-entity being queried
744 //! \param sub_index Index of sub-entity being queried
745 //! \param sub_entity_conn Connectivity of sub-entity, based on parent_conn and canonical
746 //! ordering for parent_type
747 //! \param num_sub_vertices Number of vertices in sub-entity
748 // void MBCN_SubEntityConn(const void *parent_conn, const int parent_type,
749 // const int sub_dimension,
750 // const int sub_index,
751 // void *sub_entity_conn, int &num_sub_vertices) {return CN::SubEntityConn();}
752
753 //! For a specified set of sides of given dimension, return the intersection
754 //! or union of all sides of specified target dimension adjacent to those sides.
755 //! \param this_type Type of entity for which sub-entity connectivity is being queried
756 //! \param source_indices Indices of sides being queried
757 //! \param num_source_indices Number of entries in <em>source_indices</em>
758 //! \param source_dim Dimension of source entity
759 //! \param target_dim Dimension of target entity
760 //! \param index_list Indices of target entities (returned)
761 //! \param num_indices Number of indices of target entities (returned)
762 //! \param operation_type Specify either CN::INTERSECT (0) or CN::UNION (1) to get intersection
763 //! or union of target entity lists over source entities
764 //! \param rval Error code indicating success or failure (returned)
MBCN_AdjacentSubEntities(const int this_type,const int * source_indices,const int num_source_indices,const int source_dim,const int target_dim,int * index_list,int * num_indices,const int operation_type,int * rval)765 void MBCN_AdjacentSubEntities(const int this_type,
766 const int *source_indices,
767 const int num_source_indices,
768 const int source_dim,
769 const int target_dim,
770 int *index_list,
771 int *num_indices,
772 const int operation_type, int *rval)
773 {
774 std::vector<int> tmp_index_list;
775 *rval = CN::AdjacentSubEntities((EntityType)this_type, source_indices,
776 num_source_indices, source_dim, target_dim,
777 tmp_index_list, operation_type);
778 std::copy(tmp_index_list.begin(), tmp_index_list.end(), index_list);
779 *num_indices = tmp_index_list.size();
780 }
781
782 //! return the side index represented in the input sub-entity connectivity
783 //! \param parent_type Entity type of parent entity
784 //! \param child_conn_indices Child connectivity to query, specified as indices
785 //! into the connectivity list of the parent.
786 //! \param child_num_verts Number of values in <em>child_conn_indices</em>
787 //! \param child_dim Dimension of child entity being queried
788 //! \param side_no Side number of child entity (returned)
789 //! \param sense Sense of child entity with respect to order in <em>child_conn</em> (returned)
790 //! \param offset Offset of <em>child_conn</em> with respect to canonical ordering data (returned)
791 //! \return status Returns zero if successful, -1 if not
MBCN_SideNumber(const int parent_type,const int * child_conn_indices,const int child_num_verts,const int child_dim,int * side_no,int * sense,int * offset)792 void MBCN_SideNumber(const int parent_type,
793 const int *child_conn_indices, const int child_num_verts,
794 const int child_dim,
795 int *side_no, int *sense, int *offset)
796 {
797 CN::SideNumber((EntityType)parent_type, child_conn_indices, child_num_verts, child_dim,
798 *side_no, *sense, *offset);
799 }
800
MBCN_SideNumberInt(const int * parent_conn,const EntityType parent_type,const int * child_conn,const int child_num_verts,const int child_dim,int * side_no,int * sense,int * offset)801 void MBCN_SideNumberInt(const int *parent_conn, const EntityType parent_type,
802 const int *child_conn, const int child_num_verts,
803 const int child_dim,
804 int *side_no, int *sense, int *offset)
805 {
806 moab::side_number(parent_conn, parent_type, child_conn, child_num_verts,
807 child_dim, *side_no, *sense, *offset);
808 }
809
MBCN_SideNumberUint(const unsigned int * parent_conn,const EntityType parent_type,const unsigned int * child_conn,const int child_num_verts,const int child_dim,int * side_no,int * sense,int * offset)810 void MBCN_SideNumberUint(const unsigned int *parent_conn, const EntityType parent_type,
811 const unsigned int *child_conn, const int child_num_verts,
812 const int child_dim,
813 int *side_no, int *sense, int *offset)
814 {
815 moab::side_number(parent_conn, parent_type, child_conn, child_num_verts,
816 child_dim, *side_no, *sense, *offset);
817 }
818
MBCN_SideNumberLong(const long * parent_conn,const EntityType parent_type,const long * child_conn,const int child_num_verts,const int child_dim,int * side_no,int * sense,int * offset)819 void MBCN_SideNumberLong(const long *parent_conn, const EntityType parent_type,
820 const long *child_conn, const int child_num_verts,
821 const int child_dim,
822 int *side_no, int *sense, int *offset)
823 {
824 moab::side_number(parent_conn, parent_type, child_conn, child_num_verts,
825 child_dim, *side_no, *sense, *offset);
826 }
827
MBCN_SideNumberUlong(const unsigned long * parent_conn,const EntityType parent_type,const unsigned long * child_conn,const int child_num_verts,const int child_dim,int * side_no,int * sense,int * offset)828 void MBCN_SideNumberUlong(const unsigned long *parent_conn, const EntityType parent_type,
829 const unsigned long *child_conn, const int child_num_verts,
830 const int child_dim,
831 int *side_no, int *sense, int *offset)
832 {
833 moab::side_number(parent_conn, parent_type, child_conn, child_num_verts,
834 child_dim, *side_no, *sense, *offset);
835 }
836
MBCN_SideNumberVoid(void * const * parent_conn,const EntityType parent_type,void * const * child_conn,const int child_num_verts,const int child_dim,int * side_no,int * sense,int * offset)837 void MBCN_SideNumberVoid(void * const *parent_conn, const EntityType parent_type,
838 void * const *child_conn, const int child_num_verts,
839 const int child_dim,
840 int *side_no, int *sense, int *offset)
841 {
842 moab::side_number(parent_conn, parent_type, child_conn, child_num_verts,
843 child_dim, *side_no, *sense, *offset);
844 }
845
846 //! return the dimension and index of the opposite side, given parent entity type and child
847 //! dimension and index. This function is only defined for certain types of parent/child types:
848 //! (Parent, Child dim->Opposite dim):
849 //! (Tri, 1->0), (Tri, 0->1), (Quad, 1->1), (Quad, 0->0),
850 //! (Tet, 2->0), (Tet, 1->1), (Tet, 0->2),
851 //! (Hex, 2->2), (Hex, 1->1)(diagonally across element), (Hex, 0->0) (diagonally across element)
852 //! All other parent types and child dimensions return an error.
853 //!
854 //! \param parent_type The type of parent element
855 //! \param child_type The type of child element
856 //! \param child_index The index of the child element
857 //! \param opposite_index The index of the opposite element
858 //! \return status Returns 0 if successful, -1 if not
MBCN_OppositeSide(const int parent_type,const int child_index,const int child_dim,int * opposite_index,int * opposite_dim,int * rval)859 void MBCN_OppositeSide(const int parent_type,
860 const int child_index,
861 const int child_dim,
862 int *opposite_index,
863 int *opposite_dim, int *rval)
864 {
865 *rval = CN::OppositeSide((EntityType)parent_type, child_index, child_dim,
866 *opposite_index, *opposite_dim);
867 }
868
869 //! given two connectivity arrays, determine whether or not they represent the same entity.
870 //! \param conn1 Connectivity array of first entity
871 //! \param conn2 Connectivity array of second entity
872 //! \param num_vertices Number of entries in <em>conn1</em> and <em>conn2</em>
873 //! \param direct If positive, entities have the same sense (returned)
874 //! \param offset Offset of <em>conn2</em>'s first vertex in <em>conn1</em>
875 //! \return rval Returns true if <em>conn1</em> and <em>conn2</em> match
MBCN_ConnectivityMatchInt(const int * conn1,const int * conn2,const int num_vertices,int * direct,int * offset,int * rval)876 void MBCN_ConnectivityMatchInt(const int *conn1,
877 const int *conn2,
878 const int num_vertices,
879 int *direct, int *offset, int *rval)
880 {
881 *rval = CN::ConnectivityMatch(conn1, conn2, num_vertices,
882 *direct, *offset);
883 }
884
MBCN_ConnectivityMatchUint(const unsigned int * conn1,const unsigned int * conn2,const int num_vertices,int * direct,int * offset,int * rval)885 void MBCN_ConnectivityMatchUint(const unsigned int *conn1,
886 const unsigned int *conn2,
887 const int num_vertices,
888 int *direct, int *offset, int *rval)
889 {
890 *rval = CN::ConnectivityMatch(conn1, conn2, num_vertices,
891 *direct, *offset);
892 }
893
MBCN_ConnectivityMatchLong(const long * conn1,const long * conn2,const int num_vertices,int * direct,int * offset,int * rval)894 void MBCN_ConnectivityMatchLong(const long* conn1,
895 const long* conn2,
896 const int num_vertices,
897 int* direct, int* offset , int *rval)
898 {
899 *rval = CN::ConnectivityMatch(conn1, conn2, num_vertices,
900 *direct, *offset);
901 }
902
MBCN_ConnectivityMatchUlong(const unsigned long * conn1,const unsigned long * conn2,const int num_vertices,int * direct,int * offset,int * rval)903 void MBCN_ConnectivityMatchUlong(const unsigned long* conn1,
904 const unsigned long* conn2,
905 const int num_vertices,
906 int *direct, int* offset , int *rval)
907 {
908 *rval = CN::ConnectivityMatch(conn1, conn2, num_vertices,
909 *direct, *offset);
910 }
911
MBCN_ConnectivityMatchVoid(void * const * conn1,void * const * conn2,const int num_vertices,int * direct,int * offset,int * rval)912 void MBCN_ConnectivityMatchVoid(void* const* conn1,
913 void* const* conn2,
914 const int num_vertices,
915 int* direct, int* offset , int *rval)
916 {
917 *rval = CN::ConnectivityMatch(conn1, conn2, num_vertices,
918 *direct, *offset);
919 }
920
921 //! true if entities of a given type and number of nodes indicates mid edge nodes are present.
922 //! \param this_type Type of entity for which sub-entity connectivity is being queried
923 //! \param num_verts Number of nodes defining entity
924 //! \return int Returns true if <em>this_type</em> combined with <em>num_nodes</em> indicates
925 //! mid-edge nodes are likely
MBCN_HasMidEdgeNodes(const int this_type,const int num_verts,int * rval)926 void MBCN_HasMidEdgeNodes(const int this_type,
927 const int num_verts, int *rval)
928 {
929 *rval = CN::HasMidEdgeNodes((EntityType)this_type, num_verts);
930 }
931
932 //! true if entities of a given type and number of nodes indicates mid face nodes are present.
933 //! \param this_type Type of entity for which sub-entity connectivity is being queried
934 //! \param num_verts Number of nodes defining entity
935 //! \return int Returns true if <em>this_type</em> combined with <em>num_nodes</em> indicates
936 //! mid-face nodes are likely
MBCN_HasMidFaceNodes(const int this_type,const int num_verts,int * rval)937 void MBCN_HasMidFaceNodes(const int this_type,
938 const int num_verts, int *rval)
939 {
940 *rval = CN::HasMidFaceNodes((EntityType)this_type, num_verts);
941 }
942
943 //! true if entities of a given type and number of nodes indicates mid region nodes are present.
944 //! \param this_type Type of entity for which sub-entity connectivity is being queried
945 //! \param num_verts Number of nodes defining entity
946 //! \return int Returns true if <em>this_type</em> combined with <em>num_nodes</em> indicates
947 //! mid-region nodes are likely
MBCN_HasMidRegionNodes(const int this_type,const int num_verts,int * rval)948 void MBCN_HasMidRegionNodes(const int this_type,
949 const int num_verts, int *rval)
950 {
951 *rval = CN::HasMidRegionNodes((EntityType)this_type, num_verts);
952 }
953
954 //! true if entities of a given type and number of nodes indicates mid edge/face/region nodes
955 //! are present.
956 //! \param this_type Type of entity for which sub-entity connectivity is being queried
957 //! \param num_verts Number of nodes defining entity
958 //! \param mid_nodes If <em>mid_nodes[i], i=1..3</em> is true, indicates that mid-edge
959 //! (i=1), mid-face (i=2), and/or mid-region (i=3) nodes are likely
MBCN_HasMidNodes(const int this_type,const int num_verts,int mid_nodes[4])960 void MBCN_HasMidNodes(const int this_type,
961 const int num_verts,
962 int mid_nodes[4])
963 {
964 return CN::HasMidNodes((EntityType)this_type, num_verts, mid_nodes);
965 }
966
967 //! given data about an element and a vertex in that element, return the dimension
968 //! and index of the sub-entity that the vertex resolves. If it does not resolve a
969 //! sub-entity, either because it's a corner node or it's not in the element, -1 is
970 //! returned in both return values.
971 //! \param elem_type Type of entity being queried
972 //! \param num_nodes The number of nodes in the element connectivity
973 //! \param ho_node_index The position of the HO node in the connectivity list (zero based)
974 //! \param parent_dim Dimension of sub-entity high-order node resolves (returned)
975 //! \param parent_index Index of sub-entity high-order node resolves (returned)
MBCN_HONodeParent(int elem_type,int num_nodes,int ho_node_index,int * parent_dim,int * parent_index)976 void MBCN_HONodeParent( int elem_type,
977 int num_nodes,
978 int ho_node_index,
979 int *parent_dim,
980 int *parent_index )
981 {
982 return CN::HONodeParent((EntityType)elem_type, num_nodes, ho_node_index,
983 *parent_dim, *parent_index);
984 }
985
986 //! for an entity of this type with num_verts vertices, and a specified subfacet
987 //! (dimension and index), return the index of the higher order node for that entity
988 //! in this entity's connectivity array
989 //! \param this_type Type of entity being queried
990 //! \param num_verts Number of vertices for the entity being queried
991 //! \param subfacet_dim Dimension of sub-entity being queried
992 //! \param subfacet_index Index of sub-entity being queried
993 //! \return index Index of sub-entity's higher-order node
MBCN_HONodeIndex(const int this_type,const int num_verts,const int subfacet_dim,const int subfacet_index,int * rval)994 void MBCN_HONodeIndex(const int this_type, const int num_verts,
995 const int subfacet_dim, const int subfacet_index, int *rval)
996
997 {
998
999 *rval = CN::HONodeIndex((EntityType)this_type, num_verts, subfacet_dim, subfacet_index);
1000
1001 }
1002
1003 namespace moab {
1004
1005 template <typename T>
permute_this(EntityType t,const int dim,T * conn,const int indices_per_ent,const int num_entries)1006 inline int permute_this(EntityType t,
1007 const int dim,
1008 T* conn,
1009 const int indices_per_ent,
1010 const int num_entries)
1011 {
1012 T tmp_conn[MAX_SUB_ENTITIES];
1013 assert(indices_per_ent <= CN::permuteVec[t][dim][MAX_SUB_ENTITIES]);
1014 if (indices_per_ent > CN::permuteVec[t][dim][MAX_SUB_ENTITIES]) return 1;
1015 short int *tvec = CN::permuteVec[t][dim];
1016 T *pvec = conn;
1017 for (int j = 0; j < num_entries; j++) {
1018 for (int i = 0; i < indices_per_ent; i++)
1019 tmp_conn[tvec[i]] = pvec[i];
1020 memcpy(pvec, tmp_conn, indices_per_ent*sizeof(T));
1021 pvec += indices_per_ent;
1022 }
1023
1024 return 0;
1025 }
1026
1027 template <typename T>
rev_permute_this(EntityType t,const int dim,T * conn,const int indices_per_ent,const int num_entries)1028 inline int rev_permute_this(EntityType t,
1029 const int dim,
1030 T* conn,
1031 const int indices_per_ent,
1032 const int num_entries)
1033 {
1034 T tmp_conn[MAX_SUB_ENTITIES];
1035 assert(indices_per_ent <= CN::revPermuteVec[t][dim][MAX_SUB_ENTITIES]);
1036 if (indices_per_ent > CN::revPermuteVec[t][dim][MAX_SUB_ENTITIES]) return 1;
1037 short int *tvec = CN::revPermuteVec[t][dim];
1038 T *pvec = conn;
1039 for (int j = 0; j < num_entries; j++) {
1040 for (int i = 0; i < indices_per_ent; i++)
1041 tmp_conn[i] = pvec[tvec[i]];
1042 memcpy(pvec, tmp_conn, indices_per_ent*sizeof(T));
1043 pvec += indices_per_ent;
1044 }
1045
1046 return 0;
1047 }
1048
1049 //! Permute this vector
permuteThis(const EntityType t,const int dim,int * pvec,const int num_indices,const int num_entries)1050 inline int CN::permuteThis(const EntityType t, const int dim, int *pvec,
1051 const int num_indices, const int num_entries)
1052 {return permute_this(t, dim, pvec, num_indices, num_entries);}
permuteThis(const EntityType t,const int dim,unsigned int * pvec,const int num_indices,const int num_entries)1053 inline int CN::permuteThis(const EntityType t, const int dim, unsigned int *pvec,
1054 const int num_indices, const int num_entries)
1055 {return permute_this(t, dim, pvec, num_indices, num_entries);}
permuteThis(const EntityType t,const int dim,long * pvec,const int num_indices,const int num_entries)1056 inline int CN::permuteThis(const EntityType t, const int dim, long *pvec,
1057 const int num_indices, const int num_entries)
1058 {return permute_this(t, dim, pvec, num_indices, num_entries);}
permuteThis(const EntityType t,const int dim,void ** pvec,const int num_indices,const int num_entries)1059 inline int CN::permuteThis(const EntityType t, const int dim, void **pvec,
1060 const int num_indices, const int num_entries)
1061 {return permute_this(t, dim, pvec, num_indices, num_entries);}
1062
1063 //! Reverse permute this vector
revPermuteThis(const EntityType t,const int dim,int * pvec,const int num_indices,const int num_entries)1064 inline int CN::revPermuteThis(const EntityType t, const int dim, int *pvec,
1065 const int num_indices, const int num_entries)
1066 {return rev_permute_this(t, dim, pvec, num_indices, num_entries);}
revPermuteThis(const EntityType t,const int dim,unsigned int * pvec,const int num_indices,const int num_entries)1067 inline int CN::revPermuteThis(const EntityType t, const int dim, unsigned int *pvec,
1068 const int num_indices, const int num_entries)
1069 {return rev_permute_this(t, dim, pvec, num_indices, num_entries);}
revPermuteThis(const EntityType t,const int dim,long * pvec,const int num_indices,const int num_entries)1070 inline int CN::revPermuteThis(const EntityType t, const int dim, long *pvec,
1071 const int num_indices, const int num_entries)
1072 {return rev_permute_this(t, dim, pvec, num_indices, num_entries);}
revPermuteThis(const EntityType t,const int dim,void ** pvec,const int num_indices,const int num_entries)1073 inline int CN::revPermuteThis(const EntityType t, const int dim, void **pvec,
1074 const int num_indices, const int num_entries)
1075 {return rev_permute_this(t, dim, pvec, num_indices, num_entries);}
1076
1077
1078 } // namespace moab
1079