MOAB: Mesh Oriented datABase
(version 5.2.1)
|
00001 /* ***************************************************************** 00002 MESQUITE -- The Mesh Quality Improvement Toolkit 00003 00004 Copyright 2009 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 (2009) kraftche@cae.wisc.edu 00024 00025 ***************************************************************** */ 00026 00027 /** \file TargetCalculator.cpp 00028 * \brief 00029 * \author Jason Kraftcheck 00030 */ 00031 00032 #include "Mesquite.hpp" 00033 #include "TargetCalculator.hpp" 00034 #include "MsqError.hpp" 00035 #include "PatchData.hpp" 00036 #include "ReferenceMesh.hpp" 00037 #include "MappingFunction.hpp" 00038 #include <assert.h> 00039 00040 namespace MBMesquite 00041 { 00042 00043 double TargetCalculator::size( const MsqMatrix< 3, 3 >& W ) 00044 { 00045 return MBMesquite::cbrt( fabs( det( W ) ) ); 00046 } 00047 00048 double TargetCalculator::size( const MsqMatrix< 3, 2 >& W ) 00049 { 00050 return sqrt( length( W.column( 0 ) * W.column( 1 ) ) ); 00051 } 00052 00053 double TargetCalculator::size( const MsqMatrix< 2, 2 >& W ) 00054 { 00055 return sqrt( fabs( det( W ) ) ); 00056 } 00057 00058 MsqMatrix< 3, 3 > TargetCalculator::skew( const MsqMatrix< 3, 3 >& W ) 00059 { 00060 MsqVector< 3 > a1 = W.column( 0 ) * ( 1.0 / length( W.column( 0 ) ) ); 00061 MsqVector< 3 > a2 = W.column( 1 ) * ( 1.0 / length( W.column( 1 ) ) ); 00062 MsqVector< 3 > a3 = W.column( 2 ) * ( 1.0 / length( W.column( 2 ) ) ); 00063 MsqVector< 3 > a1xa2 = a1 * a2; 00064 00065 double lenx = length( a1xa2 ); 00066 double alpha = fabs( a1xa2 % a3 ); 00067 double coeff = MBMesquite::cbrt( 1 / alpha ); 00068 double dot = a1xa2 % ( a1 * a3 ); 00069 00070 MsqMatrix< 3, 3 > q; 00071 q( 0, 0 ) = coeff; 00072 q( 0, 1 ) = coeff * ( a1 % a2 ); 00073 q( 0, 2 ) = coeff * ( a1 % a3 ); 00074 q( 1, 0 ) = 0.0; 00075 q( 1, 1 ) = coeff * lenx; 00076 q( 1, 2 ) = coeff * dot / lenx; 00077 q( 2, 0 ) = 0.0; 00078 q( 2, 1 ) = 0.0; 00079 q( 2, 2 ) = coeff * alpha / lenx; 00080 return q; 00081 } 00082 00083 MsqMatrix< 2, 2 > TargetCalculator::skew( const MsqMatrix< 3, 2 >& W ) 00084 { 00085 MsqVector< 3 > alpha = W.column( 0 ) * W.column( 1 ); 00086 double a1_sqr = W.column( 0 ) % W.column( 0 ); 00087 double a2_sqr = W.column( 1 ) % W.column( 1 ); 00088 double dot = W.column( 0 ) % W.column( 1 ); 00089 double a1a2 = sqrt( a1_sqr * a2_sqr ); 00090 double coeff = sqrt( a1a2 / length( alpha ) ); 00091 00092 MsqMatrix< 2, 2 > result; 00093 result( 0, 0 ) = coeff; 00094 result( 0, 1 ) = coeff * dot / a1a2; 00095 result( 1, 0 ) = 0.0; 00096 result( 1, 1 ) = 1 / coeff; 00097 return result; 00098 } 00099 00100 MsqMatrix< 2, 2 > TargetCalculator::skew( const MsqMatrix< 2, 2 >& W ) 00101 { 00102 double alpha = fabs( det( W ) ); 00103 double a1_sqr = W.column( 0 ) % W.column( 0 ); 00104 double a2_sqr = W.column( 1 ) % W.column( 1 ); 00105 double dot = W.column( 0 ) % W.column( 1 ); 00106 double a1a2 = sqrt( a1_sqr * a2_sqr ); 00107 double coeff = sqrt( a1a2 / alpha ); 00108 00109 MsqMatrix< 2, 2 > result; 00110 result( 0, 0 ) = coeff; 00111 result( 0, 1 ) = coeff * dot / a1a2; 00112 result( 1, 0 ) = 0.0; 00113 result( 1, 1 ) = 1 / coeff; 00114 return result; 00115 } 00116 00117 MsqMatrix< 3, 3 > TargetCalculator::aspect( const MsqMatrix< 3, 3 >& W ) 00118 { 00119 double a1 = length( W.column( 0 ) ); 00120 double a2 = length( W.column( 1 ) ); 00121 double a3 = length( W.column( 2 ) ); 00122 double coeff = 1.0 / MBMesquite::cbrt( a1 * a2 * a3 ); 00123 00124 MsqMatrix< 3, 3 > result( coeff ); 00125 result( 0, 0 ) *= a1; 00126 result( 1, 1 ) *= a2; 00127 result( 2, 2 ) *= a3; 00128 return result; 00129 } 00130 00131 MsqMatrix< 2, 2 > TargetCalculator::aspect( const MsqMatrix< 3, 2 >& W ) 00132 { 00133 double a1_sqr = W.column( 0 ) % W.column( 0 ); 00134 double a2_sqr = W.column( 1 ) % W.column( 1 ); 00135 double sqrt_rho = sqrt( sqrt( a1_sqr / a2_sqr ) ); 00136 00137 MsqMatrix< 2, 2 > result; 00138 result( 0, 0 ) = sqrt_rho; 00139 result( 0, 1 ) = 0.0; 00140 result( 1, 0 ) = 0.0; 00141 result( 1, 1 ) = 1 / sqrt_rho; 00142 return result; 00143 } 00144 00145 MsqMatrix< 2, 2 > TargetCalculator::aspect( const MsqMatrix< 2, 2 >& W ) 00146 { 00147 double a1_sqr = W.column( 0 ) % W.column( 0 ); 00148 double a2_sqr = W.column( 1 ) % W.column( 1 ); 00149 double sqrt_rho = sqrt( sqrt( a1_sqr / a2_sqr ) ); 00150 00151 MsqMatrix< 2, 2 > result; 00152 result( 0, 0 ) = sqrt_rho; 00153 result( 0, 1 ) = 0.0; 00154 result( 1, 0 ) = 0.0; 00155 result( 1, 1 ) = 1 / sqrt_rho; 00156 return result; 00157 } 00158 00159 MsqMatrix< 3, 3 > TargetCalculator::shape( const MsqMatrix< 3, 3 >& W ) 00160 { 00161 MsqVector< 3 > a1 = W.column( 0 ); 00162 MsqVector< 3 > a2 = W.column( 1 ); 00163 MsqVector< 3 > a3 = W.column( 2 ); 00164 MsqVector< 3 > a1xa2 = a1 * a2; 00165 00166 double len1 = length( a1 ); 00167 double lenx = length( a1xa2 ); 00168 double alpha = fabs( a1xa2 % a3 ); 00169 double coeff = MBMesquite::cbrt( 1 / alpha ); 00170 double inv1 = 1.0 / len1; 00171 double invx = 1.0 / lenx; 00172 00173 MsqMatrix< 3, 3 > q; 00174 q( 0, 0 ) = coeff * len1; 00175 q( 0, 1 ) = coeff * inv1 * ( a1 % a2 ); 00176 q( 0, 2 ) = coeff * inv1 * ( a1 % a3 ); 00177 q( 1, 0 ) = 0.0; 00178 q( 1, 1 ) = coeff * inv1 * lenx; 00179 q( 1, 2 ) = coeff * invx * inv1 * ( a1xa2 % ( a1 * a3 ) ); 00180 q( 2, 0 ) = 0.0; 00181 q( 2, 1 ) = 0.0; 00182 q( 2, 2 ) = coeff * alpha * invx; 00183 return q; 00184 } 00185 00186 MsqMatrix< 2, 2 > TargetCalculator::shape( const MsqMatrix< 3, 2 >& W ) 00187 { 00188 double len1 = length( W.column( 0 ) ); 00189 double inv1 = 1.0 / len1; 00190 double root_alpha = sqrt( length( W.column( 0 ) * W.column( 1 ) ) ); 00191 double coeff = 1.0 / root_alpha; 00192 00193 MsqMatrix< 2, 2 > result; 00194 result( 0, 0 ) = coeff * len1; 00195 result( 0, 1 ) = coeff * inv1 * ( W.column( 0 ) % W.column( 1 ) ); 00196 result( 1, 0 ) = 0.0; 00197 result( 1, 1 ) = root_alpha * inv1; 00198 return result; 00199 } 00200 00201 MsqMatrix< 2, 2 > TargetCalculator::shape( const MsqMatrix< 2, 2 >& W ) 00202 { 00203 double len1 = length( W.column( 0 ) ); 00204 double inv1 = 1.0 / len1; 00205 double root_alpha = sqrt( fabs( det( W ) ) ); 00206 double coeff = 1.0 / root_alpha; 00207 00208 MsqMatrix< 2, 2 > result; 00209 result( 0, 0 ) = coeff * len1; 00210 result( 0, 1 ) = coeff * inv1 * ( W.column( 0 ) % W.column( 1 ) ); 00211 result( 1, 0 ) = 0.0; 00212 result( 1, 1 ) = root_alpha * inv1; 00213 return result; 00214 } 00215 00216 bool TargetCalculator::factor_3D( const MsqMatrix< 3, 3 >& A, double& Lambda, MsqMatrix< 3, 3 >& V, 00217 MsqMatrix< 3, 3 >& Q, MsqMatrix< 3, 3 >& Delta, MsqError& /*err*/ ) 00218 { 00219 MsqVector< 3 > a1xa2 = A.column( 0 ) * A.column( 1 ); 00220 double alpha = a1xa2 % A.column( 2 ); 00221 Lambda = MBMesquite::cbrt( fabs( alpha ) ); 00222 if( Lambda < DBL_EPSILON ) return false; 00223 00224 double la1_sqr = A.column( 0 ) % A.column( 0 ); 00225 double la1 = sqrt( la1_sqr ); 00226 double la2 = length( A.column( 1 ) ); 00227 double la3 = length( A.column( 2 ) ); 00228 double lx = length( a1xa2 ); 00229 double a1dota2 = A.column( 0 ) % A.column( 1 ); 00230 00231 double inv_la1 = 1.0 / la1; 00232 double inv_lx = 1.0 / lx; 00233 V.set_column( 0, A.column( 0 ) * inv_la1 ); 00234 V.set_column( 1, ( la1_sqr * A.column( 1 ) - a1dota2 * A.column( 0 ) ) * inv_la1 * inv_lx ); 00235 V.set_column( 2, ( alpha / ( fabs( alpha ) * lx ) ) * a1xa2 ); 00236 00237 double inv_la2 = 1.0 / la2; 00238 double inv_la3 = 1.0 / la3; 00239 double len_prod_rt3 = MBMesquite::cbrt( la1 * la2 * la3 ); 00240 Q( 0, 0 ) = 1.0; 00241 Q( 1, 0 ) = Q( 2, 0 ) = 0.0; 00242 Q( 0, 1 ) = a1dota2 * inv_la1 * inv_la2; 00243 Q( 1, 1 ) = lx * inv_la1 * inv_la2; 00244 Q( 2, 1 ) = 0.0; 00245 Q( 0, 2 ) = ( A.column( 0 ) % A.column( 2 ) ) * inv_la1 * inv_la3; 00246 Q( 1, 2 ) = ( a1xa2 % ( A.column( 0 ) * A.column( 2 ) ) ) * inv_lx * inv_la1 * inv_la3; 00247 Q( 2, 2 ) = fabs( alpha ) * inv_lx * inv_la3; 00248 Q *= len_prod_rt3 / Lambda; 00249 00250 double inv_prod_rt3 = 1.0 / len_prod_rt3; 00251 ; 00252 Delta( 0, 0 ) = la1 * inv_prod_rt3; 00253 Delta( 0, 1 ) = 0.0; 00254 Delta( 0, 2 ) = 0.0; 00255 Delta( 1, 0 ) = 0.0; 00256 Delta( 1, 1 ) = la2 * inv_prod_rt3; 00257 Delta( 1, 2 ) = 0.0; 00258 Delta( 2, 0 ) = 0.0; 00259 Delta( 2, 1 ) = 0.0; 00260 Delta( 2, 2 ) = la3 * inv_prod_rt3; 00261 00262 return true; 00263 } 00264 00265 bool TargetCalculator::factor_surface( const MsqMatrix< 3, 2 >& A, double& Lambda, MsqMatrix< 3, 2 >& V, 00266 MsqMatrix< 2, 2 >& Q, MsqMatrix< 2, 2 >& Delta, MsqError& /*err*/ ) 00267 { 00268 MsqVector< 3 > cross = A.column( 0 ) * A.column( 1 ); 00269 double alpha = length( cross ); 00270 Lambda = sqrt( alpha ); 00271 if( Lambda < DBL_EPSILON ) return false; 00272 00273 double la1_sqr = A.column( 0 ) % A.column( 0 ); 00274 double la1 = sqrt( la1_sqr ); 00275 double la2 = length( A.column( 1 ) ); 00276 double inv_la1 = 1.0 / la1; 00277 double dot = A.column( 0 ) % A.column( 1 ); 00278 00279 V.set_column( 0, A.column( 0 ) * inv_la1 ); 00280 V.set_column( 1, ( la1_sqr * A.column( 1 ) - dot * A.column( 0 ) ) / ( la1 * alpha ) ); 00281 00282 double prod_rt2 = sqrt( la1 * la2 ); 00283 Q( 0, 0 ) = prod_rt2 / Lambda; 00284 Q( 0, 1 ) = dot / ( prod_rt2 * Lambda ); 00285 Q( 1, 0 ) = 0.0; 00286 Q( 1, 1 ) = 1.0 / Q( 0, 0 ); 00287 00288 double inv_prod_rt2 = 1.0 / prod_rt2; 00289 Delta( 0, 0 ) = la1 * inv_prod_rt2; 00290 Delta( 0, 1 ) = 0.0; 00291 Delta( 1, 0 ) = 0.0; 00292 Delta( 1, 1 ) = la2 * inv_prod_rt2; 00293 00294 return true; 00295 } 00296 00297 bool TargetCalculator::factor_2D( const MsqMatrix< 2, 2 >& A, double& Lambda, MsqMatrix< 2, 2 >& V, 00298 MsqMatrix< 2, 2 >& Q, MsqMatrix< 2, 2 >& Delta, MsqError& /*err*/ ) 00299 { 00300 double alpha = det( A ); 00301 Lambda = sqrt( fabs( alpha ) ); 00302 if( Lambda < DBL_EPSILON ) return false; 00303 00304 double la1_sqr = A.column( 0 ) % A.column( 0 ); 00305 double la1 = sqrt( la1_sqr ); 00306 double la2 = length( A.column( 1 ) ); 00307 double inv_la1 = 1.0 / la1; 00308 double dot = A.column( 0 ) % A.column( 1 ); 00309 00310 V.set_column( 0, A.column( 0 ) * inv_la1 ); 00311 V.set_column( 1, ( la1_sqr * A.column( 1 ) - dot * A.column( 0 ) ) / ( la1 * alpha ) ); 00312 00313 double prod_rt2 = sqrt( la1 * la2 ); 00314 Q( 0, 0 ) = prod_rt2 / Lambda; 00315 Q( 0, 1 ) = dot / ( prod_rt2 * Lambda ); 00316 Q( 1, 0 ) = 0.0; 00317 Q( 1, 1 ) = 1.0 / Q( 0, 0 ); 00318 00319 double inv_prod_rt2 = 1.0 / prod_rt2; 00320 Delta( 0, 0 ) = la1 * inv_prod_rt2; 00321 Delta( 0, 1 ) = 0.0; 00322 Delta( 1, 0 ) = 0.0; 00323 Delta( 1, 1 ) = la2 * inv_prod_rt2; 00324 00325 return true; 00326 } 00327 00328 MsqMatrix< 3, 3 > TargetCalculator::new_orientation_3D( const MsqVector< 3 >& b1, const MsqVector< 3 >& b2 ) 00329 { 00330 double lb1_sqr = b1 % b1; 00331 MsqVector< 3 > cross = b1 * b2; 00332 double lb1 = sqrt( lb1_sqr ); 00333 double inv_lb1 = 1.0 / lb1; 00334 double inv_lx = 1.0 / length( cross ); 00335 MsqMatrix< 3, 3 > V; 00336 V.set_column( 0, inv_lb1 * b1 ); 00337 V.set_column( 1, ( lb1_sqr * b2 - ( b1 % b2 ) * lb1 ) * inv_lb1 * inv_lx ); 00338 V.set_column( 2, cross * inv_lx ); 00339 return V; 00340 } 00341 00342 MsqMatrix< 3, 2 > TargetCalculator::new_orientation_2D( const MsqVector< 3 >& b1, const MsqVector< 3 >& b2 ) 00343 { 00344 double lb1_sqr = b1 % b1; 00345 double inv_lb1 = 1.0 / sqrt( lb1_sqr ); 00346 MsqMatrix< 3, 2 > V; 00347 V.set_column( 0, b1 * inv_lb1 ); 00348 V.set_column( 1, ( lb1_sqr * b2 - ( b1 % b2 ) * b1 ) * ( inv_lb1 / length( b1 * b2 ) ) ); 00349 return V; 00350 } 00351 00352 /** If, for the specified element type, the skew is constant for 00353 * an ideal element and the aspect is identity everywhere within 00354 * the element, pass back the constant skew/shape term and return 00355 * true. Otherwise return false. 00356 */ 00357 static inline bool ideal_constant_skew_I_3D( EntityTopology element_type, MsqMatrix< 3, 3 >& q ) 00358 { 00359 switch( element_type ) 00360 { 00361 case HEXAHEDRON: 00362 q = MsqMatrix< 3, 3 >( 1.0 ); // Identity 00363 return false; 00364 case TETRAHEDRON: 00365 // [ x, x/2, x/2 ] x^6 = 2 00366 // [ 0, y, y/3 ] y^2 = 3/4 x^2 00367 // [ 0, 0, z ] z^2 = 2/3 x^2 00368 q( 0, 0 ) = 1.122462048309373; 00369 q( 0, 1 ) = q( 0, 2 ) = 0.56123102415468651; 00370 q( 1, 0 ) = q( 2, 0 ) = q( 2, 1 ) = 0.0; 00371 q( 1, 1 ) = 0.97208064861983279; 00372 q( 1, 2 ) = 0.32402688287327758; 00373 q( 2, 2 ) = 0.91648642466573493; 00374 return true; 00375 case PRISM: 00376 // [ 1 0 0 ] 00377 // a^(-1/3) [ 0 1 1/2] 00378 // [ 0 0 a ] 00379 // 00380 // a = sqrt(3)/2 00381 // 00382 q( 0, 0 ) = q( 1, 1 ) = 1.0491150634216482; 00383 q( 0, 1 ) = q( 0, 2 ) = q( 1, 0 ) = 0.0; 00384 q( 1, 2 ) = 0.52455753171082409; 00385 q( 2, 0 ) = q( 2, 1 ) = 0.0; 00386 q( 2, 2 ) = 0.90856029641606972; 00387 return true; 00388 default: 00389 return false; 00390 } 00391 } 00392 00393 /** If, for the specified element type, the skew is constant for 00394 * an ideal element and the aspect is identity everywhere within 00395 * the element, pass back the constant skew/shape term and return 00396 * true. Otherwise return false. 00397 */ 00398 static inline bool ideal_constant_skew_I_2D( EntityTopology element_type, MsqMatrix< 2, 2 >& q ) 00399 { 00400 switch( element_type ) 00401 { 00402 case QUADRILATERAL: 00403 q = MsqMatrix< 2, 2 >( 1.0 ); // Identity 00404 return true; 00405 case TRIANGLE: 00406 // [ x, x/2 ] x = pow(4/3, 0.25) 00407 // [ 0, y ] y = 1/x 00408 q( 0, 0 ) = 1.074569931823542; 00409 q( 0, 1 ) = 0.537284965911771; 00410 q( 1, 0 ) = 0.0; 00411 q( 1, 1 ) = 0.93060485910209956; 00412 return true; 00413 default: 00414 return false; 00415 } 00416 } 00417 00418 void TargetCalculator::ideal_skew_3D( EntityTopology element_type, Sample s, const PatchData& pd, MsqMatrix< 3, 3 >& q, 00419 MsqError& err ) 00420 { 00421 if( !ideal_constant_skew_I_3D( element_type, q ) ) 00422 { 00423 const MappingFunction3D* map = pd.get_mapping_function_3D( element_type ); 00424 if( !map ) 00425 { 00426 MSQ_SETERR( err )( MsqError::UNSUPPORTED_ELEMENT ); 00427 return; 00428 } 00429 map->ideal( s, q, err );MSQ_ERRRTN( err ); 00430 q = TargetCalculator::skew( q ); 00431 } 00432 } 00433 00434 void TargetCalculator::ideal_skew_2D( EntityTopology element_type, Sample s, const PatchData& pd, MsqMatrix< 2, 2 >& q, 00435 MsqError& err ) 00436 { 00437 if( !ideal_constant_skew_I_2D( element_type, q ) ) 00438 { 00439 const MappingFunction2D* map = pd.get_mapping_function_2D( element_type ); 00440 if( !map ) 00441 { 00442 MSQ_SETERR( err )( MsqError::UNSUPPORTED_ELEMENT ); 00443 return; 00444 } 00445 MsqMatrix< 3, 2 > J; 00446 map->ideal( s, J, err );MSQ_ERRRTN( err ); 00447 q = TargetCalculator::skew( J ); 00448 } 00449 } 00450 00451 void TargetCalculator::ideal_shape_3D( EntityTopology element_type, Sample s, const PatchData& pd, MsqMatrix< 3, 3 >& q, 00452 MsqError& err ) 00453 { 00454 if( !ideal_constant_skew_I_3D( element_type, q ) ) 00455 { 00456 const MappingFunction3D* map = pd.get_mapping_function_3D( element_type ); 00457 if( !map ) 00458 { 00459 MSQ_SETERR( err )( MsqError::UNSUPPORTED_ELEMENT ); 00460 return; 00461 } 00462 map->ideal( s, q, err );MSQ_ERRRTN( err ); 00463 q = TargetCalculator::shape( q ); 00464 } 00465 } 00466 00467 void TargetCalculator::ideal_shape_2D( EntityTopology element_type, Sample s, const PatchData& pd, MsqMatrix< 2, 2 >& q, 00468 MsqError& err ) 00469 { 00470 if( !ideal_constant_skew_I_2D( element_type, q ) ) 00471 { 00472 const MappingFunction2D* map = pd.get_mapping_function_2D( element_type ); 00473 if( !map ) 00474 { 00475 MSQ_SETERR( err )( MsqError::UNSUPPORTED_ELEMENT ); 00476 return; 00477 } 00478 MsqMatrix< 3, 2 > J; 00479 map->ideal( s, J, err );MSQ_ERRRTN( err ); 00480 q = TargetCalculator::shape( J ); 00481 } 00482 } 00483 00484 MsqMatrix< 3, 3 > TargetCalculator::new_aspect_3D( const MsqVector< 3 >& r ) 00485 { 00486 MsqMatrix< 3, 3 > W( 0.0 ); 00487 W( 0, 0 ) = MBMesquite::cbrt( r[0] * r[0] / ( r[1] * r[2] ) ); 00488 W( 1, 1 ) = MBMesquite::cbrt( r[1] * r[1] / ( r[0] * r[2] ) ); 00489 W( 2, 2 ) = MBMesquite::cbrt( r[2] * r[2] / ( r[1] * r[0] ) ); 00490 return W; 00491 } 00492 00493 MsqMatrix< 2, 2 > TargetCalculator::new_aspect_2D( const MsqVector< 2 >& r ) 00494 { 00495 return new_aspect_2D( r[0] / r[1] ); 00496 } 00497 00498 MsqMatrix< 2, 2 > TargetCalculator::new_aspect_2D( double rho ) 00499 { 00500 MsqMatrix< 2, 2 > W( sqrt( rho ) ); 00501 W( 1, 1 ) = 1.0 / W( 0, 0 ); 00502 return W; 00503 } 00504 00505 static NodeSet get_nodeset( EntityTopology type, int num_nodes, MsqError& err ) 00506 { 00507 bool midedge, midface, midvol; 00508 TopologyInfo::higher_order( type, num_nodes, midedge, midface, midvol, err ); 00509 if( MSQ_CHKERR( err ) ) return NodeSet(); 00510 00511 NodeSet bits; 00512 bits.set_all_corner_nodes( type ); 00513 if( midedge ) bits.set_all_mid_edge_nodes( type ); 00514 if( midface ) bits.set_all_mid_face_nodes( type ); 00515 if( TopologyInfo::dimension( type ) == 3 && midvol ) bits.set_mid_region_node(); 00516 00517 return bits; 00518 } 00519 00520 void TargetCalculator::jacobian_3D( PatchData& pd, EntityTopology type, int num_nodes, Sample location, 00521 const Vector3D* coords, MsqMatrix< 3, 3 >& J, MsqError& err ) 00522 { 00523 // Get element properties 00524 NodeSet bits = get_nodeset( type, num_nodes, err );MSQ_ERRRTN( err ); 00525 const MappingFunction3D* mf = pd.get_mapping_function_3D( type ); 00526 if( !mf ) 00527 { 00528 MSQ_SETERR( err )( MsqError::UNSUPPORTED_ELEMENT ); 00529 return; 00530 } 00531 00532 // Get mapping function derivatives 00533 const int MAX_NODES = 27; 00534 assert( num_nodes <= MAX_NODES ); 00535 size_t indices[MAX_NODES], n; 00536 MsqVector< 3 > derivs[MAX_NODES]; 00537 mf->derivatives( location, bits, indices, derivs, n, err );MSQ_ERRRTN( err ); 00538 00539 // calculate Jacobian 00540 assert( sizeof( Vector3D ) == sizeof( MsqVector< 3 > ) ); 00541 const MsqVector< 3 >* verts = reinterpret_cast< const MsqVector< 3 >* >( coords ); 00542 assert( n > 0 ); 00543 J = outer( verts[indices[0]], derivs[0] ); 00544 for( size_t i = 1; i < n; ++i ) 00545 J += outer( verts[indices[i]], derivs[i] ); 00546 } 00547 00548 void TargetCalculator::jacobian_2D( PatchData& pd, EntityTopology type, int num_nodes, Sample location, 00549 const Vector3D* coords, MsqMatrix< 3, 2 >& J, MsqError& err ) 00550 { 00551 // Get element properties 00552 NodeSet bits = get_nodeset( type, num_nodes, err );MSQ_ERRRTN( err ); 00553 const MappingFunction2D* mf = pd.get_mapping_function_2D( type ); 00554 if( !mf ) 00555 { 00556 MSQ_SETERR( err )( MsqError::UNSUPPORTED_ELEMENT ); 00557 return; 00558 } 00559 00560 // Get mapping function derivatives 00561 const int MAX_NODES = 9; 00562 assert( num_nodes <= MAX_NODES ); 00563 size_t indices[MAX_NODES], n; 00564 MsqVector< 2 > derivs[MAX_NODES]; 00565 mf->derivatives( location, bits, indices, derivs, n, err );MSQ_ERRRTN( err ); 00566 00567 // calculate Jacobian 00568 assert( sizeof( Vector3D ) == sizeof( MsqVector< 3 > ) ); 00569 const MsqVector< 3 >* verts = reinterpret_cast< const MsqVector< 3 >* >( coords ); 00570 assert( n > 0 ); 00571 J = outer( verts[indices[0]], derivs[0] ); 00572 for( size_t i = 1; i < n; ++i ) 00573 J += outer( verts[indices[i]], derivs[i] ); 00574 } 00575 00576 void TargetCalculator::get_refmesh_Jacobian_3D( ReferenceMeshInterface* ref_mesh, PatchData& pd, size_t element, 00577 Sample sample, MsqMatrix< 3, 3 >& W_out, MsqError& err ) 00578 { 00579 // get element 00580 MsqMeshEntity& elem = pd.element_by_index( element ); 00581 const EntityTopology type = elem.get_element_type(); 00582 const unsigned n = elem.node_count(); 00583 00584 const unsigned MAX_NODES = 27; 00585 assert( n <= MAX_NODES ); 00586 00587 // get vertices 00588 Mesh::VertexHandle elem_verts[MAX_NODES]; 00589 const std::size_t* vtx_idx = elem.get_vertex_index_array(); 00590 const Mesh::VertexHandle* vtx_hdl = pd.get_vertex_handles_array(); 00591 for( unsigned i = 0; i < n; ++i ) 00592 elem_verts[i] = vtx_hdl[vtx_idx[i]]; 00593 00594 // get vertex coordinates 00595 Vector3D vert_coords[MAX_NODES]; 00596 ref_mesh->get_reference_vertex_coordinates( elem_verts, n, vert_coords, err );MSQ_ERRRTN( err ); 00597 00598 // calculate Jacobian 00599 jacobian_3D( pd, type, n, sample, vert_coords, W_out, err );MSQ_ERRRTN( err ); 00600 } 00601 00602 void TargetCalculator::get_refmesh_Jacobian_2D( ReferenceMeshInterface* ref_mesh, PatchData& pd, size_t element, 00603 Sample sample, MsqMatrix< 3, 2 >& W_out, MsqError& err ) 00604 { 00605 // get element 00606 MsqMeshEntity& elem = pd.element_by_index( element ); 00607 const EntityTopology type = elem.get_element_type(); 00608 const unsigned n = elem.node_count(); 00609 00610 const unsigned MAX_NODES = 9; 00611 assert( n <= MAX_NODES ); 00612 00613 // get vertices 00614 Mesh::VertexHandle elem_verts[MAX_NODES]; 00615 const std::size_t* vtx_idx = elem.get_vertex_index_array(); 00616 const Mesh::VertexHandle* vtx_hdl = pd.get_vertex_handles_array(); 00617 for( unsigned i = 0; i < n; ++i ) 00618 elem_verts[i] = vtx_hdl[vtx_idx[i]]; 00619 00620 // get vertex coordinates 00621 Vector3D vert_coords[MAX_NODES]; 00622 ref_mesh->get_reference_vertex_coordinates( elem_verts, n, vert_coords, err );MSQ_ERRRTN( err ); 00623 00624 // calculate Jacobian 00625 jacobian_2D( pd, type, n, sample, vert_coords, W_out, err );MSQ_ERRRTN( err ); 00626 } 00627 00628 TargetCalculator::~TargetCalculator() {} 00629 00630 void TargetCalculator::initialize_queue( MeshDomainAssoc*, const Settings*, MsqError& ) {} 00631 00632 } // namespace MBMesquite