MOAB: Mesh Oriented datABase  (version 5.2.1)
size_adapt_shape.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     (2009) kraftche@cae.wisc.edu
00024 
00025   ***************************************************************** */
00026 
00027 /** \file main.cpp
00028  *  \brief
00029  *  \author Jason Kraftcheck
00030  */
00031 #include "TestUtil.hpp"
00032 
00033 #include <iostream>
00034 using std::cout;
00035 using std::endl;
00036 #include <cstdlib>
00037 #include <cstdio>
00038 #include <vector>
00039 #include <algorithm>
00040 
00041 #include "Mesquite.hpp"
00042 #include "MsqError.hpp"
00043 #include "MeshImpl.hpp"
00044 #include "SizeAdaptShapeWrapper.hpp"
00045 #include "SphericalDomain.hpp"
00046 #include "MsqVertex.hpp"
00047 
00048 using namespace MBMesquite;
00049 
00050 std::string DEFAULT_INPUT = TestDir + "/2D/vtk/quads/untangled/bias-sphere-quads.vtk";
00051 
00052 void help( const char* argv0 )
00053 {
00054     std::cerr << "Usage: " << argv0 << " [<output_file>]" << std::endl
00055               << "  Input file is always: " << DEFAULT_INPUT << std::endl
00056               << "  default is no output file" << std::endl;
00057     exit( 1 );
00058 }
00059 
00060 typedef std::vector< Mesh::ElementHandle > elem_vec_t;
00061 
00062 // Assume mesh is on a sphere of radius 10 with an axis equal to Z.
00063 // Return elements near poles and elements near equator.
00064 void find_z10_extreme_elements( Mesh& mesh, elem_vec_t& polar_elems, elem_vec_t& equatorial_elems, MsqError& err );
00065 
00066 // Get areas of quads and tris
00067 void elem_areas( Mesh& mesh, const elem_vec_t& elems, double& min, double& mean, double& max, MsqError& err );
00068 
00069 int main( int argc, char* argv[] )
00070 {
00071     const char* input_file  = DEFAULT_INPUT.c_str();
00072     const char* output_file = NULL;
00073     switch( argc )
00074     {
00075         default:
00076             help( argv[0] );
00077         case 2:
00078             if( !strcmp( argv[1], "-h" ) ) help( argv[0] );
00079             output_file = argv[1];
00080         case 1:;
00081     }
00082 
00083     /* Read a VTK Mesh file */
00084     MsqPrintError err( cout );
00085     MBMesquite::MeshImpl mesh;
00086     mesh.read_vtk( input_file, err );
00087     if( err )
00088     {
00089         std::cerr << input_file << ": file read failed" << endl;
00090         return 1;
00091     }
00092 
00093     elem_vec_t polar, equatorial;
00094     find_z10_extreme_elements( mesh, polar, equatorial, err );
00095     if( err ) return 1;
00096 
00097     double eq_min, eq_max, eq_mean, pol_min, pol_max, pol_mean;
00098     elem_areas( mesh, polar, pol_min, pol_mean, pol_max, err );
00099     if( err ) return 1;
00100     elem_areas( mesh, equatorial, eq_min, eq_mean, eq_max, err );
00101     if( err ) return 1;
00102 
00103     /* Run optimizer */
00104     SphericalDomain geom( Vector3D( 0, 0, 0 ), 10.0 );
00105     SizeAdaptShapeWrapper smoother( 1e-2 );
00106     MeshDomainAssoc mesh_and_domain = MeshDomainAssoc( &mesh, &geom );
00107     smoother.run_instructions( &mesh_and_domain, err );
00108     if( err ) return 1;
00109 
00110     if( output_file )
00111     {
00112         mesh.write_vtk( output_file, err );
00113         if( err )
00114         {
00115             std::cerr << output_file << ": file write failed" << endl;
00116             return 1;
00117         }
00118         else
00119         {
00120             std::cout << "Wrote file: " << output_file << endl;
00121         }
00122     }
00123 
00124     double eq2_min, eq2_max, eq2_mean, pol2_min, pol2_max, pol2_mean;
00125     elem_areas( mesh, polar, pol2_min, pol2_mean, pol2_max, err );
00126     if( err ) return 1;
00127     elem_areas( mesh, equatorial, eq2_min, eq2_mean, eq2_max, err );
00128     if( err ) return 1;
00129 
00130     double eq_min_pct   = 100 * fabs( eq_min - eq2_min ) / eq_min;
00131     double eq_max_pct   = 100 * fabs( eq_max - eq2_max ) / eq_max;
00132     double eq_mean_pct  = 100 * fabs( eq_mean - eq2_mean ) / eq_mean;
00133     double pol_min_pct  = 100 * fabs( pol_min - pol2_min ) / pol_min;
00134     double pol_max_pct  = 100 * fabs( pol_max - pol2_max ) / pol_max;
00135     double pol_mean_pct = 100 * fabs( pol_mean - pol2_mean ) / pol_mean;
00136 
00137     std::printf( "\n\n" );
00138     std::printf( "AREAS:      Initial     Final       Difference  Change\n" );
00139     std::printf( "Polar:\n" );
00140     std::printf( "  Minimum:  %10f  %10f  %10f  %10f%%\n", pol_min, pol2_min, fabs( pol_min - pol2_min ), pol_min_pct );
00141     std::printf( "  Average:  %10f  %10f  %10f  %10f%%\n", pol_mean, pol2_mean, fabs( pol_mean - pol2_mean ),
00142                  pol_mean_pct );
00143     std::printf( "  Maximum:  %10f  %10f  %10f  %10f%%\n", pol_max, pol2_max, fabs( pol_max - pol2_max ), pol_max_pct );
00144     std::printf( "Equatorial:\n" );
00145     std::printf( "  Minimum:  %10f  %10f  %10f  %10f%%\n", eq_min, eq2_min, fabs( eq_min - eq2_min ), eq_min_pct );
00146     std::printf( "  Average:  %10f  %10f  %10f  %10f%%\n", eq_mean, eq2_mean, fabs( eq_mean - eq2_mean ), eq_mean_pct );
00147     std::printf( "  Maximum:  %10f  %10f  %10f  %10f%%\n", eq_max, eq2_max, fabs( eq_max - eq2_max ), eq_max_pct );
00148 
00149     bool success = pol_min_pct < 6.0 && pol_max_pct < 6.5 && eq_min_pct < 25.0 && eq_max_pct < 6.5;
00150     return !success;
00151 }
00152 
00153 void find_z10_extreme_elements( Mesh& mesh, elem_vec_t& polar_elems, elem_vec_t& equatorial_elems, MsqError& err )
00154 {
00155     elem_vec_t elems;
00156     mesh.get_all_elements( elems, err );MSQ_ERRRTN( err );
00157 
00158     std::vector< Mesh::VertexHandle > verts;
00159     std::vector< MsqVertex > coords;
00160     std::vector< size_t > junk;
00161     for( elem_vec_t::iterator i = elems.begin(); i != elems.end(); ++i )
00162     {
00163         verts.clear();
00164         junk.clear();
00165         mesh.elements_get_attached_vertices( &*i, 1, verts, junk, err );MSQ_ERRRTN( err );
00166         coords.resize( verts.size() );
00167         mesh.vertices_get_coordinates( arrptr( verts ), arrptr( coords ), verts.size(), err );MSQ_ERRRTN( err );
00168 
00169         for( std::vector< MsqVertex >::iterator j = coords.begin(); j != coords.end(); ++j )
00170         {
00171             double z = ( *j )[2];
00172             if( fabs( z ) < 1e-6 )
00173             {
00174                 equatorial_elems.push_back( *i );
00175                 break;
00176             }
00177             else if( fabs( z ) - 10 < 1e-6 )
00178             {
00179                 polar_elems.push_back( *i );
00180                 break;
00181             }
00182         }
00183     }
00184 }
00185 
00186 void elem_areas( Mesh& mesh, const elem_vec_t& elems, double& min, double& mean, double& max, MsqError& err )
00187 {
00188     min  = HUGE_VAL;
00189     max  = -1;
00190     mean = 0.0;
00191 
00192     std::vector< EntityTopology > types( elems.size() );
00193     mesh.elements_get_topologies( arrptr( elems ), arrptr( types ), elems.size(), err );MSQ_ERRRTN( err );
00194 
00195     std::vector< Mesh::VertexHandle > verts;
00196     std::vector< MsqVertex > coords;
00197     std::vector< size_t > junk;
00198     for( size_t i = 0; i < elems.size(); ++i )
00199     {
00200         verts.clear();
00201         junk.clear();
00202         mesh.elements_get_attached_vertices( &elems[i], 1, verts, junk, err );MSQ_ERRRTN( err );
00203         coords.resize( verts.size() );
00204         mesh.vertices_get_coordinates( arrptr( verts ), arrptr( coords ), verts.size(), err );MSQ_ERRRTN( err );
00205 
00206         Vector3D v1, v2;
00207         double area;
00208         if( types[i] == TRIANGLE )
00209         {
00210             assert( coords.size() == 3 );
00211             v1   = coords[1] - coords[0];
00212             v2   = coords[2] - coords[0];
00213             area = 0.5 * ( v1 * v2 ).length();
00214         }
00215         else if( types[i] == QUADRILATERAL )
00216         {
00217             assert( coords.size() == 4 );
00218             v1   = coords[0] + coords[1] - coords[2] - coords[3];
00219             v2   = coords[0] + coords[3] - coords[1] - coords[2];
00220             area = 0.25 * ( v1 * v2 ).length();
00221         }
00222         else
00223         {
00224             MSQ_SETERR( err )
00225             ( "Input file contains volume elements", MsqError::UNSUPPORTED_ELEMENT );
00226             return;
00227         }
00228 
00229         if( min > area ) min = area;
00230         if( max < area ) max = area;
00231         mean += area;
00232     }
00233 
00234     mean /= elems.size();
00235 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines