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