From 5695b4ddca9f59db744bba427efd91d2b9b0de05 Mon Sep 17 00:00:00 2001
From: Simon Praetorius <simon.praetorius@tu-dresden.de>
Date: Wed, 15 Jun 2011 07:52:22 +0000
Subject: [PATCH] correction of ellipt demo
---
demo/src/ellipt.cc | 304 ++-------------------------------------------
1 file changed, 11 insertions(+), 293 deletions(-)
diff --git a/demo/src/ellipt.cc b/demo/src/ellipt.cc
index a5ad31d1..8c1704a2 100644
--- a/demo/src/ellipt.cc
+++ b/demo/src/ellipt.cc
@@ -1,6 +1,4 @@
#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 =====
--
GitLab