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