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 : : \file ConditionNumberQualityMetric.cpp
29 : : \brief
30 : :
31 : : \author Michael Brewer
32 : : \date 2002-06-9
33 : : */
34 : : #include <vector>
35 : : #include "ConditionNumberQualityMetric.hpp"
36 : : #include <math.h>
37 : : #include "Vector3D.hpp"
38 : : #include "ConditionNumberFunctions.hpp"
39 : :
40 : : using namespace MBMesquite;
41 : :
42 [ + - ]: 9 : ConditionNumberQualityMetric::ConditionNumberQualityMetric() : AveragingQM( QualityMetric::LINEAR ) {}
43 : :
44 : 20 : std::string ConditionNumberQualityMetric::get_name() const
45 : : {
46 [ + - ]: 20 : return "Condition Number";
47 : : }
48 : :
49 : 361 : int ConditionNumberQualityMetric::get_negate_flag() const
50 : : {
51 : 361 : return 1;
52 : : }
53 : :
54 : 1999700 : bool ConditionNumberQualityMetric::evaluate( PatchData& pd, size_t handle, double& fval, MsqError& err )
55 : : {
56 [ + - ]: 1999700 : const MsqMeshEntity* const element = &pd.element_by_index( handle );
57 : : bool return_flag;
58 : : double met_vals[MSQ_MAX_NUM_VERT_PER_ENT];
59 : 1999700 : fval = MSQ_MAX_CAP;
60 [ + - ]: 1999700 : const size_t* v_i = element->get_vertex_index_array();
61 : : // only 3 temp_vec will be sent to cond-num calculator, but the
62 : : // additional vector3Ds may be needed during the calculations
63 [ + - ][ + + ]: 13997900 : Vector3D temp_vec[6];
64 [ + - ]: 1999700 : const MsqVertex* vertices = pd.get_vertex_array( err );
65 [ + - ]: 1999700 : EntityTopology type = element->get_element_type();
66 [ - + + + : 1999700 : switch( type )
+ + - ]
67 : : {
68 : : case TRIANGLE:
69 [ # # ][ # # ]: 0 : temp_vec[0] = vertices[v_i[1]] - vertices[v_i[0]];
70 [ # # ][ # # ]: 0 : temp_vec[2] = vertices[v_i[2]] - vertices[v_i[0]];
71 : : // make relative to equilateral
72 [ # # ][ # # ]: 0 : temp_vec[1] = ( ( 2 * temp_vec[2] ) - temp_vec[0] ) * MSQ_SQRT_THREE_INV;
[ # # ][ # # ]
73 [ # # ]: 0 : return_flag = condition_number_2d( temp_vec, handle, pd, fval, err );
74 [ # # ][ # # ]: 0 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ # # ]
75 : 0 : return return_flag;
76 : : case QUADRILATERAL:
77 [ + - ][ + - ]: 544 : temp_vec[0] = vertices[v_i[1]] - vertices[v_i[0]];
78 [ + - ][ + - ]: 544 : temp_vec[1] = vertices[v_i[3]] - vertices[v_i[0]];
79 [ + - ]: 544 : return_flag = condition_number_2d( temp_vec, handle, pd, met_vals[0], err );
80 [ + - ][ - + ]: 544 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
81 [ + + ]: 544 : if( !return_flag ) return return_flag;
82 [ + - ][ + - ]: 540 : temp_vec[0] = vertices[v_i[2]] - vertices[v_i[1]];
83 [ + - ][ + - ]: 540 : temp_vec[1] = vertices[v_i[0]] - vertices[v_i[1]];
84 [ + - ]: 540 : return_flag = condition_number_2d( temp_vec, handle, pd, met_vals[1], err );
85 [ + - ][ - + ]: 540 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
86 [ - + ]: 540 : if( !return_flag ) return return_flag;
87 [ + - ][ + - ]: 540 : temp_vec[0] = vertices[v_i[3]] - vertices[v_i[2]];
88 [ + - ][ + - ]: 540 : temp_vec[1] = vertices[v_i[1]] - vertices[v_i[2]];
89 [ + - ]: 540 : return_flag = condition_number_2d( temp_vec, handle, pd, met_vals[2], err );
90 [ + - ][ - + ]: 540 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
91 [ + + ]: 540 : if( !return_flag ) return return_flag;
92 [ + - ][ + - ]: 537 : temp_vec[0] = vertices[v_i[0]] - vertices[v_i[3]];
93 [ + - ][ + - ]: 537 : temp_vec[1] = vertices[v_i[2]] - vertices[v_i[3]];
94 [ + - ]: 537 : return_flag = condition_number_2d( temp_vec, handle, pd, met_vals[3], err );
95 [ + - ][ - + ]: 537 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
96 [ + - ]: 537 : fval = average_metrics( met_vals, 4, err );
97 : 537 : return return_flag;
98 : : case TETRAHEDRON:
99 [ + - ][ + - ]: 1991876 : temp_vec[0] = vertices[v_i[1]] - vertices[v_i[0]];
100 [ + - ][ + - ]: 1991876 : temp_vec[3] = vertices[v_i[2]] - vertices[v_i[0]];
101 [ + - ][ + - ]: 1991876 : temp_vec[4] = vertices[v_i[3]] - vertices[v_i[0]];
102 : : // transform to equilateral tet
103 [ + - ][ + - ]: 1991876 : temp_vec[1] = ( ( 2 * temp_vec[3] ) - temp_vec[0] ) / MSQ_SQRT_THREE;
[ + - ][ + - ]
104 [ + - ][ + - ]: 1991876 : temp_vec[2] = ( ( 3 * temp_vec[4] ) - temp_vec[0] - temp_vec[3] ) / ( MSQ_SQRT_THREE * MSQ_SQRT_TWO );
[ + - ][ + - ]
[ + - ]
105 [ + - ]: 1991876 : return_flag = condition_number_3d( temp_vec, pd, fval, err );
106 [ + - ][ - + ]: 1991876 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
107 : 1991876 : return return_flag;
108 : :
109 : : case HEXAHEDRON:
110 : : // transform to v_i[0]
111 [ + - ][ + - ]: 32 : temp_vec[0] = vertices[v_i[1]] - vertices[v_i[0]];
112 [ + - ][ + - ]: 32 : temp_vec[1] = vertices[v_i[3]] - vertices[v_i[0]];
113 [ + - ][ + - ]: 32 : temp_vec[2] = vertices[v_i[4]] - vertices[v_i[0]];
114 [ + - ]: 32 : return_flag = condition_number_3d( temp_vec, pd, met_vals[0], err );
115 [ + - ][ - + ]: 32 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
116 [ - + ]: 32 : if( !return_flag ) return return_flag;
117 [ + - ][ + - ]: 32 : temp_vec[0] = vertices[v_i[2]] - vertices[v_i[1]];
118 [ + - ][ + - ]: 32 : temp_vec[1] = vertices[v_i[0]] - vertices[v_i[1]];
119 [ + - ][ + - ]: 32 : temp_vec[2] = vertices[v_i[5]] - vertices[v_i[1]];
120 [ + - ]: 32 : return_flag = condition_number_3d( temp_vec, pd, met_vals[1], err );
121 [ + - ][ - + ]: 32 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
122 [ - + ]: 32 : if( !return_flag ) return return_flag;
123 [ + - ][ + - ]: 32 : temp_vec[0] = vertices[v_i[3]] - vertices[v_i[2]];
124 [ + - ][ + - ]: 32 : temp_vec[1] = vertices[v_i[1]] - vertices[v_i[2]];
125 [ + - ][ + - ]: 32 : temp_vec[2] = vertices[v_i[6]] - vertices[v_i[2]];
126 [ + - ]: 32 : return_flag = condition_number_3d( temp_vec, pd, met_vals[2], err );
127 [ + - ][ - + ]: 32 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
128 [ - + ]: 32 : if( !return_flag ) return return_flag;
129 [ + - ][ + - ]: 32 : temp_vec[0] = vertices[v_i[0]] - vertices[v_i[3]];
130 [ + - ][ + - ]: 32 : temp_vec[1] = vertices[v_i[2]] - vertices[v_i[3]];
131 [ + - ][ + - ]: 32 : temp_vec[2] = vertices[v_i[7]] - vertices[v_i[3]];
132 [ + - ]: 32 : return_flag = condition_number_3d( temp_vec, pd, met_vals[3], err );
133 [ + - ][ - + ]: 32 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
134 [ - + ]: 32 : if( !return_flag ) return return_flag;
135 [ + - ][ + - ]: 32 : temp_vec[0] = vertices[v_i[7]] - vertices[v_i[4]];
136 [ + - ][ + - ]: 32 : temp_vec[1] = vertices[v_i[5]] - vertices[v_i[4]];
137 [ + - ][ + - ]: 32 : temp_vec[2] = vertices[v_i[0]] - vertices[v_i[4]];
138 [ + - ]: 32 : return_flag = condition_number_3d( temp_vec, pd, met_vals[4], err );
139 [ + - ][ - + ]: 32 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
140 [ - + ]: 32 : if( !return_flag ) return return_flag;
141 [ + - ][ + - ]: 32 : temp_vec[0] = vertices[v_i[4]] - vertices[v_i[5]];
142 [ + - ][ + - ]: 32 : temp_vec[1] = vertices[v_i[6]] - vertices[v_i[5]];
143 [ + - ][ + - ]: 32 : temp_vec[2] = vertices[v_i[1]] - vertices[v_i[5]];
144 [ + - ]: 32 : return_flag = condition_number_3d( temp_vec, pd, met_vals[5], err );
145 [ + - ][ - + ]: 32 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
146 [ - + ]: 32 : if( !return_flag ) return return_flag;
147 [ + - ][ + - ]: 32 : temp_vec[0] = vertices[v_i[5]] - vertices[v_i[6]];
148 [ + - ][ + - ]: 32 : temp_vec[1] = vertices[v_i[7]] - vertices[v_i[6]];
149 [ + - ][ + - ]: 32 : temp_vec[2] = vertices[v_i[2]] - vertices[v_i[6]];
150 [ + - ]: 32 : return_flag = condition_number_3d( temp_vec, pd, met_vals[6], err );
151 [ + - ][ - + ]: 32 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
152 [ - + ]: 32 : if( !return_flag ) return return_flag;
153 [ + - ][ + - ]: 32 : temp_vec[0] = vertices[v_i[6]] - vertices[v_i[7]];
154 [ + - ][ + - ]: 32 : temp_vec[1] = vertices[v_i[4]] - vertices[v_i[7]];
155 [ + - ][ + - ]: 32 : temp_vec[2] = vertices[v_i[3]] - vertices[v_i[7]];
156 [ + - ]: 32 : return_flag = condition_number_3d( temp_vec, pd, met_vals[7], err );
157 [ + - ][ - + ]: 32 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
158 [ + - ]: 32 : fval = average_metrics( met_vals, 8, err );
159 [ + - ][ - + ]: 32 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
160 : 32 : return return_flag;
161 : :
162 : : case PYRAMID: {
163 : : unsigned num_adj;
164 : : const unsigned* adj_idx;
165 : 2880 : return_flag = true;
166 [ + + ]: 14400 : for( size_t j = 0; j < 4; ++j ) // skip fifth vertex (apex)
167 : : {
168 [ + - ]: 11520 : adj_idx = TopologyInfo::adjacent_vertices( type, j, num_adj );
169 [ - + ]: 11520 : assert( num_adj == 3 );
170 : :
171 [ + - ][ + - ]: 11520 : temp_vec[0] = vertices[v_i[adj_idx[0]]] - vertices[v_i[j]];
172 [ + - ][ + - ]: 11520 : temp_vec[1] = vertices[v_i[adj_idx[1]]] - vertices[v_i[j]];
173 : : // calculate last vect map to right tetrahedron
174 [ + - ][ + - ]: 11520 : temp_vec[3] = vertices[v_i[adj_idx[2]]] - vertices[v_i[adj_idx[0]]];
175 [ + - ][ + - ]: 11520 : temp_vec[4] = vertices[v_i[adj_idx[2]]] - vertices[v_i[adj_idx[1]]];
176 [ + - ][ + - ]: 11520 : temp_vec[2] = 0.5 * ( temp_vec[3] + temp_vec[4] );
[ + - ]
177 : :
178 [ + - ][ + - ]: 11520 : return_flag = return_flag && condition_number_3d( temp_vec, pd, met_vals[j], err );
[ + - ]
179 : : }
180 [ + - ]: 2880 : fval = average_metrics( met_vals, 4, err );
181 [ + - ][ - + ]: 2880 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
182 : 2880 : return return_flag;
183 : : }
184 : :
185 : : case PRISM: {
186 : : unsigned num_adj;
187 : : const unsigned* adj_idx;
188 : 4368 : return_flag = true;
189 [ + + ]: 30576 : for( size_t j = 0; j < 6; ++j )
190 : : {
191 [ + - ]: 26208 : adj_idx = TopologyInfo::adjacent_vertices( type, j, num_adj );
192 [ - + ]: 26208 : assert( num_adj == 3 );
193 : :
194 [ + - ][ + - ]: 26208 : temp_vec[0] = vertices[v_i[adj_idx[0]]] - vertices[v_i[j]];
195 [ + - ][ + - ]: 26208 : temp_vec[1] = vertices[v_i[adj_idx[1]]] - vertices[v_i[j]];
196 [ + - ][ + - ]: 26208 : temp_vec[2] = vertices[v_i[adj_idx[2]]] - vertices[v_i[j]];
197 : : // map to right tetrahedron
198 [ + - ]: 26208 : temp_vec[1] += vertices[v_i[adj_idx[1]]];
199 [ + - ]: 26208 : temp_vec[1] -= vertices[v_i[adj_idx[0]]];
200 [ + - ]: 26208 : temp_vec[1] *= MSQ_SQRT_THREE_INV;
201 : :
202 [ + - ][ + - ]: 26208 : return_flag = return_flag && condition_number_3d( temp_vec, pd, met_vals[j], err );
[ + - ]
203 : : }
204 [ + - ]: 4368 : fval = average_metrics( met_vals, 6, err );
205 [ + - ][ - + ]: 4368 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
206 : 4368 : return return_flag;
207 : : }
208 : :
209 : : default:
210 : : MSQ_SETERR( err )
211 [ # # ][ # # ]: 0 : ( MsqError::UNSUPPORTED_ELEMENT, "Unsupported cell type (%d) for Condition Number quality metric.", type );
212 : :
213 : 0 : fval = MSQ_MAX_CAP;
214 : 1999700 : return false;
215 : : } // end switch over element type
216 : : return false;
217 [ + - ][ + - ]: 36 : }
|
