MOAB: Mesh Oriented datABase  (version 5.2.1)
deforming.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) kraftche@cae.wisc.edu
00024 
00025   ***************************************************************** */
00026 
00027 /** \file deform.cpp
00028  *  \brief
00029  *  \author Jason Kraftcheck
00030  */
00031 
00032 #include "TestUtil.hpp"
00033 #include "Mesquite.hpp"
00034 #include "DeformingDomainWrapper.hpp"
00035 #include "MeshImpl.hpp"
00036 #include "MsqError.hpp"
00037 #include "MeshDomain1D.hpp"
00038 #include "PlanarDomain.hpp"
00039 #include "QualityAssessor.hpp"
00040 #include "MsqVertex.hpp"
00041 
00042 #include <iostream>
00043 #include <cstdlib>
00044 #include <algorithm>
00045 
00046 using namespace MBMesquite;
00047 
00048 std::string INPUT_FILE = std::string( STRINGIFY( SRCDIR ) ) + "/sph-10-zsquare.vtk";
00049 // const char INPUT_FILE[] = "test.vtk";
00050 const double Z = 7.0;
00051 // size of new domain
00052 const double HD = 4.5;
00053 
00054 void usage( const char* argv0 )
00055 {
00056     std::cerr << "Usage: " << argv0 << " [-t <file>] [-d <file>] [-f <file>]" << std::endl
00057               << "\t-t : write initial mesh with tags containing reference mesh" << std::endl
00058               << "\t-d : write file containing deformed mesh with smoothed curves " << std::endl
00059               << "\t-f : write file containing final mesh" << std::endl
00060               << std::endl;
00061     std::exit( 1 );
00062 }
00063 
00064 void classify_boundary( Mesh* mesh, Mesh::VertexHandle corners_out[4], std::vector< Mesh::VertexHandle > curves_out[4],
00065                         MsqError& err );
00066 
00067 void cond_write_file( MeshImpl& mesh, const char* filename )
00068 {
00069     if( filename )
00070     {
00071         MsqPrintError err( std::cerr );
00072         mesh.write_vtk( filename, err );
00073         if( MSQ_CHKERR( err ) )
00074         {
00075             std::cerr << filename << ": failed to write file" << std::endl;
00076             exit( 1 );
00077         }
00078         std::cout << "Wrote file: " << filename << std::endl;
00079     }
00080 }
00081 
00082 int main( int argc, char* argv[] )
00083 {
00084     const char* deformed_file = 0;
00085     const char* smoothed_file = 0;
00086     const char* tag_file      = 0;
00087     struct
00088     {
00089         const char* flag;
00090         const char** name_ptr;
00091     } flags[] = { { "-d", &deformed_file }, { "-t", &tag_file }, { "-f", &smoothed_file }, { 0, 0 } };
00092 
00093     for( int i = 1; i < argc; ++i )
00094     {
00095         int j;
00096         for( j = 0; flags[j].flag && strcmp( flags[j].flag, argv[i] ); ++j )
00097             ;
00098         if( !flags[j].flag )
00099         {
00100             std::cerr << "Invalid argument: \"" << argv[i] << '"' << std::endl;
00101             usage( argv[0] );
00102         }
00103         else if( ++i == argc )
00104         {
00105             std::cerr << "Expected argument following \"" << argv[i - 1] << '"' << std::endl;
00106             usage( argv[0] );
00107         }
00108         *( flags[j].name_ptr ) = argv[i];
00109     }
00110 
00111     // load mesh
00112     MsqPrintError err( std::cerr );
00113     MeshImpl mesh;
00114     mesh.read_vtk( INPUT_FILE.c_str(), err );
00115     if( MSQ_CHKERR( err ) ) return 1;
00116 
00117     // find boundary vertices
00118     std::vector< Mesh::VertexHandle > curves[4];
00119     Mesh::VertexHandle corners[4];
00120     classify_boundary( &mesh, corners, curves, err );
00121     if( MSQ_CHKERR( err ) ) return 1;
00122 
00123     // new, "deformed" domain will be an 2HDx2HD planar square
00124     const double corner_coords[][3] = { { -HD, -HD, Z }, { HD, -HD, Z }, { HD, HD, Z }, { -HD, HD, Z } };
00125     LineDomain lines[4]             = { LineDomain( Vector3D( corner_coords[0] ), Vector3D( 1, 0, 0 ) ),
00126                             LineDomain( Vector3D( corner_coords[1] ), Vector3D( 0, 1, 0 ) ),
00127                             LineDomain( Vector3D( corner_coords[2] ), Vector3D( -1, 0, 0 ) ),
00128                             LineDomain( Vector3D( corner_coords[3] ), Vector3D( 0, -1, 0 ) ) };
00129     PlanarDomain surface( PlanarDomain::XY, Z );
00130 
00131     // save initial mesh state
00132     DeformingCurveSmoother curve_tool;
00133     for( int i = 0; i < 4; ++i )
00134     {
00135         curve_tool.store_initial_mesh( &mesh, &curves[i][0], curves[i].size(), &lines[i], err );
00136         if( MSQ_CHKERR( err ) ) return 1;
00137     }
00138     DeformingDomainWrapper wrapper;
00139     wrapper.store_initial_mesh( &mesh, err );
00140     if( MSQ_CHKERR( err ) ) return 1;
00141 
00142     cond_write_file( mesh, tag_file );
00143 
00144     // move corner vertices to new location
00145     for( int i = 0; i < 4; ++i )
00146     {
00147         Vector3D vect( corner_coords[i] );
00148         mesh.vertex_set_coordinates( corners[i], vect, err );
00149         if( MSQ_CHKERR( err ) ) return 1;
00150     }
00151     std::vector< bool > fixed( 4, true );
00152     mesh.vertices_set_fixed_flag( corners, fixed, 4, err );
00153     if( MSQ_CHKERR( err ) ) return 1;
00154 
00155     // smooth curves
00156     for( int i = 0; i < 4; ++i )
00157     {
00158         curve_tool.smooth_curve( &mesh, &curves[i][0], curves[i].size(), &lines[i],
00159                                  DeformingCurveSmoother::PROPORTIONAL, err );
00160         if( MSQ_CHKERR( err ) ) return 1;
00161         fixed.resize( curves[i].size(), true );
00162         mesh.vertices_set_fixed_flag( &curves[i][0], fixed, curves[i].size(), err );
00163         if( MSQ_CHKERR( err ) ) return 1;
00164     }
00165 
00166     cond_write_file( mesh, deformed_file );
00167 
00168     // smooth surface mesh
00169     MeshDomainAssoc mesh_and_domain = MeshDomainAssoc( &mesh, &surface );
00170     wrapper.run_instructions( &mesh_and_domain, err );
00171     if( MSQ_CHKERR( err ) ) return 1;
00172 
00173     cond_write_file( mesh, smoothed_file );
00174     return wrapper.quality_assessor().invalid_elements();
00175 }
00176 
00177 #define CHKMESH( COND, ERR )                                                         \
00178     do                                                                               \
00179     {                                                                                \
00180         if( !( COND ) )                                                              \
00181         {                                                                            \
00182             MSQ_SETERR( err )( "Unexpected mesh topology", MsqError::INVALID_MESH ); \
00183             return;                                                                  \
00184         }                                                                            \
00185     } while( false )
00186 
00187 void classify_boundary( Mesh* mesh, Mesh::VertexHandle corners_out[4], std::vector< Mesh::VertexHandle > curves_out[4],
00188                         MsqError& err )
00189 {
00190     std::vector< Mesh::VertexHandle > verts;
00191     mesh->get_all_vertices( verts, err );MSQ_ERRRTN( err );
00192 
00193     // Find the corner vertex that has negative X and Y coordinates
00194     Mesh::VertexHandle start;
00195     bool have_start = false;
00196     std::vector< Mesh::ElementHandle > elems;
00197     std::vector< size_t > offsets;
00198     for( size_t i = 0; i < verts.size(); ++i )
00199     {
00200         elems.clear();
00201         offsets.clear();
00202         mesh->vertices_get_attached_elements( &verts[i], 1, elems, offsets, err );MSQ_ERRRTN( err );
00203         if( elems.size() == 1 )
00204         {
00205             MsqVertex coords;
00206             mesh->vertices_get_coordinates( &verts[i], &coords, 1, err );MSQ_ERRRTN( err );
00207             if( coords[0] < 0.0 && coords[1] < 0.0 )
00208             {
00209                 CHKMESH( !have_start, err );
00210                 have_start = true;
00211                 start      = verts[i];
00212             }
00213         }
00214     }
00215     CHKMESH( have_start, err );
00216 
00217     // starting at a a corner vertex, find skin vertices
00218     std::vector< Mesh::VertexHandle > boundary;
00219     boundary.push_back( start );
00220     elems.clear();
00221     offsets.clear();
00222     mesh->vertices_get_attached_elements( &start, 1, elems, offsets, err );MSQ_ERRRTN( err );
00223     Mesh::ElementHandle prev = elems.front();
00224     corners_out[0]           = start;
00225     int ncorner              = 1;
00226     do
00227     {
00228         verts.clear();
00229         offsets.clear();
00230         mesh->elements_get_attached_vertices( &prev, 1, verts, offsets, err );
00231         size_t idx = std::find( verts.begin(), verts.end(), boundary.back() ) - verts.begin();
00232         CHKMESH( idx < verts.size(), err );
00233 
00234         Mesh::VertexHandle next = verts[( idx + 1 ) % verts.size()];
00235         elems.clear();
00236         offsets.clear();
00237         mesh->vertices_get_attached_elements( &next, 1, elems, offsets, err );MSQ_ERRRTN( err );
00238         CHKMESH( elems.size() == 1 || elems.size() == 2, err );
00239         if( elems.size() == 2 )
00240         {
00241             idx = std::find( elems.begin(), elems.end(), prev ) - elems.begin();
00242             CHKMESH( idx < 2, err );
00243             prev = elems[1 - idx];
00244         }
00245 
00246         if( elems.size() == 1 )
00247         {
00248             CHKMESH( ncorner < 5, err );
00249             if( ncorner < 4 )  // don't add start vertx twice
00250                 corners_out[ncorner] = next;
00251             ++ncorner;
00252         }
00253 
00254         boundary.push_back( next );
00255     } while( boundary.front() != boundary.back() );
00256     CHKMESH( ncorner == 5, err );
00257 
00258     // find curve vertices
00259     size_t idx = 0;
00260     for( int c = 0; c < 4; ++c )
00261     {
00262         Mesh::VertexHandle s, e;
00263         s = corners_out[c];
00264         e = corners_out[( c + 1 ) % 4];
00265         CHKMESH( idx < boundary.size(), err );
00266         CHKMESH( boundary[idx] == s, err );
00267         for( ; boundary[idx] != e; ++idx )
00268             curves_out[c].push_back( boundary[idx] );
00269         curves_out[c].push_back( e );
00270     }
00271 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines