1 /* *****************************************************************
2 MESQUITE -- The Mesh Quality Improvement Toolkit
3
4 Copyright 2004 Sandia Corporation and Argonne National
5 Laboratory. Under the terms of Contract DE-AC04-94AL85000
6 with Sandia Corporation, the U.S. Government 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 diachin2@llnl.gov, djmelan@sandia.gov, mbrewer@sandia.gov,
24 pknupp@sandia.gov, tleurent@mcs.anl.gov, tmunson@mcs.anl.gov
25
26 ***************************************************************** */
27 //
28 // ORIG-DATE: 16-May-02 at 10:26:21
29 // LAST-MOD: 18-Jun-04 at 11:36:07 by Thomas Leurent
30 //
31 /*! \file MsqMeshEntity.cpp
32
33 \brief All elements in Mesquite are of type MsqMeshEntity. Their associated
34 functionality is implemented in this file.
35
36 \author Thomas Leurent
37 \author Michael Brewer
38 \author Darryl Melander
39 \date 2002-05-16
40 */
41
42 #include "Mesquite.hpp"
43 #include "MsqMeshEntity.hpp"
44 #include "MsqVertex.hpp"
45 #include "PatchData.hpp"
46
47 #include <vector>
48 #include <ostream>
49 using std::vector;
50 using std::ostream;
51 using std::endl;
52
53 namespace MBMesquite {
54
55 //! Gets the indices of the vertices of this element.
56 //! The indices are only valid in the PatchData from which
57 //! this element was retrieved.
58 //! The order of the vertices is the canonical order for this
59 //! element's type.
get_vertex_indices(std::vector<std::size_t> & vertices) const60 void MsqMeshEntity::get_vertex_indices(std::vector<std::size_t> &vertices) const
61 {
62 vertices.resize( vertex_count() );
63 std::copy( vertexIndices, vertexIndices + vertex_count(), vertices.begin() );
64 }
65
66 //! Gets the indices of the vertices of this element.
67 //! The indices are only valid in the PatchData from which
68 //! this element was retrieved.
69 //! The order of the vertices is the canonical order for this
70 //! element's type.
71 //! The indices are placed appended to the end of the list.
72 //! The list is not cleared before appending this entity's vertices.
append_vertex_indices(std::vector<std::size_t> & vertex_list) const73 void MsqMeshEntity::append_vertex_indices(std::vector<std::size_t> &vertex_list) const
74 {
75 vertex_list.insert(vertex_list.end(),
76 vertexIndices,
77 vertexIndices + vertex_count());
78 }
79
get_node_indices(std::vector<std::size_t> & indices) const80 void MsqMeshEntity::get_node_indices( std::vector<std::size_t>& indices ) const
81 {
82 indices.resize( node_count() );
83 std::copy( vertexIndices, vertexIndices + node_count(), indices.begin() );
84 }
85
append_node_indices(std::vector<std::size_t> & indices) const86 void MsqMeshEntity::append_node_indices( std::vector<std::size_t>& indices ) const
87 {
88 indices.insert( indices.end(), vertexIndices, vertexIndices + node_count() );
89 }
90
91
92 /*! The centroid of an element containing n vertices with equal masses is located at
93 \f[ \b{x} = \frac{ \sum_{i=1}^{n} \b{x}_i }{ n } \f]
94 where \f$ \b{x}_i ,\, i=1,...,n\f$ are the vertices coordinates.
95 */
get_centroid(Vector3D & centroid,const PatchData & pd,MsqError & err) const96 void MsqMeshEntity::get_centroid(Vector3D ¢roid, const PatchData &pd, MsqError &err) const
97 {
98 centroid = 0.0;
99 const MsqVertex* vtces = pd.get_vertex_array(err); MSQ_ERRRTN(err);
100 size_t nve = vertex_count();
101 for (size_t i=0; i<nve; ++i)
102 centroid += vtces[vertexIndices[i]];
103 centroid /= nve;
104 }
105
106
corner_volume(const Vector3D & v0,const Vector3D & v1,const Vector3D & v2,const Vector3D & v3)107 static inline double corner_volume( const Vector3D& v0,
108 const Vector3D& v1,
109 const Vector3D& v2,
110 const Vector3D& v3 )
111 {
112 return (v1 - v0) * (v2 - v0) % (v3 - v0);
113 }
114
115 /*!
116 \brief Computes the area of the given element. Returned value is
117 always non-negative. If the entity passed is not a two-dimensional
118 element, an error is set.*/
compute_unsigned_area(PatchData & pd,MsqError & err)119 double MsqMeshEntity::compute_unsigned_area(PatchData &pd, MsqError &err) {
120 const MsqVertex* verts=pd.get_vertex_array(err);MSQ_ERRZERO(err);
121 double tem=0.0;
122 switch (mType)
123 {
124
125 case TRIANGLE:
126 tem = ((verts[vertexIndices[1]]-verts[vertexIndices[0]])*
127 (verts[vertexIndices[2]]-verts[vertexIndices[0]])).length();
128 return 0.5*tem;
129
130 case QUADRILATERAL:
131 tem = ((verts[vertexIndices[1]]-verts[vertexIndices[0]])*
132 (verts[vertexIndices[3]]-verts[vertexIndices[0]])).length();
133 tem += ((verts[vertexIndices[3]]-verts[vertexIndices[2]])*
134 (verts[vertexIndices[1]]-verts[vertexIndices[2]])).length();
135 return (tem/2.0);
136
137 case POLYGON:
138 // assume convex
139 for (unsigned i = 1; i < numVertexIndices-1; ++i)
140 tem += ((verts[vertexIndices[i]] - verts[vertexIndices[0]]) *
141 (verts[vertexIndices[i+1]] - verts[vertexIndices[0]])).length();
142 return 0.5 * tem;
143
144 case TETRAHEDRON:
145 return 1.0/6.0 * fabs( corner_volume( verts[vertexIndices[0]],
146 verts[vertexIndices[1]],
147 verts[vertexIndices[2]],
148 verts[vertexIndices[3]] ) );
149
150 case PYRAMID: {
151 Vector3D m = verts[vertexIndices[0]] + verts[vertexIndices[1]]
152 + verts[vertexIndices[2]] + verts[vertexIndices[3]];
153 Vector3D t1 = verts[vertexIndices[0]] - verts[vertexIndices[2]];
154 Vector3D t2 = verts[vertexIndices[1]] - verts[vertexIndices[3]];
155 tem = ((t1 + t2) * (t1 - t2)) % (verts[vertexIndices[4]] - 0.25 * m);
156 return (1.0/12.0) * fabs(tem);
157 }
158
159 case PRISM: {
160 tem = corner_volume( verts[vertexIndices[0]],
161 verts[vertexIndices[1]],
162 verts[vertexIndices[2]],
163 verts[vertexIndices[3]] );
164
165 tem += corner_volume( verts[vertexIndices[1]],
166 verts[vertexIndices[2]],
167 verts[vertexIndices[3]],
168 verts[vertexIndices[4]] );
169
170 tem += corner_volume( verts[vertexIndices[2]],
171 verts[vertexIndices[3]],
172 verts[vertexIndices[4]],
173 verts[vertexIndices[5]] );
174
175 return 1.0/6.0 * fabs(tem);
176 }
177
178 case HEXAHEDRON: {
179
180 tem = corner_volume( verts[vertexIndices[1]],
181 verts[vertexIndices[2]],
182 verts[vertexIndices[0]],
183 verts[vertexIndices[5]] );
184
185 tem += corner_volume( verts[vertexIndices[3]],
186 verts[vertexIndices[0]],
187 verts[vertexIndices[2]],
188 verts[vertexIndices[7]] );
189
190 tem += corner_volume( verts[vertexIndices[4]],
191 verts[vertexIndices[7]],
192 verts[vertexIndices[5]],
193 verts[vertexIndices[0]] );
194
195 tem += corner_volume( verts[vertexIndices[6]],
196 verts[vertexIndices[5]],
197 verts[vertexIndices[7]],
198 verts[vertexIndices[2]] );
199
200 tem += corner_volume( verts[vertexIndices[5]],
201 verts[vertexIndices[2]],
202 verts[vertexIndices[0]],
203 verts[vertexIndices[7]] );
204
205 return (1.0/6.0) * fabs(tem);
206 }
207
208 default:
209 MSQ_SETERR(err)("Invalid type of element passed to compute unsigned area.",
210 MsqError::UNSUPPORTED_ELEMENT);
211 return 0;
212 }
213 return 0;
214 }
215
216 /*!
217 \brief Computes the area of the given element. Returned value can be
218 negative. If the entity passed is not a two-dimensional element, an
219 error is set.*/
compute_signed_area(PatchData & pd,MsqError & err)220 double MsqMeshEntity::compute_signed_area(PatchData &pd, MsqError &err) {
221 const MsqVertex* verts=pd.get_vertex_array(err);MSQ_ERRZERO(err);
222 double tem=0.0;
223 double tem2=0.0;
224 Vector3D surface_normal;
225 Vector3D cross_vec;
226 size_t element_index=pd.get_element_index(this);
227
228 switch (mType)
229 {
230
231 case TRIANGLE:
232 cross_vec=((verts[vertexIndices[1]]-verts[vertexIndices[0]])*
233 (verts[vertexIndices[2]]-verts[vertexIndices[0]]));
234 pd.get_domain_normal_at_element(element_index,surface_normal,err);
235 MSQ_ERRZERO(err);
236 tem = (cross_vec.length()/2.0);
237 //if normals do not point in same general direction, negate area
238 if(cross_vec%surface_normal<0){
239 tem *= -1;
240 }
241
242 return tem;
243
244 case QUADRILATERAL:
245 cross_vec=((verts[vertexIndices[1]]-verts[vertexIndices[0]])*
246 (verts[vertexIndices[3]]-verts[vertexIndices[0]]));
247 pd.get_domain_normal_at_element(element_index,surface_normal,err);
248 MSQ_ERRZERO(err);
249 tem = (cross_vec.length()/2.0);
250 //if normals do not point in same general direction, negate area
251 if(cross_vec%surface_normal<0){
252 tem *= -1;
253 }
254 cross_vec=((verts[vertexIndices[3]]-verts[vertexIndices[2]])*
255 (verts[vertexIndices[1]]-verts[vertexIndices[2]]));
256 tem2 = (cross_vec.length()/2.0);
257 //if normals do not point in same general direction, negate area
258 if(cross_vec%surface_normal<0){
259 tem2 *= -1;
260 //test to make sure surface normal existed
261 //if(surface_normal.length_squared()<.5){
262 //err.set_msg("compute_signed_area called without surface_normal available.");
263 //}
264 }
265 return (tem + tem2);
266
267 default:
268 MSQ_SETERR(err)("Invalid type of element passed to compute unsigned area.",
269 MsqError::UNSUPPORTED_ELEMENT);
270 return 0;
271 };
272 return 0.0;
273 }
274
275 /*!Appends the indices (in the vertex array) of the vertices to connected
276 to vertex_array[vertex_index] to the end of the vector vert_indices.
277 The connected vertices are right-hand ordered as defined by the
278 entity.
279
280 */
get_connected_vertices(std::size_t vertex_index,std::vector<std::size_t> & vert_indices,MsqError & err)281 void MsqMeshEntity::get_connected_vertices(
282 std::size_t vertex_index,
283 std::vector<std::size_t> &vert_indices,
284 MsqError &err)
285 {
286 //index is set to the index in the vertexIndices corresponding
287 //to vertex_index
288 int index;
289 for (index = vertex_count() - 1; index >= 0; --index)
290 if (vertexIndices[index] == vertex_index)
291 break;
292 if (index < 0)
293 {
294 MSQ_SETERR(err)("Invalid vertex index.", MsqError::INVALID_ARG);
295 return;
296 }
297
298 unsigned n;
299 const unsigned* indices = TopologyInfo::adjacent_vertices( mType, index, n );
300 if (!indices)
301 MSQ_SETERR(err)("Element type not available", MsqError::INVALID_ARG);
302 for (unsigned i = 0; i < n; ++i)
303 vert_indices.push_back( vertexIndices[indices[i]] );
304 }
305
306 /*! Gives the normal at the surface point corner_pt ... but if not available,
307 gives the normalized cross product of corner_vec1 and corner_vec2.
308 */
309 /*
310 void MsqMeshEntity::compute_corner_normal(size_t corner,
311 Vector3D &normal,
312 PatchData &pd,
313 MsqError &err)
314 {
315 if ( get_element_type()==TRIANGLE || get_element_type()==QUADRILATERAL )
316 {
317 // There are two cases where we cannot get a normal from the
318 // geometry that are not errors:
319 // 1) There is no domain set
320 // 2) The vertex is at a degenerate point on the geometry (e.g.
321 // tip of a cone.)
322
323 bool have_normal = false;
324
325 // Get normal from domain
326 if (pd.domain_set())
327 {
328 size_t index = pd.get_element_index(this);
329 pd.get_domain_normal_at_corner( index, corner, normal, err );
330 MSQ_ERRRTN(err);
331
332 double length = normal.length();
333 if (length > DBL_EPSILON)
334 {
335 have_normal = true;
336 normal /= length;
337 }
338 }
339
340 // If failed to get normal from domain, calculate
341 // from adjacent vertices.
342 if ( !have_normal )
343 {
344 const int num_corner = this->vertex_count();
345 const int prev_corner = (corner + num_corner - 1) % num_corner;
346 const int next_corner = (corner + 1 ) % num_corner;
347 const size_t this_idx = vertexIndices[corner];
348 const size_t prev_idx = vertexIndices[prev_corner];
349 const size_t next_idx = vertexIndices[next_corner];
350 normal = (pd.vertex_by_index(next_idx) - pd.vertex_by_index(this_idx))
351 * (pd.vertex_by_index(prev_idx) - pd.vertex_by_index(this_idx));
352 normal.normalize();
353 }
354 }
355 else
356 MSQ_SETERR(err)("Should only be used for faces (tri, quads, ...).",
357 MsqError::INVALID_ARG);
358 }
359 */
360
compute_corner_normals(Vector3D normals[],PatchData & pd,MsqError & err)361 void MsqMeshEntity::compute_corner_normals( Vector3D normals[],
362 PatchData &pd,
363 MsqError &err)
364 {
365 EntityTopology type = get_element_type();
366 if (type != TRIANGLE && type != QUADRILATERAL && type != POLYGON)
367 {
368 MSQ_SETERR(err)("Should only be used for faces (tri, quads, ...).",
369 MsqError::INVALID_ARG);
370 return;
371 }
372
373
374 // There are two cases where we cannot get a normal from the
375 // geometry that are not errors:
376 // 1) There is no domain set
377 // 2) The vertex is at a degenerate point on the geometry (e.g.
378 // tip of a cone.)
379
380 // Get normal from domain
381 if (pd.domain_set())
382 {
383 size_t index = pd.get_element_index(this);
384 pd.get_domain_normals_at_corners( index, normals, err );
385 MSQ_ERRRTN(err);
386 }
387
388 // Check if normals are valid (none are valid if !pd.domain_set())
389 const unsigned count = vertex_count();
390 for (unsigned i = 0; i < count; ++i)
391 {
392 // If got valid normal from domain,
393 // make it a unit vector and continue.
394 if (pd.domain_set())
395 {
396 double length = normals[i].length();
397 if (length > DBL_EPSILON)
398 {
399 normals[i] /= length;
400 continue;
401 }
402 }
403
404 const size_t prev_idx = vertexIndices[(i + count - 1) % count];
405 const size_t this_idx = vertexIndices[i];
406 const size_t next_idx = vertexIndices[(i + 1) % count];
407
408 // Calculate normal using edges adjacent to corner
409 normals[i] = (pd.vertex_by_index(next_idx) - pd.vertex_by_index(this_idx))
410 * (pd.vertex_by_index(prev_idx) - pd.vertex_by_index(this_idx));
411 normals[i].normalize();
412 }
413 }
414
operator <<(ostream & stream,const MsqMeshEntity & entity)415 ostream& operator<<( ostream& stream, const MsqMeshEntity& entity )
416 {
417 size_t num_vert = entity.vertex_count();
418 stream << "MsqMeshEntity " << &entity << " with vertices ";
419 for (size_t i = 0; i < num_vert; ++i)
420 stream << entity.vertexIndices[i] << " ";
421 stream << endl;
422 return stream;
423 }
424
425
426
427
428 /*! Get a array of indices that specifies for the given vertex
429 the correct matrix map. This is used by the I_DFT point
430 relaxation methods in the laplacian smoothers.
431
432 */
get_local_matrix_map_about_vertex(PatchData & pd,MsqVertex * vert,size_t local_map_size,int * local_map,MsqError & err) const433 size_t MsqMeshEntity::get_local_matrix_map_about_vertex(
434 PatchData &pd, MsqVertex* vert, size_t local_map_size,
435 int* local_map, MsqError &err) const
436 {
437 //i iterates through elem's vertices
438 int i=0;
439 //index is set to the index in the vertexIndices corresponding
440 //to vertex_index
441 int index=-1;
442 int return_val=0;
443 const MsqVertex* vertex_array = pd.get_vertex_array(err);
444 if(err)
445 return return_val;
446
447 switch (mType)
448 {
449 case TRIANGLE:
450 MSQ_SETERR(err)("Requested function not yet supported for Triangles.",
451 MsqError::NOT_IMPLEMENTED);
452
453 break;
454
455 case QUADRILATERAL:
456 MSQ_SETERR(err)("Requested function not yet supported for Quadrilaterals.",
457 MsqError::NOT_IMPLEMENTED);
458
459 break;
460
461 case TETRAHEDRON:
462 if(local_map_size<4){
463 MSQ_SETERR(err)("Array of incorrect length sent to function.",
464 MsqError::INVALID_ARG);
465 return return_val;
466 }
467 return_val = 4;
468 while(i<4)
469 {
470 if(&vertex_array[vertexIndices[i]]==vert)
471 {
472 index=i;
473 break;
474 }
475 ++i;
476 }
477 switch(index){
478 case(0):
479 local_map[0]=0;
480 local_map[1]=1;
481 local_map[2]=2;
482 local_map[3]=3;
483 break;
484 case(1):
485 local_map[0]=1;
486 local_map[1]=0;
487 local_map[2]=3;
488 local_map[3]=2;
489 break;
490 case(2):
491 local_map[0]=2;
492 local_map[1]=3;
493 local_map[2]=0;
494 local_map[3]=1;
495 break;
496 case(3):
497 local_map[0]=3;
498 local_map[1]=2;
499 local_map[2]=1;
500 local_map[3]=0;
501 break;
502 default:
503 local_map[0]=-1;
504 local_map[1]=-1;
505 local_map[2]=-1;
506 local_map[3]=-1;
507 };
508
509 break;
510
511 case HEXAHEDRON:
512 MSQ_SETERR(err)("Requested function not yet supported for Hexahedrons.",
513 MsqError::NOT_IMPLEMENTED);
514
515 break;
516 default:
517 MSQ_SETERR(err)("Element type not available", MsqError::UNSUPPORTED_ELEMENT);
518 break;
519 }
520 return return_val;
521
522 }
523
524
check_element_orientation(PatchData & pd,int & inverted,int & total,MsqError & err)525 void MsqMeshEntity::check_element_orientation(
526 PatchData &pd, int& inverted, int& total, MsqError &err)
527 {
528 NodeSet all = all_nodes( err ); MSQ_ERRRTN(err);
529 unsigned i;
530
531 if (TopologyInfo::dimension(mType) == 2) {
532 if (!pd.domain_set()) {
533 total = 0;
534 inverted = 0;
535 return;
536 }
537
538 const MappingFunction2D* mf = pd.get_mapping_function_2D( mType );
539 if (!mf) {
540 check_element_orientation_corners( pd, inverted, total, err );
541 return;
542 }
543
544 NodeSet sample = mf->sample_points( all );
545 total = sample.num_nodes();
546 inverted = 0;
547
548 if (sample.have_any_corner_node()) {
549 for (i = 0; i < TopologyInfo::corners(mType); ++i)
550 if (sample.corner_node(i))
551 inverted += inverted_jacobian_2d( pd, all, Sample(0,i), err );
552 }
553 if (sample.have_any_mid_edge_node()) {
554 for (i = 0; i < TopologyInfo::edges(mType); ++i)
555 if (sample.mid_edge_node(i))
556 inverted += inverted_jacobian_2d( pd, all, Sample(1,i), err );
557 }
558 if (sample.have_any_mid_face_node())
559 inverted += inverted_jacobian_2d( pd, all, Sample(2,0), err );
560 }
561 else {
562 const MappingFunction3D* mf = pd.get_mapping_function_3D( mType );
563 if (!mf) {
564 check_element_orientation_corners( pd, inverted, total, err );
565 return;
566 }
567
568 NodeSet sample = mf->sample_points( all );
569 total = sample.num_nodes();
570 inverted = 0;
571
572 if (sample.have_any_corner_node()) {
573 for (i = 0; i < TopologyInfo::corners(mType); ++i)
574 if (sample.corner_node(i))
575 inverted += inverted_jacobian_3d( pd, all, Sample(0,i), err );
576 }
577 if (sample.have_any_mid_edge_node()) {
578 for (i = 0; i < TopologyInfo::edges(mType); ++i)
579 if (sample.mid_edge_node(i))
580 inverted += inverted_jacobian_3d( pd, all, Sample(1,i), err );
581 }
582 if (sample.have_any_mid_face_node()) {
583 for (i = 0; i < TopologyInfo::faces(mType); ++i)
584 if (sample.mid_face_node(i))
585 inverted += inverted_jacobian_3d( pd, all, Sample(2,i), err );
586 }
587 if (sample.have_any_mid_region_node()) {
588 inverted += inverted_jacobian_3d( pd, all, Sample(3,0), err );
589 }
590 }
591 }
592
593 bool
inverted_jacobian_3d(PatchData & pd,NodeSet nodes,Sample sample,MsqError & err)594 MsqMeshEntity::inverted_jacobian_3d( PatchData& pd, NodeSet nodes, Sample sample, MsqError& err )
595 {
596 MsqMatrix<3,3> J;
597 MsqVector<3> junk[27];
598 size_t junk2[27], junk3;
599 assert(node_count() <= 27);
600
601 const MappingFunction3D* mf = pd.get_mapping_function_3D( mType );
602 mf->jacobian( pd, pd.get_element_index(this), nodes,
603 sample, junk2, junk, junk3, J, err );
604 MSQ_ERRZERO(err);
605 //const double size_eps_sqr = sqr_Frobenius( J ) * DBL_EPSILON;
606 const double d = det(J);
607 double l1 = J.column(0) % J.column(0);
608 double l2 = J.column(1) % J.column(1);
609 double l3 = J.column(2) % J.column(2);
610 return d < 0 || d*d < DBL_EPSILON*DBL_EPSILON * l1*l2*l3;
611 }
612
613 bool
inverted_jacobian_2d(PatchData & pd,NodeSet nodes,Sample sample,MsqError & err)614 MsqMeshEntity::inverted_jacobian_2d( PatchData& pd, NodeSet nodes, Sample sample, MsqError& err )
615 {
616 MsqMatrix<3,2> J;
617 MsqVector<2> junk[9];
618 size_t junk2[9], junk3;
619 assert(node_count() <= 9);
620
621 const int idx = pd.get_element_index(this);
622 const MappingFunction2D* mf = pd.get_mapping_function_2D( mType );
623 mf->jacobian( pd, idx, nodes, sample, junk2, junk, junk3, J, err );
624 MSQ_ERRZERO(err);
625 const MsqVector<3> cross = J.column(0) * J.column(1);
626
627 if (pd.domain_set()) {
628 Vector3D norm;
629 pd.get_domain_normal_at_sample( pd.get_element_index(this), sample, norm, err );
630 MSQ_ERRZERO(err);
631
632 const MsqVector<3> N(&norm[0]);
633 if (cross % N < 0.0)
634 return true;
635 }
636
637 const double l1 = J.column(0) % J.column(0);
638 const double l2 = J.column(1) % J.column(1);
639 return cross % cross < DBL_EPSILON*DBL_EPSILON * l1*l2;
640 }
641
642 NodeSet
all_nodes(MsqError & err) const643 MsqMeshEntity::all_nodes( MsqError& err ) const
644 {
645 bool mid_edge, mid_face, mid_vol;
646 TopologyInfo::higher_order( mType, node_count(), mid_edge, mid_face, mid_vol, err );
647 NodeSet result;
648 result.set_all_corner_nodes( mType );
649 if (mid_edge)
650 result.set_all_mid_edge_nodes( mType );
651 if (mid_face)
652 result.set_all_mid_face_nodes( mType );
653 if (mid_vol)
654 result.set_all_mid_region_nodes( mType );
655 return result;
656 }
657
check_element_orientation_corners(PatchData & pd,int & inverted,int & total,MsqError & err)658 void MsqMeshEntity::check_element_orientation_corners(
659 PatchData &pd,
660 int& inverted,
661 int& total,
662 MsqError &err)
663 {
664 int num_nodes = node_count();
665 total = inverted = 0;
666
667 if (node_count() > vertex_count()) {
668 MSQ_SETERR(err)("Cannot perform inversion test for higher-order element"
669 " without mapping function.", MsqError::UNSUPPORTED_ELEMENT );
670 return;
671 }
672
673 const MsqVertex *vertices = pd.get_vertex_array(err); MSQ_ERRRTN(err);
674
675 const Vector3D d_con(1.0, 1.0, 1.0);
676
677 int i;
678 Vector3D coord_vectors[3];
679 Vector3D center_vector;
680
681 switch(mType) {
682 case TRIANGLE:
683
684 if (!pd.domain_set())
685 return;
686
687 pd.get_domain_normal_at_element(this, coord_vectors[2], err); MSQ_ERRRTN(err);
688 coord_vectors[2] = coord_vectors[2] / coord_vectors[2].length();// Need unit normal
689 center_vector = vertices[vertexIndices[0]];
690 coord_vectors[0] = vertices[vertexIndices[1]]-center_vector;
691 coord_vectors[1] = vertices[vertexIndices[2]]-center_vector;
692 total = 1;
693 inverted = (coord_vectors[2]%(coord_vectors[0]*coord_vectors[1] ) <= 0.0);
694 break;
695
696 case QUADRILATERAL:
697
698 if (!pd.domain_set())
699 return;
700
701 pd.get_domain_normal_at_element(this, coord_vectors[2], err); MSQ_ERRRTN(err);
702 coord_vectors[2] = coord_vectors[2] / coord_vectors[2].length();// Need unit normal
703
704 for (i = 0; i < 4; ++i) {
705 center_vector = vertices[vertexIndices[i]];
706 coord_vectors[0] = vertices[vertexIndices[(i+1)%4]]-center_vector;
707 coord_vectors[1] = vertices[vertexIndices[(i+3)%4]]-center_vector;
708 ++total;
709 inverted += (coord_vectors[2]%(coord_vectors[0]*coord_vectors[1] ) <= 0.0);
710 }
711 break;
712
713 case TETRAHEDRON:
714 center_vector = vertices[vertexIndices[0]];
715 coord_vectors[0] = vertices[vertexIndices[1]]-center_vector;
716 coord_vectors[1] = vertices[vertexIndices[2]]-center_vector;
717 coord_vectors[2] = vertices[vertexIndices[3]]-center_vector;
718 total = 1;
719 inverted = ( coord_vectors[0]%(coord_vectors[1]*coord_vectors[2] ) <= 0.0);
720 break;
721
722 case POLYGON:
723
724 if (!pd.domain_set())
725 return;
726
727 pd.get_domain_normal_at_element(this, coord_vectors[2], err); MSQ_ERRRTN(err);
728 coord_vectors[2] = coord_vectors[2] / coord_vectors[2].length();// Need unit normal
729
730 for (i = 0; i < num_nodes; ++i) {
731 center_vector = vertices[vertexIndices[i]];
732 coord_vectors[0] = vertices[vertexIndices[(i+1)%num_nodes]]-center_vector;
733 coord_vectors[1] = vertices[vertexIndices[(i+num_nodes-1)%num_nodes]]-center_vector;
734 ++total;
735 inverted += (coord_vectors[2]%(coord_vectors[0]*coord_vectors[1] ) <= 0.0);
736 }
737 break;
738
739 default: // generic code for 3D elements
740 {
741 size_t num_corners = corner_count();
742 unsigned num_adj;
743 const unsigned* adj_idx;
744 for (unsigned j = 0; j < num_corners; ++j)
745 {
746 adj_idx = TopologyInfo::adjacent_vertices( mType, j, num_adj );
747 if (3 != num_adj)
748 {
749 MSQ_SETERR(err)("Unsupported element type.", MsqError::INVALID_ARG);
750 return;
751 }
752
753 center_vector = vertices[vertexIndices[j]];
754 coord_vectors[0] = vertices[vertexIndices[adj_idx[0]]] - center_vector;
755 coord_vectors[1] = vertices[vertexIndices[adj_idx[1]]] - center_vector;
756 coord_vectors[2] = vertices[vertexIndices[adj_idx[2]]] - center_vector;
757 ++total;
758 inverted += (coord_vectors[0]%(coord_vectors[1]*coord_vectors[2] ) <= 0.0);
759 }
760 break;
761 }
762 } // end switch over element type
763 }
764
765 } // namespace MBMesquite
766