MOAB: Mesh Oriented datABase  (version 5.4.0)
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,
00049                                                          int input_type,
00050                                                          int output_type,
00051                                                          size_t* index_list,
00052                                                          unsigned num_indices,
00053                                                          MsqError& err )
00054 {
00055     bool in_edges, in_faces, in_region, out_edges, out_faces, out_region;
00056     TopologyInfo::higher_order( topo, input_type, in_edges, in_faces, in_region, err );MSQ_ERRRTN( err );
00057     TopologyInfo::higher_order( topo, output_type, out_edges, out_faces, out_region, err );MSQ_ERRRTN( err );
00058 
00059     // We could probably use TopologyInfo to do this more forward-compatible,
00060     // but for efficiency assume the current ITAPS node ordering, where
00061     // all mid-edge nodes occur before mid-face nodes and the mid-region
00062     // node is always last.
00063 
00064     // If both have mid-region nodes and they don't have the same stuff
00065     // preceeding the mid-region node, then we need to change that index.
00066     bool region_diff = in_region && out_region && ( in_faces != out_faces || in_edges != out_edges );
00067     // If both have mid-face nodes and one has mid-edge nodes and the other
00068     // does not, then we need to change the face indices.
00069     bool face_diff = in_faces && out_faces && in_edges != out_edges;
00070     // if nothing to change, return
00071     if( !face_diff && !region_diff ) return;
00072 
00073     const unsigned corners         = TopologyInfo::corners( topo );
00074     const unsigned edges           = TopologyInfo::edges( topo );
00075     const unsigned faces           = TopologyInfo::faces( topo );
00076     const unsigned in_face_offset  = in_edges ? corners + edges : corners;
00077     const unsigned in_regn_offset  = in_faces ? in_face_offset + faces : in_face_offset;
00078     const unsigned out_face_offset = out_edges ? corners + edges : corners;
00079     const unsigned out_regn_offset = out_faces ? out_face_offset + faces : out_face_offset;
00080 
00081     // In the code below, assertions are used to validate the input
00082     // connectivity data as we assume it is an internal mesquite coding
00083     // error for it to be inconsistent.  True error checking is used
00084     // if the elements are incompatible (the index list for the input
00085     // type contains indices for which there is no correpsonding logical
00086     // node location in the connectivity list of the output element type)
00087     // because that indicates an invalid setup (the combination of element
00088     // type and slave nodes does not result in a reduced element that is
00089     // compatible with the mapping function.)  The latter should probably
00090     // have been caught by the mapping function, but to be safe we check
00091     // again here.
00092 
00093     for( size_t i = 0; i < num_indices; ++i )
00094     {
00095         if( index_list[i] < in_face_offset )
00096         {  // corner or mid-edge node
00097             // nothing to change for these, but check that if it is a mid-edge
00098             // node that the other type also has edges
00099             if( index_list[i] >= corners && !out_edges )
00100             {
00101                 MSQ_SETERR( err )( "Incompatible nodes present.", MsqError::UNSUPPORTED_ELEMENT );
00102                 return;
00103             }
00104         }
00105         else if( index_list[i] < in_regn_offset )
00106         {  // mid-face node
00107             assert( TopologyInfo::dimension( topo ) == 3 || index_list[i] == (unsigned)input_type - 1 );
00108             if( !out_faces )
00109             {
00110                 MSQ_SETERR( err )( "Incompatible nodes present.", MsqError::UNSUPPORTED_ELEMENT );
00111                 return;
00112             }
00113             // working with unsigned type (size_t), so make sure we express this
00114             // such that there are no intermediate negative values.
00115             index_list[i] = index_list[i] + out_face_offset - in_face_offset;
00116         }
00117         else
00118         {  // region
00119             assert( in_region );
00120             assert( TopologyInfo::dimension( topo ) == 3 && index_list[i] == (unsigned)input_type - 1 );
00121             if( !out_region )
00122             {
00123                 MSQ_SETERR( err )( "Incompatible nodes present.", MsqError::UNSUPPORTED_ELEMENT );
00124                 return;
00125             }
00126             // working with unsigned type (size_t), so make sure we express this
00127             // such that there are no intermediate negative values.
00128             index_list[i] = index_list[i] + out_regn_offset - in_regn_offset;
00129         }
00130     }
00131 }
00132 
00133 void MappingFunction2D::jacobian( const PatchData& pd,
00134                                   size_t element_number,
00135                                   NodeSet nodeset,
00136                                   Sample location,
00137                                   size_t* vertex_patch_indices_out,
00138                                   MsqVector< 2 >* d_coeff_d_xi_out,
00139                                   size_t& num_vtx_out,
00140                                   MsqMatrix< 3, 2 >& jacobian_out,
00141                                   MsqError& err ) const
00142 {
00143     const MsqMeshEntity& elem = pd.element_by_index( element_number );
00144     const size_t* conn        = elem.get_vertex_index_array();
00145 
00146     derivatives( location, nodeset, vertex_patch_indices_out, d_coeff_d_xi_out, num_vtx_out, err );MSQ_ERRRTN( err );
00147 
00148     convert_connectivity_indices( elem.node_count(), vertex_patch_indices_out, num_vtx_out, err );MSQ_ERRRTN( err );
00149 
00150     jacobian_out.zero();
00151     size_t w = 0;
00152     for( size_t r = 0; r < num_vtx_out; ++r )
00153     {
00154         size_t i = conn[vertex_patch_indices_out[r]];
00155         MsqMatrix< 3, 1 > coords( pd.vertex_by_index( i ).to_array() );
00156         jacobian_out += coords * transpose( d_coeff_d_xi_out[r] );
00157         if( i < pd.num_free_vertices() )
00158         {
00159             vertex_patch_indices_out[w] = i;
00160             d_coeff_d_xi_out[w]         = d_coeff_d_xi_out[r];
00161             ++w;
00162         }
00163     }
00164     num_vtx_out = w;
00165 }
00166 
00167 void MappingFunction2D::ideal( Sample location, MsqMatrix< 3, 2 >& J, MsqError& err ) const
00168 {
00169     const Vector3D* coords = unit_element( element_topology(), true );
00170     if( !coords )
00171     {
00172         MSQ_SETERR( err )( MsqError::UNSUPPORTED_ELEMENT );
00173         return;
00174     }
00175 
00176     const unsigned MAX_VERTS = 4;
00177     MsqVector< 2 > d_coeff_d_xi[MAX_VERTS];
00178     size_t indices[MAX_VERTS], num_vtx = 0;
00179     derivatives( location, NodeSet(), indices, d_coeff_d_xi, num_vtx, err );MSQ_ERRRTN( err );
00180     assert( num_vtx > 0 && num_vtx <= MAX_VERTS );
00181 
00182     J.zero();
00183     for( size_t r = 0; r < num_vtx; ++r )
00184     {
00185         MsqMatrix< 3, 1 > c( coords[indices[r]].to_array() );
00186         J += c * transpose( d_coeff_d_xi[r] );
00187     }
00188 
00189     double size = sqrt( sqrt( fabs( det( transpose( J ) * J ) ) ) );
00190     assert( size > -1e-15 );  // no negative jacobians for ideal elements!
00191     divide( 1.0, size, size );
00192     J *= size;
00193 }
00194 
00195 void MappingFunction3D::jacobian( const PatchData& pd,
00196                                   size_t element_number,
00197                                   NodeSet nodeset,
00198                                   Sample location,
00199                                   size_t* vertex_patch_indices_out,
00200                                   MsqVector< 3 >* d_coeff_d_xi_out,
00201                                   size_t& num_vtx_out,
00202                                   MsqMatrix< 3, 3 >& jacobian_out,
00203                                   MsqError& err ) const
00204 {
00205     const MsqMeshEntity& elem = pd.element_by_index( element_number );
00206     const size_t* conn        = elem.get_vertex_index_array();
00207 
00208     derivatives( location, nodeset, vertex_patch_indices_out, d_coeff_d_xi_out, num_vtx_out, err );MSQ_ERRRTN( err );
00209 
00210     convert_connectivity_indices( elem.node_count(), vertex_patch_indices_out, num_vtx_out, err );MSQ_ERRRTN( err );
00211 
00212     jacobian_out.zero();
00213     size_t w = 0;
00214     for( size_t r = 0; r < num_vtx_out; ++r )
00215     {
00216         size_t i = conn[vertex_patch_indices_out[r]];
00217         MsqMatrix< 3, 1 > coords( pd.vertex_by_index( i ).to_array() );
00218         jacobian_out += coords * transpose( d_coeff_d_xi_out[r] );
00219         if( i < pd.num_free_vertices() )
00220         {
00221             vertex_patch_indices_out[w] = i;
00222             d_coeff_d_xi_out[w]         = d_coeff_d_xi_out[r];
00223             ++w;
00224         }
00225     }
00226     num_vtx_out = w;
00227 }
00228 
00229 void MappingFunction3D::ideal( Sample location, MsqMatrix< 3, 3 >& J, MsqError& err ) const
00230 {
00231     const Vector3D* coords = unit_element( element_topology(), true );
00232     if( !coords )
00233     {
00234         MSQ_SETERR( err )( MsqError::UNSUPPORTED_ELEMENT );
00235         return;
00236     }
00237 
00238     const unsigned MAX_VERTS = 8;
00239     MsqVector< 3 > d_coeff_d_xi[MAX_VERTS];
00240     size_t indices[MAX_VERTS], num_vtx = 0;
00241     derivatives( location, NodeSet(), indices, d_coeff_d_xi, num_vtx, err );MSQ_ERRRTN( err );
00242     assert( num_vtx > 0 && num_vtx <= MAX_VERTS );
00243 
00244     J.zero();
00245     for( size_t r = 0; r < num_vtx; ++r )
00246     {
00247         MsqMatrix< 3, 1 > c( coords[indices[r]].to_array() );
00248         J += c * transpose( d_coeff_d_xi[r] );
00249     }
00250 
00251     double size = MBMesquite::cbrt( fabs( det( J ) ) );
00252     assert( size > -1e-15 );  // no negative jacobians for ideal elements!
00253     divide( 1.0, size, size );
00254     J *= size;
00255 }
00256 
00257 }  // namespace MBMesquite
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines