MOAB: Mesh Oriented datABase  (version 5.2.1)
MappingFunction.cpp
Go to the documentation of this file.
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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines