Branch data Line data Source code
1 : : /* *****************************************************************
2 : : MESQUITE -- The Mesh Quality Improvement Toolkit
3 : :
4 : : Copyright 2010 Sandia National Laboratories. Developed at the
5 : : University of Wisconsin--Madison under SNL contract number
6 : : 624796. The U.S. Government and the University of Wisconsin
7 : : retain certain rights to 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 : : (2010) [email protected]
24 : :
25 : : ***************************************************************** */
26 : :
27 : : /** \file MeshUtil.cpp
28 : : * \brief Implement MeshUtil class
29 : : * \author Jason Kraftcheck
30 : : */
31 : :
32 : : #include "Mesquite.hpp"
33 : : #include "MeshUtil.hpp"
34 : : #include "MeshInterface.hpp"
35 : : #include "EdgeIterator.hpp"
36 : : #include "SimpleStats.hpp"
37 : : #include "MsqError.hpp"
38 : : #include "TMPQualityMetric.hpp"
39 : : #include "MappingFunction.hpp"
40 : : #include "TargetCalculator.hpp"
41 : :
42 : : #include <vector>
43 : : #include <algorithm>
44 : : #include <limits>
45 : : #include <sstream>
46 : :
47 : : namespace MBMesquite
48 : : {
49 : :
50 : 33 : MeshUtil::~MeshUtil()
51 : : {
52 [ + - ]: 33 : delete globalPatch;
53 : 33 : }
54 : :
55 : 53 : PatchData* MeshUtil::get_global_patch( MsqError& err )
56 : : {
57 [ + + ]: 53 : if( !globalPatch )
58 : : {
59 [ + - ]: 33 : globalPatch = new PatchData;
60 : 33 : globalPatch->set_mesh( myMesh );
61 [ + + ]: 33 : if( mySettings ) globalPatch->attach_settings( mySettings );
62 : 33 : globalPatch->fill_global_patch( err );
63 [ - + ][ # # ]: 33 : if( MSQ_CHKERR( err ) )
[ - + ]
64 : : {
65 [ # # ]: 0 : delete globalPatch;
66 : 0 : globalPatch = 0;
67 : : }
68 : : }
69 : 53 : return globalPatch;
70 : : }
71 : :
72 : 30 : void MeshUtil::edge_length_distribution( SimpleStats& results, MsqError& err )
73 : : {
74 [ + - ][ + - ]: 60 : PatchData* pd = get_global_patch( err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
75 [ + - ][ + - ]: 30 : EdgeIterator iter( pd, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
76 [ + - ][ + + ]: 58399 : while( !iter.is_at_end() )
77 : : {
78 [ + - ][ + - ]: 58369 : Vector3D diff = iter.start() - iter.end();
[ + - ]
79 [ + - ][ + - ]: 58369 : results.add_squared( diff % diff );
80 [ + - ][ + - ]: 58369 : iter.step( err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
81 : : }
82 : :
83 [ + - ][ - + ]: 30 : if( results.empty() )
84 : : { // no mesh
85 [ # # ][ # # ]: 30 : MSQ_SETERR( err )( "Mesh contains no elements", MsqError::INVALID_MESH );
[ + - ]
86 : 30 : }
87 : : }
88 : :
89 : 23 : void MeshUtil::lambda_distribution( SimpleStats& results, MsqError& err )
90 : : {
91 [ + - ][ + - ]: 46 : PatchData& pd = *get_global_patch( err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
92 : :
93 [ + - ]: 23 : std::vector< size_t > handles;
94 [ + - ]: 23 : TMPQualityMetric::get_patch_evaluations( pd, handles, false, err );
95 : :
96 : 23 : const size_t N = 64;
97 : : size_t count;
98 : : size_t indices[N];
99 [ + - ][ + + ]: 1495 : MsqVector< 2 > derivs2D[N];
100 [ + - ][ + + ]: 1495 : MsqVector< 3 > derivs3D[N];
101 : :
102 [ + + ]: 73799 : for( size_t i = 0; i < handles.size(); ++i )
103 : : {
104 : : double lambda;
105 [ + - ][ + - ]: 73776 : const Sample s = ElemSampleQM::sample( handles[i] );
106 [ + - ][ + - ]: 73776 : const size_t e = ElemSampleQM::elem( handles[i] );
107 [ + - ][ + - ]: 73776 : EntityTopology type = pd.element_by_index( e ).get_element_type();
108 [ + - ]: 73776 : const NodeSet bits = pd.non_slave_node_set( e );
109 : :
110 [ + - ][ - + ]: 73776 : if( TopologyInfo::dimension( type ) == 3 )
111 : : {
112 [ # # ]: 0 : const MappingFunction3D* mf = pd.get_mapping_function_3D( type );
113 [ # # ]: 0 : if( !mf )
114 : : {
115 : : MSQ_SETERR( err )
116 [ # # ][ # # ]: 0 : ( "No mapping function for element type", MsqError::UNSUPPORTED_ELEMENT );
117 : 0 : return;
118 : : }
119 : :
120 [ # # ]: 0 : MsqMatrix< 3, 3 > W;
121 [ # # ][ # # ]: 0 : mf->jacobian( pd, e, bits, s, indices, derivs3D, count, W, err );MSQ_ERRRTN( err );
[ # # ][ # # ]
[ # # ][ # # ]
122 [ # # ]: 0 : assert( N >= count );
123 [ # # ]: 0 : lambda = TargetCalculator::size( W );
124 : : }
125 : : else
126 : : {
127 [ + - ][ - + ]: 73776 : assert( TopologyInfo::dimension( type ) == 2 );
128 [ + - ]: 73776 : const MappingFunction2D* mf = pd.get_mapping_function_2D( type );
129 [ - + ]: 73776 : if( !mf )
130 : : {
131 : : MSQ_SETERR( err )
132 [ # # ][ # # ]: 0 : ( "No mapping function for element type", MsqError::UNSUPPORTED_ELEMENT );
133 : 0 : return;
134 : : }
135 : :
136 [ + - ]: 73776 : MsqMatrix< 3, 2 > W;
137 [ + - ][ + - ]: 73776 : mf->jacobian( pd, e, bits, s, indices, derivs2D, count, W, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
138 [ - + ]: 73776 : assert( N >= count );
139 [ + - ]: 73776 : lambda = TargetCalculator::size( W );
140 : : }
141 : :
142 [ + - ]: 73776 : results.add_value( lambda );
143 : : }
144 : :
145 [ + - ][ - + ]: 23 : if( results.empty() )
146 : : { // no mesh
147 [ # # ][ # # ]: 23 : MSQ_SETERR( err )( "Mesh contains no elements", MsqError::INVALID_MESH );
[ + - ]
148 : 23 : }
149 : : }
150 : :
151 : 0 : bool MeshUtil::meshes_are_different( Mesh& mesh1, Mesh& mesh2, MsqError& err, double tol, bool do_print )
152 : : {
153 : : // MsqError& err = (do_print ? MsqPrintError(cout) : MsqError());
154 : :
155 : : // mesh 1
156 [ # # ]: 0 : std::vector< Mesh::ElementHandle > elements1;
157 [ # # ]: 0 : std::vector< Mesh::VertexHandle > vertices1;
158 [ # # ]: 0 : std::vector< Mesh::VertexHandle > connectivity1;
159 [ # # ]: 0 : std::vector< size_t > offsets1;
160 : :
161 [ # # ]: 0 : mesh1.get_all_elements( elements1, err );
162 [ # # ][ # # ]: 0 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ # # ]
163 [ # # ]: 0 : mesh1.get_all_vertices( vertices1, err );
164 [ # # ][ # # ]: 0 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ # # ]
165 [ # # ][ # # ]: 0 : mesh1.elements_get_attached_vertices( arrptr( elements1 ), elements1.size(), connectivity1, offsets1, err );
166 [ # # ][ # # ]: 0 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ # # ]
167 [ # # ]: 0 : std::vector< EntityTopology > types1( elements1.size() );
168 [ # # ][ # # ]: 0 : mesh1.elements_get_topologies( arrptr( elements1 ), arrptr( types1 ), elements1.size(), err );
[ # # ]
169 [ # # ][ # # ]: 0 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ # # ]
170 : :
171 : : // mesh 2
172 [ # # ]: 0 : std::vector< Mesh::ElementHandle > elements2;
173 [ # # ]: 0 : std::vector< Mesh::VertexHandle > vertices2;
174 [ # # ]: 0 : std::vector< Mesh::VertexHandle > connectivity2;
175 [ # # ]: 0 : std::vector< size_t > offsets2;
176 : :
177 [ # # ]: 0 : mesh2.get_all_elements( elements2, err );
178 [ # # ][ # # ]: 0 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ # # ]
179 [ # # ]: 0 : mesh2.get_all_vertices( vertices2, err );
180 [ # # ][ # # ]: 0 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ # # ]
181 [ # # ][ # # ]: 0 : mesh2.elements_get_attached_vertices( arrptr( elements2 ), elements2.size(), connectivity2, offsets2, err );
182 [ # # ][ # # ]: 0 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ # # ]
183 [ # # ]: 0 : std::vector< EntityTopology > types2( elements2.size() );
184 [ # # ][ # # ]: 0 : mesh2.elements_get_topologies( arrptr( elements2 ), arrptr( types2 ), elements2.size(), err );
[ # # ]
185 [ # # ][ # # ]: 0 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ # # ]
186 : :
187 : : #define MDRET( msg ) \
188 : : do \
189 : : { \
190 : : if( do_print ) std::cout << "MeshUtil::mesh_diff meshes are different: " << msg << std::endl; \
191 : : return true; \
192 : : } while( 0 )
193 [ # # ][ # # ]: 0 : if( elements1.size() != elements2.size() ) MDRET( "elements sizes differ" );
[ # # ][ # # ]
[ # # ]
194 [ # # ][ # # ]: 0 : if( vertices1.size() != vertices2.size() ) MDRET( "vertices sizes differ" );
[ # # ][ # # ]
[ # # ]
195 [ # # ][ # # ]: 0 : if( offsets1.size() != offsets2.size() ) MDRET( "offsets sizes differ" );
[ # # ][ # # ]
[ # # ]
196 [ # # ][ # # ]: 0 : if( connectivity1.size() != connectivity2.size() ) MDRET( "connectivity sizes differ" );
[ # # ][ # # ]
[ # # ]
197 : :
198 [ # # ]: 0 : for( unsigned i = 0; i < offsets1.size(); i++ )
199 : : {
200 [ # # ][ # # ]: 0 : if( offsets1[i] != offsets2[i] ) MDRET( "offets differ" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
201 : : }
202 : :
203 [ # # ]: 0 : for( unsigned i = 0; i < vertices1.size(); i++ )
204 : : {
205 [ # # ][ # # ]: 0 : MsqVertex vert1, vert2;
206 [ # # ][ # # ]: 0 : mesh1.vertices_get_coordinates( &vertices1[i], &vert1, 1, err );
207 [ # # ][ # # ]: 0 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ # # ]
208 [ # # ][ # # ]: 0 : mesh2.vertices_get_coordinates( &vertices2[i], &vert2, 1, err );
209 [ # # ][ # # ]: 0 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ # # ]
210 [ # # ]: 0 : double dist = Vector3D::distance_between( vert1, vert2 );
211 [ # # ]: 0 : double v1 = vert1.length();
212 [ # # ]: 0 : double v2 = vert2.length();
213 [ # # ]: 0 : if( dist > 0.5 * ( v1 + v2 ) * tol )
214 : : {
215 [ # # ]: 0 : std::ostringstream ost;
216 [ # # ][ # # ]: 0 : ost << "vertices coordinates differ more than tolerance [" << tol << "], vert1= " << vert1
[ # # ][ # # ]
217 [ # # ][ # # ]: 0 : << " vert2= " << vert2;
218 [ # # ][ # # ]: 0 : MDRET( ost.str() );
[ # # ][ # # ]
[ # # ]
219 : : }
220 : : }
221 : :
222 [ # # ]: 0 : for( unsigned i = 0; i < connectivity1.size(); i++ )
223 : : {
224 [ # # ][ # # ]: 0 : MsqVertex vert1, vert2;
225 [ # # ][ # # ]: 0 : mesh1.vertices_get_coordinates( &connectivity1[i], &vert1, 1, err );
226 [ # # ][ # # ]: 0 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ # # ]
227 [ # # ][ # # ]: 0 : mesh2.vertices_get_coordinates( &connectivity2[i], &vert2, 1, err );
228 [ # # ][ # # ]: 0 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ # # ]
229 [ # # ]: 0 : double dist = Vector3D::distance_between( vert1, vert2 );
230 [ # # ]: 0 : double v1 = vert1.length();
231 [ # # ]: 0 : double v2 = vert2.length();
232 [ # # ]: 0 : if( dist > 0.5 * ( v1 + v2 ) * tol )
233 : : {
234 [ # # ]: 0 : std::ostringstream ost;
235 [ # # ][ # # ]: 0 : ost << "connectivity coordinates differ more than tolerance [" << tol << "], vert1= " << vert1
[ # # ][ # # ]
236 [ # # ][ # # ]: 0 : << " vert2= " << vert2;
237 [ # # ][ # # ]: 0 : MDRET( ost.str() );
[ # # ][ # # ]
[ # # ]
238 : : }
239 : : }
240 : :
241 : 0 : return false;
242 : : }
243 : :
244 [ + - ][ + - ]: 120 : } // namespace MBMesquite
|
