Commit d906fc98 authored by Thomas Witkowski's avatar Thomas Witkowski

Code refactoring, so added some no bugs.

parent a125580e
......@@ -10,57 +10,44 @@
namespace AMDiS {
CoarseningManager* CoarseningManager::traversePtr = NULL;
/****************************************************************************/
/* sets the mark on all elements that have to be coarsend */
/****************************************************************************/
int CoarseningManager::coarsenMarkFunction(ElInfo *el_info)
{
el_info->getElement()->setMark(traversePtr->globalMark);
return 0;
}
/****************************************************************************/
/* tries to coarsen every element of mesh at least mark times */
/****************************************************************************/
Flag CoarseningManager::globalCoarsen(Mesh *aMesh, int mark)
{
if (mark >= 0) return(0);
globalMark = mark;
traversePtr = this;
aMesh->traverse(-1, Mesh::CALL_LEAF_EL, coarsenMarkFunction);
return(coarsenMesh(aMesh));
}
int CoarseningManager::spreadCoarsenMarkFunction(ElInfo* el_info)
{
Element *el = el_info->getElement();
signed char mark;
if (el->getChild(0)) {
/****************************************************************************/
/* interior node of the tree */
/****************************************************************************/
mark = max(el->getChild(0)->getMark(), el->getChild(1)->getMark());
el->setMark(std::min(mark + 1, 0));
} else {
/****************************************************************************/
/* leaf node of the tree */
/****************************************************************************/
if (el->getMark() < 0)
el->setMark(el->getMark() - 1);
if (mark >= 0)
return 0;
TraverseStack stack;
ElInfo *elInfo = stack.traverseFirst(aMesh, -1, Mesh::CALL_LEAF_EL);
while (elInfo) {
elInfo->getElement()->setMark(mark);
elInfo = stack.traverseNext(elInfo);
}
return 0;
return coarsenMesh(aMesh);
}
void CoarseningManager::spreadCoarsenMark()
{
traversePtr = this;
mesh->traverse(-1, Mesh::CALL_EVERY_EL_POSTORDER, spreadCoarsenMarkFunction);
TraverseStack stack;
ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_POSTORDER);
while (elInfo) {
Element *el = elInfo->getElement();
if (el->getChild(0)) {
// interior node of the tree
signed char mark = max(el->getChild(0)->getMark(), el->getChild(1)->getMark());
el->setMark(std::min(mark + 1, 0));
} else {
// leaf node of the tree
if (el->getMark() < 0)
el->setMark(el->getMark() - 1);
}
elInfo = stack.traverseNext(elInfo);
}
}
......@@ -69,17 +56,15 @@ namespace AMDiS {
/* resets the element marks */
/****************************************************************************/
int CoarseningManager::cleanUpAfterCoarsenFunction(ElInfo *el_info)
{
Element *el = el_info->getElement();
el->setMark(max(el->getMark(), 0));
return 0;
}
void CoarseningManager::cleanUpAfterCoarsen()
{
traversePtr = this;
mesh->traverse(-1, Mesh::CALL_LEAF_EL, cleanUpAfterCoarsenFunction);
TraverseStack stack;
ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
while (elInfo) {
Element *el = elInfo->getElement();
el->setMark(max(el->getMark(), 0));
elInfo = stack.traverseNext(elInfo);
}
}
/****************************************************************************/
......
......@@ -42,7 +42,6 @@ namespace AMDiS {
CoarseningManager()
: mesh(NULL),
stack(NULL),
globalMark(0),
doMore(0)
{}
......@@ -92,15 +91,6 @@ namespace AMDiS {
*/
void spreadCoarsenMark();
/// Used while traversal in spreadCoarsenMark
static int spreadCoarsenMarkFunction(ElInfo* el_info);
/// Used while traversal in cleanUpAfterCoarsen
static int cleanUpAfterCoarsenFunction(ElInfo* el_info);
/// Sets the mark on all elements that have to be coarsend
static int coarsenMarkFunction(ElInfo *el_info);
/// Resets the element marks
void cleanUpAfterCoarsen();
......@@ -111,18 +101,9 @@ namespace AMDiS {
/// Used for non recursive mesh traversal.
TraverseStack *stack;
/// Used by globalCoarsen to remember the given mark value
int globalMark;
/// Spezifies whether the coarsening operation is still in progress
bool doMore;
/** \brief
* Used while mesh traversal to have a pointer to this CoarseningManager
* from a static method
*/
static CoarseningManager* traversePtr;
/// Spezifies how many DOFVectors should restricted while coarsening
int callCoarseRestrict;
......
......@@ -68,13 +68,6 @@ namespace AMDiS {
double* max,
bool writeLeafData = false,
int comp = 0);
// methods for traversal
static int maxErrAtQpFct(ElInfo* elinfo);
static int H1ErrFct(ElInfo* elinfo);
static int L2ErrFct(ElInfo* elinfo);
static int relFct(ElInfo* elinfo);
public:
static T errUFct(const DimVec<double>& lambda);
......@@ -114,13 +107,6 @@ namespace AMDiS {
static const AbstractFunction<WorldVector<T>, WorldVector<double> >* pGrdU;
static const BasisFunction* basFct;
static const DOFVector<T>* errUh;
static double maxErr;
static double l2Err2;
static double l2Norm2;
static int relative;
static double relNorm2;
static double h1Err2;
static double h1Norm2;
static bool writeInLeafData;
static int component;
};
......@@ -131,13 +117,6 @@ namespace AMDiS {
template<typename T> const AbstractFunction<WorldVector<T>, WorldVector<double> >* Error<T>::pGrdU = NULL;
template<typename T> const BasisFunction* Error<T>::basFct = NULL;
template<typename T> const DOFVector<T>* Error<T>::errUh = NULL;
template<typename T> double Error<T>::maxErr = 0.0;
template<typename T> double Error<T>::l2Err2 = 0.0;
template<typename T> double Error<T>::l2Norm2 = 0.0;
template<typename T> int Error<T>::relative = 0;
template<typename T> double Error<T>::relNorm2 = 0.0;
template<typename T> double Error<T>::h1Err2 = 0.0;
template<typename T> double Error<T>::h1Norm2 = 0.0;
template<typename T> typename Error<T>::AbstrFctErrU Error<T>::errU;
template<typename T> typename Error<T>::AbstrFctGrdErrU Error<T>::grdErrU;
template<typename T> bool Error<T>::writeInLeafData = false;
......
#include "Mesh.h"
#include "Parametric.h"
#include "Quadrature.h"
#include "Traverse.h"
namespace AMDiS {
......@@ -20,150 +21,80 @@ namespace AMDiS {
return ((*pGrdU)(x));
}
template<typename T>
int Error<T>::maxErrAtQpFct(ElInfo* el_info)
{
double err = 0.0;
const double *u_vec, *uh_vec;
elinfo = el_info;
u_vec = quadFast->getQuadrature()->fAtQp(errU, NULL);
uh_vec = errUh->getVecAtQPs(el_info,
NULL,
quadFast,
NULL);
int numPoints = quadFast->getNumPoints();
for (int i = 0; i < numPoints; i++) {
err = u_vec[i] > uh_vec[i] ? u_vec[i] - uh_vec[i] : uh_vec[i] - u_vec[i];
maxErr = max(maxErr, err);
}
return;
}
template<typename T>
double Error<T>::maxErrAtQp(const AbstractFunction<T, WorldVector<double> >& u,
const DOFVector<T>& uh,
const Quadrature* q)
{
FUNCNAME("Error<T>::maxErrAtQp()");
const FiniteElemSpace *fe_space;
const FiniteElemSpace *fe_space;
if (!(pU = &u)) {
ERROR("no function u specified; doing nothing\n");
return(-1.0);
}
if (!(errUh = &uh) || !(fe_space = uh->getFESpace()))
{
ERROR("no discrete function or no fe_space for it; doing nothing\n");
return(-1.0);
}
if (!(basFct = fe_space->getBasisFcts()))
{
ERROR("no basis functions at discrete solution ; doing nothing\n");
return(-1.0);
}
int dim = fe_space->getMesh()->getDim();
if (!(errUh = &uh) || !(fe_space = uh->getFESpace())) {
ERROR("no discrete function or no fe_space for it; doing nothing\n");
return(-1.0);
}
if (!(basFct = fe_space->getBasisFcts())) {
ERROR("no basis functions at discrete solution ; doing nothing\n");
return(-1.0);
}
if (!q)
q = Quadrature::provideQuadrature(dim,
2*fe_space->getBasisFcts()->getDegree()
- 2);
quadFast = FastQuadrature::provideFastQuadrature(basFct, *q, INIT_PHI);
q = Quadrature::provideQuadrature(fe_space->getMesh()->getDim(),
2 * fe_space->getBasisFcts()->getDegree() - 2);
maxErr = 0.0;
fe_space->getMesh()->traverse(-1,
Mesh::FILL_COORDS |
Mesh::CALL_LEAF_EL,
maxErrAtQpFct);
return(maxErr);
}
template<typename T>
int Error<T>::H1ErrFct(ElInfo* el_info)
{
int i, j;
double err, err_2, h1_err_el, norm_el, norm2, det, exact;
const WorldVector<double> *grdu_vec, *grduh_vec;
elinfo = el_info;
grdu_vec = quadFast->getQuadrature()->grdFAtQp(grdErrU, NULL);
det = el_info->getDet();
grduh_vec = errUh->getGrdAtQPs(elinfo, NULL, quadFast, NULL);
int numPoints = quadFast->getNumPoints();
int dow = Global::getGeo(WORLD);
for (h1_err_el = i = 0; i < numPoints; i++)
{
for (err_2 = j = 0; j < dow; j++)
{
err = grdu_vec[i][j] - grduh_vec[i][j];
err_2 += sqr(err);
}
h1_err_el += quadFast->getWeight(i)*err_2;
quadFast = FastQuadrature::provideFastQuadrature(basFct, *q, INIT_PHI);
double maxErr = 0.0;
TraverseStack stack;
ElInfo *elInfo = stack.traverseFirst(fe_space->getMesh(), -1,
Mesh::FILL_COORDS | Mesh::CALL_LEAF_EL);
while (elInfo) {
double err = 0.0;
const double *u_vec, *uh_vec;
u_vec = quadFast->getQuadrature()->fAtQp(errU, NULL);
uh_vec = errUh->getVecAtQPs(elInfo, NULL, quadFast, NULL);
int nPoints = quadFast->getNumPoints();
for (int i = 0; i < nPoints; i++) {
err = u_vec[i] > uh_vec[i] ? u_vec[i] - uh_vec[i] : uh_vec[i] - u_vec[i];
maxErr = max(maxErr, err);
}
exact = det*h1_err_el;
h1Err2 += exact;
maxErr = max(maxErr, exact);
if (writeInLeafData)
el_info->getElement()->setEstimation(exact, component);
if (relative) {
for (norm_el = i = 0; i < numPoints; i++) {
for (norm2 = j = 0; j < dow; j++)
norm2 += sqr(grdu_vec[i][j]);
norm_el += quadFast->getWeight(i)*norm2;
}
h1Norm2 += det*norm_el;
elInfo = stack.traverseNext(elInfo);
}
return 0;
return maxErr;
}
template<typename T>
double Error<T>::H1Err(
const AbstractFunction<WorldVector<T>, WorldVector<double> >& grdU,
const DOFVector<T>& uh,
int relErr,
double* max,
double Error<T>::H1Err(const AbstractFunction<WorldVector<T>, WorldVector<double> > &grdU,
const DOFVector<T> &uh,
int relErr,
double* max,
bool writeLeafData,
int comp)
{
FUNCNAME("Error<T>::H1Err");
FUNCNAME("Error<T>::H1Err()");
const FiniteElemSpace *fe_space;
writeInLeafData = writeLeafData;
component = comp;
Quadrature *q = NULL;
pGrdU = &grdU;
errUh = &uh;
if(!(fe_space = uh.getFESpace()))
{
ERROR("no fe_space for uh; doing nothing\n");
return(0.0);
}
if (!(basFct = fe_space->getBasisFcts()))
{
ERROR("no basis functions at discrete solution ; doing nothing\n");
return(0.0);
}
if (!(fe_space = uh.getFESpace())) {
ERROR("no fe_space for uh; doing nothing\n");
return(0.0);
}
if (!(basFct = fe_space->getBasisFcts())) {
ERROR("no basis functions at discrete solution ; doing nothing\n");
return(0.0);
}
int dim = fe_space->getMesh()->getDim();
int deg = grdU.getDegree();
......@@ -174,69 +105,74 @@ namespace AMDiS {
*q,
INIT_GRD_PHI);
relative = relErr;
maxErr = h1Err2 = h1Norm2 = 0.0;
double relative = relErr;
double maxErr = 0.0, h1Err2 = 0.0, h1Norm2 = 0.0;
TraverseStack stack;
ElInfo *elInfo = stack.traverseFirst(fe_space->getMesh(), -1,
Mesh::FILL_COORDS |
Mesh::CALL_LEAF_EL |
Mesh::FILL_DET |
Mesh::FILL_GRD_LAMBDA);
while (elInfo) {
int i, j;
double err, err_2, h1_err_el, norm_el, norm2, det, exact;
const WorldVector<double> *grdu_vec, *grduh_vec;
grdu_vec = quadFast->getQuadrature()->grdFAtQp(grdErrU, NULL);
det = elInfo->getDet();
grduh_vec = errUh->getGrdAtQPs(elinfo, NULL, quadFast, NULL);
int nPoints = quadFast->getNumPoints();
int dow = Global::getGeo(WORLD);
for (h1_err_el = i = 0; i < nPoints; i++) {
for (err_2 = j = 0; j < dow; j++) {
err = grdu_vec[i][j] - grduh_vec[i][j];
err_2 += sqr(err);
}
h1_err_el += quadFast->getWeight(i)*err_2;
}
exact = det*h1_err_el;
h1Err2 += exact;
maxErr = std::max(maxErr, exact);
if (writeInLeafData)
elInfo->getElement()->setEstimation(exact, component);
if (relative) {
for (norm_el = i = 0; i < nPoints; i++) {
for (norm2 = j = 0; j < dow; j++)
norm2 += sqr(grdu_vec[i][j]);
norm_el += quadFast->getWeight(i)*norm2;
}
h1Norm2 += det*norm_el;
}
fe_space->getMesh()->traverse(-1,
Mesh::FILL_COORDS |
Mesh::CALL_LEAF_EL |
Mesh::FILL_DET |
Mesh::FILL_GRD_LAMBDA,
H1ErrFct);
elInfo = stack.traverseNext(elInfo);
}
if (relative) {
relNorm2 = h1Norm2 + 1.e-15;
fe_space->getMesh()->traverse(-1, Mesh::CALL_LEAF_EL, relFct);
double relNorm2 = h1Norm2 + 1.e-15;
elInfo = stack.traverseFirst(fe_space->getMesh(), -1, Mesh::CALL_LEAF_EL);
while (elInfo) {
double exact = elInfo->getElement()->getEstimation(component) / relNorm2;
if (writeInLeafData)
elinfo->getElement()->setEstimation(exact, component);
elInfo = stack.traverseNext(elInfo);
}
h1Err2 /= relNorm2;
maxErr /= relNorm2;
}
if (max) *max = maxErr;
return(sqrt(h1Err2));
}
template<typename T>
int Error<T>::L2ErrFct(ElInfo* el_info)
{
int i;
double err, det, l2_err_el, norm_el, exact;
const double *u_vec, *uh_vec;
elinfo = el_info;
u_vec = quadFast->getQuadrature()->fAtQp(errU, NULL);
uh_vec = errUh->getVecAtQPs(el_info,
NULL,
quadFast,
NULL);
det = el_info->getDet();
int numPoints = quadFast->getNumPoints();
for (l2_err_el = i = 0; i < numPoints; i++) {
err = u_vec[i] - uh_vec[i];
l2_err_el += quadFast->getWeight(i)*sqr(err);
}
exact = det*l2_err_el;
l2Err2 += exact;
maxErr = max(maxErr, exact);
if (writeInLeafData) {
el_info->getElement()->setEstimation(exact, component);
}
if (relative) {
for (norm_el = i = 0; i < numPoints; i++)
norm_el += quadFast->getWeight(i)*sqr(u_vec[i]);
l2Norm2 += det*norm_el;
}
if (max)
*max = maxErr;
return 0;
return sqrt(h1Err2);
}
template<typename T>
......@@ -278,20 +214,58 @@ namespace AMDiS {
q = Quadrature::provideQuadrature(dim, degree);
quadFast = FastQuadrature::provideFastQuadrature(basFct, *q, INIT_PHI);
relative = relErr;
maxErr = l2Err2 = l2Norm2 = 0.0;
double relative = relErr;
double maxErr = 0.0, l2Err2 = 0.0, l2Norm2 = 0.0;
TraverseStack stack;
ElInfo *elInfo = stack.traverseFirst(fe_space->getMesh(), -1,
Mesh::FILL_COORDS |
Mesh::CALL_LEAF_EL |
Mesh::FILL_DET |
Mesh::FILL_GRD_LAMBDA);
while (elInfo) {
int i;
double err, det, l2_err_el, norm_el, exact;
const double *u_vec, *uh_vec;
u_vec = quadFast->getQuadrature()->fAtQp(errU, NULL);
uh_vec = errUh->getVecAtQPs(elInfo, NULL, quadFast, NULL);
det = elInfo->getDet();
int nPoints = quadFast->getNumPoints();
for (l2_err_el = i = 0; i < nPoints; i++) {
err = u_vec[i] - uh_vec[i];
l2_err_el += quadFast->getWeight(i) * sqr(err);
}
exact = det * l2_err_el;
l2Err2 += exact;
maxErr = std::max(maxErr, exact);
if (writeInLeafData)
elInfo->getElement()->setEstimation(exact, component);
if (relative) {
for (norm_el = i = 0; i < nPoints; i++)
norm_el += quadFast->getWeight(i) * sqr(u_vec[i]);
l2Norm2 += det * norm_el;
}
fe_space->getMesh()->traverse(-1,
Mesh::FILL_COORDS |
Mesh::FILL_DET |
Mesh::FILL_GRD_LAMBDA |
Mesh::CALL_LEAF_EL,
L2ErrFct);
elInfo = stack.traverseNext(elInfo);
}
if (relative) {
relNorm2 = l2Norm2 + 1.e-15;
fe_space->getMesh()->traverse(-1, Mesh::CALL_LEAF_EL, relFct);
double relNorm2 = l2Norm2 + 1.e-15;
elInfo = stack.traverseFirst(fe_space->getMesh(), -1, Mesh::CALL_LEAF_EL);
while (elInfo) {
double exact = elInfo->getElement()->getEstimation(component) / relNorm2;
if (writeInLeafData)
elInfo->getElement()->setEstimation(exact, component);
elInfo = stack.traverseNext(elInfo);
}
l2Err2 /= relNorm2;
}
......@@ -301,14 +275,4 @@ namespace AMDiS {
return (sqrt(l2Err2));
}
template<typename T>
int Error<T>::relFct(ElInfo* elinfo)
{
double exact = elinfo->getElement()->getEstimation(component);
exact /= relNorm2;
if (writeInLeafData)
elinfo->getElement()->setEstimation(exact, component);
return 0;
}
}
......@@ -731,12 +731,12 @@ namespace AMDiS {
if (projector && projector->getType() == VOLUME_PROJECTION) {
mel[i]->setProjection(0, projector);
} else { // boundary projection
for(j = 0; j < mesh->getGeo(NEIGH); j++) {
for (j = 0; j < mesh->getGeo(NEIGH); j++) {
projector = Projection::getProjection((*ind)[j]);
if(projector) {
if (projector) {
mel[i]->setProjection(j, projector);
if(dim > 2) {
for(k = 0; k < numEdgesAtBoundary; k++) {
if (dim > 2) {
for (k = 0; k < numEdgesAtBoundary; k++) {
int edgeNr = Global::getReferenceElement(dim)->getEdgeOfFace(j, k);
mel[i]->setProjection(numFaces + edgeNr, projector);
}
......@@ -1091,16 +1091,14 @@ namespace AMDiS {
DimVec<bool> periodic(dim, DEFAULT_VALUE, false);
if(ed) {
if (ed