Commit 02d0cc3a authored by Naumann, Andreas's avatar Naumann, Andreas

CreatorMap: pollutionError added

AdaptInfo: test-cout added (could be deleted)
ConditionalEstimator: should run over all elements, but there's an error (adapts only one times)
parent 50c66bbc
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
......@@ -195,8 +195,11 @@ namespace AMDiS {
{
int size = scalContents.getSize();
for (int i = 0; i < size; i++)
{
std::cout<<"est_sum:"<<scalContents[i]->est_sum<<" spaceTol:"<<scalContents[i]->spaceTolerance<<std::endl;
if (!(scalContents[i]->est_sum < scalContents[i]->spaceTolerance))
return false;
}
return true;
}
......
......@@ -27,17 +27,17 @@ namespace AMDiS {
(elInfo->getElement()->getElementData(PARTITION_ED));
TEST_EXIT_DBG(elInfo->getElement()->isLeaf())("element not leaf\n");
TEST_EXIT_DBG(elData)("no partitoin data on leaf element %d (rank %d)\n",
TEST_EXIT_DBG(elData)("no partition data on leaf element %d (rank %d)\n",
elInfo->getElement()->getIndex(),
MPI::COMM_WORLD.Get_rank());
PartitionStatus status = elData->getPartitionStatus();
if(status == IN || status == OVERLAP) {
//if(status == IN || status == OVERLAP) {
decoratedEstimator_->estimateElement(elInfo);
} else {
elInfo->getElement()->setMark(0);
}
//} else {
// elInfo->getElement()->setMark(0);
//}
elInfo = stack.traverseNext(elInfo);
}
......@@ -46,10 +46,10 @@ namespace AMDiS {
PartitionElementData *elData = dynamic_cast<PartitionElementData*>
(elInfo->getElement()->getElementData(PARTITION_ED));
PartitionStatus status = elData->getPartitionStatus();
if(status == IN) {
//if(status == IN) {
elementCount_++;
partition_sum += elInfo->getElement()->getEstimation(row_);
}
//}
elInfo = stack.traverseNext(elInfo);
}
......
......@@ -8,6 +8,7 @@
#include "Estimator.h"
#include "RecoveryEstimator.h"
#include "ResidualEstimator.h"
#include "pollutionError_konvex.h"
#include "LeafData.h"
#include "SurfaceRegion_ED.h"
#include "DOFMatrix.h"
......@@ -99,6 +100,9 @@ namespace AMDiS {
creator = new RecoveryEstimator::Creator;
addCreator("recovery", creator);
creator = new PollutionErrorKonvex::Creator;
addCreator("pollution2d",creator);
}
template<>
......
......@@ -1355,7 +1355,7 @@ namespace AMDiS {
#ifdef DEBUG
//test mesh
TraverseStack testStack;
/*TraverseStack testStack;
ElInfo* testElInfo = testStack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
while(testElInfo)
{
......@@ -1367,7 +1367,7 @@ namespace AMDiS {
std::cout<<"none";
std::cout<<"\t"<<testElement->getIndex()<<"\t"<<testElInfo->getLevel()<<std::endl;
testElInfo=testStack.traverseNext(testElInfo);
}
}*/
//end test
#endif
......
#include "pollutionError_konvex.h"
#include "mpi.h"
namespace AMDiS
{
bool AMDiS::PollutionErrorKonvex::BoundaryTraverseStack::proofConditions(ElInfo *elInfo)
{
if(elInfo==NULL) return false;
Element *curEl=elInfo->getElement();
PartitionElementData *ped=static_cast<PartitionElementData*>(curEl->getElementData(PARTITION_ED));
TEST_EXIT_DBG(ped!=NULL)("No PartitionElementData found");
if(ped!=NULL && ped->getPartitionStatus()==OVERLAP)
{
//check neighbour
int nrNeigh=curEl->getGeo(NEIGH);
for(int i=0;i<nrNeigh;i++)
{
Element* curNeigh=elInfo->getNeighbour(i);
if(curNeigh!=NULL)
{
ped=static_cast<PartitionElementData*>(curNeigh->getElementData(PARTITION_ED));
if(ped!=NULL && ped->getPartitionStatus()==ps)
{
return true;
}
}
}
}
return false;
}
ElInfo * AMDiS::PollutionErrorKonvex::BoundaryTraverseStack::traverseFirst(Mesh *mesh,PartitionStatus ps_)
{
FUNCNAME("PollutionErrorKonvex::BoundaryTraverseStack::traverseFirst");
ElInfo *swap=TraverseStack::traverseFirst(mesh,-1,
Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_NEIGH);
ps=ps_;
if(swap!=NULL && !proofConditions(swap))
swap=traverseNext(swap);
return swap;
}
ElInfo * AMDiS::PollutionErrorKonvex::BoundaryTraverseStack::traverseNext(ElInfo *old)
{
FUNCNAME("PollutionErrorKonvex::BoundaryTraverseStack::traverseNext");
ElInfo *swap=TraverseStack::traverseNext(old);
while(swap!=NULL && !proofConditions(swap))
{
swap=TraverseStack::traverseNext(swap);
}
return swap;
}
PollutionErrorKonvex::PollutionErrorKonvex(std::string name, int r):ResidualEstimator(name,r)
{
FUNCNAME("PollutionErrorKonvexEstimator::PollutionErrorKonvexEstimator()");
GET_PARAMETER(0,name+"->C1p","%f",&C1p);
GET_PARAMETER(0,name+"->C2p","%f",&C2p);
/*C1=sqrt(C1);
C2=sqrt(C2);*/
}
void PollutionErrorKonvex::init(double timestep)
{
FUNCNAME("PollutionErrorKonvex::init");
ResidualEstimator::init(timestep);
//set the distance (ie. width of the overlapped area)
//approximate the overlapped area as rectangle. then set d=max(edge1, edge2)
//retrieve all elements at the boundary overlapped -> out
BoundaryTraverseStack* bts=new BoundaryTraverseStack();
ElInfo* curEl=bts->traverseFirst(mesh,OUT);
FixVec<WorldVector<double>, VERTEX> &fv=curEl->getCoords();
//std::cout<<"got fixvec"<<std::endl;
WorldVector<double>* koordinaten=fv.begin();
double xomin,xomax;
xomin=xomax=(*koordinaten)[0];
double yomin,yomax;
yomin=yomax=(*koordinaten)[1];
for(;koordinaten!=fv.end();koordinaten++)
{
xomin=min(xomin,(*koordinaten)[0]);
xomax=max(xomax,(*koordinaten)[0]);
yomin=min(yomin,(*koordinaten)[1]);
yomax=max(yomax,(*koordinaten)[1]);
}
curEl=bts->traverseNext(curEl);
//find xmax,xmin,ymax and ymin
while(curEl!=NULL)
{
fv=curEl->getCoords();
for(koordinaten=fv.begin();koordinaten!=fv.end();koordinaten++)
{
xomin=min(xomin,(*koordinaten)[0]);
xomax=max(xomax,(*koordinaten)[0]);
yomin=min(yomin,(*koordinaten)[1]);
yomax=max(yomax,(*koordinaten)[1]);
}
curEl=bts->traverseNext(curEl);
}
//retrieve all elements at the boundary overlapped -> in
curEl=bts->traverseFirst(mesh,OUT);
fv=curEl->getCoords();
koordinaten=fv.begin();
double ximin,ximax;
ximin=ximax=(*koordinaten)[0];
double yimin,yimax;
yimin=yimax=(*koordinaten)[1];
for(;koordinaten!=fv.end();koordinaten++)
{
ximin=min(ximin,(*koordinaten)[0]);
ximax=max(ximax,(*koordinaten)[0]);
yimin=min(yimin,(*koordinaten)[1]);
yimax=max(yimax,(*koordinaten)[1]);
}
curEl=bts->traverseNext(curEl);
//find xmax,xmin,ymax and ymin
while(curEl!=NULL)
{
fv=curEl->getCoords();
for(koordinaten=fv.begin();koordinaten!=fv.end();koordinaten++)
{
ximin=min(ximin,(*koordinaten)[0]);
ximax=max(ximax,(*koordinaten)[0]);
yimin=min(yimin,(*koordinaten)[1]);
yimax=max(yimax,(*koordinaten)[1]);
}
curEl=bts->traverseNext(curEl);
}
delete bts;
double m1=max(abs(ximin-xomax),abs(xomin-ximax));
double m2=max(abs(yimin-yomax),abs(yomin-yimax));
d=min(m1,m2);
#if DEBUG
std::cout<<"distance in polError:"<<d<<std::endl;
#endif
C2p=C2p/d;
}
void PollutionErrorKonvex::estimateElement(ElInfo *elInfo)
{
FUNCNAME("PollutionErrorKonvex::estimateElement()");
//switch between inner and outer
//inner part --> use H_1Norm
//outer part --> use L_2Norm
PartitionElementData* ped=static_cast<PartitionElementData*>(elInfo->getElement()->getElementData(PARTITION_ED));
TEST_EXIT_DBG(ped!=NULL)("No partition Elementdata found");
PartitionStatus ps=ped->getPartitionStatus();
TEST_EXIT_DBG(ps!=UNDEFINED)("Elementstatus undefined");
//save old values
double C0,C1;
C0=ResidualEstimator::C0;
C1=ResidualEstimator::C1;
if(ps==OUT)
{
//outer part
ResidualEstimator::norm=L2_NORM;
ResidualEstimator::C0 *=C2p;
ResidualEstimator::C1 *=C2p;
}else
{
//inner part
ResidualEstimator::norm=H1_NORM;
ResidualEstimator::C0 *=C1p;
ResidualEstimator::C1 *=C1p;
}
ResidualEstimator::estimateElement(elInfo);
//set old values
ResidualEstimator::C0=C0;
ResidualEstimator::C1=C1;
}
// void exit(bool output=true)
// {
//
// }
}
#include "pollutionError_konvex.h"
#include "mpi.h"
namespace AMDiS
{
bool AMDiS::PollutionErrorKonvex::BoundaryTraverseStack::proofConditions(ElInfo *elInfo)
{
if(elInfo==NULL) return false;
Element *curEl=elInfo->getElement();
PartitionElementData *ped=static_cast<PartitionElementData*>(curEl->getElementData(PARTITION_ED));
TEST_EXIT_DBG(ped!=NULL)("No PartitionElementData found");
if(ped!=NULL && ped->getPartitionStatus()==OVERLAP)
{
//check neighbour
int nrNeigh=curEl->getGeo(NEIGH);
for(int i=0;i<nrNeigh;i++)
{
Element* curNeigh=elInfo->getNeighbour(i);
if(curNeigh!=NULL)
{
ped=static_cast<PartitionElementData*>(curNeigh->getElementData(PARTITION_ED));
if(ped!=NULL && ped->getPartitionStatus()==ps)
{
return true;
}
}
}
}
return false;
}
ElInfo * AMDiS::PollutionErrorKonvex::BoundaryTraverseStack::traverseFirst(Mesh *mesh,PartitionStatus ps_)
{
FUNCNAME("PollutionErrorKonvex::BoundaryTraverseStack::traverseFirst");
ElInfo *swap=TraverseStack::traverseFirst(mesh,-1,
Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_NEIGH);
ps=ps_;
if(swap!=NULL && !proofConditions(swap))
swap=traverseNext(swap);
return swap;
}
ElInfo * AMDiS::PollutionErrorKonvex::BoundaryTraverseStack::traverseNext(ElInfo *old)
{
FUNCNAME("PollutionErrorKonvex::BoundaryTraverseStack::traverseNext");
ElInfo *swap=TraverseStack::traverseNext(old);
while(swap!=NULL && !proofConditions(swap))
{
swap=TraverseStack::traverseNext(swap);
}
return swap;
}
PollutionErrorKonvex::PollutionErrorKonvex(std::string name, int r):ResidualEstimator(name,r)
{
FUNCNAME("PollutionErrorKonvexEstimator::PollutionErrorKonvexEstimator()");
GET_PARAMETER(0,name+"->C1p","%f",&C1p);
GET_PARAMETER(0,name+"->C2p","%f",&C2p);
/*C1=sqrt(C1);
C2=sqrt(C2);*/
}
void PollutionErrorKonvex::init(double timestep)
{
FUNCNAME("PollutionErrorKonvex::init");
ResidualEstimator::init(timestep);
//set the distance (ie. width of the overlapped area)
//approximate the overlapped area as rectangle. then set d=max(edge1, edge2)
//retrieve all elements at the boundary overlapped -> out
BoundaryTraverseStack* bts=new BoundaryTraverseStack();
ElInfo* curEl=bts->traverseFirst(mesh,OUT);
FixVec<WorldVector<double>, VERTEX> &fv=curEl->getCoords();
//std::cout<<"got fixvec"<<std::endl;
WorldVector<double>* koordinaten=fv.begin();
double xomin,xomax;
xomin=xomax=(*koordinaten)[0];
double yomin,yomax;
yomin=yomax=(*koordinaten)[1];
for(;koordinaten!=fv.end();koordinaten++)
{
xomin=min(xomin,(*koordinaten)[0]);
xomax=max(xomax,(*koordinaten)[0]);
yomin=min(yomin,(*koordinaten)[1]);
yomax=max(yomax,(*koordinaten)[1]);
}
curEl=bts->traverseNext(curEl);
//find xmax,xmin,ymax and ymin
while(curEl!=NULL)
{
fv=curEl->getCoords();
for(koordinaten=fv.begin();koordinaten!=fv.end();koordinaten++)
{
xomin=min(xomin,(*koordinaten)[0]);
xomax=max(xomax,(*koordinaten)[0]);
yomin=min(yomin,(*koordinaten)[1]);
yomax=max(yomax,(*koordinaten)[1]);
}
curEl=bts->traverseNext(curEl);
}
//retrieve all elements at the boundary overlapped -> in
curEl=bts->traverseFirst(mesh,OUT);
fv=curEl->getCoords();
koordinaten=fv.begin();
double ximin,ximax;
ximin=ximax=(*koordinaten)[0];
double yimin,yimax;
yimin=yimax=(*koordinaten)[1];
for(;koordinaten!=fv.end();koordinaten++)
{
ximin=min(ximin,(*koordinaten)[0]);
ximax=max(ximax,(*koordinaten)[0]);
yimin=min(yimin,(*koordinaten)[1]);
yimax=max(yimax,(*koordinaten)[1]);
}
curEl=bts->traverseNext(curEl);
//find xmax,xmin,ymax and ymin
while(curEl!=NULL)
{
fv=curEl->getCoords();
for(koordinaten=fv.begin();koordinaten!=fv.end();koordinaten++)
{
ximin=min(ximin,(*koordinaten)[0]);
ximax=max(ximax,(*koordinaten)[0]);
yimin=min(yimin,(*koordinaten)[1]);
yimax=max(yimax,(*koordinaten)[1]);
}
curEl=bts->traverseNext(curEl);
}
delete bts;
double m1=max(abs(ximin-xomax),abs(xomin-ximax));
double m2=max(abs(yimin-yomax),abs(yomin-yimax));
d=min(m1,m2);
#if DEBUG
std::cout<<"distance:"<<d<<std::endl;
#endif
C2p=C2p/d;
}
void PollutionErrorKonvex::estimateElement(ElInfo *elInfo)
{
FUNCNAME("PollutionErrorKonvex::estimateElement()");
//switch between inner and outer
//inner part --> use H_1Norm
//outer part --> use L_2Norm
PartitionElementData* ped=static_cast<PartitionElementData*>(elInfo->getElement()->getElementData(PARTITION_ED));
TEST_EXIT_DBG(ped!=NULL)("No partition Elementdata found");
PartitionStatus ps=ped->getPartitionStatus();
TEST_EXIT_DBG(ps!=UNDEFINED)("Elementstatus undefined");
//save old values
double C0,C1;
C0=ResidualEstimator::C0;
C1=ResidualEstimator::C1;
if(ps==OUT)
{
//outer part
ResidualEstimator::norm=L2_NORM;
ResidualEstimator::C0 *=C2p;
ResidualEstimator::C1 *=C2p;
}else
{
//inner part
ResidualEstimator::norm=H1_NORM;
ResidualEstimator::C0 *=C1p;
ResidualEstimator::C1 *=C1p;
}
ResidualEstimator::estimateElement(elInfo);
//set old values
ResidualEstimator::C0=C0;
ResidualEstimator::C1=C1;
}
// void exit(bool output=true)
// {
//
// }
}
#include "pollutionError_konvex.h"
#include "mpi.h"
namespace AMDiS
{
bool AMDiS::PollutionErrorKonvex::BoundaryTraverseStack::proofConditions(ElInfo *elInfo)
{
if(elInfo==NULL) return false;
Element *curEl=elInfo->getElement();
PartitionElementData *ped=static_cast<PartitionElementData*>(curEl->getElementData(PARTITION_ED));
TEST_EXIT_DBG(ped!=NULL)("No PartitionElementData found");
if(ped!=NULL && ped->getPartitionStatus()==OVERLAP)
{
//check neighbour
int nrNeigh=curEl->getGeo(NEIGH);
std::cout<<"nrNeigh:"<<nrNeigh<<std::endl;
for(int i=0;i<nrNeigh;i++)
{
std::cout <<"i:"<<i<<std::endl;
Element* curNeigh=elInfo->getNeighbour(i);
if(curNeigh!=NULL)
{
std::cout<<"curNeigh:"<<curNeigh<<std::endl;
ped=static_cast<PartitionElementData*>(curNeigh->getElementData(PARTITION_ED));
std::cout<<"ped:"<<ped<<std::endl;
if(ped!=NULL && ped->getPartitionStatus()==ps)
return true;
}
}
}
std::cout<<"return false"<<std::endl;
return false;
}
ElInfo * AMDiS::PollutionErrorKonvex::BoundaryTraverseStack::traverseFirst(Mesh *mesh,PartitionStatus ps_)
{
FUNCNAME("PollutionErrorKonvex::BoundaryTraverseStack::traverseFirst");
ElInfo *swap=TraverseStack::traverseFirst(mesh,-1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
ps=ps_;
if(swap!=NULL && !proofConditions(swap))
swap=traverseNext(swap);
return swap;
}
ElInfo * AMDiS::PollutionErrorKonvex::BoundaryTraverseStack::traverseNext(ElInfo *old)
{
FUNCNAME("PollutionErrorKonvex::BoundaryTraverseStack::traverseNext");
ElInfo *swap=TraverseStack::traverseNext(old);
while(swap!=NULL && !proofConditions(swap))
{
swap=TraverseStack::traverseNext(swap);
}
return swap;
}
PollutionErrorKonvex::PollutionErrorKonvex(std::string name, int r):ResidualEstimator(name,r)
{
FUNCNAME("PollutionErrorKonvexEstimator::PollutionErrorKonvexEstimator()");
GET_PARAMETER(0,name+"->C1","%f",&C1);
GET_PARAMETER(0,name+"->C2r","%f",&C2);
/*C1=sqrt(C1);
C2=sqrt(C2);*/
}
void PollutionErrorKonvex::init(double timestep)
{
FUNCNAME("PollutionErrorKonvex::init");
ResidualEstimator::init(timestep);
//set the distance (ie. width of the overlapped area)
//approximate the overlapped area as rectangle. then set d=max(edge1, edge2)
//retrieve all elements at the boundary overlapped -> out
BoundaryTraverseStack* bts=new BoundaryTraverseStack();
ElInfo* curEl=bts->traverseFirst(mesh,OUT);
FixVec<WorldVector<double>, VERTEX> &fv=curEl->getCoords();
//std::cout<<"got fixvec"<<std::endl;
WorldVector<double>* koordinaten=fv.begin();
double xomin,xomax;
xomin=xomax=(*koordinaten)[0];
double yomin,yomax;
yomin=yomax=(*koordinaten)[1];
for(;koordinaten!=fv.end();koordinaten++)
{
xomin=min(xomin,(*koordinaten)[0]);
xomax=max(xomax,(*koordinaten)[0]);
yomin=min(yomin,(*koordinaten)[1]);
yomax=max(yomax,(*koordinaten)[1]);
}
curEl=bts->traverseNext(curEl);
//find xmax,xmin,ymax and ymin
while(curEl!=NULL)
{
fv=curEl->getCoords();
for(koordinaten=fv.begin();koordinaten!=fv.end();koordinaten++)
{
xomin=min(xomin,(*koordinaten)[0]);
xomax=max(xomax,(*koordinaten)[0]);
yomin=min(yomin,(*koordinaten)[1]);
yomax=max(yomax,(*koordinaten)[1]);
}
curEl=bts->traverseNext(curEl);
}
//retrieve all elements at the boundary overlapped -> in
curEl=bts->traverseFirst(mesh,OUT);
fv=curEl->getCoords();
koordinaten=fv.begin();
double ximin,ximax;
ximin=ximax=(*koordinaten)[0];
double yimin,yimax;
yimin=yimax=(*koordinaten)[1];
for(;koordinaten!=fv.end();koordinaten++)
{
ximin=min(ximin,(*koordinaten)[0]);
ximax=max(ximax,(*koordinaten)[0]);
yimin=min(yimin,(*koordinaten)[1]);
yimax=max(yimax,(*koordinaten)[1]);
}
curEl=bts->traverseNext(curEl);
//find xmax,xmin,ymax and ymin
while(curEl!=NULL)
{
fv=curEl->getCoords();
for(koordinaten=fv.begin();koordinaten!=fv.end();koordinaten++)
{
ximin=min(ximin,(*koordinaten)[0]);
ximax=max(ximax,(*koordinaten)[0]);
yimin=min(yimin,(*koordinaten)[1]);
yimax=max(yimax,(*koordinaten)[1]);
}
curEl=bts->traverseNext(curEl);
}
delete bts;
double m1=max(abs(ximin-xomax),abs(xomin-ximax));
double m2=max(abs(yimin-yomax),abs(yomin-yimax));
std::cout<<"m1:"<<m1<<" m2:"<<m2<<std::endl;
d=min(m1,m2);
std::cout<<"d:"<<d<<std::endl;
d=0.0000001;
C2=C2/d;
}
void PollutionErrorKonvex::estimateElement(ElInfo *elInfo)
{
FUNCNAME("PollutionErrorKonvex::estimateElement()");
//switch between inner and outer
//inner part --> use H_1Norm
//outer part --> use L_2Norm
PartitionElementData* ped=static_cast<PartitionElementData*>(elInfo->getElement()->getElementData(PARTITION_ED));
TEST_EXIT_DBG(ped!=NULL)("No partition Elementdata found");
PartitionStatus ps=ped->getPartitionStatus();
TEST_EXIT_DBG(ps!=UNDEFINED)("Elementstatus undefined");
if(ps==OUT)
{
//outer part
ResidualEstimator::norm=L2_NORM;
ResidualEstimator::C0=ResidualEstimator::C1=C2;
}else
{
//inner part
ResidualEstimator::norm=H1_NORM;
ResidualEstimator::C0=ResidualEstimator::C1=C1;
}
ResidualEstimator::estimateElement(elInfo);
}
// void exit(bool output=true)
// {
//
// }
}
/** \file pollutionError_konvex.h */