MOAB: Mesh Oriented datABase
(version 5.2.1)
|
00001 /* ***************************************************************** 00002 MESQUITE -- The Mesh Quality Improvement Toolkit 00003 00004 Copyright 2006 Sandia National Laboratories. Developed at the 00005 University of Wisconsin--Madison under SNL contract number 00006 624796. The U.S. Government and the University of Wisconsin 00007 retain certain rights to this software. 00008 00009 This library is free software; you can redistribute it and/or 00010 modify it under the terms of the GNU Lesser General Public 00011 License as published by the Free Software Foundation; either 00012 version 2.1 of the License, or (at your option) any later version. 00013 00014 This library is distributed in the hope that it will be useful, 00015 but WITHOUT ANY WARRANTY; without even the implied warranty of 00016 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00017 Lesser General Public License for more details. 00018 00019 You should have received a copy of the GNU Lesser General Public License 00020 (lgpl.txt) along with this library; if not, write to the Free Software 00021 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 00022 00023 (2006) kraftche@cae.wisc.edu 00024 00025 ***************************************************************** */ 00026 00027 /** \file TMPQualityMetric.cpp 00028 * \brief 00029 * \author Jason Kraftcheck 00030 */ 00031 00032 #undef PRINT_INFO 00033 00034 #include "Mesquite.hpp" 00035 #include "TMPQualityMetric.hpp" 00036 #include "MsqMatrix.hpp" 00037 #include "ElementQM.hpp" 00038 #include "MsqError.hpp" 00039 #include "Vector3D.hpp" 00040 #include "PatchData.hpp" 00041 #include "MappingFunction.hpp" 00042 #include "WeightCalculator.hpp" 00043 #include "TargetCalculator.hpp" 00044 #include "TargetMetricUtil.hpp" 00045 00046 #ifdef PRINT_INFO 00047 #include <iostream> 00048 #endif 00049 00050 #include <functional> 00051 #include <algorithm> 00052 00053 namespace MBMesquite 00054 { 00055 00056 int TMPQualityMetric::get_negate_flag() const 00057 { 00058 return 1; 00059 } 00060 00061 void TMPQualityMetric::get_evaluations( PatchData& pd, std::vector< size_t >& handles, bool free, MsqError& err ) 00062 { 00063 get_sample_pt_evaluations( pd, handles, free, err ); 00064 } 00065 00066 void TMPQualityMetric::get_patch_evaluations( PatchData& pd, std::vector< size_t >& handles, bool free, MsqError& err ) 00067 { 00068 get_sample_pt_evaluations( pd, handles, free, err ); 00069 } 00070 00071 void TMPQualityMetric::get_element_evaluations( PatchData& pd, size_t p_elem, std::vector< size_t >& handles, 00072 MsqError& err ) 00073 { 00074 get_elem_sample_points( pd, p_elem, handles, err ); 00075 } 00076 00077 bool TMPQualityMetric::evaluate( PatchData& pd, size_t p_handle, double& value, MsqError& err ) 00078 { 00079 size_t num_idx; 00080 bool valid = evaluate_internal( pd, p_handle, value, mIndices, num_idx, err ); 00081 if( MSQ_CHKERR( err ) || !valid ) return false; 00082 00083 // apply target weight to value 00084 if( weightCalc ) 00085 { 00086 const Sample s = ElemSampleQM::sample( p_handle ); 00087 const size_t e = ElemSampleQM::elem( p_handle ); 00088 double ck = weightCalc->get_weight( pd, e, s, err ); 00089 MSQ_ERRZERO( err ); 00090 value *= ck; 00091 } 00092 return true; 00093 } 00094 00095 bool TMPQualityMetric::evaluate_with_indices( PatchData& pd, size_t p_handle, double& value, 00096 std::vector< size_t >& indices, MsqError& err ) 00097 { 00098 indices.resize( MAX_ELEM_NODES ); 00099 size_t num_idx = 0; 00100 bool result = evaluate_internal( pd, p_handle, value, arrptr( indices ), num_idx, err ); 00101 if( MSQ_CHKERR( err ) || !result ) return false; 00102 00103 indices.resize( num_idx ); 00104 00105 // apply target weight to value 00106 if( weightCalc ) 00107 { 00108 const Sample s = ElemSampleQM::sample( p_handle ); 00109 const size_t e = ElemSampleQM::elem( p_handle ); 00110 double ck = weightCalc->get_weight( pd, e, s, err ); 00111 MSQ_ERRZERO( err ); 00112 value *= ck; 00113 } 00114 00115 return true; 00116 } 00117 00118 static void get_u_perp( const MsqVector< 3 >& u, MsqVector< 3 >& u_perp ) 00119 { 00120 double a = sqrt( u[0] * u[0] + u[1] * u[1] ); 00121 if( a < 1e-10 ) 00122 { 00123 u_perp[0] = 1.0; 00124 u_perp[1] = u_perp[2] = 0.0; 00125 } 00126 else 00127 { 00128 double b = -u[2] / a; 00129 u_perp[0] = u[0] * b; 00130 u_perp[1] = u[1] * b; 00131 u_perp[2] = a; 00132 } 00133 } 00134 00135 /* Do transform M_hat = S_a M_{3x2}, M_{2x2} Theta^-1 M_hat 00136 * where the plane into which we are projecting is orthogonal 00137 * to the passed u vector. 00138 */ 00139 static inline bool project_to_perp_plane( MsqMatrix< 3, 2 > J, const MsqVector< 3 >& u, const MsqVector< 3 >& u_perp, 00140 MsqMatrix< 2, 2 >& A, MsqMatrix< 3, 2 >& S_a_transpose_Theta ) 00141 { 00142 MsqVector< 3 > n_a = J.column( 0 ) * J.column( 1 ); 00143 double sc, len = length( n_a ); 00144 if( !divide( 1.0, len, sc ) ) return false; 00145 n_a *= sc; 00146 double ndot = n_a % u; 00147 double sigma = ( ndot < 0.0 ) ? -1 : 1; 00148 double cosphi = sigma * ndot; 00149 MsqVector< 3 > cross = n_a * u; 00150 double sinphi = length( cross ); 00151 00152 MsqMatrix< 3, 2 > Theta; 00153 Theta.set_column( 0, u_perp ); 00154 Theta.set_column( 1, u * u_perp ); 00155 00156 // If columns of J are not in plane orthogonal to u, then 00157 // rotate J such that they are. 00158 if( sinphi > 1e-12 ) 00159 { 00160 MsqVector< 3 > m = sigma * cross; 00161 MsqVector< 3 > n = ( 1 / sinphi ) * m; 00162 MsqVector< 3 > p = ( 1 - cosphi ) * n; 00163 double s_a[] = { p[0] * n[0] + cosphi, p[0] * n[1] - m[2], p[0] * n[2] + m[1], 00164 p[1] * n[0] + m[2], p[1] * n[1] + cosphi, p[1] * n[2] - m[0], 00165 p[2] * n[0] - m[1], p[2] * n[1] + m[0], p[2] * n[2] + cosphi }; 00166 MsqMatrix< 3, 3 > S_a( s_a ); 00167 J = S_a * J; 00168 S_a_transpose_Theta = transpose( S_a ) * Theta; 00169 } 00170 else 00171 { 00172 S_a_transpose_Theta = Theta; 00173 // J *= sigma; 00174 } 00175 00176 // Project to get 2x2 A from A_hat (which might be equal to J) 00177 A = transpose( Theta ) * J; 00178 return true; 00179 } 00180 00181 /* Do transform M_hat = S_a M_{3x2}, M_{2x2} Theta^-1 M_hat 00182 * where the plane into which we are projecting is the cross 00183 * product of the columns of M, such that S_a is I. Use the 00184 * first column of M as u_perp. 00185 * 00186 * Also pass back the cross product of the columns of M as u, 00187 * and the first column of M as u_perp, both normalized. 00188 */ 00189 static inline void project_to_matrix_plane( const MsqMatrix< 3, 2 >& M_in, MsqMatrix< 2, 2 >& M_out, MsqVector< 3 >& u, 00190 MsqVector< 3 >& u_perp ) 00191 { 00192 u = M_in.column( 0 ) * M_in.column( 1 ); 00193 u_perp = M_in.column( 0 ); 00194 double len0 = length( u_perp ); 00195 double u_len = length( u ); 00196 double d_perp, d_u; 00197 if( !divide( 1.0, len0, d_perp ) ) 00198 { 00199 // try the other column 00200 u_perp = M_in.column( 1 ); 00201 len0 = length( u_perp ); 00202 if( !divide( 1.0, len0, d_perp ) ) 00203 { 00204 // matrix is all zeros 00205 u[0] = 0; 00206 u[1] = 0; 00207 u[2] = 1; 00208 u_perp[0] = 1; 00209 u_perp[1] = 0; 00210 u_perp[2] = 0; 00211 M_out = MsqMatrix< 2, 2 >( 0.0 ); 00212 } 00213 else 00214 { 00215 MsqMatrix< 3, 2 > junk; 00216 get_u_perp( u_perp, u ); 00217 project_to_perp_plane( M_in, u, u_perp, M_out, junk ); 00218 } 00219 } 00220 else if( !divide( 1.0, u_len, d_u ) ) 00221 { 00222 MsqMatrix< 3, 2 > junk; 00223 get_u_perp( u_perp, u ); 00224 project_to_perp_plane( M_in, u, u_perp, M_out, junk ); 00225 } 00226 else 00227 { // the normal case (neither column is zero) 00228 u *= d_u; 00229 u_perp *= d_perp; 00230 00231 // M_out = transpose(theta)*M_in 00232 M_out( 0, 0 ) = len0; 00233 M_out( 0, 1 ) = u_perp % M_in.column( 1 ); 00234 M_out( 1, 0 ) = 0.0; 00235 M_out( 1, 1 ) = u_len / len0; 00236 } 00237 } 00238 00239 bool TMPQualityMetric::evaluate_surface_common( PatchData& pd, Sample s, size_t e, const NodeSet& bits, size_t* indices, 00240 size_t& num_indices, MsqVector< 2 >* derivs, MsqMatrix< 2, 2 >& W, 00241 MsqMatrix< 2, 2 >& A, MsqMatrix< 3, 2 >& S_a_transpose_Theta, 00242 MsqError& err ) 00243 { 00244 EntityTopology type = pd.element_by_index( e ).get_element_type(); 00245 00246 const MappingFunction2D* mf = pd.get_mapping_function_2D( type ); 00247 if( !mf ) 00248 { 00249 MSQ_SETERR( err )( "No mapping function for element type", MsqError::UNSUPPORTED_ELEMENT ); 00250 return false; 00251 } 00252 00253 MsqMatrix< 3, 2 > J; 00254 mf->jacobian( pd, e, bits, s, indices, derivs, num_indices, J, err ); 00255 00256 // If we have a 3x2 target matrix 00257 if( targetCalc->have_surface_orient() ) 00258 { 00259 MsqVector< 3 > u, u_perp; 00260 MsqMatrix< 3, 2 > W_hat; 00261 targetCalc->get_surface_target( pd, e, s, W_hat, err ); 00262 MSQ_ERRZERO( err ); 00263 // Use the cross product of the columns of W as the normal of the 00264 // plane to work in (i.e. u.). W should have been constructed such 00265 // that said cross product is in the direction of (n_s)_init. And if 00266 // for some reason it as not, then using something other than said 00267 // cross product is likely to produce very wrong results. 00268 project_to_matrix_plane( W_hat, W, u, u_perp ); 00269 // Do the transforms on A to align it with W and project into the plane. 00270 if( !project_to_perp_plane( J, u, u_perp, A, S_a_transpose_Theta ) ) return false; 00271 } 00272 // Otherwise if we have a 2x2 target matrix (i.e. the target does 00273 // not contain orientation information), project into the plane 00274 // tangent to J. 00275 else 00276 { 00277 MsqVector< 3 > u, u_perp; 00278 targetCalc->get_2D_target( pd, e, s, W, err ); 00279 MSQ_ERRZERO( err ); 00280 project_to_matrix_plane( J, A, u, u_perp ); 00281 S_a_transpose_Theta.set_column( 0, u_perp ); 00282 S_a_transpose_Theta.set_column( 1, u * u_perp ); 00283 // If the domain is set, adjust the sign of things correctly 00284 // for the case where the element is inverted with respect 00285 // to the domain. 00286 if( pd.domain_set() ) 00287 { 00288 Vector3D n; 00289 pd.get_domain_normal_at_sample( e, s, n, err ); 00290 MSQ_ERRZERO( err ); 00291 // if sigma == -1 00292 if( Vector3D( u.data() ) % n < 0.0 ) 00293 { 00294 // flip u 00295 u = -u; 00296 // S_a_transpose_Theta == Theta, because S_a == I here. 00297 // u_perp is unaffected by flipping u, so only the second 00298 // column of S_a_transpose_Theta and the second row of A 00299 // are flipped because u x u_perp will be flipped. 00300 S_a_transpose_Theta.set_column( 1, -S_a_transpose_Theta.column( 1 ) ); 00301 A.set_row( 1, -A.row( 1 ) ); 00302 } 00303 } 00304 } 00305 00306 return true; 00307 } 00308 00309 void TMPQualityMetric::weight( PatchData& pd, Sample p_sample, size_t p_elem, int num_idx, double& value, 00310 Vector3D* grad, SymMatrix3D* diag, Matrix3D* hess, MsqError& err ) 00311 { 00312 if( !weightCalc ) return; 00313 00314 double ck = weightCalc->get_weight( pd, p_elem, p_sample, err );MSQ_ERRRTN( err ); 00315 value *= ck; 00316 if( grad ) 00317 { 00318 for( int i = 0; i < num_idx; ++i ) 00319 grad[i] *= ck; 00320 } 00321 if( diag ) 00322 { 00323 for( int i = 0; i < num_idx; ++i ) 00324 diag[i] *= ck; 00325 } 00326 if( hess ) 00327 { 00328 const int n = num_idx * ( num_idx + 1 ) / 2; 00329 for( int i = 0; i < n; ++i ) 00330 hess[i] *= ck; 00331 } 00332 } 00333 00334 void TMPQualityMetric::initialize_queue( MeshDomainAssoc* mesh_and_domain, const Settings* settings, MsqError& err ) 00335 { 00336 targetCalc->initialize_queue( mesh_and_domain, settings, err );MSQ_ERRRTN( err ); 00337 if( weightCalc ) 00338 { 00339 weightCalc->initialize_queue( mesh_and_domain, settings, err );MSQ_ERRRTN( err ); 00340 } 00341 } 00342 00343 } // namespace MBMesquite