Commit 6aa136b2 authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

* On the way to independent meshes

parent f0f8d33c
......@@ -55,13 +55,7 @@ namespace AMDiS {
// call standard traverse
*elInfo1 = stack1.traverseFirst(mesh1, level1, flag1);
while (*elInfo1 != NULL && skipEl1(*elInfo1)) {
*elInfo1 = stack1.traverseNext(*elInfo1);
}
*elInfo2 = stack2.traverseFirst(mesh2, level2, flag2);
while (*elInfo2 != NULL && skipEl2(*elInfo2)) {
*elInfo2 = stack2.traverseNext(*elInfo2);
}
// finished ?
if (*elInfo1 == NULL || *elInfo2 == NULL) {
......@@ -71,8 +65,7 @@ namespace AMDiS {
rest = 1.0;
bool accepted = check(elInfo1, elInfo2,
elInfoSmall, elInfoLarge);
bool accepted = check(elInfo1, elInfo2, elInfoSmall, elInfoLarge);
// check for non domain covering level traverse
if (!accepted ||
......@@ -95,14 +88,10 @@ namespace AMDiS {
{
// call standard traverse
if (inc1) {
do {
*elInfo1 = stack1.traverseNext(*elInfo1);
} while (*elInfo1 != NULL && skipEl1(*elInfo1));
}
if (inc2) {
do {
*elInfo2 = stack2.traverseNext(*elInfo2);
} while (*elInfo2 != NULL && skipEl2(*elInfo2));
}
// finished ?
......@@ -111,8 +100,7 @@ namespace AMDiS {
return false;
}
bool accepted = check(elInfo1, elInfo2,
elInfoSmall, elInfoLarge);
bool accepted = check(elInfo1, elInfo2, elInfoSmall, elInfoLarge);
// check for non domain covering level traverse
if (!accepted ||
......@@ -144,17 +132,34 @@ namespace AMDiS {
*elInfo2;
// update rest
rest -= 1.0 / (1 << ((*elInfoSmall)->getLevel() - (*elInfoLarge)->getLevel()));
double newRest = rest - 1.0 / (1 << ((*elInfoSmall)->getLevel() - (*elInfoLarge)->getLevel()));
if (rest < 1e-32) {
if (newRest < 1e-32) {
// large element finished -> increment both elements
rest = 1.0;
inc1 = true;
inc2 = true;
stack1.stopDisplacementCalculation();
stack2.stopDisplacementCalculation();
} else {
// increment only small element
inc1 = (*elInfo1 == *elInfoSmall) ? true : false;
inc2 = (*elInfo2 == *elInfoSmall) ? true : false;
// If rest is 1, than this is the first element pair with this large
// element. Therefore, fill the displacement stack of the small element
// traverse stack.
if (rest == 1.0) {
if (*elInfo1 == *elInfoSmall) {
stack1.startDisplacementCalculation((*elInfo2)->getLevel());
// (*elInfo1)->setDisplacement(0);
} else {
stack2.startDisplacementCalculation((*elInfo1)->getLevel());
// (*elInfo2)->setDisplacement(0);
}
}
rest = newRest;
}
}
......@@ -163,12 +168,17 @@ namespace AMDiS {
ElInfo *elInfoSmall,
ElInfo *elInfoLarge)
{
if (!fillSubElemInfo)
if (!fillSubElemMat)
return;
VectorOfFixVecs<DimVec<double> > *subCoords = elInfoSmall->getSubElemCoords();
if (!subCoords) {
subCoords = NEW VectorOfFixVecs<DimVec<double> >(1, 2, NO_INIT);
VectorOfFixVecs<DimVec<double> > *subCoords =
NEW VectorOfFixVecs<DimVec<double> >(elInfo1->getMesh()->getDim(),
elInfo1->getMesh()->getDim() + 1,
NO_INIT);
DimMat<double> *subCoordsMat = elInfoSmall->getSubElemCoordsMat();
if (!subCoordsMat) {
subCoordsMat = NEW DimMat<double>(elInfo1->getMesh()->getDim());
}
if (elInfo1 == elInfoSmall) {
......@@ -176,6 +186,15 @@ namespace AMDiS {
} else {
stack2.getCoordsInElem(elInfo1, subCoords);
}
elInfoSmall->setSubElemCoords(subCoords);
for (int i = 0; i < elInfo1->getMesh()->getDim() + 1; i++) {
for (int j = 0; j < elInfo1->getMesh()->getDim() + 1; j++) {
(*subCoordsMat)[j][i] = (*subCoords)[i][j];
}
}
elInfoSmall->setSubElemCoordsMat(subCoordsMat);
DELETE subCoords;
}
}
......@@ -44,10 +44,11 @@ namespace AMDiS {
MEMORY_MANAGED(DualTraverse);
DualTraverse()
: fillSubElemInfo(false)
{};
: fillSubElemMat(false)
{}
~DualTraverse() {};
~DualTraverse()
{}
/** \brief
* Start dual traversal
......@@ -71,14 +72,6 @@ namespace AMDiS {
ElInfo **elInfoSmall,
ElInfo **elInfoLarge);
bool skipEl1(ElInfo *elInfo) {
return false;
}
bool skipEl2(ElInfo *elInfo) {
return false;
}
bool check(ElInfo **elInfo1,
ElInfo **elInfo2,
ElInfo **elInfoSmall,
......@@ -88,8 +81,8 @@ namespace AMDiS {
return true;
}
inline void setFillSubElemInfo(bool b) {
fillSubElemInfo = b;
inline void setFillSubElemMat(bool b) {
fillSubElemMat = b;
}
protected:
......@@ -154,7 +147,12 @@ namespace AMDiS {
*/
bool callLeafElLevel2_;
bool fillSubElemInfo;
/** \brief
* If true, during dual mesh traverse for the smaller element the transformation
* matrix will be computed. This matrix defines the transformation mapping for
* points defined on the larger element to the coordinates of the smaller element.
*/
bool fillSubElemMat;
};
}
......
......@@ -23,7 +23,8 @@ namespace AMDiS {
element_(NULL),
parent_(NULL),
macroElement_(NULL),
level_(0),
level(0),
displacement(0),
iChild(0),
coord_(mesh_->getDim(), NO_INIT),
boundary_(mesh_->getDim(), DEFAULT_VALUE, INTERIOR),
......@@ -33,7 +34,7 @@ namespace AMDiS {
neighbourCoord_(mesh_->getDim(), NO_INIT),
oppVertex_(mesh_->getDim(), NO_INIT),
grdLambda_(mesh_->getDim(), NO_INIT),
subElemCoords(NULL)
subElemCoordsMat(NULL)
{
projection_.set(NULL);
......@@ -46,8 +47,8 @@ namespace AMDiS {
ElInfo::~ElInfo()
{
if (subElemCoords) {
DELETE subElemCoords;
if (subElemCoordsMat) {
DELETE subElemCoordsMat;
}
}
......
......@@ -90,7 +90,7 @@ namespace AMDiS {
parent_ = rhs.parent_;
macroElement_ = rhs.macroElement_;
fillFlag_ = rhs.fillFlag_;
level_ = rhs.level_;
level = rhs.level;
iChild = rhs.iChild;
coord_ = rhs.coord_;
boundary_ = rhs.boundary_;
......@@ -99,7 +99,7 @@ namespace AMDiS {
neighbourCoord_ = rhs.neighbourCoord_;
oppVertex_ = rhs.oppVertex_;
return *this;
};
}
// ===== getting-methods ======================================================
public:
......@@ -112,49 +112,56 @@ namespace AMDiS {
*/
inline Mesh* getMesh() const {
return mesh_;
};
}
/** \brief
* Get ElInfo's \ref macroElement_
*/
inline MacroElement* getMacroElement() const {
return macroElement_;
};
}
/** \brief
* Get ElInfo's \ref element
*/
inline Element* getElement() const {
return element_;
};
}
/** \brief
* Get ElInfo's \ref parent_
*/
inline Element* getParent() const {
return parent_;
};
}
/** \brief
* Get ElInfo's \ref fillFlag_
*/
inline Flag getFillFlag() const {
return fillFlag_;
};
}
/** \brief
* Get ElInfo's \ref level_
* Get ElInfo's \ref level
*/
inline unsigned char getLevel() const {
return level_;
};
inline int getLevel() const {
return level;
}
/** \brief
* Get ElInfo's \ref displacement. Is used only during dual traverse.
*/
inline int getDisplacement() const {
return displacement;
}
/** \brief
* Get ElInfo's \ref iChild
*/
inline int getIChild() const {
return iChild;
};
}
/** \brief
* Get ElInfo's \ref coord_[i]. This is a WorldVector<double> filled with the
......@@ -162,7 +169,7 @@ namespace AMDiS {
*/
inline WorldVector<double>& getCoord(int i) {
return coord_[i];
};
}
/** \brief
* Get ElInfo's \ref coord_[i]. This is a WorldVector<double> filled with the
......@@ -170,7 +177,7 @@ namespace AMDiS {
*/
inline const WorldVector<double>& getCoord(int i) const {
return coord_[i];
};
}
/** \brief
* Get ElInfo's \ref coord_. This is a FixVec<WorldVector<double> > filled with the
......@@ -178,7 +185,7 @@ namespace AMDiS {
*/
inline FixVec<WorldVector<double>, VERTEX>& getCoords() {
return coord_;
};
}
/** \brief
* Get ElInfo's \ref coord_. This is a FixVec<WorldVector<double> > filled with the
......@@ -186,7 +193,7 @@ namespace AMDiS {
*/
inline const FixVec<WorldVector<double>, VERTEX>& getCoords() const {
return coord_;
};
}
/** \brief
* Get ElInfo's \ref oppCoord_[i]
......@@ -200,7 +207,7 @@ namespace AMDiS {
*/
inline BoundaryType getBoundary(int i) const {
return boundary_[i];
};
}
/** \brief
* Get boundary type of i-th vertex/edge/face (pos).
......@@ -260,8 +267,12 @@ namespace AMDiS {
return parametric_;
}
inline VectorOfFixVecs<DimVec<double> > *getSubElemCoords() const {
return subElemCoords;
/** \brief
* Returns the element transformation matrix \ref subElemCoordsMat .
* This value is set only during dual traverse.
*/
inline DimMat<double> *getSubElemCoordsMat() const {
return subElemCoordsMat;
}
/** \} */
......@@ -277,56 +288,63 @@ namespace AMDiS {
*/
inline void setMesh(Mesh* aMesh) {
mesh_ = aMesh;
};
}
/** \brief
* Set ElInfo's \ref macroElement_
*/
inline void setMacroElement(MacroElement* mel) {
macroElement_ = mel;
};
}
/** \brief
* Set ElInfo's \ref element
*/
inline void setElement(Element* elem) {
element_ = elem;
};
}
/** \brief
* Set ElInfo's \ref parent_
*/
inline void setParent(Element* elem) {
parent_ = elem;
};
}
/** \brief
* Set ElInfo's \ref fillFlag_
*/
inline void setFillFlag(Flag flag) {
fillFlag_ = flag;
};
}
/** \brief
* Sets ElInfo's \ref coord_[i].
*/
inline void setCoord(int i,WorldVector<double>& coord) {
coord_[i] = coord;
};
}
/** \brief
* Sets ElInfo's \ref coord.
*/
inline void setCoords(FixVec<WorldVector<double>,VERTEX >& coords) {
coord_ = coords;
};
}
/** \brief
* Set ElInfo's \ref level_
* Set ElInfo's \ref level
*/
inline void setLevel(int l) {
level_ = l;
};
level = l;
}
/** \brief
* Set ElInfo's \ref displacement. Used only during dual traverse.
*/
inline void setDisplacement(int d) {
displacement = d;
}
/** \brief
* Set ElInfo's \ref boundary_[i]
......@@ -356,8 +374,12 @@ namespace AMDiS {
parametric_ = param;
}
inline void setSubElemCoords(VectorOfFixVecs<DimVec<double> > *coords) {
subElemCoords = coords;
/** \brief
* Sets the element transformation matrix \ref subElemCoordsMat .
* This value is used only during dual traverse.
*/
inline void setSubElemCoordsMat(DimMat<double> *mat) {
subElemCoordsMat = mat;
}
/** \} */
......@@ -390,8 +412,7 @@ namespace AMDiS {
* Otherwise the function itself provides memory for this vector. In this
* case the vector is overwritten during the next call of coordToWorld.
*/
virtual const WorldVector<double>
*coordToWorld(const DimVec<double>& lambda,
virtual const WorldVector<double> *coordToWorld(const DimVec<double>& lambda,
WorldVector<double>* world) const;
......@@ -454,7 +475,7 @@ namespace AMDiS {
ERROR("virtual function not implemented in this sub-class ");
return(0.0);
};
}
virtual void getRefSimplexCoords(VectorOfFixVecs<DimVec<double> > *coords) const = 0;
......@@ -495,7 +516,34 @@ namespace AMDiS {
* of the children is (level of the parent + 1). level_ is filled always by
* the traversal routines.
*/
unsigned char level_;
unsigned char level;
/** \brief
* This value is set only during dual traverse, and only if this is the smaller
* element of an element pair. Then the value describes the displacement of
* this element to the first element in the refinement hierachy starting from
* the bigger element of the element pair.
*
* Example, if the level difference of the two elements is 1:
*
* a
* / \
* b c
*
* The element b has displacement 0, the element c has displacement 1.
*
* Example, if the level difference of the two elements is 2:
*
* a
* / \
* / \
* b c
* / \ / \
* d e f g
*
* The displacements are as follows: d = 0, e = 1, f = 2, g = 3.
*/
int displacement;
/** \brief
* This ElInfo is the iChild-th child of the parent element.
......@@ -567,7 +615,14 @@ namespace AMDiS {
*/
int dimOfWorld;
VectorOfFixVecs<DimVec<double> > *subElemCoords;
/** \brief
* This is a transformation matrix used during dual traverse. It is set, if
* the current element is the smaller element of an element pair in the traverse.
* Then this matrix defines a mapping for points defined in barycentric
* coordinates on the larger element, to the barycentric coordinates of the smaller
* element.
*/
DimMat<double> *subElemCoordsMat;
// ===== static public members ================================================
public:
......
......@@ -23,7 +23,7 @@ namespace AMDiS {
macroElement_ = const_cast<MacroElement*>( mel);
element_ = const_cast<Element*>( mel->getElement());
parent_ = NULL;
level_ = 0;
level = 0;
int vertices = mesh_->getGeo(VERTEX);
......@@ -200,7 +200,7 @@ namespace AMDiS {
macroElement_ = elInfoOld->macroElement_;
fillFlag_ = elInfoOld->fillFlag_;
parent_ = elem;
level_ = elInfoOld->level_ + 1;
level = elInfoOld->level + 1;
iChild = ichild;
int neighbours = mesh_->getGeo(NEIGH);
......@@ -274,10 +274,11 @@ namespace AMDiS {
void ElInfo1d::getRefSimplexCoords(VectorOfFixVecs<DimVec<double> > *coords) const
{
(*coords)[0][0] = 0.0;
(*coords)[0][1] = 1.0;
(*coords)[1][0] = 1.0;
(*coords)[1][1] = 0.0;
(*coords)[0][0] = 1.0;
(*coords)[0][1] = 0.0;
(*coords)[1][0] = 0.0;
(*coords)[1][1] = 1.0;
}
void ElInfo1d::getSubElementCoords(VectorOfFixVecs<DimVec<double> > *coords,
......
......@@ -43,7 +43,9 @@ namespace AMDiS {
/** \brief
* Constructor. Calls ElInfo's protected Constructor.
*/
ElInfo1d(Mesh* aMesh) : ElInfo(aMesh) {};
ElInfo1d(Mesh* aMesh)
: ElInfo(aMesh)
{}
/** \brief
* 1-dimensional realisation of ElInfo's fillElInfo method.
......@@ -77,7 +79,7 @@ namespace AMDiS {
int getSideOfNeighbour(int i) {
return (i + 1) % 2;
};
}
void getRefSimplexCoords(VectorOfFixVecs<DimVec<double> > *coords) const;
......
......@@ -36,7 +36,7 @@ namespace AMDiS {
macroElement_ = const_cast<MacroElement*>(mel);
element_ = const_cast<Element*>(mel->getElement());
parent_ = NULL;
level_ = 0;
level = 0;
if (fillFlag_.isSet(Mesh::FILL_COORDS) ||
fillFlag_.isSet(Mesh::FILL_DET) ||
......@@ -173,7 +173,7 @@ namespace AMDiS {
macroElement_ = elInfoOld->macroElement_;
fillFlag_ = fill_flag;
parent_ = elem;
level_ = elInfoOld->level_ + 1;
level = elInfoOld->level + 1;
iChild = ichild;
if (fillFlag_.isSet(Mesh::FILL_COORDS) ||
......
......@@ -16,7 +16,7 @@ namespace AMDiS {
void ElInfo3d::fillMacroInfo(const MacroElement * mel)
{
FUNCNAME("ElInfo3d::fillMacroInfo");
FUNCNAME("ElInfo3d::fillMacroInfo()");
Element *nb;
MacroElement *mnb;
Flag fill_opp_coords;
......@@ -24,7 +24,7 @@ namespace AMDiS {
macroElement_ = const_cast<MacroElement*>( mel);
element_ = const_cast<Element*>( mel->getElement());
parent_ = NULL;
level_ = 0;
level = 0;
el_type = const_cast<MacroElement*>(mel)->getElType();
int vertices = mesh_->getGeo(VERTEX);
......@@ -372,7 +372,7 @@ namespace AMDiS {
macroElement_ = elInfoOld->macroElement_;
fillFlag_ = fillFlag__local;
parent_ = el_old;
level_ = elInfoOld->level_ + 1;
level = elInfoOld->level + 1;
iChild = ichild;
int el_type_local = (dynamic_cast<const ElInfo3d*>(elInfoOld))->getType();
el_type = (el_type_local + 1) % 3;
......
......@@ -38,6 +38,7 @@ namespace AMDiS {
// ===========================================================================
class Mesh;
template<typename T> class DimMat;
template<typename T> class WorldVector;
template<typename T> class WorldMatrix;
......@@ -414,21 +415,26 @@ namespace AMDiS {
*/
DimVec(int dim, InitType initType=NO_INIT)
: FixVec<T,PARTS>(dim, initType)
{};
{}
/** \brief
* Calls the corresponding constructor of FixVec
*/
DimVec(int dim, InitType initType, T* ini)
: FixVec<T,PARTS>(dim, initType, ini)
{};
{}
/** \brief
* Calls the corresponding constructor of FixVec
*/
DimVec(int dim, InitType initType, const T& ini)
: FixVec<T,PARTS>(dim