MOAB: Mesh Oriented datABase
(version 5.2.1)
|
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 (2007) kraftche@cae.wisc.edu 00024 00025 ***************************************************************** */ 00026 00027 /** \file TargetMetricUtil.cpp 00028 * \brief 00029 * \author Jason Kraftcheck 00030 */ 00031 00032 #include "Mesquite.hpp" 00033 #include "TargetMetricUtil.hpp" 00034 #include "MsqMatrix.hpp" 00035 #include "PatchData.hpp" 00036 #include "ElementQM.hpp" 00037 #include "ElemSampleQM.hpp" 00038 00039 namespace MBMesquite 00040 { 00041 00042 void surface_to_2d( const MsqMatrix< 3, 2 >& A, const MsqMatrix< 3, 2 >& W, MsqMatrix< 2, 2 >& W_22, 00043 MsqMatrix< 3, 2 >& RZ ) 00044 { 00045 MsqMatrix< 3, 1 > W1 = W.column( 0 ); 00046 MsqMatrix< 3, 1 > W2 = W.column( 1 ); 00047 MsqMatrix< 3, 1 > nw = W1 * W2; 00048 nw *= 1.0 / length( nw ); 00049 00050 MsqMatrix< 3, 1 > z[2]; 00051 z[0] = W1 * ( 1.0 / length( W1 ) ); 00052 z[1] = nw * z[0]; 00053 MsqMatrix< 3, 2 > Z( z ); 00054 00055 MsqMatrix< 3, 1 > np = A.column( 0 ) * A.column( 1 ); 00056 np *= 1.0 / length( np ); 00057 double dot = np % nw; 00058 MsqMatrix< 3, 1 > nr = ( dot >= 0.0 ) ? nw : -nw; 00059 MsqMatrix< 3, 1 > v = nr * np; 00060 double vlen = length( v ); 00061 if( vlen < DBL_EPSILON ) 00062 { 00063 RZ = Z; // R = I 00064 } 00065 else 00066 { 00067 v *= 1.0 / vlen; 00068 MsqMatrix< 3, 1 > r1[3] = { v, np, v * np }; 00069 MsqMatrix< 3, 1 > r2[3] = { v, nr, v * nr }; 00070 MsqMatrix< 3, 3 > R1( r1 ), R2( r2 ); 00071 RZ = R1 * transpose( R2 ) * Z; 00072 } 00073 00074 W_22 = transpose( Z ) * W; 00075 } 00076 /* 00077 void surface_to_2d( const MsqMatrix<3,2>& App, 00078 const MsqMatrix<3,2>& Wp, 00079 MsqMatrix<2,2>& A, 00080 MsqMatrix<2,2>& W ) 00081 { 00082 MsqMatrix<3,1> Wp1 = Wp.column(0); 00083 MsqMatrix<3,1> Wp2 = Wp.column(1); 00084 MsqMatrix<3,1> nwp = Wp1 * Wp2; 00085 nwp *= 1.0/length(nwp); 00086 00087 MsqMatrix<3,1> z[2]; 00088 z[0] = Wp1 * (1.0 / length( Wp1 )); 00089 z[1] = nwp * z[0]; 00090 MsqMatrix<3,2> Z(z); 00091 W = transpose(Z) * Wp; 00092 00093 MsqMatrix<3,1> npp = App.column(0) * App.column(1); 00094 npp *= 1.0 / length(npp); 00095 double dot = npp % nwp; 00096 MsqMatrix<3,1> nr = (dot >= 0.0) ? nwp : -nwp; 00097 MsqMatrix<3,1> v = nr * npp; 00098 double vlen = length(v); 00099 if (vlen > DBL_EPSILON) { 00100 v *= 1.0 / vlen; 00101 MsqMatrix<3,1> r1[3] = { v, npp, v * npp }, r2[3] = { v, nr, v * nr }; 00102 MsqMatrix<3,3> R1( r1 ), R2( r2 ); 00103 MsqMatrix<3,3> RT = R2 * transpose(R1); 00104 MsqMatrix<3,2> Ap = RT * App; 00105 A = transpose(Z) * Ap; 00106 } 00107 else { 00108 A = transpose(Z) * App; 00109 } 00110 } 00111 */ 00112 00113 static inline void append_elem_samples( PatchData& pd, size_t e, std::vector< size_t >& handles ) 00114 { 00115 NodeSet samples = pd.get_samples( e ); 00116 EntityTopology type = pd.element_by_index( e ).get_element_type(); 00117 size_t curr_size = handles.size(); 00118 handles.resize( curr_size + samples.num_nodes() ); 00119 std::vector< size_t >::iterator i = handles.begin() + curr_size; 00120 for( unsigned j = 0; j < TopologyInfo::corners( type ); ++j ) 00121 if( samples.corner_node( j ) ) *( i++ ) = ElemSampleQM::handle( Sample( 0, j ), e ); 00122 for( unsigned j = 0; j < TopologyInfo::edges( type ); ++j ) 00123 if( samples.mid_edge_node( j ) ) *( i++ ) = ElemSampleQM::handle( Sample( 1, j ), e ); 00124 for( unsigned j = 0; j < TopologyInfo::faces( type ); ++j ) 00125 if( samples.mid_face_node( j ) ) *( i++ ) = ElemSampleQM::handle( Sample( 2, j ), e ); 00126 if( TopologyInfo::dimension( type ) == 3 && samples.mid_region_node() ) 00127 *( i++ ) = ElemSampleQM::handle( Sample( 3, 0 ), e ); 00128 } 00129 00130 void get_sample_pt_evaluations( PatchData& pd, std::vector< size_t >& handles, bool free, MsqError& err ) 00131 { 00132 handles.clear(); 00133 std::vector< size_t > elems; 00134 ElementQM::get_element_evaluations( pd, elems, free, err );MSQ_ERRRTN( err ); 00135 for( std::vector< size_t >::iterator i = elems.begin(); i != elems.end(); ++i ) 00136 append_elem_samples( pd, *i, handles ); 00137 } 00138 00139 void get_elem_sample_points( PatchData& pd, size_t elem, std::vector< size_t >& handles, MsqError& /*err*/ ) 00140 { 00141 handles.clear(); 00142 append_elem_samples( pd, elem, handles ); 00143 } 00144 00145 } // namespace MBMesquite