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

* Error computation

parent 0aca9309
...@@ -568,7 +568,12 @@ namespace AMDiS { ...@@ -568,7 +568,12 @@ namespace AMDiS {
/** \brief /** \brief
* Returns maximum of DOFVector * Returns maximum of DOFVector
*/ */
T max() const; T max() const;
/** \brief
* Returns absolute maximum of DOFVector
*/
T absMax() const;
/** \brief /** \brief
* Used by interpol while mesh traversal * Used by interpol while mesh traversal
......
#include <list> #include <list>
#include <algorithm>
#include "FixVec.h" #include "FixVec.h"
#include "Boundary.h" #include "Boundary.h"
...@@ -277,6 +278,12 @@ namespace AMDiS { ...@@ -277,6 +278,12 @@ namespace AMDiS {
return m; return m;
} }
template<typename T>
T DOFVector<T>::absMax() const
{
return std::max(abs(max()), abs(min()));
}
template<typename T> template<typename T>
void gemv(MatrixTranspose transpose, T alpha, void gemv(MatrixTranspose transpose, T alpha,
const DOFMatrix& a, const DOFVector<T>& x, const DOFMatrix& a, const DOFVector<T>& x,
......
...@@ -552,12 +552,13 @@ namespace AMDiS { ...@@ -552,12 +552,13 @@ namespace AMDiS {
for (int i = 0; i < nComponents; i++) { for (int i = 0; i < nComponents; i++) {
DOFVector<double> *tmp = NEW DOFVector<double>(componentSpaces[i], "tmp"); DOFVector<double> *tmp = NEW DOFVector<double>(componentSpaces[i], "tmp");
tmp->interpol(exactSolutionFcts[i]); 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)); *tmp -= *(solution_->getDOFVector(i));
double l2Error = tmp->L2Norm(); double l2Error = tmp->L2Norm();
double maxError = tmp->max() / t; double maxError = tmp->absMax() / t;
MSG("L2 error = %.8e\n", l2Error); MSG("L2 error = %.8e\n", l2Error);
MSG("Max error = %.8e\n", maxError); MSG("L-inf error = %.8e\n", maxError);
DELETE tmp; DELETE tmp;
} }
} else { } else {
......
...@@ -4,35 +4,33 @@ RecoveryStructure& RecoveryStructure::operator=(const RecoveryStructure& rhs) { ...@@ -4,35 +4,33 @@ RecoveryStructure& RecoveryStructure::operator=(const RecoveryStructure& rhs) {
if (rhs.coords) { if (rhs.coords) {
if (!coords) if (!coords) {
coords=NEW WorldVector<double>; coords = NEW WorldVector<double>;
*coords=*rhs.coords; }
} *coords = *rhs.coords;
else { } else {
if (coords) { if (coords) {
DELETE coords; DELETE coords;
coords=NULL; coords = NULL;
} }
} }
if (rhs.A) { if (rhs.A) {
if (!A) if (!A)
A=NEW Matrix<double>(rhs.A->getNumRows(), rhs.A->getNumCols()); A = NEW Matrix<double>(rhs.A->getNumRows(), rhs.A->getNumCols());
*A=*rhs.A; *A = *rhs.A;
} } else {
else {
if (A) { if (A) {
DELETE A; DELETE A;
A=NULL; A = NULL;
} }
} }
if (rhs.rec_uh) { if (rhs.rec_uh) {
if (!rec_uh) if (!rec_uh)
rec_uh=NEW Vector<double>(rhs.rec_uh->getSize()); rec_uh = NEW Vector<double>(rhs.rec_uh->getSize());
*rec_uh=*rhs.rec_uh; *rec_uh = *rhs.rec_uh;
} } else {
else {
if (rec_uh) { if (rec_uh) {
DELETE rec_uh; DELETE rec_uh;
rec_uh=NULL; rec_uh=NULL;
...@@ -42,24 +40,22 @@ RecoveryStructure& RecoveryStructure::operator=(const RecoveryStructure& rhs) { ...@@ -42,24 +40,22 @@ RecoveryStructure& RecoveryStructure::operator=(const RecoveryStructure& rhs) {
if (rhs.rec_grdUh) { if (rhs.rec_grdUh) {
if (!rec_grdUh) if (!rec_grdUh)
rec_grdUh = NEW Vector<WorldVector<double> >(rhs.rec_grdUh->getSize()); rec_grdUh = NEW Vector<WorldVector<double> >(rhs.rec_grdUh->getSize());
*rec_grdUh=*rhs.rec_grdUh; *rec_grdUh = *rhs.rec_grdUh;
} } else {
else {
if (rec_grdUh) { if (rec_grdUh) {
DELETE rec_grdUh; DELETE rec_grdUh;
rec_grdUh=NULL; rec_grdUh = NULL;
} }
} }
if (rhs.neighbors) { if (rhs.neighbors) {
if (!neighbors) if (!neighbors)
neighbors = NEW std::set<DegreeOfFreedom>; neighbors = NEW std::set<DegreeOfFreedom>;
*neighbors=*rhs.neighbors ; *neighbors = *rhs.neighbors ;
} } else {
else {
if (neighbors) { if (neighbors) {
DELETE neighbors; DELETE neighbors;
neighbors=NULL; neighbors = NULL;
} }
} }
...@@ -903,7 +899,7 @@ Recovery::recovery(DOFVector<double> *uh, const FiniteElemSpace *fe_space, ...@@ -903,7 +899,7 @@ Recovery::recovery(DOFVector<double> *uh, const FiniteElemSpace *fe_space,
AbstractFunction<double, double> *f_scal, AbstractFunction<double, double> *f_scal,
DOFVector<double> *aux_vec) DOFVector<double> *aux_vec)
{ {
FUNCNAME("Recovery::recovery"); FUNCNAME("Recovery::recovery()");
clear(); clear();
...@@ -916,55 +912,52 @@ Recovery::recovery(DOFVector<double> *uh, const FiniteElemSpace *fe_space, ...@@ -916,55 +912,52 @@ Recovery::recovery(DOFVector<double> *uh, const FiniteElemSpace *fe_space,
DOFVector<RecoveryStructure>::Iterator SV_it(struct_vec, USED_DOFS); DOFVector<RecoveryStructure>::Iterator SV_it(struct_vec, USED_DOFS);
// Solving local systems. // Solving local systems.
for (SV_it.reset(); !SV_it.end(); ++SV_it) for (SV_it.reset(); !SV_it.end(); ++SV_it) {
{ if ((*SV_it).A) {
if ((*SV_it).A) { int error = Cholesky::solve((*SV_it).A,
int error = Cholesky::solve((*SV_it).A, (*SV_it).rec_grdUh,
(*SV_it).rec_grdUh, (*SV_it).rec_grdUh);
(*SV_it).rec_grdUh); TEST_EXIT_DBG(error)
TEST_EXIT_DBG(error) ("There must be some error, matrix is not positive definite.\n");
("There must be some error, matrix is not positive definite.\n");
}
} }
}
// define result vector // define result vector
static DOFVector<WorldVector<double> > *vec = NULL; static DOFVector<WorldVector<double> > *vec = NULL;
DOFVector<WorldVector<double> > *result = NULL; DOFVector<WorldVector<double> > *result = NULL;
// Allocate memory for result vector // Allocate memory for result vector
if (vec && vec->getFESpace() != feSpace) if (vec && vec->getFESpace() != feSpace) {
{ DELETE vec;
DELETE vec; vec = NULL;
vec = NULL; }
}
if (!vec) if (!vec) {
vec = NEW DOFVector<WorldVector<double> >(feSpace, "gradient"); vec = NEW DOFVector<WorldVector<double> >(feSpace, "gradient");
}
result = vec; result = vec;
result->set(WorldVector<double>(DEFAULT_VALUE, 0.0)); result->set(WorldVector<double>(DEFAULT_VALUE, 0.0));
DOFVector<WorldVector<double> >::Iterator grdIt(result, USED_DOFS); DOFVector<WorldVector<double> >::Iterator grdIt(result, USED_DOFS);
std::set<DegreeOfFreedom>::const_iterator setIterator; std::set<DegreeOfFreedom>::const_iterator setIterator;
int i;
for (SV_it.reset(), grdIt.reset(); !grdIt.end(); ++SV_it, ++grdIt) for (SV_it.reset(), grdIt.reset(); !grdIt.end(); ++SV_it, ++grdIt) {
{ if ((*SV_it).rec_grdUh) {
if ((*SV_it).rec_grdUh) *grdIt = (*(*SV_it).rec_grdUh)[0];
*grdIt = (*(*SV_it).rec_grdUh)[0]; } else {
else for (setIterator = (*SV_it).neighbors->begin();
{ setIterator != (*SV_it).neighbors->end();
for (setIterator = (*SV_it).neighbors->begin(); ++setIterator) {
setIterator != (*SV_it).neighbors->end(); for (int i = 0; i < n_monomials; i++)
++setIterator) *grdIt = *grdIt + (*(*struct_vec)[*setIterator].rec_grdUh)[i] *
{ (*(*matrix_fcts)[0][i])(*(*SV_it).coords,
for (i=0; i<n_monomials; i++) *(*struct_vec)[*setIterator].coords);
*grdIt = *grdIt + (*(*struct_vec)[*setIterator].rec_grdUh)[i] * }
(*(*matrix_fcts)[0][i])(*(*SV_it).coords, *grdIt = *grdIt * (1.0 / (*SV_it).neighbors->size());
*(*struct_vec)[*setIterator].coords);
}
*grdIt = *grdIt * (1.0/(*SV_it).neighbors->size());
}
} }
}
return result; return result;
} }
...@@ -975,7 +968,7 @@ Recovery::recovery(DOFVector<double> *uh, ...@@ -975,7 +968,7 @@ Recovery::recovery(DOFVector<double> *uh,
AbstractFunction<double, double> *f_scal, AbstractFunction<double, double> *f_scal,
DOFVector<double> *aux_vec) DOFVector<double> *aux_vec)
{ {
FUNCNAME("Recovery::simpleAveraging"); FUNCNAME("Recovery::simpleAveraging()");
TEST_EXIT_DBG(!(f_vec && f_scal))("Only one diffusion function, please!\n"); TEST_EXIT_DBG(!(f_vec && f_scal))("Only one diffusion function, please!\n");
...@@ -983,16 +976,16 @@ Recovery::recovery(DOFVector<double> *uh, ...@@ -983,16 +976,16 @@ Recovery::recovery(DOFVector<double> *uh,
// define result vector // define result vector
static DOFVector<WorldVector<double> > *vec = NULL; static DOFVector<WorldVector<double> > *vec = NULL;
DOFVector<WorldVector<double> > *result = NULL; DOFVector<WorldVector<double> > *result = NULL;
// Allocate memory for result vector // Allocate memory for result vector
if (vec && vec->getFESpace() != fe_space) if (vec && vec->getFESpace() != fe_space) {
{ DELETE vec;
DELETE vec; vec = NULL;
vec = NULL; }
} if (!vec) {
if (!vec)
vec = NEW DOFVector<WorldVector<double> >(fe_space, "gradient"); vec = NEW DOFVector<WorldVector<double> >(fe_space, "gradient");
}
result = vec; result = vec;
result->set(WorldVector<double>(DEFAULT_VALUE, 0.0)); result->set(WorldVector<double>(DEFAULT_VALUE, 0.0));
...@@ -1000,13 +993,11 @@ Recovery::recovery(DOFVector<double> *uh, ...@@ -1000,13 +993,11 @@ Recovery::recovery(DOFVector<double> *uh,
DOFVector<double> volume(fe_space, "volume"); DOFVector<double> volume(fe_space, "volume");
volume.set(0.0); volume.set(0.0);
int i;
Mesh *mesh = fe_space->getMesh(); Mesh *mesh = fe_space->getMesh();
int dim = mesh->getDim(); int dim = mesh->getDim();
const BasisFunction *basFcts = fe_space->getBasisFcts(); const BasisFunction *basFcts = fe_space->getBasisFcts();
DOFAdmin *admin = fe_space->getAdmin(); DOFAdmin *admin = fe_space->getAdmin();
int numPreDOFs = admin->getNumberOfPreDOFs(0); int numPreDOFs = admin->getNumberOfPreDOFs(0);
...@@ -1044,7 +1035,7 @@ Recovery::recovery(DOFVector<double> *uh, ...@@ -1044,7 +1035,7 @@ Recovery::recovery(DOFVector<double> *uh,
fAtBary = (*f_scal)(fAtBary); fAtBary = (*f_scal)(fAtBary);
} }
for (i = 0; i < dim + 1; i++) { for (int i = 0; i < dim + 1; i++) {
DegreeOfFreedom dofIndex = dof[i][numPreDOFs]; DegreeOfFreedom dofIndex = dof[i][numPreDOFs];
(*result)[dofIndex] += grd * fAtBary * det; (*result)[dofIndex] += grd * fAtBary * det;
volume[dofIndex] += det; volume[dofIndex] += det;
...@@ -1065,7 +1056,7 @@ Recovery::recovery(DOFVector<double> *uh, ...@@ -1065,7 +1056,7 @@ Recovery::recovery(DOFVector<double> *uh,
void Recovery::test(DOFVector<double> *uh, const FiniteElemSpace *fe_space) void Recovery::test(DOFVector<double> *uh, const FiniteElemSpace *fe_space)
{ {
FUNCNAME("Recovery::test"); FUNCNAME("Recovery::test()");
clear(); clear();
...@@ -1078,13 +1069,12 @@ void Recovery::test(DOFVector<double> *uh, const FiniteElemSpace *fe_space) ...@@ -1078,13 +1069,12 @@ void Recovery::test(DOFVector<double> *uh, const FiniteElemSpace *fe_space)
WorldVector<double> coord; WorldVector<double> coord;
// for every DOFs // for every DOFs
for (LM_iterator.reset(); !LM_iterator.end(); ++LM_iterator) for (LM_iterator.reset(); !LM_iterator.end(); ++LM_iterator) {
{ position = LM_iterator.getDOFIndex();
position = LM_iterator.getDOFIndex(); MSG("Node: ");
MSG("Node: "); std::cout << position << std::endl;
std::cout << position << std::endl; (*struct_vec)[position].print();
(*struct_vec)[position].print(); }
}
return; return;
} }
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment