MOAB: Mesh Oriented datABase  (version 5.3.1)
main.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     (2010) kraftche@cae.wisc.edu
00024 
00025   ***************************************************************** */
00026 
00027 #include "Mesquite.hpp"
00028 #include "MsqIGeom.hpp"
00029 #include "FacetModifyEngine.hpp"
00030 #include "MeshImpl.hpp"
00031 #include "MsqError.hpp"
00032 #include "ShapeImprovementWrapper.hpp"
00033 #include "MsqVertex.hpp"
00034 #include "QualityAssessor.hpp"
00035 #include "SphericalDomain.hpp"
00036 #include "TestUtil.hpp"
00037 
00038 #include <memory>
00039 #include <iostream>
00040 
00041 #include "MsqIBase.hpp"
00042 
00043 using namespace MBMesquite;
00044 
00045 // characteristics of geometry
00046 const double SPHERE_RADIUS = 3.0;
00047 const Vector3D SPHERE_CENTER( 2.0, 2.0, 0.0 );
00048 
00049 #define CHKIGEOM \
00050     if( chk_igeom_error( ierr, __FILE__, __LINE__ ) ) return 0
00051 
00052 bool chk_igeom_error( int ierr, const char* file, int line )
00053 {
00054     if( iBase_SUCCESS == ierr ) return false;
00055     std::cerr << "iGeom call failed at " << file << ":" << line << std::endl;
00056     std::cerr << process_itaps_error( ierr ) << std::endl;
00057     return true;
00058 }
00059 
00060 std::string default_file_name = TestDir + "unittest/mesquite/2D/vtk/quads/untangled/quads_on_sphere_529.vtk";
00061 
00062 void usage()
00063 {
00064     std::cout << "main [-N] [filename] [-o <output_file>]" << std::endl;
00065     std::cout << "  -N : Use native representation instead of iGeom implementation\n" << std::endl;
00066     std::cout << "  If no file name is specified, will use \"" << default_file_name << '"' << std::endl;
00067     exit( 1 );
00068 }
00069 
00070 void run_smoother( Mesh& mesh, MeshDomain* domain, MsqError& err );
00071 bool check_results( Mesh& mesh, MeshDomain* domain, MsqError& err );
00072 MeshDomain* get_itaps_domain();
00073 MeshDomain* get_mesquite_domain();
00074 
00075 int main( int argc, char* argv[] )
00076 {
00077     // command line arguments
00078     const char* file_name   = 0;
00079     const char* output_file = 0;
00080     bool use_native = false, opts_done = false;
00081     for( int arg = 1; arg < argc; ++arg )
00082     {
00083         if( !opts_done && argv[arg][0] == '-' )
00084         {
00085             if( !strcmp( argv[arg], "-N" ) )
00086                 use_native = true;
00087             else if( !strcmp( argv[arg], "-o" ) )
00088                 output_file = argv[arg++];
00089             else if( !strcmp( argv[arg], "--" ) )
00090                 opts_done = true;
00091             else
00092                 usage();
00093         }
00094         else if( !file_name )
00095             file_name = argv[arg];
00096         else
00097             usage();
00098     }
00099     if( !file_name )
00100     {
00101         file_name = default_file_name.c_str();
00102         std::cout << "No file specified: using default: " << default_file_name << std::endl;
00103     }
00104 
00105     MsqError err;
00106     MeshImpl mesh;
00107     mesh.read_vtk( file_name, err );
00108     if( MSQ_CHKERR( err ) )
00109     {
00110         std::cerr << err << std::endl;
00111         std::cerr << "Failed to read input file: " << file_name << std::endl;
00112         return 1;
00113     }
00114 
00115     std::unique_ptr< MeshDomain > dom( use_native ? get_mesquite_domain() : get_itaps_domain() );
00116 
00117     if( !dom.get() )
00118     {
00119         MSQ_SETERR( err )( "Domain creation failed", MsqError::INTERNAL_ERROR );
00120         std::cerr << err << std::endl;
00121         return 1;
00122     }
00123 
00124     run_smoother( mesh, dom.get(), err );
00125     if( MSQ_CHKERR( err ) )
00126     {
00127         std::cerr << err << std::endl;
00128         std::cerr << "Smoother failed." << std::endl;
00129         return 1;
00130     }
00131 
00132     if( output_file )
00133     {
00134         mesh.write_vtk( output_file, err );
00135         if( MSQ_CHKERR( err ) )
00136         {
00137             std::cerr << err << std::endl;
00138             std::cerr << "Failed to write file: " << output_file << std::endl;
00139             return 1;
00140         }
00141     }
00142 
00143     bool valid = check_results( mesh, dom.get(), err );
00144     if( MSQ_CHKERR( err ) )
00145     {
00146         std::cerr << err << std::endl;
00147         std::cerr << "Error while validating results." << std::endl;
00148         return 1;
00149     }
00150 
00151     return !valid;
00152 }
00153 
00154 void run_smoother( Mesh& mesh, MeshDomain* dom, MsqError& err )
00155 {
00156     ShapeImprovementWrapper smoother;
00157     MeshDomainAssoc mesh_and_domain = MeshDomainAssoc( &mesh, dom );
00158     smoother.run_instructions( &mesh_and_domain, err );MSQ_CHKERR( err );
00159 
00160     if( smoother.quality_assessor().invalid_elements() )
00161     {
00162         MSQ_SETERR( err )( "Smoothing produced invalid elements", MsqError::INVALID_MESH );
00163         return;
00164     }
00165 }
00166 
00167 bool check_results( Mesh& mesh, MeshDomain* dom, MsqError& err )
00168 {
00169     double EPS = 1e-2;
00170 
00171     std::vector< Mesh::VertexHandle > handles;
00172     mesh.get_all_vertices( handles, err );
00173     MSQ_ERRZERO( err );
00174     std::vector< MsqVertex > coords( handles.size() );
00175     mesh.vertices_get_coordinates( arrptr( handles ), arrptr( coords ), handles.size(), err );
00176     MSQ_ERRZERO( err );
00177 
00178     bool valid = true;
00179     size_t c   = 0;
00180     for( size_t i = 0; i < coords.size(); ++i )
00181     {
00182         double d = ( coords[i] - SPHERE_CENTER ).length();
00183         if( fabs( d - SPHERE_RADIUS ) > EPS ) ++c;
00184     }
00185     if( c )
00186     {
00187         std::cerr << c << " vertices not on domain" << std::endl;
00188         valid = false;
00189     }
00190 
00191     std::vector< unsigned short > dof( handles.size() ), exp_dof( handles.size(), 2 );
00192     dom->domain_DoF( arrptr( handles ), arrptr( dof ), handles.size(), err );
00193     MSQ_ERRZERO( err );
00194     if( dof != exp_dof )
00195     {
00196         std::cerr << "Invalid domain dimension for one or more vertices" << std::endl;
00197         valid = false;
00198     }
00199 
00200     c = 0;
00201     std::vector< Vector3D > normals( coords.begin(), coords.end() );
00202     dom->vertex_normal_at( arrptr( handles ), arrptr( normals ), handles.size(), err );
00203     MSQ_ERRZERO( err );
00204     for( size_t i = 0; i < handles.size(); ++i )
00205     {
00206         Vector3D exp_norm = ~( coords[i] - SPHERE_CENTER );
00207         Vector3D act_norm = ~normals[i];
00208         Vector3D diff     = exp_norm - act_norm;
00209         if( diff.length() > EPS ) ++c;
00210     }
00211     if( c )
00212     {
00213         std::cerr << c << " invalid vertex normals" << std::endl;
00214         valid = false;
00215     }
00216 
00217     return valid;
00218 }
00219 
00220 MeshDomain* get_itaps_domain()
00221 {
00222     const double EPS = 1e-2;
00223     int ierr;
00224     iGeom_Instance igeom;
00225     iGeom_newGeom( "", &igeom, &ierr, 0 );
00226     CHKIGEOM;
00227 
00228 #ifdef MOAB_HAVE_CGM_FACET
00229     FacetModifyEngine::set_modify_enabled( CUBIT_TRUE );
00230 #endif
00231 
00232     iBase_EntityHandle sphere_vol;
00233     iGeom_createSphere( igeom, SPHERE_RADIUS, &sphere_vol, &ierr );
00234     CHKIGEOM;
00235     iGeom_moveEnt( igeom, sphere_vol, SPHERE_CENTER[0], SPHERE_CENTER[1], SPHERE_CENTER[2], &ierr );
00236     CHKIGEOM;
00237 
00238     iBase_EntityHandle sphere_surf;
00239     iBase_EntityHandle* ptr = &sphere_surf;
00240     int size = 0, alloc = 1;
00241     iGeom_getEntAdj( igeom, sphere_vol, iBase_FACE, &ptr, &alloc, &size, &ierr );
00242     CHKIGEOM;
00243     if( size != 1 )
00244     {
00245         std::cerr << "Failed to get spherical surface from iGeom instance" << std::endl;
00246         return 0;
00247     }
00248 
00249     Vector3D bmin, bmax;
00250     iGeom_getEntBoundBox( igeom, sphere_surf, &bmin[0], &bmin[1], &bmin[2], &bmax[0], &bmax[1], &bmax[2], &ierr );
00251     CHKIGEOM;
00252     Vector3D center = 0.5 * ( bmin + bmax );
00253     Vector3D rad    = 0.5 * ( bmax - bmin );
00254     if( ( center - SPHERE_CENTER ).length() > EPS || fabs( rad[0] - SPHERE_RADIUS ) > EPS ||
00255         fabs( rad[1] - SPHERE_RADIUS ) > EPS || fabs( rad[2] - SPHERE_RADIUS ) > EPS )
00256     {
00257         std::cerr << "iGeom implementation returned invalid box for sphere" << std::endl
00258                   << "  Expected Sphere Center: " << SPHERE_CENTER << std::endl
00259                   << "  Expected Sphere Radius: " << SPHERE_RADIUS << std::endl
00260                   << "  Actual Box Minimum    : " << bmin << std::endl
00261                   << "  Actual Box Maximum    : " << bmax << std::endl;
00262         return 0;
00263     }
00264 
00265     return new MsqIGeom( igeom, sphere_surf );
00266 }
00267 
00268 MeshDomain* get_mesquite_domain()
00269 {
00270     return new SphericalDomain( SPHERE_CENTER, SPHERE_RADIUS );
00271 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines