MOAB: Mesh Oriented datABase  (version 5.2.1)
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         { get_domain()->snap_to( vertexHandlesArray[vertex_index], vertexArray[vertex_index] ); }
01132         // entire domain is 2-D (all vertices have a single normal)
01133         else if( vertexNormalIndices.empty() )
01134         {
01135             get_domain()->closest_point( vertexHandlesArray[vertex_index], Vector3D( vertexArray[vertex_index] ),
01136                                          vertexArray[vertex_index], normalData[vertex_index], err );MSQ_ERRRTN( err );
01137         }
01138         else if( vertexNormalIndices[vertex_index] < normalData.size() )
01139         {  // vertex has a unique normal
01140             get_domain()->closest_point( vertexHandlesArray[vertex_index], Vector3D( vertexArray[vertex_index] ),
01141                                          vertexArray[vertex_index], normalData[vertexNormalIndices[vertex_index]],
01142                                          err );MSQ_ERRRTN( err );
01143         }
01144         else
01145         {  // vertex has multiple normals
01146             get_domain()->snap_to( vertexHandlesArray[vertex_index], vertexArray[vertex_index] );
01147         }
01148     }
01149 }
01150 
01151 void PatchData::update_cached_normals( MsqError& err )
01152 {
01153     size_t i;
01154 
01155     MeshDomain* domain = get_domain();
01156     if( !domain )
01157     {
01158         MSQ_SETERR( err )( "No domain constraint set.", MsqError::INVALID_STATE );
01159         return;
01160     }
01161 
01162     // Determine which vertices lie on surfaces
01163     vertexDomainDOF.resize( num_nodes() );
01164     domain->domain_DoF( arrptr( vertexHandlesArray ), arrptr( vertexDomainDOF ), num_nodes(), err );MSQ_ERRRTN( err );
01165 
01166     // Count how many vertices have a single normal
01167     // Sun doesn't support partial template specialization, so can't use std::count
01168     // size_t n = std::count( vertexDomainDOF.begin(), vertexDomainDOF.end(), 2 );
01169     size_t n = 0;
01170     std::vector< unsigned short >::iterator k;
01171     for( k = vertexDomainDOF.begin(); k != vertexDomainDOF.end(); ++k )
01172         if( *k == 2 ) ++n;
01173 
01174     normalData.resize( n );
01175 
01176     // If all vertices are on a surface, pass in the existing handles array
01177     // and store a single normal per vertex.
01178     if( n == num_nodes() )
01179     {
01180         std::copy( vertexArray.begin(), vertexArray.end(), normalData.begin() );
01181         domain->vertex_normal_at( arrptr( vertexHandlesArray ), arrptr( normalData ), num_nodes(), err );
01182         vertexNormalIndices.clear();
01183         vertexDomainDOF.clear();MSQ_ERRRTN( err );
01184     }
01185     else
01186     {
01187         vertexNormalIndices.resize( num_nodes() );
01188         size_t nn = 0;
01189         for( i = 0; i < num_nodes(); ++i )
01190         {
01191             if( vertexDomainDOF[i] == 2 )
01192             {
01193                 normalData[nn] = vertexArray[i];
01194                 domain->vertex_normal_at( vertexHandlesArray[i], normalData[nn] );
01195                 vertexNormalIndices[i] = nn;
01196                 ++nn;
01197             }
01198             else
01199             {
01200                 vertexNormalIndices[i] = n + 1;
01201             }
01202         }
01203         assert( nn == n );
01204     }
01205 }
01206 
01207 void PatchData::get_domain_normal_at_element( size_t elem_index, Vector3D& surf_norm, MsqError& err )
01208 {
01209     // check if element as a mid-face node
01210     const MsqMeshEntity& elem = element_by_index( elem_index );
01211     const int mid_node = TopologyInfo::higher_order_from_side( elem.get_element_type(), elem.node_count(), 2, 0, err );MSQ_ERRRTN( err );
01212     // if we have the mid-element vertex, get cached normal for it
01213     if( mid_node > 0 )
01214     {
01215         get_domain_normal_at_vertex( elem.get_vertex_index_array()[mid_node], elementHandlesArray[elem_index],
01216                                      surf_norm, err );MSQ_ERRRTN( err );
01217     }
01218     // otherwise query domain for normal at element centroid
01219     else if( domain_set() )
01220     {
01221         assert( elem_index < elementArray.size() );
01222         elementArray[elem_index].get_centroid( surf_norm, *this, err );MSQ_ERRRTN( err );
01223         get_domain()->element_normal_at( elementHandlesArray[elem_index], surf_norm );
01224     }
01225     else
01226         MSQ_SETERR( err )( "No domain constraint set.", MsqError::INVALID_STATE );
01227 }
01228 
01229 void PatchData::get_domain_normal_at_mid_edge( size_t elem_index, unsigned edge_num, Vector3D& normal, MsqError& err )
01230 {
01231     // check if element as a mid-edge node
01232     const MsqMeshEntity& elem = element_by_index( elem_index );
01233     const int mid_node =
01234         TopologyInfo::higher_order_from_side( elem.get_element_type(), elem.node_count(), 1, edge_num, err );MSQ_ERRRTN( err );
01235     // if we have the mid-edge vertex, get cached normal for it
01236     if( mid_node > 0 )
01237     {
01238         get_domain_normal_at_vertex( elem.get_vertex_index_array()[mid_node], elementHandlesArray[elem_index], normal,
01239                                      err );MSQ_ERRRTN( err );
01240     }
01241     // otherwise query domain for normal at center of edge
01242     else if( domain_set() )
01243     {
01244         const unsigned* edge = TopologyInfo::edge_vertices( elem.get_element_type(), edge_num, err );MSQ_ERRRTN( err );
01245         const MsqVertex& v1 = vertex_by_index( elem.get_vertex_index_array()[edge[0]] );
01246         const MsqVertex& v2 = vertex_by_index( elem.get_vertex_index_array()[edge[1]] );
01247         normal              = 0.5 * ( v1 + v2 );
01248         get_domain()->element_normal_at( elementHandlesArray[elem_index], normal );
01249     }
01250     else
01251     {
01252         MSQ_SETERR( err )( "No domain constraint set.", MsqError::INVALID_STATE );
01253         return;
01254     }
01255 }
01256 
01257 void PatchData::get_domain_normals_at_corners( size_t elem_index, Vector3D normals_out[], MsqError& err )
01258 {
01259     if( !domain_set() )
01260     {
01261         MSQ_SETERR( err )( "No domain constraint set.", MsqError::INVALID_STATE );
01262         return;
01263     }
01264 
01265     assert( elem_index < elementArray.size() );
01266     if( 2 != TopologyInfo::dimension( elementArray[elem_index].get_element_type() ) )
01267     {
01268         MSQ_SETERR( err )( "Attempt to get corners of non-surface element", MsqError::INVALID_ARG );
01269         return;
01270     }
01271 
01272     if( normalData.empty() )
01273     {
01274         update_cached_normals( err );MSQ_ERRRTN( err );
01275     }
01276 
01277     MsqMeshEntity& elem                = elementArray[elem_index];
01278     const unsigned count               = elem.vertex_count();
01279     const size_t* const vertex_indices = elem.get_vertex_index_array();
01280     for( size_t i = 0; i < count; ++i )
01281     {
01282         const size_t v = vertex_indices[i];
01283         if( vertexNormalIndices.empty() ) { normals_out[i] = normalData[v]; }
01284         else if( vertexNormalIndices[v] < normalData.size() )
01285         {
01286             normals_out[i] = normalData[vertexNormalIndices[v]];
01287         }
01288         else
01289         {
01290             normals_out[i] = vertexArray[v];
01291             get_domain()->element_normal_at( elementHandlesArray[elem_index], normals_out[i] );
01292         }
01293     }
01294 }
01295 
01296 void PatchData::get_domain_normal_at_vertex( size_t vert_index, Mesh::EntityHandle handle, Vector3D& normal,
01297                                              MsqError& err )
01298 {
01299     if( !domain_set() )
01300     {
01301         MSQ_SETERR( err )( "No domain constraint set.", MsqError::INVALID_STATE );
01302         return;
01303     }
01304 
01305     if( normalData.empty() )
01306     {
01307         update_cached_normals( err );MSQ_ERRRTN( err );
01308     }
01309 
01310     if( vertexNormalIndices.empty() ) { normal = normalData[vert_index]; }
01311     else if( vertexNormalIndices[vert_index] < normalData.size() )
01312     {
01313         normal = normalData[vertexNormalIndices[vert_index]];
01314     }
01315     else
01316     {
01317         normal = vertexArray[vert_index];
01318         get_domain()->element_normal_at( handle, normal );
01319     }
01320 }
01321 
01322 void PatchData::get_domain_normal_at_corner( size_t elem_index, unsigned corner, Vector3D& normal, MsqError& err )
01323 {
01324     assert( elem_index < elementArray.size() );
01325     if( 2 != TopologyInfo::dimension( elementArray[elem_index].get_element_type() ) )
01326     {
01327         MSQ_SETERR( err )( "Attempt to get corners of non-surface element", MsqError::INVALID_ARG );
01328         return;
01329     }
01330 
01331     MsqMeshEntity& elem                = elementArray[elem_index];
01332     const size_t* const vertex_indices = elem.get_vertex_index_array();
01333     get_domain_normal_at_vertex( vertex_indices[corner], elementHandlesArray[elem_index], normal, err );MSQ_CHKERR( err );
01334 }
01335 
01336 void PatchData::set_mesh( Mesh* ms )
01337 {
01338     myMesh = ms;
01339     // observers should treat this the same as if the
01340     // instance of this object wzs being deleted.
01341     notify_patch_destroyed();
01342 }
01343 
01344 void PatchData::set_domain( MeshDomain* d )
01345 {
01346     myDomain = d;
01347 
01348     // clear all cached data from the previous domain
01349     vertexNormalIndices.clear();
01350     normalData.clear();
01351     // vertexDomainDOF.clear();
01352 
01353     // observers should treat this the same as if the
01354     // instance of this object wzs being deleted.
01355     notify_patch_destroyed();
01356 }
01357 
01358 static int width( double d )
01359 {
01360     if( d == 0.0 ) return 1;
01361     const int max_precision = 6;
01362     int w                   = (int)ceil( log10( 0.001 + fabs( d ) ) );
01363     if( w < 0 ) w = 2 + std::min( max_precision, -w );
01364     if( d < 0.0 ) ++w;
01365     return w;
01366 }
01367 static int width( size_t t )
01368 {
01369     return t ? (int)ceil( log10( (double)( 1 + t ) ) ) : 1;
01370 }
01371 static int width( const void* ptr )
01372 {
01373     return width( (size_t)ptr );
01374 }
01375 
01376 ostream& operator<<( ostream& stream, const PatchData& pd )
01377 {
01378     size_t i;
01379     int fw = 5;  // width of bit flags
01380     int hw = 6;  // width of a handle
01381     int cw = 4;  // with of a coordinate value
01382     int iw = 3;  // with of an index
01383     int tw = 3;  // with of the string name of an element type
01384     int xw = cw, yw = cw, zw = cw;
01385 
01386     for( i = 0; i < pd.num_nodes(); ++i )
01387     {
01388         int w = 2 + width( pd.vertexHandlesArray[i] );
01389         if( w > hw ) hw = w;
01390         w = width( pd.vertexArray[i].x() );
01391         if( w > xw ) xw = w;
01392         w = width( pd.vertexArray[i].y() );
01393         if( w > yw ) yw = w;
01394         w = width( pd.vertexArray[i].z() );
01395         if( w > zw ) zw = w;
01396     }
01397     for( i = 0; i < pd.num_elements(); ++i )
01398     {
01399         int w = 2 + width( pd.elementHandlesArray[i] );
01400         if( w > hw ) hw = w;
01401         const char* name = TopologyInfo::short_name( pd.elementArray[i].get_element_type() );
01402         if( name && (int)strlen( name ) > tw ) tw = strlen( name );
01403     }
01404     if( iw < (int)ceil( log10( (double)( 1 + pd.num_nodes() ) ) ) )
01405         iw = (int)ceil( log10( (double)( 1 + pd.num_nodes() ) ) );
01406     if( iw < (int)ceil( log10( (double)( 1 + pd.num_elements() ) ) ) )
01407         iw = (int)ceil( log10( (double)( 1 + pd.num_elements() ) ) );
01408 
01409     stream << "Vertices: " << endl;
01410     stream << "Flags: C: culled, F: fixed, S: slave, P: patch vertex, M: marked" << endl;
01411     stream << setw( iw ) << "Idx"
01412            << " " << setw( hw ) << "Handle"
01413            << " " << setw( cw ) << "X"
01414            << "," << setw( cw ) << "Y"
01415            << "," << setw( cw ) << "Z"
01416            << " " << setw( fw ) << "Flags"
01417            << " "
01418            << "Adj.Elems" << endl
01419            << setw( iw ) << setfill( '-' ) << ""
01420            << " " << setw( hw ) << setfill( '-' ) << ""
01421            << " " << setw( cw ) << setfill( '-' ) << ""
01422            << "," << setw( cw ) << setfill( '-' ) << ""
01423            << "," << setw( cw ) << setfill( '-' ) << ""
01424            << " " << setw( fw ) << setfill( '-' ) << ""
01425            << " " << setfill( ' ' ) << "-------------" << std::endl;
01426     for( i = 0; i < pd.num_nodes(); ++i )
01427     {
01428         stream << setw( iw ) << i << " " << setw( hw ) << pd.vertexHandlesArray[i] << " " << setw( cw )
01429                << pd.vertexArray[i].x() << "," << setw( cw ) << pd.vertexArray[i].y() << "," << setw( cw )
01430                << pd.vertexArray[i].z() << " ";
01431         if( pd.vertexArray[i].is_flag_set( MsqVertex::MSQ_CULLED ) )
01432             stream << "C";
01433         else
01434             stream << " ";
01435         if( pd.vertexArray[i].is_flag_set( MsqVertex::MSQ_HARD_FIXED ) )
01436             stream << "F";
01437         else
01438             stream << " ";
01439         if( pd.vertexArray[i].is_flag_set( MsqVertex::MSQ_DEPENDENT ) )
01440             stream << "S";
01441         else
01442             stream << " ";
01443         if( pd.vertexArray[i].is_flag_set( MsqVertex::MSQ_PATCH_FIXED ) )
01444             stream << "f";
01445         else
01446             stream << " ";
01447         if( pd.vertexArray[i].is_flag_set( MsqVertex::MSQ_MARK ) )
01448             stream << "M";
01449         else
01450             stream << " ";
01451 
01452         if( pd.vertAdjacencyArray.size() )
01453         {
01454             size_t j   = pd.vertAdjacencyOffsets[i];
01455             size_t end = pd.vertAdjacencyOffsets[i + 1];
01456             if( j != end ) stream << " " << pd.vertAdjacencyArray[j++];
01457             for( ; j < end; ++j )
01458                 stream << "," << pd.vertAdjacencyArray[j];
01459         }
01460 
01461         stream << endl;
01462     }
01463 
01464     stream << "Elements: " << endl;
01465     stream << setw( iw ) << "Idx"
01466            << " " << setw( hw ) << "Handle"
01467            << " " << setw( tw + 2 ) << "Type"
01468            << " "
01469            << "Connectivity" << std::endl
01470            << setw( iw ) << setfill( '-' ) << ""
01471            << " " << setw( hw ) << setfill( '-' ) << ""
01472            << " " << setw( tw + 2 ) << setfill( '-' ) << ""
01473            << " " << setfill( ' ' ) << "--------------------------" << std::endl;
01474     for( i = 0; i < pd.num_elements(); ++i )
01475     {
01476         EntityTopology type = pd.elementArray[i].get_element_type();
01477         stream << setw( iw ) << i << " " << setw( hw ) << pd.elementHandlesArray[i] << " " << setw( tw )
01478                << TopologyInfo::short_name( type ) << left << setw( 2 ) << pd.elementArray[i].node_count() << internal
01479                << " " << setw( iw ) << pd.elementArray[i].get_vertex_index_array()[0];
01480         for( size_t j = 1; j < pd.elementArray[i].node_count(); ++j )
01481             stream << "," << setw( iw ) << pd.elementArray[i].get_vertex_index_array()[j];
01482         stream << endl;
01483     }
01484     stream << endl;
01485 
01486     stream << "Mesh: " << ( pd.myMesh ? "yes" : "no" ) << endl;
01487     stream << "Domain: " << ( pd.myDomain ? "yes" : "no" ) << endl;
01488     //   stream << "mType: " << (pd.mType==PatchData::VERTICES_ON_VERTEX_PATCH?"vert-on-vert":
01489     //                           pd.mType==PatchData::ELEMENTS_ON_VERTEX_PATCH?"elem-on-vert":
01490     //                           pd.mType==PatchData::GLOBAL_PATCH?"global":"unknown") << endl;
01491 
01492     if( pd.haveComputedInfos )
01493     {
01494         stream << "ComputedInfos:" << endl;
01495         if( pd.have_computed_info( PatchData::MIN_UNSIGNED_AREA ) )
01496             stream << "\t MIN_UNSINGED_AREA = " << pd.computedInfos[PatchData::MIN_UNSIGNED_AREA] << endl;
01497         if( pd.have_computed_info( PatchData::MAX_UNSIGNED_AREA ) )
01498             stream << "\t MAX_UNSIGNED_AREA = " << pd.computedInfos[PatchData::MAX_UNSIGNED_AREA] << endl;
01499         if( pd.have_computed_info( PatchData::MIN_EDGE_LENGTH ) )
01500             stream << "\t MIN_EDGE_LENGTH = " << pd.computedInfos[PatchData::MIN_EDGE_LENGTH] << endl;
01501         if( pd.have_computed_info( PatchData::MAX_EDGE_LENGTH ) )
01502             stream << "\t MAX_EDGE_LENGTH = " << pd.computedInfos[PatchData::MAX_EDGE_LENGTH] << endl;
01503         if( pd.have_computed_info( PatchData::MINMAX_SIGNED_DET2D ) )
01504             stream << "\t MINMAX_SIGNED_DET2D = " << pd.computedInfos[PatchData::MINMAX_SIGNED_DET2D] << endl;
01505         if( pd.have_computed_info( PatchData::MINMAX_SIGNED_DET3D ) )
01506             stream << "\t MINMAX_SIGNED_DET3D = " << pd.computedInfos[PatchData::MINMAX_SIGNED_DET3D] << endl;
01507         if( pd.have_computed_info( PatchData::AVERAGE_DET3D ) )
01508             stream << "\t AVERAGE_DET3D = " << pd.computedInfos[PatchData::AVERAGE_DET3D] << endl;
01509     }
01510 
01511     return stream << endl;
01512 }
01513 
01514 void print_patch_data( const PatchData& pd )
01515 {
01516     std::cout << pd << std::endl;
01517 }
01518 
01519 void PatchData::enslave_higher_order_nodes( const size_t* elem_offset_array, unsigned char* vertex_flags,
01520                                             MsqError& ) const
01521 {
01522     for( size_t i = 0; i < elementArray.size(); ++i )
01523     {
01524         size_t start    = elem_offset_array[i];
01525         size_t conn_len = elem_offset_array[i + 1] - start;
01526         for( size_t j = elementArray[i].vertex_count(); j < conn_len; ++j )
01527         {
01528             const size_t vert_idx = elemConnectivityArray[start + j];
01529             assert( vert_idx < vertexHandlesArray.size() );
01530             if( !( vertex_flags[vert_idx] & MsqVertex::MSQ_HARD_FIXED ) )
01531                 vertex_flags[vert_idx] |= MsqVertex::MSQ_DEPENDENT;
01532         }
01533     }
01534 }
01535 
01536 void PatchData::initialize_data( size_t* elem_offset_array, unsigned char* vertex_flags, MsqError& )
01537 {
01538     // Clear out data specific to patch
01539     vertexNormalIndices.clear();
01540     normalData.clear();
01541     // vertexDomainDOF.clear();
01542 
01543     // Clear any vertex->element adjacency data.  It
01544     // is probably invalid, and certainly will be by the time
01545     // this function completes if the mesh contains higher-order
01546     // elements.
01547     vertAdjacencyArray.clear();
01548     vertAdjacencyOffsets.clear();
01549     size_t i, j;
01550     for( i = 0; i < elementArray.size(); ++i )
01551     {
01552         size_t start    = elem_offset_array[i];
01553         size_t conn_len = elem_offset_array[i + 1] - start;
01554         assert( conn_len > 0 );
01555         elementArray[i].set_connectivity( &elemConnectivityArray[start], conn_len );
01556     }
01557 
01558     // Use vertAdjacencyOffsets array as temporary storage.
01559     vertAdjacencyOffsets.resize( vertexHandlesArray.size() + 1 );
01560     size_t* vertex_index_map = arrptr( vertAdjacencyOffsets );
01561 
01562     // Count number of free vertices and initialize vertex_index_map
01563     numFreeVertices = 0;
01564     for( i = 0; i < vertexHandlesArray.size(); ++i )
01565     {
01566         if( !( vertex_flags[i] & MsqVertex::MSQ_FIXED ) ) ++numFreeVertices;
01567         vertex_index_map[i] = i;
01568     }
01569 
01570     // Re-order vertices such that all free vertices are
01571     // first in the list.  Construct map from old to new
01572     // position in list for updating element connectivity.
01573     i = 0;
01574     j = numFreeVertices;
01575     for( ;; ++i, ++j )
01576     {
01577         // find next fixed vertex in the range [i,numFreeVertices)
01578         for( ; i < numFreeVertices && !( vertex_flags[i] & MsqVertex::MSQ_FIXED ); ++i )
01579             ;
01580         // if no more fixed vertices in the free vertex range [0,numFreeVertices)
01581         if( i == numFreeVertices ) break;
01582         // find the next free vertex in the range [j,num_nodes)
01583         for( ; vertex_flags[j] & MsqVertex::MSQ_FIXED; ++j )
01584             ;
01585         // swap fixed (i) and free (j) vertices
01586         vertex_index_map[i] = j;
01587         vertex_index_map[j] = i;
01588         std::swap( vertexHandlesArray[i], vertexHandlesArray[j] );
01589         std::swap( vertex_flags[i], vertex_flags[j] );
01590     }
01591     assert( i == numFreeVertices );
01592     assert( j <= vertexHandlesArray.size() );
01593 
01594     // Update element connectivity data for new vertex indices
01595     for( i = 0; i < elemConnectivityArray.size(); ++i )
01596         elemConnectivityArray[i] = vertex_index_map[elemConnectivityArray[i]];
01597 
01598     // Reorder vertices such that free, slave vertices
01599     // occur after free, non-slave vertices in list.
01600     numSlaveVertices = 0;
01601     for( i = 0; i < vertexHandlesArray.size(); ++i )
01602     {
01603         if( ( vertex_flags[i] & MsqVertex::MSQ_DEPENDENT ) && !( vertex_flags[i] & MsqVertex::MSQ_FIXED ) )
01604             ++numSlaveVertices;
01605     }
01606     numFreeVertices -= numSlaveVertices;
01607 
01608     if( numSlaveVertices )
01609     {
01610         // re-initialize vertex index map
01611         for( i = 0; i < vertexHandlesArray.size(); ++i )
01612             vertex_index_map[i] = i;
01613 
01614         // Re-order free vertices such that all slave vertices are
01615         // last in the list.  Construct map from old to new
01616         // position in list for updating element connectivity.
01617         i = 0;
01618         j = numFreeVertices;
01619         for( ;; ++i, ++j )
01620         {
01621             // find next slave vertex in the range [i,numFreeVertices)
01622             for( ; i < numFreeVertices && !( vertex_flags[i] & MsqVertex::MSQ_DEPENDENT ); ++i )
01623                 ;
01624             // if no more slave vertices in [0,numFreeVertices), then done.
01625             if( i == numFreeVertices ) break;
01626             // find the next free (non-slave) vertex in the range
01627             //   [numFreeVertices,numFreeVertices+numSlaveVertices)
01628             for( ; vertex_flags[j] & MsqVertex::MSQ_DEPENDENT; ++j )
01629                 ;
01630             // swap free (j) and slave (i) vertices
01631             vertex_index_map[i] = j;
01632             vertex_index_map[j] = i;
01633             std::swap( vertexHandlesArray[i], vertexHandlesArray[j] );
01634             std::swap( vertex_flags[i], vertex_flags[j] );
01635         }
01636         assert( i == numFreeVertices );
01637         assert( j <= numFreeVertices + numSlaveVertices );
01638 
01639         // Update element connectivity data for new vertex indices
01640         for( i = 0; i < elemConnectivityArray.size(); ++i )
01641             elemConnectivityArray[i] = vertex_index_map[elemConnectivityArray[i]];
01642     }
01643 
01644     // Clear temporary data
01645     vertAdjacencyOffsets.clear();
01646 
01647     notify_new_patch();
01648 }
01649 
01650 size_t PatchData::num_corners() const
01651 {
01652     size_t result = 0;
01653     for( unsigned i = 0; i < elementArray.size(); ++i )
01654         result += elementArray[i].vertex_count();
01655     return result;
01656 }
01657 
01658 void PatchData::fill( size_t num_vertex, const double* coords, size_t num_elem, EntityTopology type,
01659                       const size_t* connectivity, const bool* fixed, MsqError& err )
01660 {
01661     std::vector< EntityTopology > types( num_elem );
01662     std::fill( types.begin(), types.end(), type );
01663     const EntityTopology* type_ptr = num_elem ? arrptr( types ) : 0;
01664     this->fill( num_vertex, coords, num_elem, type_ptr, connectivity, fixed, err );MSQ_CHKERR( err );
01665 }
01666 
01667 void PatchData::fill( size_t num_vertex, const double* coords, size_t num_elem, const EntityTopology* types,
01668                       const size_t* conn, const bool* fixed, MsqError& err )
01669 {
01670     std::vector< size_t > lengths( num_elem );
01671     std::transform( types, types + num_elem, lengths.begin(), std::ptr_fun( TopologyInfo::corners ) );
01672     const size_t* len_ptr = num_elem ? arrptr( lengths ) : 0;
01673     this->fill( num_vertex, coords, num_elem, types, len_ptr, conn, fixed, err );MSQ_CHKERR( err );
01674 }
01675 
01676 void PatchData::fill( size_t num_vertex, const double* coords, size_t num_elem, const EntityTopology* types,
01677                       const size_t* lengths, const size_t* conn, const bool* fixed, MsqError& err )
01678 {
01679     size_t i;
01680 
01681     // count vertex uses
01682     size_t num_uses = std::accumulate( lengths, lengths + num_elem, 0 );
01683 
01684     // Allocate storage for data
01685     vertexArray.resize( num_vertex );
01686     vertexHandlesArray.resize( num_vertex );
01687     elementArray.resize( num_elem );
01688     elementHandlesArray.resize( num_elem );
01689     elemConnectivityArray.resize( num_uses );
01690     numFreeVertices  = 0;
01691     numSlaveVertices = 0;
01692 
01693     // Must call clear() first so that any stale values get
01694     // zero'd when we call resize.
01695     byteArray.clear();
01696     if( fixed )
01697     {
01698         byteArray.resize( num_vertex, 0 );
01699         for( i = 0; i < num_vertex; ++i )
01700             if( fixed[i] ) byteArray[i] |= ( MsqVertex::MSQ_HARD_FIXED | MsqVertex::MSQ_PATCH_FIXED );
01701     }
01702 
01703     for( i = 0; i < num_elem; ++i )
01704     {
01705         element_by_index( i ).set_element_type( types[i] );
01706         elementHandlesArray[i] = (Mesh::ElementHandle)i;
01707     }
01708     for( i = 0; i < num_vertex; ++i )
01709         vertexHandlesArray[i] = (Mesh::VertexHandle)i;
01710 
01711     memcpy( get_connectivity_array(), conn, num_uses * sizeof( size_t ) );
01712 
01713     std::vector< size_t > offsets( num_elem + 1 );
01714     size_t sum = offsets[0] = 0;
01715     for( i = 1; i <= num_elem; ++i )
01716         offsets[i] = sum += lengths[i - 1];
01717 
01718     const Settings::HigherOrderSlaveMode ho_mode =
01719         mSettings ? mSettings->get_slaved_ho_node_mode() : Settings::SLAVE_ALL;
01720     switch( ho_mode )
01721     {
01722         case Settings::SLAVE_ALL:
01723             byteArray.resize( num_vertex, 0 );
01724             enslave_higher_order_nodes( arrptr( offsets ), arrptr( byteArray ), err );MSQ_ERRRTN( err );
01725             break;
01726         case Settings::SLAVE_NONE:
01727             // Do nothing.  We clear other bits when processing the 'fixed' array above.
01728             break;
01729         default:
01730             MSQ_SETERR( err )
01731             ( "Specified higher-order noded slaving scheme not supported "
01732               "when initializind PatchData using PatchData::fill",
01733               MsqError::NOT_IMPLEMENTED );
01734             return;
01735     }
01736 
01737     this->initialize_data( arrptr( offsets ), arrptr( byteArray ), err );MSQ_ERRRTN( err );
01738 
01739     // initialize_data will re-order vertex handles and
01740     // update element connectivity accordingly.  Use
01741     // the values we stored in vertexHandlesArray to
01742     // figure out the new index of each vertex, and initialize
01743     // the vertex.
01744     for( i = 0; i < num_vertex; ++i )
01745         vertexArray[i] = coords + 3 * (size_t)vertexHandlesArray[i];
01746 
01747     for( i = 0; i < num_vertex; ++i )
01748         vertexArray[i].flags() = byteArray[i];
01749 }
01750 
01751 void PatchData::make_handles_unique( Mesh::EntityHandle* handles, size_t& count, size_t* index_map )
01752 {
01753     if( count < 2 ) { return; }
01754     // save this now, as we'll be changing count later
01755     const size_t* index_end = index_map + count;
01756 
01757     if( index_map )
01758     {
01759         // Copy input handle list into index map as a temporary
01760         // copy of the input list.
01761         assert( sizeof( Mesh::EntityHandle ) == sizeof( size_t ) );
01762         memcpy( index_map, handles, count * sizeof( size_t ) );
01763     }
01764 
01765     // Make handles a unique list
01766     std::sort( handles, handles + count );
01767     Mesh::EntityHandle* end = std::unique( handles, handles + count );
01768     count                   = end - handles;
01769 
01770     if( index_map )
01771     {
01772         // Replace each handle in index_map with the index of
01773         // its position in the handles array
01774         Mesh::EntityHandle* pos;
01775         for( size_t* iter = index_map; iter != index_end; ++iter )
01776         {
01777             pos   = std::lower_bound( handles, handles + count, (Mesh::EntityHandle)*iter );
01778             *iter = pos - handles;
01779         }
01780     }
01781 }
01782 
01783 void PatchData::fill_global_patch( MsqError& err )
01784 {
01785     GlobalPatch gp;
01786     gp.set_mesh( get_mesh() );
01787     PatchIterator iter( &gp );
01788     bool b = iter.get_next_patch( *this, err );MSQ_ERRRTN( err );
01789     if( !b ) MSQ_SETERR( err )( "Empty patch in PatchData::fill_global_patch", MsqError::INVALID_MESH );
01790     assert( b );
01791 }
01792 
01793 void PatchData::set_mesh_entities( std::vector< Mesh::ElementHandle >& elements,
01794                                    std::vector< Mesh::VertexHandle >& free_vertices, MsqError& err )
01795 {
01796     Mesh* current_mesh = get_mesh();
01797     if( !current_mesh )
01798     {
01799         MSQ_SETERR( err )( "No Mesh associated with PatchData.", MsqError::INVALID_STATE );
01800         return;
01801     }
01802 
01803     if( elements.empty() )
01804     {
01805         clear();
01806         return;
01807     }
01808 
01809     elementHandlesArray = elements;
01810     get_mesh()->elements_get_attached_vertices( arrptr( elementHandlesArray ), elementHandlesArray.size(),
01811                                                 vertexHandlesArray, offsetArray, err );MSQ_ERRRTN( err );
01812 
01813     // Construct CSR-rep element connectivity data
01814     size_t num_vert = vertexHandlesArray.size();
01815     elemConnectivityArray.resize( num_vert );
01816     make_handles_unique( arrptr( vertexHandlesArray ), num_vert, arrptr( elemConnectivityArray ) );
01817     vertexHandlesArray.resize( num_vert );
01818 
01819     // Get element topologies
01820     std::vector< EntityTopology > elem_topologies( elementHandlesArray.size() );
01821     get_mesh()->elements_get_topologies( arrptr( elementHandlesArray ), arrptr( elem_topologies ),
01822                                          elementHandlesArray.size(), err );MSQ_ERRRTN( err );
01823 
01824     // get vertex bits from mesh
01825     byteArray.resize( vertexHandlesArray.size() );
01826     get_mesh()->vertices_get_byte( arrptr( vertexHandlesArray ), arrptr( byteArray ), vertexHandlesArray.size(), err );MSQ_ERRRTN( err );
01827 
01828     // if free_vertices is not empty, then we need to mark as
01829     // fixed any vertices *not* in that list.
01830     if( free_vertices.size() == 1 )
01831     {
01832         for( size_t i = 0; i < vertexHandlesArray.size(); ++i )
01833             if( vertexHandlesArray[i] == free_vertices.front() )
01834                 byteArray[i] &= ~MsqVertex::MSQ_PATCH_FIXED;
01835             else
01836                 byteArray[i] |= MsqVertex::MSQ_PATCH_FIXED;
01837     }
01838     else if( !free_vertices.empty() )
01839     {
01840         // sort and remove duplicates from free_vertices list.
01841         std::sort( free_vertices.begin(), free_vertices.end() );
01842         free_vertices.erase( std::unique( free_vertices.begin(), free_vertices.end() ), free_vertices.end() );
01843 
01844         for( size_t i = 0; i < vertexHandlesArray.size(); ++i )
01845         {
01846             if( ( byteArray[i] & MsqVertex::MSQ_DEPENDENT ) ||
01847                 std::binary_search( free_vertices.begin(), free_vertices.end(), vertexHandlesArray[i] ) )
01848                 byteArray[i] &= ~MsqVertex::MSQ_PATCH_FIXED;
01849             else
01850                 byteArray[i] |= MsqVertex::MSQ_PATCH_FIXED;
01851         }
01852     }
01853 
01854     // set element types
01855     elementArray.resize( elementHandlesArray.size() );
01856     for( size_t i = 0; i < elementHandlesArray.size(); ++i )
01857         elementArray[i].set_element_type( elem_topologies[i] );
01858 
01859     // get element connectivity, group vertices by free/slave/fixed state
01860     initialize_data( arrptr( offsetArray ), arrptr( byteArray ), err );MSQ_ERRRTN( err );
01861 
01862     // get vertex coordinates
01863     vertexArray.resize( vertexHandlesArray.size() );
01864     get_mesh()->vertices_get_coordinates( arrptr( vertexHandlesArray ), arrptr( vertexArray ),
01865                                           vertexHandlesArray.size(), err );MSQ_ERRRTN( err );
01866 
01867     // set vertex flags
01868     for( size_t i = 0; i < vertexArray.size(); ++i )
01869         vertexArray[i].flags() = byteArray[i];
01870 }
01871 
01872 void PatchData::get_sample_location( size_t element_index, Sample sample, Vector3D& result, MsqError& err ) const
01873 {
01874     const MsqMeshEntity& elem      = element_by_index( element_index );
01875     const NodeSet ho_bits          = non_slave_node_set( element_index );
01876     const MappingFunction* const f = get_mapping_function( elem.get_element_type() );
01877     if( !f )
01878     {
01879         MSQ_SETERR( err )( "No mapping function", MsqError::UNSUPPORTED_ELEMENT );
01880         return;
01881     }
01882 
01883     double coeff[27];
01884     size_t num_coeff, index[27];
01885     f->coefficients( sample, ho_bits, coeff, index, num_coeff, err );MSQ_ERRRTN( err );
01886     f->convert_connectivity_indices( elem.node_count(), index, num_coeff, err );MSQ_ERRRTN( err );
01887 
01888     const size_t* const conn = elem.get_vertex_index_array();
01889     assert( num_coeff > 0 );
01890     result = coeff[0] * vertex_by_index( conn[index[0]] );
01891     for( unsigned i = 1; i < num_coeff; ++i )
01892         result += coeff[i] * vertex_by_index( conn[index[i]] );
01893 }
01894 
01895 NodeSet PatchData::non_slave_node_set( size_t element_index ) const
01896 {
01897     const MsqMeshEntity& elem = element_by_index( element_index );
01898     EntityTopology type       = elem.get_element_type();
01899     const size_t* conn        = elem.get_vertex_index_array();
01900     const size_t n            = elem.node_count();
01901 
01902     MsqError err;
01903     bool have_midedge, have_midface, have_midelem;
01904     unsigned num_edge = 0, num_face = 0, num_corner = TopologyInfo::corners( type );
01905     TopologyInfo::higher_order( type, n, have_midedge, have_midface, have_midelem, err );
01906     num_edge = TopologyInfo::edges( type );
01907     if( TopologyInfo::dimension( type ) == 2 )
01908         num_face = 1;
01909     else
01910         num_face = TopologyInfo::faces( type );
01911 
01912     NodeSet result;
01913     result.set_all_corner_nodes( type );
01914     if( have_midedge )
01915     {
01916         for( unsigned i = 0; i < num_edge; ++i )
01917         {
01918             if( !( vertex_by_index( conn[num_corner + i] ).get_flags() & MsqVertex::MSQ_DEPENDENT ) )
01919                 result.set_mid_edge_node( i );
01920         }
01921     }
01922     if( have_midface )
01923     {
01924         for( unsigned i = 0; i < num_face; ++i )
01925         {
01926             if( !( vertex_by_index( conn[num_corner + num_edge + i] ).get_flags() & MsqVertex::MSQ_DEPENDENT ) )
01927                 result.set_mid_face_node( i );
01928         }
01929     }
01930     if( have_midelem &&
01931         !( vertex_by_index( conn[num_corner + num_edge + num_face] ).get_flags() & MsqVertex::MSQ_DEPENDENT ) )
01932         result.set_mid_region_node();
01933 
01934     return result;
01935 }
01936 
01937 void PatchData::get_samples( size_t element, std::vector< Sample >& samples, MsqError& ) const
01938 {
01939     NodeSet ns = get_samples( element );
01940     samples.resize( ns.num_nodes() );
01941     std::vector< Sample >::iterator i = samples.begin();
01942 
01943     unsigned j;
01944     EntityTopology type = element_by_index( element ).get_element_type();
01945     for( j = 0; j < TopologyInfo::corners( type ); ++j )
01946         if( ns.corner_node( j ) ) *( i++ ) = Sample( 0, j );
01947     for( j = 0; j < TopologyInfo::edges( type ); ++j )
01948         if( ns.mid_edge_node( j ) ) *( i++ ) = Sample( 1, j );
01949     if( TopologyInfo::dimension( type ) == 3 )
01950     {
01951         for( j = 0; j < TopologyInfo::faces( type ); ++j )
01952             if( ns.mid_face_node( j ) ) *( i++ ) = Sample( 2, j );
01953         if( ns.mid_region_node() ) *( i++ ) = Sample( 3, 0 );
01954     }
01955     else if( ns.mid_face_node( 0 ) )
01956         *( i++ ) = Sample( 2, 0 );
01957 
01958     assert( i == samples.end() );
01959 }
01960 
01961 bool PatchData::attach_extra_data( ExtraData* data )
01962 {
01963     if( data->patchNext ) { return false; }
01964 
01965     if( !data->patchPtr )
01966         data->patchPtr = this;
01967     else if( data->patchPtr != this )
01968         return false;
01969 
01970     data->patchNext = dataList;
01971     dataList        = data;
01972     return true;
01973 }
01974 
01975 bool PatchData::remove_extra_data( ExtraData* data )
01976 {
01977     if( data->patchPtr != this ) return false;
01978 
01979     if( dataList == data )
01980     {
01981         dataList        = data->patchNext;
01982         data->patchNext = 0;
01983         data->patchPtr  = 0;
01984         return true;
01985     }
01986 
01987     for( ExtraData* iter = dataList; iter; iter = iter->patchNext )
01988         if( iter->patchNext == data )
01989         {
01990             iter->patchNext = data->patchNext;
01991             data->patchNext = 0;
01992             data->patchPtr  = 0;
01993             return true;
01994         }
01995 
01996     return false;
01997 }
01998 
01999 void PatchData::notify_new_patch()
02000 {
02001     for( ExtraData* iter = dataList; iter; iter = iter->patchNext )
02002         iter->notify_new_patch();
02003 }
02004 
02005 void PatchData::notify_sub_patch( PatchData& sub_patch, const size_t* vertex_map, const size_t* element_map,
02006                                   MsqError& err )
02007 {
02008     for( ExtraData* iter = dataList; iter; iter = iter->patchNext )
02009     {
02010         iter->notify_sub_patch( sub_patch, vertex_map, element_map, err );MSQ_ERRRTN( err );
02011     }
02012 }
02013 
02014 void PatchData::notify_patch_destroyed()
02015 {
02016     // Remove all ExtraDatas from list and notify them
02017     // that they are being destroyed.
02018     while( dataList )
02019     {
02020         ExtraData* dead_data = dataList;
02021         dataList             = dataList->patchNext;
02022         dead_data->patchNext = 0;
02023         dead_data->patchPtr  = 0;
02024         dead_data->notify_patch_destroyed();
02025     }
02026 }
02027 
02028 }  // namespace MBMesquite
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines