MOAB: Mesh Oriented datABase  (version 5.2.1)
DomainClassifier.cpp
Go to the documentation of this file.
00001 /* *****************************************************************
00002     MESQUITE -- The Mesh Quality Improvement Toolkit
00003 
00004     Copyright 2008 Lawrence Livermore National Laboratory.  Under
00005     the terms of Contract B545069 with the University of Wisconsin --
00006     Madison, Lawrence Livermore National Laboratory 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     (2008) kraftche@cae.wisc.edu
00024 
00025   ***************************************************************** */
00026 
00027 /** \file DomainClassifier.cpp
00028  *  \brief
00029  *  \author Jason Kraftcheck
00030  */
00031 
00032 #include "DomainClassifier.hpp"
00033 #include "MsqError.hpp"
00034 #include "TopologyInfo.hpp"
00035 #include "MsqVertex.hpp"
00036 #include <map>
00037 #include <set>
00038 #include <algorithm>
00039 #include <iterator>
00040 
00041 namespace MBMesquite
00042 {
00043 
00044 /**\brief Helper function for 'find_skin'
00045  *
00046  * Check the vertices on element sides to see if the
00047  * elements are adjacent.
00048  *\param num_vtx  Number of vertices in each element side
00049  *\param idx_set_1 Indices into the first element's vertex list
00050  *                 at which the side vertices are located.
00051  *\param vtx_set_1 Vertex of first element.
00052  *\param idx_set_2 Indices into the second element's vertex list
00053  *                 at which the side vertices are located.
00054  *\param vtx_set_2 Vertex of second element.
00055  */
00056 static bool compare_sides( unsigned num_vtx, const unsigned* idx_set_1, const Mesh::VertexHandle* vtx_set_1,
00057                            const unsigned* idx_set_2, const Mesh::VertexHandle* vtx_set_2 );
00058 
00059 /**\brief Skin mesh
00060  *
00061  * Find vertices and lower-dimension elements on the boundary
00062  * of a mesh.
00063  */
00064 static void find_skin( Mesh* mesh, std::vector< Mesh::VertexHandle >& skin_verts,
00065                        std::vector< Mesh::ElementHandle >& skin_elems, MsqError& err );
00066 
00067 /**\brief Group MeshDomain objectects by topological dimension.
00068  *
00069  * Given an array of MeshDomain pointers and and array of dimensions,
00070  * where each dimension value is the topological dimension of the
00071  * corresponding domain, construct a DomainSet vector containing
00072  * the domains, grouped from lowest dimension to highest.
00073  *\param domains  Array of domains
00074  *\param dims     Topological dimension of each domain
00075  *\param count    Length of 'domains' and 'dims' arrays
00076  *\param results  Output
00077  *\param dim_indices Output: index in output at which the first
00078  *                domain of the corresponding topological dimension
00079  *                occurs.
00080  */
00081 static void dimension_sort_domains( MeshDomain** domains, const int* dims, unsigned count,
00082                                     std::vector< DomainClassifier::DomainSet >& results, int dim_indices[4],
00083                                     MsqError& err );
00084 
00085 /**\brief Classify mesh vertices using vertex position
00086  *
00087  * Classify vertices by distance from each geometric domain.
00088  * A vertex is assigned to a domain if it lies within the domain.
00089  * If a vertex lies in multiple domains, it is assigned to the
00090  * domain with the lowest topological dimension.
00091  *\param mesh  The MBMesquite::Mesh instance
00092  *\param dim_sorted_domains  The list of MeshDomains, sorted from lowest
00093  *             to highest dimension.
00094  *\param num_domain  Length of 'dim_sorted_domains' array.
00095  *\param vertices    Vertices to classify
00096  *\param epsilon     Maximum distance a vertex may deviate from its domain.
00097  */
00098 static void geom_classify_vertices( Mesh* mesh, DomainClassifier::DomainSet* dim_sorted_domains, unsigned num_domain,
00099                                     const std::vector< Mesh::VertexHandle >& vertices, double epsilon, MsqError& err );
00100 
00101 /**\brief Classify mesh elements using vertex classification
00102  *
00103  * Classify elements using vertex classification, according to the
00104  * following rules, listed highest-precedence first.
00105  * Case 1: one or more vertices has no classification -> element has no classification
00106  * Case 2: all vertices classified to domain with dimension < 2 -> unknown
00107  * Case 3: one or more vertices has dimension-2 classification:
00108  *  Case 3a: all have same dimension-2 domain -> element assigned to same domain
00109  *  Case 3b: different dimension-2 domains -> element has no classification
00110  *
00111  *\param mesh  The MBMesquite::Mesh instance
00112  *\param dim_sorted_domains  The list of MeshDomains, sorted from lowest
00113  *             to highest dimension.
00114  *\param num_domain  Length of 'dim_sorted_domains' array.
00115  *\param elems    Elements to classify
00116  *\param uknown   Elements for which all vertices are classified to
00117  *                domains with dimension less than two.
00118  */
00119 static void vert_classify_elements( Mesh* mesh, DomainClassifier::DomainSet* dim_sorted_domains, unsigned num_domain,
00120                                     int dim_indices[4], const std::vector< Mesh::ElementHandle >& elems,
00121                                     std::vector< Mesh::ElementHandle >& unknown_elems, MsqError& err );
00122 
00123 /**\brief Classify mesh elements geometrically. */
00124 static void geom_classify_elements( Mesh* mesh, DomainClassifier::DomainSet* dim_sorted_domains, unsigned num_domain,
00125                                     int dim_indices[4], const std::vector< Mesh::ElementHandle >& elems,
00126                                     std::vector< Mesh::ElementHandle >& unknown_elems, double epsilon, MsqError& err );
00127 
00128 static bool compare_sides( unsigned num_vtx, const unsigned* idx_set_1, const Mesh::VertexHandle* vtx_set_1,
00129                            const unsigned* idx_set_2, const Mesh::VertexHandle* vtx_set_2 )
00130 {
00131     unsigned i;
00132     for( i = 0;; ++i )
00133     {
00134         if( i >= num_vtx ) return false;
00135         if( vtx_set_1[idx_set_1[0]] == vtx_set_2[idx_set_2[i]] ) break;
00136     }
00137 
00138     if( num_vtx == 1 ) return true;
00139 
00140     if( vtx_set_1[idx_set_1[1]] == vtx_set_2[idx_set_2[( i + 1 ) % num_vtx]] )
00141     {
00142         for( unsigned j = 2; j < num_vtx; ++j )
00143             if( vtx_set_1[idx_set_1[j]] != vtx_set_2[idx_set_2[( i + j ) % num_vtx]] ) return false;
00144         return true;
00145     }
00146     else if( vtx_set_1[idx_set_1[1]] == vtx_set_2[idx_set_2[( i + num_vtx - 1 ) % num_vtx]] )
00147     {
00148         for( unsigned j = 2; j < num_vtx; ++j )
00149             if( vtx_set_1[idx_set_1[j]] != vtx_set_2[idx_set_2[( i + num_vtx - j ) % num_vtx]] ) return false;
00150         return true;
00151     }
00152     else
00153     {
00154         return false;
00155     }
00156 }
00157 
00158 static void find_skin( Mesh* mesh, std::vector< Mesh::VertexHandle >& skin_verts,
00159                        std::vector< Mesh::ElementHandle >& skin_elems, MsqError& err )
00160 {
00161     std::vector< Mesh::ElementHandle > elements, adj_elem, side_elem;
00162     std::vector< Mesh::VertexHandle > vertices, adj_vtx;
00163     std::vector< size_t > junk;
00164     mesh->get_all_elements( elements, err );MSQ_ERRRTN( err );
00165     if( elements.empty() ) return;
00166     const unsigned elem_idx[] = { 0, 1, 2, 3 };
00167 
00168     std::vector< unsigned > extra_side_vtx;
00169     for( size_t e = 0; e < elements.size(); ++e )
00170     {
00171         Mesh::ElementHandle elem = elements[e];
00172         vertices.clear();
00173         mesh->elements_get_attached_vertices( &elem, 1, vertices, junk, err );MSQ_ERRRTN( err );
00174 
00175         EntityTopology type;
00176         mesh->elements_get_topologies( &elem, &type, 1, err );MSQ_ERRRTN( err );
00177         const unsigned dim    = TopologyInfo::dimension( type );
00178         const unsigned nsides = TopologyInfo::sides( type );
00179         for( unsigned s = 0; s < nsides; ++s )
00180         {
00181             unsigned num_side_vtx;
00182             const unsigned* side_vtx = TopologyInfo::side_vertices( type, dim - 1, s, num_side_vtx );
00183             adj_elem.clear();
00184             side_elem.clear();
00185             mesh->vertices_get_attached_elements( &vertices[*side_vtx], 1, adj_elem, junk, err );MSQ_ERRRTN( err );
00186 
00187             bool found_adj = false;
00188             for( size_t a = 0; a < adj_elem.size(); ++a )
00189             {
00190                 if( adj_elem[a] == elem ) continue;
00191                 EntityTopology adj_type;
00192                 mesh->elements_get_topologies( &adj_elem[a], &adj_type, 1, err );MSQ_ERRRTN( err );
00193                 adj_vtx.clear();
00194                 mesh->elements_get_attached_vertices( &adj_elem[a], 1, adj_vtx, junk, err );MSQ_ERRRTN( err );
00195                 if( TopologyInfo::dimension( adj_type ) == dim )
00196                 {
00197                     unsigned nsides2 = TopologyInfo::sides( adj_type );
00198                     for( unsigned s2 = 0; s2 < nsides2; ++s2 )
00199                     {
00200                         unsigned num_side_vtx2;
00201                         const unsigned* side_vtx2 = TopologyInfo::side_vertices( adj_type, dim - 1, s2, num_side_vtx2 );
00202                         if( num_side_vtx2 == num_side_vtx &&
00203                             compare_sides( num_side_vtx, side_vtx, arrptr( vertices ), side_vtx2, arrptr( adj_vtx ) ) )
00204                         {
00205                             found_adj = true;
00206                             break;
00207                         }
00208                     }
00209                 }
00210                 else
00211                 {
00212                     if( TopologyInfo::corners( adj_type ) == num_side_vtx )
00213                         if( compare_sides( num_side_vtx, side_vtx, arrptr( vertices ), elem_idx, arrptr( adj_vtx ) ) )
00214                             side_elem.push_back( adj_elem[a] );
00215                 }
00216             }
00217 
00218             if( !found_adj )
00219             {
00220                 int ho = TopologyInfo::higher_order( type, vertices.size(), err );MSQ_ERRRTN( err );
00221                 // mask mid-element node (because that cannot be on the skin)
00222                 ho &= ~( 1 << dim );
00223                 // if element has other higher-order nodes, we need to find the subset on the skin
00224                 // side
00225                 if( ho )
00226                 {
00227                     extra_side_vtx.clear();
00228                     extra_side_vtx.insert( extra_side_vtx.end(), side_vtx, side_vtx + num_side_vtx );
00229                     // get mid-side node
00230                     if( ho & ( 1 << ( dim - 1 ) ) )
00231                     {
00232                         int idx = TopologyInfo::higher_order_from_side( type, vertices.size(), dim - 1, s, err );MSQ_ERRRTN( err );
00233                         extra_side_vtx.push_back( idx );
00234                     }
00235                     if( dim == 3 && ( ho & 2 ) )
00236                     {  // if volume element, need mid-edge nodes too
00237                         const unsigned nedges = TopologyInfo::edges( type );
00238                         for( unsigned je = 0; je < nedges; ++je )
00239                         {
00240                             unsigned tmp;
00241                             const unsigned* edge_vtx = TopologyInfo::side_vertices( type, 1, je, tmp );
00242                             assert( edge_vtx && 2 == tmp );
00243                             unsigned idx = std::find( side_vtx, side_vtx + num_side_vtx, edge_vtx[0] ) - side_vtx;
00244                             if( idx < num_side_vtx )
00245                             {
00246                                 if( ( edge_vtx[1] == side_vtx[( idx + 1 ) % num_side_vtx] ) ||
00247                                     ( edge_vtx[1] == side_vtx[( idx + num_side_vtx - 1 ) % num_side_vtx] ) )
00248                                 {
00249                                     idx = TopologyInfo::higher_order_from_side( type, vertices.size(), 1, je, err );MSQ_ERRRTN( err );
00250                                     extra_side_vtx.push_back( idx );
00251                                 }
00252                             }
00253                         }
00254                     }
00255 
00256                     num_side_vtx = extra_side_vtx.size();
00257                     side_vtx     = arrptr( extra_side_vtx );
00258                 }
00259 
00260                 for( unsigned v = 0; v < num_side_vtx; ++v )
00261                     skin_verts.push_back( vertices[side_vtx[v]] );
00262                 for( unsigned j = 0; j < side_elem.size(); ++j )
00263                     skin_elems.push_back( side_elem[j] );
00264             }
00265         }
00266     }
00267 
00268     std::sort( skin_verts.begin(), skin_verts.end() );
00269     skin_verts.erase( std::unique( skin_verts.begin(), skin_verts.end() ), skin_verts.end() );
00270     std::sort( skin_elems.begin(), skin_elems.end() );
00271     skin_elems.erase( std::unique( skin_elems.begin(), skin_elems.end() ), skin_elems.end() );
00272 }
00273 
00274 static void dimension_sort_domains( MeshDomain** domains, const int* dims, unsigned count,
00275                                     std::vector< DomainClassifier::DomainSet >& results, int dim_indices[4],
00276                                     MsqError& err )
00277 {
00278     results.clear();
00279     results.resize( count );
00280     unsigned idx = 0;
00281     int dim;
00282     for( dim = 0; idx < count; ++dim )
00283     {
00284         if( dim > 3 )
00285         {
00286             MSQ_SETERR( err )( "Invalid domain dimension", MsqError::INVALID_ARG );
00287             return;
00288         }
00289         dim_indices[dim] = idx;
00290         for( unsigned j = 0; j < count; ++j )
00291         {
00292             if( dims[j] == dim )
00293             {
00294                 results[idx].domain = domains[j];
00295                 ++idx;
00296             }
00297         }
00298     }
00299     for( ; dim <= 3; ++dim )
00300         dim_indices[dim] = idx;
00301 }
00302 
00303 static void geom_classify_vertices( Mesh* mesh, DomainClassifier::DomainSet* dim_sorted_domains, unsigned num_domain,
00304                                     const std::vector< Mesh::VertexHandle >& vertices, double epsilon, MsqError& err )
00305 {
00306     Vector3D pt;
00307     MsqVertex coord;
00308     unsigned i;
00309     unsigned short dim;
00310     const double epssqr = epsilon * epsilon;
00311     std::vector< Mesh::VertexHandle >::const_iterator iter;
00312     for( iter = vertices.begin(); iter != vertices.end(); ++iter )
00313     {
00314         mesh->vertices_get_coordinates( &*iter, &coord, 1, err );MSQ_ERRRTN( err );
00315         for( i = 0; i < num_domain; ++i )
00316         {
00317             // Call DoF function first just to see if it fails
00318             // (e.g. if domain knows which vertices it is responsible for,
00319             //  it should fail if this vertex isn't one.)
00320             dim_sorted_domains[i].domain->domain_DoF( &*iter, &dim, 1, err );
00321             if( err || dim > 2 )
00322             {
00323                 err.clear();
00324                 continue;
00325             }
00326             pt = coord;
00327             dim_sorted_domains[i].domain->snap_to( *iter, pt );
00328             if( ( pt - coord ).length_squared() <= epssqr )
00329             {
00330                 dim_sorted_domains[i].vertices.push_back( *iter );
00331                 break;
00332             }
00333         }
00334     }
00335 }
00336 
00337 static void vert_classify_elements( Mesh* mesh, DomainClassifier::DomainSet* dim_sorted_domains, unsigned num_domain,
00338                                     int dim_indices[4], const std::vector< Mesh::ElementHandle >& elems,
00339                                     std::vector< Mesh::ElementHandle >& unknown_elems, MsqError& err )
00340 {
00341     // sort vertex lists for faster search
00342     for( unsigned i = 0; i < num_domain; ++i )
00343         std::sort( dim_sorted_domains[i].vertices.begin(), dim_sorted_domains[i].vertices.end() );
00344 
00345     std::vector< Mesh::ElementHandle >::const_iterator iter;
00346     std::vector< Mesh::VertexHandle > verts;
00347     std::vector< size_t > junk;
00348     for( iter = elems.begin(); iter != elems.end(); ++iter )
00349     {
00350         int dom = -1;
00351         verts.clear();
00352         mesh->elements_get_attached_vertices( &*iter, 1, verts, junk, err );MSQ_ERRRTN( err );
00353         for( unsigned v = 0; v < verts.size(); ++v )
00354         {
00355             int i;
00356             for( i = 0; i < dim_indices[3]; ++i )
00357             {
00358                 DomainClassifier::DomainSet* d = dim_sorted_domains + i;
00359                 if( std::binary_search( d->vertices.begin(), d->vertices.end(), verts[v] ) ) break;
00360             }
00361             // if any vertex in element has no domain, then element has no domain
00362             if( i < 0 )
00363             {
00364                 dom = -2;
00365                 break;
00366             }
00367             // if vertex is in curve or point, ignore it
00368             else if( i < dim_indices[2] )
00369                 continue;
00370             // else if we already have a vertex on a different surface,
00371             // element must be in volume (no domain)
00372             else if( dom >= 0 && dom != i )
00373             {
00374                 dom = -2;
00375                 break;
00376             }
00377             else
00378             {
00379                 dom = i;
00380             }
00381         }
00382 
00383         if( dom >= 0 )
00384             dim_sorted_domains[dom].elements.push_back( *iter );
00385         else if( dom == -1 )  // all verts on curves or points
00386             unknown_elems.push_back( *iter );
00387     }
00388 }
00389 
00390 static void geom_classify_elements( Mesh* mesh, DomainClassifier::DomainSet* dim_sorted_domains, unsigned,
00391                                     int dim_indices[4], const std::vector< Mesh::ElementHandle >& elems,
00392                                     std::vector< Mesh::ElementHandle >& unknown_elems, double epsilon, MsqError& err )
00393 {
00394     if( elems.empty() ) return;
00395 
00396     const double epssqr = epsilon * epsilon;
00397     Vector3D pt;
00398     MsqVertex coord;
00399     std::vector< Mesh::ElementHandle >::const_iterator iter;
00400     std::vector< Mesh::VertexHandle > verts;
00401     std::vector< MsqVertex > coords;
00402     std::vector< size_t > junk;
00403     for( iter = elems.begin(); iter != elems.end(); ++iter )
00404     {
00405         verts.clear();
00406         mesh->elements_get_attached_vertices( &*iter, 1, verts, junk, err );MSQ_ERRRTN( err );
00407         coords.resize( verts.size() );
00408         mesh->vertices_get_coordinates( arrptr( verts ), arrptr( coords ), verts.size(), err );MSQ_ERRRTN( err );
00409 
00410         // get all 2D domains that contain first vertex
00411         std::vector< int > doms;
00412         for( int i = dim_indices[2]; i < dim_indices[3]; ++i )
00413         {
00414             pt = coords[0];
00415             dim_sorted_domains[i].domain->snap_to( verts[0], pt );
00416             if( ( pt - coords[0] ).length_squared() <= epssqr ) doms.push_back( i );
00417         }
00418 
00419         // of the 2D domains containing the first vertex,
00420         // remove any that do not contain all other corners of the element.
00421         for( unsigned i = 1; i < verts.size(); ++i )
00422         {
00423             for( unsigned j = 0; j < doms.size(); )
00424             {
00425                 pt = coords[i];
00426                 dim_sorted_domains[doms[j]].domain->snap_to( verts[j], pt );
00427                 if( ( pt - coords[0] ).length_squared() <= epssqr )
00428                     ++j;
00429                 else
00430                     doms.erase( doms.begin() + j );
00431             }
00432         }
00433 
00434         if( doms.size() == 1 )
00435             dim_sorted_domains[doms.front()].elements.push_back( *iter );
00436         else
00437             unknown_elems.push_back( *iter );
00438     }
00439 }
00440 
00441 void DomainClassifier::classify_geometrically( DomainClassifier& result, Mesh* mesh, double tolerance,
00442                                                MeshDomain** domain_array, const int* dimension_array,
00443                                                unsigned array_length, MsqError& err )
00444 {
00445     if( !array_length ) return;
00446 
00447     std::vector< DomainSet > domains;
00448     int dim_indices[4];
00449     dimension_sort_domains( domain_array, dimension_array, array_length, domains, dim_indices, err );MSQ_ERRRTN( err );
00450 
00451     // classify vertices by position
00452     std::vector< Mesh::VertexHandle > vertices;
00453     mesh->get_all_vertices( vertices, err );MSQ_ERRRTN( err );
00454     geom_classify_vertices( mesh, arrptr( domains ), dim_indices[3], vertices, tolerance, err );MSQ_ERRRTN( err );
00455 
00456     // get elements and types
00457     std::vector< Mesh::ElementHandle > elems;
00458     mesh->get_all_elements( elems, err );MSQ_ERRRTN( err );
00459     if( elems.empty() ) return;
00460     std::vector< EntityTopology > types( elems.size() );
00461     mesh->elements_get_topologies( arrptr( elems ), arrptr( types ), elems.size(), err );MSQ_ERRRTN( err );
00462 
00463     // get rid of elements w/ dimension other than 2
00464     size_t w = 0;
00465     for( size_t r = 0; r < elems.size(); ++r )
00466     {
00467         if( TopologyInfo::dimension( types[r] ) == 2 )
00468         {
00469             elems[w] = elems[r];
00470             ++w;
00471         }
00472     }
00473     elems.resize( w );
00474 
00475     // classify elements using vertex classification
00476     std::vector< Mesh::ElementHandle > unknown;
00477     vert_classify_elements( mesh, arrptr( domains ), domains.size(), dim_indices, elems, unknown, err );MSQ_ERRRTN( err );
00478 
00479     // classify unknown elements geometrically
00480     elems.swap( unknown );
00481     unknown.clear();
00482     geom_classify_elements( mesh, arrptr( domains ), domains.size(), dim_indices, elems, unknown, tolerance, err );
00483 
00484     classify_by_handle( result, mesh, arrptr( domains ), domains.size(), err );
00485 }
00486 
00487 void DomainClassifier::classify_skin_geometrically( DomainClassifier& result, Mesh* mesh, double tolerance,
00488                                                     MeshDomain** domain_array, const int* dimension_array,
00489                                                     unsigned array_length, MsqError& err )
00490 {
00491     if( !array_length ) return;
00492 
00493     std::vector< DomainSet > domains;
00494     int dim_indices[4];
00495     dimension_sort_domains( domain_array, dimension_array, array_length, domains, dim_indices, err );MSQ_ERRRTN( err );
00496 
00497     std::vector< Mesh::VertexHandle > vertices;
00498     std::vector< Mesh::ElementHandle > elements;
00499     find_skin( mesh, vertices, elements, err );MSQ_ERRRTN( err );
00500 
00501     // classify vertices by position
00502     geom_classify_vertices( mesh, arrptr( domains ), dim_indices[3], vertices, tolerance, err );MSQ_ERRRTN( err );
00503 
00504     // classify elements using vertex classification
00505     std::vector< Mesh::ElementHandle > unknown;
00506     vert_classify_elements( mesh, arrptr( domains ), domains.size(), dim_indices, elements, unknown, err );MSQ_ERRRTN( err );
00507 
00508     // classify unknown elements geometrically
00509     elements.swap( unknown );
00510     unknown.clear();
00511     geom_classify_elements( mesh, arrptr( domains ), domains.size(), dim_indices, elements, unknown, tolerance, err );
00512 
00513     if( !unknown.empty() )
00514     {
00515         MSQ_SETERR( err )( "Count not classify all skin elements", MsqError::INVALID_MESH );
00516         return;
00517     }
00518 
00519     classify_by_handle( result, mesh, arrptr( domains ), domains.size(), err );
00520 }
00521 
00522 void DomainClassifier::classify_by_tag( DomainClassifier& result, Mesh* mesh, const char* tag_name,
00523                                         MeshDomain** domain_array, const int* id_array, unsigned array_length,
00524                                         MsqError& err )
00525 {
00526     TagHandle tag = mesh->tag_get( tag_name, err );MSQ_ERRRTN( err );
00527 
00528     std::string name2;
00529     Mesh::TagType type;
00530     unsigned size;
00531     mesh->tag_properties( tag, name2, type, size, err );MSQ_ERRRTN( err );
00532     if( type != Mesh::INT || size != 1 )
00533     {
00534         MSQ_SETERR( err )
00535         ( MsqError::TAG_NOT_FOUND, "Tag does not contain single integer value: %s\n", tag_name );
00536         return;
00537     }
00538 
00539     std::vector< DomainSet > sets( array_length );
00540     std::map< int, DomainSet* > idmap;
00541     for( unsigned i = 0; i < array_length; ++i )
00542     {
00543         sets[i].domain     = domain_array[i];
00544         idmap[id_array[i]] = &sets[i];
00545     }
00546 
00547     std::vector< Mesh::VertexHandle > vertices;
00548     std::vector< Mesh::ElementHandle > elements;
00549     mesh->get_all_vertices( vertices, err );MSQ_ERRRTN( err );
00550     mesh->get_all_elements( elements, err );MSQ_ERRRTN( err );
00551     if( vertices.empty() ) return;
00552 
00553     std::map< int, DomainSet* >::const_iterator iter;
00554     std::vector< int > ids( vertices.size() );
00555     mesh->tag_get_vertex_data( tag, vertices.size(), arrptr( vertices ), arrptr( ids ), err );MSQ_ERRRTN( err );
00556     for( size_t j = 0; j < vertices.size(); ++j )
00557     {
00558         iter = idmap.find( ids[j] );
00559         if( iter != idmap.end() ) iter->second->vertices.push_back( vertices[j] );
00560     }
00561 
00562     ids.clear();
00563     ids.resize( elements.size() );
00564     if( !elements.empty() )
00565     {
00566         mesh->tag_get_element_data( tag, elements.size(), arrptr( elements ), arrptr( ids ), err );
00567         if( err ) { err.clear(); }
00568         else
00569         {
00570             for( size_t k = 0; k < elements.size(); ++k )
00571             {
00572                 iter = idmap.find( ids[k] );
00573                 if( iter != idmap.end() ) iter->second->elements.push_back( elements[k] );
00574             }
00575         }
00576     }
00577 
00578     if( !sets.empty() ) classify_by_handle( result, mesh, arrptr( sets ), sets.size(), err );
00579 }
00580 
00581 void DomainClassifier::classify_by_handle( DomainClassifier& result, Mesh* mesh, DomainSet* sets, unsigned array_length,
00582                                            MsqError& err )
00583 {
00584     result.clear();
00585 
00586     // get all vertices and elements
00587     std::vector< Mesh::VertexHandle > vertices;
00588     std::vector< Mesh::ElementHandle > elements;
00589     mesh->get_all_vertices( vertices, err );MSQ_ERRRTN( err );
00590     mesh->get_all_elements( elements, err );MSQ_ERRRTN( err );
00591 
00592     // sort all arrays so we can merge
00593     std::sort( vertices.begin(), vertices.end() );
00594     std::sort( elements.begin(), elements.end() );
00595     for( unsigned i = 0; i < array_length; ++i )
00596     {
00597         std::sort( sets[i].vertices.begin(), sets[i].vertices.end() );
00598         std::sort( sets[i].elements.begin(), sets[i].elements.end() );
00599     }
00600 
00601     // build vertex block list
00602     std::vector< size_t > indices( array_length, 0 );
00603     size_t idx = 0;
00604     while( idx < vertices.size() )
00605     {
00606         DomainBlock block;
00607         block.firstHandle = vertices[idx];
00608         block.lastHandle  = vertices[idx];
00609 
00610         // find domain
00611         unsigned dom = 0;
00612         for( dom = 0; dom < array_length; ++dom )
00613         {
00614             const size_t i = indices[dom];
00615             if( i < sets[dom].vertices.size() && sets[dom].vertices[i] == vertices[idx] )
00616             {
00617                 ++indices[dom];
00618                 break;
00619             }
00620         }
00621         ++idx;
00622         if( dom == indices.size() )  // no domain
00623             continue;
00624 
00625         block.domain = sets[dom].domain;
00626         while( indices[dom] < sets[dom].vertices.size() && sets[dom].vertices[indices[dom]] == vertices[idx] )
00627         {
00628             block.lastHandle = vertices[idx];
00629             ++idx;
00630             ++indices[dom];
00631         }
00632         result.vertexList.push_back( block );
00633     }
00634 
00635     // build vertex block list
00636     indices.clear();
00637     indices.resize( array_length, 0 );
00638     idx = 0;
00639     while( idx < elements.size() )
00640     {
00641         DomainBlock block;
00642         block.firstHandle = elements[idx];
00643         block.lastHandle  = elements[idx];
00644 
00645         // find domain
00646         unsigned dom = 0;
00647         for( dom = 0; dom < array_length; ++dom )
00648         {
00649             const size_t i = indices[dom];
00650             if( i < sets[dom].elements.size() && sets[dom].elements[i] == elements[idx] )
00651             {
00652                 ++indices[dom];
00653                 break;
00654             }
00655         }
00656         ++idx;
00657         if( dom == indices.size() )  // no domain
00658             continue;
00659 
00660         block.domain = sets[dom].domain;
00661         while( indices[dom] < sets[dom].elements.size() && sets[dom].elements[indices[dom]] == elements[idx] )
00662         {
00663             block.lastHandle = elements[idx];
00664             ++idx;
00665             ++indices[dom];
00666         }
00667         result.elementList.push_back( block );
00668     }
00669     /*
00670       printf("\n");
00671       for (unsigned i = 0; i < result.vertexList.size(); ++i)
00672         printf("v %2lu-%2lu @ %p\n", (unsigned long)result.vertexList[i].firstHandle,
00673                                    (unsigned long)result.vertexList[i].lastHandle,
00674                                    result.vertexList[i].domain );
00675       for (unsigned i = 0; i < result.elementList.size(); ++i)
00676         printf("e %2lu-%2lu @ %p\n", (unsigned long)result.elementList[i].firstHandle,
00677                                    (unsigned long)result.elementList[i].lastHandle,
00678                                    result.elementList[i].domain );
00679       printf("\n");
00680     */
00681 }
00682 
00683 static bool next_vertex( Mesh* mesh, Mesh::VertexHandle& vtx, std::set< Mesh::VertexHandle >& unseen, MsqError& err )
00684 {
00685     std::vector< Mesh::ElementHandle > vtx_elems;
00686     std::vector< size_t > junk;
00687 
00688     mesh->vertices_get_attached_elements( &vtx, 1, vtx_elems, junk, err );
00689     MSQ_ERRZERO( err );
00690 
00691     std::vector< EntityTopology > elem_types( vtx_elems.size() );
00692     if( !vtx_elems.empty() )
00693     {
00694         mesh->elements_get_topologies( arrptr( vtx_elems ), arrptr( elem_types ), vtx_elems.size(), err );
00695         MSQ_ERRZERO( err );
00696     }
00697 
00698     std::vector< Mesh::VertexHandle > corners;
00699     for( size_t j = 0; j < vtx_elems.size(); ++j )
00700     {
00701         corners.clear();
00702         mesh->elements_get_attached_vertices( &vtx_elems[j], 1, corners, junk, err );
00703         MSQ_ERRZERO( err );
00704 
00705         unsigned nedges = TopologyInfo::edges( elem_types[j] );
00706         unsigned vidx   = std::find( corners.begin(), corners.end(), vtx ) - corners.begin();
00707 
00708         // Check mid-edge nodes first, if present
00709         int ho = TopologyInfo::higher_order( elem_types[j], corners.size(), err );
00710         MSQ_ERRZERO( err );
00711         // If element has mid-edge nodes *and* current vertex is a corner vertex
00712         if( ( ho & 2 ) && ( vidx < TopologyInfo::corners( elem_types[j] ) ) )
00713         {
00714             for( unsigned e = 0; e < nedges; ++e )
00715             {
00716                 const unsigned* edge_verts = TopologyInfo::edge_vertices( elem_types[j], e, err );
00717                 MSQ_ERRZERO( err );
00718                 if( edge_verts[0] == vidx || edge_verts[1] == vidx )
00719                 {
00720                     int idx = TopologyInfo::higher_order_from_side( elem_types[j], corners.size(), 1, e, err );
00721                     MSQ_ERRZERO( err );
00722                     std::set< Mesh::VertexHandle >::iterator f = unseen.find( corners[idx] );
00723                     if( f != unseen.end() )
00724                     {
00725                         vtx = *f;
00726                         unseen.erase( f );
00727                         return true;
00728                     }
00729                 }
00730             }
00731         }
00732 
00733         // if current vertx is a mid-edge node
00734         else if( ho & 2 )
00735         {
00736             unsigned d, e;
00737             TopologyInfo::side_from_higher_order( elem_types[j], corners.size(), vidx, d, e, err );
00738             MSQ_ERRZERO( err );
00739             if( d != 1 ) continue;
00740             const unsigned* edge_verts = TopologyInfo::edge_vertices( elem_types[j], e, err );
00741             MSQ_ERRZERO( err );
00742             for( int v = 0; v < 2; ++v )
00743             {
00744                 std::set< Mesh::VertexHandle >::iterator f = unseen.find( corners[edge_verts[v]] );
00745                 if( f != unseen.end() )
00746                 {
00747                     vtx = *f;
00748                     unseen.erase( f );
00749                     return true;
00750                 }
00751             }
00752         }
00753 
00754         else
00755         {
00756             for( unsigned e = 0; e < nedges; ++e )
00757             {
00758                 const unsigned* edge_verts = TopologyInfo::edge_vertices( elem_types[j], e, err );
00759                 MSQ_ERRZERO( err );
00760                 int idx;
00761                 if( edge_verts[0] == vidx )
00762                     idx = edge_verts[1];
00763                 else if( edge_verts[1] == vidx )
00764                     idx = edge_verts[0];
00765                 else
00766                     continue;
00767 
00768                 std::set< Mesh::VertexHandle >::iterator f = unseen.find( corners[idx] );
00769                 if( f != unseen.end() )
00770                 {
00771                     vtx = *f;
00772                     unseen.erase( f );
00773                     return true;
00774                 }
00775             }
00776         }
00777     }
00778 
00779     return false;
00780 }
00781 
00782 void DomainClassifier::test_valid_classification( Mesh* mesh, MsqError& err )
00783 {
00784     size_t i;
00785 
00786     // Get all mesh entities
00787     std::vector< Mesh::VertexHandle > vertices;
00788     std::vector< Mesh::ElementHandle > elements;
00789     mesh->get_all_vertices( vertices, err );MSQ_ERRRTN( err );
00790     mesh->get_all_elements( elements, err );MSQ_ERRRTN( err );
00791     std::sort( vertices.begin(), vertices.end() );
00792     std::sort( elements.begin(), elements.end() );
00793 
00794     // Get contents of each domain.
00795     std::map< MeshDomain*, int >::iterator iter;
00796     std::map< MeshDomain*, int > idxmap;
00797     for( i = 0; i < vertexList.size(); ++i )
00798         idxmap[vertexList[i].domain] = 0;
00799     for( i = 0; i < elementList.size(); ++i )
00800         idxmap[elementList[i].domain] = 0;
00801     std::vector< DomainSet > domains( idxmap.size() );
00802     int idx = 0;
00803     for( iter = idxmap.begin(); iter != idxmap.end(); ++iter )
00804     {
00805         iter->second        = idx;
00806         domains[idx].domain = iter->first;
00807         ++idx;
00808     }
00809     for( i = 0; i < vertexList.size(); ++i )
00810     {
00811         std::vector< Mesh::VertexHandle >::const_iterator s, e;
00812         s              = std::lower_bound( vertices.begin(), vertices.end(), vertexList[i].firstHandle );
00813         e              = std::upper_bound( vertices.begin(), vertices.end(), vertexList[i].lastHandle );
00814         DomainSet& set = domains[idxmap[vertexList[i].domain]];
00815         std::copy( s, e, std::back_inserter( set.vertices ) );
00816     }
00817     for( i = 0; i < elementList.size(); ++i )
00818     {
00819         std::vector< Mesh::ElementHandle >::const_iterator s, e;
00820         s              = std::lower_bound( elements.begin(), elements.end(), elementList[i].firstHandle );
00821         e              = std::upper_bound( elements.begin(), elements.end(), elementList[i].lastHandle );
00822         DomainSet& set = domains[idxmap[elementList[i].domain]];
00823         std::copy( s, e, std::back_inserter( set.elements ) );
00824     }
00825 
00826     // guess geometric dimension for each domain
00827     std::map< DomainSet*, int > dimmap;
00828     std::vector< unsigned short > dof;
00829     for( i = 0; i < domains.size(); ++i )
00830     {
00831         if( !domains[i].elements.empty() )
00832             dimmap[&domains[i]] = 2;
00833         else
00834         {
00835             dof.resize( domains[i].vertices.size() );
00836             domains[i].domain->domain_DoF( &( domains[i].vertices[0] ), arrptr( dof ), dof.size(), err );MSQ_ERRRTN( err );
00837             unsigned short dim  = *std::max_element( dof.begin(), dof.end() );
00838             dimmap[&domains[i]] = dim;
00839         }
00840     }
00841 
00842     // group domains by dimension
00843     std::vector< DomainSet* > points, curves, surfaces;
00844     for( std::map< DomainSet*, int >::iterator it = dimmap.begin(); it != dimmap.end(); ++it )
00845     {
00846         switch( it->second )
00847         {
00848             case 0:
00849                 points.push_back( it->first );
00850                 break;
00851             case 1:
00852                 curves.push_back( it->first );
00853                 break;
00854             case 2:
00855                 surfaces.push_back( it->first );
00856                 break;
00857             default:
00858                 MSQ_SETERR( err )
00859                 ( MsqError::INVALID_STATE, "Invalid domain dimension: %d\n", it->second );
00860                 break;
00861         }
00862     }
00863 
00864     // check that each point-domain has a single vertex
00865     for( i = 0; i < points.size(); ++i )
00866     {
00867         if( points[i]->vertices.size() != 1 || !points[i]->elements.empty() )
00868         {
00869             MSQ_SETERR( err )( "Point domain mesh not a single vertex.", MsqError::INVALID_STATE );
00870             return;
00871         }
00872     }
00873 
00874     // check that each curve domain has a chain of connected edges
00875     std::set< Mesh::VertexHandle > unseen;
00876     for( i = 0; i < curves.size(); ++i )
00877     {
00878         if( !curves[i]->elements.empty() )
00879         {
00880             MSQ_SETERR( err )( "Elements associated with 1D domain.", MsqError::INVALID_STATE );
00881             return;
00882         }
00883 
00884         unseen.clear();
00885         std::copy( curves[i]->vertices.begin(), curves[i]->vertices.end(), std::inserter( unseen, unseen.begin() ) );
00886 
00887         const Mesh::VertexHandle first_vtx = *unseen.begin();
00888         unseen.erase( unseen.begin() );
00889         Mesh::VertexHandle vtx = first_vtx;
00890         // find chain of vertices
00891         while( next_vertex( mesh, vtx, unseen, err ) )
00892             MSQ_ERRRTN( err );
00893         // search again from the starting vertex because
00894         // it probably wasn't the first one on the curve
00895         vtx = first_vtx;
00896         while( next_vertex( mesh, vtx, unseen, err ) )
00897             MSQ_ERRRTN( err );
00898         // were all vertices in a chain?
00899         if( !unseen.empty() )
00900         {
00901             MSQ_SETERR( err )
00902             ( "Curve domain contains vertices not in a simply connected chain.", MsqError::INVALID_STATE );
00903             return;
00904         }
00905     }
00906 
00907     std::set< Mesh::VertexHandle > seen;
00908     std::set< Mesh::ElementHandle > remaining;
00909     std::vector< Mesh::VertexHandle > verts, verts2;
00910     std::vector< Mesh::ElementHandle > stack, vert_elems;
00911     std::vector< EntityTopology > types;
00912     std::vector< size_t > junk;
00913     // if surface contains elements...
00914     for( i = 0; i < surfaces.size(); ++i )
00915     {
00916         if( surfaces[i]->elements.empty() ) continue;
00917 
00918         // Check that any vertices on surface are contained in an
00919         // element on the surface.
00920         verts.clear();
00921         mesh->elements_get_attached_vertices( &( surfaces[i]->elements[0] ), surfaces[i]->elements.size(), verts, junk,
00922                                               err );MSQ_ERRRTN( err );
00923         seen.clear();
00924         std::copy( verts.begin(), verts.end(), std::inserter( seen, seen.begin() ) );
00925 
00926         std::vector< Mesh::VertexHandle >::const_iterator v;
00927         for( v = surfaces[i]->vertices.begin(); v != surfaces[i]->vertices.end(); ++v )
00928         {
00929             std::set< Mesh::VertexHandle >::iterator j = seen.find( *v );
00930             if( j == seen.end() )
00931             {
00932                 MSQ_SETERR( err )
00933                 ( "Vertex on surface domain not in any element.", MsqError::INVALID_STATE );
00934                 return;
00935             }
00936         }
00937 
00938         // check that elements form 2D patch
00939         stack.clear();
00940         remaining.clear();
00941         std::copy( surfaces[i]->elements.begin(), surfaces[i]->elements.end(),
00942                    std::inserter( remaining, remaining.begin() ) );
00943         stack.push_back( *remaining.begin() );
00944         remaining.erase( remaining.begin() );
00945         while( !stack.empty() )
00946         {
00947             Mesh::ElementHandle elem = stack.back();
00948             stack.pop_back();
00949             verts.clear();
00950             mesh->elements_get_attached_vertices( &elem, 1, verts, junk, err );MSQ_ERRRTN( err );
00951             // for each edge
00952             for( size_t j = 0; j < verts.size(); ++j )
00953             {
00954                 Mesh::VertexHandle v1 = verts[j], v2 = verts[( j + 1 ) % verts.size()];
00955                 vert_elems.clear();
00956                 mesh->vertices_get_attached_elements( &v1, 1, vert_elems, junk, err );MSQ_ERRRTN( err );
00957                 types.resize( vert_elems.size() );
00958                 if( !vert_elems.empty() )
00959                 {
00960                     mesh->elements_get_topologies( arrptr( vert_elems ), arrptr( types ), vert_elems.size(), err );MSQ_ERRRTN( err );
00961                 }
00962                 while( !vert_elems.empty() )
00963                 {
00964                     Mesh::ElementHandle e2 = vert_elems.back();
00965                     EntityTopology type    = types.back();
00966                     vert_elems.pop_back();
00967                     types.pop_back();
00968                     if( TopologyInfo::dimension( type ) != 2 ) continue;
00969                     verts2.clear();
00970                     mesh->elements_get_attached_vertices( &e2, 1, verts2, junk, err );MSQ_ERRRTN( err );
00971                     size_t jdx = std::find( verts2.begin(), verts2.end(), v1 ) - verts2.begin();
00972                     if( verts2[( jdx + 1 ) % verts2.size()] != v2 &&
00973                         verts2[( jdx + verts2.size() - 1 ) % verts2.size()] != v2 )
00974                         continue;
00975                     std::set< Mesh::ElementHandle >::iterator r = remaining.find( e2 );
00976                     if( r == remaining.end() ) continue;
00977 
00978                     stack.push_back( *r );
00979                     remaining.erase( r );
00980                 }
00981             }
00982         }
00983 
00984         if( !remaining.empty() )
00985         {
00986             MSQ_SETERR( err )
00987             ( "Surface mesh not a single, simply-connected patch", MsqError::INVALID_STATE );
00988             return;
00989         }
00990     }
00991 
00992     // check that sides of volume elements that are on surface
00993     // form simply connected patch
00994     for( i = 0; i < surfaces.size(); ++i )
00995     {
00996         // build list of sides known to be on surface
00997         std::sort( surfaces[i]->vertices.begin(), surfaces[i]->vertices.end() );
00998         std::set< std::pair< Mesh::ElementHandle, int > > sides;
00999         for( size_t j = 0; j < surfaces[i]->vertices.size(); ++j )
01000         {
01001             Mesh::VertexHandle v = surfaces[i]->vertices[j];
01002             vert_elems.clear();
01003             mesh->vertices_get_attached_elements( &v, 1, vert_elems, junk, err );MSQ_ERRRTN( err );
01004             types.resize( vert_elems.size() );
01005             if( !vert_elems.empty() )
01006             {
01007                 mesh->elements_get_topologies( arrptr( vert_elems ), arrptr( types ), vert_elems.size(), err );MSQ_ERRRTN( err );
01008             }
01009             while( !vert_elems.empty() )
01010             {
01011                 Mesh::ElementHandle e = vert_elems.back();
01012                 EntityTopology type   = types.back();
01013                 vert_elems.pop_back();
01014                 types.pop_back();
01015                 if( TopologyInfo::dimension( type ) != 3 ) continue;
01016                 verts.clear();
01017                 mesh->elements_get_attached_vertices( &e, 1, verts, junk, err );MSQ_ERRRTN( err );
01018 
01019                 for( unsigned s = 0; s < TopologyInfo::faces( type ); ++s )
01020                 {
01021                     unsigned n;
01022                     const unsigned* si = TopologyInfo::face_vertices( type, s, n );
01023                     unsigned ns        = 0;
01024                     for( unsigned k = 0; k < n; ++k )
01025                     {
01026                         if( std::binary_search( surfaces[i]->vertices.begin(), surfaces[i]->vertices.end(),
01027                                                 verts[si[k]] ) )
01028                             ++ns;
01029                     }
01030                     if( ns >= 3 ) sides.insert( std::pair< Mesh::ElementHandle, int >( e, s ) );
01031                 }
01032             }
01033         }
01034 
01035         std::vector< std::pair< Mesh::ElementHandle, int > > sstack;
01036         sstack.push_back( *sides.begin() );
01037         sides.erase( sides.begin() );
01038         while( !sstack.empty() )
01039         {
01040             Mesh::ElementHandle e = sstack.back().first;
01041             int s                 = sstack.back().second;
01042             sstack.pop_back();
01043 
01044             verts.clear();
01045             mesh->elements_get_attached_vertices( &e, 1, verts, junk, err );MSQ_ERRRTN( err );
01046             EntityTopology type;
01047             mesh->elements_get_topologies( &e, &type, 1, err );MSQ_ERRRTN( err );
01048             unsigned n;
01049             const unsigned* si = TopologyInfo::face_vertices( type, s, n );
01050 
01051             // for each edge
01052             for( unsigned j = 0; j < n; ++j )
01053             {
01054                 Mesh::VertexHandle v1 = verts[si[j]], v2 = verts[si[( j + 1 ) % n]];
01055                 vert_elems.clear();
01056                 mesh->vertices_get_attached_elements( &v1, 1, vert_elems, junk, err );MSQ_ERRRTN( err );
01057                 types.resize( vert_elems.size() );
01058                 if( !vert_elems.empty() )
01059                 {
01060                     mesh->elements_get_topologies( arrptr( vert_elems ), arrptr( types ), vert_elems.size(), err );MSQ_ERRRTN( err );
01061                 }
01062                 while( !vert_elems.empty() )
01063                 {
01064                     Mesh::ElementHandle e2 = vert_elems.back();
01065                     EntityTopology type2   = types.back();
01066                     vert_elems.pop_back();
01067                     types.pop_back();
01068                     if( TopologyInfo::dimension( type ) != 3 ) continue;
01069                     verts2.clear();
01070                     mesh->elements_get_attached_vertices( &e2, 1, verts2, junk, err );MSQ_ERRRTN( err );
01071                     // for each face
01072                     for( unsigned s2 = 0; s2 < TopologyInfo::faces( type2 ); ++s2 )
01073                     {
01074                         std::set< std::pair< Mesh::ElementHandle, int > >::iterator side;
01075                         side = sides.find( std::pair< Mesh::ElementHandle, int >( e2, s2 ) );
01076                         if( side == sides.end() ) continue;
01077 
01078                         unsigned n2;
01079                         const unsigned* si2 = TopologyInfo::face_vertices( type2, s2, n2 );
01080                         unsigned jdx;
01081                         for( jdx = 0; jdx < n2; ++jdx )
01082                             if( verts2[si2[jdx]] == v1 ) break;
01083                         assert( jdx < n2 );
01084 
01085                         if( verts2[si2[( jdx + 1 ) % n2]] == v2 || verts2[si2[( jdx + n2 - 1 ) % n2]] == v2 )
01086                         {
01087                             sstack.push_back( *side );
01088                             sides.erase( side );
01089                         }
01090                     }
01091                 }
01092             }
01093         }
01094 
01095         if( !sides.empty() )
01096         {
01097             MSQ_SETERR( err )
01098             ( "Surface mesh not a single, simply-connected patch", MsqError::INVALID_STATE );
01099             return;
01100         }
01101     }
01102 }
01103 
01104 static inline bool operator<( const Mesh::EntityHandle h, const DomainClassifier::DomainBlock& b )
01105 {
01106     return h < b.firstHandle;
01107 }
01108 
01109 static inline bool operator<( const DomainClassifier::DomainBlock& b, const Mesh::EntityHandle h )
01110 {
01111     return b.lastHandle < h;
01112 }
01113 
01114 static inline bool operator<( const DomainClassifier::DomainBlock& b, const DomainClassifier::DomainBlock& c )
01115 {
01116     return b.lastHandle < c.firstHandle;
01117 }
01118 
01119 MeshDomain* DomainClassifier::find_domain( Mesh::EntityHandle handle,
01120                                            const std::vector< DomainClassifier::DomainBlock >& list )
01121 {
01122     std::vector< DomainClassifier::DomainBlock >::const_iterator i;
01123     i = std::lower_bound( list.begin(), list.end(), handle );
01124     return ( i != list.end() && i->firstHandle <= handle ) ? i->domain : NULL;
01125 }
01126 
01127 void DomainClassifier::snap_to( Mesh::EntityHandle entity_handle, Vector3D& coordinate ) const
01128 {
01129     if( const MeshDomain* dom = find_vertex_domain( entity_handle ) ) dom->snap_to( entity_handle, coordinate );
01130 }
01131 
01132 void DomainClassifier::vertex_normal_at( Mesh::VertexHandle entity_handle, Vector3D& coordinate ) const
01133 {
01134     if( const MeshDomain* dom = find_vertex_domain( entity_handle ) )
01135         dom->vertex_normal_at( entity_handle, coordinate );
01136 }
01137 
01138 void DomainClassifier::element_normal_at( Mesh::ElementHandle entity_handle, Vector3D& coordinate ) const
01139 {
01140     if( const MeshDomain* dom = find_element_domain( entity_handle ) )
01141         dom->element_normal_at( entity_handle, coordinate );
01142 }
01143 
01144 void DomainClassifier::vertex_normal_at( const Mesh::VertexHandle* handles, Vector3D coordinates[], unsigned count,
01145                                          MsqError& err ) const
01146 {
01147     for( unsigned i = 0; i < count; ++i )
01148     {
01149         const MeshDomain* dom = find_vertex_domain( handles[i] );
01150         if( !dom )
01151         {
01152             MSQ_SETERR( err )( MsqError::INVALID_ARG );
01153             return;
01154         }
01155         dom->vertex_normal_at( handles + i, coordinates + i, 1, err );MSQ_ERRRTN( err );
01156     }
01157 }
01158 
01159 void DomainClassifier::closest_point( Mesh::VertexHandle handle, const Vector3D& position, Vector3D& closest,
01160                                       Vector3D& normal, MsqError& err ) const
01161 {
01162     if( const MeshDomain* dom = find_vertex_domain( handle ) )
01163         dom->closest_point( handle, position, closest, normal, err );
01164 }
01165 
01166 void DomainClassifier::domain_DoF( const Mesh::VertexHandle* handles, unsigned short* dof_array, size_t num_handles,
01167                                    MsqError& err ) const
01168 {
01169     for( size_t i = 0; i < num_handles; ++i )
01170     {
01171         const MeshDomain* dom = find_vertex_domain( handles[i] );
01172         if( !dom )
01173             dof_array[i] = 3;
01174         else
01175         {
01176             dom->domain_DoF( handles + i, dof_array + i, 1, err );MSQ_ERRRTN( err );
01177         }
01178     }
01179 }
01180 
01181 DomainClassifier::~DomainClassifier()
01182 {
01183     if( deleteSubDomains ) delete_all_sub_domains();
01184 }
01185 
01186 void DomainClassifier::delete_all_sub_domains()
01187 {
01188     // get unique list of domains
01189     std::set< MeshDomain* > domains;
01190     std::vector< DomainBlock >::iterator i;
01191     for( i = vertexList.begin(); i != vertexList.end(); ++i )
01192         domains.insert( i->domain );
01193     for( i = elementList.begin(); i != elementList.end(); ++i )
01194         domains.insert( i->domain );
01195     std::set< MeshDomain* >::iterator j;
01196     for( j = domains.begin(); j != domains.end(); ++j )
01197         delete *j;
01198     clear();
01199 }
01200 
01201 }  // namespace MBMesquite
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines