Skip to content
Snippets Groups Projects
Commit 383e7e00 authored by Stefen Guzik's avatar Stefen Guzik
Browse files

Only 2D elements have been tested
parent a2f4ff59
No related branches found
No related tags found
No related merge requests found
#ifndef _MNEIGHBOUR_H_
#define _MNEIGHBOUR_H_
// Copyright (C) 2006 S. Guzik, C. Geuzaine, J.-F. Remacle
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
// USA.
//
// Please report all bugs and problems to <gmsh@geuz.org>.
/*******************************************************************************
*
* - The classes in this file construct a database of the elements that share
* the lower-dimensional bounding objects (e.g., vertex, edge, or face).
* These lower-dimensional objects are referred to in general as polytopes.
* - Templated traits classes are used extensively to define the characteristics
* of the elements and the entities.
*
* ****************************************************************************/
#include <algorithm>
#include <iterator>
#include <vector>
#include "MElement.h"
#include "GEdge.h"
#include "GFace.h"
#include "GRegion.h"
#include "GmshDefines.h"
// #define HAVE_HASH_MAP
#if defined(HAVE_HASH_MAP)
#include "HashMap.h"
#endif
/*==============================================================================
* File scope types
*============================================================================*/
typedef std::list<MElement*> Neighbours;
typedef Neighbours::const_iterator NeighboursConstIterator;
typedef Neighbours::iterator NeighboursIterator;
struct Range_t
{
unsigned int num;
NeighboursIterator begin;
Range_t() : num(0) { }
};
//--Use a hash map for neighbour lookup if possible, otherwise a map will do
#if defined(HAVE_HASH_MAP)
struct Hash_VertexPtr : public std::unary_function<MVertex*, size_t>
{
size_t operator()(const MVertex *const v) const
{
return HashFNV1a<sizeof(MVertex*)>::eval(v);
}
};
typedef HASH_MAP<MVertex*, Range_t, Hash_VertexPtr, std::equal_to<MVertex*> >
VertNRange;
typedef HASH_MAP<MEdge, Range_t, Hash_Edge, Equal_Edge> EdgeNRange;
typedef HASH_MAP<MFace, Range_t, Hash_Face, Equal_Face> FaceNRange;
#else
typedef std::map<MVertex*, Range_t, std::less<MVertex*> > VertNRange;
typedef std::map<MEdge, Range_t, Less_Edge> EdgeNRange;
typedef std::map<MFace, Range_t, Less_Face> FaceNRange;
#endif
typedef VertNRange::iterator VertNRangeIterator;
typedef VertNRange::const_iterator VertNRangeConstIterator;
typedef EdgeNRange::iterator EdgeNRangeIterator;
typedef EdgeNRange::const_iterator EdgeNRangeConstIterator;
typedef FaceNRange::iterator FaceNRangeIterator;
typedef FaceNRange::const_iterator FaceNRangeConstIterator;
/*==============================================================================
* Traits classes - that just perform compile-time checks for valid argument
* types
*============================================================================*/
//--This traits class checks for a pointer and strips it
// NOTE: If the compiler sent you here, you should be using a pointer.
template <typename T> struct StripPointer;
template <typename T> struct StripPointer<T*> { typedef T Type; };
//--This traits class checks for valid types of entities and iterators when the
//--entity is explicitly specified
// NOTE: If the compiler sent you here, you're using an invalid entity and/or
// invalid iterator. Valid pairs are shown below
template <typename Ent, typename Iter> struct ValidSpecEntityIterator;
template <> struct ValidSpecEntityIterator<GEdge, GEdge*>
{ typedef void Type; };
template <> struct ValidSpecEntityIterator<GEdge, GEntity*>
{ typedef void Type; };
template <> struct ValidSpecEntityIterator<GFace, GFace*>
{ typedef void Type; };
template <> struct ValidSpecEntityIterator<GFace, GEntity*>
{ typedef void Type; };
template <> struct ValidSpecEntityIterator<GRegion, GRegion*>
{ typedef void Type; };
template <> struct ValidSpecEntityIterator<GRegion, GEntity*>
{ typedef void Type; };
//--This traits class checks for valid types of entities and iterators when the
//--entity is not specified. It also strips the pointer.
// NOTE: If the compiler sent you here, you're using an invalid entity or
// iterator of entities
template <typename Iter> struct ValidEntityIterator;
template <> struct ValidEntityIterator<GEdge*>
{ typedef GEdge Type; };
template <> struct ValidEntityIterator<GFace*>
{ typedef GFace Type; };
template <> struct ValidEntityIterator<GRegion*>
{ typedef GRegion Type; };
//--This traits class checks for valid element types
// NOTE: If the compiler sent you here, you're using an invalid iterator for an
// element. What's listed here are NOT the only valid types. Any iterator with
// a value type of a pointer to an element should work.
template <typename Elem_o1> struct ValidElement;
template <> struct ValidElement<MElement> { typedef void Type; };
template <> struct ValidElement<MLine> { typedef void Type; };
template <> struct ValidElement<MTriangle> { typedef void Type; };
template <> struct ValidElement<MQuadrangle> { typedef void Type; };
template <> struct ValidElement<MTetrahedron> { typedef void Type; };
template <> struct ValidElement<MHexahedron> { typedef void Type; };
template <> struct ValidElement<MPrism> { typedef void Type; };
template <> struct ValidElement<MPyramid> { typedef void Type; };
/*==============================================================================
* Traits classes - that return information about a type
*============================================================================*/
//--This traits class gives the first-order base type for a second-order element
template <typename Order2> struct ElemBaseOrder
{ typedef Order2 Order1; };
template <> struct ElemBaseOrder<MLine3>
{ typedef MLine Order1; };
template <> struct ElemBaseOrder<MTriangle6>
{ typedef MTriangle Order1; };
template <> struct ElemBaseOrder<MQuadrangle8>
{ typedef MQuadrangle Order1; };
template <> struct ElemBaseOrder<MQuadrangle9>
{ typedef MQuadrangle Order1; };
template <> struct ElemBaseOrder<MTetrahedron10>
{ typedef MTetrahedron Order1; };
template <> struct ElemBaseOrder<MHexahedron20>
{ typedef MHexahedron Order1; };
template <> struct ElemBaseOrder<MHexahedron27>
{ typedef MHexahedron Order1; };
template <> struct ElemBaseOrder<MPrism15>
{ typedef MPrism Order1; };
template <> struct ElemBaseOrder<MPrism18>
{ typedef MPrism Order1; };
template <> struct ElemBaseOrder<MPyramid13>
{ typedef MPyramid Order1; };
template <> struct ElemBaseOrder<MPyramid14>
{ typedef MPyramid Order1; };
//--This traits class describes describes the lower dimensional bounding
//--polytopes of a first-order element
template <typename Elem> struct ElemTr;
template <> struct ElemTr<MElement>
{ enum { numVertex = 0, numEdge = 0, numFace = 0 }; };
template <> struct ElemTr<MLine>
{ enum { numVertex = 2, numEdge = 0, numFace = 0 }; };
template <> struct ElemTr<MTriangle>
{ enum { numVertex = 3, numEdge = 3, numFace = 0 }; };
template <> struct ElemTr<MQuadrangle>
{ enum { numVertex = 4, numEdge = 4, numFace = 0 }; };
template <> struct ElemTr<MTetrahedron>
{ enum { numVertex = 4, numEdge = 6, numFace = 4 }; };
template <> struct ElemTr<MHexahedron>
{ enum { numVertex = 8, numEdge = 12, numFace = 6 }; };
template <> struct ElemTr<MPrism>
{ enum { numVertex = 6, numEdge = 9, numFace = 5 }; };
template <> struct ElemTr<MPyramid>
{ enum { numVertex = 5, numEdge = 8, numFace = 5 }; };
//--This traits class gives the number of element types in entity Ent
template <typename Ent> struct EntTr;
template <> struct EntTr<GEdge> { enum { numElemTypes = 1 }; };
template <> struct EntTr<GFace> { enum { numElemTypes = 2 }; };
template <> struct EntTr<GRegion> { enum { numElemTypes = 4 }; };
//--This traits class gives iterator types and begin and end iterators for
//--element type number N in entity Ent.
template <typename Ent, unsigned int N> struct EntElemTr;
template <> struct EntElemTr<GEdge, 1> {
typedef MLine Elem;
typedef std::vector<MLine*>::const_iterator ConstElementIterator;
typedef std::vector<MLine*>::iterator ElementIterator;
static ConstElementIterator begin(const GEdge *const edge)
{
return edge->lines.begin();
}
static ConstElementIterator end(const GEdge *const edge)
{
return edge->lines.end();
}
};
template <> struct EntElemTr<GFace, 1> {
typedef MQuadrangle Elem;
typedef std::vector<MQuadrangle*>::const_iterator ConstElementIterator;
typedef std::vector<MQuadrangle*>::iterator ElementIterator;
static ConstElementIterator begin(const GFace *const face)
{
return face->quadrangles.begin();
}
static ConstElementIterator end(const GFace *const face)
{
return face->quadrangles.end();
}
};
template <> struct EntElemTr<GFace, 2> {
typedef MTriangle Elem;
typedef std::vector<MTriangle*>::const_iterator ConstElementIterator;
typedef std::vector<MTriangle*>::iterator ElementIterator;
static ConstElementIterator begin(const GFace *const face)
{
return face->triangles.begin();
}
static ConstElementIterator end(const GFace *const face)
{
return face->triangles.end();
}
};
template <> struct EntElemTr<GRegion, 1> {
typedef MPyramid Elem;
typedef std::vector<MPyramid*>::const_iterator ConstElementIterator;
typedef std::vector<MPyramid*>::iterator ElementIterator;
static ConstElementIterator begin(const GRegion *const region)
{
return region->pyramids.begin();
}
static ConstElementIterator end(const GRegion *const region)
{
return region->pyramids.end();
}
};
template <> struct EntElemTr<GRegion, 2> {
typedef MPrism Elem;
typedef std::vector<MPrism*>::const_iterator ConstElementIterator;
typedef std::vector<MPrism*>::iterator ElementIterator;
static ConstElementIterator begin(const GRegion *const region)
{
return region->prisms.begin();
}
static ConstElementIterator end(const GRegion *const region)
{
return region->prisms.end();
}
};
template <> struct EntElemTr<GRegion, 3> {
typedef MHexahedron Elem;
typedef std::vector<MHexahedron*>::const_iterator ConstElementIterator;
typedef std::vector<MHexahedron*>::iterator ElementIterator;
static ConstElementIterator begin(const GRegion *const region)
{
return region->hexahedra.begin();
}
static ConstElementIterator end(const GRegion *const region)
{
return region->hexahedra.end();
}
};
template <> struct EntElemTr<GRegion, 4> {
typedef MTetrahedron Elem;
typedef std::vector<MTetrahedron*>::const_iterator ConstElementIterator;
typedef std::vector<MTetrahedron*>::iterator ElementIterator;
static ConstElementIterator begin(const GRegion *const region)
{
return region->tetrahedra.begin();
}
static ConstElementIterator end(const GRegion *const region)
{
return region->tetrahedra.end();
}
};
//--This is a traits/policy class for the lower-dimension polytopes that bound
//--an element (i.e., vertex, edge, or face). It returns the corresponding
//--range type and extracts from the element the function to get the polytope.
//--Note that MVertex is a pointer because it exists somewhere whereas MEdge and
//--MFace are constructed using MVertex's by getEdge() and getFace().
template <typename Polytope> struct PolytopeTr;
template <> struct PolytopeTr<MVertex*> {
typedef VertNRange PolytopeNRange;
typedef VertNRangeConstIterator PNRConstIterator;
template <typename Elem>
static MVertex *getPolytope(Elem *const element,
const unsigned int nPolytope)
{
return element->getVertex(nPolytope);
}
};
template <> struct PolytopeTr<MEdge> {
typedef EdgeNRange PolytopeNRange;
typedef EdgeNRangeConstIterator PNRConstIterator;
template <typename Elem>
static MEdge getPolytope(Elem *const element,
const unsigned int nPolytope)
{
return element->getEdge(nPolytope);
}
};
template <> struct PolytopeTr<MFace> {
typedef FaceNRange PolytopeNRange;
typedef FaceNRangeConstIterator PNRConstIterator;
template <typename Elem>
static MFace getPolytope(Elem *const element,
const unsigned int nPolytope)
{
return element->getFace(nPolytope);
}
};
/*******************************************************************************
*
* class: PolytopeIterator
*
* Purpose
* =======
*
* Provides an iterator to iterate through the unique vertices, edges and,
* faces defined by the (hash_)maps. An iterator to the (hash_)map has a
* value type of a (key, value) pair. This iterator makes the container look
* like it only contains the bounding polytope.
*
* Constructors
* ============
*
* The class takes one template argument <Polytope> which should be either
* MVertex*, MEdge, or MFace.
*
* PolytopeIterator()
* -- default constructor
*
* PolytopeIterator(const BaseIterator &baseIter_in)
* -- the base iterator is the real iterator to the
* (hash_)map
*
* PolytopeIterator(const PolytopeIterator &iter)
* -- copy
*
* Destructor
* ==========
*
* ~PolytopeIterator -- synthesized
*
* Member Functions
* ================
*
* unsigned int num_neighbours()
* -- returns number of elements sharing the polytope
*
* Operators
* =========
*
* Includes typical bidirectional iterator operators {=, ==, !=, *, ->, ++(),
* --(), ()++, ()--} with * and -> dereferencing to the polytope.
*
* Notes
* =====
*
* - Only constant iterators are supported.
*
******************************************************************************/
template<typename Polytope>
class PolytopeIterator // Actually only a const_iterator
{
public:
// The base iterator iterates through the (hash_)maps defining the range of
// element neighbours for a polytope
typedef typename PolytopeTr<Polytope>::PNRConstIterator BaseIterator;
typedef PolytopeIterator Self;
// Standard traits
typedef std::bidirectional_iterator_tag iterator_category;
typedef Polytope value_type;
typedef BaseIterator difference_type;
typedef const Polytope* pointer;
typedef const Polytope& reference;
//--Constructors
PolytopeIterator() : baseIter() { }
PolytopeIterator(const BaseIterator &baseIter_in) : baseIter(baseIter_in) { }
PolytopeIterator(const PolytopeIterator &iter) : baseIter(iter.baseIter) { }
PolytopeIterator& operator=(const PolytopeIterator& iter)
{
if(&iter != this) baseIter = iter.baseIter;
return *this;
}
//--Comparison
bool operator==(const Self& iter) const { return baseIter == iter.baseIter; }
bool operator!=(const Self& iter) const { return baseIter != iter.baseIter; }
//--Dereference and arrow operator
reference operator*() const { return baseIter->first; }
pointer operator->() const { return &baseIter->first; }
//--Increment
// Prefix
Self& operator++()
{
++baseIter;
return *this;
}
// Postfix
Self operator++(int)
{
Self tmp = *this;
++*this;
return tmp;
}
//--Decrement
// Prefix
Self& operator--()
{
--baseIter;
return *this;
}
// Postfix
Self operator--(int)
{
Self tmp = *this;
--*this;
return tmp;
}
//--Special
unsigned int num_neighbours() { return baseIter->second.num; }
private:
NeighboursConstIterator begin_neighbours() { return baseIter->second.begin; }
//--Date members
BaseIterator baseIter;
//--Friends
friend class MNeighbour;
};
/*******************************************************************************
*
* class: MNeighbour
*
* Purpose
* =======
*
* Generates a database for quick lookup of elements sharing a
* MVertex*, MEdge, or MFace (referred to in general as polytopes. I.e., a
* geometrical construct bounding the element with a lower dimension than the
* element). The neighbour elements are stored in a list. Maps or hash maps
* provide an index into the list for a given polytope. In addition to the
* element neighbours, this class provides iterators to the unique vertices,
* edges, and faces in the database.
*
* Iterators
* =========
*
* Vertex_const_iterator
* -- iterator to sequence of unique vertices
* Edge_const_iterator
* -- iterator to sequence of unique edges
* Face_const_iterator
* -- iterator to sequence of unique faces
*
* Constructors
* ============
*
* MNeighbour() -- default constructor does nothing
*
* Destructor
* ==========
*
* ~MNeighbour() -- synthesized
*
* Member Functions
* ================
*
* Iterators
* ---------
*
* Vertex_const_iterator vertex_begin()
* -- first unique vertex
*
* Vertex_const_iterator vertex_end()
* -- one-past-last unique vertex
*
* Edeg_const_iterator edge_begin()
* -- first unique edge
*
* Edge_const_iterator edge_end()
* -- one-past-last unique edge
*
* Face_const_iterator face_begin()
* -- first unique face
*
* Face_const_iterator face_end()
* -- one-past-last unique face
*
* Routines for adding elements
* ----------------------------
*
* The routines that generate a database of neighbours from container of
* entities (add_elements_in_entitie) require knowledge of the specific type
* of entity. In cases where a container contains specific entities, the type
* need not be stated. In cases where the container contains GEntity*, the
* user must state the actual entity type. Extensive compile-time type
* checking is employed to make sure the routines are used correctly.
*
* template <typename Ent, typename EntIter>
* void add_elements_in_entities(EntIter begin, EntIter end)
* -- adds elements in entities to neighbours database.
* Entity type <Ent> must be explicitly specified. The
* value type of the iterators must be a pointer to an
* entity. It is expected, but not required to be a
* pointer to a GEntity.
* begin (I) iterator/pointer to first entity pointer
* end (I) iterator/pointer to one-past-last entity pointer
*
* template <typename Ent, typename EntP>
* void add_elements_in_entities(EntP **pEnt, unsigned int n)
* -- adds elements in entities to neighbours database.
* Entity type <Ent> must be explicitly specified. The
* array element must be a pointer to an element. It is
* expected, but not required to be a pointer to a
* GEntity.
* pEnt (I) first entity in the array
* n (I) number of entities in the array
*
* template <typename EntIter>
* void add_elements_in_entities(EntIter begin, EntIter end)
* -- adds elements in entities to neighbours database.
* The entity type must not be GEntity. The value type
* of the iterators must be a pointer to an entity.
* begin (I) iterator/pointer to first entity pointer
* end (I) iterator/pointer to one-past-last entity pointer
*
* template <typename Ent>
* void add_elements_in_entities(Ent **pEnt, unsigned int n)
* -- adds elements in entities to neighbours database.
* The entity type must not be GEntity. The array
* element must be a pointer to an element.
* pEnt (I) first entity in the array
* n (I) number of entities in the array
*
* template <typename EntPtr>
* void add_elements_in_entity(EntPtr entity)
* -- adds elements in a single entity to the neighbours
* database. The entity type must not be GEntity. The
* dereferenced quantity must be a pointer to an entity.
* entity (I) pointer to an entity
*
* The routines that generate a database of neighbours from containers of
* elements do not require the user to state the specific type of element.
* If the container constains a specific element pointer, compiler-time
* polymorphism is used. If the container contains the base pointers
* MElement*, run-time switching based on element introspection is used.
*
* template <typename ElemIter>
* void add_elements(ElemIter begin, ElemIter end)
* -- adds elements from iterator sequence to the neighbours
* database.
* begin (I) iterator/pointer to first element pointer
* end (I) iterator/pointer to one-past-last element pointer
*
* template <typename Elem>
* void add_elements(Elem **pElem, const unsigned int n)
* -- adds elements in an array to the neighbours database.
* pElem (I) first element in the array
* n (I) number of elements in the array
*
* Routines for querying the neighbours database
* ---------------------------------------------
*
* unsigned int max_vertex_neighbours()
* -- returns max number of elements sharing a vertex. This
* is also the max number of elements sharing any
* bounding polytope.
*
* unsigned int max_edge_neighbours()
* -- returns max number of elements sharing an edge
*
* unsigned int max_face_neighbours()
* -- returns max number of elements sharing a face
*
* int vertex_neighbours(Vertex_const_iterator &itVert,
* std::vector<MElement*> &vElem)
* -- return all elements sharing a vertex
* itVert (I) iterator to a unique vertex. Avoids a find()
* vElem (O) vector with elements loaded via push_back
* return number of neighbour elements
*
* int vertex_neighbours(Vertex_const_iterator &itVert, MElement **pElem)
* -- return all elements sharing a vertex
* itVert (I) iterator to a unique vertex. Avoids a find()
* pElem (O) array with elements loaded starting at pElem[0]
* return number of neighbour elements
*
* int vertex_neighbours(MVertex *const vertex, std::vector<MElement*> &vElem)
* -- returns all elements sharing a vertex
* vertex (I) pointer to vertex. It must be found in (hash_)map
* vElem (O) vector with elements loaded via push_back
* return number of neighbour elements. 0 if vertex not found
*
* int vertex_neighbours(MVertex *const vertex, MElement **pElem)
* -- returns all elements sharing a vertex
* vertex (I) pointer to vertex. It must be found in (hash_)map
* pElem (O) array with elements loaded starting at pElem[0]
* return number of neighbour elements. 0 if vertex not found
*
* void edge_neighbours(Edge_const_iterator &itEdge,
* std::vector<MElement*> &vElem)
* -- return all elements sharing a edge
* itEdge (I) iterator to a unique edge. Avoids a find()
* vElem (O) vector with elements loaded via push_back
* return number of neighbour elements
*
* void edge_neighbours(Edge_const_iterator &itEdge, MElement **pElem)
* -- return all elements sharing a edge
* itEdge (I) iterator to a unique edge. Avoids a find()
* pElem (O) array with elements loaded starting at pElem[0]
* return number of neighbour elements
*
* int edge_neighbours(MEdge *const edge, std::vector<MElement*> &vElem)
* -- returns all elements sharing a edge
* edge (I) pointer to edge. It must be found in (hash_)map
* vElem (O) vector with elements loaded via push_back
* return number of neighbour elements. 0 if edge not found
*
* int edge_neighbours(MEdge *const edge, MElement **pElem)
* -- returns all elements sharing a edge
* edge (I) pointer to edge. It must be found in (hash_)map
* pElem (O) array with elements loaded starting at pElem[0]
* return number of neighbour elements. 0 if edge not found
*
* void face_neighbours(Face_const_iterator &itFace,
* std::vector<MElement*> &vElem)
* -- return all elements sharing a face
* itFace (I) iterator to a unique face. Avoids a find()
* vElem (O) vector with elements loaded via push_back
* return number of neighbour elements
*
* void face_neighbours(Face_const_iterator &itFace, MElement **pElem)
* -- return all elements sharing a face
* itFace (I) iterator to a unique face. Avoids a find()
* pElem (O) array with elements loaded starting at pElem[0]
* return number of neighbour elements
*
* int face_neighbours(MFace *const face, std::vector<MElement*> &vElem)
* -- returns all elements sharing a face
* face (I) pointer to face. It must be found in (hash_)map
* vElem (O) vector with elements loaded via push_back
* return number of neighbour elements. 0 if face not found
*
* int face_neighbours(MFace *const face, MElement **pElem)
* -- returns all elements sharing a face
* face (I) pointer to face. It must be found in (hash_)map
* pElem (O) array with elements loaded starting at pElem[0]
* return number of neighbour elements. 0 if face not found
*
* unsigned int vertex_num_neighbours(MVertex *const vertex) const
* -- returns the number of elements sharing a vertex
*
* unsigned int edge_num_neighbours(const MEdge& edge) const
* -- returns the number of elements sharing an edge
*
* unsigned int edge_num_neighbours(const MFace& face) const
* -- returns the number of elements sharing a face
*
* void clear() -- clears the neighbours database. Use this before
* constructing a new database of element neighbours.
* Otherwise, new elements are added to the existing
* database.
*
* MElement *element_neighbour(MElement *const element,
* const unsigned int iPolytope) const
* -- find the neighbour element across the d-1 bounding
* polytope. E.g., for a triangle, find an element
* across an edge.
* element (I) element to find
* iPolytope (I) which polytope to find. E.g., edge 1
* return the neighbour element or 0 if polytope not found
* in neighbours database.
*
* Example
* =======
*
*------------------------------------------------------------------------------
// Examples for adding elements from containers (the type of container doesn't
// matter.
MNeighbour meshNeighbour;
std::list<GEntity*> faceEnt1; // Assume known to be GFace type
meshNeighbour.add_elements_in_entities<GFace>(faceEnt1.begin(),
faceEnt1.end());
std::vector<GFace*> faceEnt2;
meshNeighbour.add_elements_in_entities(faceEnt2.begin(), faceEnt2.end());
meshNeighbour.add_elements_in_entity(faceEnt2.begin());
GEntity *faceEnt3[2]; // Assume known to be GFace type
meshNeighbour.add_elements_in_entities<GFace>(faceEnt3, 2);
GFace *faceEnt4[10];
meshNeighbour.add_elements_in_entities(faceEnt4, 10);
std::list<MElement*> elem1;
// Run time introspection. Elements can be mixed.
meshNeighbour.add_elements(elem1.begin(), elem1.end());
std::vector<MTriangle*> elem2;
// Compile time introspection.
meshNeighbour.add_elements(elem2.begin(), elem2.end());
MElement *elem3[5];
// Run time introspection. Elements can be mixed.
meshNeighbour.add_elements(elem3, 5);
MQuadrangle *elem4[6];
// Compile time introspection.
meshNeighbour.add_elements(elem4, 6);
*------------------------------------------------------------------------------
*
******************************************************************************/
class MNeighbour
{
/*==============================================================================
* Iterators
*============================================================================*/
public:
// Iterators over polytopes
typedef PolytopeIterator<MVertex*> Vertex_const_iterator;
typedef PolytopeIterator<MEdge> Edge_const_iterator;
typedef PolytopeIterator<MFace> Face_const_iterator;
// Begin and end points for these iterators
Vertex_const_iterator vertex_begin()
{
return Vertex_const_iterator(vertNRange.begin());
}
Vertex_const_iterator vertex_end()
{
return Vertex_const_iterator(vertNRange.end());
}
Edge_const_iterator edge_begin()
{
return Edge_const_iterator(edgeNRange.begin());
}
Edge_const_iterator edge_end()
{
return Edge_const_iterator(edgeNRange.end());
}
Face_const_iterator face_begin()
{
return Face_const_iterator(faceNRange.begin());
}
Face_const_iterator face_end()
{
return Face_const_iterator(faceNRange.end());
}
/*==============================================================================
* Template meta-programming classes
*============================================================================*/
private:
template <typename Elem, typename Polytope, unsigned int NP>
struct FindNeighbours;
template <typename Ent, unsigned int N = EntTr<Ent>::numElemTypes>
struct FindNeighboursInEntity;
/*==============================================================================
* Member functions and data
*============================================================================*/
public:
//--Default constructor.
// The constructor cannot really do anything because the initialization
// functions (find_neighbours_*) cannot perform template argument deduction.
MNeighbour() : maxVertNeighbours(0), maxEdgeNeighbours(0),
maxFaceNeighbours(0) { }
/*--------------------------------------------------------------------*
* Elements added from entities.
*--------------------------------------------------------------------*/
//--Compute database of neighbour elements. The entity pointers may be of a
//--base or specific type and the call must state the specific type.
template <typename Ent, typename EntIter>
void add_elements_in_entities(EntIter begin, EntIter end) {
// Check the Ent and EntIter types
// NOTE: If the compiler sent you here, you're using an invalid entity as a
// template parameter and/or and invalid iterator. See struct
// ValidEntityIterator above for valid types
typedef typename ValidSpecEntityIterator<Ent,
typename EntIter::value_type>::Type Check;
// Find the neighbours of each vertex, edge, and face
for(EntIter itEnt = begin; itEnt != end; ++itEnt) {
Ent *entity = static_cast<Ent*>(*itEnt);
FindNeighboursInEntity<Ent>::eval(entity, vertNRange, edgeNRange,
faceNRange, neighbours);
}
maxNeighbours();
}
//--Same as the routine above except the arguments are pointers instead of
//--iterators. Hence, value_type does not exist.
template <typename Ent, typename EntP>
void add_elements_in_entities(EntP **pEnt, unsigned int n) {
// Check the Ent and EntP types
// NOTE: If the compiler sent you here, you're using an invalid entity as a
// template parameter and/or invalid pointers. See struct
// ValidEntityIterator above for valid types
typedef typename ValidSpecEntityIterator<Ent, EntP*>::Type Check;
// Find the neighbours of each vertex, edge, and face
while(n--) {
Ent *entity = static_cast<Ent*>(*pEnt++);
FindNeighboursInEntity<Ent>::eval(entity, vertNRange, edgeNRange,
faceNRange, neighbours);
}
maxNeighbours();
}
//--Compute database of neighbour elements. The entity pointers must be of a
//--specific type (not a base class). The type of entity is determined from
//--Iterator::value_type and does not need to be stated.
template <typename EntIter>
void add_elements_in_entities(EntIter begin, EntIter end) {
// Check the value type in EntIter (it must not be GEntity*)
// NOTE: If the compiler sent you here, you're using an invalid entity as a
// template parameter and/or and invalid iterator for the constructor. See
// struct ValidEntityIterator above for valid types
typedef typename ValidEntityIterator<typename EntIter::value_type>::
Type Ent;
// Find the neighbours of each vertex, edge, and face
for(EntIter itEnt = begin; itEnt != end; ++itEnt) {
FindNeighboursInEntity<Ent>::eval(*itEnt, vertNRange, edgeNRange,
faceNRange, neighbours);
}
maxNeighbours();
}
//--Same as the routine above except the arguments are pointers instead of
//--iterators. Hence, value_type does not exist.
template <typename Ent>
void add_elements_in_entities(Ent **pEnt, unsigned int n) {
// Check the type Ent (it must not be GEntity)
// NOTE: If the compiler sent you here, you're using an invalid entity as a
// pointer See struct ValidEntityIterator above for valid types
typedef typename ValidEntityIterator<Ent*>::Type Check;
// Find the neighbours of each vertex, edge, and face
while(n--) {
FindNeighboursInEntity<Ent>::eval(*pEnt++, vertNRange, edgeNRange,
faceNRange, neighbours);
}
maxNeighbours();
}
//--Find neighbours from the elements attached to a pointer to a single entity.
//--The type of entity does not need to be stated.
template <typename EntPtr>
void add_elements_in_entity(EntPtr entity) {
// Check the type of EntPtr (it must not be GEntity*)
// NOTE: If the compiler sent you here, you're using an invalid entity as a
// pointer See struct ValidEntityIterator above for valid types
typedef typename ValidEntityIterator<EntPtr>::Type Ent;
// Find the neighbours of each vertex, edge, and face
FindNeighboursInEntity<Ent>::eval(entity, vertNRange, edgeNRange,
faceNRange, neighbours);
maxNeighbours();
}
/*--------------------------------------------------------------------*
* Elements added directly
*--------------------------------------------------------------------*/
//--Add elements from a container using iterators
template <typename ElemIter>
void add_elements(ElemIter begin, ElemIter end) {
// Strip the pointer to the element type and convert it to first order
typedef typename ElemBaseOrder<
typename StripPointer<
typename ElemIter::value_type
>::Type
>::Order1 Elem_1o;
// Check for a valid type of element
// NOTE: If the compiler sent you here, you're using an invalid iterator.
// The dereferenced object needs to be a pointer to an element.
typedef typename ValidElement<Elem_1o>::Type Check;
add_elements1<Elem_1o, ElemIter>(begin, end);
}
//--Add elements from an array of pointers (value_type does not exist)
template <typename Elem>
void add_elements(Elem **pElem, const unsigned int n) {
// Convert the element type to first order
typedef typename ElemBaseOrder<Elem>::Order1 Elem_1o;
// Check for a valid type of element
// NOTE: If the compiler sent you here, you're using an invalid pointer.
// The dereferenced object needs to be a pointer to an element.
typedef typename ValidElement<Elem_1o>::Type Check;
add_elements1<Elem_1o, Elem**>(pElem, pElem+n);
}
/*--------------------------------------------------------------------*
* Query the database
*--------------------------------------------------------------------*/
//--Get the max number of elements sharing a vertex, edge, or face
unsigned int max_vertex_neighbours() const { return maxVertNeighbours; }
unsigned int max_edge_neighbours() const { return maxEdgeNeighbours; }
unsigned int max_face_neighbours() const { return maxFaceNeighbours; }
//--Return elements that share a vertex
// From a vertex iterator
int vertex_neighbours(Vertex_const_iterator &itVert,
std::vector<MElement*> &vElem) const
{
NeighboursConstIterator itElem = itVert.begin_neighbours();
for(unsigned int n = itVert.num_neighbours(); n--;)
vElem.push_back(*itElem++);
return itVert.num_neighbours();
}
int vertex_neighbours(Vertex_const_iterator &itVert, MElement **pElem) const
{
NeighboursConstIterator itElem = itVert.begin_neighbours();
for(unsigned int n = itVert.num_neighbours(); n--;)
*pElem++ = *itElem++;
return itVert.num_neighbours();
}
// From a vertex
int vertex_neighbours(MVertex *const vertex,
std::vector<MElement*> &vElem) const
{
VertNRangeConstIterator itVert = vertNRange.find(vertex);
if(itVert == vertNRange.end()) return 0; // Vertex not found
NeighboursConstIterator itElem = itVert->second.begin;
for(unsigned int n = itVert->second.num; n--;) vElem.push_back(*itElem++);
return itVert->second.num;
}
int vertex_neighbours(MVertex *const vertex, MElement **pElem) const
{
VertNRangeConstIterator itVert = vertNRange.find(vertex);
if(itVert == vertNRange.end()) return 0; // Vertex not found
NeighboursConstIterator itElem = itVert->second.begin;
for(unsigned int n = itVert->second.num; n--;) *pElem++ = *itElem++;
return itVert->second.num;
}
//--Return elements that share an edge
// From an edge iterator
int edge_neighbours(Edge_const_iterator &itEdge,
std::vector<MElement*> &vElem) const
{
NeighboursConstIterator itElem = itEdge.begin_neighbours();
for(unsigned int n = itEdge.num_neighbours(); n--;)
vElem.push_back(*itElem++);
return itEdge.num_neighbours();
}
int edge_neighbours(Edge_const_iterator &itEdge, MElement **pElem) const
{
NeighboursConstIterator itElem = itEdge.begin_neighbours();
for(unsigned int n = itEdge.num_neighbours(); n--;)
*pElem++ = *itElem++;
return itEdge.num_neighbours();
}
// From an edge
int edge_neighbours(const MEdge& edge, std::vector<MElement*> &vElem) const
{
EdgeNRangeConstIterator itEdge = edgeNRange.find(edge);
if(itEdge == edgeNRange.end()) return 0; // Edge not found
NeighboursConstIterator itElem = itEdge->second.begin;
for(unsigned int n = itEdge->second.num; n--;) vElem.push_back(*itElem++);
return itEdge->second.num;
}
int edge_neighbours(const MEdge& edge, MElement **pElem) const
{
EdgeNRangeConstIterator itEdge = edgeNRange.find(edge);
if(itEdge == edgeNRange.end()) return 0; // Edge not found
NeighboursConstIterator itElem = itEdge->second.begin;
for(unsigned int n = itEdge->second.num; n--;) *pElem++ = *itElem++;
return itEdge->second.num;
}
//--Return elements that share a face
// From a face iterator
int face_neighbours(Face_const_iterator &itFace,
std::vector<MElement*> &vElem) const
{
NeighboursConstIterator itElem = itFace.begin_neighbours();
for(unsigned int n = itFace.num_neighbours(); n--;)
vElem.push_back(*itElem++);
return itFace.num_neighbours();
}
int face_neighbours(Face_const_iterator &itFace, MElement **pElem) const
{
NeighboursConstIterator itElem = itFace.begin_neighbours();
for(unsigned int n = itFace.num_neighbours(); n--;)
*pElem++ = *itElem++;
return itFace.num_neighbours();
}
// From a face
int face_neighbours(const MFace& face, std::vector<MElement*> &vElem) const
{
FaceNRangeConstIterator itFace = faceNRange.find(face);
if(itFace == faceNRange.end()) return 0; // Face not found
NeighboursConstIterator itElem = itFace->second.begin;
for(unsigned int n = itFace->second.num; n--;) vElem.push_back(*itElem++);
return itFace->second.num;
}
int face_neighbours(const MFace& face, MElement **pElem) const
{
FaceNRangeConstIterator itFace = faceNRange.find(face);
if(itFace == faceNRange.end()) return 0; // Face not found
NeighboursConstIterator itElem = itFace->second.begin;
for(unsigned int n = itFace->second.num; n--;) *pElem++ = *itElem++;
return itFace->second.num;
}
//--Return the number of elements that share a vertex
unsigned int vertex_num_neighbours(MVertex *const vertex) const
{
VertNRangeConstIterator itVert = vertNRange.find(vertex);
if(itVert == vertNRange.end()) return 0; // Vertex not found
return itVert->second.num;
}
//--Return the number of elements that share an edge
unsigned int edge_num_neighbours(const MEdge& edge) const
{
EdgeNRangeConstIterator itEdge = edgeNRange.find(edge);
if(itEdge == edgeNRange.end()) return 0; // Edge not found
return itEdge->second.num;
}
//--Return the number of elements that share a face
unsigned int face_num_neighbours(const MFace& face) const
{
FaceNRangeConstIterator itFace = faceNRange.find(face);
if(itFace == faceNRange.end()) return 0; // Face not found
return itFace->second.num;
}
//--Clear the containers
void clear()
{
vertNRange.clear();
edgeNRange.clear();
faceNRange.clear();
neighbours.clear();
maxVertNeighbours = 0;
maxEdgeNeighbours = 0;
maxFaceNeighbours = 0;
}
//--Return the element neighbour across a d-1 polytope
MElement *element_neighbour(MElement *const element,
const unsigned int iPolytope) const
{
unsigned int num;
NeighboursConstIterator itElem;
switch(element->getTypeForMSH()) {
case MSH_LIN_2:
case MSH_LIN_3:
{
MVertex *vertex = element->getVertex(iPolytope);
VertNRangeConstIterator itVert = vertNRange.find(vertex);
if(itVert == vertNRange.end()) return 0;
num = itVert->second.num;
itElem = itVert->second.begin;
}
break;
case MSH_TRI_3:
case MSH_TRI_6:
case MSH_QUA_4:
case MSH_QUA_8:
case MSH_QUA_9:
{
MEdge edge = element->getEdge(iPolytope);
EdgeNRangeConstIterator itEdge = edgeNRange.find(edge);
if(itEdge == edgeNRange.end()) return 0;
num = itEdge->second.num;
itElem = itEdge->second.begin;
}
break;
case MSH_TET_4:
case MSH_TET_10:
case MSH_HEX_8:
case MSH_HEX_20:
case MSH_HEX_27:
case MSH_PRI_6:
case MSH_PRI_15:
case MSH_PRI_18:
case MSH_PYR_5:
case MSH_PYR_13:
case MSH_PYR_14: {
MFace face = element->getFace(iPolytope);
FaceNRangeConstIterator itFace = faceNRange.find(face);
if(itFace == faceNRange.end()) return 0;
num = itFace->second.num;
itElem = itFace->second.begin;
}
break;
}
if(*itElem != element) return *itElem;
if(num > 1) {
++itElem;
if(*itElem != element) return *itElem;
}
return 0;
}
private:
//--Data members
VertNRange vertNRange;
EdgeNRange edgeNRange;
FaceNRange faceNRange;
Neighbours neighbours;
unsigned int maxVertNeighbours;
unsigned int maxEdgeNeighbours;
unsigned int maxFaceNeighbours;
//--Working routine to determine the neighbours from elements
template <typename Elem, typename ElemIter>
void add_elements1(ElemIter begin, ElemIter end)
{
// Test if the type of element is known
if(ElemTr<Elem>::numVertex) { // 0 only for MElement
// Find the neighbours of each vertex, edge, and face
for(ElemIter itElem = begin; itElem != end; ++itElem) {
// Find neighbours for vertices
FindNeighbours<Elem, MVertex*, ElemTr<Elem>::numVertex>::
eval(*itElem, vertNRange, neighbours);
// Find neighbours for edges
FindNeighbours<Elem, MEdge, ElemTr<Elem>::numEdge>::
eval(*itElem, edgeNRange, neighbours);
// Find neighbours for faces
FindNeighbours<Elem, MFace, ElemTr<Elem>::numFace>::
eval(*itElem, faceNRange, neighbours);
}
}
else {
// We only have an iterator to an MElement* and have to check the type in
// run-time.
for(ElemIter itElem = begin; itElem != end; ++itElem) {
switch((*itElem)->getTypeForMSH()) {
case MSH_LIN_2:
case MSH_LIN_3:
find_polytope_neighbours<MLine>(*itElem);
break;
case MSH_TRI_3:
case MSH_TRI_6:
find_polytope_neighbours<MTriangle>(*itElem);
break;
case MSH_QUA_4:
case MSH_QUA_8:
case MSH_QUA_9:
find_polytope_neighbours<MQuadrangle>(*itElem);
break;
case MSH_TET_4:
case MSH_TET_10:
find_polytope_neighbours<MTetrahedron>(*itElem);
break;
case MSH_HEX_8:
case MSH_HEX_20:
case MSH_HEX_27:
find_polytope_neighbours<MHexahedron>(*itElem);
break;
case MSH_PRI_6:
case MSH_PRI_15:
case MSH_PRI_18:
find_polytope_neighbours<MPrism>(*itElem);
break;
case MSH_PYR_5:
case MSH_PYR_13:
case MSH_PYR_14:
find_polytope_neighbours<MPyramid>(*itElem);
break;
}
}
}
maxNeighbours();
}
//--Find the neighbours sharing each polytope
template <typename Elem>
void find_polytope_neighbours(MElement *element)
{
// Find neighbours for vertices
FindNeighbours<Elem, MVertex*, ElemTr<Elem>::numVertex>::
eval(static_cast<Elem*>(element), vertNRange, neighbours);
// Find neighbours for edges
FindNeighbours<Elem, MEdge, ElemTr<Elem>::numEdge>::
eval(static_cast<Elem*>(element), edgeNRange, neighbours);
// Find neighbours for faces
FindNeighbours<Elem, MFace, ElemTr<Elem>::numFace>::
eval(static_cast<Elem*>(element), faceNRange, neighbours);
}
//--Compute max neighbours
void maxNeighbours()
{
maxVertNeighbours = 0;
for(VertNRangeConstIterator itVert = vertNRange.begin();
itVert != vertNRange.end(); ++itVert)
maxVertNeighbours = std::max(maxVertNeighbours, itVert->second.num);
maxEdgeNeighbours = 0;
for(EdgeNRangeConstIterator itEdge = edgeNRange.begin();
itEdge != edgeNRange.end(); ++itEdge)
maxEdgeNeighbours = std::max(maxEdgeNeighbours, itEdge->second.num);
maxFaceNeighbours = 0;
for(FaceNRangeConstIterator itFace = faceNRange.begin();
itFace != faceNRange.end(); ++itFace)
maxFaceNeighbours = std::max(maxFaceNeighbours, itFace->second.num);
}
};
/*==============================================================================
* Template meta-programming classes
*============================================================================*/
/*--------------------------------------------------------------------
* This is the main working routine for finding the elements sharing
* a lower-dimension polytope (i.e, the structures bounding the
* element: vertex, edge, or face). It is templated for any type of
* polytope and uses template meta-programming to unroll the loop over
* the number of that type of polytope in the element.
*--------------------------------------------------------------------*/
//--Entry point
template <typename Elem, // Elem is the type of element
typename Polytope, // Polytope is a lower dimensional
// bounding object (MVertex*, MEdge, or
// MFace) of the element. We want to
// find the elements sharing a given
// Polytope.
unsigned int NP> // NP is an index through the number
// of Polytope in the element (e.g.,
// number of vertices)
struct MNeighbour::FindNeighbours
{
typedef typename PolytopeTr<Polytope>::PolytopeNRange PolytopeNRange;
static void eval(Elem *element, PolytopeNRange &polytopeNRange,
Neighbours &neighbours)
{
Polytope polytope = PolytopeTr<Polytope>::template
getPolytope<Elem>(element, NP-1);
Range_t &range = polytopeNRange[polytope];
if(range.num == 0) range.begin = neighbours.end();
range.begin = neighbours.insert(range.begin, element);
++range.num;
FindNeighbours<Elem, Polytope, NP-1>::eval(element, polytopeNRange,
neighbours);
}
};
//--Terminate loop when NP == 1
template <typename Elem, typename Polytope>
struct MNeighbour::FindNeighbours<Elem, Polytope, 1>
{
typedef typename PolytopeTr<Polytope>::PolytopeNRange PolytopeNRange;
static void eval(Elem *element, PolytopeNRange &polytopeNRange,
Neighbours &neighbours)
{
Polytope polytope = PolytopeTr<Polytope>::template
getPolytope<Elem>(element, 0);
Range_t &range = polytopeNRange[polytope];
if(range.num == 0) range.begin = neighbours.end();
range.begin = neighbours.insert(range.begin, element);
++range.num;
}
};
//--Nothing to do if the dimension of this polytope is equal to or greater than
//--that of the element (e.g., no faces (2D) in a 2D element). This is
//--indicated by NP = 0.
template <typename Elem, typename Polytope>
struct MNeighbour::FindNeighbours<Elem, Polytope, 0>
{
typedef typename PolytopeTr<Polytope>::PolytopeNRange PolytopeNRange;
static void eval(Elem *element, PolytopeNRange &polytopeNRange,
Neighbours &neighbours) { }
};
/*--------------------------------------------------------------------*
* This class uses traits to determine the element types in an entity
* and iterate over them. A template meta-programming loop is used to
* iterate over the types of elements.
*--------------------------------------------------------------------*/
//--Entry point
template <typename Ent, unsigned int N>
struct MNeighbour::FindNeighboursInEntity {
typedef typename EntElemTr<Ent, N>::Elem Elem;
static void eval(const Ent *const entity, VertNRange &vertNRange,
EdgeNRange &edgeNRange, FaceNRange &faceNRange,
Neighbours &neighbours)
{
for(typename EntElemTr<Ent, N>::ConstElementIterator itElem =
EntElemTr<Ent, N>::begin(entity);
itElem != EntElemTr<Ent, N>::end(entity); ++itElem) {
// Find neighbours for vertices
FindNeighbours<Elem, MVertex*, ElemTr<Elem>::numVertex>::
eval(*itElem, vertNRange, neighbours);
// Find neighbours for edges
FindNeighbours<Elem, MEdge, ElemTr<Elem>::numEdge>::
eval(*itElem, edgeNRange, neighbours);
// Find neighbours for faces
FindNeighbours<Elem, MFace, ElemTr<Elem>::numFace>::
eval(*itElem, faceNRange, neighbours);
}
FindNeighboursInEntity<Ent, N-1>::eval(entity, vertNRange, edgeNRange,
faceNRange, neighbours);
}
};
//--Terminate loop when N = 1
template <typename Ent>
struct MNeighbour::FindNeighboursInEntity<Ent, 1> {
typedef typename EntElemTr<Ent, 1>::Elem Elem;
static void eval(const Ent *const entity, VertNRange &vertNRange,
EdgeNRange &edgeNRange, FaceNRange &faceNRange,
Neighbours &neighbours)
{
for(typename EntElemTr<Ent, 1>::ConstElementIterator itElem =
EntElemTr<Ent, 1>::begin(entity);
itElem != EntElemTr<Ent, 1>::end(entity); ++itElem) {
// Find neighbours for vertices
FindNeighbours<Elem, MVertex*, ElemTr<Elem>::numVertex>::
eval(*itElem, vertNRange, neighbours);
// Find neighbours for edges
FindNeighbours<Elem, MEdge, ElemTr<Elem>::numEdge>::
eval(*itElem, edgeNRange, neighbours);
// Find neighbours for faces
FindNeighbours<Elem, MFace, ElemTr<Elem>::numFace>::
eval(*itElem, faceNRange, neighbours);
}
}
};
#endif
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment