Commit cc0849b7 authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

Work on residual estimator for multiple meshes in problem definition.

parent d51a4a4b
......@@ -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;
}
......
......@@ -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;
}
}
......@@ -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;
......
......@@ -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;
}
......
......@@ -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,
......
......@@ -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);
}
}
......
......@@ -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,
......
#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);
}
}
}
......@@ -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)
{