From b66281563813691317fbc473324dc98bbddd577d Mon Sep 17 00:00:00 2001
From: Thomas Witkowski <thomas.witkowski@gmx.de>
Date: Fri, 29 Aug 2008 09:09:52 +0000
Subject: [PATCH] * Error computation
---
AMDiS/src/DOFVector.h | 7 +-
AMDiS/src/DOFVector.hh | 7 ++
AMDiS/src/ProblemVec.cc | 9 +--
AMDiS/src/Recovery.cc | 148 +++++++++++++++++++---------------------
4 files changed, 87 insertions(+), 84 deletions(-)
diff --git a/AMDiS/src/DOFVector.h b/AMDiS/src/DOFVector.h
index aad09fe7..572887ab 100644
--- a/AMDiS/src/DOFVector.h
+++ b/AMDiS/src/DOFVector.h
@@ -568,7 +568,12 @@ namespace AMDiS {
/** \brief
* Returns maximum of DOFVector
*/
- T max() const;
+ T max() const;
+
+ /** \brief
+ * Returns absolute maximum of DOFVector
+ */
+ T absMax() const;
/** \brief
* Used by interpol while mesh traversal
diff --git a/AMDiS/src/DOFVector.hh b/AMDiS/src/DOFVector.hh
index 5361890c..93a2bb43 100644
--- a/AMDiS/src/DOFVector.hh
+++ b/AMDiS/src/DOFVector.hh
@@ -1,4 +1,5 @@
#include <list>
+#include <algorithm>
#include "FixVec.h"
#include "Boundary.h"
@@ -277,6 +278,12 @@ namespace AMDiS {
return m;
}
+ template<typename T>
+ T DOFVector<T>::absMax() const
+ {
+ return std::max(abs(max()), abs(min()));
+ }
+
template<typename T>
void gemv(MatrixTranspose transpose, T alpha,
const DOFMatrix& a, const DOFVector<T>& x,
diff --git a/AMDiS/src/ProblemVec.cc b/AMDiS/src/ProblemVec.cc
index f0be8cca..36104a09 100644
--- a/AMDiS/src/ProblemVec.cc
+++ b/AMDiS/src/ProblemVec.cc
@@ -552,12 +552,13 @@ namespace AMDiS {
for (int i = 0; i < nComponents; i++) {
DOFVector<double> *tmp = NEW DOFVector<double>(componentSpaces[i], "tmp");
tmp->interpol(exactSolutionFcts[i]);
- double t = tmp->max();
+ // double t = max(abs(tmp->max()), abs(tmp->min()));
+ double t = tmp->absMax();
*tmp -= *(solution_->getDOFVector(i));
double l2Error = tmp->L2Norm();
- double maxError = tmp->max() / t;
- MSG("L2 error = %.8e\n", l2Error);
- MSG("Max error = %.8e\n", maxError);
+ double maxError = tmp->absMax() / t;
+ MSG("L2 error = %.8e\n", l2Error);
+ MSG("L-inf error = %.8e\n", maxError);
DELETE tmp;
}
} else {
diff --git a/AMDiS/src/Recovery.cc b/AMDiS/src/Recovery.cc
index cdadc23d..b3f69172 100644
--- a/AMDiS/src/Recovery.cc
+++ b/AMDiS/src/Recovery.cc
@@ -4,35 +4,33 @@ RecoveryStructure& RecoveryStructure::operator=(const RecoveryStructure& rhs) {
if (rhs.coords) {
- if (!coords)
- coords=NEW WorldVector<double>;
- *coords=*rhs.coords;
- }
- else {
+ if (!coords) {
+ coords = NEW WorldVector<double>;
+ }
+ *coords = *rhs.coords;
+ } else {
if (coords) {
DELETE coords;
- coords=NULL;
+ coords = NULL;
}
}
if (rhs.A) {
if (!A)
- A=NEW Matrix<double>(rhs.A->getNumRows(), rhs.A->getNumCols());
- *A=*rhs.A;
- }
- else {
+ A = NEW Matrix<double>(rhs.A->getNumRows(), rhs.A->getNumCols());
+ *A = *rhs.A;
+ } else {
if (A) {
DELETE A;
- A=NULL;
+ A = NULL;
}
}
if (rhs.rec_uh) {
if (!rec_uh)
- rec_uh=NEW Vector<double>(rhs.rec_uh->getSize());
- *rec_uh=*rhs.rec_uh;
- }
- else {
+ rec_uh = NEW Vector<double>(rhs.rec_uh->getSize());
+ *rec_uh = *rhs.rec_uh;
+ } else {
if (rec_uh) {
DELETE rec_uh;
rec_uh=NULL;
@@ -42,24 +40,22 @@ RecoveryStructure& RecoveryStructure::operator=(const RecoveryStructure& rhs) {
if (rhs.rec_grdUh) {
if (!rec_grdUh)
rec_grdUh = NEW Vector<WorldVector<double> >(rhs.rec_grdUh->getSize());
- *rec_grdUh=*rhs.rec_grdUh;
- }
- else {
+ *rec_grdUh = *rhs.rec_grdUh;
+ } else {
if (rec_grdUh) {
DELETE rec_grdUh;
- rec_grdUh=NULL;
+ rec_grdUh = NULL;
}
}
if (rhs.neighbors) {
if (!neighbors)
neighbors = NEW std::set<DegreeOfFreedom>;
- *neighbors=*rhs.neighbors ;
- }
- else {
+ *neighbors = *rhs.neighbors ;
+ } else {
if (neighbors) {
DELETE neighbors;
- neighbors=NULL;
+ neighbors = NULL;
}
}
@@ -903,7 +899,7 @@ Recovery::recovery(DOFVector<double> *uh, const FiniteElemSpace *fe_space,
AbstractFunction<double, double> *f_scal,
DOFVector<double> *aux_vec)
{
- FUNCNAME("Recovery::recovery");
+ FUNCNAME("Recovery::recovery()");
clear();
@@ -916,55 +912,52 @@ Recovery::recovery(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) {
- int error = Cholesky::solve((*SV_it).A,
- (*SV_it).rec_grdUh,
- (*SV_it).rec_grdUh);
- TEST_EXIT_DBG(error)
- ("There must be some error, matrix is not positive definite.\n");
- }
+ for (SV_it.reset(); !SV_it.end(); ++SV_it) {
+ if ((*SV_it).A) {
+ int error = Cholesky::solve((*SV_it).A,
+ (*SV_it).rec_grdUh,
+ (*SV_it).rec_grdUh);
+ TEST_EXIT_DBG(error)
+ ("There must be some error, matrix is not positive definite.\n");
}
+ }
// define result vector
static DOFVector<WorldVector<double> > *vec = NULL;
- DOFVector<WorldVector<double> > *result = NULL;
+ DOFVector<WorldVector<double> > *result = NULL;
// Allocate memory for result vector
- if (vec && vec->getFESpace() != feSpace)
- {
- DELETE vec;
- vec = NULL;
- }
- if (!vec)
+ if (vec && vec->getFESpace() != feSpace) {
+ DELETE vec;
+ vec = NULL;
+ }
+
+ if (!vec) {
vec = NEW DOFVector<WorldVector<double> >(feSpace, "gradient");
+ }
+
result = vec;
result->set(WorldVector<double>(DEFAULT_VALUE, 0.0));
DOFVector<WorldVector<double> >::Iterator grdIt(result, USED_DOFS);
std::set<DegreeOfFreedom>::const_iterator setIterator;
- int i;
- for (SV_it.reset(), grdIt.reset(); !grdIt.end(); ++SV_it, ++grdIt)
- {
- if ((*SV_it).rec_grdUh)
- *grdIt = (*(*SV_it).rec_grdUh)[0];
- else
- {
- for (setIterator = (*SV_it).neighbors->begin();
- setIterator != (*SV_it).neighbors->end();
- ++setIterator)
- {
- for (i=0; i<n_monomials; i++)
- *grdIt = *grdIt + (*(*struct_vec)[*setIterator].rec_grdUh)[i] *
- (*(*matrix_fcts)[0][i])(*(*SV_it).coords,
- *(*struct_vec)[*setIterator].coords);
- }
- *grdIt = *grdIt * (1.0/(*SV_it).neighbors->size());
- }
+ for (SV_it.reset(), grdIt.reset(); !grdIt.end(); ++SV_it, ++grdIt) {
+ if ((*SV_it).rec_grdUh) {
+ *grdIt = (*(*SV_it).rec_grdUh)[0];
+ } else {
+ for (setIterator = (*SV_it).neighbors->begin();
+ setIterator != (*SV_it).neighbors->end();
+ ++setIterator) {
+ for (int i = 0; i < n_monomials; i++)
+ *grdIt = *grdIt + (*(*struct_vec)[*setIterator].rec_grdUh)[i] *
+ (*(*matrix_fcts)[0][i])(*(*SV_it).coords,
+ *(*struct_vec)[*setIterator].coords);
+ }
+ *grdIt = *grdIt * (1.0 / (*SV_it).neighbors->size());
}
+ }
return result;
}
@@ -975,7 +968,7 @@ Recovery::recovery(DOFVector<double> *uh,
AbstractFunction<double, double> *f_scal,
DOFVector<double> *aux_vec)
{
- FUNCNAME("Recovery::simpleAveraging");
+ FUNCNAME("Recovery::simpleAveraging()");
TEST_EXIT_DBG(!(f_vec && f_scal))("Only one diffusion function, please!\n");
@@ -983,16 +976,16 @@ Recovery::recovery(DOFVector<double> *uh,
// define result vector
static DOFVector<WorldVector<double> > *vec = NULL;
- DOFVector<WorldVector<double> > *result = NULL;
+ DOFVector<WorldVector<double> > *result = NULL;
// Allocate memory for result vector
- if (vec && vec->getFESpace() != fe_space)
- {
- DELETE vec;
- vec = NULL;
- }
- if (!vec)
+ if (vec && vec->getFESpace() != fe_space) {
+ DELETE vec;
+ vec = NULL;
+ }
+ if (!vec) {
vec = NEW DOFVector<WorldVector<double> >(fe_space, "gradient");
+ }
result = vec;
result->set(WorldVector<double>(DEFAULT_VALUE, 0.0));
@@ -1000,13 +993,11 @@ Recovery::recovery(DOFVector<double> *uh,
DOFVector<double> volume(fe_space, "volume");
volume.set(0.0);
- int i;
-
Mesh *mesh = fe_space->getMesh();
- int dim = mesh->getDim();
+ int dim = mesh->getDim();
const BasisFunction *basFcts = fe_space->getBasisFcts();
- DOFAdmin *admin = fe_space->getAdmin();
+ DOFAdmin *admin = fe_space->getAdmin();
int numPreDOFs = admin->getNumberOfPreDOFs(0);
@@ -1044,7 +1035,7 @@ Recovery::recovery(DOFVector<double> *uh,
fAtBary = (*f_scal)(fAtBary);
}
- for (i = 0; i < dim + 1; i++) {
+ for (int i = 0; i < dim + 1; i++) {
DegreeOfFreedom dofIndex = dof[i][numPreDOFs];
(*result)[dofIndex] += grd * fAtBary * det;
volume[dofIndex] += det;
@@ -1065,7 +1056,7 @@ Recovery::recovery(DOFVector<double> *uh,
void Recovery::test(DOFVector<double> *uh, const FiniteElemSpace *fe_space)
{
- FUNCNAME("Recovery::test");
+ FUNCNAME("Recovery::test()");
clear();
@@ -1078,13 +1069,12 @@ void Recovery::test(DOFVector<double> *uh, const FiniteElemSpace *fe_space)
WorldVector<double> coord;
// for every DOFs
- for (LM_iterator.reset(); !LM_iterator.end(); ++LM_iterator)
- {
- position = LM_iterator.getDOFIndex();
- MSG("Node: ");
- std::cout << position << std::endl;
- (*struct_vec)[position].print();
- }
+ for (LM_iterator.reset(); !LM_iterator.end(); ++LM_iterator) {
+ position = LM_iterator.getDOFIndex();
+ MSG("Node: ");
+ std::cout << position << std::endl;
+ (*struct_vec)[position].print();
+ }
return;
}
--
GitLab