MOAB: Mesh Oriented datABase  (version 5.2.1)
PlanarDomain.cpp
Go to the documentation of this file.
00001 /* *****************************************************************
00002     MESQUITE -- The Mesh Quality Improvement Toolkit
00003 
00004     Copyright 2004 Sandia Corporation and Argonne National
00005     Laboratory.  Under the terms of Contract DE-AC04-94AL85000
00006     with Sandia Corporation, the U.S. Government retains certain
00007     rights in 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     diachin2@llnl.gov, djmelan@sandia.gov, mbrewer@sandia.gov,
00024     pknupp@sandia.gov, tleurent@mcs.anl.gov, tmunson@mcs.anl.gov,
00025     kraftche@cae.wisc.edu
00026 
00027   ***************************************************************** */
00028 #include "PlanarDomain.hpp"
00029 #include "MsqError.hpp"
00030 #include "MsqVertex.hpp"
00031 #include "DomainUtil.hpp"
00032 
00033 #include <algorithm>
00034 
00035 MBMesquite::PlanarDomain::~PlanarDomain() {}
00036 
00037 void MBMesquite::PlanarDomain::set_plane( const MBMesquite::Vector3D& normal, const MBMesquite::Vector3D& point )
00038 {
00039     mNormal = normal;
00040     mNormal.normalize();
00041     mCoeff = -( mNormal % point );
00042 }
00043 
00044 void MBMesquite::PlanarDomain::snap_to( MBMesquite::Mesh::VertexHandle /*entity_handle*/, Vector3D& coordinate ) const
00045 {
00046     coordinate -= mNormal * ( mNormal % coordinate + mCoeff );
00047 }
00048 
00049 void MBMesquite::PlanarDomain::vertex_normal_at( MBMesquite::Mesh::VertexHandle /*entity_handle*/,
00050                                                  MBMesquite::Vector3D& coordinate ) const
00051 {
00052     coordinate = mNormal;
00053 }
00054 
00055 void MBMesquite::PlanarDomain::element_normal_at( MBMesquite::Mesh::ElementHandle /*entity_handle*/,
00056                                                   MBMesquite::Vector3D& coordinate ) const
00057 {
00058     coordinate = mNormal;
00059 }
00060 
00061 void MBMesquite::PlanarDomain::vertex_normal_at( const MBMesquite::Mesh::VertexHandle*, Vector3D coords[],
00062                                                  unsigned count, MBMesquite::MsqError& ) const
00063 {
00064     for( unsigned i = 0; i < count; ++i )
00065         coords[i] = mNormal;
00066 }
00067 
00068 void MBMesquite::PlanarDomain::closest_point( MBMesquite::Mesh::VertexHandle, const MBMesquite::Vector3D& position,
00069                                               MBMesquite::Vector3D& closest, MBMesquite::Vector3D& normal,
00070                                               MBMesquite::MsqError& ) const
00071 {
00072     normal  = mNormal;
00073     closest = position - mNormal * ( mNormal % position + mCoeff );
00074 }
00075 
00076 void MBMesquite::PlanarDomain::domain_DoF( const Mesh::VertexHandle*, unsigned short* dof_array, size_t num_vertices,
00077                                            MsqError& ) const
00078 {
00079     std::fill( dof_array, dof_array + num_vertices, 2 );
00080 }
00081 
00082 void MBMesquite::PlanarDomain::flip()
00083 {
00084     mNormal = -mNormal;
00085     mCoeff  = -mCoeff;
00086 }
00087 
00088 void MBMesquite::PlanarDomain::fit_vertices( Mesh* mesh, MsqError& err, double epsilon )
00089 {
00090     // Our goal here is to consider only the boundary (fixed) vertices
00091     // when calculating the plane.  If for some reason the user wants
00092     // to snap a not-quite-planar mesh to a plane during optimization,
00093     // if possible we want to start with the plane containing the fixed
00094     // vertices, as those won't get snapped.  If no vertices are fixed,
00095     // then treat them all as fixed for the purpose calculating the plane
00096     // (consider them all.)
00097 
00098     std::vector< Mesh::VertexHandle > verts, fixed;
00099     mesh->get_all_vertices( verts, err );MSQ_ERRRTN( err );
00100     DomainUtil::get_fixed_vertices( mesh, arrptr( verts ), verts.size(), fixed, err );MSQ_ERRRTN( err );
00101 
00102     bool do_all_verts = true;
00103     if( fixed.size() > 2 )
00104     {
00105         do_all_verts = false;
00106         fit_vertices( mesh, arrptr( fixed ), fixed.size(), err, epsilon );
00107 
00108         // if we failed with only the fixed vertices, try again with all of the
00109         // vertices
00110         if( err )
00111         {
00112             err.clear();
00113             do_all_verts = true;
00114         }
00115     }
00116 
00117     if( do_all_verts )
00118     {
00119         fit_vertices( mesh, arrptr( verts ), verts.size(), err, epsilon );MSQ_ERRRTN( err );
00120     }
00121 
00122     // now count inverted elements
00123     size_t inverted_count = 0, total_count = 0;
00124     std::vector< Mesh::ElementHandle > elems;
00125     std::vector< size_t > junk;
00126     std::vector< MsqVertex > coords;
00127     mesh->get_all_elements( elems, err );MSQ_ERRRTN( err );
00128     for( size_t i = 0; i < elems.size(); ++i )
00129     {
00130 
00131         EntityTopology type;
00132         mesh->elements_get_topologies( &elems[i], &type, 1, err );
00133         if( TopologyInfo::dimension( type ) != 2 ) continue;
00134 
00135         verts.clear();
00136         mesh->elements_get_attached_vertices( arrptr( elems ), 1, verts, junk, err );MSQ_ERRRTN( err );
00137         if( verts.size() < 3 ) continue;
00138 
00139         coords.resize( verts.size() );
00140         mesh->vertices_get_coordinates( arrptr( verts ), arrptr( coords ), 3, err );MSQ_ERRRTN( err );
00141         Vector3D n = ( coords[1] - coords[0] ) * ( coords[2] - coords[0] );
00142         ++total_count;
00143         if( n % mNormal < 0.0 ) ++inverted_count;
00144     }
00145 
00146     // if most elements are inverted, flip normal
00147     if( 2 * inverted_count > total_count ) this->flip();
00148 }
00149 
00150 void MBMesquite::PlanarDomain::fit_vertices( MBMesquite::Mesh* mesh, const MBMesquite::Mesh::VertexHandle* verts,
00151                                              size_t num_verts, MBMesquite::MsqError& err, double epsilon )
00152 {
00153     std::vector< MsqVertex > coords( num_verts );
00154     mesh->vertices_get_coordinates( verts, arrptr( coords ), num_verts, err );MSQ_ERRRTN( err );
00155 
00156     if( epsilon <= 0.0 ) epsilon = DomainUtil::default_tolerance( arrptr( coords ), num_verts );
00157 
00158     Vector3D pts[3];
00159     if( !DomainUtil::non_colinear_vertices( arrptr( coords ), num_verts, pts, epsilon ) )
00160     {
00161         MSQ_SETERR( err )( "All vertices are colinear", MsqError::INVALID_MESH );
00162         return;
00163     }
00164 
00165     this->set_plane( ( pts[1] - pts[0] ) * ( pts[2] - pts[0] ), pts[0] );
00166 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines