MOAB: Mesh Oriented datABase  (version 5.2.1)
DeformingDomainWrapper.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 DeformingDomainWrapper.cpp
00028  *  \brief Implement DeformingDomainWrapper class
00029  *  \author Jason Kraftcheck
00030  */
00031 
00032 #include "DeformingDomainWrapper.hpp"
00033 #include "TagVertexMesh.hpp"
00034 #include "ReferenceMesh.hpp"
00035 #include "RefMeshTargetCalculator.hpp"
00036 #include "InstructionQueue.hpp"
00037 #include "QualityImprover.hpp"
00038 #include "SteepestDescent.hpp"
00039 #include "TerminationCriterion.hpp"
00040 #include "PMeanPTemplate.hpp"
00041 #include "ElementPMeanP.hpp"
00042 #include "TQualityMetric.hpp"
00043 #include "TMixed.hpp"
00044 #include "TShapeNB1.hpp"
00045 #include "TShapeSize2DNB1.hpp"
00046 #include "TShapeSize3DNB1.hpp"
00047 #include "TShapeSizeOrientNB1.hpp"
00048 #include "MeshInterface.hpp"
00049 #include "MsqError.hpp"
00050 #include "MsqVertex.hpp"
00051 #include "QualityAssessor.hpp"
00052 #include "CurveDomain.hpp"
00053 
00054 #include <numeric>  // for std::accumulate
00055 
00056 namespace MBMesquite
00057 {
00058 
00059 const DeformingDomainWrapper::MeshCharacteristic DEFAULT_METRIC_TYPE = DeformingDomainWrapper::SHAPE;
00060 
00061 const bool DEFAULT_CULLING                              = true;
00062 const double DEFAULT_CPU_TIME                           = 0.0;
00063 const double DEFAULT_MOVEMENT_FACTOR                    = 0.01;
00064 const int DEFAULT_INNER_ITERATIONS                      = 2;
00065 const char DEFAULT_CURVE_TAG[]                          = "MesquiteCurveFraction";
00066 const DeformingCurveSmoother::Scheme DEFAULT_CURVE_TYPE = DeformingCurveSmoother::PROPORTIONAL;
00067 
00068 DeformingDomainWrapper::DeformingDomainWrapper()
00069     : metricType( DEFAULT_METRIC_TYPE ), initVertexCoords( "" ), doCulling( DEFAULT_CULLING ),
00070       cpuTime( DEFAULT_CPU_TIME ), movementFactor( DEFAULT_MOVEMENT_FACTOR )
00071 {
00072 }
00073 
00074 DeformingDomainWrapper::~DeformingDomainWrapper() {}
00075 
00076 void DeformingDomainWrapper::store_initial_mesh( Mesh* mesh, MsqError& err )
00077 {
00078     TagVertexMesh tool( err, mesh, false, initVertexCoords );MSQ_ERRRTN( err );
00079     InstructionQueue q;
00080     q.add_tag_vertex_mesh( &tool, err );MSQ_ERRRTN( err );
00081     q.run_instructions( mesh, err );MSQ_ERRRTN( err );
00082     if( initVertexCoords.empty() ) initVertexCoords = tool.get_tag_name();
00083 }
00084 
00085 void DeformingDomainWrapper::set_vertex_movement_limit_factor( double f )
00086 {
00087     assert( f > 0.0 );
00088     movementFactor = f;
00089 }
00090 
00091 void DeformingDomainWrapper::run_wrapper( MeshDomainAssoc* mesh_and_domain, ParallelMesh* pmesh, Settings* settings,
00092                                           QualityAssessor* qa, MsqError& err )
00093 {
00094     if( movementFactor <= 0 )
00095     {
00096         MSQ_SETERR( err )
00097         ( MsqError::INVALID_STATE, "Optimization will not terminate with non-positive movement factor %f",
00098           movementFactor );
00099         return;
00100     }
00101 
00102     Mesh* mesh       = mesh_and_domain->get_mesh();
00103     MeshDomain* geom = mesh_and_domain->get_domain();
00104 
00105     // Move initial mesh to domain in case caller did not
00106     move_to_domain( mesh, geom, err );MSQ_ERRRTN( err );
00107 
00108     const double P = 1.0;
00109     TagVertexMesh init_mesh( err, mesh );MSQ_ERRRTN( err );
00110     ReferenceMesh ref_mesh( &init_mesh );
00111     RefMeshTargetCalculator W( &ref_mesh );
00112 
00113     TShapeNB1 mu_s;
00114     TShapeSize2DNB1 mu_2d_ss;
00115     TShapeSize3DNB1 mu_3d_ss;
00116     TMixed mu_ss( &mu_2d_ss, &mu_3d_ss );
00117     TShapeSizeOrientNB1 mu_sso;
00118     TMetric* mu = 0;
00119     switch( metricType )
00120     {
00121         case SHAPE:
00122             mu = &mu_s;
00123             break;
00124         case SHAPE_SIZE:
00125             mu = &mu_ss;
00126             break;
00127         case SHAPE_SIZE_ORIENT:
00128             mu = &mu_sso;
00129             break;
00130     }
00131 
00132     TQualityMetric sample_metric( &W, mu );
00133     ElementPMeanP elem_metric( P, &sample_metric );
00134     PMeanPTemplate obj_func( P, &elem_metric );
00135 
00136     TerminationCriterion inner, outer;
00137     SteepestDescent improver( &obj_func );
00138     improver.use_element_on_vertex_patch();  // Nash optimization
00139     inner.add_iteration_limit( DEFAULT_INNER_ITERATIONS );
00140     if( doCulling )
00141         inner.cull_on_absolute_vertex_movement_edge_length( movementFactor );
00142     else
00143         outer.add_absolute_vertex_movement_edge_length( movementFactor );
00144     if( cpuTime > 0.0 ) outer.add_cpu_time( cpuTime );
00145     improver.set_inner_termination_criterion( &inner );
00146     improver.set_outer_termination_criterion( &outer );
00147 
00148     qa->add_quality_assessment( &elem_metric );
00149     InstructionQueue q;
00150     q.add_quality_assessor( qa, err );MSQ_ERRRTN( err );
00151     q.set_master_quality_improver( &improver, err );MSQ_ERRRTN( err );
00152     q.add_quality_assessor( qa, err );MSQ_ERRRTN( err );
00153     q.run_common( mesh_and_domain, pmesh, settings, err );MSQ_ERRRTN( err );
00154 }
00155 
00156 void DeformingDomainWrapper::move_to_domain( Mesh* mesh, MeshDomain* geom, MsqError& err )
00157 {
00158     std::vector< Mesh::VertexHandle > verts;
00159     mesh->get_all_vertices( verts, err );MSQ_ERRRTN( err );
00160 
00161     MsqVertex coords;
00162     std::vector< Mesh::VertexHandle >::const_iterator i;
00163     for( i = verts.begin(); i != verts.end(); ++i )
00164     {
00165         mesh->vertices_get_coordinates( &*i, &coords, 1, err );MSQ_ERRRTN( err );
00166         geom->snap_to( *i, coords );
00167         mesh->vertex_set_coordinates( *i, coords, err );MSQ_ERRRTN( err );
00168     }
00169 }
00170 
00171 DeformingCurveSmoother::DeformingCurveSmoother() : metricType( DEFAULT_CURVE_TYPE ), initFractTag( DEFAULT_CURVE_TAG )
00172 {
00173 }
00174 
00175 DeformingCurveSmoother::~DeformingCurveSmoother() {}
00176 
00177 void DeformingCurveSmoother::store_initial_mesh( Mesh* mesh, const Mesh::VertexHandle* verts, int nverts, CurveDomain*,
00178                                                  MsqError& err )
00179 {
00180     if( nverts < 2 )
00181     {
00182         MSQ_SETERR( err )
00183         ( "Invalid curve mesh.  Need at least two end vertices.", MsqError::INVALID_MESH );
00184         return;
00185     }
00186 
00187     // get edge lengths
00188     std::vector< double > vals( nverts - 1 );
00189     MsqVertex coords[2];
00190     int prev = 0;
00191     mesh->vertices_get_coordinates( verts, coords + prev, 1, err );MSQ_ERRRTN( err );
00192     for( int i = 1; i < nverts; ++i )
00193     {
00194         int next = 1 - prev;
00195         mesh->vertices_get_coordinates( verts + i, coords + next, 1, err );MSQ_ERRRTN( err );
00196         vals[i - 1] = ( coords[0] - coords[1] ).length();
00197         prev        = next;
00198     }
00199 
00200     // convert to length fraction before each iterior vertex
00201     // (sum of lengths of adjacent edges over total curve length)
00202     const double total = std::accumulate( vals.begin(), vals.end(), 0.0 );
00203     for( int i = 1; i < nverts - 1; ++i )
00204         vals[i - 1] = vals[i - 1] / total;
00205     vals.resize( nverts - 2 );
00206 
00207     // create tag
00208     TagHandle tag = mesh->tag_create( initFractTag, Mesh::DOUBLE, 1, 0, err );
00209     if( err.error_code() == MsqError::TAG_ALREADY_EXISTS )
00210     {
00211         err.clear();
00212         tag = get_tag( mesh, err );
00213     }
00214     MSQ_ERRRTN( err );
00215 
00216     // store tag data on interior vertices
00217     mesh->tag_set_vertex_data( tag, nverts - 2, verts + 1, &vals[0], err );MSQ_ERRRTN( err );
00218 }
00219 
00220 void DeformingCurveSmoother::smooth_curve( Mesh* mesh, const Mesh::VertexHandle* verts, int nverts, CurveDomain* geom,
00221                                            Scheme, MsqError& err )
00222 {
00223     if( nverts < 2 )
00224     {
00225         MSQ_SETERR( err )
00226         ( "Invalid curve mesh.  Need at least two end vertices.", MsqError::INVALID_MESH );
00227         return;
00228     }
00229 
00230     // Verify that end vertices are on curve.  We cannot snap end vertices
00231     // to ends of curve for application because we don't know where the
00232     // ends of the curve are.
00233     MsqVertex coords[2], coords2;
00234     Mesh::VertexHandle ends[2] = { verts[0], verts[nverts - 1] };
00235     for( int i = 0; i < 2; ++i )
00236     {
00237         mesh->vertices_get_coordinates( ends + i, coords + i, 1, err );MSQ_ERRRTN( err );
00238         geom->position_from_length( coords[i].to_array(), 0, coords2.to_array(), err );MSQ_ERRRTN( err );
00239         if( ( coords[i] - coords2 ).length_squared() > DBL_EPSILON )
00240         {
00241             MSQ_SETERR( err )
00242             ( "Curve end vertices do not line on curve.  Move ends to curve end points first", MsqError::INVALID_MESH );
00243             return;
00244         }
00245     }
00246     const double total = geom->arc_length( coords[0].to_array(), coords[1].to_array(), err );MSQ_ERRRTN( err );
00247 
00248     std::vector< double > vals( nverts - 1 );
00249     if( metricType == EQUAL )
00250     {
00251         std::fill( vals.begin(), vals.end(), 1.0 / ( nverts - 1 ) );
00252         // fracsum = 1.0;
00253     }
00254     else
00255     {  // metricType == PROPORTIONAL
00256         TagHandle tag = get_tag( mesh, err );MSQ_ERRRTN( err );
00257         mesh->tag_get_vertex_data( tag, nverts - 2, verts + 1, &vals[0], err );MSQ_ERRRTN( err );
00258         double sum = std::accumulate( vals.begin(), vals.end() - 1, 0.0 );
00259         if( 1.0 - sum > 1e-8 )
00260             vals.back() = 1.0 - sum;
00261         else
00262         {
00263             vals.back() = *std::min_element( vals.begin(), vals.end() - 1 );
00264             sum += vals.back();
00265             for( size_t i = 0; i < vals.size(); ++i )
00266                 vals[i] /= sum;
00267         }
00268     }
00269 
00270     double frac_sum = 0.0;
00271     for( int i = 1; i < nverts - 1; ++i )
00272     {
00273         frac_sum += vals[i - 1];
00274         geom->position_from_length( coords[0].to_array(), total * frac_sum, coords[1].to_array(), err );MSQ_ERRRTN( err );
00275         mesh->vertex_set_coordinates( verts[i], coords[1], err );MSQ_ERRRTN( err );
00276     }
00277 }
00278 
00279 TagHandle DeformingCurveSmoother::get_tag( Mesh* mesh, MsqError& err )
00280 {
00281     TagHandle h = mesh->tag_get( initFractTag, err );
00282     MSQ_ERRZERO( err );
00283     std::string name;
00284     Mesh::TagType type;
00285     unsigned len;
00286     mesh->tag_properties( h, name, type, len, err );
00287     MSQ_ERRZERO( err );
00288     if( type != Mesh::DOUBLE || len != 1 )
00289     {
00290         MSQ_SETERR( err )
00291         ( MsqError::INVALID_MESH, "Tag \"%s\" exists but is not 1 double value per vertex.", initFractTag.c_str() );
00292     }
00293     return h;
00294 }
00295 
00296 }  // namespace MBMesquite
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines