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