Commit 196751d4 authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

* Dirichlet boundary problem for different FE spaces fixed.

parent 3d291c93
...@@ -45,116 +45,6 @@ namespace AMDiS { ...@@ -45,116 +45,6 @@ namespace AMDiS {
*/ */
typedef signed int BoundaryType; typedef signed int BoundaryType;
// /** \ingroup Triangulation
// * \brief
// * Holds information about the type of boundary associated to an edge/face,
// * and how new vertices are projected to the boundary in the case of curved
// * boundaries.
//
// class Boundary
// {
// public:
// MEMORY_MANAGED(Boundary);
// /** \brief
// * constructor
//
// Boundary(BoundaryType type=0) {
// bound = type;
// };
// /** \brief
// * copy constructor
//
// Boundary(const Boundary& old) { bound = old.getBound(); };
// /** \brief
// * destructor
//
// virtual ~Boundary() {};
// /** \brief
// * assignment operator
//
// Boundary& operator=(const Boundary& old) {
// if (this!=&old) bound = old.getBound();
// return *this;
// };
// /** \brief
// * Returns
// * -true: if a new vertex should be projected to a curved boundary
// * -false: otherwise
//
// virtual bool interpolateBoundary() { return false; };
// /** \brief
// * Projection to the curved boundary
//
// virtual void interpolateBoundary(WorldVector<double>& ) {};
// /** \brief
// * Returns \ref bound
//
// inline const BoundaryType getBound() const { return bound; };
// /** \brief
// * Returns
// * -true: is \ref bound is INTERIOR
// * -false: otherwise
//
// inline const bool isInterior() const {return (bound == INTERIOR);};
// /** \brief
// * Returns
// * -true: is \ref bound is DIRICHLET
// * -false: otherwise
//
// inline const bool isDirichlet() const {return (bound == DIRICHLET);};
// /** \brief
// * Returns
// * -true: is \ref bound is NEUMANN
// * -false: otherwise
//
// inline const bool isNeumann() const {return (bound == NEUMANN);};
// /** \brief
// * Returns
// * -true: is \ref bound is ROBIN
// * -false: otherwise
//
// inline const bool isRobin() const {return (bound == ROBIN);};
// /** \brief
// * Returns the new value of \ref bound with respect to its old value and
// * the value of bound_.
//
// BoundaryType newVal(const BoundaryType bound_);
// /** \brief
// * Returns the Boundary with given type from \ref boundaryMap.
//
// static Boundary* getBoundary(BoundaryType type);
// /** \brief
// * Adds Boundary b to \ref boundaryMap.
//
// //static void addBoundary(Boundary *b);
// protected:
// /** \brief
// * type of this boundary
//
// BoundaryType bound;
// protected:
// /** \brief
// * stl map of all existing boundaries.
//
// static ::std::map<BoundaryType, Boundary*> boundaryMap;
// };
BoundaryType newBound(BoundaryType oldBound, BoundaryType newBound); BoundaryType newBound(BoundaryType oldBound, BoundaryType newBound);
} }
......
...@@ -43,7 +43,7 @@ namespace AMDiS { ...@@ -43,7 +43,7 @@ namespace AMDiS {
* Sub class of BoundaryCondition. Local boundary conditions are filled * Sub class of BoundaryCondition. Local boundary conditions are filled
* while mesh traversal. * while mesh traversal.
*/ */
class BoundaryCondition //: public BoundaryCondition class BoundaryCondition
{ {
public: public:
/** \brief /** \brief
...@@ -121,7 +121,9 @@ namespace AMDiS { ...@@ -121,7 +121,9 @@ namespace AMDiS {
*/ */
virtual double boundResidual(ElInfo *elInfo, virtual double boundResidual(ElInfo *elInfo,
DOFMatrix *matrix, DOFMatrix *matrix,
const DOFVectorBase<double> *dv) { return 0.0; }; const DOFVectorBase<double> *dv) {
return 0.0;
};
/** \brief /** \brief
* Returns whether the condition must be treated as dirichlet condition * Returns whether the condition must be treated as dirichlet condition
......
...@@ -15,8 +15,8 @@ namespace AMDiS { ...@@ -15,8 +15,8 @@ namespace AMDiS {
{ {
double result = 0; double result = 0;
::std::map<BoundaryType, BoundaryCondition*>::iterator it; ::std::map<BoundaryType, BoundaryCondition*>::iterator it;
for(it = localBCs.begin(); it != localBCs.end(); ++it) { for (it = localBCs.begin(); it != localBCs.end(); ++it) {
if((*it).second) if ((*it).second)
result += (*it).second->boundResidual(elInfo, matrix, dv); result += (*it).second->boundResidual(elInfo, matrix, dv);
} }
return result; return result;
...@@ -35,7 +35,7 @@ namespace AMDiS { ...@@ -35,7 +35,7 @@ namespace AMDiS {
::std::map<BoundaryType, BoundaryCondition*>::iterator it; ::std::map<BoundaryType, BoundaryCondition*>::iterator it;
if(localBCs.size() > 0) { if (localBCs.size() > 0) {
// get boundaries of all DOFs // get boundaries of all DOFs
localBound = basisFcts->getBound(elInfo, NULL); localBound = basisFcts->getBound(elInfo, NULL);
...@@ -45,18 +45,18 @@ namespace AMDiS { ...@@ -45,18 +45,18 @@ namespace AMDiS {
admin, NULL); admin, NULL);
// apply non dirichlet boundary conditions // apply non dirichlet boundary conditions
for(it = localBCs.begin(); it != localBCs.end(); ++it) { for (it = localBCs.begin(); it != localBCs.end(); ++it) {
if((*it).second) { if ((*it).second) {
if(!(*it).second->isDirichlet()) { if (!(*it).second->isDirichlet()) {
(*it).second->fillBoundaryCondition(vec, elInfo, dofIndices, localBound, nBasFcts); (*it).second->fillBoundaryCondition(vec, elInfo, dofIndices, localBound, nBasFcts);
} }
} }
} }
// apply dirichlet boundary conditions // apply dirichlet boundary conditions
for(it = localBCs.begin(); it != localBCs.end(); ++it) { for (it = localBCs.begin(); it != localBCs.end(); ++it) {
if((*it).second) { if ((*it).second) {
if((*it).second->isDirichlet()) { if ((*it).second->isDirichlet()) {
(*it).second->fillBoundaryCondition(vec, elInfo, dofIndices, localBound, nBasFcts); (*it).second->fillBoundaryCondition(vec, elInfo, dofIndices, localBound, nBasFcts);
} }
} }
...@@ -97,7 +97,7 @@ namespace AMDiS { ...@@ -97,7 +97,7 @@ namespace AMDiS {
for (it = localBCs.begin(); it != localBCs.end(); ++it) { for (it = localBCs.begin(); it != localBCs.end(); ++it) {
if ((*it).second) { if ((*it).second) {
if ((*it).second->isDirichlet()) { if ((*it).second->isDirichlet()) {
(*it).second->fillBoundaryCondition(mat, elInfo, dofIndices, localBound, nBasFcts); (*it).second->fillBoundaryCondition(mat, elInfo, dofIndices, localBound, nBasFcts);
} }
} }
} }
...@@ -165,16 +165,16 @@ namespace AMDiS { ...@@ -165,16 +165,16 @@ namespace AMDiS {
void BoundaryManager::exitVector(DOFVectorBase<double> *vector) void BoundaryManager::exitVector(DOFVectorBase<double> *vector)
{ {
::std::map<BoundaryType, BoundaryCondition*>::iterator it; ::std::map<BoundaryType, BoundaryCondition*>::iterator it;
for(it = localBCs.begin(); it != localBCs.end(); ++it) { for (it = localBCs.begin(); it != localBCs.end(); ++it) {
if((*it).second) { if ((*it).second) {
if(!(*it).second->isDirichlet()) { if (!(*it).second->isDirichlet()) {
(*it).second->exitVector(vector); (*it).second->exitVector(vector);
} }
} }
} }
for(it = localBCs.begin(); it != localBCs.end(); ++it) { for (it = localBCs.begin(); it != localBCs.end(); ++it) {
if((*it).second) { if ((*it).second) {
if((*it).second->isDirichlet()) { if ((*it).second->isDirichlet()) {
(*it).second->exitVector(vector); (*it).second->exitVector(vector);
} }
} }
......
...@@ -250,7 +250,7 @@ namespace AMDiS { ...@@ -250,7 +250,7 @@ namespace AMDiS {
if (condition && condition->isDirichlet()) { if (condition && condition->isDirichlet()) {
MatrixRow *matrixRow = &(matrix[row]); MatrixRow *matrixRow = &(matrix[row]);
if(coupleMatrix) { if (coupleMatrix) {
matrixRow->resize(0); matrixRow->resize(0);
} else { } else {
matrixRow->resize(1); matrixRow->resize(1);
...@@ -447,9 +447,9 @@ namespace AMDiS { ...@@ -447,9 +447,9 @@ namespace AMDiS {
ElementMatrix *DOFMatrix::assemble(double factor, ElInfo *elInfo, ElementMatrix *DOFMatrix::assemble(double factor, ElInfo *elInfo,
const BoundaryType *bound, Operator *op) const BoundaryType *bound, Operator *op)
{ {
FUNCNAME("DOFMatrix::assemble"); FUNCNAME("DOFMatrix::assemble()");
if(!op && operators.size() == 0) { if (!op && operators.size() == 0) {
//WARNING("no operator\n"); //WARNING("no operator\n");
return NULL; return NULL;
} }
...@@ -459,12 +459,12 @@ namespace AMDiS { ...@@ -459,12 +459,12 @@ namespace AMDiS {
elementMatrix = elementMatrix =
operat->getAssembler()->initElementMatrix(elementMatrix, elInfo); operat->getAssembler()->initElementMatrix(elementMatrix, elInfo);
if(op) { if (op) {
op->getElementMatrix(elInfo, elementMatrix); op->getElementMatrix(elInfo, elementMatrix);
} else { } else {
::std::vector<Operator*>::iterator it; ::std::vector<Operator*>::iterator it;
::std::vector<double*>::iterator factorIt; ::std::vector<double*>::iterator factorIt;
for(it = operators.begin(), factorIt = operatorFactor.begin(); for (it = operators.begin(), factorIt = operatorFactor.begin();
it != operators.end(); it != operators.end();
++it, ++factorIt) ++it, ++factorIt)
{ {
......
...@@ -8,7 +8,9 @@ namespace AMDiS { ...@@ -8,7 +8,9 @@ namespace AMDiS {
DirichletBC::DirichletBC(BoundaryType type, DirichletBC::DirichletBC(BoundaryType type,
DOFVectorBase<double> *vec) DOFVectorBase<double> *vec)
: BoundaryCondition(type, vec->getFESpace()), f(NULL), dofVec(vec) : BoundaryCondition(type, vec->getFESpace()),
f(NULL),
dofVec(vec)
{} {}
void DirichletBC::fillBoundaryCondition(DOFMatrix* matrix, void DirichletBC::fillBoundaryCondition(DOFMatrix* matrix,
...@@ -18,10 +20,9 @@ namespace AMDiS { ...@@ -18,10 +20,9 @@ namespace AMDiS {
int nBasFcts) int nBasFcts)
{ {
FUNCNAME("DirichletBC::fillBoundaryCondition()"); FUNCNAME("DirichletBC::fillBoundaryCondition()");
TEST_EXIT(matrix->getRowFESpace() == rowFESpace) TEST_EXIT(matrix->getRowFESpace() == rowFESpace)
("invalid row fe space\n"); ("invalid row fe space\n");
TEST_EXIT(matrix->getColFESpace() == colFESpace)
("invalid col fe space\n");
} }
void DirichletBC::fillBoundaryCondition(DOFVectorBase<double>* vector, void DirichletBC::fillBoundaryCondition(DOFVectorBase<double>* vector,
...@@ -31,20 +32,20 @@ namespace AMDiS { ...@@ -31,20 +32,20 @@ namespace AMDiS {
int nBasFcts) int nBasFcts)
{ {
FUNCNAME("DirichletBC::fillBoundaryCondition()"); FUNCNAME("DirichletBC::fillBoundaryCondition()");
TEST_EXIT(vector->getFESpace() == rowFESpace) TEST_EXIT(vector->getFESpace() == rowFESpace)
("invalid row fe space\n"); ("invalid row fe space\n");
const BasisFunction *basFcts = rowFESpace->getBasisFcts(); const BasisFunction *basFcts = rowFESpace->getBasisFcts();
int i; for (int i = 0; i < nBasFcts; i++) {
for(i=0; i < nBasFcts; i++) { if (localBound[i] == boundaryType) {
if(localBound[i] == boundaryType) { if (f) {
if(f) {
DimVec<double> *coords = basFcts->getCoords(i); DimVec<double> *coords = basFcts->getCoords(i);
const WorldVector<double> *worldCoords = elInfo->coordToWorld(*coords, NULL); const WorldVector<double> *worldCoords = elInfo->coordToWorld(*coords, NULL);
double fAtCoords = f ? (*f)(*worldCoords) : 0; double fAtCoords = f ? (*f)(*worldCoords) : 0;
(*vector)[dofIndices[i]] = fAtCoords; (*vector)[dofIndices[i]] = fAtCoords;
} }
if(dofVec) { if (dofVec) {
(*vector)[dofIndices[i]] = (*dofVec)[dofIndices[i]]; (*vector)[dofIndices[i]] = (*dofVec)[dofIndices[i]];
} }
} }
......
...@@ -51,14 +51,17 @@ namespace AMDiS { ...@@ -51,14 +51,17 @@ namespace AMDiS {
*/ */
DirichletBC(BoundaryType type, DirichletBC(BoundaryType type,
AbstractFunction<double, WorldVector<double> > *fct, AbstractFunction<double, WorldVector<double> > *fct,
FiniteElemSpace *rowFESpace_) FiniteElemSpace *rowFESpace,
: BoundaryCondition(type, rowFESpace_, NULL), f(fct), dofVec(NULL) FiniteElemSpace *colFESpace = NULL)
: BoundaryCondition(type, rowFESpace, colFESpace),
f(fct),
dofVec(NULL)
{}; {};
/** \brief /** \brief
* Constructor. * Constructor.
*/ */
DirichletBC(BoundaryType type, DirichletBC(BoundaryType type,
DOFVectorBase<double> *vec); DOFVectorBase<double> *vec);
/** \brief /** \brief
...@@ -84,9 +87,13 @@ namespace AMDiS { ...@@ -84,9 +87,13 @@ namespace AMDiS {
*/ */
double boundResidual(ElInfo*, double boundResidual(ElInfo*,
DOFMatrix *, DOFMatrix *,
const DOFVectorBase<double>*) { return 0.0; }; const DOFVectorBase<double>*) {
return 0.0;
};
bool isDirichlet() { return true; }; bool isDirichlet() {
return true;
};
inline AbstractFunction<double, WorldVector<double> > *getF() { inline AbstractFunction<double, WorldVector<double> > *getF() {
return f; return f;
......
...@@ -134,17 +134,23 @@ namespace AMDiS { ...@@ -134,17 +134,23 @@ namespace AMDiS {
* Returns true if Element is a leaf element (\ref child[0] == NULL), returns * Returns true if Element is a leaf element (\ref child[0] == NULL), returns
* false otherwise. * false otherwise.
*/ */
inline const bool isLeaf() const { return (child[0]==NULL); }; inline const bool isLeaf() const {
return (child[0]==NULL);
};
/** \brief /** \brief
* Returns \ref dof[i][j] which is the j-th DOF of the i-th node of Element. * Returns \ref dof[i][j] which is the j-th DOF of the i-th node of Element.
*/ */
const DegreeOfFreedom getDOF(int i,int j) const { return dof[i][j];}; const DegreeOfFreedom getDOF(int i, int j) const {
return dof[i][j];
};
/** \brief /** \brief
* Returns \ref dof[i] which is a pointer to the DOFs of the i-th node. * Returns \ref dof[i] which is a pointer to the DOFs of the i-th node.
*/ */
const DegreeOfFreedom* getDOF(int i) const {return dof[i];}; const DegreeOfFreedom* getDOF(int i) const {
return dof[i];
};
/** \brief /** \brief
* Returns a pointer to the DOFs of this Element * Returns a pointer to the DOFs of this Element
...@@ -156,7 +162,9 @@ namespace AMDiS { ...@@ -156,7 +162,9 @@ namespace AMDiS {
/** \brief /** \brief
* Returns \ref mesh of Element * Returns \ref mesh of Element
*/ */
inline Mesh* getMesh() const { return mesh; }; inline Mesh* getMesh() const {
return mesh;
};
/** \brief /** \brief
* Returns \ref elementData's error estimation, if Element is a leaf element * Returns \ref elementData's error estimation, if Element is a leaf element
......
...@@ -311,6 +311,14 @@ namespace AMDiS { ...@@ -311,6 +311,14 @@ namespace AMDiS {
double old_res = -1.0; double old_res = -1.0;
// ::std::cout << "TEST!" << ::std::endl;
// matVec->matVec(NoTranspose, *x, *r);
// *r -= *b;
// r->print();
if (norm(b) < TOL) { if (norm(b) < TOL) {
INFO(this->info, 2)("b == 0, x = 0 is the solution of the linear system\n"); INFO(this->info, 2)("b == 0, x = 0 is the solution of the linear system\n");
setValue(*x, 0.0); setValue(*x, 0.0);
......
...@@ -216,7 +216,7 @@ namespace AMDiS { ...@@ -216,7 +216,7 @@ namespace AMDiS {
ElementMatrix *userMat, ElementMatrix *userMat,
double factor) double factor)
{ {
if(!assembler) { if (!assembler) {
initAssembler(NULL, NULL, NULL, NULL); initAssembler(NULL, NULL, NULL, NULL);
} }
...@@ -227,7 +227,7 @@ namespace AMDiS { ...@@ -227,7 +227,7 @@ namespace AMDiS {
ElementVector *userVec, ElementVector *userVec,
double factor) double factor)
{ {
if(!assembler) { if (!assembler) {
initAssembler(NULL, NULL, NULL, NULL); initAssembler(NULL, NULL, NULL, NULL);
} }
......
...@@ -809,7 +809,7 @@ namespace AMDiS { ...@@ -809,7 +809,7 @@ namespace AMDiS {
{ {
FUNCNAME("ProblemVec::addMatrixOperator()"); FUNCNAME("ProblemVec::addMatrixOperator()");
if(!(*systemMatrix_)[i][j]) { if (!(*systemMatrix_)[i][j]) {
TEST_EXIT(i != j)("should have been created already\n"); TEST_EXIT(i != j)("should have been created already\n");
(*systemMatrix_)[i][j] = NEW DOFMatrix(componentSpaces_[i], (*systemMatrix_)[i][j] = NEW DOFMatrix(componentSpaces_[i],
componentSpaces_[j], componentSpaces_[j],
...@@ -845,8 +845,10 @@ namespace AMDiS { ...@@ -845,8 +845,10 @@ namespace AMDiS {
(*systemMatrix_)[system][i]->getBoundaryManager()->addBoundaryCondition(dirichlet); (*systemMatrix_)[system][i]->getBoundaryManager()->addBoundaryCondition(dirichlet);
} }
} }
if (rhs_) if (rhs_)
rhs_->getDOFVector(system)->getBoundaryManager()->addBoundaryCondition(dirichlet); rhs_->getDOFVector(system)->getBoundaryManager()->addBoundaryCondition(dirichlet);
if (solution_) if (solution_)
solution_->getDOFVector(system)->getBoundaryManager()->addBoundaryCondition(dirichlet); solution_->getDOFVector(system)->getBoundaryManager()->addBoundaryCondition(dirichlet);
} }
...@@ -860,7 +862,7 @@ namespace AMDiS { ...@@ -860,7 +862,7 @@ namespace AMDiS {
new NeumannBC(type, n, new NeumannBC(type, n,
componentSpaces_[row], componentSpaces_[row],
componentSpaces_[col]); componentSpaces_[col]);
if(rhs_) if (rhs_)
rhs_->getDOFVector(row)->getBoundaryManager()->addBoundaryCondition(neumann); rhs_->getDOFVector(row)->getBoundaryManager()->addBoundaryCondition(neumann);
} }
...@@ -874,9 +876,10 @@ namespace AMDiS { ...@@ -874,9 +876,10 @@ namespace AMDiS {
new RobinBC(type, n, r, new RobinBC(type, n, r,
componentSpaces_[row], componentSpaces_[row],
componentSpaces_[col]); componentSpaces_[col]);
if(rhs_) if (rhs_)
rhs_->getDOFVector(row)->getBoundaryManager()->addBoundaryCondition(robin); rhs_->getDOFVector(row)->getBoundaryManager()->addBoundaryCondition(robin);
if(systemMatrix_ && (*systemMatrix_)[row][col]) {
if (systemMatrix_ && (*systemMatrix_)[row][col]) {
(*systemMatrix_)[row][col]->getBoundaryManager()->addBoundaryCondition(robin); (*systemMatrix_)[row][col]->getBoundaryManager()->addBoundaryCondition(robin);
} }
} }
......
...@@ -446,7 +446,9 @@ namespace AMDiS { ...@@ -446,7 +446,9 @@ namespace AMDiS {
/** \brief /** \brief
* Sets \ref solver_. * Sets \ref solver_.
*/ */
inline void setSolver(OEMSolver<SystemVector>* sol) { solver_ = sol; }; inline void setSolver(OEMSolver<SystemVector>* sol) {
solver_ = sol;
};
/** \brief /** \brief
* Sets \ref leftPrecon_. * Sets \ref leftPrecon_.
......