PatchData.cpp

Go to the documentation of this file.

/* *****************************************************************
    MESQUITE -- The Mesh Quality Improvement Toolkit

    Copyright 2004 Sandia Corporation and Argonne National
    Laboratory.  Under the terms of Contract DE-AC04-94AL85000
    with Sandia Corporation, the U.S. Government retains certain
    rights in this software.

    This library is free software; you can redistribute it and/or
    modify it under the terms of the GNU Lesser General Public
    License as published by the Free Software Foundation; either
    version 2.1 of the License, or (at your option) any later version.

    This library is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
    Lesser General Public License for more details.

    You should have received a copy of the GNU Lesser General Public License
    (lgpl.txt) along with this library; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

    [email protected], [email protected], [email protected],
    [email protected], [email protected], [email protected]

  ***************************************************************** */
/*!
  \file   PatchData.cpp

  \author Thomas Leurent
  \author Michael Brewer
  \date   2002-01-17
*/

#include "PatchData.hpp"
#include "MsqMeshEntity.hpp"
#include "MsqFreeVertexIndexIterator.hpp"
#include "MeshInterface.hpp"
#include "MsqTimer.hpp"
#include "MsqDebug.hpp"
#include "GlobalPatch.hpp"
#include "PatchIterator.hpp"
#include "ExtraData.hpp"
#include "Settings.hpp"
#include "MappingFunction.hpp"

#include <list>
#include <vector>
#include <map>
#include <algorithm>
#include <numeric>
#include <functional>
#include <utility>
#include <iostream>
#include <iomanip>
using std::endl;
using std::internal;
using std::left;
using std::list;
using std::map;
using std::ostream;
using std::setfill;
using std::setw;
using std::vector;

namespace MBMesquite
{

const Settings PatchData::defaultSettings;

PatchData::PatchData()
    : myMesh( 0 ), myDomain( 0 ), numFreeVertices( 0 ), numSlaveVertices( 0 ), haveComputedInfos( 0 ), dataList( 0 ),
      mSettings( &defaultSettings )
{
}

// Destructor
PatchData::~PatchData()
{
    notify_patch_destroyed();
}

void PatchData::get_minmax_element_unsigned_area( double& min, double& max, MsqError& err )
{
    if( !have_computed_info( MAX_UNSIGNED_AREA ) || !have_computed_info( MIN_UNSIGNED_AREA ) )
    {
        max          = 0;
        min          = MSQ_DBL_MAX;
        size_t count = num_elements();
        for( size_t i = 0; i < count; ++i )
        {
            double vol;
            assert( i < elementArray.size() );
            vol = elementArray[i].compute_unsigned_area( *this, err );MSQ_ERRRTN( err );
            if( vol > max ) max = vol;
            if( vol < min ) min = vol;
        }
        note_have_info( MAX_UNSIGNED_AREA );
        note_have_info( MIN_UNSIGNED_AREA );
        computedInfos[MAX_UNSIGNED_AREA] = max;
        computedInfos[MIN_UNSIGNED_AREA] = min;
    }
    else
    {
        max = computedInfos[MAX_UNSIGNED_AREA];
        min = computedInfos[MIN_UNSIGNED_AREA];
    }

    if( max <= 0 || min < 0 || min == MSQ_DBL_MAX ) MSQ_SETERR( err )( MsqError::INTERNAL_ERROR );
}

void PatchData::get_minmax_edge_length( double& min, double& max ) const
{
    min = HUGE_VAL;
    max = -HUGE_VAL;

    for( size_t i = 0; i < num_elements(); ++i )
    {
        const MsqMeshEntity& elem = element_by_index( i );
        const size_t* conn        = elem.get_vertex_index_array();
        const unsigned num_edges  = TopologyInfo::edges( elem.get_element_type() );
        for( unsigned e = 0; e < num_edges; ++e )
        {
            const unsigned* edge = TopologyInfo::edge_vertices( elem.get_element_type(), e );
            const double len_sqr =
                ( vertex_by_index( conn[edge[0]] ) - vertex_by_index( conn[edge[1]] ) ).length_squared();
            if( len_sqr < min ) min = len_sqr;
            if( len_sqr > max ) max = len_sqr;
        }
    }
    min = sqrt( min );
    max = sqrt( max );
}

/*
double PatchData::get_average_Lambda_3d( MsqError &err)
{
  double avg;
  if (have_computed_info(AVERAGE_DET3D))
  {
    avg = computedInfos[AVERAGE_DET3D];
  }
  else
  {
    avg =0.;
    int total_num_corners =0;
    Matrix3D A[MSQ_MAX_NUM_VERT_PER_ENT];
    for (size_t i=0; i<elementArray.size(); ++i) {
      int nve = elementArray[i].corner_count();
      elementArray[i].compute_corner_matrices(*this, A, nve, err);
      MSQ_ERRZERO(err);
      total_num_corners += nve;
      for (int c=0; c<nve; ++c) {
        avg += TargetCalculator::compute_Lambda(A[c], err);
        MSQ_ERRZERO(err);
      }
    }

    avg = avg / total_num_corners;
    computedInfos[AVERAGE_DET3D] = avg;
    note_have_info(AVERAGE_DET3D);
  }
  return avg;
}
*/

/*! \fn PatchData::reorder()
   Physically reorder the vertices and elements in the PatchData to improve
   the locality of reference.  This method implements a Reverse Breadth First
   Search order starting with the vertex furthest from the origin.  Other
   orderings can also be implemented.
*/
void PatchData::reorder()
{
    size_t i, j;
    const size_t num_vertex = num_nodes();
    const size_t num_elem   = num_elements();

    // Step 1: Clear any cached data that will be invalidated by this
    vertexNormalIndices.clear();
    normalData.clear();
    // vertexDomainDOF.clear();

    // Step 2: Make sure we have vertex-to-element adjacencies
    if( !vertAdjacencyArray.size() ) generate_vertex_to_element_data();

    // Step 3: Do breadth-first search
    std::vector< bool > visited( num_vertex, false );
    std::vector< size_t > vertex_order( num_vertex );
    std::vector< size_t >::iterator q1_beg, q1_end, q2_end;
    q1_beg = q1_end = q2_end = vertex_order.begin();
    // Outer loop will be done once for each disconnected chunk of mesh.
    while( q1_beg != vertex_order.end() )
    {
        // Find vertex furthest from the origin
        double max     = -1.0;
        size_t vtx_idx = num_vertex;
        for( i = 0; i < num_vertex; ++i )
            if( !visited[i] )
            {
                double dist = vertexArray[i].length_squared();
                if( dist > max )
                {
                    max     = dist;
                    vtx_idx = i;
                }
            }
        assert( vtx_idx < num_vertex );

        *q2_end++ = vtx_idx;
        ;
        visited[vtx_idx] = true;
        do
        {
            q1_end = q2_end;
            for( ; q1_beg != q1_end; ++q1_beg )
            {
                size_t vtx_adj_offset = vertAdjacencyOffsets[*q1_beg];
                size_t vtx_adj_end    = vertAdjacencyOffsets[*q1_beg + 1];
                for( i = vtx_adj_offset; i < vtx_adj_end; ++i )
                {
                    size_t elem = vertAdjacencyArray[i];
                    assert( elem < elementArray.size() );
                    size_t num_elem_verts = elementArray[elem].node_count();
                    size_t* elem_verts    = elementArray[elem].get_vertex_index_array();
                    for( j = 0; j < num_elem_verts; ++j )
                    {
                        size_t elem_vert = elem_verts[j];
                        if( !visited[elem_vert] )
                        {
                            *q2_end++ = elem_vert;
                            ;
                            visited[elem_vert] = true;
                        }
                    }
                }
            }
        } while( q2_end != q1_end );
    }

    // Step 4: vertex_order contains the list of current vertex indices
    //         in the opposite of the order that they will occur in the
    //         reorderd patch.  The following code will construct veretx_map
    //         from vertex_order with the following properties
    //         - vertex_map will be indexed by the current vertex index and
    //           contain the new index of that vertex (inverse of vertex_order)
    //         - the vertices will be grouped by their free/slave/fixed flag.
    std::vector< size_t > vertex_map( num_vertex );
    const size_t fixed_vtx_offset = numFreeVertices + numSlaveVertices;
    size_t free_idx = 0, slave_idx = numFreeVertices, fixed_idx = fixed_vtx_offset;
    for( i = 1; i <= num_vertex; ++i )
    {
        size_t vtx_idx = vertex_order[num_vertex - i];
        if( vtx_idx < numFreeVertices )
            vertex_map[vtx_idx] = free_idx++;
        else if( vtx_idx < fixed_vtx_offset )
            vertex_map[vtx_idx] = slave_idx++;
        else
            vertex_map[vtx_idx] = fixed_idx++;
    }
    // make sure everything adds up
    assert( free_idx == numFreeVertices );
    assert( slave_idx == fixed_vtx_offset );
    assert( fixed_idx == num_vertex );

    // Step 5: compute element permutation
    // initialize all to "num_elem" to indicate unvisited
    std::vector< size_t > element_map( num_elem, num_elem );
    size_t elem_idx = 0;
    for( i = 1; i <= num_vertex; ++i )
    {
        size_t vtx_idx        = vertex_order[num_vertex - i];
        size_t vtx_adj_offset = vertAdjacencyOffsets[vtx_idx];
        size_t vtx_adj_end    = vertAdjacencyOffsets[vtx_idx + 1];
        for( j = vtx_adj_offset; j < vtx_adj_end; ++j )
        {
            size_t elem = vertAdjacencyArray[j];
            if( element_map[elem] == num_elem ) element_map[elem] = elem_idx++;
        }
    }
    // make sure everything adds up
    assert( elem_idx == num_elem );

    // Step 6:  Permute the vertices
    std::vector< MsqVertex > new_vertex_array( num_vertex );
    std::vector< Mesh::VertexHandle > new_vtx_handle_array( num_vertex );
    for( i = 0; i < num_vertex; ++i )
    {
        size_t new_idx                = vertex_map[i];
        new_vertex_array[new_idx]     = vertexArray[i];
        new_vtx_handle_array[new_idx] = vertexHandlesArray[i];
    }
    vertexArray.swap( new_vertex_array );
    vertexHandlesArray.swap( new_vtx_handle_array );

    // Step 7: Permute the elements and vertex indices for the elements
    std::vector< MsqMeshEntity > new_elem_array( num_elem );
    std::vector< Mesh::ElementHandle > new_elem_handle_array( num_elem );
    for( i = 0; i < num_elem; ++i )
    {
        assert( i < elementArray.size() );
        size_t vert_count  = elementArray[i].node_count();
        size_t* conn_array = elementArray[i].get_vertex_index_array();
        for( j = 0; j < vert_count; ++j )
        {
            conn_array[j] = vertex_map[conn_array[j]];
        }

        size_t new_idx = element_map[i];
        assert( new_idx < num_elem );
        new_elem_array[new_idx]        = elementArray[i];
        new_elem_handle_array[new_idx] = elementHandlesArray[i];
    }
    elementArray.swap( new_elem_array );
    elementHandlesArray.swap( new_elem_handle_array );

    // Step 8: Clear no-longer-valid vertex-to-element adjacency info.
    if( vertAdjacencyOffsets.size() )
    {
        vertAdjacencyOffsets.clear();
        vertAdjacencyArray.clear();
        generate_vertex_to_element_data();
    }

    notify_new_patch();
}

/*!
   PatchData::move_free_vertices_constrained() moves the free vertices
   (see MsqVertex::is_free() ) as specified by the search direction (dk)
   and scale factor (step_size). After being moved, the vertices are
   projected onto the appropriate geometry.  Fixed vertices are not moved
   regardless of their corresponding dk direction.
   It is often useful to use the create_coords_momento() function before
   calling this function.
   Compile with -DMSQ_DBG3 to check that fixed vertices
   are not moved with that call.

   \param dk: must be a [nb_vtx] array of Vector3D that contains
   the direction in which to move each vertex. Fixed vertices moving
   direction should be zero, although fixed vertices will not be
   moved regardless of their corresponding dk value.
   \param nb_vtx is the number of vertices to move. must corresponds
   to the number of vertices in the PatchData.
   \param step_size will multiply the moving direction given in dk
   for each vertex.
  */
void PatchData::move_free_vertices_constrained( Vector3D dk[], size_t nb_vtx, double step_size, MsqError& err )
{
    if( nb_vtx != num_free_vertices() )
    {
        MSQ_SETERR( err )
        ( "The directional vector must be of length numVertices.", MsqError::INVALID_ARG );
        return;
    }

    size_t i;
    for( i = 0; i < num_free_vertices(); ++i )
    {
        vertexArray[i] += ( step_size * dk[i] );
        snap_vertex_to_domain( i, err );MSQ_ERRRTN( err );
    }

    if( numSlaveVertices )
    {
        update_slave_node_coordinates( err );MSQ_CHKERR( err );
    }
}

/*! set_free_vertices_constrained is similar to
PatchData::move_free_vertices_constrained() except the original vertex positions
are those stored in the PatchDataVerticesMemento instead of the actual vertex
position stored in the PatchData Vertex array.  The final location is stored
in the PatchData; the PatchDataVerticesMemento is unchanged.

   \param dk: must be a [nb_vtx] array of Vector3D that contains
   the direction in which to move each vertex. Fixed vertices moving
   direction should be zero, although fixed vertices will not be
   moved regardless of their corresponding dk value.
   \param nb_vtx is the number of vertices to move. must corresponds
   to the number of vertices in the PatchData.
   \param step_size will multiply the moving direction given in dk
   for each vertex.
  */
void PatchData::set_free_vertices_constrained( PatchDataVerticesMemento* memento,
                                               Vector3D dk[],
                                               size_t nb_vtx,
                                               double step_size,
                                               MsqError& err )
{
    if( memento->originator != this || nb_vtx != num_free_vertices() )
    {
        MSQ_SETERR( err )( MsqError::INVALID_ARG );
        return;
    }

    size_t i;
    for( i = 0; i < num_free_vertices(); ++i )
    {
        vertexArray[i] = memento->vertices[i] + ( step_size * dk[i] );
        snap_vertex_to_domain( i, err );MSQ_ERRRTN( err );
    }

    if( numSlaveVertices )
    {
        update_slave_node_coordinates( err );MSQ_CHKERR( err );
    }

    // Checks that moving direction is zero for fixed vertices.
    if( MSQ_DBG( 3 ) )
    {
        for( i = 0; i < num_nodes(); ++i )
        {
            if( dk[i].length_squared() != 0.0 )
            {
                MSQ_DBGOUT( 3 ) << "dk[" << i << "]: " << dk[i] << endl;
                MSQ_DBGOUT( 3 ) << "moving a fixed vertex." << endl;
            }
        }
    }
}

static void project_to_plane( Vector3D& vect, const Vector3D& norm )
{
    double len_sqr = norm % norm;
    if( len_sqr > DBL_EPSILON ) vect -= norm * ( ( norm % vect ) / len_sqr );
}

void PatchData::project_gradient( std::vector< Vector3D >& gradient, MsqError& err )
{
    if( !domain_set() ) return;

    if( gradient.size() != num_free_vertices() )
    {
        MSQ_SETERR( err )( MsqError::INVALID_ARG );
        return;
    }

    if( normalData.empty() )
    {
        update_cached_normals( err );MSQ_ERRRTN( err );
    }

    Vector3D norm;
    for( size_t i = 0; i < num_free_vertices(); ++i )
    {
        if( vertexNormalIndices.empty() )
        {  // whole mesh on single 2D domain
            norm = normalData[i];
            project_to_plane( gradient[i], norm );
        }
        else if( vertexDomainDOF[i] == 2 )
        {  // vertex on surface
            norm = normalData[vertexNormalIndices[i]];
            project_to_plane( gradient[i], norm );
        }
        else if( vertexDomainDOF[i] == 1 )
        {
            size_t j, num_elem;
            const size_t* elems = get_vertex_element_adjacencies( i, num_elem, err );MSQ_ERRRTN( err );
            for( j = 0; j < num_elem; ++j )
            {
                if( 2 == TopologyInfo::dimension( element_by_index( elems[j] ).get_element_type() ) )
                {
                    norm = vertexArray[i];
                    get_domain()->element_normal_at( elementHandlesArray[elems[j]], norm );
                    project_to_plane( gradient[i], norm );
                }
            }
        }
    }
}

/*! Finds the maximum movement (in the distance norm) of the vertices in a
  patch.  The previous vertex positions are givena as a
  PatchDataVerticesMemento (memento).  The distance squared which each
  vertex has moved is then calculated, and the largest of those distances
  is returned.  This function only considers the movement of vertices
  that are currently 'free'.
  \param memento  a memento of this patch's vertex position at some
  (prior) time in the optimization.
      */
double PatchData::get_max_vertex_movement_squared( PatchDataVerticesMemento* memento, MsqError& )
{
    double max_dist = 0.0;
    for( size_t i = 0; i < memento->vertices.size(); ++i )
    {
        double temp_dist = ( vertexArray[i] - memento->vertices[i] ).length_squared();
        if( temp_dist > max_dist ) max_dist = temp_dist;
    }
    return max_dist;
}

/*!
 */
void PatchData::set_all_vertices_soft_fixed( MsqError& /*err*/ )
{
    for( size_t i = 0; i < num_free_vertices(); ++i )
        vertexArray[i].set_soft_fixed_flag();
}

/*!
 */
void PatchData::set_free_vertices_soft_fixed( MsqError& /*err*/ )
{
    for( size_t i = 0; i < num_free_vertices(); ++i )
    {
        if( vertexArray[i].is_free_vertex() ) vertexArray[i].set_soft_fixed_flag();
    }
}

/*!
 */
void PatchData::set_all_vertices_soft_free( MsqError& /*err*/ )
{
    for( size_t i = 0; i < num_nodes(); ++i )
        vertexArray[i].remove_soft_fixed_flag();
}

/*! Get coordinates of element vertices, in canonical order.

    \param elem_index The element index in the Patch
    \param coords This vector will have the coordinates appended to it.
    If necessary, make sure to clear the vector before calling the function.
  */
void PatchData::get_element_vertex_coordinates( size_t elem_index, std::vector< Vector3D >& coords, MsqError& /*err*/ )
{
    // Check index
    if( elem_index >= num_elements() ) return;

    // Ask the element for its vertex indices
    assert( elem_index < elementArray.size() );
    const size_t* vertex_indices = elementArray[elem_index].get_vertex_index_array();

    // Get the coords for each indicated vertex
    size_t num_verts = elementArray[elem_index].vertex_count();
    coords.reserve( coords.size() + num_verts );
    for( size_t i = 0; i < num_verts; i++ )
        coords.push_back( Vector3D( vertexArray[vertex_indices[i]] ) );
}

/*! This is an inefficient way of retrieving vertex_indices.
    Use PatchData::get_element_array followed by
    MsqMeshEntity::get_vertex_index_array() if you don't need
    to fill an STL vector.
*/
void PatchData::get_element_vertex_indices( size_t elem_index,
                                            std::vector< size_t >& vertex_indices,
                                            MsqError& /*err*/ )
{
    // Ask the element for its vertex indices
    assert( elem_index < elementArray.size() );
    elementArray[elem_index].get_vertex_indices( vertex_indices );
}

void PatchData::get_vertex_element_indices( size_t vertex_index, std::vector< size_t >& elem_indices, MsqError& err )
{
    size_t count;
    const size_t* ptr;
    ptr = get_vertex_element_adjacencies( vertex_index, count, err );
    elem_indices.resize( count );
    std::copy( ptr, ptr + count, elem_indices.begin() );
}

void PatchData::get_vertex_element_indices( size_t vertex_index,
                                            unsigned element_dimension,
                                            std::vector< size_t >& elem_indices,
                                            MsqError& err )
{
    elem_indices.clear();
    size_t count;
    const size_t* ptr;
    ptr = get_vertex_element_adjacencies( vertex_index, count, err );
    for( const size_t* const end = ptr + count; ptr != end; ++ptr )
    {
        assert( *ptr < elementArray.size() );
        const EntityTopology type = elementArray[*ptr].get_element_type();
        const unsigned dim        = TopologyInfo::dimension( type );
        if( dim == element_dimension ) elem_indices.push_back( *ptr );
    }
}

const size_t* PatchData::get_vertex_element_adjacencies( size_t vertex_index, size_t& array_len_out, MsqError& )
{
    // Make sure we've got the data
    if( vertAdjacencyArray.empty() )
    {
        generate_vertex_to_element_data();
    }

    const size_t begin = vertAdjacencyOffsets[vertex_index];
    const size_t end   = vertAdjacencyOffsets[vertex_index + 1];
    array_len_out      = end - begin;
    return &vertAdjacencyArray[begin];
}

/*!
    \brief This function fills a vector<size_t> with the indices
    to vertices connected to the given vertex by an edge.  If vert_indices
    is not initially empty, the function will not delete the current
    contents.  Instead, it will append the new indices at the end of
    the vector.

*/
void PatchData::get_adjacent_vertex_indices( size_t vertex_index,
                                             std::vector< size_t >& vert_indices,
                                             MsqError& err ) const
{
    bitMap.clear();
    bitMap.resize( num_nodes(), false );

    const size_t* conn;
    size_t conn_idx, curr_vtx_idx;
    const unsigned* adj;
    unsigned num_adj, i;
    std::vector< MsqMeshEntity >::const_iterator e;
    for( e = elementArray.begin(); e != elementArray.end(); ++e )
    {
        conn     = e->get_vertex_index_array();
        conn_idx = std::find( conn, conn + e->node_count(), vertex_index ) - conn;
        if( conn_idx == e->node_count() ) continue;

        // If a higher-order node, return corners of side/face
        // that node is in the center of.
        EntityTopology type = e->get_element_type();
        if( conn_idx >= TopologyInfo::corners( type ) && type != POLYGON )
        {
            unsigned dim, id;
            TopologyInfo::side_from_higher_order( type, e->node_count(), conn_idx, dim, id, err );MSQ_ERRRTN( err );
            adj = TopologyInfo::side_vertices( type, dim, id, num_adj );
        }
        else
        {
            EntityTopology topo = e->get_element_type();
            if( topo == POLYGON )
            {
                unsigned number_of_nodes = e->node_count();
                num_adj                  = 2;  // always 2 for a polygon
                unsigned vert_adj[2];
                vert_adj[0] = ( conn_idx + 1 ) % number_of_nodes;
                vert_adj[1] = ( conn_idx + number_of_nodes - 1 ) % number_of_nodes;
                for( i = 0; i < num_adj; ++i )
                {
                    curr_vtx_idx = conn[vert_adj[i]];  // get index into patch vertex list
                    if( !bitMap[curr_vtx_idx] )
                    {
                        vert_indices.push_back( curr_vtx_idx );
                        bitMap[curr_vtx_idx] = true;
                    }
                }
            }
            else
            {
                adj = TopologyInfo::adjacent_vertices( topo, conn_idx, num_adj );
                for( i = 0; i < num_adj; ++i )
                {
                    curr_vtx_idx = conn[adj[i]];  // get index into patch vertex list
                    if( !bitMap[curr_vtx_idx] )
                    {
                        vert_indices.push_back( curr_vtx_idx );
                        bitMap[curr_vtx_idx] = true;
                    }
                }
            }
        }
    }
}

/*! Fills a vector of indices into the entities array. The entities
    in the vector are connected the given entity (ent_ind) via an
    n-diminsional entity (where 'n' is a given integer).
    Thus, if n = 0, the entities must be connected via a vertex.
    If n = 1, the entities must be connected via an edge.
    If n = 2, the entities must be connected via a two-dimensional element.
    NOTE:  if n is 2 and the elements in the entity array are
    two-dimensional, no entities should meet this criterion.
    The adj_ents vector is cleared at the beginning of the call.

*/
void PatchData::get_adjacent_entities_via_n_dim( int n, size_t ent_ind, std::vector< size_t >& adj_ents, MsqError& err )
{
    // reset the vector
    adj_ents.clear();
    // vertices of this entity (given by ent_ind)
    vector< size_t > verts;
    // vector to store elements attached to the vertices in verts
    vector< size_t > elem_on_vert[MSQ_MAX_NUM_VERT_PER_ENT];
    // length of above vectos
    int length_elem_on_vert[MSQ_MAX_NUM_VERT_PER_ENT];
    // get verts on this element
    get_element_vertex_indices( ent_ind, verts, err );
    int num_vert = verts.size();
    int i        = 0;
    int j        = 0;
    for( i = 0; i < num_vert; ++i )
    {
        // get elements on the vertices in verts and the number of vertices
        get_vertex_element_indices( verts[i], elem_on_vert[i], err );MSQ_ERRRTN( err );
        length_elem_on_vert[i] = elem_on_vert[i].size();
    }
    // this_ent is the index for an entity which is a candidate to be placed
    // into adj_ents
    size_t this_ent;
    // num of times this_ent has been found in the vectors of entity indices
    int counter = 0;
    i           = 0;
    // loop of each vert on ent_ind
    while( i < num_vert )
    {
        // loop over each ent connected to vert
        j = 0;
        while( j < length_elem_on_vert[i] )
        {
            // get candidate element
            this_ent = elem_on_vert[i][j];
            // if we haven't already consider this ent
            if( this_ent != ent_ind )
            {
                // if this_ent occurred earlier we would have already considered it
                // so start at i and j+1
                int k1 = i;
                int k2 = j + 1;
                // this_ent has occured once so far
                counter = 1;
                // while k1 < num_vert
                while( k1 < num_vert )
                {
                    // loop over entries in the elem on vert vector
                    while( k2 < length_elem_on_vert[k1] )
                    {
                        // if it matches this_ent
                        if( elem_on_vert[k1][k2] == this_ent )
                        {
                            // mark it as 'seen', by making it the same as ent_ind
                            // i.e., the entity  passed to us.
                            elem_on_vert[k1][k2] = ent_ind;
                            ++counter;
                            // do not look at remaining elems in this vector
                            k2 += length_elem_on_vert[k1];
                        }
                        else
                            ++k2;
                    }
                    ++k1;
                    k2 = 0;
                }
                // if this_ent occured enough times and isn't ent_ind
                if( counter > n && this_ent != ent_ind )
                {
                    adj_ents.push_back( this_ent );
                }
            }
            ++j;
        }
        ++i;
    }
}

/*! \fn PatchData::update_mesh(MsqError &err)

    \brief This function copies to the TSTT mesh  the changes made to the
    free vertices / elements of the PatchData object.

*/
void PatchData::update_mesh( MsqError& err, const TagHandle* tag )
{
    if( !myMesh )
    {
        MSQ_SETERR( err )( "Cannot update mesh on temporary patch.\n", MsqError::INTERNAL_ERROR );
        return;
    }

    const size_t not_fixed = numFreeVertices + numSlaveVertices;
    if( tag )
    {
        for( size_t i = 0; i < not_fixed; ++i )
        {
            myMesh->tag_set_vertex_data( *tag, 1, &vertexHandlesArray[i], vertexArray[i].to_array(), err );MSQ_ERRRTN( err );
        }
    }
    else
    {
        for( size_t i = 0; i < not_fixed; ++i )
        {
            myMesh->vertex_set_coordinates( vertexHandlesArray[i], vertexArray[i], err );MSQ_ERRRTN( err );
        }
    }

    for( size_t i = 0; i < vertexArray.size(); ++i )
    {
        myMesh->vertex_set_byte( vertexHandlesArray[i], vertexArray[i].get_flags(), err );MSQ_ERRRTN( err );
    }
}

void PatchData::update_slave_node_coordinates( MsqError& err )
{
    // update slave vertices
    if( 0 == num_slave_vertices() ) return;

    // Set a mark on every slave vertex.  We'll clear the marks as we
    // set the vertex coordinates.  This way we can check that all
    // vertices got updated.
    const size_t vert_end = num_free_vertices() + num_slave_vertices();
    for( size_t i = num_free_vertices(); i < vert_end; ++i )
        vertexArray[i].flags() |= MsqVertex::MSQ_MARK;

    // For each element, calculate slave vertex coordinates from
    // mapping function.
    const int ARR_SIZE = 27;
    double coeff[ARR_SIZE];
    size_t index[ARR_SIZE];
    for( size_t i = 0; i < num_elements(); ++i )
    {
        MsqMeshEntity& elem  = element_by_index( i );
        const int num_corner = elem.vertex_count();
        const int num_node   = elem.node_count();
        assert( num_node < ARR_SIZE );

        const EntityTopology type       = elem.get_element_type();
        const MappingFunction* const mf = get_mapping_function( type );
        if( 0 == mf || num_node == num_corner ) continue;

        const size_t* conn = elem.get_vertex_index_array();

        // for each higher-order non-slave node, set bit indicating
        // that mapping function is a function of the non-slave node
        // coordinates
        NodeSet ho_bits = non_slave_node_set( i );

        // for each higher-order slave node
        for( int k = num_corner; k < num_node; ++k )
        {
            if( !is_vertex_slave( conn[k] ) ) continue;

            // check if we already did this one for an adjacent element
            MsqVertex& vert = vertexArray[conn[k]];
            if( !vert.is_flag_set( MsqVertex::MSQ_MARK ) ) continue;

            // what is this node a mid node of (i.e. face 1, edge 2, etc.)
            Sample loc = TopologyInfo::sample_from_node( type, elem.node_count(), k, err );MSQ_ERRRTN( err );

            // evaluate mapping function at logical loction of HO node.
            size_t num_coeff;
            mf->coefficients( loc, ho_bits, coeff, index, num_coeff, err );MSQ_ERRRTN( err );
            mf->convert_connectivity_indices( num_node, index, num_coeff, err );MSQ_ERRRTN( err );

            // calulate new coordinates for slave node
            assert( num_coeff > 0 );
            vert = coeff[0] * vertex_by_index( conn[index[0]] );
            for( size_t j = 1; j < num_coeff; ++j )
                vert += coeff[j] * vertex_by_index( conn[index[j]] );

            // clear mark
            vert.flags() &= ~MsqVertex::MSQ_MARK;
        }
    }

    // make sure we set the coordinates for every slave node
    for( size_t i = num_free_vertices(); i < vert_end; ++i )
    {
        if( vertex_by_index( i ).is_flag_set( MsqVertex::MSQ_MARK ) )
        {
            MSQ_SETERR( err )
            ( MsqError::INVALID_MESH, "No element with mapping function adjacent to slave vertex %lu (%lu)\n",
              (unsigned long)i, (unsigned long)get_vertex_handles_array()[i] );
            // make sure we finish with all marks cleared
            vertexArray[i].flags() &= ~MsqVertex::MSQ_MARK;
        }
    }

    // snap vertices to domain
    if( domain_set() )
    {
        for( size_t i = num_free_vertices(); i < vert_end; ++i )
        {
            snap_vertex_to_domain( i, err );MSQ_ERRRTN( err );
        }
    }
}

void PatchData::update_slave_node_coordinates( const size_t* elements, size_t num_elems, MsqError& err )
{
    // update slave vertices
    if( 0 == num_slave_vertices() ) return;

    // set a mark on each vertex so we don't update shared
    // vertices more than once.
    for( size_t i = 0; i < num_elems; ++i )
    {
        MsqMeshEntity& elem  = element_by_index( elements[i] );
        const int num_corner = elem.vertex_count();
        const int num_node   = elem.node_count();
        const size_t* conn   = elem.get_vertex_index_array();
        for( int j = num_corner; j < num_node; ++j )
            vertexArray[conn[j]].flags() |= MsqVertex::MSQ_MARK;
    }

    // For each element, calculate slave vertex coordinates from
    // mapping function.
    const int ARR_SIZE = 27;
    double coeff[ARR_SIZE];
    size_t index[ARR_SIZE];
    for( size_t i = 0; i < num_elems; ++i )
    {
        MsqMeshEntity& elem  = element_by_index( elements[i] );
        const int num_corner = elem.vertex_count();
        const int num_node   = elem.node_count();
        assert( num_node < ARR_SIZE );

        const EntityTopology type       = elem.get_element_type();
        const MappingFunction* const mf = get_mapping_function( type );
        if( 0 == mf || num_node == num_corner ) continue;

        const size_t* conn = elem.get_vertex_index_array();

        // for each higher-order non-slave node, set bit indicating
        // that mapping function is a function of the non-slave node
        // coordinates
        NodeSet ho_bits = non_slave_node_set( i );

        // for each higher-order slave node
        for( int k = num_corner; k < num_node; ++k )
        {
            if( !is_vertex_slave( conn[k] ) ) continue;

            // check if we already did this one for an adjacent element
            MsqVertex& vert = vertexArray[conn[k]];
            if( !vert.is_flag_set( MsqVertex::MSQ_MARK ) ) continue;

            // what is this node a mid node of (i.e. face 1, edge 2, etc.)
            Sample loc = TopologyInfo::sample_from_node( type, elem.node_count(), k, err );MSQ_ERRRTN( err );

            // evaluate mapping function at logical loction of HO node.
            size_t num_coeff;
            mf->coefficients( loc, ho_bits, coeff, index, num_coeff, err );MSQ_ERRRTN( err );
            mf->convert_connectivity_indices( num_node, index, num_coeff, err );MSQ_ERRRTN( err );

            // calulate new coordinates for slave node
            assert( num_coeff > 0 );
            vert = coeff[0] * vertex_by_index( conn[index[0]] );
            for( size_t j = 1; j < num_coeff; ++j )
                vert += coeff[j] * vertex_by_index( conn[index[j]] );

            // snap vertices to domain
            if( domain_set() )
            {
                snap_vertex_to_domain( conn[k], err );MSQ_ERRRTN( err );
            }

            // clear mark
            vert.flags() &= ~MsqVertex::MSQ_MARK;
        }
    }
}

void PatchData::generate_vertex_to_element_data()
{
    MSQ_FUNCTION_TIMER( "PatchData::generate_vertex_to_element_data" );

    // Skip if data already exists
    if( !vertAdjacencyArray.empty() ) return;

    // Skip if patch is empty
    if( 0 == num_nodes() ) return;

    // Allocate offset array
    vertAdjacencyOffsets.clear();
    vertAdjacencyOffsets.resize( num_nodes() + 1, 0 );

    // Temporarily use offsets array to hold per-vertex element count
    std::vector< MsqMeshEntity >::iterator elem_iter;
    const std::vector< MsqMeshEntity >::iterator elem_end = elementArray.end();
    for( elem_iter = elementArray.begin(); elem_iter != elem_end; ++elem_iter )
    {
        size_t* conn_iter      = elem_iter->get_vertex_index_array();
        const size_t* conn_end = conn_iter + elem_iter->node_count();
        for( ; conn_iter != conn_end; ++conn_iter )
            ++vertAdjacencyOffsets[*conn_iter];
    }

    // Convert counts to end indices.
    // When done, vertAdjacencyOffsets will contain, for each vertex,
    // one more than the *last* index for that vertex's data in the
    // adjacency array.  This is *not* the final state for this data.
    // See comments for next loop.
    std::vector< size_t >::iterator off_iter      = vertAdjacencyOffsets.begin();
    const std::vector< size_t >::iterator off_end = vertAdjacencyOffsets.end();
    size_t prev                                   = *off_iter;
    ++off_iter;
    for( ; off_iter != off_end; ++off_iter )
    {
        prev += *off_iter;
        *off_iter = prev;
    }

    // Allocate space for element numbers
    const size_t num_vert_uses = vertAdjacencyOffsets[num_nodes() - 1];
    assert( num_vert_uses == elemConnectivityArray.size() );
    vertAdjacencyArray.resize( num_vert_uses );

    // Fill vertAdjacencyArray, using the indices in vertAdjacencyOffsets
    // as the location to insert the next element number in
    // vertAdjacencyArray.  When done, vertAdjacenyOffsets will contain
    // the start index for each vertex, rather than one past the last
    // index.
    for( size_t i = 0; i < elementArray.size(); ++i )
    {
        size_t* conn_iter      = elementArray[i].get_vertex_index_array();
        const size_t* conn_end = conn_iter + elementArray[i].node_count();
        for( ; conn_iter != conn_end; ++conn_iter )
        {
            const size_t array_index        = --vertAdjacencyOffsets[*conn_iter];
            vertAdjacencyArray[array_index] = i;
        }
    }

    // Last entry should be number of vertex uses (one past the
    // last index of the last vertex.)
    vertAdjacencyOffsets[num_nodes()] = num_vert_uses;
}

void PatchData::get_subpatch( size_t center_vertex_index,
                              unsigned num_adj_elem_layers,
                              PatchData& subpatch,
                              MsqError& err )
{
    // Make sure we're in range
    if( center_vertex_index >= num_free_vertices() )
    {
        MSQ_SETERR( err )( "Invalid index for center vertex", MsqError::INVALID_ARG );
        return;
    }

    // Notify any observers of the existing subpatch that the mesh
    // in the patch is to be changed.
    subpatch.notify_new_patch();

    // Get list of vertices and elements in subpatch.
    // Ultimately, end up with arrays of unique, sorted indices.
    // It is important that the vertex indices be sorted so later
    // a reverse lookup can be done using a binary search (std::lower_bound).
    std::vector< size_t > elements, vertices, offsets;
    vertices.push_back( center_vertex_index );
    for( unsigned i = 0; i < num_adj_elem_layers; ++i )
    {
        elements.clear();
        for( unsigned v = 0; v < vertices.size(); ++v )
        {
            size_t num_elem;
            const size_t* vert_elems = get_vertex_element_adjacencies( vertices[v], num_elem, err );MSQ_ERRRTN( err );
            elements.insert( elements.end(), vert_elems, vert_elems + num_elem );
        }
        std::sort( elements.begin(), elements.end() );
        elements.erase( std::unique( elements.begin(), elements.end() ), elements.end() );

        vertices.clear();
        for( unsigned e = 0; e < elements.size(); ++e )
        {
            MsqMeshEntity& elem      = element_by_index( elements[e] );
            size_t num_vert          = elem.node_count();
            const size_t* elem_verts = elem.get_vertex_index_array();
            vertices.insert( vertices.end(), elem_verts, elem_verts + num_vert );
        }
        std::sort( vertices.begin(), vertices.end() );
        vertices.erase( std::unique( vertices.begin(), vertices.end() ), vertices.end() );
    }

    // Allocate space for element connectivity info.
    size_t num_vert_uses = 0;
    for( unsigned i = 0; i < elements.size(); ++i )
        num_vert_uses += element_by_index( elements[i] ).node_count();
    subpatch.elementArray.resize( elements.size() );
    subpatch.elementHandlesArray.resize( elements.size() );
    subpatch.elemConnectivityArray.resize( num_vert_uses );
    offsets.resize( elements.size() + 1 );

    // Construct element connectivity data in new patch,
    // and copy element type into new patch
    size_t curr_offset = 0;
    for( unsigned i = 0; i < elements.size(); ++i )
    {
        MsqMeshEntity& elem = element_by_index( elements[i] );
        assert( i < elementArray.size() );
        subpatch.elementArray[i].set_element_type( elem.get_element_type() );
        subpatch.elementHandlesArray[i] = elementHandlesArray[elements[i]];
        const size_t* verts             = elem.get_vertex_index_array();
        offsets[i]                      = curr_offset;
        for( unsigned j = 0; j < elem.node_count(); ++j )
        {
            subpatch.elemConnectivityArray[curr_offset++] =
                std::lower_bound( vertices.begin(), vertices.end(), verts[j] ) - vertices.begin();
        }
    }
    offsets[elements.size()] = curr_offset;

    // Store index in this patch in vertex handle array of subpatch
    // so we can determine how vertices were reordered when setting
    // vertex coordinates.
    assert( sizeof( size_t ) == sizeof( void* ) );
    subpatch.vertexHandlesArray.resize( vertices.size() );
    size_t* vert_handles = reinterpret_cast< size_t* >( &subpatch.vertexHandlesArray[0] );
    std::copy( vertices.begin(), vertices.end(), vert_handles );

    // All vertices except vertex at center_vertex_index are fixed.
    subpatch.byteArray.resize( vertices.size() );
    for( size_t pi = 0; pi < vertices.size(); ++pi )
    {
        if( vertices[pi] == center_vertex_index )
            subpatch.byteArray[pi] = vertexArray[vertices[pi]].get_flags() & ~MsqVertex::MSQ_PATCH_FIXED;
        else
            subpatch.byteArray[pi] = vertexArray[vertices[pi]].get_flags() | MsqVertex::MSQ_PATCH_FIXED;
    }

    // Re-order vertices and initialize other data in subpatch
    subpatch.initialize_data( arrptr( offsets ), &subpatch.byteArray[0], err );MSQ_ERRRTN( err );

    // Copy vertex data into subpatch.  subpatch.vertexHandlesArray contains
    // the indices into this PatchData for each vertex, as reordered by the
    // call to initialize_data.
    subpatch.vertexArray.resize( vertices.size() );
    for( unsigned i = 0; i < vertices.size(); ++i )
    {
        size_t vert_index               = (size_t)( subpatch.vertexHandlesArray[i] );
        vertices[i]                     = vert_index;
        subpatch.vertexHandlesArray[i]  = vertexHandlesArray[vert_index];
        subpatch.vertexArray[i]         = vertexArray[vert_index];
        subpatch.vertexArray[i].flags() = subpatch.byteArray[i];
    }

    subpatch.myMesh    = myMesh;
    subpatch.myDomain  = myDomain;
    subpatch.mSettings = mSettings;

    notify_sub_patch( subpatch, arrptr( vertices ), elements.empty() ? 0 : arrptr( elements ), err );MSQ_CHKERR( err );
}

//! Adjust the position of the specified vertex so that it
//! lies on its constraining domain.  The actual domain constraint
//! is managed by the TSTT mesh implementation
void PatchData::snap_vertex_to_domain( size_t vertex_index, MsqError& err )
{
    if( domain_set() )
    {
        // if not doing normal caching or vertex is not on a single surface
        if( normalData.empty() )
        {
            get_domain()->snap_to( vertexHandlesArray[vertex_index], vertexArray[vertex_index] );
        }
        // entire domain is 2-D (all vertices have a single normal)
        else if( vertexNormalIndices.empty() )
        {
            get_domain()->closest_point( vertexHandlesArray[vertex_index], Vector3D( vertexArray[vertex_index] ),
                                         vertexArray[vertex_index], normalData[vertex_index], err );MSQ_ERRRTN( err );
        }
        else if( vertexNormalIndices[vertex_index] < normalData.size() )
        {  // vertex has a unique normal
            get_domain()->closest_point( vertexHandlesArray[vertex_index], Vector3D( vertexArray[vertex_index] ),
                                         vertexArray[vertex_index], normalData[vertexNormalIndices[vertex_index]],
                                         err );MSQ_ERRRTN( err );
        }
        else
        {  // vertex has multiple normals
            get_domain()->snap_to( vertexHandlesArray[vertex_index], vertexArray[vertex_index] );
        }
    }
}

void PatchData::update_cached_normals( MsqError& err )
{
    size_t i;

    MeshDomain* domain = get_domain();
    if( !domain )
    {
        MSQ_SETERR( err )( "No domain constraint set.", MsqError::INVALID_STATE );
        return;
    }

    // Determine which vertices lie on surfaces
    vertexDomainDOF.resize( num_nodes() );
    domain->domain_DoF( arrptr( vertexHandlesArray ), arrptr( vertexDomainDOF ), num_nodes(), err );MSQ_ERRRTN( err );

    // Count how many vertices have a single normal
    // Sun doesn't support partial template specialization, so can't use std::count
    // size_t n = std::count( vertexDomainDOF.begin(), vertexDomainDOF.end(), 2 );
    size_t n = 0;
    std::vector< unsigned short >::iterator k;
    for( k = vertexDomainDOF.begin(); k != vertexDomainDOF.end(); ++k )
        if( *k == 2 ) ++n;

    normalData.resize( n );

    // If all vertices are on a surface, pass in the existing handles array
    // and store a single normal per vertex.
    if( n == num_nodes() )
    {
        std::copy( vertexArray.begin(), vertexArray.end(), normalData.begin() );
        domain->vertex_normal_at( arrptr( vertexHandlesArray ), arrptr( normalData ), num_nodes(), err );
        vertexNormalIndices.clear();
        vertexDomainDOF.clear();MSQ_ERRRTN( err );
    }
    else
    {
        vertexNormalIndices.resize( num_nodes() );
        size_t nn = 0;
        for( i = 0; i < num_nodes(); ++i )
        {
            if( vertexDomainDOF[i] == 2 )
            {
                normalData[nn] = vertexArray[i];
                domain->vertex_normal_at( vertexHandlesArray[i], normalData[nn] );
                vertexNormalIndices[i] = nn;
                ++nn;
            }
            else
            {
                vertexNormalIndices[i] = n + 1;
            }
        }
        assert( nn == n );
    }
}

void PatchData::get_domain_normal_at_element( size_t elem_index, Vector3D& surf_norm, MsqError& err )
{
    // check if element as a mid-face node
    const MsqMeshEntity& elem = element_by_index( elem_index );
    const int mid_node = TopologyInfo::higher_order_from_side( elem.get_element_type(), elem.node_count(), 2, 0, err );MSQ_ERRRTN( err );
    // if we have the mid-element vertex, get cached normal for it
    if( mid_node > 0 )
    {
        get_domain_normal_at_vertex( elem.get_vertex_index_array()[mid_node], elementHandlesArray[elem_index],
                                     surf_norm, err );MSQ_ERRRTN( err );
    }
    // otherwise query domain for normal at element centroid
    else if( domain_set() )
    {
        assert( elem_index < elementArray.size() );
        elementArray[elem_index].get_centroid( surf_norm, *this, err );MSQ_ERRRTN( err );
        get_domain()->element_normal_at( elementHandlesArray[elem_index], surf_norm );
    }
    else
        MSQ_SETERR( err )( "No domain constraint set.", MsqError::INVALID_STATE );
}

void PatchData::get_domain_normal_at_mid_edge( size_t elem_index, unsigned edge_num, Vector3D& normal, MsqError& err )
{
    // check if element as a mid-edge node
    const MsqMeshEntity& elem = element_by_index( elem_index );
    const int mid_node =
        TopologyInfo::higher_order_from_side( elem.get_element_type(), elem.node_count(), 1, edge_num, err );MSQ_ERRRTN( err );
    // if we have the mid-edge vertex, get cached normal for it
    if( mid_node > 0 )
    {
        get_domain_normal_at_vertex( elem.get_vertex_index_array()[mid_node], elementHandlesArray[elem_index], normal,
                                     err );MSQ_ERRRTN( err );
    }
    // otherwise query domain for normal at center of edge
    else if( domain_set() )
    {
        const unsigned* edge = TopologyInfo::edge_vertices( elem.get_element_type(), edge_num, err );MSQ_ERRRTN( err );
        const MsqVertex& v1 = vertex_by_index( elem.get_vertex_index_array()[edge[0]] );
        const MsqVertex& v2 = vertex_by_index( elem.get_vertex_index_array()[edge[1]] );
        normal              = 0.5 * ( v1 + v2 );
        get_domain()->element_normal_at( elementHandlesArray[elem_index], normal );
    }
    else
    {
        MSQ_SETERR( err )( "No domain constraint set.", MsqError::INVALID_STATE );
        return;
    }
}

void PatchData::get_domain_normals_at_corners( size_t elem_index, Vector3D normals_out[], MsqError& err )
{
    if( !domain_set() )
    {
        MSQ_SETERR( err )( "No domain constraint set.", MsqError::INVALID_STATE );
        return;
    }

    assert( elem_index < elementArray.size() );
    if( 2 != TopologyInfo::dimension( elementArray[elem_index].get_element_type() ) )
    {
        MSQ_SETERR( err )( "Attempt to get corners of non-surface element", MsqError::INVALID_ARG );
        return;
    }

    if( normalData.empty() )
    {
        update_cached_normals( err );MSQ_ERRRTN( err );
    }

    MsqMeshEntity& elem                = elementArray[elem_index];
    const unsigned count               = elem.vertex_count();
    const size_t* const vertex_indices = elem.get_vertex_index_array();
    for( size_t i = 0; i < count; ++i )
    {
        const size_t v = vertex_indices[i];
        if( vertexNormalIndices.empty() )
        {
            normals_out[i] = normalData[v];
        }
        else if( vertexNormalIndices[v] < normalData.size() )
        {
            normals_out[i] = normalData[vertexNormalIndices[v]];
        }
        else
        {
            normals_out[i] = vertexArray[v];
            get_domain()->element_normal_at( elementHandlesArray[elem_index], normals_out[i] );
        }
    }
}

void PatchData::get_domain_normal_at_vertex( size_t vert_index,
                                             Mesh::EntityHandle handle,
                                             Vector3D& normal,
                                             MsqError& err )
{
    if( !domain_set() )
    {
        MSQ_SETERR( err )( "No domain constraint set.", MsqError::INVALID_STATE );
        return;
    }

    if( normalData.empty() )
    {
        update_cached_normals( err );MSQ_ERRRTN( err );
    }

    if( vertexNormalIndices.empty() )
    {
        normal = normalData[vert_index];
    }
    else if( vertexNormalIndices[vert_index] < normalData.size() )
    {
        normal = normalData[vertexNormalIndices[vert_index]];
    }
    else
    {
        normal = vertexArray[vert_index];
        get_domain()->element_normal_at( handle, normal );
    }
}

void PatchData::get_domain_normal_at_corner( size_t elem_index, unsigned corner, Vector3D& normal, MsqError& err )
{
    assert( elem_index < elementArray.size() );
    if( 2 != TopologyInfo::dimension( elementArray[elem_index].get_element_type() ) )
    {
        MSQ_SETERR( err )( "Attempt to get corners of non-surface element", MsqError::INVALID_ARG );
        return;
    }

    MsqMeshEntity& elem                = elementArray[elem_index];
    const size_t* const vertex_indices = elem.get_vertex_index_array();
    get_domain_normal_at_vertex( vertex_indices[corner], elementHandlesArray[elem_index], normal, err );MSQ_CHKERR( err );
}

void PatchData::set_mesh( Mesh* ms )
{
    myMesh = ms;
    // observers should treat this the same as if the
    // instance of this object wzs being deleted.
    notify_patch_destroyed();
}

void PatchData::set_domain( MeshDomain* d )
{
    myDomain = d;

    // clear all cached data from the previous domain
    vertexNormalIndices.clear();
    normalData.clear();
    // vertexDomainDOF.clear();

    // observers should treat this the same as if the
    // instance of this object wzs being deleted.
    notify_patch_destroyed();
}

static int width( double d )
{
    if( d == 0.0 ) return 1;
    const int max_precision = 6;
    int w                   = (int)ceil( log10( 0.001 + fabs( d ) ) );
    if( w < 0 ) w = 2 + std::min( max_precision, -w );
    if( d < 0.0 ) ++w;
    return w;
}
static int width( size_t t )
{
    return t ? (int)ceil( log10( (double)( 1 + t ) ) ) : 1;
}
static int width( const void* ptr )
{
    return width( (size_t)ptr );
}

ostream& operator<<( ostream& stream, const PatchData& pd )
{
    size_t i;
    int fw = 5;  // width of bit flags
    int hw = 6;  // width of a handle
    int cw = 4;  // with of a coordinate value
    int iw = 3;  // with of an index
    int tw = 3;  // with of the string name of an element type
    int xw = cw, yw = cw, zw = cw;

    for( i = 0; i < pd.num_nodes(); ++i )
    {
        int w = 2 + width( pd.vertexHandlesArray[i] );
        if( w > hw ) hw = w;
        w = width( pd.vertexArray[i].x() );
        if( w > xw ) xw = w;
        w = width( pd.vertexArray[i].y() );
        if( w > yw ) yw = w;
        w = width( pd.vertexArray[i].z() );
        if( w > zw ) zw = w;
    }
    for( i = 0; i < pd.num_elements(); ++i )
    {
        int w = 2 + width( pd.elementHandlesArray[i] );
        if( w > hw ) hw = w;
        const char* name = TopologyInfo::short_name( pd.elementArray[i].get_element_type() );
        if( name && (int)strlen( name ) > tw ) tw = strlen( name );
    }
    if( iw < (int)ceil( log10( (double)( 1 + pd.num_nodes() ) ) ) )
        iw = (int)ceil( log10( (double)( 1 + pd.num_nodes() ) ) );
    if( iw < (int)ceil( log10( (double)( 1 + pd.num_elements() ) ) ) )
        iw = (int)ceil( log10( (double)( 1 + pd.num_elements() ) ) );

    stream << "Vertices: " << endl;
    stream << "Flags: C: culled, F: fixed, S: slave, P: patch vertex, M: marked" << endl;
    stream << setw( iw ) << "Idx"
           << " " << setw( hw ) << "Handle"
           << " " << setw( cw ) << "X"
           << "," << setw( cw ) << "Y"
           << "," << setw( cw ) << "Z"
           << " " << setw( fw ) << "Flags"
           << " "
           << "Adj.Elems" << endl
           << setw( iw ) << setfill( '-' ) << ""
           << " " << setw( hw ) << setfill( '-' ) << ""
           << " " << setw( cw ) << setfill( '-' ) << ""
           << "," << setw( cw ) << setfill( '-' ) << ""
           << "," << setw( cw ) << setfill( '-' ) << ""
           << " " << setw( fw ) << setfill( '-' ) << ""
           << " " << setfill( ' ' ) << "-------------" << std::endl;
    for( i = 0; i < pd.num_nodes(); ++i )
    {
        stream << setw( iw ) << i << " " << setw( hw ) << pd.vertexHandlesArray[i] << " " << setw( cw )
               << pd.vertexArray[i].x() << "," << setw( cw ) << pd.vertexArray[i].y() << "," << setw( cw )
               << pd.vertexArray[i].z() << " ";
        if( pd.vertexArray[i].is_flag_set( MsqVertex::MSQ_CULLED ) )
            stream << "C";
        else
            stream << " ";
        if( pd.vertexArray[i].is_flag_set( MsqVertex::MSQ_HARD_FIXED ) )
            stream << "F";
        else
            stream << " ";
        if( pd.vertexArray[i].is_flag_set( MsqVertex::MSQ_DEPENDENT ) )
            stream << "S";
        else
            stream << " ";
        if( pd.vertexArray[i].is_flag_set( MsqVertex::MSQ_PATCH_FIXED ) )
            stream << "f";
        else
            stream << " ";
        if( pd.vertexArray[i].is_flag_set( MsqVertex::MSQ_MARK ) )
            stream << "M";
        else
            stream << " ";

        if( pd.vertAdjacencyArray.size() )
        {
            size_t j   = pd.vertAdjacencyOffsets[i];
            size_t end = pd.vertAdjacencyOffsets[i + 1];
            if( j != end ) stream << " " << pd.vertAdjacencyArray[j++];
            for( ; j < end; ++j )
                stream << "," << pd.vertAdjacencyArray[j];
        }

        stream << endl;
    }

    stream << "Elements: " << endl;
    stream << setw( iw ) << "Idx"
           << " " << setw( hw ) << "Handle"
           << " " << setw( tw + 2 ) << "Type"
           << " "
           << "Connectivity" << std::endl
           << setw( iw ) << setfill( '-' ) << ""
           << " " << setw( hw ) << setfill( '-' ) << ""
           << " " << setw( tw + 2 ) << setfill( '-' ) << ""
           << " " << setfill( ' ' ) << "--------------------------" << std::endl;
    for( i = 0; i < pd.num_elements(); ++i )
    {
        EntityTopology type = pd.elementArray[i].get_element_type();
        stream << setw( iw ) << i << " " << setw( hw ) << pd.elementHandlesArray[i] << " " << setw( tw )
               << TopologyInfo::short_name( type ) << left << setw( 2 ) << pd.elementArray[i].node_count() << internal
               << " " << setw( iw ) << pd.elementArray[i].get_vertex_index_array()[0];
        for( size_t j = 1; j < pd.elementArray[i].node_count(); ++j )
            stream << "," << setw( iw ) << pd.elementArray[i].get_vertex_index_array()[j];
        stream << endl;
    }
    stream << endl;

    stream << "Mesh: " << ( pd.myMesh ? "yes" : "no" ) << endl;
    stream << "Domain: " << ( pd.myDomain ? "yes" : "no" ) << endl;
    //   stream << "mType: " << (pd.mType==PatchData::VERTICES_ON_VERTEX_PATCH?"vert-on-vert":
    //                           pd.mType==PatchData::ELEMENTS_ON_VERTEX_PATCH?"elem-on-vert":
    //                           pd.mType==PatchData::GLOBAL_PATCH?"global":"unknown") << endl;

    if( pd.haveComputedInfos )
    {
        stream << "ComputedInfos:" << endl;
        if( pd.have_computed_info( PatchData::MIN_UNSIGNED_AREA ) )
            stream << "\t MIN_UNSINGED_AREA = " << pd.computedInfos[PatchData::MIN_UNSIGNED_AREA] << endl;
        if( pd.have_computed_info( PatchData::MAX_UNSIGNED_AREA ) )
            stream << "\t MAX_UNSIGNED_AREA = " << pd.computedInfos[PatchData::MAX_UNSIGNED_AREA] << endl;
        if( pd.have_computed_info( PatchData::MIN_EDGE_LENGTH ) )
            stream << "\t MIN_EDGE_LENGTH = " << pd.computedInfos[PatchData::MIN_EDGE_LENGTH] << endl;
        if( pd.have_computed_info( PatchData::MAX_EDGE_LENGTH ) )
            stream << "\t MAX_EDGE_LENGTH = " << pd.computedInfos[PatchData::MAX_EDGE_LENGTH] << endl;
        if( pd.have_computed_info( PatchData::MINMAX_SIGNED_DET2D ) )
            stream << "\t MINMAX_SIGNED_DET2D = " << pd.computedInfos[PatchData::MINMAX_SIGNED_DET2D] << endl;
        if( pd.have_computed_info( PatchData::MINMAX_SIGNED_DET3D ) )
            stream << "\t MINMAX_SIGNED_DET3D = " << pd.computedInfos[PatchData::MINMAX_SIGNED_DET3D] << endl;
        if( pd.have_computed_info( PatchData::AVERAGE_DET3D ) )
            stream << "\t AVERAGE_DET3D = " << pd.computedInfos[PatchData::AVERAGE_DET3D] << endl;
    }

    return stream << endl;
}

void print_patch_data( const PatchData& pd )
{
    std::cout << pd << std::endl;
}

void PatchData::enslave_higher_order_nodes( const size_t* elem_offset_array,
                                            unsigned char* vertex_flags,
                                            MsqError& ) const
{
    for( size_t i = 0; i < elementArray.size(); ++i )
    {
        size_t start    = elem_offset_array[i];
        size_t conn_len = elem_offset_array[i + 1] - start;
        for( size_t j = elementArray[i].vertex_count(); j < conn_len; ++j )
        {
            const size_t vert_idx = elemConnectivityArray[start + j];
            assert( vert_idx < vertexHandlesArray.size() );
            if( !( vertex_flags[vert_idx] & MsqVertex::MSQ_HARD_FIXED ) )
                vertex_flags[vert_idx] |= MsqVertex::MSQ_DEPENDENT;
        }
    }
}

void PatchData::initialize_data( size_t* elem_offset_array, unsigned char* vertex_flags, MsqError& )
{
    // Clear out data specific to patch
    vertexNormalIndices.clear();
    normalData.clear();
    // vertexDomainDOF.clear();

    // Clear any vertex->element adjacency data.  It
    // is probably invalid, and certainly will be by the time
    // this function completes if the mesh contains higher-order
    // elements.
    vertAdjacencyArray.clear();
    vertAdjacencyOffsets.clear();
    size_t i, j;
    for( i = 0; i < elementArray.size(); ++i )
    {
        size_t start    = elem_offset_array[i];
        size_t conn_len = elem_offset_array[i + 1] - start;
        assert( conn_len > 0 );
        elementArray[i].set_connectivity( &elemConnectivityArray[start], conn_len );
    }

    // Use vertAdjacencyOffsets array as temporary storage.
    vertAdjacencyOffsets.resize( vertexHandlesArray.size() + 1 );
    size_t* vertex_index_map = arrptr( vertAdjacencyOffsets );

    // Count number of free vertices and initialize vertex_index_map
    numFreeVertices = 0;
    for( i = 0; i < vertexHandlesArray.size(); ++i )
    {
        if( !( vertex_flags[i] & MsqVertex::MSQ_FIXED ) ) ++numFreeVertices;
        vertex_index_map[i] = i;
    }

    // Re-order vertices such that all free vertices are
    // first in the list.  Construct map from old to new
    // position in list for updating element connectivity.
    i = 0;
    j = numFreeVertices;
    for( ;; ++i, ++j )
    {
        // find next fixed vertex in the range [i,numFreeVertices)
        for( ; i < numFreeVertices && !( vertex_flags[i] & MsqVertex::MSQ_FIXED ); ++i )
            ;
        // if no more fixed vertices in the free vertex range [0,numFreeVertices)
        if( i == numFreeVertices ) break;
        // find the next free vertex in the range [j,num_nodes)
        for( ; vertex_flags[j] & MsqVertex::MSQ_FIXED; ++j )
            ;
        // swap fixed (i) and free (j) vertices
        vertex_index_map[i] = j;
        vertex_index_map[j] = i;
        std::swap( vertexHandlesArray[i], vertexHandlesArray[j] );
        std::swap( vertex_flags[i], vertex_flags[j] );
    }
    assert( i == numFreeVertices );
    assert( j <= vertexHandlesArray.size() );

    // Update element connectivity data for new vertex indices
    for( i = 0; i < elemConnectivityArray.size(); ++i )
        elemConnectivityArray[i] = vertex_index_map[elemConnectivityArray[i]];

    // Reorder vertices such that free, slave vertices
    // occur after free, non-slave vertices in list.
    numSlaveVertices = 0;
    for( i = 0; i < vertexHandlesArray.size(); ++i )
    {
        if( ( vertex_flags[i] & MsqVertex::MSQ_DEPENDENT ) && !( vertex_flags[i] & MsqVertex::MSQ_FIXED ) )
            ++numSlaveVertices;
    }
    numFreeVertices -= numSlaveVertices;

    if( numSlaveVertices )
    {
        // re-initialize vertex index map
        for( i = 0; i < vertexHandlesArray.size(); ++i )
            vertex_index_map[i] = i;

        // Re-order free vertices such that all slave vertices are
        // last in the list.  Construct map from old to new
        // position in list for updating element connectivity.
        i = 0;
        j = numFreeVertices;
        for( ;; ++i, ++j )
        {
            // find next slave vertex in the range [i,numFreeVertices)
            for( ; i < numFreeVertices && !( vertex_flags[i] & MsqVertex::MSQ_DEPENDENT ); ++i )
                ;
            // if no more slave vertices in [0,numFreeVertices), then done.
            if( i == numFreeVertices ) break;
            // find the next free (non-slave) vertex in the range
            //   [numFreeVertices,numFreeVertices+numSlaveVertices)
            for( ; vertex_flags[j] & MsqVertex::MSQ_DEPENDENT; ++j )
                ;
            // swap free (j) and slave (i) vertices
            vertex_index_map[i] = j;
            vertex_index_map[j] = i;
            std::swap( vertexHandlesArray[i], vertexHandlesArray[j] );
            std::swap( vertex_flags[i], vertex_flags[j] );
        }
        assert( i == numFreeVertices );
        assert( j <= numFreeVertices + numSlaveVertices );

        // Update element connectivity data for new vertex indices
        for( i = 0; i < elemConnectivityArray.size(); ++i )
            elemConnectivityArray[i] = vertex_index_map[elemConnectivityArray[i]];
    }

    // Clear temporary data
    vertAdjacencyOffsets.clear();

    notify_new_patch();
}

size_t PatchData::num_corners() const
{
    size_t result = 0;
    for( unsigned i = 0; i < elementArray.size(); ++i )
        result += elementArray[i].vertex_count();
    return result;
}

void PatchData::fill( size_t num_vertex,
                      const double* coords,
                      size_t num_elem,
                      EntityTopology type,
                      const size_t* connectivity,
                      const bool* fixed,
                      MsqError& err )
{
    std::vector< EntityTopology > types( num_elem );
    std::fill( types.begin(), types.end(), type );
    const EntityTopology* type_ptr = num_elem ? arrptr( types ) : 0;
    this->fill( num_vertex, coords, num_elem, type_ptr, connectivity, fixed, err );MSQ_CHKERR( err );
}

void PatchData::fill( size_t num_vertex,
                      const double* coords,
                      size_t num_elem,
                      const EntityTopology* types,
                      const size_t* conn,
                      const bool* fixed,
                      MsqError& err )
{
    std::vector< size_t > lengths( num_elem );
    std::transform( types, types + num_elem, lengths.begin(), std::ptr_fun( TopologyInfo::corners ) );
    const size_t* len_ptr = num_elem ? arrptr( lengths ) : 0;
    this->fill( num_vertex, coords, num_elem, types, len_ptr, conn, fixed, err );MSQ_CHKERR( err );
}

void PatchData::fill( size_t num_vertex,
                      const double* coords,
                      size_t num_elem,
                      const EntityTopology* types,
                      const size_t* lengths,
                      const size_t* conn,
                      const bool* fixed,
                      MsqError& err )
{
    size_t i;

    // count vertex uses
    size_t num_uses = std::accumulate( lengths, lengths + num_elem, 0 );

    // Allocate storage for data
    vertexArray.resize( num_vertex );
    vertexHandlesArray.resize( num_vertex );
    elementArray.resize( num_elem );
    elementHandlesArray.resize( num_elem );
    elemConnectivityArray.resize( num_uses );
    numFreeVertices  = 0;
    numSlaveVertices = 0;

    // Must call clear() first so that any stale values get
    // zero'd when we call resize.
    byteArray.clear();
    if( fixed )
    {
        byteArray.resize( num_vertex, 0 );
        for( i = 0; i < num_vertex; ++i )
            if( fixed[i] ) byteArray[i] |= ( MsqVertex::MSQ_HARD_FIXED | MsqVertex::MSQ_PATCH_FIXED );
    }

    for( i = 0; i < num_elem; ++i )
    {
        element_by_index( i ).set_element_type( types[i] );
        elementHandlesArray[i] = (Mesh::ElementHandle)i;
    }
    for( i = 0; i < num_vertex; ++i )
        vertexHandlesArray[i] = (Mesh::VertexHandle)i;

    memcpy( get_connectivity_array(), conn, num_uses * sizeof( size_t ) );

    std::vector< size_t > offsets( num_elem + 1 );
    size_t sum = offsets[0] = 0;
    for( i = 1; i <= num_elem; ++i )
        offsets[i] = sum += lengths[i - 1];

    const Settings::HigherOrderSlaveMode ho_mode =
        mSettings ? mSettings->get_slaved_ho_node_mode() : Settings::SLAVE_ALL;
    switch( ho_mode )
    {
        case Settings::SLAVE_ALL:
            byteArray.resize( num_vertex, 0 );
            enslave_higher_order_nodes( arrptr( offsets ), arrptr( byteArray ), err );MSQ_ERRRTN( err );
            break;
        case Settings::SLAVE_NONE:
            // Do nothing.  We clear other bits when processing the 'fixed' array above.
            break;
        default:
            MSQ_SETERR( err )
            ( "Specified higher-order noded slaving scheme not supported "
              "when initializind PatchData using PatchData::fill",
              MsqError::NOT_IMPLEMENTED );
            return;
    }

    this->initialize_data( arrptr( offsets ), arrptr( byteArray ), err );MSQ_ERRRTN( err );

    // initialize_data will re-order vertex handles and
    // update element connectivity accordingly.  Use
    // the values we stored in vertexHandlesArray to
    // figure out the new index of each vertex, and initialize
    // the vertex.
    for( i = 0; i < num_vertex; ++i )
        vertexArray[i] = coords + 3 * (size_t)vertexHandlesArray[i];

    for( i = 0; i < num_vertex; ++i )
        vertexArray[i].flags() = byteArray[i];
}

void PatchData::make_handles_unique( Mesh::EntityHandle* handles, size_t& count, size_t* index_map )
{
    if( count < 2 )
    {
        return;
    }
    // save this now, as we'll be changing count later
    const size_t* index_end = index_map + count;

    if( index_map )
    {
        // Copy input handle list into index map as a temporary
        // copy of the input list.
        assert( sizeof( Mesh::EntityHandle ) == sizeof( size_t ) );
        memcpy( index_map, handles, count * sizeof( size_t ) );
    }

    // Make handles a unique list
    std::sort( handles, handles + count );
    Mesh::EntityHandle* end = std::unique( handles, handles + count );
    count                   = end - handles;

    if( index_map )
    {
        // Replace each handle in index_map with the index of
        // its position in the handles array
        Mesh::EntityHandle* pos;
        for( size_t* iter = index_map; iter != index_end; ++iter )
        {
            pos   = std::lower_bound( handles, handles + count, (Mesh::EntityHandle)*iter );
            *iter = pos - handles;
        }
    }
}

void PatchData::fill_global_patch( MsqError& err )
{
    GlobalPatch gp;
    gp.set_mesh( get_mesh() );
    PatchIterator iter( &gp );
    bool b = iter.get_next_patch( *this, err );MSQ_ERRRTN( err );
    if( !b ) MSQ_SETERR( err )( "Empty patch in PatchData::fill_global_patch", MsqError::INVALID_MESH );
    assert( b );
}

void PatchData::set_mesh_entities( std::vector< Mesh::ElementHandle >& elements,
                                   std::vector< Mesh::VertexHandle >& free_vertices,
                                   MsqError& err )
{
    Mesh* current_mesh = get_mesh();
    if( !current_mesh )
    {
        MSQ_SETERR( err )( "No Mesh associated with PatchData.", MsqError::INVALID_STATE );
        return;
    }

    if( elements.empty() )
    {
        clear();
        return;
    }

    elementHandlesArray = elements;
    get_mesh()->elements_get_attached_vertices( arrptr( elementHandlesArray ), elementHandlesArray.size(),
                                                vertexHandlesArray, offsetArray, err );MSQ_ERRRTN( err );

    // Construct CSR-rep element connectivity data
    size_t num_vert = vertexHandlesArray.size();
    elemConnectivityArray.resize( num_vert );
    make_handles_unique( arrptr( vertexHandlesArray ), num_vert, arrptr( elemConnectivityArray ) );
    vertexHandlesArray.resize( num_vert );

    // Get element topologies
    std::vector< EntityTopology > elem_topologies( elementHandlesArray.size() );
    get_mesh()->elements_get_topologies( arrptr( elementHandlesArray ), arrptr( elem_topologies ),
                                         elementHandlesArray.size(), err );MSQ_ERRRTN( err );

    // get vertex bits from mesh
    byteArray.resize( vertexHandlesArray.size() );
    get_mesh()->vertices_get_byte( arrptr( vertexHandlesArray ), arrptr( byteArray ), vertexHandlesArray.size(), err );MSQ_ERRRTN( err );

    // if free_vertices is not empty, then we need to mark as
    // fixed any vertices *not* in that list.
    if( free_vertices.size() == 1 )
    {
        for( size_t i = 0; i < vertexHandlesArray.size(); ++i )
            if( vertexHandlesArray[i] == free_vertices.front() )
                byteArray[i] &= ~MsqVertex::MSQ_PATCH_FIXED;
            else
                byteArray[i] |= MsqVertex::MSQ_PATCH_FIXED;
    }
    else if( !free_vertices.empty() )
    {
        // sort and remove duplicates from free_vertices list.
        std::sort( free_vertices.begin(), free_vertices.end() );
        free_vertices.erase( std::unique( free_vertices.begin(), free_vertices.end() ), free_vertices.end() );

        for( size_t i = 0; i < vertexHandlesArray.size(); ++i )
        {
            if( ( byteArray[i] & MsqVertex::MSQ_DEPENDENT ) ||
                std::binary_search( free_vertices.begin(), free_vertices.end(), vertexHandlesArray[i] ) )
                byteArray[i] &= ~MsqVertex::MSQ_PATCH_FIXED;
            else
                byteArray[i] |= MsqVertex::MSQ_PATCH_FIXED;
        }
    }

    // set element types
    elementArray.resize( elementHandlesArray.size() );
    for( size_t i = 0; i < elementHandlesArray.size(); ++i )
        elementArray[i].set_element_type( elem_topologies[i] );

    // get element connectivity, group vertices by free/slave/fixed state
    initialize_data( arrptr( offsetArray ), arrptr( byteArray ), err );MSQ_ERRRTN( err );

    // get vertex coordinates
    vertexArray.resize( vertexHandlesArray.size() );
    get_mesh()->vertices_get_coordinates( arrptr( vertexHandlesArray ), arrptr( vertexArray ),
                                          vertexHandlesArray.size(), err );MSQ_ERRRTN( err );

    // set vertex flags
    for( size_t i = 0; i < vertexArray.size(); ++i )
        vertexArray[i].flags() = byteArray[i];
}

void PatchData::get_sample_location( size_t element_index, Sample sample, Vector3D& result, MsqError& err ) const
{
    const MsqMeshEntity& elem      = element_by_index( element_index );
    const NodeSet ho_bits          = non_slave_node_set( element_index );
    const MappingFunction* const f = get_mapping_function( elem.get_element_type() );
    if( !f )
    {
        MSQ_SETERR( err )( "No mapping function", MsqError::UNSUPPORTED_ELEMENT );
        return;
    }

    double coeff[27];
    size_t num_coeff, index[27];
    f->coefficients( sample, ho_bits, coeff, index, num_coeff, err );MSQ_ERRRTN( err );
    f->convert_connectivity_indices( elem.node_count(), index, num_coeff, err );MSQ_ERRRTN( err );

    const size_t* const conn = elem.get_vertex_index_array();
    assert( num_coeff > 0 );
    result = coeff[0] * vertex_by_index( conn[index[0]] );
    for( unsigned i = 1; i < num_coeff; ++i )
        result += coeff[i] * vertex_by_index( conn[index[i]] );
}

NodeSet PatchData::non_slave_node_set( size_t element_index ) const
{
    const MsqMeshEntity& elem = element_by_index( element_index );
    EntityTopology type       = elem.get_element_type();
    const size_t* conn        = elem.get_vertex_index_array();
    const size_t n            = elem.node_count();

    MsqError err;
    bool have_midedge, have_midface, have_midelem;
    unsigned num_edge = 0, num_face = 0, num_corner = TopologyInfo::corners( type );
    TopologyInfo::higher_order( type, n, have_midedge, have_midface, have_midelem, err );
    num_edge = TopologyInfo::edges( type );
    if( TopologyInfo::dimension( type ) == 2 )
        num_face = 1;
    else
        num_face = TopologyInfo::faces( type );

    NodeSet result;
    result.set_all_corner_nodes( type );
    if( have_midedge )
    {
        for( unsigned i = 0; i < num_edge; ++i )
        {
            if( !( vertex_by_index( conn[num_corner + i] ).get_flags() & MsqVertex::MSQ_DEPENDENT ) )
                result.set_mid_edge_node( i );
        }
    }
    if( have_midface )
    {
        for( unsigned i = 0; i < num_face; ++i )
        {
            if( !( vertex_by_index( conn[num_corner + num_edge + i] ).get_flags() & MsqVertex::MSQ_DEPENDENT ) )
                result.set_mid_face_node( i );
        }
    }
    if( have_midelem &&
        !( vertex_by_index( conn[num_corner + num_edge + num_face] ).get_flags() & MsqVertex::MSQ_DEPENDENT ) )
        result.set_mid_region_node();

    return result;
}

void PatchData::get_samples( size_t element, std::vector< Sample >& samples, MsqError& ) const
{
    NodeSet ns = get_samples( element );
    samples.resize( ns.num_nodes() );
    std::vector< Sample >::iterator i = samples.begin();

    unsigned j;
    EntityTopology type = element_by_index( element ).get_element_type();
    for( j = 0; j < TopologyInfo::corners( type ); ++j )
        if( ns.corner_node( j ) ) *( i++ ) = Sample( 0, j );
    for( j = 0; j < TopologyInfo::edges( type ); ++j )
        if( ns.mid_edge_node( j ) ) *( i++ ) = Sample( 1, j );
    if( TopologyInfo::dimension( type ) == 3 )
    {
        for( j = 0; j < TopologyInfo::faces( type ); ++j )
            if( ns.mid_face_node( j ) ) *( i++ ) = Sample( 2, j );
        if( ns.mid_region_node() ) *( i++ ) = Sample( 3, 0 );
    }
    else if( ns.mid_face_node( 0 ) )
        *( i++ ) = Sample( 2, 0 );

    assert( i == samples.end() );
}

bool PatchData::attach_extra_data( ExtraData* data )
{
    if( data->patchNext )
    {
        return false;
    }

    if( !data->patchPtr )
        data->patchPtr = this;
    else if( data->patchPtr != this )
        return false;

    data->patchNext = dataList;
    dataList        = data;
    return true;
}

bool PatchData::remove_extra_data( ExtraData* data )
{
    if( data->patchPtr != this ) return false;

    if( dataList == data )
    {
        dataList        = data->patchNext;
        data->patchNext = 0;
        data->patchPtr  = 0;
        return true;
    }

    for( ExtraData* iter = dataList; iter; iter = iter->patchNext )
        if( iter->patchNext == data )
        {
            iter->patchNext = data->patchNext;
            data->patchNext = 0;
            data->patchPtr  = 0;
            return true;
        }

    return false;
}

void PatchData::notify_new_patch()
{
    for( ExtraData* iter = dataList; iter; iter = iter->patchNext )
        iter->notify_new_patch();
}

void PatchData::notify_sub_patch( PatchData& sub_patch,
                                  const size_t* vertex_map,
                                  const size_t* element_map,
                                  MsqError& err )
{
    for( ExtraData* iter = dataList; iter; iter = iter->patchNext )
    {
        iter->notify_sub_patch( sub_patch, vertex_map, element_map, err );MSQ_ERRRTN( err );
    }
}

void PatchData::notify_patch_destroyed()
{
    // Remove all ExtraDatas from list and notify them
    // that they are being destroyed.
    while( dataList )
    {
        ExtraData* dead_data = dataList;
        dataList             = dataList->patchNext;
        dead_data->patchNext = 0;
        dead_data->patchPtr  = 0;
        dead_data->notify_patch_destroyed();
    }
}

}  // namespace MBMesquite