diff --git a/AMDiS/src/Assembler.cc b/AMDiS/src/Assembler.cc
index 297ba661b12e244e4905173ee5a389848fa42ff0..adc1b6423b249b472d5e7d76074d39dd2d7fd41b 100644
--- a/AMDiS/src/Assembler.cc
+++ b/AMDiS/src/Assembler.cc
@@ -133,16 +133,13 @@ namespace AMDiS {
ElementMatrix m(nRow, nRow);
smallElInfo->getSubElemCoordsMat(m, rowFESpace->getBasisFcts()->getDegree());
-#if 0
- std::cout << "ASM MAT:" << std::endl;
- std::cout << mat << std::endl;
- std::cout << "MULT MAT:" << std::endl;
- std::cout << m << std::endl;
-#endif
-
ElementMatrix tmpMat(nRow, nRow);
- //tmpMat = m * mat;
- tmpMat = mat * m;
+
+ if (smallElInfo == colElInfo)
+ tmpMat = m * mat;
+ else
+ tmpMat = mat * trans(m);
+
mat = tmpMat;
}
diff --git a/AMDiS/src/ComponentTraverseInfo.cc b/AMDiS/src/ComponentTraverseInfo.cc
index a830f81a16efcccdf7d050cbe91a4a4e534c6636..772b39fc2fd7f81f2d0274061d74d111e34d5198 100644
--- a/AMDiS/src/ComponentTraverseInfo.cc
+++ b/AMDiS/src/ComponentTraverseInfo.cc
@@ -11,32 +11,32 @@ namespace AMDiS {
void SingleComponentInfo::updateStatus()
{
- if (rowFESpace == NULL) {
+ if (rowFeSpace == NULL) {
status = SingleComponentInfo::EMPTY;
return;
}
- if (colFESpace == NULL ||
- (colFESpace != NULL && rowFESpace->getMesh() == colFESpace->getMesh())) {
- if (auxFESpaces.size() == 0) {
+ if (colFeSpace == NULL ||
+ (colFeSpace != NULL && rowFeSpace->getMesh() == colFeSpace->getMesh())) {
+ if (auxFeSpaces.size() == 0) {
status = SingleComponentInfo::EQ_SPACES_NO_AUX;
} else {
status = SingleComponentInfo::EQ_SPACES_WITH_AUX;
- for (int i = 0; i < static_cast<int>(auxFESpaces.size()); i++) {
- if (auxFESpaces[i]->getMesh() != rowFESpace->getMesh()) {
+ for (int i = 0; i < static_cast<int>(auxFeSpaces.size()); i++) {
+ if (auxFeSpaces[i]->getMesh() != rowFeSpace->getMesh()) {
status = SingleComponentInfo::EQ_SPACES_WITH_DIF_AUX;
break;
}
}
}
} else {
- if (auxFESpaces.size() == 0) {
+ if (auxFeSpaces.size() == 0) {
status = SingleComponentInfo::DIF_SPACES_NO_AUX;
} else {
status = SingleComponentInfo::DIF_SPACES_WITH_AUX;
- for (int i = 0; i < static_cast<int>(auxFESpaces.size()); i++) {
- if (auxFESpaces[i]->getMesh() != rowFESpace->getMesh() &&
- auxFESpaces[i]->getMesh() != colFESpace->getMesh()) {
+ for (int i = 0; i < static_cast<int>(auxFeSpaces.size()); i++) {
+ if (auxFeSpaces[i]->getMesh() != rowFeSpace->getMesh() &&
+ auxFeSpaces[i]->getMesh() != colFeSpace->getMesh()) {
status = SingleComponentInfo::DIF_SPACES_WITH_DIF_AUX;
break;
}
@@ -45,4 +45,41 @@ namespace AMDiS {
}
}
+
+ const FiniteElemSpace* ComponentTraverseInfo::getRowFeSpace(int row)
+ {
+ FUNCNAME("ComponentTraverseInfo::getRowFeSpace()");
+
+ TEST_EXIT_DBG(row < nComponents)("No component traverse info for this row!\n");
+ TEST_EXIT_DBG(matrixComponents[row][row].getRowFeSpace() ==
+ matrixComponents[row][row].getColFeSpace())
+ ("Should not happen!\n");
+
+ return matrixComponents[row][row].getRowFeSpace();
+ }
+
+
+ const FiniteElemSpace* ComponentTraverseInfo::getNonRowFeSpace(int row)
+ {
+ FUNCNAME("ComponentTraverseInfo::getNonRowFeSpace()");
+
+ TEST_EXIT_DBG(row < nComponents)("No component traverse info for this row!\n");
+
+ const FiniteElemSpace* rowFeSpace = getRowFeSpace(row);
+
+ TEST_EXIT_DBG(rowFeSpace != NULL)("No row FE space!\n");
+
+ for (int i = 0; i < nComponents; i++) {
+ if (matrixComponents[row][i].getColFeSpace() != rowFeSpace)
+ return matrixComponents[row][i].getColFeSpace();
+ if (matrixComponents[row][i].getAuxFeSpace(0) != rowFeSpace)
+ return matrixComponents[row][i].getAuxFeSpace(0);
+ }
+
+ if (vectorComponents[row].getAuxFeSpace(0) != rowFeSpace)
+ return vectorComponents[row].getAuxFeSpace(0);
+
+ return NULL;
+ }
+
}
diff --git a/AMDiS/src/ComponentTraverseInfo.h b/AMDiS/src/ComponentTraverseInfo.h
index 4174599244e50029a90013d13d166c619872cb4d..0a7d8a27bf7c086319c9ae9a0c33a1d93db3e9b7 100644
--- a/AMDiS/src/ComponentTraverseInfo.h
+++ b/AMDiS/src/ComponentTraverseInfo.h
@@ -23,7 +23,7 @@
#define AMDIS_COMPONENTTRAVERSEINFO_H
#include <vector>
-
+#include "Global.h"
#include "FiniteElemSpace.h"
namespace AMDiS {
@@ -32,58 +32,58 @@ namespace AMDiS {
{
public:
SingleComponentInfo()
- : rowFESpace(NULL),
- colFESpace(NULL),
- auxFESpaces(0),
+ : rowFeSpace(NULL),
+ colFeSpace(NULL),
+ auxFeSpaces(0),
status(0)
{}
- void setFESpace(FiniteElemSpace *row, FiniteElemSpace *col = NULL)
+ void setFeSpace(FiniteElemSpace *row, FiniteElemSpace *col = NULL)
{
- rowFESpace = row;
- colFESpace = col;
+ rowFeSpace = row;
+ colFeSpace = col;
}
- void setAuxFESpaces(std::vector<const FiniteElemSpace*> feSpaces)
+ void setAuxFeSpaces(std::vector<const FiniteElemSpace*> feSpaces)
{
- auxFESpaces = feSpaces;
+ auxFeSpaces = feSpaces;
}
- void addAuxFESpace(const FiniteElemSpace *fe)
+ void addAuxFeSpace(const FiniteElemSpace *fe)
{
- auxFESpaces.push_back(fe);
+ auxFeSpaces.push_back(fe);
}
- bool hasFESpace()
+ bool hasFeSpace()
{
- return rowFESpace != NULL;
+ return rowFeSpace != NULL;
}
void updateStatus();
- int getNumAuxFESpaces()
+ int getNumAuxFeSpaces()
{
- return auxFESpaces.size();
+ return auxFeSpaces.size();
}
- FiniteElemSpace *getRowFESpace()
+ FiniteElemSpace *getRowFeSpace()
{
- return rowFESpace;
+ return rowFeSpace;
}
- FiniteElemSpace *getColFESpace()
+ FiniteElemSpace *getColFeSpace()
{
- return colFESpace;
+ return colFeSpace;
}
- const FiniteElemSpace *getAuxFESpace(int i)
+ const FiniteElemSpace *getAuxFeSpace(int i)
{
- return ((i < static_cast<int>(auxFESpaces.size())) ? auxFESpaces[i] : NULL);
+ return ((i < static_cast<int>(auxFeSpaces.size())) ? auxFeSpaces[i] : NULL);
}
- void setAuxFESpace(const FiniteElemSpace* fe, int pos)
+ void setAuxFeSpace(const FiniteElemSpace* fe, int pos)
{
- auxFESpaces[pos] = fe;
+ auxFeSpaces[pos] = fe;
}
int getStatus()
@@ -92,11 +92,11 @@ namespace AMDiS {
}
protected:
- FiniteElemSpace *rowFESpace;
+ FiniteElemSpace *rowFeSpace;
- FiniteElemSpace *colFESpace;
+ FiniteElemSpace *colFeSpace;
- std::vector<const FiniteElemSpace*> auxFESpaces;
+ std::vector<const FiniteElemSpace*> auxFeSpaces;
/// Status of the component, see the possible status flags below.
int status;
@@ -175,15 +175,36 @@ namespace AMDiS {
return vectorComponents[row].getStatus();
}
- const FiniteElemSpace* getAuxFESpace(int row, int col)
+ const FiniteElemSpace* getAuxFeSpace(int row, int col)
{
- return matrixComponents[row][col].getAuxFESpace(0);
+ return matrixComponents[row][col].getAuxFeSpace(0);
}
- const FiniteElemSpace* getAuxFESpace(int row)
+ const FiniteElemSpace* getAuxFeSpace(int row)
{
- return vectorComponents[row].getAuxFESpace(0);
+ return vectorComponents[row].getAuxFeSpace(0);
}
+
+ /** \brief
+ * Returns the row FE space for a given row number, i.e., the FE space
+ * of the diagonal matrix.
+ *
+ * \param[in] row Row number of the matrix line for which the FE space
+ * should be returned.
+ */
+ const FiniteElemSpace* getRowFeSpace(int row);
+
+
+ /** \brief
+ * Returns the non row FE space for a given row number. This is either the
+ * col FE space of an off diagonal matrix or the aux fe space of another
+ * matrix in the row or of the right hand side vector. If there is no non row
+ * FE space, this function returns a null pointer.
+ *
+ * \param[in] row Row number of the matrix line for which the non FE space
+ * should be returned.
+ */
+ const FiniteElemSpace* getNonRowFeSpace(int row);
protected:
int nComponents;
diff --git a/AMDiS/src/DOFMatrix.h b/AMDiS/src/DOFMatrix.h
index 11fd8c8ca070fdd94f75178f17d5171cc265335b..f297d46ca38a56eea4130dea7fbf2a267425be0e 100644
--- a/AMDiS/src/DOFMatrix.h
+++ b/AMDiS/src/DOFMatrix.h
@@ -346,17 +346,17 @@ namespace AMDiS {
bool symmetric();
- inline std::vector<Operator*> getOperators()
+ inline std::vector<Operator*>& getOperators()
{
return operators;
}
- inline std::vector<double*> getOperatorFactor()
+ inline std::vector<double*>& getOperatorFactor()
{
return operatorFactor;
}
- inline std::vector<double*> getOperatorEstFactor()
+ inline std::vector<double*>& getOperatorEstFactor()
{
return operatorEstFactor;
}
diff --git a/AMDiS/src/DOFVector.hh b/AMDiS/src/DOFVector.hh
index ef41f136f907885c1d9abb1825f21ab6011637cc..d97dede0e2fb5138877ea9e3ea87db160548121c 100644
--- a/AMDiS/src/DOFVector.hh
+++ b/AMDiS/src/DOFVector.hh
@@ -923,8 +923,9 @@ namespace AMDiS {
return result;
}
+
template<typename T>
- inline const DOFVector<T>& add(const DOFVector<T>& v1,
+ inline const DOFVector<T>& add(const DOFVector<T>& v1,
const DOFVector<T>& v2,
DOFVector<T>& result)
{
@@ -938,9 +939,15 @@ namespace AMDiS {
return result;
}
+
template<typename T>
const T *DOFVectorBase<T>::getLocalVector(const Element *el, T *d) const
{
+ FUNCNAME("DOFVectorBase<T>::getLocalVector()");
+
+ TEST_EXIT_DBG(feSpace->getMesh() == el->getMesh())
+ ("Element is defined on a different mesh than the DOF vector!\n");
+
static T* localVec = NULL;
static int localVecSize = 0;
const DOFAdmin* admin = feSpace->getAdmin();
@@ -956,6 +963,7 @@ namespace AMDiS {
delete [] localVec;
localVec = new T[nBasFcts];
}
+
if (!localVec)
localVec = new T[nBasFcts];
@@ -977,6 +985,7 @@ namespace AMDiS {
return result;
}
+
template<typename T>
const T *DOFVectorBase<T>::getVecAtQPs(const ElInfo *elInfo,
const Quadrature *quad,
diff --git a/AMDiS/src/DualTraverse.cc b/AMDiS/src/DualTraverse.cc
index 958112792b3d4459f83670cea89f242b83b67dc7..de4ccd81d86b2e012f4adc93ebe04e385d39d069 100644
--- a/AMDiS/src/DualTraverse.cc
+++ b/AMDiS/src/DualTraverse.cc
@@ -56,11 +56,12 @@ namespace AMDiS {
// call standard traverse
*elInfo1 = stack1.traverseFirst(mesh1, level1, flag1);
while (*elInfo1 != NULL && skipEl1(*elInfo1)) {
- *elInfo1 = stack1.traverseNext(*elInfo1);
+ *elInfo1 = stack1.traverseNext(*elInfo1);
}
+
*elInfo2 = stack2.traverseFirst(mesh2, level2, flag2);
while (*elInfo2 != NULL && skipEl2(*elInfo2)) {
- *elInfo2 = stack2.traverseNext(*elInfo2);
+ *elInfo2 = stack2.traverseNext(*elInfo2);
}
// finished ?
@@ -87,28 +88,29 @@ namespace AMDiS {
return true;
}
+
bool DualTraverse::traverseNext(ElInfo **elInfo1,
ElInfo **elInfo2,
ElInfo **elInfoSmall,
ElInfo **elInfoLarge)
{
// 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));
- }
+ 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 ?
- if (*elInfo1 == NULL || *elInfo2 == NULL) {
- TEST_EXIT(*elInfo1 == *elInfo2)("invalid dual traverse\n");
- return false;
- }
+ if (*elInfo1 == NULL || *elInfo2 == NULL) {
+ TEST_EXIT(*elInfo1 == *elInfo2)("invalid dual traverse\n");
+ return false;
+ }
// finished ?
if (*elInfo1 == NULL || *elInfo2 == NULL) {
@@ -132,6 +134,7 @@ namespace AMDiS {
return true;
}
+
void DualTraverse::prepareNextStep(ElInfo **elInfo1,
ElInfo **elInfo2,
ElInfo **elInfoSmall,
@@ -170,18 +173,18 @@ namespace AMDiS {
if (!fillSubElemMat)
return;
-// mtl::dense2D<double>& subCoordsMat =
-// const_cast<mtl::dense2D<double>&>(elInfoSmall->getSubElemCoordsMat());
-// mtl::dense2D<double>& subCoordsMat_so =
-// const_cast<mtl::dense2D<double>&>(elInfoSmall->getSubElemCoordsMat_so());
+ // mtl::dense2D<double>& subCoordsMat =
+ // const_cast<mtl::dense2D<double>&>(elInfoSmall->getSubElemCoordsMat());
+ // mtl::dense2D<double>& subCoordsMat_so =
+ // const_cast<mtl::dense2D<double>&>(elInfoSmall->getSubElemCoordsMat_so());
if (elInfo1 == elInfoSmall) {
-// stack1.getCoordsInElem(elInfo2, basisFcts, subCoordsMat);
-// stack1.getCoordsInElem_so(elInfo2, basisFcts, subCoordsMat_so);
+ // stack1.getCoordsInElem(elInfo2, basisFcts, subCoordsMat);
+ // stack1.getCoordsInElem_so(elInfo2, basisFcts, subCoordsMat_so);
stack1.fillRefinementPath(*elInfoSmall, *elInfo2);
} else {
-// stack2.getCoordsInElem(elInfo1, basisFcts, subCoordsMat);
-// stack2.getCoordsInElem_so(elInfo1, basisFcts, subCoordsMat_so);
+ // stack2.getCoordsInElem(elInfo1, basisFcts, subCoordsMat);
+ // stack2.getCoordsInElem_so(elInfo1, basisFcts, subCoordsMat_so);
stack2.fillRefinementPath(*elInfoSmall, *elInfo1);
}
}
diff --git a/AMDiS/src/DualTraverse.h b/AMDiS/src/DualTraverse.h
index bce6ccf6871abc2de92bf6ba015e66d21f8a0230..0c4059b2039bf2fbe7873325a4b58e9d65fb4c31 100644
--- a/AMDiS/src/DualTraverse.h
+++ b/AMDiS/src/DualTraverse.h
@@ -28,6 +28,18 @@
namespace AMDiS {
+ /** \brief
+ * Stores the four pointers to element info structures, that are required for the
+ * dual mesh traverse.
+ */
+ struct DualElInfo
+ {
+ ElInfo *rowElInfo;
+ ElInfo *colElInfo;
+ ElInfo *smallElInfo;
+ ElInfo *largeElInfo;
+ };
+
/// Parallel traversal of two meshes.
class DualTraverse
{
@@ -51,12 +63,34 @@ namespace AMDiS {
ElInfo **elInfoSmall,
ElInfo **elInfoLarge);
+ /// Alternative use for starting dual traversal.
+ inline bool traverseFirst(Mesh *mesh1, Mesh *mesh2,
+ int level1, int level2,
+ Flag flag1, Flag flag2,
+ DualElInfo &dualElInfo)
+ {
+ return traverseFirst(mesh1, mesh2, level1, level2, flag1, flag2,
+ &(dualElInfo.rowElInfo),
+ &(dualElInfo.colElInfo),
+ &(dualElInfo.smallElInfo),
+ &(dualElInfo.largeElInfo));
+ }
+
/// Get next ElInfo combination
bool traverseNext(ElInfo **elInfoNext1,
ElInfo **elInfoNext2,
ElInfo **elInfoSmall,
ElInfo **elInfoLarge);
+ /// Alternative use for getting the next elements in the dual traversal.
+ inline bool traverseNext(DualElInfo &dualElInfo)
+ {
+ return traverseNext(&(dualElInfo.rowElInfo),
+ &(dualElInfo.colElInfo),
+ &(dualElInfo.smallElInfo),
+ &(dualElInfo.largeElInfo));
+ }
+
bool check(ElInfo **elInfo1,
ElInfo **elInfo2,
ElInfo **elInfoSmall,
diff --git a/AMDiS/src/Estimator.cc b/AMDiS/src/Estimator.cc
index 5550158d397422adc7e1f4798f66c0093ef6ad76..c534c2227697107a1835d77e85b81a5e07b8bce9 100644
--- a/AMDiS/src/Estimator.cc
+++ b/AMDiS/src/Estimator.cc
@@ -1,32 +1,96 @@
#include "Estimator.h"
#include "Traverse.h"
#include "Parameters.h"
+#include "DualTraverse.h"
namespace AMDiS {
Estimator::Estimator(std::string name_, int r)
: name(name_),
norm(NO_NORM),
- row(r)
+ row(r),
+ mesh(NULL),
+ auxMesh(NULL),
+ traverseInfo(0)
{
+ FUNCNAME("Estimator::Estimator()");
+
GET_PARAMETER(0, name + "->error norm", "%d", &norm);
}
+
double Estimator::estimate(double ts)
{
FUNCNAME("Estimator::estimate()");
-
+
+ bool dualTraverse = false;
+ for (unsigned int i = 0; i < matrix.size(); i++) {
+ TEST_EXIT(traverseInfo.getStatus(row, i) != SingleComponentInfo::DIF_SPACES_WITH_DIF_AUX)
+ ("Not yet implemented!\n");
+
+ if (traverseInfo.getStatus(row, i) == SingleComponentInfo::EQ_SPACES_WITH_DIF_AUX ||
+ traverseInfo.getStatus(row, i) == SingleComponentInfo::DIF_SPACES_NO_AUX ||
+ traverseInfo.getStatus(row, i) == SingleComponentInfo::DIF_SPACES_WITH_AUX)
+ dualTraverse = true;
+ }
+
+ if (!dualTraverse) {
+ mesh = uh[row == -1 ? 0 : row]->getFESpace()->getMesh();
+ auxMesh = NULL;
+ } else {
+ const FiniteElemSpace *mainFeSpace = traverseInfo.getRowFeSpace(row);
+ const FiniteElemSpace *auxFeSpace = traverseInfo.getNonRowFeSpace(row);
+
+ TEST_EXIT(mainFeSpace)("No main FE space!\n");
+ TEST_EXIT(auxFeSpace)("No aux FE space!\n");
+
+ mesh = mainFeSpace->getMesh();
+ auxMesh = auxFeSpace->getMesh();
+
+ TEST_EXIT_DBG(mainFeSpace->getBasisFcts()->getDegree() ==
+ auxFeSpace->getBasisFcts()->getDegree())
+ ("Mh, do you really want to do this? Think about it ...\n");
+ }
+
init(ts);
+ if (!dualTraverse)
+ singleMeshTraverse();
+ else
+ dualMeshTraverse();
+
+ exit();
+
+ return est_sum;
+ }
+
+
+ void Estimator::singleMeshTraverse()
+ {
+ FUNCNAME("Estimator::singleMeshTraverse()");
+
TraverseStack stack;
ElInfo *elInfo = stack.traverseFirst(mesh, -1, traverseFlag);
while (elInfo) {
estimateElement(elInfo);
elInfo = stack.traverseNext(elInfo);
}
+ }
- exit();
- return est_sum;
+ void Estimator::dualMeshTraverse()
+ {
+ FUNCNAME("Estimator::dualMeshTraverse()");
+
+ DualTraverse dualTraverse;
+ DualElInfo dualElInfo;
+
+ bool cont = dualTraverse.traverseFirst(mesh, auxMesh, -1, -1,
+ traverseFlag, traverseFlag,
+ dualElInfo);
+ while (cont) {
+ estimateElement(dualElInfo.rowElInfo, &dualElInfo);
+ cont = dualTraverse.traverseNext(dualElInfo);
+ }
}
}
diff --git a/AMDiS/src/Estimator.h b/AMDiS/src/Estimator.h
index 1a0e9cedec23852cfe929ebffb745667e7f61480..e8a4e937e5e8f437017ef97260ad297f968d7543 100644
--- a/AMDiS/src/Estimator.h
+++ b/AMDiS/src/Estimator.h
@@ -32,6 +32,8 @@
#include "FixVec.h"
#include "SystemVector.h"
#include "CreatorInterface.h"
+#include "ComponentTraverseInfo.h"
+#include "DualTraverse.h"
namespace AMDiS {
@@ -44,7 +46,9 @@ namespace AMDiS {
class Estimator
{
public:
- Estimator() {}
+ Estimator()
+ : traverseInfo(0)
+ {}
/// Constructor.
Estimator(std::string name_, int r);
@@ -64,11 +68,17 @@ namespace AMDiS {
///
virtual void init(double timestep) {}
- ///
- virtual void estimateElement(ElInfo *elInfo) {}
+ /** \brief
+ * Estimates the error on an element. If there is more than one mesh used in the
+ * problem description, it may be necessary to used the dual mesh traverse. In this
+ * case elInfo is the element of the main mesh, i.e., the mesh of the row FE space,
+ * and dualElInfo contains all elInfo informations about the main mesh element and
+ * the col (or aux) mesh element.
+ */
+ virtual void estimateElement(ElInfo *elInfo, DualElInfo *dualElInfo = NULL) {}
///
- virtual void exit(bool output=true) {}
+ virtual void exit(bool output = true) {}
/// Returns \ref est_sum of the Estimator
inline double getErrorSum() const
@@ -161,6 +171,19 @@ namespace AMDiS {
return row;
}
+ /// Sets \ref traverseInfo.
+ void setTraverseInfo(const ComponentTraverseInfo &ti)
+ {
+ traverseInfo = ti;
+ }
+
+ protected:
+ /// Traverse one mesh to estimate the error.
+ void singleMeshTraverse();
+
+ /// Traverses two meshes to estimate the error.
+ void dualMeshTraverse();
+
protected:
/// Name of the Estimator
std::string name;
@@ -209,9 +232,26 @@ namespace AMDiS {
Flag traverseFlag;
+ /** \brief
+ * The mesh on which the error must be estimated. If there is more than one mesh
+ * used, here the main, i.e., the row mesh, is stored.
+ */
Mesh *mesh;
+ /** \brief
+ * If there is only one mesh used at all, this variable is not used. In the case
+ * that the error must be estimated on a system row with more than one mesh, here
+ * either the column mesh or the auxiliary mesh is stored.
+ */
+ Mesh *auxMesh;
+
double timestep;
+
+ /** \brief
+ * Stores information about which mesh(es) must be travesed to estimate
+ * the error on the component matrices.
+ */
+ ComponentTraverseInfo traverseInfo;
};
diff --git a/AMDiS/src/Mesh.cc b/AMDiS/src/Mesh.cc
index 059afb148ea2b8795aac3671f7ebef0071694c0a..9e421e01de9fab37477d0559583c61236275496b 100644
--- a/AMDiS/src/Mesh.cc
+++ b/AMDiS/src/Mesh.cc
@@ -631,6 +631,8 @@ namespace AMDiS {
ElInfo* Mesh::createNewElInfo()
{
+ FUNCNAME("Mesh::createNewElInfo()");
+
switch(dim) {
case 1:
return new ElInfo1d(this);
@@ -647,6 +649,7 @@ namespace AMDiS {
}
}
+
bool Mesh::findElInfoAtPoint(const WorldVector<double>& xy,
ElInfo *el_info,
DimVec<double>& bary,
diff --git a/AMDiS/src/Operator.h b/AMDiS/src/Operator.h
index f64dfa51fde6774004bf7db3d701b4151d37b59b..0c3f4cfcac368d402b1017715a8c2c7f6ab3b258 100644
--- a/AMDiS/src/Operator.h
+++ b/AMDiS/src/Operator.h
@@ -76,7 +76,7 @@ namespace AMDiS {
{}
/// Returs \auxFESpaces, the list of all aux fe spaces the operator makes use off.
- std::vector<const FiniteElemSpace*>& getAuxFESpaces()
+ std::vector<const FiniteElemSpace*>& getAuxFeSpaces()
{
return auxFESpaces;
}
@@ -3237,7 +3237,7 @@ namespace AMDiS {
}
/// Returns \ref auxFESpaces.
- inline std::vector<const FiniteElemSpace*> getAuxFESpaces()
+ inline std::vector<const FiniteElemSpace*> getAuxFeSpaces()
{
return auxFESpaces;
}
diff --git a/AMDiS/src/Operator.hh b/AMDiS/src/Operator.hh
index 0752563f7b7ac142681605f1ad584eda73deab0f..c1b36002be752e22acc5bd80f1a557037bfa63f9 100644
--- a/AMDiS/src/Operator.hh
+++ b/AMDiS/src/Operator.hh
@@ -6,8 +6,8 @@ namespace AMDiS {
secondOrder[0].push_back(term);
term->operat = this;
auxFESpaces.insert(auxFESpaces.end(),
- term->getAuxFESpaces().begin(),
- term->getAuxFESpaces().end());
+ term->getAuxFeSpaces().begin(),
+ term->getAuxFeSpaces().end());
for (int i = 1; i < omp_get_overall_max_threads(); i++) {
T *newTerm = new T(static_cast<const T>(*term));
@@ -26,8 +26,8 @@ namespace AMDiS {
}
term->operat = this;
auxFESpaces.insert(auxFESpaces.end(),
- term->getAuxFESpaces().begin(),
- term->getAuxFESpaces().end());
+ term->getAuxFeSpaces().begin(),
+ term->getAuxFeSpaces().end());
for (int i = 1; i < omp_get_overall_max_threads(); i++) {
T *newTerm = new T(static_cast<const T>(*term));
@@ -45,8 +45,8 @@ namespace AMDiS {
zeroOrder[0].push_back(term);
term->operat = this;
auxFESpaces.insert(auxFESpaces.end(),
- term->getAuxFESpaces().begin(),
- term->getAuxFESpaces().end());
+ term->getAuxFeSpaces().begin(),
+ term->getAuxFeSpaces().end());
for (int i = 1; i < omp_get_overall_max_threads(); i++) {
T *newTerm = new T(static_cast<const T>(*term));
diff --git a/AMDiS/src/ProblemVec.cc b/AMDiS/src/ProblemVec.cc
index f4335ee01aa57b291031a5a8481d47068b8c72ab..5668230f6145569e497299e1198a48f0325ab41b 100644
--- a/AMDiS/src/ProblemVec.cc
+++ b/AMDiS/src/ProblemVec.cc
@@ -306,9 +306,9 @@ namespace AMDiS {
for (int i = 0; i < nComponents; i++) {
for (int j = 0; j < nComponents; j++)
- traverseInfo.getMatrix(i, j).setFESpace(componentSpaces[i], componentSpaces[j]);
+ traverseInfo.getMatrix(i, j).setFeSpace(componentSpaces[i], componentSpaces[j]);
- traverseInfo.getVector(i).setFESpace(componentSpaces[i]);
+ traverseInfo.getVector(i).setFeSpace(componentSpaces[i]);
}
// create dof admin for vertex dofs if neccessary
@@ -531,7 +531,10 @@ namespace AMDiS {
Estimator *scalEstimator = estimator[i];
if (scalEstimator) {
+ traverseInfo.updateStatus();
+ scalEstimator->setTraverseInfo(traverseInfo);
scalEstimator->estimate(adaptInfo->getTimestep());
+
adaptInfo->setEstSum(scalEstimator->getErrorSum(), i);
adaptInfo->setEstMax(scalEstimator->getErrorMax(), i);
adaptInfo->setTimeEstSum(scalEstimator->getTimeEst(), i);
@@ -659,7 +662,7 @@ namespace AMDiS {
for (int j = 0; j < nComponents; j++) {
- std::cout << "-------" << i << " " << j << "----------------" << std::endl;
+ // std::cout << "-------" << i << " " << j << "----------------" << std::endl;
// Only if this variable is true, the current matrix will be assembled.
bool assembleMatrix = true;
@@ -709,8 +712,6 @@ namespace AMDiS {
// The simplest case: either the right hand side has no operaters, no aux
// fe spaces, or all aux fe spaces are equal to the row and col fe space.
-
- // std::cout << "ASM 1" << std::endl;
assembleOnOneMesh(componentSpaces[i],
assembleFlag,
assembleMatrix ? matrix : NULL,
@@ -720,9 +721,6 @@ namespace AMDiS {
// Row fe space and col fe space are both equal, but right hand side has at
// least one another aux fe space.
-
- // std::cout << "ASM 2" << std::endl;
-
if (assembleMatrix)
assembleOnOneMesh(componentSpaces[i],
assembleFlag,
@@ -731,7 +729,7 @@ namespace AMDiS {
if (i == j && asmVector)
assembleOnDifMeshes2(componentSpaces[i],
- traverseInfo.getAuxFESpace(i),
+ traverseInfo.getAuxFeSpace(i),
assembleFlag,
NULL,
rhs->getDOFVector(i));
@@ -740,21 +738,18 @@ namespace AMDiS {
}
} else if (traverseInfo.getStatus(i, j) == SingleComponentInfo::EQ_SPACES_WITH_DIF_AUX) {
- // std::cout << "ASM 4" << std::endl;
assembleOnOneMesh(componentSpaces[i],
assembleFlag,
assembleMatrix ? matrix : NULL,
((i == j) && asmVector) ? rhs->getDOFVector(i) : NULL);
- // std::cout << "ASM 5" << std::endl;
assembleOnDifMeshes2(componentSpaces[i],
- traverseInfo.getAuxFESpace(i, j),
+ traverseInfo.getAuxFeSpace(i, j),
assembleFlag,
assembleMatrix ? matrix : NULL,
((i == j) && asmVector) ? rhs->getDOFVector(i) : NULL);
} else if (traverseInfo.getStatus(i, j) == SingleComponentInfo::DIF_SPACES_NO_AUX ||
traverseInfo.getStatus(i, j) == SingleComponentInfo::DIF_SPACES_WITH_AUX) {
- // std::cout << "ASM 6" << std::endl;
assembleOnDifMeshes(componentSpaces[i], componentSpaces[j],
assembleFlag,
assembleMatrix ? matrix : NULL,
@@ -782,8 +777,6 @@ namespace AMDiS {
assembleFlag);
}
- // exit(0);
-
solverMatrix.setMatrix(*systemMatrix);
createPrecon();
@@ -884,11 +877,11 @@ namespace AMDiS {
(*systemMatrix)[i][j]->addOperator(op, factor, estFactor);
- traverseInfo.getMatrix(i, j).setAuxFESpaces(op->getAuxFESpaces());
+ traverseInfo.getMatrix(i, j).setAuxFeSpaces(op->getAuxFeSpaces());
- for (int k = 0; k < static_cast<int>(op->getAuxFESpaces().size()); k++) {
- if ((op->getAuxFESpaces())[k]->getMesh() != componentSpaces[i]->getMesh() ||
- (op->getAuxFESpaces())[k]->getMesh() != componentSpaces[j]->getMesh()) {
+ for (int k = 0; k < static_cast<int>(op->getAuxFeSpaces().size()); k++) {
+ if ((op->getAuxFeSpaces())[k]->getMesh() != componentSpaces[i]->getMesh() ||
+ (op->getAuxFeSpaces())[k]->getMesh() != componentSpaces[j]->getMesh()) {
op->setNeedDualTraverse(true);
break;
}
@@ -917,10 +910,10 @@ namespace AMDiS {
rhs->getDOFVector(i)->addOperator(op, factor, estFactor);
- traverseInfo.getVector(i).setAuxFESpaces(op->getAuxFESpaces());
+ traverseInfo.getVector(i).setAuxFeSpaces(op->getAuxFeSpaces());
- for (int j = 0; j < static_cast<int>(op->getAuxFESpaces().size()); j++) {
- if ((op->getAuxFESpaces())[j]->getMesh() != componentSpaces[i]->getMesh()) {
+ for (int j = 0; j < static_cast<int>(op->getAuxFeSpaces().size()); j++) {
+ if ((op->getAuxFeSpaces())[j]->getMesh() != componentSpaces[i]->getMesh()) {
op->setNeedDualTraverse(true);
break;
}
@@ -1196,10 +1189,6 @@ namespace AMDiS {
&rowElInfo, &colElInfo,
&smallElInfo, &largeElInfo);
while (cont) {
-// std::cout << "ROW = " << rowElInfo->getElement()->getIndex() << " "
-// << "COL = " << colElInfo->getElement()->getIndex() << " "
-// << "SMA = " << smallElInfo->getElement()->getIndex() << " "
-// << "LAR = " << largeElInfo->getElement()->getIndex() << std::endl;
if (useGetBound)
basisFcts->getBound(rowElInfo, bound);
@@ -1396,6 +1385,7 @@ namespace AMDiS {
solution->serialize(out);
}
+
void ProblemVec::deserialize(std::istream &in)
{
FUNCNAME("ProblemVec::deserialize()");
diff --git a/AMDiS/src/Recovery.cc b/AMDiS/src/Recovery.cc
index ccff6ec1268333c391e36db37e219bee67e65762..f9a669033456438b51e0a4c0d4cb87dca2c3062c 100644
--- a/AMDiS/src/Recovery.cc
+++ b/AMDiS/src/Recovery.cc
@@ -1,12 +1,10 @@
#include "Recovery.h"
-RecoveryStructure& RecoveryStructure::operator=(const RecoveryStructure& rhs) {
-
-
+RecoveryStructure& RecoveryStructure::operator=(const RecoveryStructure& rhs)
+{
if (rhs.coords) {
- if (!coords) {
- coords = new WorldVector<double>;
- }
+ if (!coords)
+ coords = new WorldVector<double>;
*coords = *rhs.coords;
} else {
if (coords) {
@@ -33,7 +31,7 @@ RecoveryStructure& RecoveryStructure::operator=(const RecoveryStructure& rhs) {
} else {
if (rec_uh) {
delete rec_uh;
- rec_uh=NULL;
+ rec_uh = NULL;
}
}
@@ -60,211 +58,183 @@ RecoveryStructure& RecoveryStructure::operator=(const RecoveryStructure& rhs) {
}
return *this;
-};
-
+}
void RecoveryStructure::print()
{
- FUNCNAME("RecoveryStructure::print");
-
- int i, j;
+ FUNCNAME("RecoveryStructure::print()");
std::cout << std::endl;
MSG("Coordinates of the node: ");
std::cout << *coords << std::endl;
- if (A)
- {
- MSG("Interior vertex: printing system information.\n\n");
-
- int n = A->getNumRows();
-
- MSG("System matrix:\n");
- for (i=0; i<n; i++)
- {
- MSG("( ");
- for (j=0; j<i; j++)
- std::cout << "* ";
- for (j=i; j<n; j++)
- std::cout << (*A)[i][j] << " ";
- std::cout << ")" << std::endl;
- }
-
- MSG("Right hand side:\n");
- for (i=0; i<n; i++)
- {
- MSG("( ");
- if (rec_grdUh)
- std::cout << (*rec_grdUh)[i];
- else
- std::cout << (*rec_uh)[i];
- std::cout << " )" << std::endl;
- }
-
- if (neighbors)
- {
- MSG("Printing neighbors vertices\n\n");
-
- MSG("Number of neighbors: ");
- std::cout << neighbors->size() << std::endl << std::endl;
-
- MSG("List of neighbors: ");
- std::set<DegreeOfFreedom>::const_iterator setIterator;
-
- for (setIterator = neighbors->begin();
- setIterator != neighbors->end();
- ++setIterator)
- std::cout << " " << *setIterator;
+ if (A) {
+ MSG("Interior vertex: printing system information.\n\n");
+
+ int n = A->getNumRows();
+
+ MSG("System matrix:\n");
+ for (int i = 0; i < n; i++) {
+ MSG("( ");
+ for (int j = 0; j < i; j++)
+ std::cout << "* ";
+ for (int j = i; j < n; j++)
+ std::cout << (*A)[i][j] << " ";
+ std::cout << ")" << std::endl;
+ }
- std::cout << std::endl << std::endl;
- }
+ MSG("Right hand side:\n");
+ for (int i = 0; i < n; i++) {
+ MSG("( ");
+ if (rec_grdUh)
+ std::cout << (*rec_grdUh)[i];
+ else
+ std::cout << (*rec_uh)[i];
+ std::cout << " )" << std::endl;
}
- else
- {
- MSG("Boundary vertex or not a vertex node: printing interior neighbors\n\n");
+
+ if (neighbors) {
+ MSG("Printing neighbors vertices\n\n");
MSG("Number of neighbors: ");
std::cout << neighbors->size() << std::endl << std::endl;
-
+
MSG("List of neighbors: ");
std::set<DegreeOfFreedom>::const_iterator setIterator;
-
- for (setIterator = neighbors->begin();
- setIterator != neighbors->end();
+
+ for (setIterator = neighbors->begin(); setIterator != neighbors->end();
++setIterator)
std::cout << " " << *setIterator;
-
+
std::cout << std::endl << std::endl;
}
+ } else {
+ MSG("Boundary vertex or not a vertex node: printing interior neighbors\n\n");
- WAIT;
+ MSG("Number of neighbors: ");
+ std::cout << neighbors->size() << std::endl << std::endl;
+
+ MSG("List of neighbors: ");
+ std::set<DegreeOfFreedom>::const_iterator setIterator;
+
+ for (setIterator = neighbors->begin(); setIterator != neighbors->end();
+ ++setIterator)
+ std::cout << " " << *setIterator;
- return;
+ std::cout << std::endl << std::endl;
+ }
+
+ WAIT;
}
-/*****************************************************************************/
void Recovery::set_feSpace(const FiniteElemSpace *fe_space)
{
- FUNCNAME("Recovery::set_feSpace");
-
- if (!feSpace || (feSpace != fe_space))
- {
- if (struct_vec)
- {
- delete struct_vec;
- struct_vec = NULL;
- }
-
- feSpace = fe_space;
+ FUNCNAME("Recovery::set_feSpace()");
- // create new structure vector
- struct_vec = new DOFVector<RecoveryStructure>(feSpace, "struct vec");
+ if (!feSpace || feSpace != fe_space) {
+ if (struct_vec) {
+ delete struct_vec;
+ struct_vec = NULL;
}
- return;
+ feSpace = fe_space;
+
+ // create new structure vector
+ struct_vec = new DOFVector<RecoveryStructure>(feSpace, "struct vec");
+ }
}
+
int Recovery::set_exponents(int degree)
{
- FUNCNAME("Recovery::set_exponents");
+ FUNCNAME("Recovery::set_exponents()");
int dow = Global::getGeo(WORLD);
int number_monomials = degree + 1;
// Computing number of monomials.
- if (dow > 1)
- {
- number_monomials *= (degree + 2);
- number_monomials /= 2;
- }
+ if (dow > 1) {
+ number_monomials *= (degree + 2);
+ number_monomials /= 2;
+ }
- if (dow == 3)
- {
- number_monomials *= (degree + 3);
- number_monomials /= 3;
- }
+ if (dow == 3) {
+ number_monomials *= (degree + 3);
+ number_monomials /= 3;
+ }
// Allocating memory.
- if (n_monomials != number_monomials)
- {
- n_monomials = number_monomials;
- exponents.resize(n_monomials);
-
- if (matrix_fcts)
- matrix_fcts->resize(n_monomials, n_monomials);
- else
- matrix_fcts = new Matrix<Monomial*>(n_monomials, n_monomials);
- }
+ if (n_monomials != number_monomials) {
+ n_monomials = number_monomials;
+ exponents.resize(n_monomials);
+
+ if (matrix_fcts)
+ matrix_fcts->resize(n_monomials, n_monomials);
+ else
+ matrix_fcts = new Matrix<Monomial*>(n_monomials, n_monomials);
+ }
// Setting vector of exponents.
- int i, j, k, count = 0;
-
- switch (dow)
- {
- case 1: // 1D monomials.
- for (i=0; i<=degree; i++)
- exponents[count++][0] = i;
- break;
-
- case 2: // 2D monomials.
- for (i=0; i<=degree; i++)
- for (j=0; j<=i; j++)
- {
- exponents[count][0] = i - j;
- exponents[count++][1] = j;
- }
- break;
-
- case 3: // 3D monomials.
- for (i=0; i<=degree; i++)
- for (j=0; j<=i; j++)
- for (k=0; k<=j; k++)
- {
- exponents[count][0] = i - j;
- exponents[count][1] = j - k;
- exponents[count++][2] = k;
- }
- break;
-
- default:
- ERROR_EXIT("Which dimension have your world???\n");
- }
+ int count = 0;
+
+ switch (dow) {
+ case 1: // 1D monomials.
+ for (int i = 0; i <= degree; i++)
+ exponents[count++][0] = i;
+ break;
+
+ case 2: // 2D monomials.
+ for (int i = 0; i <= degree; i++)
+ for (int j = 0; j <= i; j++) {
+ exponents[count][0] = i - j;
+ exponents[count++][1] = j;
+ }
+ break;
+
+ case 3: // 3D monomials.
+ for (int i = 0; i <= degree; i++)
+ for (int j = 0; j <= i; j++)
+ for (int k = 0; k <= j; k++) {
+ exponents[count][0] = i - j;
+ exponents[count][1] = j - k;
+ exponents[count++][2] = k;
+ }
+ break;
- TEST_EXIT_DBG(count == n_monomials)("There must be an error!\n");
+ default:
+ ERROR_EXIT("Which dimension have your world???\n");
+ }
+ TEST_EXIT_DBG(count == n_monomials)("There must be an error!\n");
+
// Setting matrix of monomials.
WorldVector<int> sum;
- for (i=0; i<n_monomials; i++)
- for (j=i; j<n_monomials; j++)
- {
- // Computing exponent vector of monomial.
- sum = exponents[i] + exponents[j];
- (*matrix_fcts)[i][j] = new Monomial(sum);
- }
+ for (int i = 0; i < n_monomials; i++)
+ for (int j = i; j < n_monomials; j++) {
+ // Computing exponent vector of monomial.
+ sum = exponents[i] + exponents[j];
+ (*matrix_fcts)[i][j] = new Monomial(sum);
+ }
return n_monomials;
}
+
void Recovery::compute_integrals(DOFVector<double> *uh, ElInfo *elInfo,
RecoveryStructure *rec_struct,
AbstractFunction<double, WorldVector<double> > *f_vec,
AbstractFunction<double, double> *f_scal,
DOFVector<double> *aux_vec)
{
- FUNCNAME("Recovery::compute_integrals");
+ FUNCNAME("Recovery::compute_integrals()");
TEST_EXIT_DBG(!(f_vec && f_scal))("Only one diffusion function, please!\n");
- int i, j, k;
- double sum;
WorldVector<double> vec_sum;
-
- Quadrature *quad;
- int n_points;
WorldVector<double> quad_pts; // For world coordinates of quadrature points.
int deg_f = 0;
@@ -278,278 +248,238 @@ void Recovery::compute_integrals(DOFVector<double> *uh, ElInfo *elInfo,
else
deg_f += feSpace->getBasisFcts()->getDegree();
- for (i=0; i<n_monomials; i++)
- {
- // Computing contributions to system matrix.
- for (j=i; j<n_monomials; j++)
- {
- sum = 0.0;
- quad = Quadrature::provideQuadrature(Global::getGeo(WORLD),
- (*matrix_fcts)[i][j]->getDegree());
- n_points = quad->getNumPoints();
-
- for (k=0; k<n_points; k++)
- {
- elInfo->coordToWorld(quad->getLambda(k), quad_pts);
- sum += quad->getWeight(k) *
- (*(*matrix_fcts)[i][j])(quad_pts, *rec_struct->coords);
- }
- (*(rec_struct->A))[i][j] += sum * elInfo->getDet();
- }
+ for (int i = 0; i < n_monomials; i++) {
+ // Computing contributions to system matrix.
+ for (int j = i; j < n_monomials; j++) {
+ double sum = 0.0;
+ Quadrature *quad = Quadrature::provideQuadrature(Global::getGeo(WORLD),
+ (*matrix_fcts)[i][j]->getDegree());
+ int n_points = quad->getNumPoints();
+
+ for (int k = 0; k < n_points; k++) {
+ elInfo->coordToWorld(quad->getLambda(k), quad_pts);
+ sum += quad->getWeight(k) *
+ (*(*matrix_fcts)[i][j])(quad_pts, *rec_struct->coords);
+ }
+ (*(rec_struct->A))[i][j] += sum * elInfo->getDet();
+ }
- quad = Quadrature::
- provideQuadrature(Global::getGeo(WORLD),
+ Quadrature *quad = Quadrature::
+ provideQuadrature(Global::getGeo(WORLD),
(*matrix_fcts)[0][i]->getDegree() + deg_f);
- n_points = quad->getNumPoints();
-
- double *uhAtQP = new double[n_points];
-
- // Computing contributions to right hand side.
- if (gradient) // For gradient recovery.
- {
- double fAtQP = 1.0;
- if (f_scal)
- {
- if (aux_vec)
- aux_vec->getVecAtQPs(elInfo, quad, NULL, uhAtQP);
- else
- uh->getVecAtQPs(elInfo, quad, NULL, uhAtQP);
- }
-
- // Get gradient at quadrature points
- WorldVector<double> *grdAtQP = new WorldVector<double>[n_points];
- uh->getGrdAtQPs(elInfo, quad, NULL, grdAtQP);
- vec_sum = 0.0;
- for (k=0; k<n_points; k++)
- {
- elInfo->coordToWorld(quad->getLambda(k), quad_pts);
- if (f_vec)
- fAtQP = (*f_vec)(quad_pts);
- if (f_scal)
- fAtQP = (*f_scal)(uhAtQP[k]);
-
- vec_sum = vec_sum + grdAtQP[k] * fAtQP * quad->getWeight(k)
- * (*(*matrix_fcts)[0][i])(quad_pts, *rec_struct->coords);
- }
- (*rec_struct->rec_grdUh)[i] = (*rec_struct->rec_grdUh)[i]
- + vec_sum * elInfo->getDet();
-
- delete [] grdAtQP;
- }
- else // For recovery of DOFVector.
- {
- // Get uh at quadrature points
+ int n_points = quad->getNumPoints();
+ double *uhAtQP = new double[n_points];
+
+ // Computing contributions to right hand side.
+ if (gradient) { // For gradient recovery.
+ double fAtQP = 1.0;
+ if (f_scal) {
+ if (aux_vec)
+ aux_vec->getVecAtQPs(elInfo, quad, NULL, uhAtQP);
+ else
uh->getVecAtQPs(elInfo, quad, NULL, uhAtQP);
- sum = 0.0;
- for (k=0; k<n_points; k++)
- {
- elInfo->coordToWorld(quad->getLambda(k), quad_pts);
- sum += uhAtQP[k] * quad->getWeight(k)
- * (*(*matrix_fcts)[0][i])(quad_pts, *rec_struct->coords);
- }
- (*rec_struct->rec_uh)[i] += sum * elInfo->getDet();
- }
- delete [] uhAtQP;
+ }
+
+ // Get gradient at quadrature points
+ WorldVector<double> *grdAtQP = new WorldVector<double>[n_points];
+ uh->getGrdAtQPs(elInfo, quad, NULL, grdAtQP);
+ vec_sum = 0.0;
+ for (int k = 0; k < n_points; k++) {
+ elInfo->coordToWorld(quad->getLambda(k), quad_pts);
+ if (f_vec)
+ fAtQP = (*f_vec)(quad_pts);
+ if (f_scal)
+ fAtQP = (*f_scal)(uhAtQP[k]);
+
+ vec_sum = vec_sum + grdAtQP[k] * fAtQP * quad->getWeight(k)
+ * (*(*matrix_fcts)[0][i])(quad_pts, *rec_struct->coords);
+ }
+ (*rec_struct->rec_grdUh)[i] = (*rec_struct->rec_grdUh)[i]
+ + vec_sum * elInfo->getDet();
+
+ delete [] grdAtQP;
+ } else { // For recovery of DOFVector.
+ // Get uh at quadrature points
+ uh->getVecAtQPs(elInfo, quad, NULL, uhAtQP);
+ double sum = 0.0;
+ for (int k = 0; k < n_points; k++) {
+ elInfo->coordToWorld(quad->getLambda(k), quad_pts);
+ sum += uhAtQP[k] * quad->getWeight(k)
+ * (*(*matrix_fcts)[0][i])(quad_pts, *rec_struct->coords);
+ }
+ (*rec_struct->rec_uh)[i] += sum * elInfo->getDet();
}
+ delete [] uhAtQP;
+ }
}
+
void Recovery::compute_interior_sums(DOFVector<double> *uh, ElInfo *elInfo,
RecoveryStructure *rec_struct, Quadrature *quad,
AbstractFunction<double, WorldVector<double> > *f_vec,
AbstractFunction<double, double> *f_scal,
DOFVector<double> *aux_vec)
{
- FUNCNAME("Recovery::compute_sums");
+ FUNCNAME("Recovery::compute_sums()");
TEST_EXIT_DBG(gradient)("SPR of solution need computing node sums.\n");
-
TEST_EXIT_DBG(!(f_vec && f_scal))("Only one diffusion function, please!\n");
- int i, j, k;
- double sum;
WorldVector<double> vec_sum;
-
int n_points = quad->getNumPoints();
WorldVector<double> quad_pts; // For world coordinates of quadrature points.
-
double *uhAtQP = new double[n_points];
WorldVector<double> *grdAtQP = new WorldVector<double>[n_points];
- for (i=0; i<n_monomials; i++)
- {
- // Computing contributions to system matrix.
- for (j=i; j<n_monomials; j++)
- {
- sum = 0.0;
- for (k=0; k<n_points; k++)
- {
- elInfo->coordToWorld(quad->getLambda(k), quad_pts);
- sum += (*(*matrix_fcts)[i][j])(quad_pts, *rec_struct->coords);
- }
- (*(rec_struct->A))[i][j] += sum;
- }
+ for (int i = 0; i < n_monomials; i++) {
+ // Computing contributions to system matrix.
+ for (int j = i; j < n_monomials; j++) {
+ double sum = 0.0;
+ for (int k = 0; k < n_points; k++) {
+ elInfo->coordToWorld(quad->getLambda(k), quad_pts);
+ sum += (*(*matrix_fcts)[i][j])(quad_pts, *rec_struct->coords);
+ }
+ (*(rec_struct->A))[i][j] += sum;
+ }
- // Computing contributions to right hand side.
- double fAtQP = 1.0;
+ // Computing contributions to right hand side.
+ double fAtQP = 1.0;
+ if (f_scal) {
+ if (aux_vec)
+ aux_vec->getVecAtQPs(elInfo, quad, NULL, uhAtQP);
+ else
+ uh->getVecAtQPs(elInfo, quad, NULL, uhAtQP);
+ }
+
+ // Get gradient at quadrature points
+ uh->getGrdAtQPs(elInfo, quad, NULL, grdAtQP);
+ vec_sum = 0.0;
+ for (int k = 0; k < n_points; k++) {
+ elInfo->coordToWorld(quad->getLambda(k), quad_pts);
+ if (f_vec)
+ fAtQP = (*f_vec)(quad_pts);
if (f_scal)
- {
- if (aux_vec)
- aux_vec->getVecAtQPs(elInfo, quad, NULL, uhAtQP);
- else
- uh->getVecAtQPs(elInfo, quad, NULL, uhAtQP);
- }
+ fAtQP = (*f_scal)(uhAtQP[k]);
- // Get gradient at quadrature points
- uh->getGrdAtQPs(elInfo, quad, NULL, grdAtQP);
- vec_sum = 0.0;
- for (k=0; k<n_points; k++)
- {
- elInfo->coordToWorld(quad->getLambda(k), quad_pts);
- if (f_vec)
- fAtQP = (*f_vec)(quad_pts);
- if (f_scal)
- fAtQP = (*f_scal)(uhAtQP[k]);
-
- vec_sum = vec_sum + grdAtQP[k] * fAtQP
- * (*(*matrix_fcts)[0][i])(quad_pts, *rec_struct->coords);
- }
- (*rec_struct->rec_grdUh)[i] = (*rec_struct->rec_grdUh)[i] + vec_sum;
+ vec_sum = vec_sum + grdAtQP[k] * fAtQP
+ * (*(*matrix_fcts)[0][i])(quad_pts, *rec_struct->coords);
}
+ (*rec_struct->rec_grdUh)[i] = (*rec_struct->rec_grdUh)[i] + vec_sum;
+ }
delete [] uhAtQP;
delete [] grdAtQP;
}
+
void Recovery::compute_node_sums(DOFVector<double> *uh, ElInfo *elInfo,
RecoveryStructure *rec_struct, DimVec<int> preDOFs,
int n_vertices, int n_edges, int n_faces)
{
- FUNCNAME("Recovery::compute_sums");
+ FUNCNAME("Recovery::compute_sums()");
TEST_EXIT_DBG(!gradient)
("SPR of flux or gradient need computing interior sums\n");
-
TEST_EXIT_DBG(feSpace->getMesh()->getDim()==1)
("At the moment only for linear finite elements.\n");
- int i, j, l;
-
WorldVector<double> node; // For world coordinates at nodes.
-
const DegreeOfFreedom **dof = elInfo->getElement()->getDOF();
const double *uh_loc = uh->getLocalVector(elInfo->getElement(), NULL);
- for (l=0; l<n_vertices; l++)
- {
- // Computing contributions of vertex nodes
- if (rec_struct->neighbors->insert(dof[l][preDOFs[VERTEX]]).second)
- {
- node = elInfo->getCoord(l);
- for (i=0; i<n_monomials; i++)
- {
- // Computing contributions to system matrix.
- for (j=i; j<n_monomials; j++)
- (*(rec_struct->A))[i][j] += (*(*matrix_fcts)[i][j])(node,
- *rec_struct->coords);
-
- // Computing contributions to right hand side.
- (*rec_struct->rec_uh)[i] += uh_loc[l]
- * (*(*matrix_fcts)[0][i])(node, *rec_struct->coords);
- }
- }
+ for (int l = 0; l < n_vertices; l++) {
+ // Computing contributions of vertex nodes
+ if (rec_struct->neighbors->insert(dof[l][preDOFs[VERTEX]]).second) {
+ node = elInfo->getCoord(l);
+ for (int i = 0; i < n_monomials; i++) {
+ // Computing contributions to system matrix.
+ for (int j = i; j < n_monomials; j++)
+ (*(rec_struct->A))[i][j] += (*(*matrix_fcts)[i][j])(node,
+ *rec_struct->coords);
+
+ // Computing contributions to right hand side.
+ (*rec_struct->rec_uh)[i] += uh_loc[l]
+ * (*(*matrix_fcts)[0][i])(node, *rec_struct->coords);
+ }
}
-
- return;
+ }
}
+
void Recovery::compute_sums_linear(DOFVector<double> *uh, ElInfo *elInfo,
RecoveryStructure *rec_struct,
int vertex, DimVec<int> preDOFs,
int n_vertices)
{
- FUNCNAME("Recovery::compute_sums_linear");
+ FUNCNAME("Recovery::compute_sums_linear()");
TEST_EXIT_DBG(!gradient)
("SPR of flux or gradient need computing interior sums\n");
- int i, j, l;
- DegreeOfFreedom k;
-
WorldVector<double> node; // For world coordinates at nodes.
-
const DegreeOfFreedom **dof = elInfo->getElement()->getDOF();
const double *uh_loc = uh->getLocalVector(elInfo->getElement(), NULL);
- for (l=0; l<n_vertices; l++)
- {
- // Computing contributions of vertex nodes
- k = dof[l][preDOFs[VERTEX]];
- if (rec_struct->neighbors->insert(k).second)
- {
- node = elInfo->getCoord(l);
- for (i=0; i<n_monomials; i++)
- {
- // Computing contributions to system matrix.
- for (j=i; j<n_monomials; j++)
- (*(rec_struct->A))[i][j] += (*(*matrix_fcts)[i][j])(node,
- *rec_struct->coords);
-
- // Computing contributions to right hand side.
- (*rec_struct->rec_uh)[i] += uh_loc[l]
- * (*(*matrix_fcts)[0][i])(node, *rec_struct->coords);
- }
- }
+ for (int l = 0; l < n_vertices; l++) {
+ // Computing contributions of vertex nodes
+ DegreeOfFreedom k = dof[l][preDOFs[VERTEX]];
+ if (rec_struct->neighbors->insert(k).second) {
+ node = elInfo->getCoord(l);
+ for (int i = 0; i < n_monomials; i++) {
+ // Computing contributions to system matrix.
+ for (int j = i; j < n_monomials; j++)
+ (*(rec_struct->A))[i][j] += (*(*matrix_fcts)[i][j])(node,
+ *rec_struct->coords);
+
+ // Computing contributions to right hand side.
+ (*rec_struct->rec_uh)[i] += uh_loc[l]
+ * (*(*matrix_fcts)[0][i])(node, *rec_struct->coords);
+ }
}
+ }
- if (vertex>1)
- if (elInfo->getNeighbour(vertex))
- {
- int oppVertex = elInfo->getOppVertex(vertex);
- k = elInfo->getNeighbour(vertex)->getDOF(oppVertex)[preDOFs[VERTEX]];
-
- if (rec_struct->neighbors->insert(k).second)
- {
- node = elInfo->getOppCoord(vertex);
- for (i=0; i<n_monomials; i++)
- {
- // Computing contributions to system matrix.
- for (j=i; j<n_monomials; j++)
- (*(rec_struct->A))[i][j] += (*(*matrix_fcts)[i][j])(node,
- *rec_struct->coords);
-
- // Computing contributions to right hand side.
- (*rec_struct->rec_uh)[i] += (*uh)[k]
- * (*(*matrix_fcts)[0][i])(node, *rec_struct->coords);
- }
- }
+ if (vertex > 1 && elInfo->getNeighbour(vertex)) {
+ int oppVertex = elInfo->getOppVertex(vertex);
+ DegreeOfFreedom k = elInfo->getNeighbour(vertex)->getDOF(oppVertex)[preDOFs[VERTEX]];
+
+ if (rec_struct->neighbors->insert(k).second) {
+ node = elInfo->getOppCoord(vertex);
+ for (int i = 0; i < n_monomials; i++) {
+ // Computing contributions to system matrix.
+ for (int j = i; j < n_monomials; j++)
+ (*(rec_struct->A))[i][j] += (*(*matrix_fcts)[i][j])(node, *rec_struct->coords);
+
+ // Computing contributions to right hand side.
+ (*rec_struct->rec_uh)[i] += (*uh)[k] * (*(*matrix_fcts)[0][i])(node, *rec_struct->coords);
}
-
- return;
+ }
+ }
}
+
void Recovery::fill_struct_vec(DOFVector<double> *uh,
AbstractFunction<double, WorldVector<double> > *f_vec,
AbstractFunction<double, double> *f_scal,
DOFVector<double> *aux_vec)
{
- FUNCNAME("Recovery::fill_struct_vec");
+ FUNCNAME("Recovery::fill_struct_vec()");
// Information on the mesh.
Mesh *mesh = feSpace->getMesh();
- int dim = mesh->getDim();
+ int dim = mesh->getDim();
// Geometric information.
int n_vertices = Global::getGeo(VERTEX, dim);
- int n_edges = Global::getGeo(EDGE, dim);
- int n_faces = Global::getGeo(FACE, dim);
+ int n_edges = Global::getGeo(EDGE, dim);
+ int n_faces = Global::getGeo(FACE, dim);
// Information concerning the finite element space.
const BasisFunction *basis_fcts = feSpace->getBasisFcts();
- DimVec<int> *nDOFs = basis_fcts->getNumberOfDOFs();
+ DimVec<int> *nDOFs = basis_fcts->getNumberOfDOFs();
// Information from DOFAdmin.
const DOFAdmin *admin = feSpace->getAdmin();
- DimVec<int> preDOFs = admin->getNumberOfPreDOFs();
+ DimVec<int> preDOFs = admin->getNumberOfPreDOFs();
// Variables for storing temporary information.
DimVec<DegreeOfFreedom> interior_vertices(dim);
@@ -564,201 +494,171 @@ void Recovery::fill_struct_vec(DOFVector<double> *uh,
DimVec<int> pre_dofs(dim, NO_INIT);
if (!gradient)
pre_dofs = uh->getFESpace()->getAdmin()->getNumberOfPreDOFs();
-
- // Auxiliar variables.
- int i, j, l, m, n_neighbors, n_count;
- DegreeOfFreedom k;
-
+
// Variables for traversing the mesh.
- TraverseStack stack;
- ElInfo *el_info;
- Flag fill_flag;
-
- fill_flag =
+ Flag fill_flag =
Mesh::CALL_LEAF_EL |
Mesh::FILL_COORDS |
Mesh::FILL_BOUND |
Mesh::FILL_DET |
Mesh::FILL_GRD_LAMBDA;
- if (degree==2 && dim>1)
+ if (degree == 2 && dim > 1)
fill_flag |= Mesh::FILL_NEIGH | Mesh::FILL_OPP_COORDS;
- el_info = stack.traverseFirst(mesh, -1, fill_flag);
+ TraverseStack stack;
+ ElInfo *el_info = stack.traverseFirst(mesh, -1, fill_flag);
- while (el_info) // traversing the mesh.
- {
- const DegreeOfFreedom **dof = el_info->getElement()->getDOF();
+ while (el_info) { // traversing the mesh.
+ const DegreeOfFreedom **dof = el_info->getElement()->getDOF();
+
+ int n_neighbors = 0; // counting interior vertices of element
+ for (int i = 0; i < n_vertices; i++) {
+ DegreeOfFreedom k = dof[i][preDOFs[VERTEX]];
+ if (el_info->getBoundary(VERTEX, i) == INTERIOR)
+ interior_vertices[n_neighbors++] = k;
+ }
+
+ TEST_EXIT_DBG(n_neighbors)
+ ("Each element should have a least one interior vertex!\n");
+
+ for (int i = 0; i < n_vertices; i++) { // Handling nodes on vertices.
+ DegreeOfFreedom k = dof[i][preDOFs[VERTEX]];
+
+ // Setting world coordinates of node.
+ if (!(*struct_vec)[k].coords) {
+ (*struct_vec)[k].coords = new WorldVector<double>;
+ *(*struct_vec)[k].coords = el_info->getCoord(i);
+ }
- n_neighbors = 0; // counting interior vertices of element
- for (i=0; i<n_vertices; i++)
- {
- k = dof[i][preDOFs[VERTEX]];
- if (el_info->getBoundary(VERTEX, i) == INTERIOR)
- interior_vertices[n_neighbors++] = k;
+ if (el_info->getBoundary(VERTEX, i) == INTERIOR) {
+ // Allocating memory for matrix and right hand side.
+ if (!(*struct_vec)[k].A) {
+ (*struct_vec)[k].A = new Matrix<double>(n_monomials, n_monomials);
+ *(*struct_vec)[k].A = 0.0;
+
+ if (gradient) {
+ (*struct_vec)[k].rec_grdUh = new Vector<WorldVector<double> >(n_monomials);
+ for (int j = 0; j < n_monomials; j++)
+ (*(*struct_vec)[k].rec_grdUh)[j] = 0.0;
+ } else {
+ (*struct_vec)[k].rec_uh = new Vector<double>(n_monomials);
+ *(*struct_vec)[k].rec_uh = 0.0;
+ }
}
- TEST_EXIT_DBG(n_neighbors)
- ("Each element should have a least one interior vertex!\n");
-
- for (i=0; i<n_vertices; i++) // Handling nodes on vertices.
- {
- k = dof[i][preDOFs[VERTEX]];
-
- // Setting world coordinates of node.
- if (!(*struct_vec)[k].coords)
- {
- (*struct_vec)[k].coords = new WorldVector<double>;
- *(*struct_vec)[k].coords = el_info->getCoord(i);
- }
-
- if (el_info->getBoundary(VERTEX, i) == INTERIOR)
- {
- // Allocating memory for matrix and right hand side.
- if (!(*struct_vec)[k].A)
- {
- (*struct_vec)[k].A = new Matrix<double>(n_monomials, n_monomials);
- *(*struct_vec)[k].A = 0.0;
-
- if (gradient)
- {
- (*struct_vec)[k].rec_grdUh
- = new Vector<WorldVector<double> >(n_monomials);
- for (j=0; j<n_monomials; j++)
- (*(*struct_vec)[k].rec_grdUh)[j] = 0.0;
- }
- else
- {
- (*struct_vec)[k].rec_uh = new Vector<double>(n_monomials);
- *(*struct_vec)[k].rec_uh = 0.0;
- }
- }
-
- // Computing the integrals.
- if (method)
- compute_integrals(uh, el_info, &(*struct_vec)[k],
- f_vec, f_scal, aux_vec);
- else if (gradient)
- compute_interior_sums(uh, el_info, &(*struct_vec)[k], quad,
- f_vec, f_scal, aux_vec);
- else
- {
- if (!(*struct_vec)[k].neighbors)
- (*struct_vec)[k].neighbors = new std::set<DegreeOfFreedom>;
-
- if (degree==2 && dim>1)
- compute_sums_linear(uh, el_info, &(*struct_vec)[k],
- i, pre_dofs, n_vertices);
- else
- compute_node_sums(uh, el_info, &(*struct_vec)[k],
- pre_dofs, n_vertices, n_edges, n_faces);
- }
- }
- else // Setting list of adjacent interior vertices.
- {
- if (!(*struct_vec)[k].neighbors)
- (*struct_vec)[k].neighbors = new std::set<DegreeOfFreedom>;
-
- for (j=0; j<n_neighbors; j++)
- (*struct_vec)[k].neighbors->insert(interior_vertices[j]);
- }
+ // Computing the integrals.
+ if (method)
+ compute_integrals(uh, el_info, &(*struct_vec)[k],
+ f_vec, f_scal, aux_vec);
+ else if (gradient)
+ compute_interior_sums(uh, el_info, &(*struct_vec)[k], quad,
+ f_vec, f_scal, aux_vec);
+ else {
+ if (!(*struct_vec)[k].neighbors)
+ (*struct_vec)[k].neighbors = new std::set<DegreeOfFreedom>;
+
+ if (degree == 2 && dim > 1)
+ compute_sums_linear(uh, el_info, &(*struct_vec)[k],
+ i, pre_dofs, n_vertices);
+ else
+ compute_node_sums(uh, el_info, &(*struct_vec)[k],
+ pre_dofs, n_vertices, n_edges, n_faces);
}
+ } else { // Setting list of adjacent interior vertices.
+ if (!(*struct_vec)[k].neighbors)
+ (*struct_vec)[k].neighbors = new std::set<DegreeOfFreedom>;
+
+ for (int j = 0; j < n_neighbors; j++)
+ (*struct_vec)[k].neighbors->insert(interior_vertices[j]);
+ }
+ }
- n_count = n_vertices;
-
- if (dim > 1) // Handling nodes on edges.
- for (i=0; i<n_edges; i++)
- for (j=0; j<(*nDOFs)[EDGE]; j++)
- {
- k = dof[n_vertices+i][preDOFs[EDGE]+j];
-
- if (!(*struct_vec)[k].coords)
- {
- // Setting world coordinates of node.
- el_info->coordToWorld(*basis_fcts->getCoords(n_count),
- coordinates);
- (*struct_vec)[k].coords = new WorldVector<double>;
- *(*struct_vec)[k].coords = coordinates;
-
- // Setting list of adjacent interior vertices.
- (*struct_vec)[k].neighbors = new std::set<DegreeOfFreedom>;
-
- if (el_info->getBoundary(EDGE, i) == INTERIOR)
- for (m=0; m<2; m++)
- {
- l = Global::getReferenceElement(dim)->getVertexOfEdge(i, m);
- if (el_info->getBoundary(VERTEX, l) == INTERIOR)
- (*struct_vec)[k].neighbors->insert(dof[l][preDOFs[VERTEX]]);
- }
- else
- for (m=0; m<n_neighbors; m++)
- (*struct_vec)[k].neighbors->insert(interior_vertices[m]);
- }
-
- n_count++;
- }
-
- if (dim == 3) // Handling nodes on faces.
- for (i=0; i<n_faces; i++)
- for (j=0; j<(*nDOFs)[FACE]; j++)
- {
- k = dof[n_vertices+n_edges+i][preDOFs[FACE]+j];
-
- if (!(*struct_vec)[k].coords)
- {
- // Setting world coordinates of node.
- el_info->coordToWorld(*basis_fcts->getCoords(n_count),
- coordinates);
- (*struct_vec)[k].coords = new WorldVector<double>;
- *(*struct_vec)[k].coords = coordinates;
-
- // Setting list of adjacent interior vertices.
- (*struct_vec)[k].neighbors = new std::set<DegreeOfFreedom>;
-
- if (el_info->getBoundary(FACE, i) == INTERIOR)
- for (m=0; m<3; m++)
- {
- l = Global::getReferenceElement(dim)
- ->getVertexOfPosition(FACE, i, m);
- if (el_info->getBoundary(VERTEX, l) == INTERIOR)
- (*struct_vec)[k].neighbors->insert(dof[l][preDOFs[VERTEX]]);
- }
- else
- for (m=0; m<n_neighbors; m++)
- (*struct_vec)[k].neighbors->insert(interior_vertices[m]);
- }
-
- n_count++;
- }
-
- if ((*nDOFs)[CENTER]) // Handling nodes on center of element.
- for (j=0; j<(*nDOFs)[CENTER]; j++)
- {
- k = dof[n_vertices+n_edges+n_faces][preDOFs[CENTER]+j];
+ int n_count = n_vertices;
+ if (dim > 1) // Handling nodes on edges.
+ for (int i = 0; i < n_edges; i++)
+ for (int j = 0; j < (*nDOFs)[EDGE]; j++) {
+ DegreeOfFreedom k = dof[n_vertices+i][preDOFs[EDGE] + j];
+
+ if (!(*struct_vec)[k].coords) {
// Setting world coordinates of node.
- el_info->coordToWorld(*basis_fcts->getCoords(n_count), coordinates);
- (*struct_vec)[k].coords = new WorldVector<double>;
+ el_info->coordToWorld(*basis_fcts->getCoords(n_count),
+ coordinates);
+ (*struct_vec)[k].coords = new WorldVector<double>;
*(*struct_vec)[k].coords = coordinates;
-
+
// Setting list of adjacent interior vertices.
(*struct_vec)[k].neighbors = new std::set<DegreeOfFreedom>;
-
- for (m=0; m<n_neighbors; m++)
- (*struct_vec)[k].neighbors->insert(interior_vertices[m]);
-
- n_count++;
+
+ if (el_info->getBoundary(EDGE, i) == INTERIOR)
+ for (int m = 0; m < 2; m++) {
+ int l = Global::getReferenceElement(dim)->getVertexOfEdge(i, m);
+ if (el_info->getBoundary(VERTEX, l) == INTERIOR)
+ (*struct_vec)[k].neighbors->insert(dof[l][preDOFs[VERTEX]]);
+ } else
+ for (int m = 0; m < n_neighbors; m++)
+ (*struct_vec)[k].neighbors->insert(interior_vertices[m]);
}
+
+ n_count++;
+ }
- el_info = stack.traverseNext(el_info);
- }
-
- return;
+ if (dim == 3) // Handling nodes on faces.
+ for (int i = 0; i < n_faces; i++)
+ for (int j = 0; j < (*nDOFs)[FACE]; j++) {
+ DegreeOfFreedom k = dof[n_vertices+n_edges+i][preDOFs[FACE] + j];
+
+ if (!(*struct_vec)[k].coords) {
+ // Setting world coordinates of node.
+ el_info->coordToWorld(*basis_fcts->getCoords(n_count),
+ coordinates);
+ (*struct_vec)[k].coords = new WorldVector<double>;
+ *(*struct_vec)[k].coords = coordinates;
+
+ // Setting list of adjacent interior vertices.
+ (*struct_vec)[k].neighbors = new std::set<DegreeOfFreedom>;
+
+ if (el_info->getBoundary(FACE, i) == INTERIOR)
+ for (int m = 0; m < 3; m++) {
+ int l = Global::getReferenceElement(dim)->getVertexOfPosition(FACE, i, m);
+ if (el_info->getBoundary(VERTEX, l) == INTERIOR)
+ (*struct_vec)[k].neighbors->insert(dof[l][preDOFs[VERTEX]]);
+ }
+ else
+ for (int m = 0; m < n_neighbors; m++)
+ (*struct_vec)[k].neighbors->insert(interior_vertices[m]);
+ }
+
+ n_count++;
+ }
+
+ if ((*nDOFs)[CENTER]) // Handling nodes on center of element.
+ for (int j = 0; j < (*nDOFs)[CENTER]; j++) {
+ DegreeOfFreedom k = dof[n_vertices+n_edges+n_faces][preDOFs[CENTER] + j];
+
+ // Setting world coordinates of node.
+ el_info->coordToWorld(*basis_fcts->getCoords(n_count), coordinates);
+ (*struct_vec)[k].coords = new WorldVector<double>;
+ *(*struct_vec)[k].coords = coordinates;
+
+ // Setting list of adjacent interior vertices.
+ (*struct_vec)[k].neighbors = new std::set<DegreeOfFreedom>;
+
+ for (int m = 0; m < n_neighbors; m++)
+ (*struct_vec)[k].neighbors->insert(interior_vertices[m]);
+
+ n_count++;
+ }
+
+ el_info = stack.traverseNext(el_info);
+ }
}
+
void Recovery::recoveryUh(DOFVector<double> *uh, DOFVector<double> &rec_vec)
{
- FUNCNAME("Recovery::recoveryUh");
+ FUNCNAME("Recovery::recoveryUh()");
clear();
@@ -772,55 +672,46 @@ void Recovery::recoveryUh(DOFVector<double> *uh, DOFVector<double> &rec_vec)
DOFVector<RecoveryStructure>::Iterator SV_it(struct_vec, USED_DOFS);
// Solving local systems.
- for (SV_it.reset(); !SV_it.end(); ++SV_it)
- {
- if ((*SV_it).A)
- {
- TEST(Cholesky::solve((*SV_it).A, (*SV_it).rec_uh, (*SV_it).rec_uh))
- ("There must be an error, matrix is not positive definite.\n");
- }
+ for (SV_it.reset(); !SV_it.end(); ++SV_it) {
+ if ((*SV_it).A) {
+ TEST(Cholesky::solve((*SV_it).A, (*SV_it).rec_uh, (*SV_it).rec_uh))
+ ("There must be an error, matrix is not positive definite.\n");
}
+ }
// define result vector
- DOFVector<double> *result = &rec_vec;
+ DOFVector<double> *result = &rec_vec;
result->set(0.0);
DOFVector<double>::Iterator result_it(result, USED_DOFS);
std::set<DegreeOfFreedom>::const_iterator setIterator;
- int i;
- for (SV_it.reset(), result_it.reset(); !result_it.end();
- ++SV_it, ++result_it)
- {
- if ((*SV_it).rec_uh)
- *result_it = (*(*SV_it).rec_uh)[0];
- else
- {
- if ((*SV_it).neighbors) {
- for (setIterator = (*SV_it).neighbors->begin();
- setIterator != (*SV_it).neighbors->end();
- ++setIterator)
- {
- for (i=0; i<n_monomials; i++)
- *result_it = *result_it + (*(*struct_vec)[*setIterator].rec_uh)[i] *
- (*(*matrix_fcts)[0][i])(*(*SV_it).coords,
- *(*struct_vec)[*setIterator].coords);
- }
- *result_it /= (*SV_it).neighbors->size();
- }
- else
- *result_it=0.0;
+ for (SV_it.reset(), result_it.reset(); !result_it.end(); ++SV_it, ++result_it) {
+ if ((*SV_it).rec_uh) {
+ *result_it = (*(*SV_it).rec_uh)[0];
+ } else {
+ if ((*SV_it).neighbors) {
+ for (setIterator = (*SV_it).neighbors->begin();
+ setIterator != (*SV_it).neighbors->end();
+ ++setIterator) {
+ for (int i = 0; i < n_monomials; i++)
+ *result_it = *result_it + (*(*struct_vec)[*setIterator].rec_uh)[i] *
+ (*(*matrix_fcts)[0][i])(*(*SV_it).coords,
+ *(*struct_vec)[*setIterator].coords);
}
+ *result_it /= (*SV_it).neighbors->size();
+ } else
+ *result_it=0.0;
}
-
- return;
+ }
}
+
DOFVector<double>*
Recovery::recoveryUh(DOFVector<double> *uh, const FiniteElemSpace *fe_space)
{
- FUNCNAME("Recovery::recoveryUh");
+ FUNCNAME("Recovery::recoveryUh()");
clear();
@@ -833,25 +724,23 @@ Recovery::recoveryUh(DOFVector<double> *uh, const FiniteElemSpace *fe_space)
DOFVector<RecoveryStructure>::Iterator SV_it(struct_vec, USED_DOFS);
// Solving local systems.
- for (SV_it.reset(); !SV_it.end(); ++SV_it)
- {
- if ((*SV_it).A)
- {
- TEST(Cholesky::solve((*SV_it).A, (*SV_it).rec_uh, (*SV_it).rec_uh))
- ("There must be an error, matrix is not positive definite.\n");
- }
+ for (SV_it.reset(); !SV_it.end(); ++SV_it) {
+ if ((*SV_it).A) {
+ TEST(Cholesky::solve((*SV_it).A, (*SV_it).rec_uh, (*SV_it).rec_uh))
+ ("There must be an error, matrix is not positive definite.\n");
}
-
+ }
+
// define result vector
static DOFVector<double> *vec = NULL;
- DOFVector<double> *result = NULL;
+ DOFVector<double> *result = NULL;
// Allocate memory for result vector
- if (vec && vec->getFESpace() != feSpace)
- {
- delete vec;
- vec = NULL;
- }
+ if (vec && vec->getFESpace() != feSpace) {
+ delete vec;
+ vec = NULL;
+ }
+
if (!vec)
vec = new DOFVector<double>(feSpace, "gradient");
result = vec;
@@ -860,35 +749,31 @@ Recovery::recoveryUh(DOFVector<double> *uh, const FiniteElemSpace *fe_space)
DOFVector<double>::Iterator result_it(result, USED_DOFS);
std::set<DegreeOfFreedom>::const_iterator setIterator;
- int i;
for (SV_it.reset(), result_it.reset(); !result_it.end();
- ++SV_it, ++result_it)
- {
- if ((*SV_it).rec_uh)
- *result_it = (*(*SV_it).rec_uh)[0];
- else
- {
- if ((*SV_it).neighbors) {
- for (setIterator = (*SV_it).neighbors->begin();
- setIterator != (*SV_it).neighbors->end();
- ++setIterator)
- {
- for (i=0; i<n_monomials; i++)
- *result_it = *result_it + (*(*struct_vec)[*setIterator].rec_uh)[i] *
- (*(*matrix_fcts)[0][i])(*(*SV_it).coords,
- *(*struct_vec)[*setIterator].coords);
- }
- *result_it /= (*SV_it).neighbors->size();
- }
- else
- *result_it=0.0;
+ ++SV_it, ++result_it) {
+ if ((*SV_it).rec_uh) {
+ *result_it = (*(*SV_it).rec_uh)[0];
+ } else {
+ if ((*SV_it).neighbors) {
+ for (setIterator = (*SV_it).neighbors->begin();
+ setIterator != (*SV_it).neighbors->end();
+ ++setIterator) {
+ for (int i = 0; i < n_monomials; i++)
+ *result_it = *result_it + (*(*struct_vec)[*setIterator].rec_uh)[i] *
+ (*(*matrix_fcts)[0][i])(*(*SV_it).coords,
+ *(*struct_vec)[*setIterator].coords);
}
+ *result_it /= (*SV_it).neighbors->size();
+ } else
+ *result_it=0.0;
}
+ }
return result;
}
+
DOFVector<WorldVector<double> >*
Recovery::recovery(DOFVector<double> *uh, const FiniteElemSpace *fe_space,
AbstractFunction<double, WorldVector<double> > *f_vec,
@@ -958,6 +843,7 @@ Recovery::recovery(DOFVector<double> *uh, const FiniteElemSpace *fe_space,
return result;
}
+
DOFVector<WorldVector<double> >*
Recovery::recovery(DOFVector<double> *uh,
AbstractFunction<double, WorldVector<double> > *f_vec,
@@ -979,11 +865,10 @@ Recovery::recovery(DOFVector<double> *uh,
delete vec;
vec = NULL;
}
- if (!vec) {
- vec = new DOFVector<WorldVector<double> >(fe_space, "gradient");
- }
- result = vec;
+ if (!vec)
+ vec = new DOFVector<WorldVector<double> >(fe_space, "gradient");
+ result = vec;
result->set(WorldVector<double>(DEFAULT_VALUE, 0.0));
DOFVector<double> volume(fe_space, "volume");
@@ -991,21 +876,16 @@ Recovery::recovery(DOFVector<double> *uh,
Mesh *mesh = fe_space->getMesh();
int dim = mesh->getDim();
-
const BasisFunction *basFcts = fe_space->getBasisFcts();
DOFAdmin *admin = fe_space->getAdmin();
-
int numPreDOFs = admin->getNumberOfPreDOFs(0);
-
DimVec<double> bary(dim, DEFAULT_VALUE, (1.0 / (dim + 1.0)));
WorldVector<double> barycenter; // For world coordinates at barycenter
// traverse mesh
TraverseStack stack;
-
Flag fillFlag =
Mesh::CALL_LEAF_EL | Mesh::FILL_DET | Mesh::FILL_GRD_LAMBDA | Mesh::FILL_COORDS;
-
ElInfo *elInfo = stack.traverseFirst(mesh, -1, fillFlag);
while (elInfo) {
@@ -1071,6 +951,4 @@ void Recovery::test(DOFVector<double> *uh, const FiniteElemSpace *fe_space)
std::cout << position << std::endl;
(*struct_vec)[position].print();
}
-
- return;
}
diff --git a/AMDiS/src/Recovery.h b/AMDiS/src/Recovery.h
index 498189d75c6daf2bd5ecd6475174a5181667a10c..7fa0a0feec5dcfd997a40b22d6a8d4639d646d17 100644
--- a/AMDiS/src/Recovery.h
+++ b/AMDiS/src/Recovery.h
@@ -22,250 +22,237 @@
#ifndef AMDIS_RECOVERY_H
#define AMDIS_RECOVERY_H
-#include "set"
+#include <set>
#include "Lagrange.h"
#include "Traverse.h"
#include "DOFVector.h"
#include "Cholesky.h"
-using namespace std;
-using namespace AMDiS;
+namespace AMDiS {
-class Monomial : public
-BinaryAbstractFunction<double, WorldVector<double>, WorldVector<double> >
-{
-public:
-
- // Constructor.
- Monomial(WorldVector<int> expon)
- : BinaryAbstractFunction<double, WorldVector<double>, WorldVector<double> >(),
- exponent(expon)
- {
- degree_ = exponent[0];
- for (int i = 1; i<exponent.size(); i++)
- degree_ += exponent[i];
- }
-
- virtual ~Monomial() {};
-
- double operator()(const WorldVector<double>& y,
- const WorldVector<double>& z) const
+ class Monomial : public
+ BinaryAbstractFunction<double, WorldVector<double>, WorldVector<double> >
{
- double result = pow(y[0] - z[0], exponent[0]);
- for (int i = 1; i < exponent.size(); i++)
- result *= pow(y[i] - z[i], exponent[i]);
-
- return result;
- }
-
-private:
- WorldVector<int> exponent;
-};
-
-
-class RecoveryStructure
-{
-public: // construtor, destructor
- // constructor
- RecoveryStructure() :
- coords(NULL),
- A(NULL),
- rec_uh(NULL),
- rec_grdUh(NULL),
- neighbors(NULL)
- {}
-
- // Destructor?
- ~RecoveryStructure()
- {
- clear();
- }
-
- RecoveryStructure(const RecoveryStructure& rhs)
- : coords(NULL),
- A(NULL),
- rec_uh(NULL),
- rec_grdUh(NULL),
- neighbors(NULL)
- {
- *this=rhs;
- }
-
- RecoveryStructure& operator=(const RecoveryStructure& rhs);
-
-
- // Clear recovery structure
- inline void clear()
+ public:
+ Monomial(WorldVector<int> expon)
+ : BinaryAbstractFunction<double, WorldVector<double>, WorldVector<double> >(),
+ exponent(expon)
+ {
+ degree_ = exponent[0];
+ for (int i = 1; i<exponent.size(); i++)
+ degree_ += exponent[i];
+ }
+
+ virtual ~Monomial() {}
+
+ double operator()(const WorldVector<double>& y,
+ const WorldVector<double>& z) const
+ {
+ double result = pow(y[0] - z[0], exponent[0]);
+ for (int i = 1; i < exponent.size(); i++)
+ result *= pow(y[i] - z[i], exponent[i]);
+
+ return result;
+ }
+
+ private:
+ WorldVector<int> exponent;
+ };
+
+
+ class RecoveryStructure
{
- if (coords)
- {
+ public:
+ RecoveryStructure()
+ : coords(NULL),
+ A(NULL),
+ rec_uh(NULL),
+ rec_grdUh(NULL),
+ neighbors(NULL)
+ {}
+
+ ~RecoveryStructure()
+ {
+ clear();
+ }
+
+ RecoveryStructure(const RecoveryStructure& rhs)
+ : coords(NULL),
+ A(NULL),
+ rec_uh(NULL),
+ rec_grdUh(NULL),
+ neighbors(NULL)
+ {
+ *this = rhs;
+ }
+
+ RecoveryStructure& operator=(const RecoveryStructure& rhs);
+
+ /// Clear recovery structure
+ inline void clear()
+ {
+ if (coords) {
delete coords;
- coords = NULL;
+ coords = NULL;
}
-
- if (A)
- {
+
+ if (A) {
delete A;
- A = NULL;
+ A = NULL;
}
-
- if (rec_uh)
- {
+
+ if (rec_uh) {
delete rec_uh;
- rec_uh = NULL;
+ rec_uh = NULL;
}
-
- if (rec_grdUh)
- {
+
+ if (rec_grdUh) {
delete rec_grdUh;
rec_grdUh = NULL;
}
- if (neighbors)
- {
+ if (neighbors) {
delete neighbors;
neighbors = NULL;
}
+ }
+
+ /// Prints recovery structure (for test purposes)
+ void print();
+
+
+ private:
+ /// World coordinates of the node.
+ WorldVector<double> *coords;
+
+ /// For interior nodes.
+ Matrix<double> *A;
+ Vector<double> *rec_uh;
+ Vector<WorldVector<double> > *rec_grdUh;
+
+ /// For boundary, edge, face, of center nodes: interior neighbors nodes.
+ /// For interior vertices in uh-recovery: neighbors nodes.
+ std::set<DegreeOfFreedom> *neighbors;
+
+ friend class Recovery;
};
-
- // Prints recovery structure (for test purposes)
- void print();
-
-
-private: // members
-
- // World coordinates of the node.
- WorldVector<double> *coords;
-
- // For interior nodes.
- Matrix<double> *A;
- Vector<double> *rec_uh;
- Vector<WorldVector<double> > *rec_grdUh;
-
- // For boundary, edge, face, of center nodes: interior neighbors nodes.
- // For interior vertices in uh-recovery: neighbors nodes.
- std::set<DegreeOfFreedom> *neighbors;
-
- friend class Recovery;
-};
-
-
-// ============================================================================
-// ===== class Recovery =======================================================
-// ============================================================================
-
-/**
- * \ingroup Recovery
- *
- * \brief
- * Recovering the gradient of a finite element function.
- */
-
-class Recovery
-{
-public:
- // constructor
- Recovery(int norm, int method_)
- : struct_vec(NULL), feSpace(NULL), matrix_fcts(NULL), method(method_)
- {
- n_monomials = 0;
- gradient = (norm == H1_NORM);
- };
-
- // Destructor?
-
- // Clear vector of recovery structures
- inline void clear()
+
+
+ /**
+ * \ingroup Recovery
+ *
+ * \brief
+ * Recovering the gradient of a finite element function.
+ */
+ class Recovery
{
- if (struct_vec)
- {
+ public:
+ Recovery(int norm, int method_)
+ : struct_vec(NULL),
+ feSpace(NULL),
+ matrix_fcts(NULL),
+ method(method_)
+ {
+ n_monomials = 0;
+ gradient = (norm == H1_NORM);
+ }
+
+ /// Clear vector of recovery structures
+ inline void clear()
+ {
+ if (struct_vec) {
DOFVector<RecoveryStructure>::Iterator SV_it(struct_vec, ALL_DOFS);
for (SV_it.reset(); !SV_it.end(); ++SV_it)
(*SV_it).clear();
}
- };
-
- /** \brief
- * Recovers flux or gradient of given DOFVector.
- */
- DOFVector<WorldVector<double> >*
- recovery(DOFVector<double> *uh,
- AbstractFunction<double, WorldVector<double> > *f_vec = NULL,
- AbstractFunction<double, double> *f_scal = NULL,
- DOFVector<double> *aux_vec = NULL);
-
- DOFVector<WorldVector<double> >*
- recovery(DOFVector<double> *uh, const FiniteElemSpace *fe_space,
- AbstractFunction<double, WorldVector<double> > *f_vec = NULL,
- AbstractFunction<double, double> *f_scal = NULL,
- DOFVector<double> *aux_vec = NULL);
-
- /** \brief
- * Computes higher order approximation of given DOFVector.
- */
- void recoveryUh(DOFVector<double> *uh, DOFVector<double> &rec_vec);
- DOFVector<double>* recoveryUh(DOFVector<double> *uh,
- const FiniteElemSpace *fe_space);
-
- // for test purposes
- void test(DOFVector<double> *uh, const FiniteElemSpace *fe_space);
-
-
-private: // Private functions.
-
- // Defines a new finite element space if necessary.
- void set_feSpace(const FiniteElemSpace *fe_space);
-
- // Fills vector of exponents of the monomials.
- int set_exponents(int degree);
-
- // Fills vector of recovery structures.
- void fill_struct_vec(DOFVector<double> *uh,
- AbstractFunction<double, WorldVector<double> > *f_vec = NULL,
- AbstractFunction<double, double> *f = NULL,
- DOFVector<double> *aux_vec = NULL);
-
- // Compute integrals defining matrix and vector on element
- // (continuous ZZ-recovery)
- void compute_integrals(DOFVector<double> *uh, ElInfo *elInfo,
- RecoveryStructure *rec_struct,
+ }
+
+ /// Recovers flux or gradient of given DOFVector.
+ DOFVector<WorldVector<double> >*
+ recovery(DOFVector<double> *uh,
+ AbstractFunction<double, WorldVector<double> > *f_vec = NULL,
+ AbstractFunction<double, double> *f_scal = NULL,
+ DOFVector<double> *aux_vec = NULL);
+
+ DOFVector<WorldVector<double> >*
+ recovery(DOFVector<double> *uh, const FiniteElemSpace *fe_space,
+ AbstractFunction<double, WorldVector<double> > *f_vec = NULL,
+ AbstractFunction<double, double> *f_scal = NULL,
+ DOFVector<double> *aux_vec = NULL);
+
+ /// Computes higher order approximation of given DOFVector.
+ void recoveryUh(DOFVector<double> *uh, DOFVector<double> &rec_vec);
+
+ DOFVector<double>* recoveryUh(DOFVector<double> *uh,
+ const FiniteElemSpace *fe_space);
+
+ /// For test purposes
+ void test(DOFVector<double> *uh, const FiniteElemSpace *fe_space);
+
+ private:
+ /// Defines a new finite element space if necessary.
+ void set_feSpace(const FiniteElemSpace *fe_space);
+
+ /// Fills vector of exponents of the monomials.
+ int set_exponents(int degree);
+
+ /// Fills vector of recovery structures.
+ void fill_struct_vec(DOFVector<double> *uh,
AbstractFunction<double, WorldVector<double> > *f_vec = NULL,
- AbstractFunction<double, double> *f_scal = NULL,
+ AbstractFunction<double, double> *f = NULL,
DOFVector<double> *aux_vec = NULL);
-
- // Compute integrals defining matrix and vector on element
- // (superconvergent patch recovery)
- void compute_interior_sums(DOFVector<double> *uh, ElInfo *elInfo,
- RecoveryStructure *rec_struct, Quadrature *quad,
- AbstractFunction<double, WorldVector<double> > *f_vec = NULL,
- AbstractFunction<double, double> *f_scal = NULL,
- DOFVector<double> *aux_vec = NULL);
-
- void compute_node_sums(DOFVector<double> *uh, ElInfo *elInfo,
- RecoveryStructure *rec_struct, DimVec<int> preDOFs,
- int n_vertices, int n_edges, int n_faces);
-
- void compute_sums_linear(DOFVector<double> *uh, ElInfo *elInfo,
+
+ /// Compute integrals defining matrix and vector on elemen (continuous ZZ-recovery)
+ void compute_integrals(DOFVector<double> *uh, ElInfo *elInfo,
RecoveryStructure *rec_struct,
- int vertex, DimVec<int> preDOFs, int n_vertices);
-
-private: // Private members.
-
- // Structure storing needed information.
- DOFVector<RecoveryStructure> *struct_vec;
-
- const FiniteElemSpace *feSpace; // Working finite element space.
-
- int n_monomials; // Number of monomials.
- Vector<WorldVector<int> > exponents; // Exponents of the monomials.
- Matrix<Monomial*> *matrix_fcts; // Functions for system matrix.
-
- bool gradient; // True if gradient or flux should be recovered;
- // false if seeking for higher order approximation of uh.
-
- int method; // 0: superconvergent patch recovery (discrete ZZ)
- // 1: local L2-averaging (continuous ZZ-recovery)
- // 2: simple averaging
-};
-
+ AbstractFunction<double, WorldVector<double> > *f_vec = NULL,
+ AbstractFunction<double, double> *f_scal = NULL,
+ DOFVector<double> *aux_vec = NULL);
+
+ /// Compute integrals defining matrix and vector on element (superconvergent patch recovery)
+ void compute_interior_sums(DOFVector<double> *uh, ElInfo *elInfo,
+ RecoveryStructure *rec_struct, Quadrature *quad,
+ AbstractFunction<double, WorldVector<double> > *f_vec = NULL,
+ AbstractFunction<double, double> *f_scal = NULL,
+ DOFVector<double> *aux_vec = NULL);
+
+ void compute_node_sums(DOFVector<double> *uh, ElInfo *elInfo,
+ RecoveryStructure *rec_struct, DimVec<int> preDOFs,
+ int n_vertices, int n_edges, int n_faces);
+
+ void compute_sums_linear(DOFVector<double> *uh, ElInfo *elInfo,
+ RecoveryStructure *rec_struct,
+ int vertex, DimVec<int> preDOFs, int n_vertices);
+
+ private:
+ /// Structure storing needed information.
+ DOFVector<RecoveryStructure> *struct_vec;
+
+ /// Working finite element space.
+ const FiniteElemSpace *feSpace;
+
+ /// Number of monomials.
+ int n_monomials;
+
+ /// Exponents of the monomials.
+ Vector<WorldVector<int> > exponents;
+
+ /// Functions for system matrix.
+ Matrix<Monomial*> *matrix_fcts;
+
+ /** \brief
+ * True if gradient or flux should be recovered.
+ * False if seeking for higher order approximation of uh.
+ */
+ bool gradient;
+
+ /** \brief
+ * 0: superconvergent patch recovery (discrete ZZ)
+ * 1: local L2-averaging (continuous ZZ-recovery)
+ * 2: simple averaging
+ */
+ int method;
+ };
+
+}
#endif
diff --git a/AMDiS/src/RecoveryEstimator.cc b/AMDiS/src/RecoveryEstimator.cc
index df129567d575e6b090b9303feefa845cc17574f3..7eac009fda4e268d96f9ce58cd122f7499899a2e 100644
--- a/AMDiS/src/RecoveryEstimator.cc
+++ b/AMDiS/src/RecoveryEstimator.cc
@@ -1,222 +1,202 @@
#include "RecoveryEstimator.h"
#include "Parameters.h"
-RecoveryEstimator::RecoveryEstimator(std::string name, DOFVector<double> *uh, int r)
- : Estimator(name, r),
- uh_(uh),
- relative_(0),
- C(1.0),
- method(0),
- feSpace(NULL),
- f_vec(NULL),
- f_scal(NULL),
- aux_vec(NULL),
- rec_struct(NULL)
-{
- FUNCNAME("RecoveryEstimator::constructor()");
-
- GET_PARAMETER(0, name + "->rec method", "%d", &method);
- GET_PARAMETER(0, name + "->rel error", "%d", &relative_);
-
- GET_PARAMETER(0, name + "->C", "%f", &C);
- C = C > 1e-25 ? sqr(C) : 1.0;
-
- if (norm == H1_NORM) {
- feSpace = uh->getFESpace();
- degree = feSpace->getBasisFcts()->getDegree();
-
- if ((degree <= 2) && (C != 1.0)) {
- WARNING("Recovery estimators in the H1_NORM usually very good for linear and quadratic finite element; normally you do not need to overwrite the default value of C\n");
- WAIT;
- }
- }
- else {
- degree = uh->getFESpace()->getBasisFcts()->getDegree() + 1;
+namespace AMDiS {
+
+ RecoveryEstimator::RecoveryEstimator(std::string name, DOFVector<double> *uh, int r)
+ : Estimator(name, r),
+ uh_(uh),
+ relative_(0),
+ C(1.0),
+ method(0),
+ feSpace(NULL),
+ f_vec(NULL),
+ f_scal(NULL),
+ aux_vec(NULL),
+ rec_struct(NULL)
+ {
+ FUNCNAME("RecoveryEstimator::constructor()");
- feSpace = FiniteElemSpace::provideFESpace(NULL,
- Lagrange::getLagrange(uh->getFESpace()->getMesh()->getDim(), degree),
- uh->getFESpace()->getMesh(),
- name + "->feSpace");
-
+ GET_PARAMETER(0, name + "->rec method", "%d", &method);
+ GET_PARAMETER(0, name + "->rel error", "%d", &relative_);
+
+ GET_PARAMETER(0, name + "->C", "%f", &C);
+ C = C > 1e-25 ? sqr(C) : 1.0;
+
+ if (norm == H1_NORM) {
+ feSpace = uh->getFESpace();
+ degree = feSpace->getBasisFcts()->getDegree();
+
+ if (degree <= 2 && C != 1.0) {
+ WARNING("Recovery estimators in the H1_NORM usually very good for linear and quadratic finite element; normally you do not need to overwrite the default value of C\n");
+ WAIT;
+ }
+ } else {
+ degree = uh->getFESpace()->getBasisFcts()->getDegree() + 1;
+ feSpace = FiniteElemSpace::provideFESpace(NULL,
+ Lagrange::getLagrange(uh->getFESpace()->getMesh()->getDim(), degree),
+ uh->getFESpace()->getMesh(),
+ name + "->feSpace");
+
+ if (method == 2) {
+ ERROR("Simple averaging only for the H1_NORM; using SPR instead\n");
+ WAIT;
+ method = 0;
+ }
+ }
- if (method == 2) {
- ERROR("Simple averaging only for the H1_NORM; using SPR instead\n");
+ if (method == 2 && degree !=1) {
+ ERROR("Simple averaging only for linear elements; using SPR instead\n");
WAIT;
method = 0;
}
- }
- if ((method == 2) && (degree !=1)) {
- ERROR("Simple averaging only for linear elements; using SPR instead\n");
- WAIT;
- method = 0;
+ rec_struct = new Recovery(norm, method);
}
- rec_struct = new Recovery(norm, method);
-}
-
-double RecoveryEstimator::estimate(double timestep)
-{
- FUNCNAME("RecoveryEstimator::estimate()");
-
- const BasisFunction *basFcts = uh_->getFESpace()->getBasisFcts();
- Mesh *mesh = uh_->getFESpace()->getMesh();
- int dim = mesh->getDim();
- int dow = Global::getGeo(WORLD);
- double h1Norm2 = 0.0;
- Flag flag =
- Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET | Mesh::FILL_GRD_LAMBDA;
- TraverseStack stack;
- ElInfo *elInfo = NULL;
-
- int i, j;
- double det, estEl, errEl, normEl;
- Element *el = NULL;
-
- if (norm == H1_NORM) // sets recovery gradient.
- {
+ double RecoveryEstimator::estimate(double timestep)
+ {
+ FUNCNAME("RecoveryEstimator::estimate()");
+
+ const BasisFunction *basFcts = uh_->getFESpace()->getBasisFcts();
+ Mesh *mesh = uh_->getFESpace()->getMesh();
+ int dim = mesh->getDim();
+ int dow = Global::getGeo(WORLD);
+ double h1Norm2 = 0.0;
+
+ if (norm == H1_NORM) { // sets recovery gradient.
if (method == 2)
- rec_grd = rec_struct->recovery(uh_, f_vec, f_scal, aux_vec);
+ rec_grd = rec_struct->recovery(uh_, f_vec, f_scal, aux_vec);
else
- rec_grd = rec_struct->recovery(uh_, feSpace, f_vec, f_scal, aux_vec);
+ rec_grd = rec_struct->recovery(uh_, feSpace, f_vec, f_scal, aux_vec);
+
rec_basFcts = rec_grd->getFESpace()->getBasisFcts();
- }
- else // sets higher-order recovery solution.
- {
- rec_uh = rec_struct->recoveryUh(uh_, feSpace);
+ } else { // sets higher-order recovery solution.
+ rec_uh = rec_struct->recoveryUh(uh_, feSpace);
rec_basFcts = rec_uh->getFESpace()->getBasisFcts();
}
- int deg = 2 * std::max(basFcts->getDegree(),
- rec_basFcts->getDegree());
- Quadrature *quad = Quadrature::provideQuadrature(dim, deg);
- int numPoints = quad->getNumPoints();
-
- WorldVector<double> quad_pt;
- WorldVector<double> *grdAtQP = new WorldVector<double>[numPoints];
- WorldVector<double> *recoveryGrdAtQP = new WorldVector<double>[numPoints];
-
- double *uhAtQP = new double[numPoints];
- double *recoveryUhAtQP = new double[numPoints];
-
- FastQuadrature *quadFast =
- FastQuadrature::provideFastQuadrature(basFcts, *quad,
- INIT_PHI | INIT_GRD_PHI);
- FastQuadrature *rec_quadFast =
- FastQuadrature::provideFastQuadrature(rec_basFcts, *quad,
- INIT_PHI | INIT_GRD_PHI);
-
- // clear error indicators
- elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
- while (elInfo)
- {
+ int deg = 2 * std::max(basFcts->getDegree(), rec_basFcts->getDegree());
+ Quadrature *quad = Quadrature::provideQuadrature(dim, deg);
+ int nPoints = quad->getNumPoints();
+
+ WorldVector<double> quad_pt;
+ WorldVector<double> *grdAtQP = new WorldVector<double>[nPoints];
+ WorldVector<double> *recoveryGrdAtQP = new WorldVector<double>[nPoints];
+
+ double *uhAtQP = new double[nPoints];
+ double *recoveryUhAtQP = new double[nPoints];
+
+ FastQuadrature *quadFast =
+ FastQuadrature::provideFastQuadrature(basFcts, *quad, INIT_PHI | INIT_GRD_PHI);
+ FastQuadrature *rec_quadFast =
+ FastQuadrature::provideFastQuadrature(rec_basFcts, *quad, INIT_PHI | INIT_GRD_PHI);
+
+ // clear error indicators
+ TraverseStack stack;
+ ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
+ while (elInfo) {
elInfo->getElement()->setEstimation(0.0, row);
elInfo = stack.traverseNext(elInfo);
}
- // traverse mesh
- est_sum = est_max = est_t_sum = 0.0;
-
- elInfo = stack.traverseFirst(mesh, -1, flag);
- while (elInfo)
- {
- el = elInfo->getElement();
- det = elInfo->getDet();
- errEl = 0.0;
-
- if (norm == H1_NORM)
- {
- // get gradient and recovery gradient at quadrature points
- uh_->getGrdAtQPs(elInfo, NULL, quadFast, grdAtQP);
- rec_grd->getVecAtQPs(elInfo, NULL, rec_quadFast, recoveryGrdAtQP);
- if (f_scal)
- {
- if (aux_vec)
- aux_vec->getVecAtQPs(elInfo, NULL, quadFast, uhAtQP);
- else
- uh_->getVecAtQPs(elInfo, NULL, quadFast, uhAtQP);
- }
-
- // calc h1 error
- for (i=0; i<numPoints; i++)
- {
- double err2 = 0.0;
- double fAtQP = 1.0;
- if (f_scal)
- fAtQP = (*f_scal)(uhAtQP[i]);
- if (f_vec)
- {
- elInfo->coordToWorld(quad->getLambda(i), quad_pt);
- fAtQP = (*f_vec)(quad_pt);
- }
-
- for (j=0; j<dow; j++)
- err2 += sqr(recoveryGrdAtQP[i][j] - fAtQP * grdAtQP[i][j]);
- errEl += quad->getWeight(i) * err2;
- }
- }
- else
- {
- // get vector and recovery vector at quadrature points
- uh_->getVecAtQPs(elInfo, NULL, quadFast, uhAtQP);
- rec_uh->getVecAtQPs(elInfo, NULL, rec_quadFast, recoveryUhAtQP);
-
- // calc l2 error
- for (i=0; i<numPoints; i++)
- errEl += quad->getWeight(i) * sqr(recoveryUhAtQP[i] - uhAtQP[i]);
+ // traverse mesh
+ est_sum = est_max = est_t_sum = 0.0;
+
+ elInfo = stack.traverseFirst(mesh, -1,
+ Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS |
+ Mesh::FILL_DET | Mesh::FILL_GRD_LAMBDA);
+ while (elInfo) {
+ Element *el = elInfo->getElement();
+ double det = elInfo->getDet();
+ double errEl = 0.0;
+
+ if (norm == H1_NORM) {
+ // get gradient and recovery gradient at quadrature points
+ uh_->getGrdAtQPs(elInfo, NULL, quadFast, grdAtQP);
+ rec_grd->getVecAtQPs(elInfo, NULL, rec_quadFast, recoveryGrdAtQP);
+ if (f_scal) {
+ if (aux_vec)
+ aux_vec->getVecAtQPs(elInfo, NULL, quadFast, uhAtQP);
+ else
+ uh_->getVecAtQPs(elInfo, NULL, quadFast, uhAtQP);
}
- estEl = C * det * errEl;
+ // calc h1 error
+ for (int i = 0; i < nPoints; i++) {
+ double err2 = 0.0;
+ double fAtQP = 1.0;
+ if (f_scal)
+ fAtQP = (*f_scal)(uhAtQP[i]);
+ if (f_vec) {
+ elInfo->coordToWorld(quad->getLambda(i), quad_pt);
+ fAtQP = (*f_vec)(quad_pt);
+ }
+
+ for (int j = 0; j < dow; j++)
+ err2 += sqr(recoveryGrdAtQP[i][j] - fAtQP * grdAtQP[i][j]);
+ errEl += quad->getWeight(i) * err2;
+ }
+ } else {
+ // get vector and recovery vector at quadrature points
+ uh_->getVecAtQPs(elInfo, NULL, quadFast, uhAtQP);
+ rec_uh->getVecAtQPs(elInfo, NULL, rec_quadFast, recoveryUhAtQP);
+
+ // calc l2 error
+ for (int i = 0; i < nPoints; i++)
+ errEl += quad->getWeight(i) * sqr(recoveryUhAtQP[i] - uhAtQP[i]);
+ }
+
+ double estEl = C * det * errEl;
el->setEstimation(estEl, row);
est_sum += estEl;
est_max = std::max(est_max, estEl);
- if (relative_)
- {
- normEl = 0.0;
-
- if (norm == H1_NORM)
- for (i=0; i<numPoints; i++)
- {
- double norm2 = 0.0;
- for (j = 0; j < dow; j++)
- norm2 += sqr(recoveryGrdAtQP[i][j]);
- normEl += quad->getWeight(i) * norm2;
- }
- else
- for (i=0; i<numPoints; i++)
- normEl += quad->getWeight(i) * sqr(recoveryUhAtQP[i]);
-
- h1Norm2 += det * normEl;
+ if (relative_) {
+ double normEl = 0.0;
+
+ if (norm == H1_NORM) {
+ for (int i = 0; i < nPoints; i++) {
+ double norm2 = 0.0;
+ for (int j = 0; j < dow; j++)
+ norm2 += sqr(recoveryGrdAtQP[i][j]);
+ normEl += quad->getWeight(i) * norm2;
+ }
+ } else {
+ for (int i = 0; i < nPoints; i++)
+ normEl += quad->getWeight(i) * sqr(recoveryUhAtQP[i]);
}
+ h1Norm2 += det * normEl;
+ }
+
elInfo = stack.traverseNext(elInfo);
}
- // Computing relative errors
- if (relative_)
- {
+ // Computing relative errors
+ if (relative_) {
elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
- while (elInfo)
- {
- estEl = elInfo->getElement()->getEstimation(row);
- estEl /= h1Norm2;
- elInfo->getElement()->setEstimation(estEl, row);
- elInfo = stack.traverseNext(elInfo);
- }
-
+ while (elInfo) {
+ double estEl = elInfo->getElement()->getEstimation(row);
+ estEl /= h1Norm2;
+ elInfo->getElement()->setEstimation(estEl, row);
+ elInfo = stack.traverseNext(elInfo);
+ }
+
est_max /= h1Norm2;
est_sum /= h1Norm2;
}
- est_sum = sqrt(est_sum);
-
- MSG("estimate = %.8e\n", est_sum);
-
- delete [] grdAtQP;
- delete [] recoveryGrdAtQP;
- delete [] uhAtQP;
- delete [] recoveryUhAtQP;
+ est_sum = sqrt(est_sum);
+
+ MSG("estimate = %.8e\n", est_sum);
+
+ delete [] grdAtQP;
+ delete [] recoveryGrdAtQP;
+ delete [] uhAtQP;
+ delete [] recoveryUhAtQP;
+
+ return est_sum;
+ }
- return(est_sum);
}
diff --git a/AMDiS/src/ResidualEstimator.cc b/AMDiS/src/ResidualEstimator.cc
index 18e9cb025d41ca9e281a6636ef9e90df3c99eaf3..13ae3f63f22a80eb1097282b5ab88cc90c0f25cc 100644
--- a/AMDiS/src/ResidualEstimator.cc
+++ b/AMDiS/src/ResidualEstimator.cc
@@ -40,9 +40,6 @@ namespace AMDiS {
FUNCNAME("ResidualEstimator::init()");
timestep = ts;
-
- mesh = uh[row == -1 ? 0 : row]->getFESpace()->getMesh();
-
nSystems = static_cast<int>(uh.size());
TEST_EXIT_DBG(nSystems > 0)("no system set\n");
@@ -111,7 +108,7 @@ namespace AMDiS {
Mesh::CALL_LEAF_EL;
neighInfo = mesh->createNewElInfo();
- // prepare date for computing jump residual
+ // === Prepare date for computing jump residual. ===
if (C1 > 0.0 && dim > 1) {
surfaceQuad = Quadrature::provideQuadrature(dim - 1, degree);
nPointsSurface = surfaceQuad->getNumPoints();
@@ -122,6 +119,15 @@ namespace AMDiS {
nNeighbours = Global::getGeo(NEIGH, dim);
lambdaNeigh = new DimVec<WorldVector<double> >(dim, NO_INIT);
lambda = new DimVec<double>(dim, NO_INIT);
+
+ secondOrderTerms.resize(nSystems);
+ for (int system = 0; system < nSystems; system++) {
+ secondOrderTerms[system] = false;
+
+ for (std::vector<Operator*>::iterator it = matrix[system]->getOperators().begin();
+ it != matrix[system]->getOperators().end(); ++it)
+ secondOrderTerms[system] = secondOrderTerms[system] || (*it)->secondOrderTerms();
+ }
}
}
@@ -188,24 +194,18 @@ namespace AMDiS {
}
- void ResidualEstimator::estimateElement(ElInfo *elInfo)
+ void ResidualEstimator::estimateElement(ElInfo *elInfo, DualElInfo *dualElInfo)
{
FUNCNAME("ResidualEstimator::estimateElement()");
TEST_EXIT_DBG(nSystems > 0)("no system set\n");
- double val = 0.0;
- std::vector<Operator*>::iterator it;
- std::vector<double*>::iterator itfac;
Element *el = elInfo->getElement();
- double det = elInfo->getDet();
- const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
double est_el = el->getEstimation(row);
- double h2 = h2_from_det(det, dim);
-
- for (int iq = 0; iq < nPoints; iq++)
- riq[iq] = 0.0;
+ std::vector<Operator*>::iterator it;
+ std::vector<double*>::iterator itfac;
+ // === Init assemblers. ===
for (int system = 0; system < nSystems; system++) {
if (matrix[system] == NULL)
continue;
@@ -213,8 +213,6 @@ namespace AMDiS {
DOFMatrix *dofMat = const_cast<DOFMatrix*>(matrix[system]);
DOFVector<double> *dofVec = const_cast<DOFVector<double>*>(fh[system]);
- // === init assemblers ===
-
for (it = dofMat->getOperatorsBegin(), itfac = dofMat->getOperatorEstFactorBegin();
it != dofMat->getOperatorsEnd(); ++it, ++itfac)
if (*itfac == NULL || **itfac != 0.0)
@@ -223,10 +221,47 @@ namespace AMDiS {
if (C0 > 0.0)
for (it = dofVec->getOperatorsBegin(); it != dofVec->getOperatorsEnd(); ++it)
(*it)->getAssembler(omp_get_thread_num())->initElement(elInfo, NULL, quad);
-
+ }
+
+
+ // === Compute element residuals and time error estimation. ===
+ if (C0 || C3)
+ est_el += computeElementResidual(elInfo, dualElInfo);
+
+ // === Compute jump residuals. ===
+ if (C1 && dim > 1)
+ est_el += computeJumpResidual(elInfo, dualElInfo);
+
+ // === Update global residual variables. ===
+ el->setEstimation(est_el, row);
+ el->setMark(0);
+ est_sum += est_el;
+ est_max = max(est_max, est_el);
+ }
+
+
+ double ResidualEstimator::computeElementResidual(ElInfo *elInfo,
+ DualElInfo *dualElInfo)
+ {
+ FUNCNAME("ResidualEstimator::computeElementResidual()");
+
+ TEST_EXIT(dualElInfo)("Not yet implemented!\n");
+
+ std::vector<Operator*>::iterator it;
+ std::vector<double*>::iterator itfac;
+ double det = elInfo->getDet();
+ double h2 = h2_from_det(det, dim);
+
+ for (int iq = 0; iq < nPoints; iq++)
+ riq[iq] = 0.0;
+
+ for (int system = 0; system < nSystems; system++) {
+ if (matrix[system] == NULL)
+ continue;
+
if (timestep && uhOld[system]) {
TEST_EXIT_DBG(uhOld[system])("no uhOld\n");
- uhOld[system]->getLocalVector(el, uhOldEl[system]);
+ uhOld[system]->getLocalVector(elInfo->getElement(), uhOldEl[system]);
// === Compute time error. ===
@@ -235,12 +270,12 @@ namespace AMDiS {
uhOld[system]->getVecAtQPs(elInfo, NULL, quadFast[system], uhOldQP);
if (C3 > 0.0 && uhOldQP && system == std::max(row, 0)) {
- val = 0.0;
+ double result = 0.0;
for (int iq = 0; iq < nPoints; iq++) {
double tiq = (uhQP[iq] - uhOldQP[iq]);
- val += quad->getWeight(iq) * tiq * tiq;
+ result += quad->getWeight(iq) * tiq * tiq;
}
- double v = C3 * det * val;
+ double v = C3 * det * result;
est_t_sum += v;
est_t_max = max(est_t_max, v);
}
@@ -248,7 +283,10 @@ namespace AMDiS {
}
// === Compute element residual. ===
- if (C0 > 0.0) {
+ if (C0 > 0.0) {
+ DOFMatrix *dofMat = const_cast<DOFMatrix*>(matrix[system]);
+ DOFVector<double> *dofVec = const_cast<DOFVector<double>*>(fh[system]);
+
for (it = dofMat->getOperatorsBegin(), itfac = dofMat->getOperatorEstFactorBegin();
it != dofMat->getOperatorsEnd(); ++it, ++itfac) {
if (*itfac == NULL || **itfac != 0.0) {
@@ -286,182 +324,167 @@ namespace AMDiS {
}
// add integral over r square
- val = 0.0;
+ double result = 0.0;
for (int iq = 0; iq < nPoints; iq++)
- val += quad->getWeight(iq) * riq[iq] * riq[iq];
+ result += quad->getWeight(iq) * riq[iq] * riq[iq];
if (timestep != 0.0 || norm == NO_NORM || norm == L2_NORM)
- val = C0 * h2 * h2 * det * val;
+ result = C0 * h2 * h2 * det * result;
else
- val = C0 * h2 * det * val;
+ result = C0 * h2 * det * result;
- est_el += val;
+ return result;
+ }
- // === Compute jump residuals. ===
+ double ResidualEstimator::computeJumpResidual(ElInfo *elInfo,
+ DualElInfo *dualElInfo)
+ {
+ FUNCNAME("ResidualEstimator::computeJumpResidual()");
+
+ double result = 0.0;
+ int dow = Global::getGeo(WORLD);
+ Element *el = elInfo->getElement();
+ const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
+ double det = elInfo->getDet();
+ double h2 = h2_from_det(det, dim);
- if (C1 && dim > 1) {
- int dow = Global::getGeo(WORLD);
+ for (int face = 0; face < nNeighbours; face++) {
+ Element *neigh = const_cast<Element*>(elInfo->getNeighbour(face));
- for (int face = 0; face < nNeighbours; face++) {
- Element *neigh = const_cast<Element*>(elInfo->getNeighbour(face));
- if (neigh && neigh->getMark()) {
- int oppV = elInfo->getOppVertex(face);
-
- el->sortFaceIndices(face, &faceIndEl);
- neigh->sortFaceIndices(oppV, &faceIndNeigh);
-
- neighInfo->setElement(const_cast<Element*>(neigh));
- neighInfo->setFillFlag(Mesh::FILL_COORDS);
-
- for (int i = 0; i < dow; i++)
- neighInfo->getCoord(oppV)[i] = elInfo->getOppCoord(face)[i];
-
- // periodic leaf data ?
- ElementData *ldp = el->getElementData()->getElementData(PERIODIC);
-
- bool periodicCoords = false;
-
- if (ldp) {
- std::list<LeafDataPeriodic::PeriodicInfo>::iterator it;
- std::list<LeafDataPeriodic::PeriodicInfo>& infoList =
- dynamic_cast<LeafDataPeriodic*>(ldp)->getInfoList();
-
- for (it = infoList.begin(); it != infoList.end(); ++it) {
- if (it->elementSide == face) {
- for (int i = 0; i < dim; i++) {
- int i1 = faceIndEl[i];
- int i2 = faceIndNeigh[i];
-
- int j = 0;
- for (; j < dim; j++) {
- if (i1 == el->getVertexOfPosition(INDEX_OF_DIM(dim - 1,
- dim),
- face,
- j)) {
- break;
- }
- }
-
- TEST_EXIT_DBG(j != dim)("vertex i1 not on face ???\n");
-
- neighInfo->getCoord(i2) = (*(it->periodicCoords))[j];
- }
- periodicCoords = true;
- break;
- }
- }
- }
+ if (!(neigh && neigh->getMark()))
+ continue;
- if (!periodicCoords) {
+ int oppV = elInfo->getOppVertex(face);
+
+ el->sortFaceIndices(face, &faceIndEl);
+ neigh->sortFaceIndices(oppV, &faceIndNeigh);
+ neighInfo->setElement(const_cast<Element*>(neigh));
+ neighInfo->setFillFlag(Mesh::FILL_COORDS);
+
+ for (int i = 0; i < dow; i++)
+ neighInfo->getCoord(oppV)[i] = elInfo->getOppCoord(face)[i];
+
+ // periodic leaf data ?
+ ElementData *ldp = el->getElementData()->getElementData(PERIODIC);
+ bool periodicCoords = false;
+
+ if (ldp) {
+ typedef std::list<LeafDataPeriodic::PeriodicInfo> PerInfList;
+ PerInfList& infoList = dynamic_cast<LeafDataPeriodic*>(ldp)->getInfoList();
+
+ for (PerInfList::iterator it = infoList.begin(); it != infoList.end(); ++it) {
+ if (it->elementSide == face) {
for (int i = 0; i < dim; i++) {
int i1 = faceIndEl[i];
int i2 = faceIndNeigh[i];
- for (int j = 0; j < dow; j++)
- neighInfo->getCoord(i2)[j] = elInfo->getCoord(i1)[j];
- }
- }
-
- Parametric *parametric = mesh->getParametric();
- if (parametric)
- neighInfo = parametric->addParametricInfo(neighInfo);
- double detNeigh = abs(neighInfo->calcGrdLambda(*lambdaNeigh));
+ int j = 0;
+ for (; j < dim; j++)
+ if (i1 == el->getVertexOfPosition(INDEX_OF_DIM(dim - 1, dim), face, j))
+ break;
- for (int iq = 0; iq < nPointsSurface; iq++)
- jump[iq].set(0.0);
-
- for (int system = 0; system < nSystems; system++) {
- if (matrix[system] == NULL)
- continue;
+ TEST_EXIT_DBG(j != dim)("vertex i1 not on face ???\n");
- uh[system]->getLocalVector(el, uhEl[system]);
- uh[system]->getLocalVector(neigh, uhNeigh[system]);
-
- for (int iq = 0; iq < nPointsSurface; iq++) {
- (*lambda)[face] = 0.0;
- for (int i = 0; i < dim; i++)
- (*lambda)[faceIndEl[i]] = surfaceQuad->getLambda(iq, i);
-
- basFcts[system]->evalGrdUh(*lambda,
- grdLambda,
- uhEl[system],
- &grdUhEl[iq]);
-
- (*lambda)[oppV] = 0.0;
- for (int i = 0; i < dim; i++)
- (*lambda)[faceIndNeigh[i]] = surfaceQuad->getLambda(iq, i);
-
- basFcts[system]->evalGrdUh(*lambda,
- *lambdaNeigh,
- uhNeigh[system],
- &grdUhNeigh[iq]);
-
- grdUhEl[iq] -= grdUhNeigh[iq];
- }
-
- std::vector<double*>::iterator fac;
-
- for (it = const_cast<DOFMatrix*>(matrix[system])->getOperatorsBegin(),
- fac = const_cast<DOFMatrix*>(matrix[system])->getOperatorEstFactorBegin();
- it != const_cast<DOFMatrix*>(matrix[system])->getOperatorsEnd();
- ++it, ++fac) {
-
- if (*fac == NULL || **fac != 0.0) {
- for (int iq = 0; iq < nPointsSurface; iq++)
- localJump[iq].set(0.0);
-
- (*it)->weakEvalSecondOrder(nPointsSurface,
- grdUhEl.getValArray(),
- localJump.getValArray());
- double factor = *fac ? **fac : 1.0;
- if (factor != 1.0)
- for (int i = 0; i < nPointsSurface; i++)
- localJump[i] *= factor;
-
- for (int i = 0; i < nPointsSurface; i++)
- jump[i] += localJump[i];
- }
+ neighInfo->getCoord(i2) = (*(it->periodicCoords))[j];
}
+ periodicCoords = true;
+ break;
}
+ }
+ } // if (ldp)
+
+ if (!periodicCoords) {
+ for (int i = 0; i < dim; i++) {
+ int i1 = faceIndEl[i];
+ int i2 = faceIndNeigh[i];
+ for (int j = 0; j < dow; j++)
+ neighInfo->getCoord(i2)[j] = elInfo->getCoord(i1)[j];
+ }
+ }
+
+ Parametric *parametric = mesh->getParametric();
+ if (parametric)
+ neighInfo = parametric->addParametricInfo(neighInfo);
- val = 0.0;
- for (int iq = 0; iq < nPointsSurface; iq++)
- val += surfaceQuad->getWeight(iq) * (jump[iq] * jump[iq]);
-
- double d = 0.5 * (det + detNeigh);
-
- if (norm == NO_NORM || norm == L2_NORM)
- val *= C1 * h2_from_det(d, dim) * d;
- else
- val *= C1 * d;
+ double detNeigh = abs(neighInfo->calcGrdLambda(*lambdaNeigh));
+
+ for (int iq = 0; iq < nPointsSurface; iq++)
+ jump[iq].set(0.0);
+
+ for (int system = 0; system < nSystems; system++) {
+ if (matrix[system] == NULL || secondOrderTerms[system] == false)
+ continue;
- if (parametric)
- neighInfo = parametric->removeParametricInfo(neighInfo);
+ uh[system]->getLocalVector(el, uhEl[system]);
+ uh[system]->getLocalVector(neigh, uhNeigh[system]);
+
+ for (int iq = 0; iq < nPointsSurface; iq++) {
+ (*lambda)[face] = 0.0;
+ for (int i = 0; i < dim; i++)
+ (*lambda)[faceIndEl[i]] = surfaceQuad->getLambda(iq, i);
+
+ basFcts[system]->evalGrdUh(*lambda, grdLambda, uhEl[system], &grdUhEl[iq]);
+
+ (*lambda)[oppV] = 0.0;
+ for (int i = 0; i < dim; i++)
+ (*lambda)[faceIndNeigh[i]] = surfaceQuad->getLambda(iq, i);
+
+ basFcts[system]->evalGrdUh(*lambda, *lambdaNeigh, uhNeigh[system], &grdUhNeigh[iq]);
+
+ grdUhEl[iq] -= grdUhNeigh[iq];
+ } // for iq
+
+ std::vector<double*>::iterator fac;
+ std::vector<Operator*>::iterator it;
+ DOFMatrix *mat = const_cast<DOFMatrix*>(matrix[system]);
+ for (it = mat->getOperatorsBegin(), fac = mat->getOperatorEstFactorBegin();
+ it != mat->getOperatorsEnd(); ++it, ++fac) {
+
+ if (*fac == NULL || **fac != 0.0) {
+ for (int iq = 0; iq < nPointsSurface; iq++)
+ localJump[iq].set(0.0);
+
+ (*it)->weakEvalSecondOrder(nPointsSurface, &(grdUhEl[0]), &(localJump[0]));
- neigh->setEstimation(neigh->getEstimation(row) + val, row);
- est_el += val;
- }
- }
-
- val = fh[std::max(row, 0)]->getBoundaryManager()->
- boundResidual(elInfo, matrix[std::max(row, 0)], uh[std::max(row, 0)]);
+ double factor = *fac ? **fac : 1.0;
+ if (factor != 1.0)
+ for (int i = 0; i < nPointsSurface; i++)
+ localJump[i] *= factor;
+
+ for (int i = 0; i < nPointsSurface; i++)
+ jump[i] += localJump[i];
+ }
+ } // for (it = ...
+ } // for system
+
+ double val = 0.0;
+ for (int iq = 0; iq < nPointsSurface; iq++)
+ val += surfaceQuad->getWeight(iq) * (jump[iq] * jump[iq]);
+ double d = 0.5 * (det + detNeigh);
+
if (norm == NO_NORM || norm == L2_NORM)
- val *= C1 * h2;
+ val *= C1 * h2_from_det(d, dim) * d;
else
- val *= C1;
-
- est_el += val;
- }
-
-
- el->setEstimation(est_el, row);
-
- est_sum += est_el;
- est_max = max(est_max, est_el);
+ val *= C1 * d;
+
+ if (parametric)
+ neighInfo = parametric->removeParametricInfo(neighInfo);
+
+ neigh->setEstimation(neigh->getEstimation(row) + val, row);
+ result += val;
+ } // for face
+
+ double val = fh[std::max(row, 0)]->getBoundaryManager()->
+ boundResidual(elInfo, matrix[std::max(row, 0)], uh[std::max(row, 0)]);
+ if (norm == NO_NORM || norm == L2_NORM)
+ val *= C1 * h2;
+ else
+ val *= C1;
+ result += val;
- elInfo->getElement()->setMark(0);
+ return result;
}
diff --git a/AMDiS/src/ResidualEstimator.h b/AMDiS/src/ResidualEstimator.h
index e21d18af473a775a80990db180b4dff3c9f5e601..ba737d6db20010e60734e0d9c0972ef630ff36e7 100644
--- a/AMDiS/src/ResidualEstimator.h
+++ b/AMDiS/src/ResidualEstimator.h
@@ -81,11 +81,22 @@ namespace AMDiS {
/// Constructor.
ResidualEstimator(std::string name, int r);
- virtual void init(double timestep);
+ void init(double timestep);
- virtual void estimateElement(ElInfo *elInfo);
+ /** \brief
+ * Estimates the error on an element. For more information about the parameter,
+ * see the description \ref Estimator::estimateElement.
+ */
+ void estimateElement(ElInfo *elInfo, DualElInfo *dualElInfo = NULL);
- virtual void exit(bool output = true);
+ void exit(bool output = true);
+
+ protected:
+ /// Computes the element residual for a given element.
+ double computeElementResidual(ElInfo *elInfo, DualElInfo *dualElInfo);
+
+ /// Computes the jump residual for a given element.
+ double computeJumpResidual(ElInfo *elInfo, DualElInfo *dualElInfo);
protected:
/// Constant in front of element residual
@@ -139,13 +150,13 @@ namespace AMDiS {
int nPointsSurface;
- Vector<WorldVector<double> > grdUhEl;
+ std::vector<WorldVector<double> > grdUhEl;
- Vector<WorldVector<double> > grdUhNeigh;
+ std::vector<WorldVector<double> > grdUhNeigh;
- Vector<WorldVector<double> > jump;
+ std::vector<WorldVector<double> > jump;
- Vector<WorldVector<double> > localJump;
+ std::vector<WorldVector<double> > localJump;
WorldVector<int> faceIndEl;
@@ -157,6 +168,13 @@ namespace AMDiS {
/// Maximal number of neighbours an element may have in the used dimension.
int nNeighbours;
+
+ /** \brief
+ * Defines for every system if there are second order terms. These values
+ * are used to ommit computations of the jump residual that is defined
+ * only on second order terms.
+ */
+ std::vector<bool> secondOrderTerms;
};
}