MOAB: Mesh Oriented datABase  (version 5.2.1)
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, 0.0, sqrt( 3.0 / 2.0 )
00056 };
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, size_t p_elem, std::vector< size_t >& handles,
00085                                                MsqError& err )
00086 {
00087     get_elem_sample_points( pd, p_elem, handles, err );
00088 }
00089 
00090 bool AffineMapMetric::evaluate( PatchData& pd, size_t p_handle, double& value, MsqError& err )
00091 {
00092     Sample s              = ElemSampleQM::sample( p_handle );
00093     size_t e              = ElemSampleQM::elem( p_handle );
00094     MsqMeshEntity& p_elem = pd.element_by_index( e );
00095     EntityTopology type   = p_elem.get_element_type();
00096     unsigned edim         = TopologyInfo::dimension( type );
00097     const size_t* conn    = p_elem.get_vertex_index_array();
00098 
00099     // This metric only supports sampling at corners, except for simplices.
00100     // If element is a simpex, then the Jacobian is constant over a linear
00101     // element.  In this case, always evaluate at any vertex.
00102     // unsigned corner = s.number;
00103     if( s.dimension != 0 )
00104     {
00105         if( type == TRIANGLE || type == TETRAHEDRON )
00106             /*corner = 0*/;
00107         else
00108         {
00109             MSQ_SETERR( err )
00110             ( "Invalid sample point for AffineMapMetric", MsqError::UNSUPPORTED_ELEMENT );
00111             return false;
00112         }
00113     }
00114 
00115     bool rval;
00116     if( edim == 3 )
00117     {  // 3x3 or 3x2 targets ?
00118         Vector3D c[3] = { Vector3D( 0, 0, 0 ), Vector3D( 0, 0, 0 ), Vector3D( 0, 0, 0 ) };
00119         unsigned n;
00120         const unsigned* adj = TopologyInfo::adjacent_vertices( type, s.number, n );
00121         c[0]                = pd.vertex_by_index( conn[adj[0]] ) - pd.vertex_by_index( conn[s.number] );
00122         c[1]                = pd.vertex_by_index( conn[adj[1]] ) - pd.vertex_by_index( conn[s.number] );
00123         c[2]                = pd.vertex_by_index( conn[adj[2]] ) - pd.vertex_by_index( conn[s.number] );
00124         MsqMatrix< 3, 3 > A;
00125         A.set_column( 0, MsqMatrix< 3, 1 >( c[0].to_array() ) );
00126         A.set_column( 1, MsqMatrix< 3, 1 >( c[1].to_array() ) );
00127         A.set_column( 2, MsqMatrix< 3, 1 >( c[2].to_array() ) );
00128         if( type == TETRAHEDRON ) A = A * TET_XFORM;
00129 
00130         MsqMatrix< 3, 3 > W;
00131         targetCalc->get_3D_target( pd, e, s, W, err );
00132         MSQ_ERRZERO( err );
00133         rval = targetMetric->evaluate( A * inverse( W ), value, err );
00134         MSQ_ERRZERO( err );
00135     }
00136     else
00137     {
00138         Vector3D c[2] = { Vector3D( 0, 0, 0 ), Vector3D( 0, 0, 0 ) };
00139         unsigned n;
00140         const unsigned* adj = TopologyInfo::adjacent_vertices( type, s.number, n );
00141         c[0]                = pd.vertex_by_index( conn[adj[0]] ) - pd.vertex_by_index( conn[s.number] );
00142         c[1]                = pd.vertex_by_index( conn[adj[1]] ) - pd.vertex_by_index( conn[s.number] );
00143         MsqMatrix< 3, 2 > App;
00144         App.set_column( 0, MsqMatrix< 3, 1 >( c[0].to_array() ) );
00145         App.set_column( 1, MsqMatrix< 3, 1 >( c[1].to_array() ) );
00146 
00147         MsqMatrix< 3, 2 > Wp;
00148         targetCalc->get_surface_target( pd, e, s, Wp, err );
00149         MSQ_ERRZERO( err );
00150 
00151         MsqMatrix< 2, 2 > A, W;
00152         MsqMatrix< 3, 2 > RZ;
00153         surface_to_2d( App, Wp, W, RZ );
00154         A = transpose( RZ ) * App;
00155         if( type == TRIANGLE ) A = A * TRI_XFORM;
00156 
00157         rval = targetMetric->evaluate( A * inverse( W ), value, err );
00158         MSQ_ERRZERO( err );
00159     }
00160 
00161     // apply target weight to value
00162     if( weightCalc )
00163     {
00164         double ck = weightCalc->get_weight( pd, e, s, err );
00165         MSQ_ERRZERO( err );
00166         value *= ck;
00167     }
00168     return rval;
00169 }
00170 
00171 bool AffineMapMetric::evaluate_with_indices( PatchData& pd, size_t p_handle, double& value,
00172                                              std::vector< size_t >& indices, MsqError& err )
00173 {
00174     Sample s              = ElemSampleQM::sample( p_handle );
00175     size_t e              = ElemSampleQM::elem( p_handle );
00176     MsqMeshEntity& p_elem = pd.element_by_index( e );
00177     EntityTopology type   = p_elem.get_element_type();
00178     const size_t* conn    = p_elem.get_vertex_index_array();
00179 
00180     // this metric only supports sampling at corners
00181     if( s.dimension != 0 )
00182     {
00183         if( type != TRIANGLE && type != TETRAHEDRON )
00184         {
00185             MSQ_SETERR( err )
00186             ( "Invalid sample point for AffineMapMetric", MsqError::UNSUPPORTED_ELEMENT );
00187             return false;
00188         }
00189         s.dimension = 0;
00190         s.number    = 0;
00191     }
00192 
00193     unsigned n;
00194     const unsigned* adj = TopologyInfo::adjacent_vertices( type, s.number, n );
00195     indices.clear();
00196     if( conn[s.number] < pd.num_free_vertices() ) indices.push_back( conn[s.number] );
00197     for( unsigned i = 0; i < n; ++i )
00198         if( conn[adj[i]] < pd.num_free_vertices() ) indices.push_back( conn[adj[i]] );
00199 
00200     return evaluate( pd, p_handle, value, err );
00201 }
00202 
00203 }  // namespace MBMesquite
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines