MOAB: Mesh Oriented datABase
(version 5.2.1)
|
00001 // test comment, remove when done 00002 /* ***************************************************************** 00003 MESQUITE -- The Mesh Quality Improvement Toolkit 00004 00005 Copyright 2007 Sandia National Laboratories. Developed at the 00006 University of Wisconsin--Madison under SNL contract number 00007 624796. The U.S. Government and the University of Wisconsin 00008 retain certain rights to this software. 00009 00010 This library is free software; you can redistribute it and/or 00011 modify it under the terms of the GNU Lesser General Public 00012 License as published by the Free Software Foundation; either 00013 version 2.1 of the License, or (at your option) any later version. 00014 00015 This library is distributed in the hope that it will be useful, 00016 but WITHOUT ANY WARRANTY; without even the implied warranty of 00017 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00018 Lesser General Public License for more details. 00019 00020 You should have received a copy of the GNU Lesser General Public License 00021 (lgpl.txt) along with this library; if not, write to the Free Software 00022 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 00023 00024 (2008) kraftche@cae.wisc.edu 00025 00026 ***************************************************************** */ 00027 00028 /** \file MappingFunction.cpp 00029 * \brief 00030 * \author Jason Kraftcheck 00031 */ 00032 00033 #include "Mesquite.hpp" 00034 #include "MappingFunction.hpp" 00035 #include "MsqError.hpp" 00036 #include "PatchData.hpp" 00037 #include "IdealElements.hpp" 00038 00039 namespace MBMesquite 00040 { 00041 00042 NodeSet MappingFunction::sample_points( NodeSet higher_order ) const 00043 { 00044 higher_order.set_all_corner_nodes( element_topology() ); 00045 return higher_order; 00046 } 00047 00048 void MappingFunction::convert_connectivity_indices_impl( EntityTopology topo, int input_type, int output_type, 00049 size_t* index_list, unsigned num_indices, MsqError& err ) 00050 { 00051 bool in_edges, in_faces, in_region, out_edges, out_faces, out_region; 00052 TopologyInfo::higher_order( topo, input_type, in_edges, in_faces, in_region, err );MSQ_ERRRTN( err ); 00053 TopologyInfo::higher_order( topo, output_type, out_edges, out_faces, out_region, err );MSQ_ERRRTN( err ); 00054 00055 // We could probably use TopologyInfo to do this more forward-compatible, 00056 // but for efficiency assume the current ITAPS node ordering, where 00057 // all mid-edge nodes occur before mid-face nodes and the mid-region 00058 // node is always last. 00059 00060 // If both have mid-region nodes and they don't have the same stuff 00061 // preceeding the mid-region node, then we need to change that index. 00062 bool region_diff = in_region && out_region && ( in_faces != out_faces || in_edges != out_edges ); 00063 // If both have mid-face nodes and one has mid-edge nodes and the other 00064 // does not, then we need to change the face indices. 00065 bool face_diff = in_faces && out_faces && in_edges != out_edges; 00066 // if nothing to change, return 00067 if( !face_diff && !region_diff ) return; 00068 00069 const unsigned corners = TopologyInfo::corners( topo ); 00070 const unsigned edges = TopologyInfo::edges( topo ); 00071 const unsigned faces = TopologyInfo::faces( topo ); 00072 const unsigned in_face_offset = in_edges ? corners + edges : corners; 00073 const unsigned in_regn_offset = in_faces ? in_face_offset + faces : in_face_offset; 00074 const unsigned out_face_offset = out_edges ? corners + edges : corners; 00075 const unsigned out_regn_offset = out_faces ? out_face_offset + faces : out_face_offset; 00076 00077 // In the code below, assertions are used to validate the input 00078 // connectivity data as we assume it is an internal mesquite coding 00079 // error for it to be inconsistent. True error checking is used 00080 // if the elements are incompatible (the index list for the input 00081 // type contains indices for which there is no correpsonding logical 00082 // node location in the connectivity list of the output element type) 00083 // because that indicates an invalid setup (the combination of element 00084 // type and slave nodes does not result in a reduced element that is 00085 // compatible with the mapping function.) The latter should probably 00086 // have been caught by the mapping function, but to be safe we check 00087 // again here. 00088 00089 for( size_t i = 0; i < num_indices; ++i ) 00090 { 00091 if( index_list[i] < in_face_offset ) 00092 { // corner or mid-edge node 00093 // nothing to change for these, but check that if it is a mid-edge 00094 // node that the other type also has edges 00095 if( index_list[i] >= corners && !out_edges ) 00096 { 00097 MSQ_SETERR( err )( "Incompatible nodes present.", MsqError::UNSUPPORTED_ELEMENT ); 00098 return; 00099 } 00100 } 00101 else if( index_list[i] < in_regn_offset ) 00102 { // mid-face node 00103 assert( TopologyInfo::dimension( topo ) == 3 || index_list[i] == (unsigned)input_type - 1 ); 00104 if( !out_faces ) 00105 { 00106 MSQ_SETERR( err )( "Incompatible nodes present.", MsqError::UNSUPPORTED_ELEMENT ); 00107 return; 00108 } 00109 // working with unsigned type (size_t), so make sure we express this 00110 // such that there are no intermediate negative values. 00111 index_list[i] = index_list[i] + out_face_offset - in_face_offset; 00112 } 00113 else 00114 { // region 00115 assert( in_region ); 00116 assert( TopologyInfo::dimension( topo ) == 3 && index_list[i] == (unsigned)input_type - 1 ); 00117 if( !out_region ) 00118 { 00119 MSQ_SETERR( err )( "Incompatible nodes present.", MsqError::UNSUPPORTED_ELEMENT ); 00120 return; 00121 } 00122 // working with unsigned type (size_t), so make sure we express this 00123 // such that there are no intermediate negative values. 00124 index_list[i] = index_list[i] + out_regn_offset - in_regn_offset; 00125 } 00126 } 00127 } 00128 00129 void MappingFunction2D::jacobian( const PatchData& pd, size_t element_number, NodeSet nodeset, Sample location, 00130 size_t* vertex_patch_indices_out, MsqVector< 2 >* d_coeff_d_xi_out, 00131 size_t& num_vtx_out, MsqMatrix< 3, 2 >& jacobian_out, MsqError& err ) const 00132 { 00133 const MsqMeshEntity& elem = pd.element_by_index( element_number ); 00134 const size_t* conn = elem.get_vertex_index_array(); 00135 00136 derivatives( location, nodeset, vertex_patch_indices_out, d_coeff_d_xi_out, num_vtx_out, err );MSQ_ERRRTN( err ); 00137 00138 convert_connectivity_indices( elem.node_count(), vertex_patch_indices_out, num_vtx_out, err );MSQ_ERRRTN( err ); 00139 00140 jacobian_out.zero(); 00141 size_t w = 0; 00142 for( size_t r = 0; r < num_vtx_out; ++r ) 00143 { 00144 size_t i = conn[vertex_patch_indices_out[r]]; 00145 MsqMatrix< 3, 1 > coords( pd.vertex_by_index( i ).to_array() ); 00146 jacobian_out += coords * transpose( d_coeff_d_xi_out[r] ); 00147 if( i < pd.num_free_vertices() ) 00148 { 00149 vertex_patch_indices_out[w] = i; 00150 d_coeff_d_xi_out[w] = d_coeff_d_xi_out[r]; 00151 ++w; 00152 } 00153 } 00154 num_vtx_out = w; 00155 } 00156 00157 void MappingFunction2D::ideal( Sample location, MsqMatrix< 3, 2 >& J, MsqError& err ) const 00158 { 00159 const Vector3D* coords = unit_element( element_topology(), true ); 00160 if( !coords ) 00161 { 00162 MSQ_SETERR( err )( MsqError::UNSUPPORTED_ELEMENT ); 00163 return; 00164 } 00165 00166 const unsigned MAX_VERTS = 4; 00167 MsqVector< 2 > d_coeff_d_xi[MAX_VERTS]; 00168 size_t indices[MAX_VERTS], num_vtx = 0; 00169 derivatives( location, NodeSet(), indices, d_coeff_d_xi, num_vtx, err );MSQ_ERRRTN( err ); 00170 assert( num_vtx > 0 && num_vtx <= MAX_VERTS ); 00171 00172 J.zero(); 00173 for( size_t r = 0; r < num_vtx; ++r ) 00174 { 00175 MsqMatrix< 3, 1 > c( coords[indices[r]].to_array() ); 00176 J += c * transpose( d_coeff_d_xi[r] ); 00177 } 00178 00179 double size = sqrt( sqrt( fabs( det( transpose( J ) * J ) ) ) ); 00180 assert( size > -1e-15 ); // no negative jacobians for ideal elements! 00181 divide( 1.0, size, size ); 00182 J *= size; 00183 } 00184 00185 void MappingFunction3D::jacobian( const PatchData& pd, size_t element_number, NodeSet nodeset, Sample location, 00186 size_t* vertex_patch_indices_out, MsqVector< 3 >* d_coeff_d_xi_out, 00187 size_t& num_vtx_out, MsqMatrix< 3, 3 >& jacobian_out, MsqError& err ) const 00188 { 00189 const MsqMeshEntity& elem = pd.element_by_index( element_number ); 00190 const size_t* conn = elem.get_vertex_index_array(); 00191 00192 derivatives( location, nodeset, vertex_patch_indices_out, d_coeff_d_xi_out, num_vtx_out, err );MSQ_ERRRTN( err ); 00193 00194 convert_connectivity_indices( elem.node_count(), vertex_patch_indices_out, num_vtx_out, err );MSQ_ERRRTN( err ); 00195 00196 jacobian_out.zero(); 00197 size_t w = 0; 00198 for( size_t r = 0; r < num_vtx_out; ++r ) 00199 { 00200 size_t i = conn[vertex_patch_indices_out[r]]; 00201 MsqMatrix< 3, 1 > coords( pd.vertex_by_index( i ).to_array() ); 00202 jacobian_out += coords * transpose( d_coeff_d_xi_out[r] ); 00203 if( i < pd.num_free_vertices() ) 00204 { 00205 vertex_patch_indices_out[w] = i; 00206 d_coeff_d_xi_out[w] = d_coeff_d_xi_out[r]; 00207 ++w; 00208 } 00209 } 00210 num_vtx_out = w; 00211 } 00212 00213 void MappingFunction3D::ideal( Sample location, MsqMatrix< 3, 3 >& J, MsqError& err ) const 00214 { 00215 const Vector3D* coords = unit_element( element_topology(), true ); 00216 if( !coords ) 00217 { 00218 MSQ_SETERR( err )( MsqError::UNSUPPORTED_ELEMENT ); 00219 return; 00220 } 00221 00222 const unsigned MAX_VERTS = 8; 00223 MsqVector< 3 > d_coeff_d_xi[MAX_VERTS]; 00224 size_t indices[MAX_VERTS], num_vtx = 0; 00225 derivatives( location, NodeSet(), indices, d_coeff_d_xi, num_vtx, err );MSQ_ERRRTN( err ); 00226 assert( num_vtx > 0 && num_vtx <= MAX_VERTS ); 00227 00228 J.zero(); 00229 for( size_t r = 0; r < num_vtx; ++r ) 00230 { 00231 MsqMatrix< 3, 1 > c( coords[indices[r]].to_array() ); 00232 J += c * transpose( d_coeff_d_xi[r] ); 00233 } 00234 00235 double size = MBMesquite::cbrt( fabs( det( J ) ) ); 00236 assert( size > -1e-15 ); // no negative jacobians for ideal elements! 00237 divide( 1.0, size, size ); 00238 J *= size; 00239 } 00240 00241 } // namespace MBMesquite