From ecbb27791f67c1971510debb3aecd0ff0524c932 Mon Sep 17 00:00:00 2001
From: Thomas Witkowski <thomas.witkowski@gmx.de>
Date: Fri, 29 Aug 2008 11:20:33 +0000
Subject: [PATCH] * Exact error computation
---
AMDiS/src/ProblemNonLin.cc | 2 +-
AMDiS/src/ProblemVec.cc | 96 ++++++++++++++++++++++++++++----------
AMDiS/src/ProblemVec.h | 13 ++++--
3 files changed, 82 insertions(+), 29 deletions(-)
diff --git a/AMDiS/src/ProblemNonLin.cc b/AMDiS/src/ProblemNonLin.cc
index 5f4bb972..6a8309da 100644
--- a/AMDiS/src/ProblemNonLin.cc
+++ b/AMDiS/src/ProblemNonLin.cc
@@ -188,7 +188,7 @@ namespace AMDiS {
// for all elements ...
for (int i = 0; i < nComponents; i++) {
- elInfo = stack.traverseFirst(componentMeshes_[i], -1,
+ elInfo = stack.traverseFirst(componentMeshes[i], -1,
Mesh::CALL_LEAF_EL |
Mesh::FILL_BOUND |
Mesh::FILL_COORDS |
diff --git a/AMDiS/src/ProblemVec.cc b/AMDiS/src/ProblemVec.cc
index 36104a09..ea9ad47a 100644
--- a/AMDiS/src/ProblemVec.cc
+++ b/AMDiS/src/ProblemVec.cc
@@ -20,6 +20,7 @@
#include "RobinBC.h"
#include "PeriodicBC.h"
#include "Lagrange.h"
+#include "Flag.h"
namespace AMDiS {
@@ -45,7 +46,7 @@ namespace AMDiS {
adoptFlag.isSet(INIT_SYSTEM) ||
adoptFlag.isSet(INIT_FE_SPACE))) {
meshes_ = adoptProblem->getMeshes();
- componentMeshes_ = adoptProblem->componentMeshes_;
+ componentMeshes = adoptProblem->componentMeshes;
refinementManager_ = adoptProblem->refinementManager_;
coarseningManager_ = adoptProblem->coarseningManager_;
}
@@ -180,7 +181,7 @@ namespace AMDiS {
{
FUNCNAME("ProblemVec::createMesh()");
- componentMeshes_.resize(nComponents);
+ componentMeshes.resize(nComponents);
std::map<int, Mesh*> meshForRefinementSet;
char number[3];
@@ -203,7 +204,7 @@ namespace AMDiS {
meshForRefinementSet[refSet] = newMesh;
meshes_.push_back(newMesh);
}
- componentMeshes_[i] = meshForRefinementSet[refSet];
+ componentMeshes[i] = meshForRefinementSet[refSet];
}
switch(dim) {
case 1:
@@ -243,17 +244,17 @@ namespace AMDiS {
TEST_EXIT(componentSpaces[i] == NULL)("feSpace already created\n");
- if (feSpaceMap[std::pair<Mesh*, int>(componentMeshes_[i], degree)] == NULL) {
+ if (feSpaceMap[std::pair<Mesh*, int>(componentMeshes[i], degree)] == NULL) {
FiniteElemSpace *newFESpace =
FiniteElemSpace::provideFESpace(NULL,
Lagrange::getLagrange(dim, degree),
- componentMeshes_[i],
+ componentMeshes[i],
name_ + "->feSpace");
- feSpaceMap[std::pair<Mesh*, int>(componentMeshes_[i], degree)] = newFESpace;
+ feSpaceMap[std::pair<Mesh*, int>(componentMeshes[i], degree)] = newFESpace;
feSpaces.push_back(newFESpace);
}
componentSpaces[i] =
- feSpaceMap[std::pair<Mesh*, int>(componentMeshes_[i], degree)];
+ feSpaceMap[std::pair<Mesh*, int>(componentMeshes[i], degree)];
}
// create dof admin for vertex dofs if neccessary
@@ -459,14 +460,14 @@ namespace AMDiS {
std::vector< DOFVector<double>* > solutionList(nComponents);
for (int i = 0; i < nComponents; i++) {
- TEST_EXIT(componentMeshes_[0] == componentMeshes_[i])
+ TEST_EXIT(componentMeshes[0] == componentMeshes[i])
("All Meshes have to be equal to write a vector file.\n");
solutionList[i] = solution_->getDOFVector(i);
}
fileWriters_.push_back(NEW FileWriter(numberedName,
- componentMeshes_[0],
+ componentMeshes[0],
solutionList));
}
@@ -481,7 +482,7 @@ namespace AMDiS {
if (filename != "") {
fileWriters_.push_back(NEW FileWriter(numberedName,
- componentMeshes_[i],
+ componentMeshes[i],
solution_->getDOFVector(i)));
}
}
@@ -549,18 +550,7 @@ namespace AMDiS {
#endif
if (computeExactError) {
- for (int i = 0; i < nComponents; i++) {
- DOFVector<double> *tmp = NEW DOFVector<double>(componentSpaces[i], "tmp");
- tmp->interpol(exactSolutionFcts[i]);
- // double t = max(abs(tmp->max()), abs(tmp->min()));
- double t = tmp->absMax();
- *tmp -= *(solution_->getDOFVector(i));
- double l2Error = tmp->L2Norm();
- double maxError = tmp->absMax() / t;
- MSG("L2 error = %.8e\n", l2Error);
- MSG("L-inf error = %.8e\n", maxError);
- DELETE tmp;
- }
+ computeError(adaptInfo);
} else {
for (int i = 0; i < nComponents; i++) {
Estimator *scalEstimator = estimator_[i];
@@ -600,7 +590,7 @@ namespace AMDiS {
Flag markFlag = 0;
for (int i = 0; i < nComponents; i++) {
if (marker[i]) {
- markFlag |= marker[i]->markMesh(adaptInfo, componentMeshes_[i]);
+ markFlag |= marker[i]->markMesh(adaptInfo, componentMeshes[i]);
} else {
WARNING("no marker for component %d\n", i);
}
@@ -745,7 +735,7 @@ namespace AMDiS {
}
TraverseStack stack;
- ElInfo *elInfo = stack.traverseFirst(componentMeshes_[i], -1, assembleFlag);
+ ElInfo *elInfo = stack.traverseFirst(componentMeshes[i], -1, assembleFlag);
while (elInfo) {
if (useGetBound_) {
@@ -787,7 +777,7 @@ namespace AMDiS {
solution_->getDOFVector(i)->getBoundaryManager()->initVector(solution_->getDOFVector(i));
TraverseStack stack;
- ElInfo *elInfo = stack.traverseFirst(componentMeshes_[i], -1, assembleFlag);
+ ElInfo *elInfo = stack.traverseFirst(componentMeshes[i], -1, assembleFlag);
while (elInfo) {
if (rhs_->getDOFVector(i)->getBoundaryManager())
rhs_->getDOFVector(i)->getBoundaryManager()->
@@ -1018,5 +1008,61 @@ namespace AMDiS {
solution_->deserialize(in);
}
+
+ void ProblemVec::computeError(AdaptInfo *adaptInfo)
+ {
+ FUNCNAME("ProblemVec::computeError()");
+
+ for (int i = 0; i < nComponents; i++) {
+ TEST_EXIT(exactSolutionFcts[i])("No solution function given!\n");
+
+ // Compute the difference between exact and computed solution
+ DOFVector<double> *tmp = NEW DOFVector<double>(componentSpaces[i], "tmp");
+ tmp->interpol(exactSolutionFcts[i]);
+ double solMax = tmp->absMax();
+ *tmp -= *(solution_->getDOFVector(i));
+
+ MSG("L2 error = %.8e\n", tmp->L2Norm());
+ MSG("L-inf error = %.8e\n", tmp->absMax() / solMax);
+
+ adaptInfo->setEstSum(tmp->absMax() / solMax, i);
+ adaptInfo->setEstMax(tmp->absMax() / solMax, i);
+
+ // To set element estimates, compute a vector with the difference
+ // between exact and computed solution for each DOF.
+ DOFVector<double> *sol = NEW DOFVector<double>(componentSpaces[i], "tmp");
+ sol->interpol(exactSolutionFcts[i]);
+ DOFVector<double>::Iterator it1(sol, USED_DOFS);
+ DOFVector<double>::Iterator it2(tmp, USED_DOFS);
+ for (it1.reset(), it2.reset(); !it1.end(); ++it1, ++it2) {
+ if ((abs(*it1) <= DBL_TOL) || (abs(*it2) <= DBL_TOL)) {
+ *it2 = 0.0;
+ } else {
+ *it2 = abs(*it2 / *it1);
+ }
+ }
+
+ // Compute estimate for every mesh element
+ Vector<DegreeOfFreedom> locInd(componentSpaces[i]->getBasisFcts()->getNumber());
+ TraverseStack stack;
+ ElInfo *elInfo = stack.traverseFirst(componentMeshes[i], -1, Mesh::CALL_LEAF_EL);
+ while (elInfo) {
+ componentSpaces[i]->getBasisFcts()->getLocalIndicesVec(elInfo->getElement(),
+ componentSpaces[i]->getAdmin(),
+ &locInd);
+ double estimate = 0.0;
+ for (int j = 0; j < componentSpaces[i]->getBasisFcts()->getNumber(); j++) {
+ estimate += (*tmp)[locInd[j]];
+ }
+ elInfo->getElement()->setEstimation(estimate, i);
+ elInfo->getElement()->setMark(0);
+
+ elInfo = stack.traverseNext(elInfo);
+ }
+
+ DELETE tmp;
+ DELETE sol;
+ }
+ }
}
diff --git a/AMDiS/src/ProblemVec.h b/AMDiS/src/ProblemVec.h
index 39a27af8..6dca6598 100644
--- a/AMDiS/src/ProblemVec.h
+++ b/AMDiS/src/ProblemVec.h
@@ -322,9 +322,9 @@ namespace AMDiS {
*/
inline Mesh* getMesh(int comp) {
FUNCNAME("ProblemVec::getMesh()");
- TEST_EXIT(comp < static_cast<int>(componentMeshes_.size()) && comp >= 0)
+ TEST_EXIT(comp < static_cast<int>(componentMeshes.size()) && comp >= 0)
("invalid component number\n");
- return componentMeshes_[comp];
+ return componentMeshes[comp];
};
/** \brief
@@ -523,6 +523,13 @@ namespace AMDiS {
return fileWriters_;
};
+ protected:
+ /** \brief
+ * If the exact solution is known, the problem can compute the exact
+ * error instead of the error estimation. This is done in this function.
+ */
+ void computeError(AdaptInfo *adaptInfo);
+
protected:
/** \brief
* Name of this problem.
@@ -552,7 +559,7 @@ namespace AMDiS {
/** \brief
* Pointer to the meshes for the different problem components
*/
- std::vector<Mesh*> componentMeshes_;
+ std::vector<Mesh*> componentMeshes;
/** \brief
* Responsible for element marking.
--
GitLab