MOAB: Mesh Oriented datABase  (version 5.4.0)
AffineMapMetric.cpp
Go to the documentation of this file.
00001 /* *****************************************************************
00002     MESQUITE -- The Mesh Quality Improvement Toolkit
00003 
00004     Copyright 2007 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     (2007) kraftche@cae.wisc.edu
00024 
00025   ***************************************************************** */
00026 
00027 /** \file AffineMapMetric.cpp
00028  *  \brief
00029  *  \author Jason Kraftcheck
00030  */
00031 
00032 #include "Mesquite.hpp"
00033 #include "AffineMapMetric.hpp"
00034 #include "MsqMatrix.hpp"
00035 #include "ElementQM.hpp"
00036 #include "MsqError.hpp"
00037 #include "Vector3D.hpp"
00038 #include "PatchData.hpp"
00039 #include "MappingFunction.hpp"
00040 #include "WeightCalculator.hpp"
00041 #include "TargetCalculator.hpp"
00042 #include "TMetric.hpp"
00043 #include "TargetMetricUtil.hpp"
00044 
00045 #include <functional>
00046 #include <algorithm>
00047 
00048 namespace MBMesquite
00049 {
00050 
00051 const double TRI_XFORM_VALS[] = { 1.0, -1.0 / sqrt( 3.0 ), 0.0, 2.0 / sqrt( 3.0 ) };
00052 MsqMatrix< 2, 2 > TRI_XFORM( TRI_XFORM_VALS );
00053 
00054 const double TET_XFORM_VALS[] = {
00055     1.0, -1.0 / sqrt( 3.0 ), -1.0 / sqrt( 6.0 ), 0.0, 2.0 / sqrt( 3.0 ), -1.0 / sqrt( 6.0 ), 0.0,
00056     0.0, sqrt( 3.0 / 2.0 ) };
00057 MsqMatrix< 3, 3 > TET_XFORM( TET_XFORM_VALS );
00058 
00059 AffineMapMetric::AffineMapMetric( TargetCalculator* tc, WeightCalculator* wc, TMetric* target_metric )
00060     : targetCalc( tc ), weightCalc( wc ), targetMetric( target_metric )
00061 {
00062 }
00063 
00064 AffineMapMetric::AffineMapMetric( TargetCalculator* tc, TMetric* target_metric )
00065     : targetCalc( tc ), weightCalc( 0 ), targetMetric( target_metric )
00066 {
00067 }
00068 
00069 int AffineMapMetric::get_negate_flag() const
00070 {
00071     return 1;
00072 }
00073 
00074 std::string AffineMapMetric::get_name() const
00075 {
00076     return std::string( "AffineMap(" ) + targetMetric->get_name() + ')';
00077 }
00078 
00079 void AffineMapMetric::get_evaluations( PatchData& pd, std::vector< size_t >& handles, bool free, MsqError& err )
00080 {
00081     get_sample_pt_evaluations( pd, handles, free, err );
00082 }
00083 
00084 void AffineMapMetric::get_element_evaluations( PatchData& pd,
00085                                                size_t p_elem,
00086                                                std::vector< size_t >& handles,
00087                                                MsqError& err )
00088 {
00089     get_elem_sample_points( pd, p_elem, handles, err );
00090 }
00091 
00092 bool AffineMapMetric::evaluate( PatchData& pd, size_t p_handle, double& value, MsqError& err )
00093 {
00094     Sample s              = ElemSampleQM::sample( p_handle );
00095     size_t e              = ElemSampleQM::elem( p_handle );
00096     MsqMeshEntity& p_elem = pd.element_by_index( e );
00097     EntityTopology type   = p_elem.get_element_type();
00098     unsigned edim         = TopologyInfo::dimension( type );
00099     const size_t* conn    = p_elem.get_vertex_index_array();
00100 
00101     // This metric only supports sampling at corners, except for simplices.
00102     // If element is a simpex, then the Jacobian is constant over a linear
00103     // element.  In this case, always evaluate at any vertex.
00104     // unsigned corner = s.number;
00105     if( s.dimension != 0 )
00106     {
00107         if( type == TRIANGLE || type == TETRAHEDRON )
00108             /*corner = 0*/;
00109         else
00110         {
00111             MSQ_SETERR( err )
00112             ( "Invalid sample point for AffineMapMetric", MsqError::UNSUPPORTED_ELEMENT );
00113             return false;
00114         }
00115     }
00116 
00117     bool rval;
00118     if( edim == 3 )
00119     {  // 3x3 or 3x2 targets ?
00120         Vector3D c[3] = { Vector3D( 0, 0, 0 ), Vector3D( 0, 0, 0 ), Vector3D( 0, 0, 0 ) };
00121         unsigned n;
00122         const unsigned* adj = TopologyInfo::adjacent_vertices( type, s.number, n );
00123         c[0]                = pd.vertex_by_index( conn[adj[0]] ) - pd.vertex_by_index( conn[s.number] );
00124         c[1]                = pd.vertex_by_index( conn[adj[1]] ) - pd.vertex_by_index( conn[s.number] );
00125         c[2]                = pd.vertex_by_index( conn[adj[2]] ) - pd.vertex_by_index( conn[s.number] );
00126         MsqMatrix< 3, 3 > A;
00127         A.set_column( 0, MsqMatrix< 3, 1 >( c[0].to_array() ) );
00128         A.set_column( 1, MsqMatrix< 3, 1 >( c[1].to_array() ) );
00129         A.set_column( 2, MsqMatrix< 3, 1 >( c[2].to_array() ) );
00130         if( type == TETRAHEDRON ) A = A * TET_XFORM;
00131 
00132         MsqMatrix< 3, 3 > W;
00133         targetCalc->get_3D_target( pd, e, s, W, err );
00134         MSQ_ERRZERO( err );
00135         rval = targetMetric->evaluate( A * inverse( W ), value, err );
00136         MSQ_ERRZERO( err );
00137     }
00138     else
00139     {
00140         Vector3D c[2] = { Vector3D( 0, 0, 0 ), Vector3D( 0, 0, 0 ) };
00141         unsigned n;
00142         const unsigned* adj = TopologyInfo::adjacent_vertices( type, s.number, n );
00143         c[0]                = pd.vertex_by_index( conn[adj[0]] ) - pd.vertex_by_index( conn[s.number] );
00144         c[1]                = pd.vertex_by_index( conn[adj[1]] ) - pd.vertex_by_index( conn[s.number] );
00145         MsqMatrix< 3, 2 > App;
00146         App.set_column( 0, MsqMatrix< 3, 1 >( c[0].to_array() ) );
00147         App.set_column( 1, MsqMatrix< 3, 1 >( c[1].to_array() ) );
00148 
00149         MsqMatrix< 3, 2 > Wp;
00150         targetCalc->get_surface_target( pd, e, s, Wp, err );
00151         MSQ_ERRZERO( err );
00152 
00153         MsqMatrix< 2, 2 > A, W;
00154         MsqMatrix< 3, 2 > RZ;
00155         surface_to_2d( App, Wp, W, RZ );
00156         A = transpose( RZ ) * App;
00157         if( type == TRIANGLE ) A = A * TRI_XFORM;
00158 
00159         rval = targetMetric->evaluate( A * inverse( W ), value, err );
00160         MSQ_ERRZERO( err );
00161     }
00162 
00163     // apply target weight to value
00164     if( weightCalc )
00165     {
00166         double ck = weightCalc->get_weight( pd, e, s, err );
00167         MSQ_ERRZERO( err );
00168         value *= ck;
00169     }
00170     return rval;
00171 }
00172 
00173 bool AffineMapMetric::evaluate_with_indices( PatchData& pd,
00174                                              size_t p_handle,
00175                                              double& value,
00176                                              std::vector< size_t >& indices,
00177                                              MsqError& err )
00178 {
00179     Sample s              = ElemSampleQM::sample( p_handle );
00180     size_t e              = ElemSampleQM::elem( p_handle );
00181     MsqMeshEntity& p_elem = pd.element_by_index( e );
00182     EntityTopology type   = p_elem.get_element_type();
00183     const size_t* conn    = p_elem.get_vertex_index_array();
00184 
00185     // this metric only supports sampling at corners
00186     if( s.dimension != 0 )
00187     {
00188         if( type != TRIANGLE && type != TETRAHEDRON )
00189         {
00190             MSQ_SETERR( err )
00191             ( "Invalid sample point for AffineMapMetric", MsqError::UNSUPPORTED_ELEMENT );
00192             return false;
00193         }
00194         s.dimension = 0;
00195         s.number    = 0;
00196     }
00197 
00198     unsigned n;
00199     const unsigned* adj = TopologyInfo::adjacent_vertices( type, s.number, n );
00200     indices.clear();
00201     if( conn[s.number] < pd.num_free_vertices() ) indices.push_back( conn[s.number] );
00202     for( unsigned i = 0; i < n; ++i )
00203         if( conn[adj[i]] < pd.num_free_vertices() ) indices.push_back( conn[adj[i]] );
00204 
00205     return evaluate( pd, p_handle, value, err );
00206 }
00207 
00208 }  // namespace MBMesquite
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines