// ============================================================================
// == ==
// == AMDiS - Adaptive multidimensional simulations ==
// == ==
// == http://www.amdis-fem.org ==
// == ==
// ============================================================================
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.
/** \file Element.h */
#ifndef AMDIS_ELEMENT_H
#define AMDIS_ELEMENT_H
#include "AMDiS_fwd.h"
#include "Global.h"
#include "RefinementManager.h"
#include "Serializable.h"
#include "ElementData.h"
#include "LeafData.h"
namespace AMDiS {
template<typename T, GeoIndex d> class FixVec;
#define AMDIS_UNDEFINED 5
/** \ingroup Triangulation
* \brief
* Base class for Line, Triangle, Tetrahedron
*
* Elements in AMDiS are always simplices (a simplex is a Line in 1d, a
* Triangle in 2d and a Tetrahedron in 3d).
* We restrict ourselves here to simplicial meshes, for several reasons:
* -# A simplex is one of the most simple geometric types and complex domains
* may be approximated by a set of simplices quite easily.
* -# Simplicial meshes allow local refinement without the need of
* nonconforming meshes (hanging nodes), parametric elements, or mixture of
* element types (which is the case for quadrilateral meshes).
* -# Polynomials of any degree are easily represented on a simplex using
* local (barycentric) coordinates.
*
* A Line element and its refinement:
*
* <img src="line.png">
*
* A Triangle element and its refinement:
*
* <img src="triangle.png">
*
* A Tetrahedron element and its refinements:
*
* <img src="tetrahedron.png">
*/
class Element : public Serializable
{
private:
/// private standard constructor because an Element must know his Mesh
Element() {}
public:
/// constructs an Element which belongs to Mesh
Element(Mesh *);
/// copy constructor
Element(const Element& old);
/// destructor
virtual ~Element();
///
void deleteElementDOFs();
/** \brief
* Clone this Element and return a reference to it. Because also the DOFs
* are cloned, \ref Mesh::serializedDOfs must be used.
*/
Element* cloneWithDOFs();
/** \name getting methods
* \{
*/
/// Returns \ref child[0]
inline Element* getFirstChild() const
{
return child[0];
}
/// Returns \ref child[1]
inline Element* getSecondChild() const
{
return child[1];
}
/// Returns \ref child[i], i=0,1
inline Element* getChild(int i) const
{
TEST_EXIT_DBG(i == 0 || i == 1)("There is only child 0 or 1! (i = %d)\n", i);
return child[i];
}
/** \brief
* Returns true if Element is a leaf element (\ref child[0] == NULL), returns
* false otherwise.
*/
inline const bool isLeaf() const
{
return (child[0] == NULL);
}
/// Returns \ref dof[i][j] which is the j-th DOF of the i-th node of Element.
inline DegreeOfFreedom getDof(int i, int j) const
{
TEST_EXIT_DBG(dofValid)("DOFs are not valid in element %d!\n", index);
return dof[i][j];
}
/// Returns \ref dof[i] which is a pointer to the DOFs of the i-th node.
inline const DegreeOfFreedom* getDof(int i) const
{
TEST_EXIT_DBG(dofValid)("DOFs are not valid in element %d!\n", index);
return dof[i];
}
/// Returns a pointer to the DOFs of this Element
inline const DegreeOfFreedom** getDof() const
{
TEST_EXIT_DBG(dofValid)("DOFs are not valid in element %d!\n", index);
return const_cast<const DegreeOfFreedom**>(dof);
}
/// Returns \ref mesh of Element
inline Mesh* getMesh() const
{
return mesh;
}
inline void setDofValid(bool b)
{
dofValid = b;
}
/** \brief
* Returns \ref elementData's error estimation, if Element is a leaf element
* and has leaf data.
*/
inline double getEstimation(int row) const
{
if (isLeaf()) {
TEST_EXIT_DBG(elementData)("leaf element without leaf data\n");
ElementData *ld = elementData->getElementData(ESTIMATABLE);
TEST_EXIT_DBG(ld)("leaf data not estimatable!\n");
return dynamic_cast<LeafDataEstimatableInterface*>(ld)->getErrorEstimate(row);
}
return 0.0;
}
/** \brief
* Returns Element's coarsening error estimation, if Element is a leaf
* element and if it has leaf data and if this leaf data are coarsenable.
*/
inline double getCoarseningEstimation(int row)
{
if (isLeaf()) {
TEST_EXIT_DBG(elementData)("leaf element without leaf data\n");
ElementData *ld = elementData->getElementData(COARSENABLE);
TEST_EXIT_DBG(ld)("element data not coarsenable!\n");
return dynamic_cast<LeafDataCoarsenableInterface*>(ld)->getCoarseningErrorEstimate(row);
}
return 0.0;
}
/// Returns region of element if defined, -1 else.
int getRegion() const;
/// Returns local vertex number of the j-th vertex of the i-th edge
virtual int getVertexOfEdge(int i, int j) const = 0;
/** \brief
* Returns local vertex number of the vertexIndex-th vertex of the
* positionIndex-th part of type position (vertex, edge, face)
*/
virtual int getVertexOfPosition(GeoIndex position,
int positionIndex,
int vertexIndex) const = 0;
///
virtual int getPositionOfVertex(int side, int vertex) const = 0;
///
virtual int getEdgeOfFace(int face, int edge) const = 0;
///
virtual DofEdge getEdge(int localEdgeIndex) const = 0;
///
virtual DofFace getFace(int localFaceIndex) const = 0;
/// Returns either vertex, edge or face of the element.
/* template<typename T> */
/* T getObject(int index); */
/// Returns the number of parts of type i in this element
virtual int getGeo(GeoIndex i) const = 0;
/// Returns Element's \ref mark
inline const int getMark() const
{
return mark;
}
/// Returns \ref newCoord[i]
inline double getNewCoord(int i) const
{
TEST_EXIT_DBG(newCoord)("newCoord = NULL\n");
return (*newCoord)[i];
}
/// Returns Element's \ref index
inline int getIndex() const
{
return index;
}
/// Returns \ref newCoord
inline WorldVector<double>* getNewCoord() const
{
return newCoord;
}
/** \} */
/** \name setting methods
* \{
*/
/// Sets \ref child[0]
virtual void setFirstChild(Element *aChild)
{
child[0] = aChild;
}
/// Sets \ref child[1]
virtual void setSecondChild(Element *aChild)
{
child[1] = aChild;
}
/// Sets \ref elementData of Element
void setElementData(ElementData* ed)
{
elementData = ed;
}
/** \brief
* Sets \ref newCoord of Element. Needed by refinement, if Element has a
* boundary edge on a curved boundary.
*/
inline void setNewCoord(WorldVector<double>* coord)
{
newCoord = coord;
}
/// Sets \ref mesh.
inline void setMesh(Mesh *m)
{
mesh = m;
}
/// Sets the pointer to the DOFs of the i-th node of Element
inline void setDof(int pos, DegreeOfFreedom* p)
{
dof[pos] = p;
}
/** \brief
* Checks whether Element is a leaf element and whether it has leaf data.
* If the checks don't fail, leaf data's error estimation is set to est.
*/
inline void setEstimation(double est, int row)
{
FUNCNAME("Element::setEstimation()");
if (isLeaf()) {
TEST_EXIT_DBG(elementData)("Leaf element %d without leaf data!\n", index);
ElementData *ld = elementData->getElementData(ESTIMATABLE);
TEST_EXIT_DBG(ld)("Leaf data %d not estimatable!\n", index);
dynamic_cast<LeafDataEstimatableInterface*>(ld)->
setErrorEstimate(row, est);
} else {
ERROR_EXIT("setEstimation only for leaf elements!\n");
}
}
/** \brief
* Sets Element's coarsening error estimation, if Element is a leaf element
* and if it has leaf data and if this leaf data are coarsenable.
*/
inline void setCoarseningEstimation(double est, int row)
{
if (isLeaf()) {
TEST_EXIT_DBG(elementData)("leaf element without leaf data\n");
ElementData *ld = elementData->getElementData(COARSENABLE);
TEST_EXIT_DBG(ld)("leaf data not coarsenable\n");
dynamic_cast<LeafDataCoarsenableInterface*>(ld)->
setCoarseningErrorEstimate(row, est);
} else {
ERROR_EXIT("setEstimation only for leaf elements!\n");
}
}
/// Sets Elements \ref mark = mark + 1;
inline void incrementMark()
{
mark++;
}
/// Sets Elements \ref mark = mark - 1;
inline void decrementMark()
{
if (0 < mark)
mark--;
}
/// Sets Element's \ref mark
inline void setMark(int m)
{
mark = m;
}
/** \} */
/** \name pure virtual methods
* \{
*/
/** \brief
* orient the vertices of edges/faces.
* Used by Estimator for the jumps => same quadrature nodes from both sides!
*/
virtual const FixVec<int,WORLD>&
sortFaceIndices(int face, FixVec<int, WORLD> *vec) const = 0;
/// Returns a copy of itself. Needed by Mesh to create Elements by a prototype.
virtual Element *clone() = 0;
/** \brief
* Returns which side of child[childnr] corresponds to side sidenr of this
* Element. If the child has no corresponding side, the return value is negative.
*/
virtual int getSideOfChild(int childnr, int sidenr, int elType = 0) const = 0;
/** \brief
* Generalization of \ref getSideOfChild to arbitrary subObject. Thus, e.g., in 3d
* we can ask for the local id of a verte, edge or face on the elements children.
*
* \param[in] childnr Either 0 or 1 for the left or right children.
* \param[in] subObj Defines whether we ask for VERTEX, EDGE or FACE.
* \param[in] ithObj Number of the object on the parent.
* \param[in] elType Type of the element. Important only in 3D.
*/
virtual int getSubObjOfChild(int childnr, GeoIndex subObj, int ithObj,
int elType = 0) const = 0;
/** \brief
* Returns which vertex of elements parent corresponds to the vertexnr of
* the element, if the element is the childnr-th child of the parent.
* If the vertex is the ner vertex at the refinement edge, -1 is returned.
*/
virtual int getVertexOfParent(int childnr, int vertexnr, int elType = 0) const = 0;
/// Returns whether Element is a Line
virtual bool isLine() const = 0;
/// Returns whether Element is a Triangle
virtual bool isTriangle() const = 0;
/// Returns whether Element is a Tetrahedron
virtual bool isTetrahedron() const = 0;
/// Returns whether Element has sideElem as one of its sides.
virtual bool hasSide(Element *sideElem) const = 0;
/** \brief
* Returns for a given element type number the element type number of the children.
* For 1d and 2d this is always 0, because element type number are used in the
* 3d case only.
*/
virtual int getChildType(int elType) const = 0;
/** \brief
* Traverses an edge/face of a given element (this includes also all children of
* the element having the same edge/face). All vertex dofs alonge this edge/face
* are assembled and put together to a list.
*
* \param[in] feSpace FE space which is used to get the dofs.
* \param[in] bound Defines the edge/face of the element on which
* all vertex dofs are assembled.
* \param[out] dofs List of dofs, where the result is stored.
*/
virtual void getVertexDofs(FiniteElemSpace* feSpace,
BoundaryObject bound,
DofContainer& dofs) const = 0;
/** \brief
* Traverses an edge/face of a given element (this includes also all children of
* the element having the same edge/face). All non vertex dofs alonge this edge/face
* are assembled and put together to a list.
*
* \param[in] feSpace FE space which is used to get the dofs.
* \param[in] bound Defines the edge/face of the element on which
* all non vertex dofs are assembled.
* \param[out] dofs All dofs are put to this dof list.
*/
virtual void getNonVertexDofs(FiniteElemSpace* feSpace,
BoundaryObject bound,
DofContainer& dofs) const = 0;
/** \} */
// ===== other public methods =================================================
/// assignment operator
Element& operator=(const Element& el);
/** \brief
* Checks whether the face with vertices dof[0],..,dof[DIM-1] is
* part of mel's boundary. returns the opposite vertex if true, -1 else
*/
int oppVertex(FixVec<DegreeOfFreedom*, DIMEN> pdof) const;
/// Refines Element's leaf data
inline void refineElementData(Element* child1, Element* child2, int elType = 0)
{
if (elementData) {
bool remove = elementData->refineElementData(this, child1, child2, elType);
if (remove) {
ElementData *tmp = elementData->getDecorated();
delete elementData;
elementData = tmp;
}
}
}
/// Coarsens Element's leaf data
inline void coarsenElementData(Element* child1, Element* child2, int elType = 0)
{
ElementData *childData;
childData = child1->getElementData();
if (childData) {
childData->coarsenElementData(this, child1, child2, elType);
delete childData;
child1->setElementData(NULL);
}
childData = child2->getElementData();
if (childData) {
childData->coarsenElementData(this, child2, child1, elType);
delete childData;
child2->setElementData(NULL);
}
}
/// Returns pointer to \ref elementData
inline ElementData* getElementData() const
{
return elementData;
}
///
inline ElementData* getElementData(int typeID) const
{
if (elementData)
return elementData->getElementData(typeID);
return NULL;
}
/// Deletes the \ref elementData with a specific typeID.
bool deleteElementData(int typeID);
/** \brief
* Returns whether element is refined at side side
* el1, el2 are the corresponding children.
* (not neccessarly the direct children!)
* elementTyp is the type of this element (comes from ElInfo)
*/
bool isRefinedAtSide(int side, Element *el1, Element *el2,
unsigned char elementTyp = 255);
/// Returns whether Element's \ref newCoord is set
inline bool isNewCoordSet() const
{
return (newCoord != NULL);
}
/// Frees memory for \ref newCoord
void eraseNewCoord();
/// Serialize the element to a file.
void serialize(std::ostream &out);
/// Deserialize an element from a file.
void deserialize(std::istream &in);
/// Sets Element's \ref dof pointer.
void createNewDofPtrs();
protected:
/// Sets Element's \ref index. Used by friend class Mesh.
inline void setIndex(int i)
{
index = i;
}
/// Used by friend class Mesh while dofCompress
void newDofFct1(const DOFAdmin*, std::vector<int>&);
/// Used by friend class Mesh while dofCompress
void newDofFct2(const DOFAdmin*);
/// Changes old dofs to negative new dofs
void changeDofs1(const DOFAdmin* admin, std::vector<int>& newDofIndex,
int n0, int nd0, int nd, int pos);
/// Changes negative new dofs to positive
void changeDofs2(int n0, int nd0, int nd, int pos);
protected:
/** \brief
* Pointers to the two children of interior elements of the tree. Pointers
* to NULL for leaf elements.
*/
Element *child[2];
/** \brief
* Vector of pointers to DOFs. These pointers must be available for elements
* vertices (for the geometric description of the mesh). There my be pointers
* for the edges, for faces and for the center of an element. They are
* ordered the following way: The first N_VERTICES entries correspond to the
* DOFs at the vertices of the element. The next ones are those at the edges,
* if present, then those at the faces, if present, and then those at the
* barycenter, if present.
*/
DegreeOfFreedom **dof;
/** \brief
* Unique global index of the element. these indices are not strictly ordered
* and may be larger than the number of elements in the binary tree (the list
* of indices may have holes after coarsening).
*/
int index;
/** \brief
* Marker for refinement and coarsening. if mark is positive for a leaf
* element, this element is refined mark times. if mark is negative for
* a leaf element, this element is coarsened -mark times.
*/
int mark;
/** \brief
* If the element has a boundary edge on a curved boundary, this is a pointer
* to the coordinates of the new vertex that is created due to the refinement
* of the element, otherwise it is a NULL pointer. Thus coordinate
* information can be also produced by the traversal routines in the case of
* curved boundary.
*/
WorldVector<double> *newCoord;
/// Pointer to the Mesh this element belongs to
Mesh* mesh;
/// Pointer to Element's leaf data
ElementData* elementData;
/** \brief
* Is true, if the DOF pointers are valid. In sequential computations this
* is always the case. In parallel computations all macro elements are stored
* in memory though not all of them are part of the mesh (because they are owned
* by another rank). In this case, there are no valid DOFs on the element.
*/
bool dofValid;
/** \brief
* This map is used for deletion of all DOFs of all elements of a mesh. Once
* a DOF-vector (all DOFS at a node, edge, etc.) is deleted, its address is
* added to this map to note not to delete it a second time.
*/
static std::map<DegreeOfFreedom*, bool> deletedDOFs;
friend class Mesh;
};
/// Writes the element hierarchie to a Graphviz dot-file. Using the dot-tool from
/// Graphvis, this dot-file can be converted to a ps-file. Useful for debugging!
void writeDotFile(Element *el, std::string filename, int maxLevels = -1);
}
#include "Element.hh"
#endif // AMDIS_ELEMENT_H