MOAB: Mesh Oriented datABase  (version 5.2.1)
SphericalDomain.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 "Mesquite.hpp"
00029 #include "SphericalDomain.hpp"
00030 #include "Vector3D.hpp"
00031 #include "MsqError.hpp"
00032 #include "MsqVertex.hpp"
00033 #include "DomainUtil.hpp"
00034 #include "MsqMatrix.hpp"
00035 #include "moab/Util.hpp"
00036 
00037 #ifdef MSQ_HAVE_IEEEFP_H
00038 #include <ieeefp.h>
00039 #endif
00040 
00041 #include <algorithm>
00042 
00043 MBMesquite::SphericalDomain::~SphericalDomain() {}
00044 
00045 void MBMesquite::SphericalDomain::snap_to( Mesh::VertexHandle /*entity_handle*/, Vector3D& coordinate ) const
00046 {
00047     // Get vector center to coordinate, store in coordinate.
00048     coordinate -= mCenter;
00049     // Get distance from center of sphere
00050     double len = coordinate.length();
00051     // Scale vector to have length of radius
00052     coordinate *= mRadius / len;
00053     // If was at center, return arbitrary position on sphere
00054     // (all possitions are equally close)
00055     if( !moab::Util::is_finite( coordinate.x() ) ) coordinate.set( mRadius, 0.0, 0.0 );
00056     // Get position from vector
00057     coordinate += mCenter;
00058 }
00059 
00060 void MBMesquite::SphericalDomain::vertex_normal_at( Mesh::VertexHandle /*entity_handle*/, Vector3D& coordinate ) const
00061 {
00062     // normal is vector from center to input position
00063     coordinate -= mCenter;
00064     // make it a unit vector
00065     double length = coordinate.length();
00066     coordinate /= length;
00067     // if input position was at center, choose same position
00068     // on sphere as snap_to.
00069     if( !moab::Util::is_finite( coordinate.x() ) ) coordinate.set( 1.0, 0.0, 0.0 );
00070 }
00071 void MBMesquite::SphericalDomain::element_normal_at( Mesh::ElementHandle h, Vector3D& coordinate ) const
00072 {
00073     SphericalDomain::vertex_normal_at( h, coordinate );
00074 }
00075 
00076 void MBMesquite::SphericalDomain::vertex_normal_at( const MBMesquite::Mesh::VertexHandle* handle,
00077                                                     MBMesquite::Vector3D coords[], unsigned count,
00078                                                     MBMesquite::MsqError& ) const
00079 {
00080     for( unsigned i = 0; i < count; ++i )
00081         vertex_normal_at( handle[i], coords[i] );
00082 }
00083 
00084 void MBMesquite::SphericalDomain::closest_point( MBMesquite::Mesh::VertexHandle, const MBMesquite::Vector3D& position,
00085                                                  MBMesquite::Vector3D& closest, MBMesquite::Vector3D& normal,
00086                                                  MBMesquite::MsqError& ) const
00087 {
00088     normal = position - mCenter;
00089     normal.normalize();
00090     if( !moab::Util::is_finite( normal.x() ) ) normal.set( 1.0, 0.0, 0.0 );
00091     closest = mCenter + mRadius * normal;
00092 }
00093 
00094 void MBMesquite::SphericalDomain::domain_DoF( const Mesh::VertexHandle*, unsigned short* dof_array, size_t num_vertices,
00095                                               MsqError& ) const
00096 {
00097     std::fill( dof_array, dof_array + num_vertices, 2 );
00098 }
00099 
00100 void MBMesquite::SphericalDomain::fit_vertices( Mesh* mesh, MsqError& err, double epsilon )
00101 {
00102     std::vector< Mesh::VertexHandle > verts;
00103     mesh->get_all_vertices( verts, err );
00104     if( !MSQ_CHKERR( err ) ) fit_vertices( mesh, arrptr( verts ), verts.size(), err, epsilon );
00105 }
00106 
00107 void MBMesquite::SphericalDomain::fit_vertices( Mesh* mesh, const Mesh::VertexHandle* verts, size_t num_verts,
00108                                                 MsqError& err, double epsilon )
00109 {
00110     std::vector< MsqVertex > coords( num_verts );
00111     mesh->vertices_get_coordinates( verts, arrptr( coords ), num_verts, err );MSQ_ERRRTN( err );
00112 
00113     if( epsilon <= 0.0 ) epsilon = DomainUtil::default_tolerance( arrptr( coords ), num_verts );
00114 
00115     Vector3D pts[4];
00116     if( !DomainUtil::non_coplanar_vertices( arrptr( coords ), num_verts, pts, epsilon ) )
00117     {
00118         MSQ_SETERR( err )( "All vertices are co-planar", MsqError::INVALID_MESH );
00119         return;
00120     }
00121 
00122     // solve deterinant form of four-point sphere
00123 
00124     // Define the bottom 4 rows of a 5x5 matrix.  The top
00125     // row contains the variables we are solving for, so just
00126     // fill it with ones.
00127     const double M_vals[25] = { 1,
00128                                 1,
00129                                 1,
00130                                 1,
00131                                 1,
00132                                 pts[0] % pts[0],
00133                                 pts[0][0],
00134                                 pts[0][1],
00135                                 pts[0][2],
00136                                 1,
00137                                 pts[1] % pts[1],
00138                                 pts[1][0],
00139                                 pts[1][1],
00140                                 pts[1][2],
00141                                 1,
00142                                 pts[2] % pts[2],
00143                                 pts[2][0],
00144                                 pts[2][1],
00145                                 pts[2][2],
00146                                 1,
00147                                 pts[3] % pts[3],
00148                                 pts[3][0],
00149                                 pts[3][1],
00150                                 pts[3][2],
00151                                 1 };
00152     MsqMatrix< 5, 5 > M( M_vals );
00153     double M11 = det( MsqMatrix< 4, 4 >( M, 0, 0 ) );
00154     double M12 = det( MsqMatrix< 4, 4 >( M, 0, 1 ) );
00155     double M13 = det( MsqMatrix< 4, 4 >( M, 0, 2 ) );
00156     double M14 = det( MsqMatrix< 4, 4 >( M, 0, 3 ) );
00157     double M15 = det( MsqMatrix< 4, 4 >( M, 0, 4 ) );
00158 
00159     // define the sphere
00160     Vector3D cent( 0.5 * M12 / M11, -0.5 * M13 / M11, 0.5 * M14 / M11 );
00161     this->set_sphere( cent, sqrt( cent % cent - M15 / M11 ) );
00162 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines