MOAB: Mesh Oriented datABase
(version 5.2.1)
|
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 }