Assign subsets of a mesh do different domains. More...

#include <DomainClassifier.hpp>

Inheritance diagram for MBMesquite::DomainClassifier:

Collaboration diagram for MBMesquite::DomainClassifier:

Classes
struct	DomainBlock
struct	DomainSet
Public Member Functions
MESQUITE_EXPORT void	test_valid_classification (Mesh *mesh, MsqError &err)
	Check that classification maps to B-Rep topology.
MESQUITE_EXPORT	DomainClassifier ()
virtual MESQUITE_EXPORT	~DomainClassifier ()
virtual MESQUITE_EXPORT void	snap_to (Mesh::VertexHandle entity_handle, Vector3D &coordinate) const
virtual MESQUITE_EXPORT void	vertex_normal_at (Mesh::VertexHandle entity_handle, Vector3D &coordinate) const
virtual MESQUITE_EXPORT void	element_normal_at (Mesh::ElementHandle entity_handle, Vector3D &coordinate) const
virtual MESQUITE_EXPORT void	vertex_normal_at (const Mesh::VertexHandle *handles, Vector3D coordinates[], unsigned count, MsqError &err) const
	evaluate surface normals
virtual MESQUITE_EXPORT void	closest_point (Mesh::VertexHandle handle, const Vector3D &position, Vector3D &closest, Vector3D &normal, MsqError &err) const
	evaluate closest point and normal
virtual MESQUITE_EXPORT void	domain_DoF (const Mesh::VertexHandle handle_array, unsigned short dof_array, size_t num_handles, MsqError &err) const
MESQUITE_EXPORT void	clear ()
	Clear all data, including MeshDomain list.
MESQUITE_EXPORT MeshDomain *	find_vertex_domain (Mesh::VertexHandle vertex) const
MESQUITE_EXPORT MeshDomain *	find_element_domain (Mesh::ElementHandle element) const
MESQUITE_EXPORT void	delete_sub_domains (bool yesno)
MESQUITE_EXPORT void	delete_all_sub_domains ()
Static Public Member Functions
static MESQUITE_EXPORT void	classify_by_tag (DomainClassifier &result, Mesh mesh, const char tag_name, MeshDomain *domain_array, const int id_array, unsigned array_length, MsqError &err)
	Classify mesh entities using tag values.
static MESQUITE_EXPORT void	classify_skin_geometrically (DomainClassifier &result, Mesh mesh, double tolerance, MeshDomain domain_array, const int dimension_array, unsigned array_length, MsqError &err)
	Skin mesh and classify skin entities geometrically.
static MESQUITE_EXPORT void	classify_geometrically (DomainClassifier &result, Mesh mesh, double tolerance, MeshDomain domain_array, const int dimension_array, unsigned array_length, MsqError &err)
	Classify all mesh entities geometrically.
static MESQUITE_EXPORT void	classify_by_handle (DomainClassifier &result, Mesh mesh, DomainSet domain_set_array, unsigned array_length, MsqError &err)
	Specify classification explicitly for each entity.
Static Private Member Functions
static MeshDomain *	find_domain (Mesh::EntityHandle handle, const std::vector< DomainBlock > &list)
Private Attributes
bool	deleteSubDomains
std::vector< DomainBlock >	vertexList
std::vector< DomainBlock >	elementList

Detailed Description

Assign subsets of a mesh do different domains.

Provide classification of mesh entities to domain(s). For example, given a domain defined by multiple surfaces and curves of intersection, associate the appropriate elements and vertices to the appropriate geometric entities.

Definition at line 48 of file DomainClassifier.hpp.

Constructor & Destructor Documentation

MESQUITE_EXPORT MBMesquite::DomainClassifier::DomainClassifier ( ) [inline]

Definition at line 179 of file DomainClassifier.hpp.

: deleteSubDomains( false ) {}

MBMesquite::DomainClassifier::~DomainClassifier ( ) [virtual]

Definition at line 1248 of file DomainClassifier.cpp.

References delete_all_sub_domains(), and deleteSubDomains.

{
    if( deleteSubDomains ) delete_all_sub_domains();
}

Member Function Documentation

void MBMesquite::DomainClassifier::classify_by_handle	(	DomainClassifier &	result,
		Mesh *	mesh,
		DomainSet *	domain_set_array,
		unsigned	array_length,
		MsqError &	err
	)		`[static]`

Specify classification explicitly for each entity.

For each mesh element and vertex that is classified to a domain, specify that classification explicitly.

Parameters:

result	DomainClassifier to popupate
mesh	The MBMesquite::Mesh instance
domain_set_array	Array of DomainClassifier::DomainSet structs specifying for each MeshDomain: the domain instance and the elements and vertices associated with the domain.
array_length	Length of 'domain_set_array'

Definition at line 638 of file DomainClassifier.cpp.

References clear(), MBMesquite::DomainClassifier::DomainSet::domain, MBMesquite::DomainClassifier::DomainBlock::domain, elementList, MBMesquite::DomainClassifier::DomainSet::elements, MBMesquite::DomainClassifier::DomainBlock::firstHandle, MBMesquite::Mesh::get_all_elements(), MBMesquite::Mesh::get_all_vertices(), MBMesquite::DomainClassifier::DomainBlock::lastHandle, MSQ_ERRRTN, vertexList, and MBMesquite::DomainClassifier::DomainSet::vertices.

Referenced by classify_by_tag(), classify_geometrically(), classify_skin_geometrically(), get_planar_example(), DomainClassifierTest::test_classify_by_handle(), HigherOrderTest::test_tri_open_domain(), and DomainClassifierTest::test_valid_classification().

{
    result.clear();

    // get all vertices and elements
    std::vector< Mesh::VertexHandle > vertices;
    std::vector< Mesh::ElementHandle > elements;
    mesh->get_all_vertices( vertices, err );MSQ_ERRRTN( err );
    mesh->get_all_elements( elements, err );MSQ_ERRRTN( err );

    // sort all arrays so we can merge
    std::sort( vertices.begin(), vertices.end() );
    std::sort( elements.begin(), elements.end() );
    for( unsigned i = 0; i < array_length; ++i )
    {
        std::sort( sets[i].vertices.begin(), sets[i].vertices.end() );
        std::sort( sets[i].elements.begin(), sets[i].elements.end() );
    }

    // build vertex block list
    std::vector< size_t > indices( array_length, 0 );
    size_t idx = 0;
    while( idx < vertices.size() )
    {
        DomainBlock block;
        block.firstHandle = vertices[idx];
        block.lastHandle  = vertices[idx];

        // find domain
        unsigned dom = 0;
        for( dom = 0; dom < array_length; ++dom )
        {
            const size_t i = indices[dom];
            if( i < sets[dom].vertices.size() && sets[dom].vertices[i] == vertices[idx] )
            {
                ++indices[dom];
                break;
            }
        }
        ++idx;
        if( dom == indices.size() )  // no domain
            continue;

        block.domain = sets[dom].domain;
        while( indices[dom] < sets[dom].vertices.size() && sets[dom].vertices[indices[dom]] == vertices[idx] )
        {
            block.lastHandle = vertices[idx];
            ++idx;
            ++indices[dom];
        }
        result.vertexList.push_back( block );
    }

    // build vertex block list
    indices.clear();
    indices.resize( array_length, 0 );
    idx = 0;
    while( idx < elements.size() )
    {
        DomainBlock block;
        block.firstHandle = elements[idx];
        block.lastHandle  = elements[idx];

        // find domain
        unsigned dom = 0;
        for( dom = 0; dom < array_length; ++dom )
        {
            const size_t i = indices[dom];
            if( i < sets[dom].elements.size() && sets[dom].elements[i] == elements[idx] )
            {
                ++indices[dom];
                break;
            }
        }
        ++idx;
        if( dom == indices.size() )  // no domain
            continue;

        block.domain = sets[dom].domain;
        while( indices[dom] < sets[dom].elements.size() && sets[dom].elements[indices[dom]] == elements[idx] )
        {
            block.lastHandle = elements[idx];
            ++idx;
            ++indices[dom];
        }
        result.elementList.push_back( block );
    }
    /*
      printf("\n");
      for (unsigned i = 0; i < result.vertexList.size(); ++i)
        printf("v %2lu-%2lu @ %p\n", (unsigned long)result.vertexList[i].firstHandle,
                                   (unsigned long)result.vertexList[i].lastHandle,
                                   result.vertexList[i].domain );
      for (unsigned i = 0; i < result.elementList.size(); ++i)
        printf("e %2lu-%2lu @ %p\n", (unsigned long)result.elementList[i].firstHandle,
                                   (unsigned long)result.elementList[i].lastHandle,
                                   result.elementList[i].domain );
      printf("\n");
    */
}

void MBMesquite::DomainClassifier::classify_by_tag	(	DomainClassifier &	result,
		Mesh *	mesh,
		const char *	tag_name,
		MeshDomain **	domain_array,
		const int *	id_array,
		unsigned	array_length,
		MsqError &	err
	)		`[static]`

Classify mesh entities using tag values.

Given a list of MeshDomain instances, a unique integer ID for each MeshDomain, and a tag name: classify mesh entities by matching the tag value on each mesh entity to the list of MeshDomain IDs.

Parameters:

result	DomainClassifier to popupate
mesh	The MBMesquite::Mesh instance
tag_name	Tag containing integer domain ID for each mesh entity.
domain_array	Array of MeshDomain instances.
id_array	Array of integer MeshDomain IDs.
array_length	Length of 'domain_array' and 'id_array'

Definition at line 572 of file DomainClassifier.cpp.

References MBMesquite::arrptr(), classify_by_handle(), MBMesquite::MsqError::clear(), MBMesquite::Mesh::get_all_elements(), MBMesquite::Mesh::get_all_vertices(), MBMesquite::Mesh::INT, MSQ_ERRRTN, MSQ_SETERR, size, MBMesquite::Mesh::tag_get(), MBMesquite::Mesh::tag_get_element_data(), MBMesquite::Mesh::tag_get_vertex_data(), MBMesquite::MsqError::TAG_NOT_FOUND, and MBMesquite::Mesh::tag_properties().

Referenced by DomainClassifierTest::test_classify_by_tag().

{
    TagHandle tag = mesh->tag_get( tag_name, err );MSQ_ERRRTN( err );

    std::string name2;
    Mesh::TagType type;
    unsigned size;
    mesh->tag_properties( tag, name2, type, size, err );MSQ_ERRRTN( err );
    if( type != Mesh::INT || size != 1 )
    {
        MSQ_SETERR( err )
        ( MsqError::TAG_NOT_FOUND, "Tag does not contain single integer value: %s\n", tag_name );
        return;
    }

    std::vector< DomainSet > sets( array_length );
    std::map< int, DomainSet* > idmap;
    for( unsigned i = 0; i < array_length; ++i )
    {
        sets[i].domain     = domain_array[i];
        idmap[id_array[i]] = &sets[i];
    }

    std::vector< Mesh::VertexHandle > vertices;
    std::vector< Mesh::ElementHandle > elements;
    mesh->get_all_vertices( vertices, err );MSQ_ERRRTN( err );
    mesh->get_all_elements( elements, err );MSQ_ERRRTN( err );
    if( vertices.empty() ) return;

    std::map< int, DomainSet* >::const_iterator iter;
    std::vector< int > ids( vertices.size() );
    mesh->tag_get_vertex_data( tag, vertices.size(), arrptr( vertices ), arrptr( ids ), err );MSQ_ERRRTN( err );
    for( size_t j = 0; j < vertices.size(); ++j )
    {
        iter = idmap.find( ids[j] );
        if( iter != idmap.end() ) iter->second->vertices.push_back( vertices[j] );
    }

    ids.clear();
    ids.resize( elements.size() );
    if( !elements.empty() )
    {
        mesh->tag_get_element_data( tag, elements.size(), arrptr( elements ), arrptr( ids ), err );
        if( err )
        {
            err.clear();
        }
        else
        {
            for( size_t k = 0; k < elements.size(); ++k )
            {
                iter = idmap.find( ids[k] );
                if( iter != idmap.end() ) iter->second->elements.push_back( elements[k] );
            }
        }
    }

    if( !sets.empty() ) classify_by_handle( result, mesh, arrptr( sets ), sets.size(), err );
}

void MBMesquite::DomainClassifier::classify_geometrically	(	DomainClassifier &	result,
		Mesh *	mesh,
		double	tolerance,
		MeshDomain **	domain_array,
		const int *	dimension_array,
		unsigned	array_length,
		MsqError &	err
	)		`[static]`

Classify all mesh entities geometrically.

Classify entities by distance from vertex coordiantes to mesh domain instances.

NOTE: Mesquite's limited MeshDomain callback interface is not sufficient for robust geometric classification. If the mesh is not sufficiently refined, this method may produce invalid results.

NOTE: This method is likely to fail for many domains for volume meshes. Use classify_skin_geometrically for volume meshes.

Parameters:

result	DomainClassifier to popupate
mesh	The MBMesquite::Mesh instance
tolerance	Maximum distance a vertex may deviate from its domain.
domain_array	Array of MeshDomain instances.
dimension_array	Topological dimensiono of each MeshDomain
array_length	Length of 'domain_array' and 'dimension_array'

Definition at line 483 of file DomainClassifier.cpp.

References MBMesquite::arrptr(), classify_by_handle(), MBMesquite::dimension_sort_domains(), domains, MBMesquite::Mesh::elements_get_topologies(), MBMesquite::geom_classify_elements(), MBMesquite::geom_classify_vertices(), MBMesquite::Mesh::get_all_elements(), MBMesquite::Mesh::get_all_vertices(), MSQ_ERRRTN, tolerance, and MBMesquite::vert_classify_elements().

Referenced by PatchDataTest::get_quad8_mesh_and_domain(), process_domain_args(), and DomainClassifierTest::test_classify_by_geometry().

{
    if( !array_length ) return;

    std::vector< DomainSet > domains;
    int dim_indices[4];
    dimension_sort_domains( domain_array, dimension_array, array_length, domains, dim_indices, err );MSQ_ERRRTN( err );

    // classify vertices by position
    std::vector< Mesh::VertexHandle > vertices;
    mesh->get_all_vertices( vertices, err );MSQ_ERRRTN( err );
    geom_classify_vertices( mesh, arrptr( domains ), dim_indices[3], vertices, tolerance, err );MSQ_ERRRTN( err );

    // get elements and types
    std::vector< Mesh::ElementHandle > elems;
    mesh->get_all_elements( elems, err );MSQ_ERRRTN( err );
    if( elems.empty() ) return;
    std::vector< EntityTopology > types( elems.size() );
    mesh->elements_get_topologies( arrptr( elems ), arrptr( types ), elems.size(), err );MSQ_ERRRTN( err );

    // get rid of elements w/ dimension other than 2
    size_t w = 0;
    for( size_t r = 0; r < elems.size(); ++r )
    {
        if( TopologyInfo::dimension( types[r] ) == 2 )
        {
            elems[w] = elems[r];
            ++w;
        }
    }
    elems.resize( w );

    // classify elements using vertex classification
    std::vector< Mesh::ElementHandle > unknown;
    vert_classify_elements( mesh, arrptr( domains ), domains.size(), dim_indices, elems, unknown, err );MSQ_ERRRTN( err );

    // classify unknown elements geometrically
    elems.swap( unknown );
    unknown.clear();
    geom_classify_elements( mesh, arrptr( domains ), domains.size(), dim_indices, elems, unknown, tolerance, err );

    classify_by_handle( result, mesh, arrptr( domains ), domains.size(), err );
}

void MBMesquite::DomainClassifier::classify_skin_geometrically	(	DomainClassifier &	result,
		Mesh *	mesh,
		double	tolerance,
		MeshDomain **	domain_array,
		const int *	dimension_array,
		unsigned	array_length,
		MsqError &	err
	)		`[static]`

Skin mesh and classify skin entities geometrically.

Calculate the boundary of the mesh, and classify boundary entities by calculating the distance from the mesh domain to each vertex.

NOTE: Mesquite's limited MeshDomain callback interface is not sufficient for robust geometric classification. If the mesh is not sufficiently refined, this method may produce invalid results.

NOTE: This method should not be used for surface meshes. Everything in a surface mesh should be classified to a domain. not just the skin. Use classify_geometrically instead.

Parameters:

result	DomainClassifier to popupate
mesh	The MBMesquite::Mesh instance
tolerance	Maximum distance a vertex may deviate from its domain.
domain_array	Array of MeshDomain instances.
dimension_array	Topological dimensiono of each MeshDomain
array_length	Length of 'domain_array' and 'dimension_array'

Definition at line 533 of file DomainClassifier.cpp.

References MBMesquite::arrptr(), classify_by_handle(), MBMesquite::dimension_sort_domains(), domains, MBMesquite::find_skin(), MBMesquite::geom_classify_elements(), MBMesquite::geom_classify_vertices(), MBMesquite::MsqError::INVALID_MESH, MSQ_ERRRTN, MSQ_SETERR, tolerance, and MBMesquite::vert_classify_elements().

Referenced by get_cut_cube_example(), get_hex_3d_part_example(), get_sphere_cube_example(), get_sphere_cylinder_example(), SlaveBoundaryVerticesTest::make_mesh(), process_domain_args(), and DomainClassifierTest::test_classify_skin().

{
    if( !array_length ) return;

    std::vector< DomainSet > domains;
    int dim_indices[4];
    dimension_sort_domains( domain_array, dimension_array, array_length, domains, dim_indices, err );MSQ_ERRRTN( err );

    std::vector< Mesh::VertexHandle > vertices;
    std::vector< Mesh::ElementHandle > elements;
    find_skin( mesh, vertices, elements, err );MSQ_ERRRTN( err );

    // classify vertices by position
    geom_classify_vertices( mesh, arrptr( domains ), dim_indices[3], vertices, tolerance, err );MSQ_ERRRTN( err );

    // classify elements using vertex classification
    std::vector< Mesh::ElementHandle > unknown;
    vert_classify_elements( mesh, arrptr( domains ), domains.size(), dim_indices, elements, unknown, err );MSQ_ERRRTN( err );

    // classify unknown elements geometrically
    elements.swap( unknown );
    unknown.clear();
    geom_classify_elements( mesh, arrptr( domains ), domains.size(), dim_indices, elements, unknown, tolerance, err );

    if( !unknown.empty() )
    {
        MSQ_SETERR( err )( "Count not classify all skin elements", MsqError::INVALID_MESH );
        return;
    }

    classify_by_handle( result, mesh, arrptr( domains ), domains.size(), err );
}

MESQUITE_EXPORT void MBMesquite::DomainClassifier::clear ( ) [inline]

Clear all data, including MeshDomain list.

Definition at line 212 of file DomainClassifier.hpp.

References elementList, and vertexList.

Referenced by classify_by_handle(), and delete_all_sub_domains().

    {
        vertexList.clear();
        elementList.clear();
    }

void MBMesquite::DomainClassifier::closest_point	(	Mesh::VertexHandle	handle,
		const Vector3D &	position,
		Vector3D &	closest,
		Vector3D &	normal,
		MsqError &	err
	)		const `[virtual]`

evaluate closest point and normal

Given a position in space, return the closest position in the domain and the domain normal at that point.

Parameters:

entity_handle	Evaluate the subset of the domain contianing this entity
position	Input position for which to evaluate
closest	Closest position in the domain.
normal	Domain normal at the location of 'closest'

Implements MBMesquite::MeshDomain.

Definition at line 1221 of file DomainClassifier.cpp.

References find_vertex_domain().

{
    if( const MeshDomain* dom = find_vertex_domain( handle ) )
        dom->closest_point( handle, position, closest, normal, err );
}

void MBMesquite::DomainClassifier::delete_all_sub_domains ( )

Definition at line 1253 of file DomainClassifier.cpp.

References clear(), domains, elementList, and vertexList.

Referenced by ~DomainClassifier().

{
    // get unique list of domains
    std::set< MeshDomain* > domains;
    std::vector< DomainBlock >::iterator i;
    for( i = vertexList.begin(); i != vertexList.end(); ++i )
        domains.insert( i->domain );
    for( i = elementList.begin(); i != elementList.end(); ++i )
        domains.insert( i->domain );
    std::set< MeshDomain* >::iterator j;
    for( j = domains.begin(); j != domains.end(); ++j )
        delete *j;
    clear();
}

MESQUITE_EXPORT void MBMesquite::DomainClassifier::delete_sub_domains ( bool yesno ) [inline]

Definition at line 237 of file DomainClassifier.hpp.

References deleteSubDomains.

Referenced by PatchDataTest::get_quad8_mesh_and_domain(), and SlaveBoundaryVerticesTest::make_mesh().

    {
        deleteSubDomains = yesno;
    }

void MBMesquite::DomainClassifier::domain_DoF	(	const Mesh::VertexHandle *	handle_array,
		unsigned short *	dof_array,
		size_t	num_handles,
		MsqError &	err
	)		const `[virtual]`

Definition at line 1231 of file DomainClassifier.cpp.

References MBMesquite::MeshDomain::domain_DoF(), find_vertex_domain(), and MSQ_ERRRTN.

Referenced by SlaveBoundaryVerticesTest::test_slaved_common().

{
    for( size_t i = 0; i < num_handles; ++i )
    {
        const MeshDomain* dom = find_vertex_domain( handles[i] );
        if( !dom )
            dof_array[i] = 3;
        else
        {
            dom->domain_DoF( handles + i, dof_array + i, 1, err );MSQ_ERRRTN( err );
        }
    }
}

void MBMesquite::DomainClassifier::element_normal_at	(	Mesh::ElementHandle	entity_handle,
		Vector3D &	coordinate
	)		const `[virtual]`

Implements MBMesquite::MeshDomain.

Definition at line 1198 of file DomainClassifier.cpp.

References find_element_domain().

{
    if( const MeshDomain* dom = find_element_domain( entity_handle ) )
        dom->element_normal_at( entity_handle, coordinate );
}

MeshDomain * MBMesquite::DomainClassifier::find_domain	(	Mesh::EntityHandle	handle,
		const std::vector< DomainBlock > &	list
	)		`[static, private]`

Definition at line 1179 of file DomainClassifier.cpp.

Referenced by find_element_domain(), and find_vertex_domain().

{
    std::vector< DomainClassifier::DomainBlock >::const_iterator i;
    i = std::lower_bound( list.begin(), list.end(), handle );
    return ( i != list.end() && i->firstHandle <= handle ) ? i->domain : NULL;
}

MESQUITE_EXPORT MeshDomain* MBMesquite::DomainClassifier::find_element_domain ( Mesh::ElementHandle element ) const [inline]

Definition at line 231 of file DomainClassifier.hpp.

References elementList, and find_domain().

Referenced by DomainClassifierTest::check_domain(), and element_normal_at().

    {
        return find_domain( element, elementList );
    }

MESQUITE_EXPORT MeshDomain* MBMesquite::DomainClassifier::find_vertex_domain ( Mesh::VertexHandle vertex ) const [inline]

Definition at line 226 of file DomainClassifier.hpp.

References find_domain(), and vertexList.

Referenced by DomainClassifierTest::check_domain(), closest_point(), domain_DoF(), snap_to(), and vertex_normal_at().

    {
        return find_domain( vertex, vertexList );
    }

void MBMesquite::DomainClassifier::snap_to	(	Mesh::VertexHandle	entity_handle,
		Vector3D &	coordinate
	)		const `[virtual]`

Modifies "coordinate" so that it lies on the domain to which "entity_handle" is constrained. The handle determines the domain. The coordinate is the proposed new position on that domain.

Implements MBMesquite::MeshDomain.

Definition at line 1187 of file DomainClassifier.cpp.

References find_vertex_domain().

{
    if( const MeshDomain* dom = find_vertex_domain( entity_handle ) ) dom->snap_to( entity_handle, coordinate );
}

void MBMesquite::DomainClassifier::test_valid_classification	(	Mesh *	mesh,
		MsqError &	err
	)

Check that classification maps to B-Rep topology.

Definition at line 842 of file DomainClassifier.cpp.

References MBMesquite::arrptr(), dim, dof(), domains, elementList, MBMesquite::Mesh::elements_get_attached_vertices(), MBMesquite::Mesh::elements_get_topologies(), MBMesquite::faces, MBMesquite::Mesh::get_all_elements(), MBMesquite::Mesh::get_all_vertices(), MBMesquite::MsqError::INVALID_STATE, MSQ_ERRRTN, MSQ_SETERR, n, MBMesquite::next_vertex(), vertexList, MBMesquite::Mesh::vertices_get_attached_elements(), and vtx().

Referenced by DomainClassifierTest::test_valid_classification().

{
    size_t i;

    // Get all mesh entities
    std::vector< Mesh::VertexHandle > vertices;
    std::vector< Mesh::ElementHandle > elements;
    mesh->get_all_vertices( vertices, err );MSQ_ERRRTN( err );
    mesh->get_all_elements( elements, err );MSQ_ERRRTN( err );
    std::sort( vertices.begin(), vertices.end() );
    std::sort( elements.begin(), elements.end() );

    // Get contents of each domain.
    std::map< MeshDomain*, int >::iterator iter;
    std::map< MeshDomain*, int > idxmap;
    for( i = 0; i < vertexList.size(); ++i )
        idxmap[vertexList[i].domain] = 0;
    for( i = 0; i < elementList.size(); ++i )
        idxmap[elementList[i].domain] = 0;
    std::vector< DomainSet > domains( idxmap.size() );
    int idx = 0;
    for( iter = idxmap.begin(); iter != idxmap.end(); ++iter )
    {
        iter->second        = idx;
        domains[idx].domain = iter->first;
        ++idx;
    }
    for( i = 0; i < vertexList.size(); ++i )
    {
        std::vector< Mesh::VertexHandle >::const_iterator s, e;
        s              = std::lower_bound( vertices.begin(), vertices.end(), vertexList[i].firstHandle );
        e              = std::upper_bound( vertices.begin(), vertices.end(), vertexList[i].lastHandle );
        DomainSet& set = domains[idxmap[vertexList[i].domain]];
        std::copy( s, e, std::back_inserter( set.vertices ) );
    }
    for( i = 0; i < elementList.size(); ++i )
    {
        std::vector< Mesh::ElementHandle >::const_iterator s, e;
        s              = std::lower_bound( elements.begin(), elements.end(), elementList[i].firstHandle );
        e              = std::upper_bound( elements.begin(), elements.end(), elementList[i].lastHandle );
        DomainSet& set = domains[idxmap[elementList[i].domain]];
        std::copy( s, e, std::back_inserter( set.elements ) );
    }

    // guess geometric dimension for each domain
    std::map< DomainSet*, int > dimmap;
    std::vector< unsigned short > dof;
    for( i = 0; i < domains.size(); ++i )
    {
        if( !domains[i].elements.empty() )
            dimmap[&domains[i]] = 2;
        else
        {
            dof.resize( domains[i].vertices.size() );
            domains[i].domain->domain_DoF( &( domains[i].vertices[0] ), arrptr( dof ), dof.size(), err );MSQ_ERRRTN( err );
            unsigned short dim  = *std::max_element( dof.begin(), dof.end() );
            dimmap[&domains[i]] = dim;
        }
    }

    // group domains by dimension
    std::vector< DomainSet* > points, curves, surfaces;
    for( std::map< DomainSet*, int >::iterator it = dimmap.begin(); it != dimmap.end(); ++it )
    {
        switch( it->second )
        {
            case 0:
                points.push_back( it->first );
                break;
            case 1:
                curves.push_back( it->first );
                break;
            case 2:
                surfaces.push_back( it->first );
                break;
            default:
                MSQ_SETERR( err )
                ( MsqError::INVALID_STATE, "Invalid domain dimension: %d\n", it->second );
                break;
        }
    }

    // check that each point-domain has a single vertex
    for( i = 0; i < points.size(); ++i )
    {
        if( points[i]->vertices.size() != 1 || !points[i]->elements.empty() )
        {
            MSQ_SETERR( err )( "Point domain mesh not a single vertex.", MsqError::INVALID_STATE );
            return;
        }
    }

    // check that each curve domain has a chain of connected edges
    std::set< Mesh::VertexHandle > unseen;
    for( i = 0; i < curves.size(); ++i )
    {
        if( !curves[i]->elements.empty() )
        {
            MSQ_SETERR( err )( "Elements associated with 1D domain.", MsqError::INVALID_STATE );
            return;
        }

        unseen.clear();
        std::copy( curves[i]->vertices.begin(), curves[i]->vertices.end(), std::inserter( unseen, unseen.begin() ) );

        const Mesh::VertexHandle first_vtx = *unseen.begin();
        unseen.erase( unseen.begin() );
        Mesh::VertexHandle vtx = first_vtx;
        // find chain of vertices
        while( next_vertex( mesh, vtx, unseen, err ) )
            MSQ_ERRRTN( err );
        // search again from the starting vertex because
        // it probably wasn't the first one on the curve
        vtx = first_vtx;
        while( next_vertex( mesh, vtx, unseen, err ) )
            MSQ_ERRRTN( err );
        // were all vertices in a chain?
        if( !unseen.empty() )
        {
            MSQ_SETERR( err )
            ( "Curve domain contains vertices not in a simply connected chain.", MsqError::INVALID_STATE );
            return;
        }
    }

    std::set< Mesh::VertexHandle > seen;
    std::set< Mesh::ElementHandle > remaining;
    std::vector< Mesh::VertexHandle > verts, verts2;
    std::vector< Mesh::ElementHandle > stack, vert_elems;
    std::vector< EntityTopology > types;
    std::vector< size_t > junk;
    // if surface contains elements...
    for( i = 0; i < surfaces.size(); ++i )
    {
        if( surfaces[i]->elements.empty() ) continue;

        // Check that any vertices on surface are contained in an
        // element on the surface.
        verts.clear();
        mesh->elements_get_attached_vertices( &( surfaces[i]->elements[0] ), surfaces[i]->elements.size(), verts, junk,
                                              err );MSQ_ERRRTN( err );
        seen.clear();
        std::copy( verts.begin(), verts.end(), std::inserter( seen, seen.begin() ) );

        std::vector< Mesh::VertexHandle >::const_iterator v;
        for( v = surfaces[i]->vertices.begin(); v != surfaces[i]->vertices.end(); ++v )
        {
            std::set< Mesh::VertexHandle >::iterator j = seen.find( *v );
            if( j == seen.end() )
            {
                MSQ_SETERR( err )
                ( "Vertex on surface domain not in any element.", MsqError::INVALID_STATE );
                return;
            }
        }

        // check that elements form 2D patch
        stack.clear();
        remaining.clear();
        std::copy( surfaces[i]->elements.begin(), surfaces[i]->elements.end(),
                   std::inserter( remaining, remaining.begin() ) );
        stack.push_back( *remaining.begin() );
        remaining.erase( remaining.begin() );
        while( !stack.empty() )
        {
            Mesh::ElementHandle elem = stack.back();
            stack.pop_back();
            verts.clear();
            mesh->elements_get_attached_vertices( &elem, 1, verts, junk, err );MSQ_ERRRTN( err );
            // for each edge
            for( size_t j = 0; j < verts.size(); ++j )
            {
                Mesh::VertexHandle v1 = verts[j], v2 = verts[( j + 1 ) % verts.size()];
                vert_elems.clear();
                mesh->vertices_get_attached_elements( &v1, 1, vert_elems, junk, err );MSQ_ERRRTN( err );
                types.resize( vert_elems.size() );
                if( !vert_elems.empty() )
                {
                    mesh->elements_get_topologies( arrptr( vert_elems ), arrptr( types ), vert_elems.size(), err );MSQ_ERRRTN( err );
                }
                while( !vert_elems.empty() )
                {
                    Mesh::ElementHandle e2 = vert_elems.back();
                    EntityTopology type    = types.back();
                    vert_elems.pop_back();
                    types.pop_back();
                    if( TopologyInfo::dimension( type ) != 2 ) continue;
                    verts2.clear();
                    mesh->elements_get_attached_vertices( &e2, 1, verts2, junk, err );MSQ_ERRRTN( err );
                    size_t jdx = std::find( verts2.begin(), verts2.end(), v1 ) - verts2.begin();
                    if( verts2[( jdx + 1 ) % verts2.size()] != v2 &&
                        verts2[( jdx + verts2.size() - 1 ) % verts2.size()] != v2 )
                        continue;
                    std::set< Mesh::ElementHandle >::iterator r = remaining.find( e2 );
                    if( r == remaining.end() ) continue;

                    stack.push_back( *r );
                    remaining.erase( r );
                }
            }
        }

        if( !remaining.empty() )
        {
            MSQ_SETERR( err )
            ( "Surface mesh not a single, simply-connected patch", MsqError::INVALID_STATE );
            return;
        }
    }

    // check that sides of volume elements that are on surface
    // form simply connected patch
    for( i = 0; i < surfaces.size(); ++i )
    {
        // build list of sides known to be on surface
        std::sort( surfaces[i]->vertices.begin(), surfaces[i]->vertices.end() );
        std::set< std::pair< Mesh::ElementHandle, int > > sides;
        for( size_t j = 0; j < surfaces[i]->vertices.size(); ++j )
        {
            Mesh::VertexHandle v = surfaces[i]->vertices[j];
            vert_elems.clear();
            mesh->vertices_get_attached_elements( &v, 1, vert_elems, junk, err );MSQ_ERRRTN( err );
            types.resize( vert_elems.size() );
            if( !vert_elems.empty() )
            {
                mesh->elements_get_topologies( arrptr( vert_elems ), arrptr( types ), vert_elems.size(), err );MSQ_ERRRTN( err );
            }
            while( !vert_elems.empty() )
            {
                Mesh::ElementHandle e = vert_elems.back();
                EntityTopology type   = types.back();
                vert_elems.pop_back();
                types.pop_back();
                if( TopologyInfo::dimension( type ) != 3 ) continue;
                verts.clear();
                mesh->elements_get_attached_vertices( &e, 1, verts, junk, err );MSQ_ERRRTN( err );

                for( unsigned s = 0; s < TopologyInfo::faces( type ); ++s )
                {
                    unsigned n;
                    const unsigned* si = TopologyInfo::face_vertices( type, s, n );
                    unsigned ns        = 0;
                    for( unsigned k = 0; k < n; ++k )
                    {
                        if( std::binary_search( surfaces[i]->vertices.begin(), surfaces[i]->vertices.end(),
                                                verts[si[k]] ) )
                            ++ns;
                    }
                    if( ns >= 3 ) sides.insert( std::pair< Mesh::ElementHandle, int >( e, s ) );
                }
            }
        }

        std::vector< std::pair< Mesh::ElementHandle, int > > sstack;
        sstack.push_back( *sides.begin() );
        sides.erase( sides.begin() );
        while( !sstack.empty() )
        {
            Mesh::ElementHandle e = sstack.back().first;
            int s                 = sstack.back().second;
            sstack.pop_back();

            verts.clear();
            mesh->elements_get_attached_vertices( &e, 1, verts, junk, err );MSQ_ERRRTN( err );
            EntityTopology type;
            mesh->elements_get_topologies( &e, &type, 1, err );MSQ_ERRRTN( err );
            unsigned n;
            const unsigned* si = TopologyInfo::face_vertices( type, s, n );

            // for each edge
            for( unsigned j = 0; j < n; ++j )
            {
                Mesh::VertexHandle v1 = verts[si[j]], v2 = verts[si[( j + 1 ) % n]];
                vert_elems.clear();
                mesh->vertices_get_attached_elements( &v1, 1, vert_elems, junk, err );MSQ_ERRRTN( err );
                types.resize( vert_elems.size() );
                if( !vert_elems.empty() )
                {
                    mesh->elements_get_topologies( arrptr( vert_elems ), arrptr( types ), vert_elems.size(), err );MSQ_ERRRTN( err );
                }
                while( !vert_elems.empty() )
                {
                    Mesh::ElementHandle e2 = vert_elems.back();
                    EntityTopology type2   = types.back();
                    vert_elems.pop_back();
                    types.pop_back();
                    if( TopologyInfo::dimension( type ) != 3 ) continue;
                    verts2.clear();
                    mesh->elements_get_attached_vertices( &e2, 1, verts2, junk, err );MSQ_ERRRTN( err );
                    // for each face
                    for( unsigned s2 = 0; s2 < TopologyInfo::faces( type2 ); ++s2 )
                    {
                        std::set< std::pair< Mesh::ElementHandle, int > >::iterator side;
                        side = sides.find( std::pair< Mesh::ElementHandle, int >( e2, s2 ) );
                        if( side == sides.end() ) continue;

                        unsigned n2;
                        const unsigned* si2 = TopologyInfo::face_vertices( type2, s2, n2 );
                        unsigned jdx;
                        for( jdx = 0; jdx < n2; ++jdx )
                            if( verts2[si2[jdx]] == v1 ) break;
                        assert( jdx < n2 );

                        if( verts2[si2[( jdx + 1 ) % n2]] == v2 || verts2[si2[( jdx + n2 - 1 ) % n2]] == v2 )
                        {
                            sstack.push_back( *side );
                            sides.erase( side );
                        }
                    }
                }
            }
        }

        if( !sides.empty() )
        {
            MSQ_SETERR( err )
            ( "Surface mesh not a single, simply-connected patch", MsqError::INVALID_STATE );
            return;
        }
    }
}

void MBMesquite::DomainClassifier::vertex_normal_at	(	Mesh::VertexHandle	entity_handle,
		Vector3D &	coordinate
	)		const `[virtual]`

Returns the normal of the domain to which "entity_handle" is constrained. For non-planar surfaces, the normal is calculated at the point on the domain that is closest to the passed in value of "coordinate". If the domain does not have a normal, or the normal cannot be determined, "coordinate" is set to (0,0,0). Otherwise, "coordinate" is set to the domain's normal at the appropriate point. In summary, the handle determines the domain. The coordinate determines the point of interest on that domain.

Implements MBMesquite::MeshDomain.

Definition at line 1192 of file DomainClassifier.cpp.

References find_vertex_domain().

{
    if( const MeshDomain* dom = find_vertex_domain( entity_handle ) )
        dom->vertex_normal_at( entity_handle, coordinate );
}

void MBMesquite::DomainClassifier::vertex_normal_at	(	const Mesh::VertexHandle *	handles,
		Vector3D	coordinates[],
		unsigned	count,
		MsqError &	err
	)		const `[virtual]`

evaluate surface normals

Returns normals for a domain.

Parameters:

handles	The domain evaluated is the one in which this mesh entity is constrained.
coordinates	As input, a list of positions at which to evaluate the domain. As output, the resulting domain normals.
count	The length of the coordinates array.

Implements MBMesquite::MeshDomain.

Definition at line 1204 of file DomainClassifier.cpp.

References find_vertex_domain(), MBMesquite::MsqError::INVALID_ARG, MSQ_ERRRTN, MSQ_SETERR, and MBMesquite::MeshDomain::vertex_normal_at().

{
    for( unsigned i = 0; i < count; ++i )
    {
        const MeshDomain* dom = find_vertex_domain( handles[i] );
        if( !dom )
        {
            MSQ_SETERR( err )( MsqError::INVALID_ARG );
            return;
        }
        dom->vertex_normal_at( handles + i, coordinates + i, 1, err );MSQ_ERRRTN( err );
    }
}