MOAB: Mesh Oriented datABase  (version 5.4.1)
MeshUtil.cpp
Go to the documentation of this file.
00001 /* *****************************************************************
00002     MESQUITE -- The Mesh Quality Improvement Toolkit
00003 
00004     Copyright 2010 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     (2010) [email protected]
00024 
00025   ***************************************************************** */
00026 
00027 /** \file MeshUtil.cpp
00028  *  \brief Implement MeshUtil class
00029  *  \author Jason Kraftcheck
00030  */
00031 
00032 #include "Mesquite.hpp"
00033 #include "MeshUtil.hpp"
00034 #include "MeshInterface.hpp"
00035 #include "EdgeIterator.hpp"
00036 #include "SimpleStats.hpp"
00037 #include "MsqError.hpp"
00038 #include "TMPQualityMetric.hpp"
00039 #include "MappingFunction.hpp"
00040 #include "TargetCalculator.hpp"
00041 
00042 #include <vector>
00043 #include <algorithm>
00044 #include <limits>
00045 #include <sstream>
00046 
00047 namespace MBMesquite
00048 {
00049 
00050 MeshUtil::~MeshUtil()
00051 {
00052     delete globalPatch;
00053 }
00054 
00055 PatchData* MeshUtil::get_global_patch( MsqError& err )
00056 {
00057     if( !globalPatch )
00058     {
00059         globalPatch = new PatchData;
00060         globalPatch->set_mesh( myMesh );
00061         if( mySettings ) globalPatch->attach_settings( mySettings );
00062         globalPatch->fill_global_patch( err );
00063         if( MSQ_CHKERR( err ) )
00064         {
00065             delete globalPatch;
00066             globalPatch = 0;
00067         }
00068     }
00069     return globalPatch;
00070 }
00071 
00072 void MeshUtil::edge_length_distribution( SimpleStats& results, MsqError& err )
00073 {
00074     PatchData* pd = get_global_patch( err );MSQ_ERRRTN( err );
00075     EdgeIterator iter( pd, err );MSQ_ERRRTN( err );
00076     while( !iter.is_at_end() )
00077     {
00078         Vector3D diff = iter.start() - iter.end();
00079         results.add_squared( diff % diff );
00080         iter.step( err );MSQ_ERRRTN( err );
00081     }
00082 
00083     if( results.empty() )
00084     {  // no mesh
00085         MSQ_SETERR( err )( "Mesh contains no elements", MsqError::INVALID_MESH );
00086     }
00087 }
00088 
00089 void MeshUtil::lambda_distribution( SimpleStats& results, MsqError& err )
00090 {
00091     PatchData& pd = *get_global_patch( err );MSQ_ERRRTN( err );
00092 
00093     std::vector< size_t > handles;
00094     TMPQualityMetric::get_patch_evaluations( pd, handles, false, err );
00095 
00096     const size_t N = 64;
00097     size_t count;
00098     size_t indices[N];
00099     MsqVector< 2 > derivs2D[N];
00100     MsqVector< 3 > derivs3D[N];
00101 
00102     for( size_t i = 0; i < handles.size(); ++i )
00103     {
00104         double lambda;
00105         const Sample s      = ElemSampleQM::sample( handles[i] );
00106         const size_t e      = ElemSampleQM::elem( handles[i] );
00107         EntityTopology type = pd.element_by_index( e ).get_element_type();
00108         const NodeSet bits  = pd.non_slave_node_set( e );
00109 
00110         if( TopologyInfo::dimension( type ) == 3 )
00111         {
00112             const MappingFunction3D* mf = pd.get_mapping_function_3D( type );
00113             if( !mf )
00114             {
00115                 MSQ_SETERR( err )
00116                 ( "No mapping function for element type", MsqError::UNSUPPORTED_ELEMENT );
00117                 return;
00118             }
00119 
00120             MsqMatrix< 3, 3 > W;
00121             mf->jacobian( pd, e, bits, s, indices, derivs3D, count, W, err );MSQ_ERRRTN( err );
00122             assert( N >= count );
00123             lambda = TargetCalculator::size( W );
00124         }
00125         else
00126         {
00127             assert( TopologyInfo::dimension( type ) == 2 );
00128             const MappingFunction2D* mf = pd.get_mapping_function_2D( type );
00129             if( !mf )
00130             {
00131                 MSQ_SETERR( err )
00132                 ( "No mapping function for element type", MsqError::UNSUPPORTED_ELEMENT );
00133                 return;
00134             }
00135 
00136             MsqMatrix< 3, 2 > W;
00137             mf->jacobian( pd, e, bits, s, indices, derivs2D, count, W, err );MSQ_ERRRTN( err );
00138             assert( N >= count );
00139             lambda = TargetCalculator::size( W );
00140         }
00141 
00142         results.add_value( lambda );
00143     }
00144 
00145     if( results.empty() )
00146     {  // no mesh
00147         MSQ_SETERR( err )( "Mesh contains no elements", MsqError::INVALID_MESH );
00148     }
00149 }
00150 
00151 bool MeshUtil::meshes_are_different( Mesh& mesh1, Mesh& mesh2, MsqError& err, double tol, bool do_print )
00152 {
00153     // MsqError& err = (do_print ? MsqPrintError(cout) : MsqError());
00154 
00155     // mesh 1
00156     std::vector< Mesh::ElementHandle > elements1;
00157     std::vector< Mesh::VertexHandle > vertices1;
00158     std::vector< Mesh::VertexHandle > connectivity1;
00159     std::vector< size_t > offsets1;
00160 
00161     mesh1.get_all_elements( elements1, err );
00162     MSQ_ERRZERO( err );
00163     mesh1.get_all_vertices( vertices1, err );
00164     MSQ_ERRZERO( err );
00165     mesh1.elements_get_attached_vertices( arrptr( elements1 ), elements1.size(), connectivity1, offsets1, err );
00166     MSQ_ERRZERO( err );
00167     std::vector< EntityTopology > types1( elements1.size() );
00168     mesh1.elements_get_topologies( arrptr( elements1 ), arrptr( types1 ), elements1.size(), err );
00169     MSQ_ERRZERO( err );
00170 
00171     // mesh 2
00172     std::vector< Mesh::ElementHandle > elements2;
00173     std::vector< Mesh::VertexHandle > vertices2;
00174     std::vector< Mesh::VertexHandle > connectivity2;
00175     std::vector< size_t > offsets2;
00176 
00177     mesh2.get_all_elements( elements2, err );
00178     MSQ_ERRZERO( err );
00179     mesh2.get_all_vertices( vertices2, err );
00180     MSQ_ERRZERO( err );
00181     mesh2.elements_get_attached_vertices( arrptr( elements2 ), elements2.size(), connectivity2, offsets2, err );
00182     MSQ_ERRZERO( err );
00183     std::vector< EntityTopology > types2( elements2.size() );
00184     mesh2.elements_get_topologies( arrptr( elements2 ), arrptr( types2 ), elements2.size(), err );
00185     MSQ_ERRZERO( err );
00186 
00187 #define MDRET( msg )                                                                                      \
00188     do                                                                                                    \
00189     {                                                                                                     \
00190         if( do_print ) std::cout << "MeshUtil::mesh_diff meshes are different: " << ( msg ) << std::endl; \
00191         return true;                                                                                      \
00192     } while( 0 )
00193     if( elements1.size() != elements2.size() ) MDRET( "elements sizes differ" );
00194     if( vertices1.size() != vertices2.size() ) MDRET( "vertices sizes differ" );
00195     if( offsets1.size() != offsets2.size() ) MDRET( "offsets sizes differ" );
00196     if( connectivity1.size() != connectivity2.size() ) MDRET( "connectivity sizes differ" );
00197 
00198     for( unsigned i = 0; i < offsets1.size(); i++ )
00199     {
00200         if( offsets1[i] != offsets2[i] ) MDRET( "offets differ" );
00201     }
00202 
00203     for( unsigned i = 0; i < vertices1.size(); i++ )
00204     {
00205         MsqVertex vert1, vert2;
00206         mesh1.vertices_get_coordinates( &vertices1[i], &vert1, 1, err );
00207         MSQ_ERRZERO( err );
00208         mesh2.vertices_get_coordinates( &vertices2[i], &vert2, 1, err );
00209         MSQ_ERRZERO( err );
00210         double dist = Vector3D::distance_between( vert1, vert2 );
00211         double v1   = vert1.length();
00212         double v2   = vert2.length();
00213         if( dist > 0.5 * ( v1 + v2 ) * tol )
00214         {
00215             std::ostringstream ost;
00216             ost << "vertices coordinates differ more than tolerance [" << tol << "], vert1= " << vert1
00217                 << " vert2= " << vert2;
00218             MDRET( ost.str() );
00219         }
00220     }
00221 
00222     for( unsigned i = 0; i < connectivity1.size(); i++ )
00223     {
00224         MsqVertex vert1, vert2;
00225         mesh1.vertices_get_coordinates( &connectivity1[i], &vert1, 1, err );
00226         MSQ_ERRZERO( err );
00227         mesh2.vertices_get_coordinates( &connectivity2[i], &vert2, 1, err );
00228         MSQ_ERRZERO( err );
00229         double dist = Vector3D::distance_between( vert1, vert2 );
00230         double v1   = vert1.length();
00231         double v2   = vert2.length();
00232         if( dist > 0.5 * ( v1 + v2 ) * tol )
00233         {
00234             std::ostringstream ost;
00235             ost << "connectivity coordinates differ more than tolerance [" << tol << "], vert1= " << vert1
00236                 << " vert2= " << vert2;
00237             MDRET( ost.str() );
00238         }
00239     }
00240 
00241     return false;
00242 }
00243 
00244 }  // namespace MBMesquite
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines