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