MOAB: Mesh Oriented datABase  (version 5.3.0)
PatchData.cpp
Go to the documentation of this file.
00001 /* *****************************************************************
00002     MESQUITE -- The Mesh Quality Improvement Toolkit
00003 
00004     Copyright 2004 Sandia Corporation and Argonne National
00005     Laboratory.  Under the terms of Contract DE-AC04-94AL85000
00006     with Sandia Corporation, the U.S. Government retains certain
00007     rights in 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     diachin2@llnl.gov, djmelan@sandia.gov, mbrewer@sandia.gov,
00024     pknupp@sandia.gov, tleurent@mcs.anl.gov, tmunson@mcs.anl.gov
00025 
00026   ***************************************************************** */
00027 /*!
00028   \file   PatchData.cpp
00029 
00030   \author Thomas Leurent
00031   \author Michael Brewer
00032   \date   2002-01-17
00033 */
00034 
00035 #include "PatchData.hpp"
00036 #include "MsqMeshEntity.hpp"
00037 #include "MsqFreeVertexIndexIterator.hpp"
00038 #include "MeshInterface.hpp"
00039 #include "MsqTimer.hpp"
00040 #include "MsqDebug.hpp"
00041 #include "GlobalPatch.hpp"
00042 #include "PatchIterator.hpp"
00043 #include "ExtraData.hpp"
00044 #include "Settings.hpp"
00045 #include "MappingFunction.hpp"
00046 
00047 #include <list>
00048 #include <vector>
00049 #include <map>
00050 #include <algorithm>
00051 #include <numeric>
00052 #include <functional>
00053 #include <utility>
00054 #include <iostream>
00055 #include <iomanip>
00056 using std::endl;
00057 using std::internal;
00058 using std::left;
00059 using std::list;
00060 using std::map;
00061 using std::ostream;
00062 using std::setfill;
00063 using std::setw;
00064 using std::vector;
00065 
00066 namespace MBMesquite
00067 {
00068 
00069 const Settings PatchData::defaultSettings;
00070 
00071 PatchData::PatchData()
00072     : myMesh( 0 ), myDomain( 0 ), numFreeVertices( 0 ), numSlaveVertices( 0 ), haveComputedInfos( 0 ), dataList( 0 ),
00073       mSettings( &defaultSettings )
00074 {
00075 }
00076 
00077 // Destructor
00078 PatchData::~PatchData()
00079 {
00080     notify_patch_destroyed();
00081 }
00082 
00083 void PatchData::get_minmax_element_unsigned_area( double& min, double& max, MsqError& err )
00084 {
00085     if( !have_computed_info( MAX_UNSIGNED_AREA ) || !have_computed_info( MIN_UNSIGNED_AREA ) )
00086     {
00087         max          = 0;
00088         min          = MSQ_DBL_MAX;
00089         size_t count = num_elements();
00090         for( size_t i = 0; i < count; ++i )
00091         {
00092             double vol;
00093             assert( i < elementArray.size() );
00094             vol = elementArray[i].compute_unsigned_area( *this, err );MSQ_ERRRTN( err );
00095             if( vol > max ) max = vol;
00096             if( vol < min ) min = vol;
00097         }
00098         note_have_info( MAX_UNSIGNED_AREA );
00099         note_have_info( MIN_UNSIGNED_AREA );
00100         computedInfos[MAX_UNSIGNED_AREA] = max;
00101         computedInfos[MIN_UNSIGNED_AREA] = min;
00102     }
00103     else
00104     {
00105         max = computedInfos[MAX_UNSIGNED_AREA];
00106         min = computedInfos[MIN_UNSIGNED_AREA];
00107     }
00108 
00109     if( max <= 0 || min < 0 || min == MSQ_DBL_MAX ) MSQ_SETERR( err )( MsqError::INTERNAL_ERROR );
00110 }
00111 
00112 void PatchData::get_minmax_edge_length( double& min, double& max ) const
00113 {
00114     min = HUGE_VAL;
00115     max = -HUGE_VAL;
00116 
00117     for( size_t i = 0; i < num_elements(); ++i )
00118     {
00119         const MsqMeshEntity& elem = element_by_index( i );
00120         const size_t* conn        = elem.get_vertex_index_array();
00121         const unsigned num_edges  = TopologyInfo::edges( elem.get_element_type() );
00122         for( unsigned e = 0; e < num_edges; ++e )
00123         {
00124             const unsigned* edge = TopologyInfo::edge_vertices( elem.get_element_type(), e );
00125             const double len_sqr =
00126                 ( vertex_by_index( conn[edge[0]] ) - vertex_by_index( conn[edge[1]] ) ).length_squared();
00127             if( len_sqr < min ) min = len_sqr;
00128             if( len_sqr > max ) max = len_sqr;
00129         }
00130     }
00131     min = sqrt( min );
00132     max = sqrt( max );
00133 }
00134 
00135 /*
00136 double PatchData::get_average_Lambda_3d( MsqError &err)
00137 {
00138   double avg;
00139   if (have_computed_info(AVERAGE_DET3D))
00140   {
00141     avg = computedInfos[AVERAGE_DET3D];
00142   }
00143   else
00144   {
00145     avg =0.;
00146     int total_num_corners =0;
00147     Matrix3D A[MSQ_MAX_NUM_VERT_PER_ENT];
00148     for (size_t i=0; i<elementArray.size(); ++i) {
00149       int nve = elementArray[i].corner_count();
00150       elementArray[i].compute_corner_matrices(*this, A, nve, err);
00151       MSQ_ERRZERO(err);
00152       total_num_corners += nve;
00153       for (int c=0; c<nve; ++c) {
00154         avg += TargetCalculator::compute_Lambda(A[c], err);
00155         MSQ_ERRZERO(err);
00156       }
00157     }
00158 
00159     avg = avg / total_num_corners;
00160     computedInfos[AVERAGE_DET3D] = avg;
00161     note_have_info(AVERAGE_DET3D);
00162   }
00163   return avg;
00164 }
00165 */
00166 
00167 /*! \fn PatchData::reorder()
00168    Physically reorder the vertices and elements in the PatchData to improve
00169    the locality of reference.  This method implements a Reverse Breadth First
00170    Search order starting with the vertex furthest from the origin.  Other
00171    orderings can also be implemented.
00172 */
00173 void PatchData::reorder()
00174 {
00175     size_t i, j;
00176     const size_t num_vertex = num_nodes();
00177     const size_t num_elem   = num_elements();
00178 
00179     // Step 1: Clear any cached data that will be invalidated by this
00180     vertexNormalIndices.clear();
00181     normalData.clear();
00182     // vertexDomainDOF.clear();
00183 
00184     // Step 2: Make sure we have vertex-to-element adjacencies
00185     if( !vertAdjacencyArray.size() ) generate_vertex_to_element_data();
00186 
00187     // Step 3: Do breadth-first search
00188     std::vector< bool > visited( num_vertex, false );
00189     std::vector< size_t > vertex_order( num_vertex );
00190     std::vector< size_t >::iterator q1_beg, q1_end, q2_end;
00191     q1_beg = q1_end = q2_end = vertex_order.begin();
00192     // Outer loop will be done once for each disconnected chunk of mesh.
00193     while( q1_beg != vertex_order.end() )
00194     {
00195         // Find vertex furthest from the origin
00196         double max     = -1.0;
00197         size_t vtx_idx = num_vertex;
00198         for( i = 0; i < num_vertex; ++i )
00199             if( !visited[i] )
00200             {
00201                 double dist = vertexArray[i].length_squared();
00202                 if( dist > max )
00203                 {
00204                     max     = dist;
00205                     vtx_idx = i;
00206                 }
00207             }
00208         assert( vtx_idx < num_vertex );
00209 
00210         *q2_end++ = vtx_idx;
00211         ;
00212         visited[vtx_idx] = true;
00213         do
00214         {
00215             q1_end = q2_end;
00216             for( ; q1_beg != q1_end; ++q1_beg )
00217             {
00218                 size_t vtx_adj_offset = vertAdjacencyOffsets[*q1_beg];
00219                 size_t vtx_adj_end    = vertAdjacencyOffsets[*q1_beg + 1];
00220                 for( i = vtx_adj_offset; i < vtx_adj_end; ++i )
00221                 {
00222                     size_t elem = vertAdjacencyArray[i];
00223                     assert( elem < elementArray.size() );
00224                     size_t num_elem_verts = elementArray[elem].node_count();
00225                     size_t* elem_verts    = elementArray[elem].get_vertex_index_array();
00226                     for( j = 0; j < num_elem_verts; ++j )
00227                     {
00228                         size_t elem_vert = elem_verts[j];
00229                         if( !visited[elem_vert] )
00230                         {
00231                             *q2_end++ = elem_vert;
00232                             ;
00233                             visited[elem_vert] = true;
00234                         }
00235                     }
00236                 }
00237             }
00238         } while( q2_end != q1_end );
00239     }
00240 
00241     // Step 4: vertex_order contains the list of current vertex indices
00242     //         in the opposite of the order that they will occur in the
00243     //         reorderd patch.  The following code will construct veretx_map
00244     //         from vertex_order with the following properties
00245     //         - vertex_map will be indexed by the current vertex index and
00246     //           contain the new index of that vertex (inverse of vertex_order)
00247     //         - the vertices will be grouped by their free/slave/fixed flag.
00248     std::vector< size_t > vertex_map( num_vertex );
00249     const size_t fixed_vtx_offset = numFreeVertices + numSlaveVertices;
00250     size_t free_idx = 0, slave_idx = numFreeVertices, fixed_idx = fixed_vtx_offset;
00251     for( i = 1; i <= num_vertex; ++i )
00252     {
00253         size_t vtx_idx = vertex_order[num_vertex - i];
00254         if( vtx_idx < numFreeVertices )
00255             vertex_map[vtx_idx] = free_idx++;
00256         else if( vtx_idx < fixed_vtx_offset )
00257             vertex_map[vtx_idx] = slave_idx++;
00258         else
00259             vertex_map[vtx_idx] = fixed_idx++;
00260     }
00261     // make sure everything adds up
00262     assert( free_idx == numFreeVertices );
00263     assert( slave_idx == fixed_vtx_offset );
00264     assert( fixed_idx == num_vertex );
00265 
00266     // Step 5: compute element permutation
00267     // initialize all to "num_elem" to indicate unvisited
00268     std::vector< size_t > element_map( num_elem, num_elem );
00269     size_t elem_idx = 0;
00270     for( i = 1; i <= num_vertex; ++i )
00271     {
00272         size_t vtx_idx        = vertex_order[num_vertex - i];
00273         size_t vtx_adj_offset = vertAdjacencyOffsets[vtx_idx];
00274         size_t vtx_adj_end    = vertAdjacencyOffsets[vtx_idx + 1];
00275         for( j = vtx_adj_offset; j < vtx_adj_end; ++j )
00276         {
00277             size_t elem = vertAdjacencyArray[j];
00278             if( element_map[elem] == num_elem ) element_map[elem] = elem_idx++;
00279         }
00280     }
00281     // make sure everything adds up
00282     assert( elem_idx == num_elem );
00283 
00284     // Step 6:  Permute the vertices
00285     std::vector< MsqVertex > new_vertex_array( num_vertex );
00286     std::vector< Mesh::VertexHandle > new_vtx_handle_array( num_vertex );
00287     for( i = 0; i < num_vertex; ++i )
00288     {
00289         size_t new_idx                = vertex_map[i];
00290         new_vertex_array[new_idx]     = vertexArray[i];
00291         new_vtx_handle_array[new_idx] = vertexHandlesArray[i];
00292     }
00293     vertexArray.swap( new_vertex_array );
00294     vertexHandlesArray.swap( new_vtx_handle_array );
00295 
00296     // Step 7: Permute the elements and vertex indices for the elements
00297     std::vector< MsqMeshEntity > new_elem_array( num_elem );
00298     std::vector< Mesh::ElementHandle > new_elem_handle_array( num_elem );
00299     for( i = 0; i < num_elem; ++i )
00300     {
00301         assert( i < elementArray.size() );
00302         size_t vert_count  = elementArray[i].node_count();
00303         size_t* conn_array = elementArray[i].get_vertex_index_array();
00304         for( j = 0; j < vert_count; ++j )
00305         {
00306             conn_array[j] = vertex_map[conn_array[j]];
00307         }
00308 
00309         size_t new_idx = element_map[i];
00310         assert( new_idx < num_elem );
00311         new_elem_array[new_idx]        = elementArray[i];
00312         new_elem_handle_array[new_idx] = elementHandlesArray[i];
00313     }
00314     elementArray.swap( new_elem_array );
00315     elementHandlesArray.swap( new_elem_handle_array );
00316 
00317     // Step 8: Clear no-longer-valid vertex-to-element adjacency info.
00318     if( vertAdjacencyOffsets.size() )
00319     {
00320         vertAdjacencyOffsets.clear();
00321         vertAdjacencyArray.clear();
00322         generate_vertex_to_element_data();
00323     }
00324 
00325     notify_new_patch();
00326 }
00327 
00328 /*!
00329    PatchData::move_free_vertices_constrained() moves the free vertices
00330    (see MsqVertex::is_free() ) as specified by the search direction (dk)
00331    and scale factor (step_size). After being moved, the vertices are
00332    projected onto the appropriate geometry.  Fixed vertices are not moved
00333    regardless of their corresponding dk direction.
00334    It is often useful to use the create_coords_momento() function before
00335    calling this function.
00336    Compile with -DMSQ_DBG3 to check that fixed vertices
00337    are not moved with that call.
00338 
00339    \param dk: must be a [nb_vtx] array of Vector3D that contains
00340    the direction in which to move each vertex. Fixed vertices moving
00341    direction should be zero, although fixed vertices will not be
00342    moved regardless of their corresponding dk value.
00343    \param nb_vtx is the number of vertices to move. must corresponds
00344    to the number of vertices in the PatchData.
00345    \param step_size will multiply the moving direction given in dk
00346    for each vertex.
00347   */
00348 void PatchData::move_free_vertices_constrained( Vector3D dk[], size_t nb_vtx, double step_size, MsqError& err )
00349 {
00350     if( nb_vtx != num_free_vertices() )
00351     {
00352         MSQ_SETERR( err )
00353         ( "The directional vector must be of length numVertices.", MsqError::INVALID_ARG );
00354         return;
00355     }
00356 
00357     size_t i;
00358     for( i = 0; i < num_free_vertices(); ++i )
00359     {
00360         vertexArray[i] += ( step_size * dk[i] );
00361         snap_vertex_to_domain( i, err );MSQ_ERRRTN( err );
00362     }
00363 
00364     if( numSlaveVertices )
00365     {
00366         update_slave_node_coordinates( err );MSQ_CHKERR( err );
00367     }
00368 }
00369 
00370 /*! set_free_vertices_constrained is similar to
00371 PatchData::move_free_vertices_constrained() except the original vertex positions
00372 are those stored in the PatchDataVerticesMemento instead of the actual vertex
00373 position stored in the PatchData Vertex array.  The final location is stored
00374 in the PatchData; the PatchDataVerticesMemento is unchanged.
00375 
00376    \param dk: must be a [nb_vtx] array of Vector3D that contains
00377    the direction in which to move each vertex. Fixed vertices moving
00378    direction should be zero, although fixed vertices will not be
00379    moved regardless of their corresponding dk value.
00380    \param nb_vtx is the number of vertices to move. must corresponds
00381    to the number of vertices in the PatchData.
00382    \param step_size will multiply the moving direction given in dk
00383    for each vertex.
00384   */
00385 void PatchData::set_free_vertices_constrained( PatchDataVerticesMemento* memento, Vector3D dk[], size_t nb_vtx,
00386                                                double step_size, MsqError& err )
00387 {
00388     if( memento->originator != this || nb_vtx != num_free_vertices() )
00389     {
00390         MSQ_SETERR( err )( MsqError::INVALID_ARG );
00391         return;
00392     }
00393 
00394     size_t i;
00395     for( i = 0; i < num_free_vertices(); ++i )
00396     {
00397         vertexArray[i] = memento->vertices[i] + ( step_size * dk[i] );
00398         snap_vertex_to_domain( i, err );MSQ_ERRRTN( err );
00399     }
00400 
00401     if( numSlaveVertices )
00402     {
00403         update_slave_node_coordinates( err );MSQ_CHKERR( err );
00404     }
00405 
00406     // Checks that moving direction is zero for fixed vertices.
00407     if( MSQ_DBG( 3 ) )
00408     {
00409         for( i = 0; i < num_nodes(); ++i )
00410         {
00411             if( dk[i].length_squared() != 0.0 )
00412             {
00413                 MSQ_DBGOUT( 3 ) << "dk[" << i << "]: " << dk[i] << endl;
00414                 MSQ_DBGOUT( 3 ) << "moving a fixed vertex." << endl;
00415             }
00416         }
00417     }
00418 }
00419 
00420 static void project_to_plane( Vector3D& vect, const Vector3D& norm )
00421 {
00422     double len_sqr = norm % norm;
00423     if( len_sqr > DBL_EPSILON ) vect -= norm * ( ( norm % vect ) / len_sqr );
00424 }
00425 
00426 void PatchData::project_gradient( std::vector< Vector3D >& gradient, MsqError& err )
00427 {
00428     if( !domain_set() ) return;
00429 
00430     if( gradient.size() != num_free_vertices() )
00431     {
00432         MSQ_SETERR( err )( MsqError::INVALID_ARG );
00433         return;
00434     }
00435 
00436     if( normalData.empty() )
00437     {
00438         update_cached_normals( err );MSQ_ERRRTN( err );
00439     }
00440 
00441     Vector3D norm;
00442     for( size_t i = 0; i < num_free_vertices(); ++i )
00443     {
00444         if( vertexNormalIndices.empty() )
00445         {  // whole mesh on single 2D domain
00446             norm = normalData[i];
00447             project_to_plane( gradient[i], norm );
00448         }
00449         else if( vertexDomainDOF[i] == 2 )
00450         {  // vertex on surface
00451             norm = normalData[vertexNormalIndices[i]];
00452             project_to_plane( gradient[i], norm );
00453         }
00454         else if( vertexDomainDOF[i] == 1 )
00455         {
00456             size_t j, num_elem;
00457             const size_t* elems = get_vertex_element_adjacencies( i, num_elem, err );MSQ_ERRRTN( err );
00458             for( j = 0; j < num_elem; ++j )
00459             {
00460                 if( 2 == TopologyInfo::dimension( element_by_index( elems[j] ).get_element_type() ) )
00461                 {
00462                     norm = vertexArray[i];
00463                     get_domain()->element_normal_at( elementHandlesArray[elems[j]], norm );
00464                     project_to_plane( gradient[i], norm );
00465                 }
00466             }
00467         }
00468     }
00469 }
00470 
00471 /*! Finds the maximum movement (in the distance norm) of the vertices in a
00472   patch.  The previous vertex positions are givena as a
00473   PatchDataVerticesMemento (memento).  The distance squared which each
00474   vertex has moved is then calculated, and the largest of those distances
00475   is returned.  This function only considers the movement of vertices
00476   that are currently 'free'.
00477   \param memento  a memento of this patch's vertex position at some
00478   (prior) time in the optimization.
00479       */
00480 double PatchData::get_max_vertex_movement_squared( PatchDataVerticesMemento* memento, MsqError& )
00481 {
00482     double max_dist = 0.0;
00483     for( size_t i = 0; i < memento->vertices.size(); ++i )
00484     {
00485         double temp_dist = ( vertexArray[i] - memento->vertices[i] ).length_squared();
00486         if( temp_dist > max_dist ) max_dist = temp_dist;
00487     }
00488     return max_dist;
00489 }
00490 
00491 /*!
00492  */
00493 void PatchData::set_all_vertices_soft_fixed( MsqError& /*err*/ )
00494 {
00495     for( size_t i = 0; i < num_free_vertices(); ++i )
00496         vertexArray[i].set_soft_fixed_flag();
00497 }
00498 
00499 /*!
00500  */
00501 void PatchData::set_free_vertices_soft_fixed( MsqError& /*err*/ )
00502 {
00503     for( size_t i = 0; i < num_free_vertices(); ++i )
00504     {
00505         if( vertexArray[i].is_free_vertex() ) vertexArray[i].set_soft_fixed_flag();
00506     }
00507 }
00508 
00509 /*!
00510  */
00511 void PatchData::set_all_vertices_soft_free( MsqError& /*err*/ )
00512 {
00513     for( size_t i = 0; i < num_nodes(); ++i )
00514         vertexArray[i].remove_soft_fixed_flag();
00515 }
00516 
00517 /*! Get coordinates of element vertices, in canonical order.
00518 
00519     \param elem_index The element index in the Patch
00520     \param coords This vector will have the coordinates appended to it.
00521     If necessary, make sure to clear the vector before calling the function.
00522   */
00523 void PatchData::get_element_vertex_coordinates( size_t elem_index, std::vector< Vector3D >& coords, MsqError& /*err*/ )
00524 {
00525     // Check index
00526     if( elem_index >= num_elements() ) return;
00527 
00528     // Ask the element for its vertex indices
00529     assert( elem_index < elementArray.size() );
00530     const size_t* vertex_indices = elementArray[elem_index].get_vertex_index_array();
00531 
00532     // Get the coords for each indicated vertex
00533     size_t num_verts = elementArray[elem_index].vertex_count();
00534     coords.reserve( coords.size() + num_verts );
00535     for( size_t i = 0; i < num_verts; i++ )
00536         coords.push_back( Vector3D( vertexArray[vertex_indices[i]] ) );
00537 }
00538 
00539 /*! This is an inefficient way of retrieving vertex_indices.
00540     Use PatchData::get_element_array followed by
00541     MsqMeshEntity::get_vertex_index_array() if you don't need
00542     to fill an STL vector.
00543 */
00544 void PatchData::get_element_vertex_indices( size_t elem_index, std::vector< size_t >& vertex_indices,
00545                                             MsqError& /*err*/ )
00546 {
00547     // Ask the element for its vertex indices
00548     assert( elem_index < elementArray.size() );
00549     elementArray[elem_index].get_vertex_indices( vertex_indices );
00550 }
00551 
00552 void PatchData::get_vertex_element_indices( size_t vertex_index, std::vector< size_t >& elem_indices, MsqError& err )
00553 {
00554     size_t count;
00555     const size_t* ptr;
00556     ptr = get_vertex_element_adjacencies( vertex_index, count, err );
00557     elem_indices.resize( count );
00558     std::copy( ptr, ptr + count, elem_indices.begin() );
00559 }
00560 
00561 void PatchData::get_vertex_element_indices( size_t vertex_index, unsigned element_dimension,
00562                                             std::vector< size_t >& elem_indices, MsqError& err )
00563 {
00564     elem_indices.clear();
00565     size_t count;
00566     const size_t* ptr;
00567     ptr = get_vertex_element_adjacencies( vertex_index, count, err );
00568     for( const size_t* const end = ptr + count; ptr != end; ++ptr )
00569     {
00570         assert( *ptr < elementArray.size() );
00571         const EntityTopology type = elementArray[*ptr].get_element_type();
00572         const unsigned dim        = TopologyInfo::dimension( type );
00573         if( dim == element_dimension ) elem_indices.push_back( *ptr );
00574     }
00575 }
00576 
00577 const size_t* PatchData::get_vertex_element_adjacencies( size_t vertex_index, size_t& array_len_out, MsqError& )
00578 {
00579     // Make sure we've got the data
00580     if( vertAdjacencyArray.empty() ) { generate_vertex_to_element_data(); }
00581 
00582     const size_t begin = vertAdjacencyOffsets[vertex_index];
00583     const size_t end   = vertAdjacencyOffsets[vertex_index + 1];
00584     array_len_out      = end - begin;
00585     return &vertAdjacencyArray[begin];
00586 }
00587 
00588 /*!
00589     \brief This function fills a vector<size_t> with the indices
00590     to vertices connected to the given vertex by an edge.  If vert_indices
00591     is not initially empty, the function will not delete the current
00592     contents.  Instead, it will append the new indices at the end of
00593     the vector.
00594 
00595 */
00596 void PatchData::get_adjacent_vertex_indices( size_t vertex_index, std::vector< size_t >& vert_indices,
00597                                              MsqError& err ) const
00598 {
00599     bitMap.clear();
00600     bitMap.resize( num_nodes(), false );
00601 
00602     const size_t* conn;
00603     size_t conn_idx, curr_vtx_idx;
00604     const unsigned* adj;
00605     unsigned num_adj, i;
00606     std::vector< MsqMeshEntity >::const_iterator e;
00607     for( e = elementArray.begin(); e != elementArray.end(); ++e )
00608     {
00609         conn     = e->get_vertex_index_array();
00610         conn_idx = std::find( conn, conn + e->node_count(), vertex_index ) - conn;
00611         if( conn_idx == e->node_count() ) continue;
00612 
00613         // If a higher-order node, return corners of side/face
00614         // that node is in the center of.
00615         EntityTopology type = e->get_element_type();
00616         if( conn_idx >= TopologyInfo::corners( type ) && type != POLYGON )
00617         {
00618             unsigned dim, id;
00619             TopologyInfo::side_from_higher_order( type, e->node_count(), conn_idx, dim, id, err );MSQ_ERRRTN( err );
00620             adj = TopologyInfo::side_vertices( type, dim, id, num_adj );
00621         }
00622         else
00623         {
00624             EntityTopology topo = e->get_element_type();
00625             if( topo == POLYGON )
00626             {
00627                 unsigned number_of_nodes = e->node_count();
00628                 num_adj                  = 2;  // always 2 for a polygon
00629                 unsigned vert_adj[2];
00630                 vert_adj[0] = ( conn_idx + 1 ) % number_of_nodes;
00631                 vert_adj[1] = ( conn_idx + number_of_nodes - 1 ) % number_of_nodes;
00632                 for( i = 0; i < num_adj; ++i )
00633                 {
00634                     curr_vtx_idx = conn[vert_adj[i]];  // get index into patch vertex list
00635                     if( !bitMap[curr_vtx_idx] )
00636                     {
00637                         vert_indices.push_back( curr_vtx_idx );
00638                         bitMap[curr_vtx_idx] = true;
00639                     }
00640                 }
00641             }
00642             else
00643             {
00644                 adj = TopologyInfo::adjacent_vertices( topo, conn_idx, num_adj );
00645                 for( i = 0; i < num_adj; ++i )
00646                 {
00647                     curr_vtx_idx = conn[adj[i]];  // get index into patch vertex list
00648                     if( !bitMap[curr_vtx_idx] )
00649                     {
00650                         vert_indices.push_back( curr_vtx_idx );
00651                         bitMap[curr_vtx_idx] = true;
00652                     }
00653                 }
00654             }
00655         }
00656     }
00657 }
00658 
00659 /*! Fills a vector of indices into the entities array. The entities
00660     in the vector are connected the given entity (ent_ind) via an
00661     n-diminsional entity (where 'n' is a given integer).
00662     Thus, if n = 0, the entities must be connected via a vertex.
00663     If n = 1, the entities must be connected via an edge.
00664     If n = 2, the entities must be connected via a two-dimensional element.
00665     NOTE:  if n is 2 and the elements in the entity array are
00666     two-dimensional, no entities should meet this criterion.
00667     The adj_ents vector is cleared at the beginning of the call.
00668 
00669 */
00670 void PatchData::get_adjacent_entities_via_n_dim( int n, size_t ent_ind, std::vector< size_t >& adj_ents, MsqError& err )
00671 {
00672     // reset the vector
00673     adj_ents.clear();
00674     // vertices of this entity (given by ent_ind)
00675     vector< size_t > verts;
00676     // vector to store elements attached to the vertices in verts
00677     vector< size_t > elem_on_vert[MSQ_MAX_NUM_VERT_PER_ENT];
00678     // length of above vectos
00679     int length_elem_on_vert[MSQ_MAX_NUM_VERT_PER_ENT];
00680     // get verts on this element
00681     get_element_vertex_indices( ent_ind, verts, err );
00682     int num_vert = verts.size();
00683     int i        = 0;
00684     int j        = 0;
00685     for( i = 0; i < num_vert; ++i )
00686     {
00687         // get elements on the vertices in verts and the number of vertices
00688         get_vertex_element_indices( verts[i], elem_on_vert[i], err );MSQ_ERRRTN( err );
00689         length_elem_on_vert[i] = elem_on_vert[i].size();
00690     }
00691     // this_ent is the index for an entity which is a candidate to be placed
00692     // into adj_ents
00693     size_t this_ent;
00694     // num of times this_ent has been found in the vectors of entity indices
00695     int counter = 0;
00696     i           = 0;
00697     // loop of each vert on ent_ind
00698     while( i < num_vert )
00699     {
00700         // loop over each ent connected to vert
00701         j = 0;
00702         while( j < length_elem_on_vert[i] )
00703         {
00704             // get candidate element
00705             this_ent = elem_on_vert[i][j];
00706             // if we haven't already consider this ent
00707             if( this_ent != ent_ind )
00708             {
00709                 // if this_ent occurred earlier we would have already considered it
00710                 // so start at i and j+1
00711                 int k1 = i;
00712                 int k2 = j + 1;
00713                 // this_ent has occured once so far
00714                 counter = 1;
00715                 // while k1 < num_vert
00716                 while( k1 < num_vert )
00717                 {
00718                     // loop over entries in the elem on vert vector
00719                     while( k2 < length_elem_on_vert[k1] )
00720                     {
00721                         // if it matches this_ent
00722                         if( elem_on_vert[k1][k2] == this_ent )
00723                         {
00724                             // mark it as 'seen', by making it the same as ent_ind
00725                             // i.e., the entity  passed to us.
00726                             elem_on_vert[k1][k2] = ent_ind;
00727                             ++counter;
00728                             // do not look at remaining elems in this vector
00729                             k2 += length_elem_on_vert[k1];
00730                         }
00731                         else
00732                             ++k2;
00733                     }
00734                     ++k1;
00735                     k2 = 0;
00736                 }
00737                 // if this_ent occured enough times and isn't ent_ind
00738                 if( counter > n && this_ent != ent_ind ) { adj_ents.push_back( this_ent ); }
00739             }
00740             ++j;
00741         }
00742         ++i;
00743     }
00744 }
00745 
00746 /*! \fn PatchData::update_mesh(MsqError &err)
00747 
00748     \brief This function copies to the TSTT mesh  the changes made to the
00749     free vertices / elements of the PatchData object.
00750 
00751 */
00752 void PatchData::update_mesh( MsqError& err, const TagHandle* tag )
00753 {
00754     if( !myMesh )
00755     {
00756         MSQ_SETERR( err )( "Cannot update mesh on temporary patch.\n", MsqError::INTERNAL_ERROR );
00757         return;
00758     }
00759 
00760     const size_t not_fixed = numFreeVertices + numSlaveVertices;
00761     if( tag )
00762     {
00763         for( size_t i = 0; i < not_fixed; ++i )
00764         {
00765             myMesh->tag_set_vertex_data( *tag, 1, &vertexHandlesArray[i], vertexArray[i].to_array(), err );MSQ_ERRRTN( err );
00766         }
00767     }
00768     else
00769     {
00770         for( size_t i = 0; i < not_fixed; ++i )
00771         {
00772             myMesh->vertex_set_coordinates( vertexHandlesArray[i], vertexArray[i], err );MSQ_ERRRTN( err );
00773         }
00774     }
00775 
00776     for( size_t i = 0; i < vertexArray.size(); ++i )
00777     {
00778         myMesh->vertex_set_byte( vertexHandlesArray[i], vertexArray[i].get_flags(), err );MSQ_ERRRTN( err );
00779     }
00780 }
00781 
00782 void PatchData::update_slave_node_coordinates( MsqError& err )
00783 {
00784     // update slave vertices
00785     if( 0 == num_slave_vertices() ) return;
00786 
00787     // Set a mark on every slave vertex.  We'll clear the marks as we
00788     // set the vertex coordinates.  This way we can check that all
00789     // vertices got updated.
00790     const size_t vert_end = num_free_vertices() + num_slave_vertices();
00791     for( size_t i = num_free_vertices(); i < vert_end; ++i )
00792         vertexArray[i].flags() |= MsqVertex::MSQ_MARK;
00793 
00794     // For each element, calculate slave vertex coordinates from
00795     // mapping function.
00796     const int ARR_SIZE = 27;
00797     double coeff[ARR_SIZE];
00798     size_t index[ARR_SIZE];
00799     for( size_t i = 0; i < num_elements(); ++i )
00800     {
00801         MsqMeshEntity& elem  = element_by_index( i );
00802         const int num_corner = elem.vertex_count();
00803         const int num_node   = elem.node_count();
00804         assert( num_node < ARR_SIZE );
00805 
00806         const EntityTopology type       = elem.get_element_type();
00807         const MappingFunction* const mf = get_mapping_function( type );
00808         if( 0 == mf || num_node == num_corner ) continue;
00809 
00810         const size_t* conn = elem.get_vertex_index_array();
00811 
00812         // for each higher-order non-slave node, set bit indicating
00813         // that mapping function is a function of the non-slave node
00814         // coordinates
00815         NodeSet ho_bits = non_slave_node_set( i );
00816 
00817         // for each higher-order slave node
00818         for( int k = num_corner; k < num_node; ++k )
00819         {
00820             if( !is_vertex_slave( conn[k] ) ) continue;
00821 
00822             // check if we already did this one for an adjacent element
00823             MsqVertex& vert = vertexArray[conn[k]];
00824             if( !vert.is_flag_set( MsqVertex::MSQ_MARK ) ) continue;
00825 
00826             // what is this node a mid node of (i.e. face 1, edge 2, etc.)
00827             Sample loc = TopologyInfo::sample_from_node( type, elem.node_count(), k, err );MSQ_ERRRTN( err );
00828 
00829             // evaluate mapping function at logical loction of HO node.
00830             size_t num_coeff;
00831             mf->coefficients( loc, ho_bits, coeff, index, num_coeff, err );MSQ_ERRRTN( err );
00832             mf->convert_connectivity_indices( num_node, index, num_coeff, err );MSQ_ERRRTN( err );
00833 
00834             // calulate new coordinates for slave node
00835             assert( num_coeff > 0 );
00836             vert = coeff[0] * vertex_by_index( conn[index[0]] );
00837             for( size_t j = 1; j < num_coeff; ++j )
00838                 vert += coeff[j] * vertex_by_index( conn[index[j]] );
00839 
00840             // clear mark
00841             vert.flags() &= ~MsqVertex::MSQ_MARK;
00842         }
00843     }
00844 
00845     // make sure we set the coordinates for every slave node
00846     for( size_t i = num_free_vertices(); i < vert_end; ++i )
00847     {
00848         if( vertex_by_index( i ).is_flag_set( MsqVertex::MSQ_MARK ) )
00849         {
00850             MSQ_SETERR( err )
00851             ( MsqError::INVALID_MESH, "No element with mapping function adjacent to slave vertex %lu (%lu)\n",
00852               (unsigned long)i, (unsigned long)get_vertex_handles_array()[i] );
00853             // make sure we finish with all marks cleared
00854             vertexArray[i].flags() &= ~MsqVertex::MSQ_MARK;
00855         }
00856     }
00857 
00858     // snap vertices to domain
00859     if( domain_set() )
00860     {
00861         for( size_t i = num_free_vertices(); i < vert_end; ++i )
00862         {
00863             snap_vertex_to_domain( i, err );MSQ_ERRRTN( err );
00864         }
00865     }
00866 }
00867 
00868 void PatchData::update_slave_node_coordinates( const size_t* elements, size_t num_elems, MsqError& err )
00869 {
00870     // update slave vertices
00871     if( 0 == num_slave_vertices() ) return;
00872 
00873     // set a mark on each vertex so we don't update shared
00874     // vertices more than once.
00875     for( size_t i = 0; i < num_elems; ++i )
00876     {
00877         MsqMeshEntity& elem  = element_by_index( elements[i] );
00878         const int num_corner = elem.vertex_count();
00879         const int num_node   = elem.node_count();
00880         const size_t* conn   = elem.get_vertex_index_array();
00881         for( int j = num_corner; j < num_node; ++j )
00882             vertexArray[conn[j]].flags() |= MsqVertex::MSQ_MARK;
00883     }
00884 
00885     // For each element, calculate slave vertex coordinates from
00886     // mapping function.
00887     const int ARR_SIZE = 27;
00888     double coeff[ARR_SIZE];
00889     size_t index[ARR_SIZE];
00890     for( size_t i = 0; i < num_elems; ++i )
00891     {
00892         MsqMeshEntity& elem  = element_by_index( elements[i] );
00893         const int num_corner = elem.vertex_count();
00894         const int num_node   = elem.node_count();
00895         assert( num_node < ARR_SIZE );
00896 
00897         const EntityTopology type       = elem.get_element_type();
00898         const MappingFunction* const mf = get_mapping_function( type );
00899         if( 0 == mf || num_node == num_corner ) continue;
00900 
00901         const size_t* conn = elem.get_vertex_index_array();
00902 
00903         // for each higher-order non-slave node, set bit indicating
00904         // that mapping function is a function of the non-slave node
00905         // coordinates
00906         NodeSet ho_bits = non_slave_node_set( i );
00907 
00908         // for each higher-order slave node
00909         for( int k = num_corner; k < num_node; ++k )
00910         {
00911             if( !is_vertex_slave( conn[k] ) ) continue;
00912 
00913             // check if we already did this one for an adjacent element
00914             MsqVertex& vert = vertexArray[conn[k]];
00915             if( !vert.is_flag_set( MsqVertex::MSQ_MARK ) ) continue;
00916 
00917             // what is this node a mid node of (i.e. face 1, edge 2, etc.)
00918             Sample loc = TopologyInfo::sample_from_node( type, elem.node_count(), k, err );MSQ_ERRRTN( err );
00919 
00920             // evaluate mapping function at logical loction of HO node.
00921             size_t num_coeff;
00922             mf->coefficients( loc, ho_bits, coeff, index, num_coeff, err );MSQ_ERRRTN( err );
00923             mf->convert_connectivity_indices( num_node, index, num_coeff, err );MSQ_ERRRTN( err );
00924 
00925             // calulate new coordinates for slave node
00926             assert( num_coeff > 0 );
00927             vert = coeff[0] * vertex_by_index( conn[index[0]] );
00928             for( size_t j = 1; j < num_coeff; ++j )
00929                 vert += coeff[j] * vertex_by_index( conn[index[j]] );
00930 
00931             // snap vertices to domain
00932             if( domain_set() )
00933             {
00934                 snap_vertex_to_domain( conn[k], err );MSQ_ERRRTN( err );
00935             }
00936 
00937             // clear mark
00938             vert.flags() &= ~MsqVertex::MSQ_MARK;
00939         }
00940     }
00941 }
00942 
00943 void PatchData::generate_vertex_to_element_data()
00944 {
00945     MSQ_FUNCTION_TIMER( "PatchData::generate_vertex_to_element_data" );
00946 
00947     // Skip if data already exists
00948     if( !vertAdjacencyArray.empty() ) return;
00949 
00950     // Skip if patch is empty
00951     if( 0 == num_nodes() ) return;
00952 
00953     // Allocate offset array
00954     vertAdjacencyOffsets.clear();
00955     vertAdjacencyOffsets.resize( num_nodes() + 1, 0 );
00956 
00957     // Temporarily use offsets array to hold per-vertex element count
00958     std::vector< MsqMeshEntity >::iterator elem_iter;
00959     const std::vector< MsqMeshEntity >::iterator elem_end = elementArray.end();
00960     for( elem_iter = elementArray.begin(); elem_iter != elem_end; ++elem_iter )
00961     {
00962         size_t* conn_iter      = elem_iter->get_vertex_index_array();
00963         const size_t* conn_end = conn_iter + elem_iter->node_count();
00964         for( ; conn_iter != conn_end; ++conn_iter )
00965             ++vertAdjacencyOffsets[*conn_iter];
00966     }
00967 
00968     // Convert counts to end indices.
00969     // When done, vertAdjacencyOffsets will contain, for each vertex,
00970     // one more than the *last* index for that vertex's data in the
00971     // adjacency array.  This is *not* the final state for this data.
00972     // See comments for next loop.
00973     std::vector< size_t >::iterator off_iter      = vertAdjacencyOffsets.begin();
00974     const std::vector< size_t >::iterator off_end = vertAdjacencyOffsets.end();
00975     size_t prev                                   = *off_iter;
00976     ++off_iter;
00977     for( ; off_iter != off_end; ++off_iter )
00978     {
00979         prev += *off_iter;
00980         *off_iter = prev;
00981     }
00982 
00983     // Allocate space for element numbers
00984     const size_t num_vert_uses = vertAdjacencyOffsets[num_nodes() - 1];
00985     assert( num_vert_uses == elemConnectivityArray.size() );
00986     vertAdjacencyArray.resize( num_vert_uses );
00987 
00988     // Fill vertAdjacencyArray, using the indices in vertAdjacencyOffsets
00989     // as the location to insert the next element number in
00990     // vertAdjacencyArray.  When done, vertAdjacenyOffsets will contain
00991     // the start index for each vertex, rather than one past the last
00992     // index.
00993     for( size_t i = 0; i < elementArray.size(); ++i )
00994     {
00995         size_t* conn_iter      = elementArray[i].get_vertex_index_array();
00996         const size_t* conn_end = conn_iter + elementArray[i].node_count();
00997         for( ; conn_iter != conn_end; ++conn_iter )
00998         {
00999             const size_t array_index        = --vertAdjacencyOffsets[*conn_iter];
01000             vertAdjacencyArray[array_index] = i;
01001         }
01002     }
01003 
01004     // Last entry should be number of vertex uses (one past the
01005     // last index of the last vertex.)
01006     vertAdjacencyOffsets[num_nodes()] = num_vert_uses;
01007 }
01008 
01009 void PatchData::get_subpatch( size_t center_vertex_index, unsigned num_adj_elem_layers, PatchData& subpatch,
01010                               MsqError& err )
01011 {
01012     // Make sure we're in range
01013     if( center_vertex_index >= num_free_vertices() )
01014     {
01015         MSQ_SETERR( err )( "Invalid index for center vertex", MsqError::INVALID_ARG );
01016         return;
01017     }
01018 
01019     // Notify any observers of the existing subpatch that the mesh
01020     // in the patch is to be changed.
01021     subpatch.notify_new_patch();
01022 
01023     // Get list of vertices and elements in subpatch.
01024     // Ultimately, end up with arrays of unique, sorted indices.
01025     // It is important that the vertex indices be sorted so later
01026     // a reverse lookup can be done using a binary search (std::lower_bound).
01027     std::vector< size_t > elements, vertices, offsets;
01028     vertices.push_back( center_vertex_index );
01029     for( unsigned i = 0; i < num_adj_elem_layers; ++i )
01030     {
01031         elements.clear();
01032         for( unsigned v = 0; v < vertices.size(); ++v )
01033         {
01034             size_t num_elem;
01035             const size_t* vert_elems = get_vertex_element_adjacencies( vertices[v], num_elem, err );MSQ_ERRRTN( err );
01036             elements.insert( elements.end(), vert_elems, vert_elems + num_elem );
01037         }
01038         std::sort( elements.begin(), elements.end() );
01039         elements.erase( std::unique( elements.begin(), elements.end() ), elements.end() );
01040 
01041         vertices.clear();
01042         for( unsigned e = 0; e < elements.size(); ++e )
01043         {
01044             MsqMeshEntity& elem      = element_by_index( elements[e] );
01045             size_t num_vert          = elem.node_count();
01046             const size_t* elem_verts = elem.get_vertex_index_array();
01047             vertices.insert( vertices.end(), elem_verts, elem_verts + num_vert );
01048         }
01049         std::sort( vertices.begin(), vertices.end() );
01050         vertices.erase( std::unique( vertices.begin(), vertices.end() ), vertices.end() );
01051     }
01052 
01053     // Allocate space for element connectivity info.
01054     size_t num_vert_uses = 0;
01055     for( unsigned i = 0; i < elements.size(); ++i )
01056         num_vert_uses += element_by_index( elements[i] ).node_count();
01057     subpatch.elementArray.resize( elements.size() );
01058     subpatch.elementHandlesArray.resize( elements.size() );
01059     subpatch.elemConnectivityArray.resize( num_vert_uses );
01060     offsets.resize( elements.size() + 1 );
01061 
01062     // Construct element connectivity data in new patch,
01063     // and copy element type into new patch
01064     size_t curr_offset = 0;
01065     for( unsigned i = 0; i < elements.size(); ++i )
01066     {
01067         MsqMeshEntity& elem = element_by_index( elements[i] );
01068         assert( i < elementArray.size() );
01069         subpatch.elementArray[i].set_element_type( elem.get_element_type() );
01070         subpatch.elementHandlesArray[i] = elementHandlesArray[elements[i]];
01071         const size_t* verts             = elem.get_vertex_index_array();
01072         offsets[i]                      = curr_offset;
01073         for( unsigned j = 0; j < elem.node_count(); ++j )
01074         {
01075             subpatch.elemConnectivityArray[curr_offset++] =
01076                 std::lower_bound( vertices.begin(), vertices.end(), verts[j] ) - vertices.begin();
01077         }
01078     }
01079     offsets[elements.size()] = curr_offset;
01080 
01081     // Store index in this patch in vertex handle array of subpatch
01082     // so we can determine how vertices were reordered when setting
01083     // vertex coordinates.
01084     assert( sizeof( size_t ) == sizeof( void* ) );
01085     subpatch.vertexHandlesArray.resize( vertices.size() );
01086     size_t* vert_handles = reinterpret_cast< size_t* >( &subpatch.vertexHandlesArray[0] );
01087     std::copy( vertices.begin(), vertices.end(), vert_handles );
01088 
01089     // All vertices except vertex at center_vertex_index are fixed.
01090     subpatch.byteArray.resize( vertices.size() );
01091     for( size_t pi = 0; pi < vertices.size(); ++pi )
01092     {
01093         if( vertices[pi] == center_vertex_index )
01094             subpatch.byteArray[pi] = vertexArray[vertices[pi]].get_flags() & ~MsqVertex::MSQ_PATCH_FIXED;
01095         else
01096             subpatch.byteArray[pi] = vertexArray[vertices[pi]].get_flags() | MsqVertex::MSQ_PATCH_FIXED;
01097     }
01098 
01099     // Re-order vertices and initialize other data in subpatch
01100     subpatch.initialize_data( arrptr( offsets ), &subpatch.byteArray[0], err );MSQ_ERRRTN( err );
01101 
01102     // Copy vertex data into subpatch.  subpatch.vertexHandlesArray contains
01103     // the indices into this PatchData for each vertex, as reordered by the
01104     // call to initialize_data.
01105     subpatch.vertexArray.resize( vertices.size() );
01106     for( unsigned i = 0; i < vertices.size(); ++i )
01107     {
01108         size_t vert_index               = ( size_t )( subpatch.vertexHandlesArray[i] );
01109         vertices[i]                     = vert_index;
01110         subpatch.vertexHandlesArray[i]  = vertexHandlesArray[vert_index];
01111         subpatch.vertexArray[i]         = vertexArray[vert_index];
01112         subpatch.vertexArray[i].flags() = subpatch.byteArray[i];
01113     }
01114 
01115     subpatch.myMesh    = myMesh;
01116     subpatch.myDomain  = myDomain;
01117     subpatch.mSettings = mSettings;
01118 
01119     notify_sub_patch( subpatch, arrptr( vertices ), elements.empty() ? 0 : arrptr( elements ), err );MSQ_CHKERR( err );
01120 }
01121 
01122 //! Adjust the position of the specified vertex so that it
01123 //! lies on its constraining domain.  The actual domain constraint
01124 //! is managed by the TSTT mesh implementation
01125 void PatchData::snap_vertex_to_domain( size_t vertex_index, MsqError& err )
01126 {
01127     if( domain_set() )
01128     {
01129         // if not doing normal caching or vertex is not on a single surface
01130         if( normalData.empty() )
01131         {
01132             get_domain()->snap_to( vertexHandlesArray[vertex_index], vertexArray[vertex_index] );
01133         }
01134         // entire domain is 2-D (all vertices have a single normal)
01135         else if( vertexNormalIndices.empty() )
01136         {
01137             get_domain()->closest_point( vertexHandlesArray[vertex_index], Vector3D( vertexArray[vertex_index] ),
01138                                          vertexArray[vertex_index], normalData[vertex_index], err );MSQ_ERRRTN( err );
01139         }
01140         else if( vertexNormalIndices[vertex_index] < normalData.size() )
01141         {  // vertex has a unique normal
01142             get_domain()->closest_point( vertexHandlesArray[vertex_index], Vector3D( vertexArray[vertex_index] ),
01143                                          vertexArray[vertex_index], normalData[vertexNormalIndices[vertex_index]],
01144                                          err );MSQ_ERRRTN( err );
01145         }
01146         else
01147         {  // vertex has multiple normals
01148             get_domain()->snap_to( vertexHandlesArray[vertex_index], vertexArray[vertex_index] );
01149         }
01150     }
01151 }
01152 
01153 void PatchData::update_cached_normals( MsqError& err )
01154 {
01155     size_t i;
01156 
01157     MeshDomain* domain = get_domain();
01158     if( !domain )
01159     {
01160         MSQ_SETERR( err )( "No domain constraint set.", MsqError::INVALID_STATE );
01161         return;
01162     }
01163 
01164     // Determine which vertices lie on surfaces
01165     vertexDomainDOF.resize( num_nodes() );
01166     domain->domain_DoF( arrptr( vertexHandlesArray ), arrptr( vertexDomainDOF ), num_nodes(), err );MSQ_ERRRTN( err );
01167 
01168     // Count how many vertices have a single normal
01169     // Sun doesn't support partial template specialization, so can't use std::count
01170     // size_t n = std::count( vertexDomainDOF.begin(), vertexDomainDOF.end(), 2 );
01171     size_t n = 0;
01172     std::vector< unsigned short >::iterator k;
01173     for( k = vertexDomainDOF.begin(); k != vertexDomainDOF.end(); ++k )
01174         if( *k == 2 ) ++n;
01175 
01176     normalData.resize( n );
01177 
01178     // If all vertices are on a surface, pass in the existing handles array
01179     // and store a single normal per vertex.
01180     if( n == num_nodes() )
01181     {
01182         std::copy( vertexArray.begin(), vertexArray.end(), normalData.begin() );
01183         domain->vertex_normal_at( arrptr( vertexHandlesArray ), arrptr( normalData ), num_nodes(), err );
01184         vertexNormalIndices.clear();
01185         vertexDomainDOF.clear();MSQ_ERRRTN( err );
01186     }
01187     else
01188     {
01189         vertexNormalIndices.resize( num_nodes() );
01190         size_t nn = 0;
01191         for( i = 0; i < num_nodes(); ++i )
01192         {
01193             if( vertexDomainDOF[i] == 2 )
01194             {
01195                 normalData[nn] = vertexArray[i];
01196                 domain->vertex_normal_at( vertexHandlesArray[i], normalData[nn] );
01197                 vertexNormalIndices[i] = nn;
01198                 ++nn;
01199             }
01200             else
01201             {
01202                 vertexNormalIndices[i] = n + 1;
01203             }
01204         }
01205         assert( nn == n );
01206     }
01207 }
01208 
01209 void PatchData::get_domain_normal_at_element( size_t elem_index, Vector3D& surf_norm, MsqError& err )
01210 {
01211     // check if element as a mid-face node
01212     const MsqMeshEntity& elem = element_by_index( elem_index );
01213     const int mid_node = TopologyInfo::higher_order_from_side( elem.get_element_type(), elem.node_count(), 2, 0, err );MSQ_ERRRTN( err );
01214     // if we have the mid-element vertex, get cached normal for it
01215     if( mid_node > 0 )
01216     {
01217         get_domain_normal_at_vertex( elem.get_vertex_index_array()[mid_node], elementHandlesArray[elem_index],
01218                                      surf_norm, err );MSQ_ERRRTN( err );
01219     }
01220     // otherwise query domain for normal at element centroid
01221     else if( domain_set() )
01222     {
01223         assert( elem_index < elementArray.size() );
01224         elementArray[elem_index].get_centroid( surf_norm, *this, err );MSQ_ERRRTN( err );
01225         get_domain()->element_normal_at( elementHandlesArray[elem_index], surf_norm );
01226     }
01227     else
01228         MSQ_SETERR( err )( "No domain constraint set.", MsqError::INVALID_STATE );
01229 }
01230 
01231 void PatchData::get_domain_normal_at_mid_edge( size_t elem_index, unsigned edge_num, Vector3D& normal, MsqError& err )
01232 {
01233     // check if element as a mid-edge node
01234     const MsqMeshEntity& elem = element_by_index( elem_index );
01235     const int mid_node =
01236         TopologyInfo::higher_order_from_side( elem.get_element_type(), elem.node_count(), 1, edge_num, err );MSQ_ERRRTN( err );
01237     // if we have the mid-edge vertex, get cached normal for it
01238     if( mid_node > 0 )
01239     {
01240         get_domain_normal_at_vertex( elem.get_vertex_index_array()[mid_node], elementHandlesArray[elem_index], normal,
01241                                      err );MSQ_ERRRTN( err );
01242     }
01243     // otherwise query domain for normal at center of edge
01244     else if( domain_set() )
01245     {
01246         const unsigned* edge = TopologyInfo::edge_vertices( elem.get_element_type(), edge_num, err );MSQ_ERRRTN( err );
01247         const MsqVertex& v1 = vertex_by_index( elem.get_vertex_index_array()[edge[0]] );
01248         const MsqVertex& v2 = vertex_by_index( elem.get_vertex_index_array()[edge[1]] );
01249         normal              = 0.5 * ( v1 + v2 );
01250         get_domain()->element_normal_at( elementHandlesArray[elem_index], normal );
01251     }
01252     else
01253     {
01254         MSQ_SETERR( err )( "No domain constraint set.", MsqError::INVALID_STATE );
01255         return;
01256     }
01257 }
01258 
01259 void PatchData::get_domain_normals_at_corners( size_t elem_index, Vector3D normals_out[], MsqError& err )
01260 {
01261     if( !domain_set() )
01262     {
01263         MSQ_SETERR( err )( "No domain constraint set.", MsqError::INVALID_STATE );
01264         return;
01265     }
01266 
01267     assert( elem_index < elementArray.size() );
01268     if( 2 != TopologyInfo::dimension( elementArray[elem_index].get_element_type() ) )
01269     {
01270         MSQ_SETERR( err )( "Attempt to get corners of non-surface element", MsqError::INVALID_ARG );
01271         return;
01272     }
01273 
01274     if( normalData.empty() )
01275     {
01276         update_cached_normals( err );MSQ_ERRRTN( err );
01277     }
01278 
01279     MsqMeshEntity& elem                = elementArray[elem_index];
01280     const unsigned count               = elem.vertex_count();
01281     const size_t* const vertex_indices = elem.get_vertex_index_array();
01282     for( size_t i = 0; i < count; ++i )
01283     {
01284         const size_t v = vertex_indices[i];
01285         if( vertexNormalIndices.empty() ) { normals_out[i] = normalData[v]; }
01286         else if( vertexNormalIndices[v] < normalData.size() )
01287         {
01288             normals_out[i] = normalData[vertexNormalIndices[v]];
01289         }
01290         else
01291         {
01292             normals_out[i] = vertexArray[v];
01293             get_domain()->element_normal_at( elementHandlesArray[elem_index], normals_out[i] );
01294         }
01295     }
01296 }
01297 
01298 void PatchData::get_domain_normal_at_vertex( size_t vert_index, Mesh::EntityHandle handle, Vector3D& normal,
01299                                              MsqError& err )
01300 {
01301     if( !domain_set() )
01302     {
01303         MSQ_SETERR( err )( "No domain constraint set.", MsqError::INVALID_STATE );
01304         return;
01305     }
01306 
01307     if( normalData.empty() )
01308     {
01309         update_cached_normals( err );MSQ_ERRRTN( err );
01310     }
01311 
01312     if( vertexNormalIndices.empty() ) { normal = normalData[vert_index]; }
01313     else if( vertexNormalIndices[vert_index] < normalData.size() )
01314     {
01315         normal = normalData[vertexNormalIndices[vert_index]];
01316     }
01317     else
01318     {
01319         normal = vertexArray[vert_index];
01320         get_domain()->element_normal_at( handle, normal );
01321     }
01322 }
01323 
01324 void PatchData::get_domain_normal_at_corner( size_t elem_index, unsigned corner, Vector3D& normal, MsqError& err )
01325 {
01326     assert( elem_index < elementArray.size() );
01327     if( 2 != TopologyInfo::dimension( elementArray[elem_index].get_element_type() ) )
01328     {
01329         MSQ_SETERR( err )( "Attempt to get corners of non-surface element", MsqError::INVALID_ARG );
01330         return;
01331     }
01332 
01333     MsqMeshEntity& elem                = elementArray[elem_index];
01334     const size_t* const vertex_indices = elem.get_vertex_index_array();
01335     get_domain_normal_at_vertex( vertex_indices[corner], elementHandlesArray[elem_index], normal, err );MSQ_CHKERR( err );
01336 }
01337 
01338 void PatchData::set_mesh( Mesh* ms )
01339 {
01340     myMesh = ms;
01341     // observers should treat this the same as if the
01342     // instance of this object wzs being deleted.
01343     notify_patch_destroyed();
01344 }
01345 
01346 void PatchData::set_domain( MeshDomain* d )
01347 {
01348     myDomain = d;
01349 
01350     // clear all cached data from the previous domain
01351     vertexNormalIndices.clear();
01352     normalData.clear();
01353     // vertexDomainDOF.clear();
01354 
01355     // observers should treat this the same as if the
01356     // instance of this object wzs being deleted.
01357     notify_patch_destroyed();
01358 }
01359 
01360 static int width( double d )
01361 {
01362     if( d == 0.0 ) return 1;
01363     const int max_precision = 6;
01364     int w                   = (int)ceil( log10( 0.001 + fabs( d ) ) );
01365     if( w < 0 ) w = 2 + std::min( max_precision, -w );
01366     if( d < 0.0 ) ++w;
01367     return w;
01368 }
01369 static int width( size_t t )
01370 {
01371     return t ? (int)ceil( log10( (double)( 1 + t ) ) ) : 1;
01372 }
01373 static int width( const void* ptr )
01374 {
01375     return width( (size_t)ptr );
01376 }
01377 
01378 ostream& operator<<( ostream& stream, const PatchData& pd )
01379 {
01380     size_t i;
01381     int fw = 5;  // width of bit flags
01382     int hw = 6;  // width of a handle
01383     int cw = 4;  // with of a coordinate value
01384     int iw = 3;  // with of an index
01385     int tw = 3;  // with of the string name of an element type
01386     int xw = cw, yw = cw, zw = cw;
01387 
01388     for( i = 0; i < pd.num_nodes(); ++i )
01389     {
01390         int w = 2 + width( pd.vertexHandlesArray[i] );
01391         if( w > hw ) hw = w;
01392         w = width( pd.vertexArray[i].x() );
01393         if( w > xw ) xw = w;
01394         w = width( pd.vertexArray[i].y() );
01395         if( w > yw ) yw = w;
01396         w = width( pd.vertexArray[i].z() );
01397         if( w > zw ) zw = w;
01398     }
01399     for( i = 0; i < pd.num_elements(); ++i )
01400     {
01401         int w = 2 + width( pd.elementHandlesArray[i] );
01402         if( w > hw ) hw = w;
01403         const char* name = TopologyInfo::short_name( pd.elementArray[i].get_element_type() );
01404         if( name && (int)strlen( name ) > tw ) tw = strlen( name );
01405     }
01406     if( iw < (int)ceil( log10( (double)( 1 + pd.num_nodes() ) ) ) )
01407         iw = (int)ceil( log10( (double)( 1 + pd.num_nodes() ) ) );
01408     if( iw < (int)ceil( log10( (double)( 1 + pd.num_elements() ) ) ) )
01409         iw = (int)ceil( log10( (double)( 1 + pd.num_elements() ) ) );
01410 
01411     stream << "Vertices: " << endl;
01412     stream << "Flags: C: culled, F: fixed, S: slave, P: patch vertex, M: marked" << endl;
01413     stream << setw( iw ) << "Idx"
01414            << " " << setw( hw ) << "Handle"
01415            << " " << setw( cw ) << "X"
01416            << "," << setw( cw ) << "Y"
01417            << "," << setw( cw ) << "Z"
01418            << " " << setw( fw ) << "Flags"
01419            << " "
01420            << "Adj.Elems" << endl
01421            << setw( iw ) << setfill( '-' ) << ""
01422            << " " << setw( hw ) << setfill( '-' ) << ""
01423            << " " << setw( cw ) << setfill( '-' ) << ""
01424            << "," << setw( cw ) << setfill( '-' ) << ""
01425            << "," << setw( cw ) << setfill( '-' ) << ""
01426            << " " << setw( fw ) << setfill( '-' ) << ""
01427            << " " << setfill( ' ' ) << "-------------" << std::endl;
01428     for( i = 0; i < pd.num_nodes(); ++i )
01429     {
01430         stream << setw( iw ) << i << " " << setw( hw ) << pd.vertexHandlesArray[i] << " " << setw( cw )
01431                << pd.vertexArray[i].x() << "," << setw( cw ) << pd.vertexArray[i].y() << "," << setw( cw )
01432                << pd.vertexArray[i].z() << " ";
01433         if( pd.vertexArray[i].is_flag_set( MsqVertex::MSQ_CULLED ) )
01434             stream << "C";
01435         else
01436             stream << " ";
01437         if( pd.vertexArray[i].is_flag_set( MsqVertex::MSQ_HARD_FIXED ) )
01438             stream << "F";
01439         else
01440             stream << " ";
01441         if( pd.vertexArray[i].is_flag_set( MsqVertex::MSQ_DEPENDENT ) )
01442             stream << "S";
01443         else
01444             stream << " ";
01445         if( pd.vertexArray[i].is_flag_set( MsqVertex::MSQ_PATCH_FIXED ) )
01446             stream << "f";
01447         else
01448             stream << " ";
01449         if( pd.vertexArray[i].is_flag_set( MsqVertex::MSQ_MARK ) )
01450             stream << "M";
01451         else
01452             stream << " ";
01453 
01454         if( pd.vertAdjacencyArray.size() )
01455         {
01456             size_t j   = pd.vertAdjacencyOffsets[i];
01457             size_t end = pd.vertAdjacencyOffsets[i + 1];
01458             if( j != end ) stream << " " << pd.vertAdjacencyArray[j++];
01459             for( ; j < end; ++j )
01460                 stream << "," << pd.vertAdjacencyArray[j];
01461         }
01462 
01463         stream << endl;
01464     }
01465 
01466     stream << "Elements: " << endl;
01467     stream << setw( iw ) << "Idx"
01468            << " " << setw( hw ) << "Handle"
01469            << " " << setw( tw + 2 ) << "Type"
01470            << " "
01471            << "Connectivity" << std::endl
01472            << setw( iw ) << setfill( '-' ) << ""
01473            << " " << setw( hw ) << setfill( '-' ) << ""
01474            << " " << setw( tw + 2 ) << setfill( '-' ) << ""
01475            << " " << setfill( ' ' ) << "--------------------------" << std::endl;
01476     for( i = 0; i < pd.num_elements(); ++i )
01477     {
01478         EntityTopology type = pd.elementArray[i].get_element_type();
01479         stream << setw( iw ) << i << " " << setw( hw ) << pd.elementHandlesArray[i] << " " << setw( tw )
01480                << TopologyInfo::short_name( type ) << left << setw( 2 ) << pd.elementArray[i].node_count() << internal
01481                << " " << setw( iw ) << pd.elementArray[i].get_vertex_index_array()[0];
01482         for( size_t j = 1; j < pd.elementArray[i].node_count(); ++j )
01483             stream << "," << setw( iw ) << pd.elementArray[i].get_vertex_index_array()[j];
01484         stream << endl;
01485     }
01486     stream << endl;
01487 
01488     stream << "Mesh: " << ( pd.myMesh ? "yes" : "no" ) << endl;
01489     stream << "Domain: " << ( pd.myDomain ? "yes" : "no" ) << endl;
01490     //   stream << "mType: " << (pd.mType==PatchData::VERTICES_ON_VERTEX_PATCH?"vert-on-vert":
01491     //                           pd.mType==PatchData::ELEMENTS_ON_VERTEX_PATCH?"elem-on-vert":
01492     //                           pd.mType==PatchData::GLOBAL_PATCH?"global":"unknown") << endl;
01493 
01494     if( pd.haveComputedInfos )
01495     {
01496         stream << "ComputedInfos:" << endl;
01497         if( pd.have_computed_info( PatchData::MIN_UNSIGNED_AREA ) )
01498             stream << "\t MIN_UNSINGED_AREA = " << pd.computedInfos[PatchData::MIN_UNSIGNED_AREA] << endl;
01499         if( pd.have_computed_info( PatchData::MAX_UNSIGNED_AREA ) )
01500             stream << "\t MAX_UNSIGNED_AREA = " << pd.computedInfos[PatchData::MAX_UNSIGNED_AREA] << endl;
01501         if( pd.have_computed_info( PatchData::MIN_EDGE_LENGTH ) )
01502             stream << "\t MIN_EDGE_LENGTH = " << pd.computedInfos[PatchData::MIN_EDGE_LENGTH] << endl;
01503         if( pd.have_computed_info( PatchData::MAX_EDGE_LENGTH ) )
01504             stream << "\t MAX_EDGE_LENGTH = " << pd.computedInfos[PatchData::MAX_EDGE_LENGTH] << endl;
01505         if( pd.have_computed_info( PatchData::MINMAX_SIGNED_DET2D ) )
01506             stream << "\t MINMAX_SIGNED_DET2D = " << pd.computedInfos[PatchData::MINMAX_SIGNED_DET2D] << endl;
01507         if( pd.have_computed_info( PatchData::MINMAX_SIGNED_DET3D ) )
01508             stream << "\t MINMAX_SIGNED_DET3D = " << pd.computedInfos[PatchData::MINMAX_SIGNED_DET3D] << endl;
01509         if( pd.have_computed_info( PatchData::AVERAGE_DET3D ) )
01510             stream << "\t AVERAGE_DET3D = " << pd.computedInfos[PatchData::AVERAGE_DET3D] << endl;
01511     }
01512 
01513     return stream << endl;
01514 }
01515 
01516 void print_patch_data( const PatchData& pd )
01517 {
01518     std::cout << pd << std::endl;
01519 }
01520 
01521 void PatchData::enslave_higher_order_nodes( const size_t* elem_offset_array, unsigned char* vertex_flags,
01522                                             MsqError& ) const
01523 {
01524     for( size_t i = 0; i < elementArray.size(); ++i )
01525     {
01526         size_t start    = elem_offset_array[i];
01527         size_t conn_len = elem_offset_array[i + 1] - start;
01528         for( size_t j = elementArray[i].vertex_count(); j < conn_len; ++j )
01529         {
01530             const size_t vert_idx = elemConnectivityArray[start + j];
01531             assert( vert_idx < vertexHandlesArray.size() );
01532             if( !( vertex_flags[vert_idx] & MsqVertex::MSQ_HARD_FIXED ) )
01533                 vertex_flags[vert_idx] |= MsqVertex::MSQ_DEPENDENT;
01534         }
01535     }
01536 }
01537 
01538 void PatchData::initialize_data( size_t* elem_offset_array, unsigned char* vertex_flags, MsqError& )
01539 {
01540     // Clear out data specific to patch
01541     vertexNormalIndices.clear();
01542     normalData.clear();
01543     // vertexDomainDOF.clear();
01544 
01545     // Clear any vertex->element adjacency data.  It
01546     // is probably invalid, and certainly will be by the time
01547     // this function completes if the mesh contains higher-order
01548     // elements.
01549     vertAdjacencyArray.clear();
01550     vertAdjacencyOffsets.clear();
01551     size_t i, j;
01552     for( i = 0; i < elementArray.size(); ++i )
01553     {
01554         size_t start    = elem_offset_array[i];
01555         size_t conn_len = elem_offset_array[i + 1] - start;
01556         assert( conn_len > 0 );
01557         elementArray[i].set_connectivity( &elemConnectivityArray[start], conn_len );
01558     }
01559 
01560     // Use vertAdjacencyOffsets array as temporary storage.
01561     vertAdjacencyOffsets.resize( vertexHandlesArray.size() + 1 );
01562     size_t* vertex_index_map = arrptr( vertAdjacencyOffsets );
01563 
01564     // Count number of free vertices and initialize vertex_index_map
01565     numFreeVertices = 0;
01566     for( i = 0; i < vertexHandlesArray.size(); ++i )
01567     {
01568         if( !( vertex_flags[i] & MsqVertex::MSQ_FIXED ) ) ++numFreeVertices;
01569         vertex_index_map[i] = i;
01570     }
01571 
01572     // Re-order vertices such that all free vertices are
01573     // first in the list.  Construct map from old to new
01574     // position in list for updating element connectivity.
01575     i = 0;
01576     j = numFreeVertices;
01577     for( ;; ++i, ++j )
01578     {
01579         // find next fixed vertex in the range [i,numFreeVertices)
01580         for( ; i < numFreeVertices && !( vertex_flags[i] & MsqVertex::MSQ_FIXED ); ++i )
01581             ;
01582         // if no more fixed vertices in the free vertex range [0,numFreeVertices)
01583         if( i == numFreeVertices ) break;
01584         // find the next free vertex in the range [j,num_nodes)
01585         for( ; vertex_flags[j] & MsqVertex::MSQ_FIXED; ++j )
01586             ;
01587         // swap fixed (i) and free (j) vertices
01588         vertex_index_map[i] = j;
01589         vertex_index_map[j] = i;
01590         std::swap( vertexHandlesArray[i], vertexHandlesArray[j] );
01591         std::swap( vertex_flags[i], vertex_flags[j] );
01592     }
01593     assert( i == numFreeVertices );
01594     assert( j <= vertexHandlesArray.size() );
01595 
01596     // Update element connectivity data for new vertex indices
01597     for( i = 0; i < elemConnectivityArray.size(); ++i )
01598         elemConnectivityArray[i] = vertex_index_map[elemConnectivityArray[i]];
01599 
01600     // Reorder vertices such that free, slave vertices
01601     // occur after free, non-slave vertices in list.
01602     numSlaveVertices = 0;
01603     for( i = 0; i < vertexHandlesArray.size(); ++i )
01604     {
01605         if( ( vertex_flags[i] & MsqVertex::MSQ_DEPENDENT ) && !( vertex_flags[i] & MsqVertex::MSQ_FIXED ) )
01606             ++numSlaveVertices;
01607     }
01608     numFreeVertices -= numSlaveVertices;
01609 
01610     if( numSlaveVertices )
01611     {
01612         // re-initialize vertex index map
01613         for( i = 0; i < vertexHandlesArray.size(); ++i )
01614             vertex_index_map[i] = i;
01615 
01616         // Re-order free vertices such that all slave vertices are
01617         // last in the list.  Construct map from old to new
01618         // position in list for updating element connectivity.
01619         i = 0;
01620         j = numFreeVertices;
01621         for( ;; ++i, ++j )
01622         {
01623             // find next slave vertex in the range [i,numFreeVertices)
01624             for( ; i < numFreeVertices && !( vertex_flags[i] & MsqVertex::MSQ_DEPENDENT ); ++i )
01625                 ;
01626             // if no more slave vertices in [0,numFreeVertices), then done.
01627             if( i == numFreeVertices ) break;
01628             // find the next free (non-slave) vertex in the range
01629             //   [numFreeVertices,numFreeVertices+numSlaveVertices)
01630             for( ; vertex_flags[j] & MsqVertex::MSQ_DEPENDENT; ++j )
01631                 ;
01632             // swap free (j) and slave (i) vertices
01633             vertex_index_map[i] = j;
01634             vertex_index_map[j] = i;
01635             std::swap( vertexHandlesArray[i], vertexHandlesArray[j] );
01636             std::swap( vertex_flags[i], vertex_flags[j] );
01637         }
01638         assert( i == numFreeVertices );
01639         assert( j <= numFreeVertices + numSlaveVertices );
01640 
01641         // Update element connectivity data for new vertex indices
01642         for( i = 0; i < elemConnectivityArray.size(); ++i )
01643             elemConnectivityArray[i] = vertex_index_map[elemConnectivityArray[i]];
01644     }
01645 
01646     // Clear temporary data
01647     vertAdjacencyOffsets.clear();
01648 
01649     notify_new_patch();
01650 }
01651 
01652 size_t PatchData::num_corners() const
01653 {
01654     size_t result = 0;
01655     for( unsigned i = 0; i < elementArray.size(); ++i )
01656         result += elementArray[i].vertex_count();
01657     return result;
01658 }
01659 
01660 void PatchData::fill( size_t num_vertex, const double* coords, size_t num_elem, EntityTopology type,
01661                       const size_t* connectivity, const bool* fixed, MsqError& err )
01662 {
01663     std::vector< EntityTopology > types( num_elem );
01664     std::fill( types.begin(), types.end(), type );
01665     const EntityTopology* type_ptr = num_elem ? arrptr( types ) : 0;
01666     this->fill( num_vertex, coords, num_elem, type_ptr, connectivity, fixed, err );MSQ_CHKERR( err );
01667 }
01668 
01669 void PatchData::fill( size_t num_vertex, const double* coords, size_t num_elem, const EntityTopology* types,
01670                       const size_t* conn, const bool* fixed, MsqError& err )
01671 {
01672     std::vector< size_t > lengths( num_elem );
01673     std::transform( types, types + num_elem, lengths.begin(), std::ptr_fun( TopologyInfo::corners ) );
01674     const size_t* len_ptr = num_elem ? arrptr( lengths ) : 0;
01675     this->fill( num_vertex, coords, num_elem, types, len_ptr, conn, fixed, err );MSQ_CHKERR( err );
01676 }
01677 
01678 void PatchData::fill( size_t num_vertex, const double* coords, size_t num_elem, const EntityTopology* types,
01679                       const size_t* lengths, const size_t* conn, const bool* fixed, MsqError& err )
01680 {
01681     size_t i;
01682 
01683     // count vertex uses
01684     size_t num_uses = std::accumulate( lengths, lengths + num_elem, 0 );
01685 
01686     // Allocate storage for data
01687     vertexArray.resize( num_vertex );
01688     vertexHandlesArray.resize( num_vertex );
01689     elementArray.resize( num_elem );
01690     elementHandlesArray.resize( num_elem );
01691     elemConnectivityArray.resize( num_uses );
01692     numFreeVertices  = 0;
01693     numSlaveVertices = 0;
01694 
01695     // Must call clear() first so that any stale values get
01696     // zero'd when we call resize.
01697     byteArray.clear();
01698     if( fixed )
01699     {
01700         byteArray.resize( num_vertex, 0 );
01701         for( i = 0; i < num_vertex; ++i )
01702             if( fixed[i] ) byteArray[i] |= ( MsqVertex::MSQ_HARD_FIXED | MsqVertex::MSQ_PATCH_FIXED );
01703     }
01704 
01705     for( i = 0; i < num_elem; ++i )
01706     {
01707         element_by_index( i ).set_element_type( types[i] );
01708         elementHandlesArray[i] = (Mesh::ElementHandle)i;
01709     }
01710     for( i = 0; i < num_vertex; ++i )
01711         vertexHandlesArray[i] = (Mesh::VertexHandle)i;
01712 
01713     memcpy( get_connectivity_array(), conn, num_uses * sizeof( size_t ) );
01714 
01715     std::vector< size_t > offsets( num_elem + 1 );
01716     size_t sum = offsets[0] = 0;
01717     for( i = 1; i <= num_elem; ++i )
01718         offsets[i] = sum += lengths[i - 1];
01719 
01720     const Settings::HigherOrderSlaveMode ho_mode =
01721         mSettings ? mSettings->get_slaved_ho_node_mode() : Settings::SLAVE_ALL;
01722     switch( ho_mode )
01723     {
01724         case Settings::SLAVE_ALL:
01725             byteArray.resize( num_vertex, 0 );
01726             enslave_higher_order_nodes( arrptr( offsets ), arrptr( byteArray ), err );MSQ_ERRRTN( err );
01727             break;
01728         case Settings::SLAVE_NONE:
01729             // Do nothing.  We clear other bits when processing the 'fixed' array above.
01730             break;
01731         default:
01732             MSQ_SETERR( err )
01733             ( "Specified higher-order noded slaving scheme not supported "
01734               "when initializind PatchData using PatchData::fill",
01735               MsqError::NOT_IMPLEMENTED );
01736             return;
01737     }
01738 
01739     this->initialize_data( arrptr( offsets ), arrptr( byteArray ), err );MSQ_ERRRTN( err );
01740 
01741     // initialize_data will re-order vertex handles and
01742     // update element connectivity accordingly.  Use
01743     // the values we stored in vertexHandlesArray to
01744     // figure out the new index of each vertex, and initialize
01745     // the vertex.
01746     for( i = 0; i < num_vertex; ++i )
01747         vertexArray[i] = coords + 3 * (size_t)vertexHandlesArray[i];
01748 
01749     for( i = 0; i < num_vertex; ++i )
01750         vertexArray[i].flags() = byteArray[i];
01751 }
01752 
01753 void PatchData::make_handles_unique( Mesh::EntityHandle* handles, size_t& count, size_t* index_map )
01754 {
01755     if( count < 2 ) { return; }
01756     // save this now, as we'll be changing count later
01757     const size_t* index_end = index_map + count;
01758 
01759     if( index_map )
01760     {
01761         // Copy input handle list into index map as a temporary
01762         // copy of the input list.
01763         assert( sizeof( Mesh::EntityHandle ) == sizeof( size_t ) );
01764         memcpy( index_map, handles, count * sizeof( size_t ) );
01765     }
01766 
01767     // Make handles a unique list
01768     std::sort( handles, handles + count );
01769     Mesh::EntityHandle* end = std::unique( handles, handles + count );
01770     count                   = end - handles;
01771 
01772     if( index_map )
01773     {
01774         // Replace each handle in index_map with the index of
01775         // its position in the handles array
01776         Mesh::EntityHandle* pos;
01777         for( size_t* iter = index_map; iter != index_end; ++iter )
01778         {
01779             pos   = std::lower_bound( handles, handles + count, (Mesh::EntityHandle)*iter );
01780             *iter = pos - handles;
01781         }
01782     }
01783 }
01784 
01785 void PatchData::fill_global_patch( MsqError& err )
01786 {
01787     GlobalPatch gp;
01788     gp.set_mesh( get_mesh() );
01789     PatchIterator iter( &gp );
01790     bool b = iter.get_next_patch( *this, err );MSQ_ERRRTN( err );
01791     if( !b ) MSQ_SETERR( err )( "Empty patch in PatchData::fill_global_patch", MsqError::INVALID_MESH );
01792     assert( b );
01793 }
01794 
01795 void PatchData::set_mesh_entities( std::vector< Mesh::ElementHandle >& elements,
01796                                    std::vector< Mesh::VertexHandle >& free_vertices, MsqError& err )
01797 {
01798     Mesh* current_mesh = get_mesh();
01799     if( !current_mesh )
01800     {
01801         MSQ_SETERR( err )( "No Mesh associated with PatchData.", MsqError::INVALID_STATE );
01802         return;
01803     }
01804 
01805     if( elements.empty() )
01806     {
01807         clear();
01808         return;
01809     }
01810 
01811     elementHandlesArray = elements;
01812     get_mesh()->elements_get_attached_vertices( arrptr( elementHandlesArray ), elementHandlesArray.size(),
01813                                                 vertexHandlesArray, offsetArray, err );MSQ_ERRRTN( err );
01814 
01815     // Construct CSR-rep element connectivity data
01816     size_t num_vert = vertexHandlesArray.size();
01817     elemConnectivityArray.resize( num_vert );
01818     make_handles_unique( arrptr( vertexHandlesArray ), num_vert, arrptr( elemConnectivityArray ) );
01819     vertexHandlesArray.resize( num_vert );
01820 
01821     // Get element topologies
01822     std::vector< EntityTopology > elem_topologies( elementHandlesArray.size() );
01823     get_mesh()->elements_get_topologies( arrptr( elementHandlesArray ), arrptr( elem_topologies ),
01824                                          elementHandlesArray.size(), err );MSQ_ERRRTN( err );
01825 
01826     // get vertex bits from mesh
01827     byteArray.resize( vertexHandlesArray.size() );
01828     get_mesh()->vertices_get_byte( arrptr( vertexHandlesArray ), arrptr( byteArray ), vertexHandlesArray.size(), err );MSQ_ERRRTN( err );
01829 
01830     // if free_vertices is not empty, then we need to mark as
01831     // fixed any vertices *not* in that list.
01832     if( free_vertices.size() == 1 )
01833     {
01834         for( size_t i = 0; i < vertexHandlesArray.size(); ++i )
01835             if( vertexHandlesArray[i] == free_vertices.front() )
01836                 byteArray[i] &= ~MsqVertex::MSQ_PATCH_FIXED;
01837             else
01838                 byteArray[i] |= MsqVertex::MSQ_PATCH_FIXED;
01839     }
01840     else if( !free_vertices.empty() )
01841     {
01842         // sort and remove duplicates from free_vertices list.
01843         std::sort( free_vertices.begin(), free_vertices.end() );
01844         free_vertices.erase( std::unique( free_vertices.begin(), free_vertices.end() ), free_vertices.end() );
01845 
01846         for( size_t i = 0; i < vertexHandlesArray.size(); ++i )
01847         {
01848             if( ( byteArray[i] & MsqVertex::MSQ_DEPENDENT ) ||
01849                 std::binary_search( free_vertices.begin(), free_vertices.end(), vertexHandlesArray[i] ) )
01850                 byteArray[i] &= ~MsqVertex::MSQ_PATCH_FIXED;
01851             else
01852                 byteArray[i] |= MsqVertex::MSQ_PATCH_FIXED;
01853         }
01854     }
01855 
01856     // set element types
01857     elementArray.resize( elementHandlesArray.size() );
01858     for( size_t i = 0; i < elementHandlesArray.size(); ++i )
01859         elementArray[i].set_element_type( elem_topologies[i] );
01860 
01861     // get element connectivity, group vertices by free/slave/fixed state
01862     initialize_data( arrptr( offsetArray ), arrptr( byteArray ), err );MSQ_ERRRTN( err );
01863 
01864     // get vertex coordinates
01865     vertexArray.resize( vertexHandlesArray.size() );
01866     get_mesh()->vertices_get_coordinates( arrptr( vertexHandlesArray ), arrptr( vertexArray ),
01867                                           vertexHandlesArray.size(), err );MSQ_ERRRTN( err );
01868 
01869     // set vertex flags
01870     for( size_t i = 0; i < vertexArray.size(); ++i )
01871         vertexArray[i].flags() = byteArray[i];
01872 }
01873 
01874 void PatchData::get_sample_location( size_t element_index, Sample sample, Vector3D& result, MsqError& err ) const
01875 {
01876     const MsqMeshEntity& elem      = element_by_index( element_index );
01877     const NodeSet ho_bits          = non_slave_node_set( element_index );
01878     const MappingFunction* const f = get_mapping_function( elem.get_element_type() );
01879     if( !f )
01880     {
01881         MSQ_SETERR( err )( "No mapping function", MsqError::UNSUPPORTED_ELEMENT );
01882         return;
01883     }
01884 
01885     double coeff[27];
01886     size_t num_coeff, index[27];
01887     f->coefficients( sample, ho_bits, coeff, index, num_coeff, err );MSQ_ERRRTN( err );
01888     f->convert_connectivity_indices( elem.node_count(), index, num_coeff, err );MSQ_ERRRTN( err );
01889 
01890     const size_t* const conn = elem.get_vertex_index_array();
01891     assert( num_coeff > 0 );
01892     result = coeff[0] * vertex_by_index( conn[index[0]] );
01893     for( unsigned i = 1; i < num_coeff; ++i )
01894         result += coeff[i] * vertex_by_index( conn[index[i]] );
01895 }
01896 
01897 NodeSet PatchData::non_slave_node_set( size_t element_index ) const
01898 {
01899     const MsqMeshEntity& elem = element_by_index( element_index );
01900     EntityTopology type       = elem.get_element_type();
01901     const size_t* conn        = elem.get_vertex_index_array();
01902     const size_t n            = elem.node_count();
01903 
01904     MsqError err;
01905     bool have_midedge, have_midface, have_midelem;
01906     unsigned num_edge = 0, num_face = 0, num_corner = TopologyInfo::corners( type );
01907     TopologyInfo::higher_order( type, n, have_midedge, have_midface, have_midelem, err );
01908     num_edge = TopologyInfo::edges( type );
01909     if( TopologyInfo::dimension( type ) == 2 )
01910         num_face = 1;
01911     else
01912         num_face = TopologyInfo::faces( type );
01913 
01914     NodeSet result;
01915     result.set_all_corner_nodes( type );
01916     if( have_midedge )
01917     {
01918         for( unsigned i = 0; i < num_edge; ++i )
01919         {
01920             if( !( vertex_by_index( conn[num_corner + i] ).get_flags() & MsqVertex::MSQ_DEPENDENT ) )
01921                 result.set_mid_edge_node( i );
01922         }
01923     }
01924     if( have_midface )
01925     {
01926         for( unsigned i = 0; i < num_face; ++i )
01927         {
01928             if( !( vertex_by_index( conn[num_corner + num_edge + i] ).get_flags() & MsqVertex::MSQ_DEPENDENT ) )
01929                 result.set_mid_face_node( i );
01930         }
01931     }
01932     if( have_midelem &&
01933         !( vertex_by_index( conn[num_corner + num_edge + num_face] ).get_flags() & MsqVertex::MSQ_DEPENDENT ) )
01934         result.set_mid_region_node();
01935 
01936     return result;
01937 }
01938 
01939 void PatchData::get_samples( size_t element, std::vector< Sample >& samples, MsqError& ) const
01940 {
01941     NodeSet ns = get_samples( element );
01942     samples.resize( ns.num_nodes() );
01943     std::vector< Sample >::iterator i = samples.begin();
01944 
01945     unsigned j;
01946     EntityTopology type = element_by_index( element ).get_element_type();
01947     for( j = 0; j < TopologyInfo::corners( type ); ++j )
01948         if( ns.corner_node( j ) ) *( i++ ) = Sample( 0, j );
01949     for( j = 0; j < TopologyInfo::edges( type ); ++j )
01950         if( ns.mid_edge_node( j ) ) *( i++ ) = Sample( 1, j );
01951     if( TopologyInfo::dimension( type ) == 3 )
01952     {
01953         for( j = 0; j < TopologyInfo::faces( type ); ++j )
01954             if( ns.mid_face_node( j ) ) *( i++ ) = Sample( 2, j );
01955         if( ns.mid_region_node() ) *( i++ ) = Sample( 3, 0 );
01956     }
01957     else if( ns.mid_face_node( 0 ) )
01958         *( i++ ) = Sample( 2, 0 );
01959 
01960     assert( i == samples.end() );
01961 }
01962 
01963 bool PatchData::attach_extra_data( ExtraData* data )
01964 {
01965     if( data->patchNext ) { return false; }
01966 
01967     if( !data->patchPtr )
01968         data->patchPtr = this;
01969     else if( data->patchPtr != this )
01970         return false;
01971 
01972     data->patchNext = dataList;
01973     dataList        = data;
01974     return true;
01975 }
01976 
01977 bool PatchData::remove_extra_data( ExtraData* data )
01978 {
01979     if( data->patchPtr != this ) return false;
01980 
01981     if( dataList == data )
01982     {
01983         dataList        = data->patchNext;
01984         data->patchNext = 0;
01985         data->patchPtr  = 0;
01986         return true;
01987     }
01988 
01989     for( ExtraData* iter = dataList; iter; iter = iter->patchNext )
01990         if( iter->patchNext == data )
01991         {
01992             iter->patchNext = data->patchNext;
01993             data->patchNext = 0;
01994             data->patchPtr  = 0;
01995             return true;
01996         }
01997 
01998     return false;
01999 }
02000 
02001 void PatchData::notify_new_patch()
02002 {
02003     for( ExtraData* iter = dataList; iter; iter = iter->patchNext )
02004         iter->notify_new_patch();
02005 }
02006 
02007 void PatchData::notify_sub_patch( PatchData& sub_patch, const size_t* vertex_map, const size_t* element_map,
02008                                   MsqError& err )
02009 {
02010     for( ExtraData* iter = dataList; iter; iter = iter->patchNext )
02011     {
02012         iter->notify_sub_patch( sub_patch, vertex_map, element_map, err );MSQ_ERRRTN( err );
02013     }
02014 }
02015 
02016 void PatchData::notify_patch_destroyed()
02017 {
02018     // Remove all ExtraDatas from list and notify them
02019     // that they are being destroyed.
02020     while( dataList )
02021     {
02022         ExtraData* dead_data = dataList;
02023         dataList             = dataList->patchNext;
02024         dead_data->patchNext = 0;
02025         dead_data->patchPtr  = 0;
02026         dead_data->notify_patch_destroyed();
02027     }
02028 }
02029 
02030 }  // namespace MBMesquite
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines