MOAB: Mesh Oriented datABase  (version 5.2.1)
ConicDomain.cpp
Go to the documentation of this file.
00001 /* *****************************************************************
00002     MESQUITE -- The Mesh Quality Improvement Toolkit
00003 
00004     Copyright 2007 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 ConicDomain.cpp
00028  *  \brief
00029  *  \author Jason Kraftcheck
00030  */
00031 
00032 #include "Mesquite.hpp"
00033 #include "ConicDomain.hpp"
00034 
00035 namespace MBMesquite
00036 {
00037 
00038 ConicDomain::~ConicDomain() {}
00039 
00040 void ConicDomain::evaluate( Mesh::VertexHandle, const Vector3D& p_point, Vector3D& closest, Vector3D& normal ) const
00041 {
00042     // translate such that cone point (mPoint) is at origin
00043     Vector3D pt = p_point - mPoint;
00044 
00045     // find the plane containing both the input point an the axis
00046     Vector3D n = mAxis * pt;
00047     double len = n.length();
00048     if( len < 1e-12 )
00049     {
00050         // point on axis
00051         // choose any plane that does't contain the axis
00052         if( fabs( mAxis[0] ) <= fabs( mAxis[1] ) && fabs( mAxis[0] ) < fabs( mAxis[2] ) )
00053             n = mAxis * Vector3D( 1, 0, 0 );
00054         else if( fabs( mAxis[1] ) <= fabs( mAxis[2] ) )
00055             n = mAxis * Vector3D( 0, 1, 0 );
00056         else
00057             n = mAxis * Vector3D( 0, 0, 1 );
00058     }
00059     else
00060     {
00061         n /= len;
00062     }
00063     // Find two points that are the intersection of the plane with the
00064     // circular cross-section of the cone centered at mPoint
00065     Vector3D p1 = mRadius * ( n * mAxis );
00066     Vector3D p2 = -p1;
00067     // Define two lines of intersect between the cone and the plane
00068     // as the two lines from the apex to each of p1 and p2.
00069     Vector3D apex = mHeight * mAxis;
00070     Vector3D v1   = p1 - apex;
00071     Vector3D v2   = p2 - apex;
00072     // Find closest point on each line to input position
00073     double t1 = v1 % ( p_point - apex ) / ( v1 % v1 );
00074     double t2 = v2 % ( p_point - apex ) / ( v2 % v2 );
00075     // Select the closest of the two
00076     p1 = apex + t1 * v1;
00077     p2 = apex + t2 * v2;
00078     double t;
00079     if( ( p1 - p_point ).length_squared() < ( p2 - p_point ).length_squared() )
00080     {
00081         normal  = v1 * n;
00082         closest = p1;
00083         t       = t1;
00084     }
00085     else
00086     {
00087         normal  = v2 * n;
00088         closest = p2;
00089         t       = t2;
00090     }
00091     // If we're below the apex (t > 0), then the normal
00092     // should be in the same direction as the axis.  Otherwise
00093     // it should be in the opposite direction.
00094     if( t * ( normal % mAxis ) < 0.0 ) normal = -normal;
00095     // normalize and translate
00096     normal *= outwardSign / normal.length();
00097     closest += mPoint;
00098 }
00099 
00100 void ConicDomain::snap_to( Mesh::VertexHandle h, Vector3D& v ) const
00101 {
00102     Vector3D p( v ), n;
00103     evaluate( h, p, v, n );
00104 }
00105 
00106 void ConicDomain::vertex_normal_at( Mesh::VertexHandle h, Vector3D& v ) const
00107 {
00108     Vector3D p( v ), l;
00109     evaluate( h, p, l, v );
00110 }
00111 
00112 void ConicDomain::element_normal_at( Mesh::ElementHandle h, Vector3D& v ) const
00113 {
00114     Vector3D p( v ), l;
00115     // NOTE: Explicitly invoke this class's evaluate method for elements.
00116     //       BoundedCylindarDomain overrides evaluate for vertices only.
00117     ConicDomain::evaluate( h, p, l, v );
00118 }
00119 
00120 void ConicDomain::vertex_normal_at( const Mesh::VertexHandle* h, Vector3D coords[], unsigned count, MsqError& ) const
00121 {
00122     for( unsigned i = 0; i < count; ++i )
00123         vertex_normal_at( h[i], coords[i] );
00124 }
00125 
00126 void ConicDomain::closest_point( Mesh::VertexHandle handle, const Vector3D& position, Vector3D& closest,
00127                                  Vector3D& normal, MsqError& ) const
00128 {
00129     evaluate( handle, position, closest, normal );
00130 }
00131 
00132 void ConicDomain::domain_DoF( const Mesh::VertexHandle*, unsigned short* dof_array, size_t count, MsqError& ) const
00133 {
00134     std::fill( dof_array, dof_array + count, (unsigned short)2 );
00135 }
00136 
00137 }  // namespace MBMesquite
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines