MOAB: Mesh Oriented datABase  (version 5.2.1)
TargetCalculator.cpp
Go to the documentation of this file.
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 <cassert>
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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines