MOAB: Mesh Oriented datABase  (version 5.3.1)
DomainUtil.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 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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines