MOAB: Mesh Oriented datABase  (version 5.4.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) [email protected]
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,
00217                                   double& Lambda,
00218                                   MsqMatrix< 3, 3 >& V,
00219                                   MsqMatrix< 3, 3 >& Q,
00220                                   MsqMatrix< 3, 3 >& Delta,
00221                                   MsqError& /*err*/ )
00222 {
00223     MsqVector< 3 > a1xa2 = A.column( 0 ) * A.column( 1 );
00224     double alpha         = a1xa2 % A.column( 2 );
00225     Lambda               = MBMesquite::cbrt( fabs( alpha ) );
00226     if( Lambda < DBL_EPSILON ) return false;
00227 
00228     double la1_sqr = A.column( 0 ) % A.column( 0 );
00229     double la1     = sqrt( la1_sqr );
00230     double la2     = length( A.column( 1 ) );
00231     double la3     = length( A.column( 2 ) );
00232     double lx      = length( a1xa2 );
00233     double a1dota2 = A.column( 0 ) % A.column( 1 );
00234 
00235     double inv_la1 = 1.0 / la1;
00236     double inv_lx  = 1.0 / lx;
00237     V.set_column( 0, A.column( 0 ) * inv_la1 );
00238     V.set_column( 1, ( la1_sqr * A.column( 1 ) - a1dota2 * A.column( 0 ) ) * inv_la1 * inv_lx );
00239     V.set_column( 2, ( alpha / ( fabs( alpha ) * lx ) ) * a1xa2 );
00240 
00241     double inv_la2      = 1.0 / la2;
00242     double inv_la3      = 1.0 / la3;
00243     double len_prod_rt3 = MBMesquite::cbrt( la1 * la2 * la3 );
00244     Q( 0, 0 )           = 1.0;
00245     Q( 1, 0 ) = Q( 2, 0 ) = 0.0;
00246     Q( 0, 1 )             = a1dota2 * inv_la1 * inv_la2;
00247     Q( 1, 1 )             = lx * inv_la1 * inv_la2;
00248     Q( 2, 1 )             = 0.0;
00249     Q( 0, 2 )             = ( A.column( 0 ) % A.column( 2 ) ) * inv_la1 * inv_la3;
00250     Q( 1, 2 )             = ( a1xa2 % ( A.column( 0 ) * A.column( 2 ) ) ) * inv_lx * inv_la1 * inv_la3;
00251     Q( 2, 2 )             = fabs( alpha ) * inv_lx * inv_la3;
00252     Q *= len_prod_rt3 / Lambda;
00253 
00254     double inv_prod_rt3 = 1.0 / len_prod_rt3;
00255     ;
00256     Delta( 0, 0 ) = la1 * inv_prod_rt3;
00257     Delta( 0, 1 ) = 0.0;
00258     Delta( 0, 2 ) = 0.0;
00259     Delta( 1, 0 ) = 0.0;
00260     Delta( 1, 1 ) = la2 * inv_prod_rt3;
00261     Delta( 1, 2 ) = 0.0;
00262     Delta( 2, 0 ) = 0.0;
00263     Delta( 2, 1 ) = 0.0;
00264     Delta( 2, 2 ) = la3 * inv_prod_rt3;
00265 
00266     return true;
00267 }
00268 
00269 bool TargetCalculator::factor_surface( const MsqMatrix< 3, 2 >& A,
00270                                        double& Lambda,
00271                                        MsqMatrix< 3, 2 >& V,
00272                                        MsqMatrix< 2, 2 >& Q,
00273                                        MsqMatrix< 2, 2 >& Delta,
00274                                        MsqError& /*err*/ )
00275 {
00276     MsqVector< 3 > cross = A.column( 0 ) * A.column( 1 );
00277     double alpha         = length( cross );
00278     Lambda               = sqrt( alpha );
00279     if( Lambda < DBL_EPSILON ) return false;
00280 
00281     double la1_sqr = A.column( 0 ) % A.column( 0 );
00282     double la1     = sqrt( la1_sqr );
00283     double la2     = length( A.column( 1 ) );
00284     double inv_la1 = 1.0 / la1;
00285     double dot     = A.column( 0 ) % A.column( 1 );
00286 
00287     V.set_column( 0, A.column( 0 ) * inv_la1 );
00288     V.set_column( 1, ( la1_sqr * A.column( 1 ) - dot * A.column( 0 ) ) / ( la1 * alpha ) );
00289 
00290     double prod_rt2 = sqrt( la1 * la2 );
00291     Q( 0, 0 )       = prod_rt2 / Lambda;
00292     Q( 0, 1 )       = dot / ( prod_rt2 * Lambda );
00293     Q( 1, 0 )       = 0.0;
00294     Q( 1, 1 )       = 1.0 / Q( 0, 0 );
00295 
00296     double inv_prod_rt2 = 1.0 / prod_rt2;
00297     Delta( 0, 0 )       = la1 * inv_prod_rt2;
00298     Delta( 0, 1 )       = 0.0;
00299     Delta( 1, 0 )       = 0.0;
00300     Delta( 1, 1 )       = la2 * inv_prod_rt2;
00301 
00302     return true;
00303 }
00304 
00305 bool TargetCalculator::factor_2D( const MsqMatrix< 2, 2 >& A,
00306                                   double& Lambda,
00307                                   MsqMatrix< 2, 2 >& V,
00308                                   MsqMatrix< 2, 2 >& Q,
00309                                   MsqMatrix< 2, 2 >& Delta,
00310                                   MsqError& /*err*/ )
00311 {
00312     double alpha = det( A );
00313     Lambda       = sqrt( fabs( alpha ) );
00314     if( Lambda < DBL_EPSILON ) return false;
00315 
00316     double la1_sqr = A.column( 0 ) % A.column( 0 );
00317     double la1     = sqrt( la1_sqr );
00318     double la2     = length( A.column( 1 ) );
00319     double inv_la1 = 1.0 / la1;
00320     double dot     = A.column( 0 ) % A.column( 1 );
00321 
00322     V.set_column( 0, A.column( 0 ) * inv_la1 );
00323     V.set_column( 1, ( la1_sqr * A.column( 1 ) - dot * A.column( 0 ) ) / ( la1 * alpha ) );
00324 
00325     double prod_rt2 = sqrt( la1 * la2 );
00326     Q( 0, 0 )       = prod_rt2 / Lambda;
00327     Q( 0, 1 )       = dot / ( prod_rt2 * Lambda );
00328     Q( 1, 0 )       = 0.0;
00329     Q( 1, 1 )       = 1.0 / Q( 0, 0 );
00330 
00331     double inv_prod_rt2 = 1.0 / prod_rt2;
00332     Delta( 0, 0 )       = la1 * inv_prod_rt2;
00333     Delta( 0, 1 )       = 0.0;
00334     Delta( 1, 0 )       = 0.0;
00335     Delta( 1, 1 )       = la2 * inv_prod_rt2;
00336 
00337     return true;
00338 }
00339 
00340 MsqMatrix< 3, 3 > TargetCalculator::new_orientation_3D( const MsqVector< 3 >& b1, const MsqVector< 3 >& b2 )
00341 {
00342     double lb1_sqr       = b1 % b1;
00343     MsqVector< 3 > cross = b1 * b2;
00344     double lb1           = sqrt( lb1_sqr );
00345     double inv_lb1       = 1.0 / lb1;
00346     double inv_lx        = 1.0 / length( cross );
00347     MsqMatrix< 3, 3 > V;
00348     V.set_column( 0, inv_lb1 * b1 );
00349     V.set_column( 1, ( lb1_sqr * b2 - ( b1 % b2 ) * lb1 ) * inv_lb1 * inv_lx );
00350     V.set_column( 2, cross * inv_lx );
00351     return V;
00352 }
00353 
00354 MsqMatrix< 3, 2 > TargetCalculator::new_orientation_2D( const MsqVector< 3 >& b1, const MsqVector< 3 >& b2 )
00355 {
00356     double lb1_sqr = b1 % b1;
00357     double inv_lb1 = 1.0 / sqrt( lb1_sqr );
00358     MsqMatrix< 3, 2 > V;
00359     V.set_column( 0, b1 * inv_lb1 );
00360     V.set_column( 1, ( lb1_sqr * b2 - ( b1 % b2 ) * b1 ) * ( inv_lb1 / length( b1 * b2 ) ) );
00361     return V;
00362 }
00363 
00364 /** If, for the specified element type, the skew is constant for
00365  *  an ideal element and the aspect is identity everywhere within
00366  *  the element, pass back the constant skew/shape term and return
00367  *  true.  Otherwise return false.
00368  */
00369 static inline bool ideal_constant_skew_I_3D( EntityTopology element_type, MsqMatrix< 3, 3 >& q )
00370 {
00371     switch( element_type )
00372     {
00373         case HEXAHEDRON:
00374             q = MsqMatrix< 3, 3 >( 1.0 );  // Identity
00375             return false;
00376         case TETRAHEDRON:
00377             // [ x,   x/2, x/2 ] x^6 = 2
00378             // [ 0,   y,   y/3 ] y^2 = 3/4 x^2
00379             // [ 0,   0,   z   ] z^2 = 2/3 x^2
00380             q( 0, 0 ) = 1.122462048309373;
00381             q( 0, 1 ) = q( 0, 2 ) = 0.56123102415468651;
00382             q( 1, 0 ) = q( 2, 0 ) = q( 2, 1 ) = 0.0;
00383             q( 1, 1 )                         = 0.97208064861983279;
00384             q( 1, 2 )                         = 0.32402688287327758;
00385             q( 2, 2 )                         = 0.91648642466573493;
00386             return true;
00387         case PRISM:
00388             //            [ 1   0   0 ]
00389             //  a^(-1/3)  [ 0   1  1/2]
00390             //            [ 0   0   a ]
00391             //
00392             // a = sqrt(3)/2
00393             //
00394             q( 0, 0 ) = q( 1, 1 ) = 1.0491150634216482;
00395             q( 0, 1 ) = q( 0, 2 ) = q( 1, 0 ) = 0.0;
00396             q( 1, 2 )                         = 0.52455753171082409;
00397             q( 2, 0 ) = q( 2, 1 ) = 0.0;
00398             q( 2, 2 )             = 0.90856029641606972;
00399             return true;
00400         default:
00401             return false;
00402     }
00403 }
00404 
00405 /** If, for the specified element type, the skew is constant for
00406  *  an ideal element and the aspect is identity everywhere within
00407  *  the element, pass back the constant skew/shape term and return
00408  *  true.  Otherwise return false.
00409  */
00410 static inline bool ideal_constant_skew_I_2D( EntityTopology element_type, MsqMatrix< 2, 2 >& q )
00411 {
00412     switch( element_type )
00413     {
00414         case QUADRILATERAL:
00415             q = MsqMatrix< 2, 2 >( 1.0 );  // Identity
00416             return true;
00417         case TRIANGLE:
00418             // [ x, x/2 ]  x = pow(4/3, 0.25)
00419             // [ 0, y   ]  y = 1/x
00420             q( 0, 0 ) = 1.074569931823542;
00421             q( 0, 1 ) = 0.537284965911771;
00422             q( 1, 0 ) = 0.0;
00423             q( 1, 1 ) = 0.93060485910209956;
00424             return true;
00425         default:
00426             return false;
00427     }
00428 }
00429 
00430 void TargetCalculator::ideal_skew_3D( EntityTopology element_type,
00431                                       Sample s,
00432                                       const PatchData& pd,
00433                                       MsqMatrix< 3, 3 >& q,
00434                                       MsqError& err )
00435 {
00436     if( !ideal_constant_skew_I_3D( element_type, q ) )
00437     {
00438         const MappingFunction3D* map = pd.get_mapping_function_3D( element_type );
00439         if( !map )
00440         {
00441             MSQ_SETERR( err )( MsqError::UNSUPPORTED_ELEMENT );
00442             return;
00443         }
00444         map->ideal( s, q, err );MSQ_ERRRTN( err );
00445         q = TargetCalculator::skew( q );
00446     }
00447 }
00448 
00449 void TargetCalculator::ideal_skew_2D( EntityTopology element_type,
00450                                       Sample s,
00451                                       const PatchData& pd,
00452                                       MsqMatrix< 2, 2 >& q,
00453                                       MsqError& err )
00454 {
00455     if( !ideal_constant_skew_I_2D( element_type, q ) )
00456     {
00457         const MappingFunction2D* map = pd.get_mapping_function_2D( element_type );
00458         if( !map )
00459         {
00460             MSQ_SETERR( err )( MsqError::UNSUPPORTED_ELEMENT );
00461             return;
00462         }
00463         MsqMatrix< 3, 2 > J;
00464         map->ideal( s, J, err );MSQ_ERRRTN( err );
00465         q = TargetCalculator::skew( J );
00466     }
00467 }
00468 
00469 void TargetCalculator::ideal_shape_3D( EntityTopology element_type,
00470                                        Sample s,
00471                                        const PatchData& pd,
00472                                        MsqMatrix< 3, 3 >& q,
00473                                        MsqError& err )
00474 {
00475     if( !ideal_constant_skew_I_3D( element_type, q ) )
00476     {
00477         const MappingFunction3D* map = pd.get_mapping_function_3D( element_type );
00478         if( !map )
00479         {
00480             MSQ_SETERR( err )( MsqError::UNSUPPORTED_ELEMENT );
00481             return;
00482         }
00483         map->ideal( s, q, err );MSQ_ERRRTN( err );
00484         q = TargetCalculator::shape( q );
00485     }
00486 }
00487 
00488 void TargetCalculator::ideal_shape_2D( EntityTopology element_type,
00489                                        Sample s,
00490                                        const PatchData& pd,
00491                                        MsqMatrix< 2, 2 >& q,
00492                                        MsqError& err )
00493 {
00494     if( !ideal_constant_skew_I_2D( element_type, q ) )
00495     {
00496         const MappingFunction2D* map = pd.get_mapping_function_2D( element_type );
00497         if( !map )
00498         {
00499             MSQ_SETERR( err )( MsqError::UNSUPPORTED_ELEMENT );
00500             return;
00501         }
00502         MsqMatrix< 3, 2 > J;
00503         map->ideal( s, J, err );MSQ_ERRRTN( err );
00504         q = TargetCalculator::shape( J );
00505     }
00506 }
00507 
00508 MsqMatrix< 3, 3 > TargetCalculator::new_aspect_3D( const MsqVector< 3 >& r )
00509 {
00510     MsqMatrix< 3, 3 > W( 0.0 );
00511     W( 0, 0 ) = MBMesquite::cbrt( r[0] * r[0] / ( r[1] * r[2] ) );
00512     W( 1, 1 ) = MBMesquite::cbrt( r[1] * r[1] / ( r[0] * r[2] ) );
00513     W( 2, 2 ) = MBMesquite::cbrt( r[2] * r[2] / ( r[1] * r[0] ) );
00514     return W;
00515 }
00516 
00517 MsqMatrix< 2, 2 > TargetCalculator::new_aspect_2D( const MsqVector< 2 >& r )
00518 {
00519     return new_aspect_2D( r[0] / r[1] );
00520 }
00521 
00522 MsqMatrix< 2, 2 > TargetCalculator::new_aspect_2D( double rho )
00523 {
00524     MsqMatrix< 2, 2 > W( sqrt( rho ) );
00525     W( 1, 1 ) = 1.0 / W( 0, 0 );
00526     return W;
00527 }
00528 
00529 static NodeSet get_nodeset( EntityTopology type, int num_nodes, MsqError& err )
00530 {
00531     bool midedge, midface, midvol;
00532     TopologyInfo::higher_order( type, num_nodes, midedge, midface, midvol, err );
00533     if( MSQ_CHKERR( err ) ) return NodeSet();
00534 
00535     NodeSet bits;
00536     bits.set_all_corner_nodes( type );
00537     if( midedge ) bits.set_all_mid_edge_nodes( type );
00538     if( midface ) bits.set_all_mid_face_nodes( type );
00539     if( TopologyInfo::dimension( type ) == 3 && midvol ) bits.set_mid_region_node();
00540 
00541     return bits;
00542 }
00543 
00544 void TargetCalculator::jacobian_3D( PatchData& pd,
00545                                     EntityTopology type,
00546                                     int num_nodes,
00547                                     Sample location,
00548                                     const Vector3D* coords,
00549                                     MsqMatrix< 3, 3 >& J,
00550                                     MsqError& err )
00551 {
00552     // Get element properties
00553     NodeSet bits = get_nodeset( type, num_nodes, err );MSQ_ERRRTN( err );
00554     const MappingFunction3D* mf = pd.get_mapping_function_3D( type );
00555     if( !mf )
00556     {
00557         MSQ_SETERR( err )( MsqError::UNSUPPORTED_ELEMENT );
00558         return;
00559     }
00560 
00561     // Get mapping function derivatives
00562     const int MAX_NODES = 27;
00563     assert( num_nodes <= MAX_NODES );
00564     size_t indices[MAX_NODES], n;
00565     MsqVector< 3 > derivs[MAX_NODES];
00566     mf->derivatives( location, bits, indices, derivs, n, err );MSQ_ERRRTN( err );
00567 
00568     // calculate Jacobian
00569     assert( sizeof( Vector3D ) == sizeof( MsqVector< 3 > ) );
00570     const MsqVector< 3 >* verts = reinterpret_cast< const MsqVector< 3 >* >( coords );
00571     assert( n > 0 );
00572     J = outer( verts[indices[0]], derivs[0] );
00573     for( size_t i = 1; i < n; ++i )
00574         J += outer( verts[indices[i]], derivs[i] );
00575 }
00576 
00577 void TargetCalculator::jacobian_2D( PatchData& pd,
00578                                     EntityTopology type,
00579                                     int num_nodes,
00580                                     Sample location,
00581                                     const Vector3D* coords,
00582                                     MsqMatrix< 3, 2 >& J,
00583                                     MsqError& err )
00584 {
00585     // Get element properties
00586     NodeSet bits = get_nodeset( type, num_nodes, err );MSQ_ERRRTN( err );
00587     const MappingFunction2D* mf = pd.get_mapping_function_2D( type );
00588     if( !mf )
00589     {
00590         MSQ_SETERR( err )( MsqError::UNSUPPORTED_ELEMENT );
00591         return;
00592     }
00593 
00594     // Get mapping function derivatives
00595     const int MAX_NODES = 9;
00596     assert( num_nodes <= MAX_NODES );
00597     size_t indices[MAX_NODES], n;
00598     MsqVector< 2 > derivs[MAX_NODES];
00599     mf->derivatives( location, bits, indices, derivs, n, err );MSQ_ERRRTN( err );
00600 
00601     // calculate Jacobian
00602     assert( sizeof( Vector3D ) == sizeof( MsqVector< 3 > ) );
00603     const MsqVector< 3 >* verts = reinterpret_cast< const MsqVector< 3 >* >( coords );
00604     assert( n > 0 );
00605     J = outer( verts[indices[0]], derivs[0] );
00606     for( size_t i = 1; i < n; ++i )
00607         J += outer( verts[indices[i]], derivs[i] );
00608 }
00609 
00610 void TargetCalculator::get_refmesh_Jacobian_3D( ReferenceMeshInterface* ref_mesh,
00611                                                 PatchData& pd,
00612                                                 size_t element,
00613                                                 Sample sample,
00614                                                 MsqMatrix< 3, 3 >& W_out,
00615                                                 MsqError& err )
00616 {
00617     // get element
00618     MsqMeshEntity& elem       = pd.element_by_index( element );
00619     const EntityTopology type = elem.get_element_type();
00620     const unsigned n          = elem.node_count();
00621 
00622     const unsigned MAX_NODES = 27;
00623     assert( n <= MAX_NODES );
00624 
00625     // get vertices
00626     Mesh::VertexHandle elem_verts[MAX_NODES];
00627     const std::size_t* vtx_idx        = elem.get_vertex_index_array();
00628     const Mesh::VertexHandle* vtx_hdl = pd.get_vertex_handles_array();
00629     for( unsigned i = 0; i < n; ++i )
00630         elem_verts[i] = vtx_hdl[vtx_idx[i]];
00631 
00632     // get vertex coordinates
00633     Vector3D vert_coords[MAX_NODES];
00634     ref_mesh->get_reference_vertex_coordinates( elem_verts, n, vert_coords, err );MSQ_ERRRTN( err );
00635 
00636     // calculate Jacobian
00637     jacobian_3D( pd, type, n, sample, vert_coords, W_out, err );MSQ_ERRRTN( err );
00638 }
00639 
00640 void TargetCalculator::get_refmesh_Jacobian_2D( ReferenceMeshInterface* ref_mesh,
00641                                                 PatchData& pd,
00642                                                 size_t element,
00643                                                 Sample sample,
00644                                                 MsqMatrix< 3, 2 >& W_out,
00645                                                 MsqError& err )
00646 {
00647     // get element
00648     MsqMeshEntity& elem       = pd.element_by_index( element );
00649     const EntityTopology type = elem.get_element_type();
00650     const unsigned n          = elem.node_count();
00651 
00652     const unsigned MAX_NODES = 9;
00653     assert( n <= MAX_NODES );
00654 
00655     // get vertices
00656     Mesh::VertexHandle elem_verts[MAX_NODES];
00657     const std::size_t* vtx_idx        = elem.get_vertex_index_array();
00658     const Mesh::VertexHandle* vtx_hdl = pd.get_vertex_handles_array();
00659     for( unsigned i = 0; i < n; ++i )
00660         elem_verts[i] = vtx_hdl[vtx_idx[i]];
00661 
00662     // get vertex coordinates
00663     Vector3D vert_coords[MAX_NODES];
00664     ref_mesh->get_reference_vertex_coordinates( elem_verts, n, vert_coords, err );MSQ_ERRRTN( err );
00665 
00666     // calculate Jacobian
00667     jacobian_2D( pd, type, n, sample, vert_coords, W_out, err );MSQ_ERRRTN( err );
00668 }
00669 
00670 TargetCalculator::~TargetCalculator() {}
00671 
00672 void TargetCalculator::initialize_queue( MeshDomainAssoc*, const Settings*, MsqError& ) {}
00673 
00674 }  // namespace MBMesquite
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines