MOAB: Mesh Oriented datABase
(version 5.2.1)
|
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 DomainUtil.cpp 00028 * \brief 00029 * \author Jason Kraftcheck 00030 */ 00031 00032 #include "Mesquite.hpp" 00033 #include "DomainUtil.hpp" 00034 #include "MeshInterface.hpp" 00035 #include "MsqVertex.hpp" 00036 #include "MsqError.hpp" 00037 00038 namespace MBMesquite 00039 { 00040 namespace DomainUtil 00041 { 00042 00043 void bounding_box( const MsqVertex* coords, size_t num_coords, Vector3D& min, Vector3D& max ) 00044 { 00045 min = max = coords[0]; 00046 for( size_t i = 1; i < num_coords; ++i ) 00047 { 00048 for( int j = 0; j < 3; ++j ) 00049 { 00050 if( coords[i][j] < min[j] ) min[j] = coords[i][j]; 00051 if( coords[i][j] > max[j] ) max[j] = coords[i][j]; 00052 } 00053 } 00054 } 00055 00056 double max_box_extent( const MsqVertex* vertex_array, size_t num_vertices ) 00057 { 00058 Vector3D min, max; 00059 bounding_box( vertex_array, num_vertices, min, max ); 00060 max -= min; 00061 return ( max[0] >= max[1] && max[0] >= max[2] ) ? max[0] : ( max[1] >= max[2] ) ? max[1] : max[2]; 00062 } 00063 00064 void get_fixed_vertices( Mesh* mesh, const Mesh::VertexHandle* verts, size_t num_verts, 00065 std::vector< Mesh::VertexHandle >& fixed_verts, MsqError& err ) 00066 { 00067 std::vector< bool > fixed( num_verts ); 00068 mesh->vertices_get_fixed_flag( verts, fixed, num_verts, err );MSQ_ERRRTN( err ); 00069 for( size_t i = 0; i < num_verts; ++i ) 00070 if( fixed[i] ) fixed_verts.push_back( verts[i] ); 00071 } 00072 00073 bool non_colinear_vertices( const MsqVertex* verts, size_t num_verts, Vector3D coords_out[3], double epsilon ) 00074 { 00075 // This function will attempt to find trhee non-colinear 00076 // vertices from the input list. Further, it will attempt 00077 // to select three such vertices that are relatively far 00078 // apart so as to minimize rounding error in any calculation 00079 // using the results of this function. 00080 00081 // Non-colinear, by definition, must be at least trhee unique points. 00082 if( num_verts < 3 ) return false; 00083 00084 // Begin with the first input vertex 00085 size_t first_idx = 0; 00086 00087 // Choose the vertex furthest from the initial one 00088 size_t second_idx = 1; 00089 double dist_sqr = ( verts[first_idx] - verts[second_idx] ).length_squared(); 00090 for( size_t i = 2; i < num_verts; ++i ) 00091 { 00092 double ds = ( verts[second_idx] - verts[i] ).length_squared(); 00093 if( ds > dist_sqr ) 00094 { 00095 dist_sqr = ds; 00096 second_idx = i; 00097 } 00098 } 00099 // fail if all vertices are coincident 00100 if( dist_sqr <= epsilon * epsilon ) return false; 00101 00102 // re-select the first vertex as the one furthest from the second 00103 for( size_t i = 1; i < num_verts; ++i ) 00104 { 00105 double ds = ( verts[second_idx] - verts[i] ).length_squared(); 00106 if( ds > dist_sqr ) 00107 { 00108 dist_sqr = ds; 00109 first_idx = i; 00110 } 00111 } 00112 00113 // select the third vertex as the one furthest from the line formed 00114 // by the first two vertices 00115 Vector3D b = verts[first_idx]; 00116 Vector3D m = verts[second_idx] - b; 00117 Vector3D mx = m * ( 1.0 / m.length_squared() ); 00118 dist_sqr = -1.0; 00119 size_t third_idx = 0; 00120 for( size_t i = 0; i < num_verts; ++i ) 00121 { 00122 double t = mx % ( verts[i] - b ); 00123 double ds = ( ( b + t * m ) - verts[i] ).length_squared(); 00124 if( ds > dist_sqr ) 00125 { 00126 third_idx = i; 00127 dist_sqr = ds; 00128 } 00129 } 00130 // fail if all vertices are colinear 00131 if( dist_sqr <= epsilon * epsilon ) return false; 00132 00133 coords_out[0] = verts[first_idx]; 00134 coords_out[1] = verts[second_idx]; 00135 coords_out[2] = verts[third_idx]; 00136 return true; 00137 } 00138 00139 bool non_coplanar_vertices( const MsqVertex* verts, size_t num_verts, Vector3D coords_out[4], double epsilon ) 00140 { 00141 // This function will attempt to find four non-coplanar 00142 // vertices from the input list. Further, it will attempt 00143 // to select four such vertices that are relatively far 00144 // apart so as to minimize rounding error in any calculation 00145 // using the results of this function. 00146 00147 // Non-coplanar, by definition, must be at least four unique points. 00148 if( num_verts < 4 ) return false; 00149 00150 // Get three non-colinear vertices 00151 if( !non_colinear_vertices( verts, num_verts, coords_out, epsilon ) ) return false; 00152 00153 // The plane of the first three vertices: 00154 Vector3D norm = ( coords_out[1] - coords_out[0] ) * ( coords_out[2] - coords_out[0] ); 00155 norm /= norm.length(); 00156 double d = -( norm % coords_out[0] ); 00157 00158 // Search for the fourth vertex that is furthest from the plane 00159 // of the first three 00160 double dist = -1.0; 00161 for( size_t i = 0; i < num_verts; ++i ) 00162 { 00163 double disti = fabs( norm % verts[i] + d ); 00164 if( disti > dist ) 00165 { 00166 dist = disti; 00167 coords_out[3] = verts[i]; 00168 } 00169 } 00170 00171 // fail if all vertices are colinear 00172 return ( dist > epsilon ); 00173 } 00174 00175 } // namespace DomainUtil 00176 } // namespace MBMesquite