diff --git a/AMDiS/src/DualTraverse.cc b/AMDiS/src/DualTraverse.cc
index ad1f1d43b094e40a4a17a6c7164d46b83b4859dd..1a9585b3ef50309826087df1003f8a651e6f904f 100644
--- a/AMDiS/src/DualTraverse.cc
+++ b/AMDiS/src/DualTraverse.cc
@@ -55,14 +55,8 @@ 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) {
TEST_EXIT(*elInfo1 == *elInfo2)("invalid dual traverse\n");
@@ -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));
+ *elInfo1 = stack1.traverseNext(*elInfo1);
}
if (inc2) {
- do {
- *elInfo2 = stack2.traverseNext(*elInfo2);
- } while (*elInfo2 != NULL && skipEl2(*elInfo2));
+ *elInfo2 = stack2.traverseNext(*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;
}
}
diff --git a/AMDiS/src/DualTraverse.h b/AMDiS/src/DualTraverse.h
index edd62d995c27314f88dd34f4ba692ddaa224bd02..38cc3832f2a495616c0c64640c635b828c923e1e 100644
--- a/AMDiS/src/DualTraverse.h
+++ b/AMDiS/src/DualTraverse.h
@@ -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;
};
}
diff --git a/AMDiS/src/ElInfo.cc b/AMDiS/src/ElInfo.cc
index 4e064fe7e5427829a100c80e00f10321197a263c..9470d9d21cf4867b41266aef3f49f4ae55319c7b 100644
--- a/AMDiS/src/ElInfo.cc
+++ b/AMDiS/src/ElInfo.cc
@@ -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;
}
}
diff --git a/AMDiS/src/ElInfo.h b/AMDiS/src/ElInfo.h
index f72bf4c5ad3e7f1b2171f9c124c5b136915614e0..3a9dc0ed3b416da979e172be5c7eb1dd626c9565 100644
--- a/AMDiS/src/ElInfo.h
+++ b/AMDiS/src/ElInfo.h
@@ -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,9 +412,8 @@ 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,
- WorldVector<double>* world) const;
+ virtual const WorldVector<double> *coordToWorld(const DimVec<double>& lambda,
+ WorldVector<double>* world) const;
/** \brief
@@ -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:
diff --git a/AMDiS/src/ElInfo1d.cc b/AMDiS/src/ElInfo1d.cc
index 9559b1a6881222b79826dc9283a273a58e083f8d..9b6fbeb6635751cd87b78f7918c1106eccbfc87c 100644
--- a/AMDiS/src/ElInfo1d.cc
+++ b/AMDiS/src/ElInfo1d.cc
@@ -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,
diff --git a/AMDiS/src/ElInfo1d.h b/AMDiS/src/ElInfo1d.h
index 57b75cbe1d3303bed72f228348976e0d10801b18..728118ff7474724a9acb8d3e467b1a70322609ce 100644
--- a/AMDiS/src/ElInfo1d.h
+++ b/AMDiS/src/ElInfo1d.h
@@ -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;
diff --git a/AMDiS/src/ElInfo2d.cc b/AMDiS/src/ElInfo2d.cc
index a8ced02b9f856f6aeb95f368ca4d67b095b82d0c..f3f33d18cd344159077ba2fb0cf932241008a803 100644
--- a/AMDiS/src/ElInfo2d.cc
+++ b/AMDiS/src/ElInfo2d.cc
@@ -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) ||
diff --git a/AMDiS/src/ElInfo3d.cc b/AMDiS/src/ElInfo3d.cc
index 2e41a93e13c93caac1602fe733a2994826776fe0..4b8589abbd834bc701771c938bc8c2f6ea4dc80f 100644
--- a/AMDiS/src/ElInfo3d.cc
+++ b/AMDiS/src/ElInfo3d.cc
@@ -16,15 +16,15 @@ namespace AMDiS {
void ElInfo3d::fillMacroInfo(const MacroElement * mel)
{
- FUNCNAME("ElInfo3d::fillMacroInfo");
- Element *nb;
+ FUNCNAME("ElInfo3d::fillMacroInfo()");
+ Element *nb;
MacroElement *mnb;
- Flag fill_opp_coords;
+ Flag fill_opp_coords;
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;
diff --git a/AMDiS/src/FixVec.h b/AMDiS/src/FixVec.h
index 975f2e70d583eb956a7c182285dad8780c5aa645..0ccaf6c4c721263a4cc0ab4f92494ea808725bd5 100644
--- a/AMDiS/src/FixVec.h
+++ b/AMDiS/src/FixVec.h
@@ -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, initType, ini)
- {};
+ {}
+
+ /** \brief
+ * Multplication of a matrix with a vector.
+ */
+ void multMatrixVec(DimMat<T> &m, DimVec<T> &v);
};
// ===========================================================================
@@ -450,7 +456,7 @@ namespace AMDiS {
*/
DimMat(int dim, InitType initType = NO_INIT)
: Matrix<T>(dim + 1, dim + 1)
- {};
+ {}
/** \brief
* Calls the corresponding constructor of VectorOfFixVecs
@@ -461,7 +467,7 @@ namespace AMDiS {
TEST_EXIT_DBG(initType == DEFAULT_VALUE)
("wrong initType or wrong initializer\n");
this->set(ini);
- };
+ }
/** \brief
* Calls the corresponding constructor of VectorOfFixVecs
@@ -471,7 +477,7 @@ namespace AMDiS {
{
TEST_EXIT_DBG(initType == VALUE_LIST)("wrong initType or wrong initializer\n");
setValues(ini);
- };
+ }
};
// ===========================================================================
diff --git a/AMDiS/src/FixVec.hh b/AMDiS/src/FixVec.hh
index 6536f374e170cf668133ae455d57bf027cc6c801..c46166b2cfb85df4ad9ea19408300be4b0ed749a 100644
--- a/AMDiS/src/FixVec.hh
+++ b/AMDiS/src/FixVec.hh
@@ -1,5 +1,19 @@
namespace AMDiS {
+ template<typename T>
+ void DimVec<T>::multMatrixVec(DimMat<T> &m, DimVec<T> &v)
+ {
+ T *mIt, *thisIt;
+ for (thisIt = this->begin(), mIt = m.begin();
+ thisIt != this->end();
+ thisIt++) {
+ *thisIt = 0;
+ for (T* vIt = v.begin(); vIt != v.end(); vIt++, mIt++) {
+ *thisIt += *vIt * *mIt;
+ }
+ }
+ }
+
template<typename T, GeoIndex d>
double absteukl(const FixVec<T,d>& a,const FixVec<T,d>& b)
{
@@ -22,14 +36,13 @@ namespace AMDiS {
T *mIt, *thisIt;
for (thisIt = this->begin(), mIt = m.begin();
thisIt != this->end();
- thisIt++)
- {
- *thisIt = 0;
-
- for (T* vIt = v.begin(); vIt != v.end(); vIt++, mIt++) {
- *thisIt += *vIt * *mIt;
- }
+ thisIt++) {
+ *thisIt = 0;
+
+ for (T* vIt = v.begin(); vIt != v.end(); vIt++, mIt++) {
+ *thisIt += *vIt * *mIt;
}
+ }
}
template<typename T>
@@ -39,6 +52,7 @@ namespace AMDiS {
for (int j = i + 1; j < this->getSize(); j++)
if (abs((*this)[i][j]) > DBL_TOL || abs((*this)[j][i]) > DBL_TOL)
return(false);
+
return(true);
}
@@ -94,10 +108,9 @@ namespace AMDiS {
T *m1It, *m2It;
for (m1It = m1.begin(), m2It = m2.begin();
m1It != m1.end();
- m1It++, m2It++)
- {
- *m1It -= *m2It;
- }
+ m1It++, m2It++) {
+ *m1It -= *m2It;
+ }
return m1;
}
@@ -108,10 +121,9 @@ namespace AMDiS {
T* m1It, *m2It;
for (m1It = m1.begin(), m2It = m2.begin();
m1It != m1.end();
- m1It++, m2It++)
- {
- *m1It += *m2It;
- }
+ m1It++, m2It++) {
+ *m1It += *m2It;
+ }
return m1;
}
diff --git a/AMDiS/src/Marker.cc b/AMDiS/src/Marker.cc
index c817a90a27122392e80e541fe13488732ff4a3f1..c0e055736ce3a6bb387bfb786544777f215b8476 100644
--- a/AMDiS/src/Marker.cc
+++ b/AMDiS/src/Marker.cc
@@ -49,7 +49,7 @@ namespace AMDiS {
double lError = el->getEstimation(row);
if (adaptInfo->isRefinementAllowed(row == -1 ? 0 : row) && lError > markRLimit) {
- if ((maxRefineLevel == -1) || (static_cast<int>(elInfo->getLevel()) < maxRefineLevel)) {
+ if ((maxRefineLevel == -1) || (elInfo->getLevel() < maxRefineLevel)) {
setMark(el, adaptInfo->getRefineBisections(row == -1 ? 0 : row));
}
} else {
diff --git a/AMDiS/src/MatrixVector.cc b/AMDiS/src/MatrixVector.cc
index 7cbaadb260d3bfd06073122e8875f1968b84604a..a2b9c54fe43ab331363d78ba84660df3c7fbe0e4 100644
--- a/AMDiS/src/MatrixVector.cc
+++ b/AMDiS/src/MatrixVector.cc
@@ -2,7 +2,20 @@
#include "DOFVector.h"
namespace AMDiS {
+
+ template<>
+ void Matrix<double>::invert2x2()
+ {
+ TEST_EXIT(rows == 2)("WRONG\n");
+ TEST_EXIT(cols == 2)("WRONG\n");
+ double det = 1 / (valArray[0] * valArray[3] - valArray[1] * valArray[2]);
+ double tmp = valArray[0];
+ valArray[0] = det * valArray[3];
+ valArray[3] = det * tmp;
+ valArray[1] *= -det;
+ valArray[2] *= -det;
+ }
}
diff --git a/AMDiS/src/MatrixVector.h b/AMDiS/src/MatrixVector.h
index db9bea122a9b7031d0ba6592d830ab44dbde6843..a2a04d23a08c51190dd6f489cb955109dd41470f 100644
--- a/AMDiS/src/MatrixVector.h
+++ b/AMDiS/src/MatrixVector.h
@@ -54,12 +54,11 @@ namespace AMDiS {
valArray = NULL;
else
valArray = new T[size];
- };
+ }
- inline bool used() const
- {
+ inline bool used() const {
return (valArray != NULL);
- };
+ }
/** \brief
* Change the size of the vector to newSize.
@@ -71,7 +70,7 @@ namespace AMDiS {
valArray = new T[newSize];
size = newSize;
}
- };
+ }
/** \brief
* Copy constructor.
@@ -81,14 +80,14 @@ namespace AMDiS {
{
valArray = new T[rhs.size];
*this = rhs; // uses operator=()
- };
+ }
/** \brief
* Destructor.
*/
virtual ~Vector() {
delete [] valArray;
- };
+ }
/** \brief
* Assignement operator
@@ -96,14 +95,14 @@ namespace AMDiS {
inline const Vector<T>& operator=(const Vector<T>& rhs) {
TEST_EXIT_DBG(rhs.size == size)("invalid size\n");
T *rhsIt, *thisIt;
- for(rhsIt = rhs.begin(), thisIt = this->begin();
- rhsIt != rhs.end();
- ++rhsIt, ++thisIt)
- {
- *thisIt = *rhsIt;
- }
+ for (rhsIt = rhs.begin(), thisIt = this->begin();
+ rhsIt != rhs.end();
+ ++rhsIt, ++thisIt) {
+ *thisIt = *rhsIt;
+ }
+
return *this;
- };
+ }
/** \brief
* Assignement operator
@@ -117,7 +116,7 @@ namespace AMDiS {
*thisIt = scal;
}
return *this;
- };
+ }
/** \brief
* Assignement operator
@@ -132,14 +131,14 @@ namespace AMDiS {
*thisIt = *vecIt;
}
return *this;
- };
+ }
/** \brief
* Sets all entries to scal.
*/
inline const Vector<T>& set(const T& scal) {
return *this = scal;
- };
+ }
/** \brief
* Sets all entries.
@@ -154,7 +153,7 @@ namespace AMDiS {
*thisIt = *valuesIt;
}
return *this;
- };
+ }
/** \brief
* Sets all entries.
@@ -178,7 +177,7 @@ namespace AMDiS {
if(*thisIt != *rhsIt) return false;
}
return true;
- };
+ }
/** \brief
* Comparison operator.
@@ -193,7 +192,7 @@ namespace AMDiS {
inline T& operator[](int i) {
TEST_EXIT_DBG(i < size && i >= 0)("invalid index\n");
return valArray[i];
- };
+ }
/** \brief
* Access to the i-th vector element for const vectors.
@@ -201,35 +200,35 @@ namespace AMDiS {
inline const T& operator[] (int i) const {
TEST_EXIT_DBG(i < size && i >= 0)("invalid index\n");
return valArray[i];
- };
+ }
/** \brief
* Returns pointer to the first vector element.
*/
inline T *begin() const {
return valArray;
- };
+ }
/** \brief
* Returns pointer after the last vector element.
*/
inline T *end() const {
return valArray + size;
- };
+ }
/** \brief
* Returns \ref size.
*/
virtual int getSize() const {
return size;
- };
+ }
/** \brief
* Returns \ref valArray as T-array
*/
inline T *getValArray() {
return valArray;
- };
+ }
void print() const {
std::cout << this->size << " vector: " << std::endl;
@@ -237,23 +236,23 @@ namespace AMDiS {
std::cout << this->valArray[i] << " ";
}
std::cout << std::endl;
- };
+ }
// ===== Serializable implementation =====
void serialize(std::ostream &out) {
out.write(reinterpret_cast<const char*>(&size), sizeof(int));
out.write(reinterpret_cast<const char*>(valArray), size * sizeof(T));
- };
+ }
void deserialize(std::istream &in) {
in.read(reinterpret_cast<char*>(&size), sizeof(int));
in.read(reinterpret_cast<char*>(valArray), size * sizeof(T));
- };
+ }
std::string getTypeName() const {
return "Vector";
- };
+ }
protected:
/** \brief
@@ -289,7 +288,7 @@ namespace AMDiS {
: Vector<T>(r*c),
rows(r),
cols(c)
- {};
+ {}
/** \brief
* Changes the size of the matrix to newRows x newCols.
@@ -300,20 +299,20 @@ namespace AMDiS {
rows = newRows;
cols = newCols;
}
- };
+ }
/** \brief
* Assignement operator.
*/
inline const Matrix<T>& operator=(const T& scal) {
return static_cast<const Matrix<T>&>(Vector<T>::operator=(scal));
- };
+ }
inline bool operator==(const Matrix<T>& rhs) const {
if (rows != rhs.getNumRows()) return false;
if (cols != rhs.getNumCols()) return false;
return Vector<T>::operator == (rhs);
- };
+ }
/** \brief
* Comparison operator.
@@ -327,21 +326,21 @@ namespace AMDiS {
*/
inline T *operator[](int i) {
return this->valArray + cols * i;
- };
+ }
/** \brief
* Acces to i-th matrix row for constant matrices.
*/
inline const T *operator[](int i) const {
return this->valArray + cols * i;
- };
+ }
/** \brief
* Returns \ref rows.
*/
inline int getNumRows() const {
return rows;
- };
+ }
/** \brief
* Return \ref cols.
@@ -362,7 +361,7 @@ namespace AMDiS {
*/
inline T *end() const {
return this->valArray + (cols * rows);
- };
+ }
void print() const {
std::cout << this->rows << " x " << this->cols << " matrix: " << std::endl;
@@ -372,7 +371,9 @@ namespace AMDiS {
}
std::cout << std::endl;
}
- };
+ }
+
+ void invert2x2();
protected:
/** \brief
diff --git a/AMDiS/src/ProblemVec.cc b/AMDiS/src/ProblemVec.cc
index 89f7b4da61f213014449ecd5a2cab817bdc7cc32..65ddfab53658b68e6bda7793cc8ad1dd47ff53fa 100644
--- a/AMDiS/src/ProblemVec.cc
+++ b/AMDiS/src/ProblemVec.cc
@@ -729,21 +729,19 @@ namespace AMDiS {
matrix->getBoundaryManager()->initMatrix(matrix);
if (componentSpaces[i] == componentSpaces[j]) {
- // std::cout << "--- assembleOnOneMesh ---\n";
assembleOnOneMesh(componentSpaces[i],
assembleFlag,
assembleMatrix ? matrix : NULL,
((i == j) && asmVector) ? rhs_->getDOFVector(i) : NULL);
- // std::cout << "--- Finished ---\n";
- // std::cout << "--- Dim: " << matrix->getUsedSize() << "x" << matrix->getNumCols() << "---\n";
} else {
- // std::cout << "--- assembleOnDifMesh ---\n";
assembleOnDifMeshes(componentSpaces[i], componentSpaces[j],
assembleFlag,
assembleMatrix ? matrix : NULL,
((i == j) && asmVector) ? rhs_->getDOFVector(i) : NULL);
- // std::cout << "--- Finished ---\n";
- // std::cout << "--- Dim: " << matrix->getUsedSize() << "x" << matrix->getNumCols() << "---\n";
+ TEST_EXIT(matrix->getUsedSize() == componentSpaces[i]->getAdmin()->getUsedSize())
+ ("Assembled matrix has wrong dimension!\n");
+ TEST_EXIT(matrix->getNumCols() == componentSpaces[j]->getAdmin()->getUsedSize())
+ ("Assembled matrix has wrong dimension!\n");
}
if (assembleMatrix && matrix->getBoundaryManager())
@@ -1035,7 +1033,6 @@ namespace AMDiS {
ElInfo *rowElInfo, *colElInfo;
ElInfo *largeElInfo, *smallElInfo;
- dualTraverse.setFillSubElemInfo(true);
bool cont = dualTraverse.traverseFirst(rowMesh, colMesh, -1, -1,
assembleFlag, assembleFlag,
&rowElInfo, &colElInfo,
@@ -1048,14 +1045,7 @@ namespace AMDiS {
if (smallElInfo == rowElInfo) {
std::cout << "Row = small\n";
ERROR_EXIT("NOCH EIN PAAR GEDANKEN MACHEN!\n");
- } else {
- std::cout << "Row = large\n";
- }
-
- if (largeElInfo == colElInfo)
- std::cout << "Col = large\n";
- else
- std::cout << "Col = small\n";
+ }
if (useGetBound_) {
basisFcts->getBound(rowElInfo, bound);
diff --git a/AMDiS/src/Quadrature.h b/AMDiS/src/Quadrature.h
index b44f5eaf7910ddc560370dd23f07e08617ea41c4..0f96518566e19f495fc150cc5ddfd535eec532d4 100644
--- a/AMDiS/src/Quadrature.h
+++ b/AMDiS/src/Quadrature.h
@@ -81,7 +81,7 @@ namespace AMDiS {
n_points(n_points_),
lambda(lambda_),
w(w_)
- {};
+ {}
public:
/** \brief
diff --git a/AMDiS/src/ScalableQuadrature.cc b/AMDiS/src/ScalableQuadrature.cc
index 6db04c06df65cd80982a760cd5706e6c495a32b3..aa2f6e18a7c2545091a4d5b507405663d8c44fa3 100644
--- a/AMDiS/src/ScalableQuadrature.cc
+++ b/AMDiS/src/ScalableQuadrature.cc
@@ -6,15 +6,17 @@ namespace AMDiS {
ScalableQuadrature::ScalableQuadrature(Quadrature *quadrature)
: Quadrature(*quadrature)
{
- /**
- * Change name of quadrature.
- */
+ // Change name of quadrature.
name = "Scalable" + getName();
- /**
- * copy quadrature->lambda to oldLambda
- */
+ // copy quadrature->lambda to oldLambda
oldLambda = NEW VectorOfFixVecs<DimVec<double> >(*(quadrature->getLambda()));
+
+ // First assume all quadrature points to be valid.
+ valid.resize(n_points);
+ for (int i = 0; i < n_points; i++) {
+ valid[i] = true;
+ }
}
void ScalableQuadrature::scaleQuadrature(const SubElInfo& subElInfo)
@@ -36,7 +38,6 @@ namespace AMDiS {
* (n_0, n_1, n_2) = l_0 * x_0 + l_1 * x_1 + l_2 * x_2 .
*/
for (int iq = 0; iq < n_points; iq++) {
-
for (int i = 0; i <= dim; i++) {
/**
@@ -49,7 +50,30 @@ namespace AMDiS {
}
(*lambda)[iq][i] = l;
+
+ // If one component is less than zero, the scaled quadrature point
+ // is not inside the target element and is therefore set to be invalid.
+ if (l < 0.0) {
+ valid[iq] = false;
+ }
}
}
}
+
+ void ScalableQuadrature::scaleQuadrature(DimMat<double> *scalMat)
+ {
+ for (int iq = 0; iq < n_points; iq++) {
+ (*lambda)[iq].multMatrixVec(*scalMat, *getOldLambda(iq));
+
+ for (int i = 0; i <= dim; i++) {
+ // If one component is less than zero, the scaled quadrature point
+ // is not inside the target element and is therefore set to be invalid.
+ if ((*lambda)[iq][i] < 0.0) {
+ valid[iq] = false;
+ break;
+ }
+ }
+ }
+ }
+
}
diff --git a/AMDiS/src/ScalableQuadrature.h b/AMDiS/src/ScalableQuadrature.h
index 7797f6b52272bc88ca64da4827d0b3e46b7fe666..5c21a5ab6495df9d1975bffeddfe4996bda0ed69 100644
--- a/AMDiS/src/ScalableQuadrature.h
+++ b/AMDiS/src/ScalableQuadrature.h
@@ -45,24 +45,30 @@ namespace AMDiS {
public:
MEMORY_MANAGED(ScalableQuadrature);
- /**
+ /** \brief
* Constructor
*/
ScalableQuadrature(Quadrature *quadrature);
- /**
+ /** \brief
* Destructor
*/
~ScalableQuadrature() {
DELETE oldLambda;
}
- /**
+ /** \brief
* Manipulates the quadrature points for the assemblage of a subelement.
*/
void scaleQuadrature(const SubElInfo& subElInfo);
- /**
+ /** \brief
+ * Scales the quadrature points using a matrix, that defines the transformation
+ * from one element to another one.
+ */
+ void scaleQuadrature(DimMat<double> *scalMat);
+
+ /** \brief
* Get b-th coordinate of the a-th original quadrature point.
*/
inline const double getOldLambda(int a, int b) const {
@@ -70,13 +76,38 @@ namespace AMDiS {
return (*oldLambda)[a][b];
else
return 0.0;
- };
+ }
+
+ /** \brief
+ * Get the a-th original quadrature point.
+ */
+ inline DimVec<double>* getOldLambda(int a) const {
+ if (oldLambda)
+ return &((*oldLambda)[a]);
+ else
+ ERROR_EXIT("Should not happen!\n");
+ }
+
+ /** \brief
+ * Returns, if the scaled quadrature point is within the target element (true),
+ * or outside of the element (false).
+ */
+ inline bool isValid(int iq) const {
+ return valid[iq];
+ }
protected:
- /**
+ /** \brief
* Original quadrature points.
*/
VectorOfFixVecs<DimVec<double> > *oldLambda;
+
+ /** \brief
+ * If quadrature points are scaled from on element to another one, they may
+ * become invalide, i.e., there are outside of the target element. In this
+ * case, the corresponding element of this field is false, otherwise true.
+ */
+ std::vector<bool> valid;
};
}
diff --git a/AMDiS/src/SubElInfo.h b/AMDiS/src/SubElInfo.h
index f50723f929e3e3d636d4dfd6cc1cee36ca3fb07e..22c7e7ffb0afa1a75bdcc7d85e9be21fdc5e5a4c 100644
--- a/AMDiS/src/SubElInfo.h
+++ b/AMDiS/src/SubElInfo.h
@@ -41,7 +41,7 @@ namespace AMDiS {
*/
~SubElInfo() {
DELETE lambda;
- };
+ }
/**
* Get b-th coordinate of the a-th vertex of subelement (barycentric
@@ -52,28 +52,28 @@ namespace AMDiS {
return (*lambda)[a][b];
else
return 0.0;
- };
+ }
/**
* Get coordinates of a-th vertex of subelement (barycentric coordinates).
*/
inline const DimVec<double>& getLambda(int a) const {
return (*lambda)[a];
- };
+ }
/**
* Get coordinates of all vertices of subelement (barycentric coordinates).
*/
inline VectorOfFixVecs<DimVec<double> > *getLambda() const {
return lambda;
- };
+ }
/**
* Get determinant corresponding to subelement.
*/
inline const double getDet() const {
return det;
- };
+ }
protected:
diff --git a/AMDiS/src/Traverse.cc b/AMDiS/src/Traverse.cc
index ae40b69a0341c9617276aa7c215e84f81f7e6d14..2b55cbfabfecc7d8a3c7311cb9bedc843215caa1 100644
--- a/AMDiS/src/Traverse.cc
+++ b/AMDiS/src/Traverse.cc
@@ -40,7 +40,7 @@ namespace AMDiS {
traverse_mel = NULL;
stack_used = 0;
-
+
return(traverseNext(NULL));
}
@@ -196,10 +196,9 @@ namespace AMDiS {
if (flag.isSet(Mesh::CALL_MG_LEVEL)) {
- mg_level = (static_cast<int>(elinfo->getLevel()) + mesh->getDim() - 1) /
- mesh->getDim();
+ mg_level = (elinfo->getLevel() + mesh->getDim() - 1) / mesh->getDim();
- if (mg_level > static_cast<int>(level)) {
+ if (mg_level > level) {
return 0;
}
@@ -333,7 +332,7 @@ namespace AMDiS {
{
FUNCNAME("TraverseStack::traverseLeafElement()");
- Element *el;
+ Element *el = NULL;
if (stack_used == 0) { /* first call */
currentMacro = traverse_mesh->firstMacroElement();
@@ -343,11 +342,11 @@ namespace AMDiS {
currentMacro++;
}
- if (currentMacro == traverse_mesh->endOfMacroElements())
+ if (currentMacro == traverse_mesh->endOfMacroElements()) {
return NULL;
+ }
traverse_mel = *currentMacro;
-
stack_used = 1;
elinfo_stack[stack_used]->fillMacroInfo(traverse_mel);
info_stack[stack_used] = 0;
@@ -364,6 +363,10 @@ namespace AMDiS {
((info_stack[stack_used] >= 2) || (el->getFirstChild() == NULL))) {
stack_used--;
el = elinfo_stack[stack_used]->getElement();
+
+ if (calcDisplacement) {
+ displacementStackPos--;
+ }
}
/* goto next macro element */
@@ -373,8 +376,9 @@ namespace AMDiS {
} while ((currentMacro != traverse_mesh->endOfMacroElements()) &&
((*currentMacro)->getIndex() % maxThreads_ != myThreadId_));
- if (currentMacro == traverse_mesh->endOfMacroElements())
+ if (currentMacro == traverse_mesh->endOfMacroElements()) {
return NULL;
+ }
traverse_mel = *currentMacro;
stack_used = 1;
@@ -390,14 +394,23 @@ namespace AMDiS {
/* go down tree until leaf */
while (el->getFirstChild()) {
- if (stack_used >= stack_size-1)
+ if (stack_used >= stack_size - 1)
enlargeTraverseStack();
int i = info_stack[stack_used];
- el = const_cast<Element*>( ((i == 0) ? el->getFirstChild() : el->getSecondChild()));
+ el = const_cast<Element*>(((i == 0) ? el->getFirstChild() : el->getSecondChild()));
info_stack[stack_used]++;
elinfo_stack[stack_used + 1]->fillElInfo(i, elinfo_stack[stack_used]);
stack_used++;
+
+ if (calcDisplacement) {
+ displacementStackPos++;
+ if (displacementStackPos < static_cast<int>(displacementStack.size())) {
+ displacementStack[displacementStackPos]++;
+ } else {
+ displacementStack.push_back(2 * (displacementStack[displacementStackPos - 1] - 1) + 1);
+ }
+ }
TEST_EXIT_DBG(stack_used < stack_size)
@@ -406,6 +419,10 @@ namespace AMDiS {
info_stack[stack_used] = 0;
}
+
+ if (calcDisplacement) {
+ elinfo_stack[stack_used]->setDisplacement(displacementStack[displacementStackPos]);
+ }
return elinfo_stack[stack_used];
}
@@ -432,58 +449,57 @@ namespace AMDiS {
ElInfo* TraverseStack::traverseMultiGridLevel()
{
- FUNCNAME("TraverseStack::traverseMultiGridLevel");
- Element *el;
- int i;
-
- if (stack_used == 0) /* first call */
- {
- currentMacro = traverse_mesh->firstMacroElement();
- traverse_mel = *currentMacro;
- if (traverse_mel == NULL) return NULL;
+ FUNCNAME("TraverseStack::traverseMultiGridLevel()");
- stack_used = 1;
- elinfo_stack[stack_used]->fillMacroInfo(traverse_mel);
- info_stack[stack_used] = 0;
-
- if((elinfo_stack[stack_used]->getLevel() == traverse_level) ||
- (elinfo_stack[stack_used]->getLevel() < traverse_level &&
- elinfo_stack[stack_used]->getElement()->isLeaf()))
- return(elinfo_stack[stack_used]);
- }
+ if (stack_used == 0) { /* first call */
+ currentMacro = traverse_mesh->firstMacroElement();
+ traverse_mel = *currentMacro;
+ if (traverse_mel == NULL) return NULL;
+
+ stack_used = 1;
+ elinfo_stack[stack_used]->fillMacroInfo(traverse_mel);
+ info_stack[stack_used] = 0;
+
+ if ((elinfo_stack[stack_used]->getLevel() == traverse_level) ||
+ (elinfo_stack[stack_used]->getLevel() < traverse_level &&
+ elinfo_stack[stack_used]->getElement()->isLeaf()))
+ return(elinfo_stack[stack_used]);
+ }
- el = elinfo_stack[stack_used]->getElement();
+ Element *el = elinfo_stack[stack_used]->getElement();
/* go up in tree until we can go down again */
- while((stack_used > 0) &&
- ((info_stack[stack_used] >= 2) || (el->getFirstChild()==NULL)))
- {
- stack_used--;
- el = elinfo_stack[stack_used]->getElement();
- }
+ while ((stack_used > 0) &&
+ ((info_stack[stack_used] >= 2) || (el->getFirstChild()==NULL))) {
+ stack_used--;
+ el = elinfo_stack[stack_used]->getElement();
+ }
/* goto next macro element */
if (stack_used < 1) {
currentMacro++;
- if (currentMacro == traverse_mesh->endOfMacroElements()) return(NULL);
- traverse_mel = *currentMacro;
+ if (currentMacro == traverse_mesh->endOfMacroElements())
+ return(NULL);
+ traverse_mel = *currentMacro;
stack_used = 1;
elinfo_stack[stack_used]->fillMacroInfo(traverse_mel);
info_stack[stack_used] = 0;
- if((elinfo_stack[stack_used]->getLevel() == traverse_level) ||
- (elinfo_stack[stack_used]->getLevel() < traverse_level &&
- elinfo_stack[stack_used]->getElement()->isLeaf()))
+ if ((elinfo_stack[stack_used]->getLevel() == traverse_level) ||
+ (elinfo_stack[stack_used]->getLevel() < traverse_level &&
+ elinfo_stack[stack_used]->getElement()->isLeaf()))
return(elinfo_stack[stack_used]);
}
/* go down tree */
- if (stack_used >= stack_size-1) enlargeTraverseStack();
- i = info_stack[stack_used];
+ if (stack_used >= stack_size - 1)
+ enlargeTraverseStack();
+
+ int i = info_stack[stack_used];
info_stack[stack_used]++;
elinfo_stack[stack_used + 1]->fillElInfo(i, elinfo_stack[stack_used]);
@@ -495,9 +511,9 @@ namespace AMDiS {
info_stack[stack_used] = 0;
- if((elinfo_stack[stack_used]->getLevel() == traverse_level) ||
- (elinfo_stack[stack_used]->getLevel() < traverse_level &&
- elinfo_stack[stack_used]->getElement()->isLeaf()))
+ if ((elinfo_stack[stack_used]->getLevel() == traverse_level) ||
+ (elinfo_stack[stack_used]->getLevel() < traverse_level &&
+ elinfo_stack[stack_used]->getElement()->isLeaf()))
return(elinfo_stack[stack_used]);
return traverseMultiGridLevel();
@@ -624,7 +640,9 @@ namespace AMDiS {
while(elinfo_stack[stack_used]->getElement()->getFirstChild()
&& (info_stack[stack_used] < 2)) {
- if (stack_used >= stack_size-1) enlargeTraverseStack();
+ if (stack_used >= stack_size-1)
+ enlargeTraverseStack();
+
i = info_stack[stack_used];
info_stack[stack_used]++;
elinfo_stack[stack_used + 1]->fillElInfo(i, elinfo_stack[stack_used]);
@@ -893,7 +911,7 @@ namespace AMDiS {
elinfo_stack[stack_used]->testFlag(Mesh::FILL_NEIGH);
el = dynamic_cast<Triangle*>(const_cast<Element*>( elinfo_stack[stack_used]->getElement()));
sav_index = el->getIndex();
- sav_el = el;
+ sav_el = el;
/* first, goto to leaf level, if necessary... */
if (!(el->isLeaf()) && (neighbour < 2)) {
@@ -925,13 +943,14 @@ namespace AMDiS {
nb = neighbour;
- while (stack_used > 1)
- {
- stack_used--;
- nb = coarse_nb[info_stack[stack_used]][nb];
- if (nb == -1) break;
- TEST_EXIT_DBG(nb >= 0)("invalid coarse_nb %d\n",nb);
- }
+ while (stack_used > 1) {
+ stack_used--;
+ nb = coarse_nb[info_stack[stack_used]][nb];
+ if (nb == -1)
+ break;
+
+ TEST_EXIT_DBG(nb >= 0)("invalid coarse_nb %d\n",nb);
+ }
/****************************************************************************/
/* Now, goto neighbouring element at the local hierarchy entry */
@@ -969,8 +988,9 @@ namespace AMDiS {
el2 = dynamic_cast<Triangle*>(const_cast<Element*>( elinfo2->getElement()));
- if (stack_used >= stack_size-1)
+ if (stack_used >= stack_size-1) {
enlargeTraverseStack();
+ }
i = 2 - info_stack[stack_used];
info_stack[stack_used] = i+1;
elinfo_stack[stack_used+1]->fillElInfo(i, elinfo_stack[stack_used]);
@@ -987,7 +1007,7 @@ namespace AMDiS {
elinfo = elinfo_stack[stack_used];
el = dynamic_cast<Triangle*>(const_cast<Element*>( elinfo->getElement()));
- while(el->getFirstChild()) {
+ while (el->getFirstChild()) {
if (nb < 2) { /* go down one level in hierarchy */
if (stack_used >= stack_size-1)
@@ -1038,6 +1058,7 @@ namespace AMDiS {
TEST_EXIT_DBG(elinfo->getNeighbour(opp_vertex) == old_elinfo->getElement())
("didn't succeed !?!?!?");
}
+
if (elinfo->getElement()->getFirstChild()) {
MSG(" looking for neighbour %d of element %d at %8X\n",
neighbour, old_elinfo->getElement()->getIndex(), old_elinfo->getElement());
@@ -1067,7 +1088,7 @@ namespace AMDiS {
TEST_EXIT_DBG(traverse_mesh->getDim() == 3)("update only in 3d\n");
for (int i = stack_used; i > 0; i--) {
- dynamic_cast<ElInfo3d*>( elinfo_stack[i])->update();
+ dynamic_cast<ElInfo3d*>(elinfo_stack[i])->update();
}
}
@@ -1085,4 +1106,18 @@ namespace AMDiS {
lowerLevel--;
}
}
+
+ void TraverseStack::startDisplacementCalculation(int level)
+ {
+ calcDisplacement = true;
+
+ displacementStack.resize(elinfo_stack[stack_used]->getLevel() - level + 1);
+ for (int i = 0; i < static_cast<int>(displacementStack.size()); i++) {
+ displacementStack[i] = 0;
+ }
+
+ displacementStackPos = displacementStack.size() - 1;
+ elinfo_stack[stack_used]->setDisplacement(0);
+ }
+
}
diff --git a/AMDiS/src/Traverse.h b/AMDiS/src/Traverse.h
index 00b764cf52a46d09e172cef1fa6e447f227a9fb5..e798cb3fba9acc375367cc9bca7a2c4335046f7c 100644
--- a/AMDiS/src/Traverse.h
+++ b/AMDiS/src/Traverse.h
@@ -70,6 +70,9 @@ namespace AMDiS {
stack_size(0),
stack_used(0),
save_stack_used(0),
+ displacementStack(0),
+ calcDisplacement(false),
+ displacementStackPos(-1),
myThreadId_(0),
maxThreads_(1)
{
@@ -77,7 +80,7 @@ namespace AMDiS {
#pragma omp critical
#endif
id = ElInfo::traverseIdCounter++;
- };
+ }
/** \brief
* Destructor
@@ -90,7 +93,7 @@ namespace AMDiS {
for (int i = 0; i < static_cast<int>(save_elinfo_stack.size()); i++) {
DELETE save_elinfo_stack[i];
}
- };
+ }
public:
/** \brief
@@ -150,6 +153,18 @@ namespace AMDiS {
void getCoordsInElem(const ElInfo *upperElInfo,
VectorOfFixVecs<DimVec<double> > *coords);
+ /** \brief
+ * Starts the calculation of the displacement relative to given level.
+ */
+ void startDisplacementCalculation(int level);
+
+ /** \brief
+ * Stops the calculation of displacement information.
+ */
+ void stopDisplacementCalculation() {
+ calcDisplacement = false;
+ }
+
/** \brief
* Is used for parallel mesh traverse.
*/
@@ -190,18 +205,18 @@ namespace AMDiS {
* the operator must be implemented (deep copy not flat copy!)
*/
void operator=(const TraverseStack& /*rhs*/) {
- FUNCNAME("TraverseStack::operator=");
+ FUNCNAME("TraverseStack::operator=()");
ERROR_EXIT("not implemented");
- };
+ }
/** \brief
* Avoids copy of a traverse stack. If copy should be possible
* the operator must be implemented (deep copy not flat copy!)
*/
TraverseStack(const TraverseStack&) {
- FUNCNAME("TraverseStack::TraverseStack(const TraverseStack&)");
+ FUNCNAME("TraverseStack::TraverseStack()");
ERROR_EXIT("not implemented");
- };
+ }
private:
/** \brief
@@ -247,9 +262,13 @@ namespace AMDiS {
std::vector<ElInfo*> elinfo_stack;
/** \brief
- *
+ * Stores for all level, which children of this element was already
+ * visited. If 0, no children were visited (this status occurs only
+ * in intermediate computations). If 1, only the first was vistied.
+ * If 2, both children of the element on this level were already
+ * visited.
*/
- std::vector<unsigned char> info_stack;
+ std::vector<int> info_stack;
/** \brief
*
@@ -276,6 +295,21 @@ namespace AMDiS {
*/
int id;
+ /** \brief
+ * Stack for counting the deplacement information for dual traverse.
+ */
+ std::vector<int> displacementStack;
+
+ /** \brief
+ * True, if \ref displacementStack should be calculated.
+ */
+ bool calcDisplacement;
+
+ /** \brief
+ * Current position within the \ref displacementStack.
+ */
+ int displacementStackPos;
+
/** \brief
* Is used for parallel mesh traverse. The thread with the id
* myThreadId is only allowed to access coarse elements, which id
@@ -320,12 +354,11 @@ namespace AMDiS {
{
TEST_EXIT(m)("No traverse without mesh!\n");
id = ElInfo::traverseIdCounter++;
- };
+ }
/** \brief
* Performs the recursive traversal
*/
- // int recursive(ElInfo*);
int recursive(ElInfoStack*);
private:
@@ -333,7 +366,7 @@ namespace AMDiS {
Flag flag;
- unsigned char level;
+ int level;
int (*el_fct)(ElInfo*);
diff --git a/AMDiS/src/ZeroOrderAssembler.cc b/AMDiS/src/ZeroOrderAssembler.cc
index f51fe81e0332fa5fcbc5add449752a97a2135468..5361c772ded2ca742ed0845b5c45a511adb325f9 100644
--- a/AMDiS/src/ZeroOrderAssembler.cc
+++ b/AMDiS/src/ZeroOrderAssembler.cc
@@ -8,8 +8,6 @@
#include "DOFVector.h"
#include "ElementMatrix.h"
#include "OpenMP.h"
-#include "SubElInfo.h"
-#include "ScalableQuadrature.h"
namespace AMDiS {
@@ -151,50 +149,48 @@ namespace AMDiS {
const ElInfo *colElInfo,
ElementMatrix *mat)
{
- const BasisFunction *psi = owner->getRowFESpace()->getBasisFcts();
- const BasisFunction *phi = owner->getColFESpace()->getBasisFcts();
+ const BasisFunction *psi = owner->getRowFESpace()->getBasisFcts(); // grosses Element
+ const BasisFunction *phi = owner->getColFESpace()->getBasisFcts(); // kleines Element
- double *phival = GET_MEMORY(double, nCol);
int nPoints = quadrature->getNumPoints();
double *c = GET_MEMORY(double, nPoints);
+ for (int iq = 0; iq < nPoints; iq++) {
+ c[iq] = 0.0;
+ }
+
int myRank = omp_get_thread_num();
std::vector<OperatorTerm*>::iterator termIt;
for (termIt = terms[myRank].begin(); termIt != terms[myRank].end(); ++termIt) {
(static_cast<ZeroOrderTerm*>((*termIt)))->getC(rowElInfo, nPoints, c);
}
-
- SubElInfo *subElInfo = NEW SubElInfo(colElInfo->getSubElemCoords(), rowElInfo);
- Quadrature psiQuadrature(*quadrature);
- ScalableQuadrature *scaledQuadrature = NEW ScalableQuadrature(&psiQuadrature);
- scaledQuadrature->scaleQuadrature(*subElInfo);
-
+
for (int iq = 0; iq < nPoints; iq++) {
- c[iq] *= rowElInfo->getDet();
+ c[iq] *= colElInfo->getDet();
+
+ for (int i = 0; i < 2; i++) { // iteration ueber col-phi basis funktionen
+ for (int j = 0; j < 2; j++) { // iteration ueber row-psi basis funktionen
+ if (j == 0) {
+ double val = quadrature->getWeight(iq) * c[iq];
+ val *= (*(phi->getPhi(i)))(quadrature->getLambda(iq)) * (*(phi->getPhi(0)))(quadrature->getLambda(iq));
+
+ (*mat)[j][i] += val;
- // calculate phi at QPs only once!
- for (int i = 0; i < nCol; i++) {
- phival[i] = (*(phi->getPhi(i)))(quadrature->getLambda(iq));
- }
+ val = quadrature->getWeight(iq) * c[iq];
+ val *= (*(phi->getPhi(i)))(quadrature->getLambda(iq)) * 0.5 * (*(phi->getPhi(1)))(quadrature->getLambda(iq));
- for (int i = 0; i < nRow; i++) {
- double psival = (*(psi->getPhi(i)))(psiQuadrature.getLambda(iq));
- for (int j = 0; j < nCol; j++) {
- (*mat)[i][j] += quadrature->getWeight(iq) * c[iq] * psival * phival[j];
+ (*mat)[j][i] += val;
+ } else {
+ double val = quadrature->getWeight(iq) * c[iq];
+ val *= (*(phi->getPhi(i)))(quadrature->getLambda(iq)) * 0.5 * (*(phi->getPhi(1)))(quadrature->getLambda(iq));
+
+ (*mat)[j][i] += val;
+ }
}
}
}
-
- DELETE scaledQuadrature;
- DELETE subElInfo;
-
- FREE_MEMORY(phival, double, nCol);
- FREE_MEMORY(c, double, nPoints);
-
- // ERROR_EXIT("SO, HIER GEHTS WEITER\n");
}
-
void StandardZOA::calculateElementVector(const ElInfo *elInfo, ElementVector *vec)
{
int nPoints = quadrature->getNumPoints();