Skip to content

GitLab

Commit 5695b4dd authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

correction of ellipt demo

parent 6374f084
Hide whitespace changes
Inline Side-by-side
demo/src/ellipt.cc
View file @ 5695b4dd
#include "AMDiS.h"
#include "ProblemStatDbg.h"
#include "Debug.h"
using namespace AMDiS;
using namespace std;
......@@ -9,284 +7,6 @@ using namespace std;
// ===== function definitions ================================================
// ===========================================================================
class MyProblemVec : public ProblemStatDbg
{
public:
MyProblemVec(std::string nameStr, ProblemIterationInterface *problemIteration = NULL)
: ProblemStatDbg(nameStr,problemIteration) {}
SolverMatrix<Matrix<DOFMatrix*> >* getSolverMatrix() { return &solverMatrix; }
void setSolverMatrix(SolverMatrix<Matrix<DOFMatrix*> >& solverMatrix_) { solverMatrix= solverMatrix_; }
inline DOFVector<double>* getRhsVector(int i = 0)
{
return rhs->getDOFVector(i);
}
};
class MyIteration : public StandardProblemIteration
{
public:
MyIteration(MyProblemVec *prob_) : StandardProblemIteration(prob_), prob(prob_) {}
~MyIteration() {}
Flag oneIteration(AdaptInfo *adaptInfo, Flag toDo = FULL_ITERATION)
{ FUNCNAME("PIteration::oneIteration()");
Flag flag, markFlag;
if (toDo.isSet(MARK))
markFlag = prob->markElements(adaptInfo);
else
markFlag = 3;
// refine
if (toDo.isSet(ADAPT) && markFlag.isSet(MESH_REFINED))
flag= prob->refineMesh(adaptInfo);
// coarsen
if (toDo.isSet(ADAPT) && markFlag.isSet(MESH_COARSENED))
flag|= prob->coarsenMesh(adaptInfo);
if (toDo.isSet(SOLVE)) {
prob->buildAfterCoarsen(adaptInfo, markFlag, true, true);
MSG("apply Dirichlet BC...\n");
bool inside= false;
WorldVector<double> p1;
p1[0]= 1.0;
p1[1]= -1.0;
for (int i= 0; i<prob->getSolution()->getSize(); ++i) {
DegreeOfFreedom idx1= getDOFidxAtCoords(prob->getSolution()->getDOFVector(i), p1, inside);
if(inside) {
MSG("set DBC to row %d...\n",idx1);
applyDbc(prob->getSystemMatrix(i,i)->getBaseMatrix(), *prob->getRhsVector(i), idx1, 0.0);
} else {
MSG("idx not found!!!\n");
}
}
// prob->writeDbgMatrix("ellipt_mat.dat");
// applySingularPertubation(prob->getSystemMatrix(0,0)->getBaseMatrix());
// applySingularDBC(prob->getSystemMatrix(0,0)->getBaseMatrix());
SolverMatrix<Matrix<DOFMatrix*> > solverMatrix;
solverMatrix.setMatrix(*prob->getSystemMatrix());
prob->setSolverMatrix(solverMatrix);
debug::writeMatlabMatrix(*(prob->getSystemMatrix()), "ellipt_mat.m");
prob->solve(adaptInfo, false);
}
if (toDo.isSet(ESTIMATE)) {
prob->estimate(adaptInfo);
}
return flag;
};
template <typename Matrix>
void applySingularPertubation(Matrix& m)
{
using namespace mtl;
// Type of m's elements
typedef typename mtl::Collection<Matrix>::value_type value_type;
int nnz= m.nnz_local(m.num_rows()-1);
// Create inserter for matrix m
// Existing values are not overwritten but inserted
matrix::inserter<Matrix, update_times<value_type> > ins(m, nnz);
ins[m.num_rows()-1][m.num_rows()-1] << 1.0+1.e-5;
}
template <typename Matrix>
void applySingularDBC(Matrix& m)
{
using namespace mtl;
typename traits::row<Matrix>::type r(m);
typename traits::col<Matrix>::type c(m);
typename traits::value<Matrix>::type v(m);
typedef typename traits::range_generator<tag::row, Matrix>::type c_type;
typedef typename traits::range_generator<tag::nz, c_type>::type ic_type;
// typedef typename Collection<Matrix>::value_type value_type;
for (c_type cursor(begin<tag::row>(m)), cend(end<tag::row>(m)); cursor != cend; ++cursor) {
for (ic_type icursor(begin<tag::nz>(cursor)), icend(end<tag::nz>(cursor)); icursor != icend; ++icursor) {
v(*icursor, (r(*icursor)==c(*icursor)?1.0:0.0));
}
break;
}
}
template <typename Matrix, typename Vector>
void applyDbc(Matrix& m, Vector& vec, DegreeOfFreedom idx, double value)
{
using namespace mtl;
typename Matrix::size_type idx_= idx;
typename traits::row<Matrix>::type r(m);
typename traits::col<Matrix>::type c(m);
typename traits::value<Matrix>::type v(m);
typedef typename traits::range_generator<tag::row, Matrix>::type c_type;
typedef typename traits::range_generator<tag::nz, c_type>::type ic_type;
// typedef typename Collection<Matrix>::value_type value_type;
for (c_type cursor(begin<tag::row>(m)+idx_), cend(begin<tag::row>(m)+idx_+1); cursor != cend; ++cursor) {
for (ic_type icursor(begin<tag::nz>(cursor)), icend(end<tag::nz>(cursor)); icursor != icend; ++icursor) {
v(*icursor, (r(*icursor)==c(*icursor)?1.0:0.0));
if(r(*icursor)==c(*icursor)) {
MSG("DBC set to (%d, %d)\n", r(*icursor), c(*icursor));
}
}
break;
}
vec[idx]= value;
}
DegreeOfFreedom getDOFidxAtCoords(DOFVector<double> *vec, WorldVector<double> &p, bool &inside, ElInfo *oldElInfo=NULL, bool useOldElInfo=false)
{
const FiniteElemSpace *feSpace = vec->getFeSpace();
Mesh *mesh = vec->getFeSpace()->getMesh();
const BasisFunction *basFcts = feSpace->getBasisFcts();
int dim = mesh->getDim();
int numBasFcts = basFcts->getNumber();
DegreeOfFreedom *localIndices = new DegreeOfFreedom[numBasFcts];
DimVec<double> lambda(dim, NO_INIT);
ElInfo *elInfo = mesh->createNewElInfo();
int dofIdx=0;
inside=false;
if(oldElInfo && useOldElInfo && oldElInfo->getMacroElement()) {
inside=mesh->findElInfoAtPoint(p,elInfo,lambda,oldElInfo->getMacroElement(),NULL,NULL);
delete oldElInfo;
} else {
inside=mesh->findElInfoAtPoint(p,elInfo,lambda,NULL,NULL,NULL);
}
if(oldElInfo) oldElInfo=elInfo;
if(inside) {
basFcts->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(), localIndices);
FixVec<WorldVector<double>, VERTEX> coords = elInfo->getMacroElement()->getCoord();
int minIdx=-1;
double minDist=1.e15;
for(int i=0;i<coords.size();++i) {
WorldVector<double> dist = coords[i]-p;
double newDist = AMDiS::norm(dist);
if(newDist<minDist) {
minIdx=i;
minDist=newDist;
}
}
if(minIdx>=0) {
dofIdx = localIndices[minIdx];
}
}
delete [] localIndices;
if(!oldElInfo) delete elInfo;
return dofIdx;
};
protected:
MyProblemVec *prob;
};
/// < factor*phase*grad(u) , grad(psi) >
class Phase_SOT : public SecondOrderTerm
{
public:
/// Constructor.
Phase_SOT(DOFVectorBase<double>* phaseDV_, double fac_=1.0) :
SecondOrderTerm(6), phaseDV(phaseDV_), fac(fac_)
{
TEST_EXIT(phaseDV_)("no vector phase!\n");
auxFeSpaces.insert(phaseDV_->getFeSpace());
};
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
Quadrature *quad = NULL) {
getVectorAtQPs(phaseDV, elInfo, subAssembler, quad, phase);
};
void initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
SubAssembler* subAssembler,
Quadrature *quad = NULL) {
getVectorAtQPs(phaseDV, smallElInfo, largeElInfo, subAssembler, quad, phase);
};
inline void getLALt(const ElInfo *elInfo, std::vector<mtl::dense2D<double> > &LALt) const
{
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
const int nPoints = static_cast<int>(LALt.size());
for (int iq = 0; iq < nPoints; iq++)
l1lt(Lambda, LALt[iq], f(iq) * phase[iq] * fac);
};
inline void eval(int nPoints,
const AMDiS::ElementVector&,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
ElementVector& result,
double opFactor)
{
if (D2UhAtQP) {
for (int iq = 0; iq < nPoints; iq++) {
double feval = f(iq) * phase[iq] * opFactor * fac;
double resultQP = 0.0;
for (int i = 0; i < dimOfWorld; i++)
resultQP += D2UhAtQP[iq][i][i];
result[iq] += feval * resultQP;
}
}
};
void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
std::vector<WorldVector<double> > &result)
{
int nPoints = grdUhAtQP.size();
for (int iq = 0; iq < nPoints; iq++) {
double factor = f(iq) * phase[iq] * fac;
axpy(factor, grdUhAtQP[iq], result[iq]);
}
};
virtual double f(const int iq) const
{
return 1.0;
};
void setFactor(const double fac_) { fac=fac_; }
protected:
DOFVectorBase<double> *phaseDV;
mtl::dense_vector<double> phase;
double fac;
};
/// simple circle: phi_0(x,y) := sqrt(x^2+y^2)-r
class Circle : public AbstractFunction<double, WorldVector<double> >
{
public:
Circle(double radius_, WorldVector<double> midPoint_) :
radius(radius_), midPoint(midPoint_) {}
Circle(double radius_) :
radius(radius_) { midPoint.set(0.0); }
double operator()(const WorldVector<double>& x) const
{
double result = 0.0;
for(int k=0; k<x.getSize(); k++) {
result += sqr(x[k]-midPoint[k]);
}
result = sqrt(result)-radius;
return result;
};
private:
double radius;
WorldVector<double> midPoint;
};
/// Dirichlet boundary function
class G : public AbstractFunction<double, WorldVector<double> >
{
......@@ -295,7 +15,7 @@ public:
/// Implementation of AbstractFunction::operator().
double operator()(const WorldVector<double>& x) const
{
return 0.0;
return exp(-10.0 * (x * x));
}
};
......@@ -309,7 +29,10 @@ public:
/// Implementation of AbstractFunction::operator().
double operator()(const WorldVector<double>& x) const
{
return sin(M_PI*x[0]);
int dow = Global::getGeo(WORLD);
double r2 = (x * x);
double ux = exp(-10.0 * r2);
return -(400.0 * r2 - 20.0 * dow) * ux;
}
};
......@@ -325,30 +48,27 @@ int main(int argc, char* argv[])
TEST_EXIT(argc >= 2)("usage: ellipt initfile\n");
// ===== init parameters =====
// ===== init parameters =====
Parameters::init(argv[1]);
// ===== create and init the scalar problem =====
MyProblemVec ellipt("ellipt");
ProblemStat ellipt("ellipt");
ellipt.initialize(INIT_ALL);
MyIteration elIter(&ellipt);
// === create adapt info ===
AdaptInfo adaptInfo("ellipt->adapt", 2);
AdaptInfo adaptInfo("ellipt->adapt");
// === create adapt ===
AdaptStationary adapt("ellipt->adapt", elIter, adaptInfo);
AdaptStationary adapt("ellipt->adapt", ellipt, adaptInfo);
// DOFVector<double> phase(ellipt.getFeSpace(), "phase");
// phase.interpol(new Circle(0.25));
// ===== create matrix operator =====
Operator matrixOperator(ellipt.getFeSpace());
matrixOperator.addSecondOrderTerm(new Simple_SOT);
ellipt.addMatrixOperator(matrixOperator, 0, 0);
ellipt.addMatrixOperator(matrixOperator, 1, 1);
// ===== create rhs operator =====
......@@ -356,12 +76,10 @@ int main(int argc, char* argv[])
Operator rhsOperator(ellipt.getFeSpace());
rhsOperator.addZeroOrderTerm(new CoordsAtQP_ZOT(new F(degree)));
ellipt.addVectorOperator(rhsOperator, 0);
ellipt.addVectorOperator(rhsOperator, 1);
// ===== add boundary conditions =====
ellipt.addDirichletBC(3, 0, 0, new G);
ellipt.addDirichletBC(5, 1, 1, new G);
ellipt.addDirichletBC(1, 0, 0, new G);
// ===== start adaption loop =====
......
Markdown is supported
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