From 5cede7ac9053a4b6ea2f5511d36636a010e8fa88 Mon Sep 17 00:00:00 2001
From: Thomas Witkowski <thomas.witkowski@gmx.de>
Date: Mon, 25 May 2009 09:15:53 +0000
Subject: [PATCH] Dirichlet boundary conditions fixed.
---
AMDiS/src/BoundaryCondition.h | 30 +++++++++---
AMDiS/src/BoundaryManager.h | 15 +++---
AMDiS/src/DOFMatrix.cc | 91 +++++------------------------------
AMDiS/src/DOFMatrix.h | 2 -
AMDiS/src/DirichletBC.cc | 18 +++----
AMDiS/src/DirichletBC.h | 26 ++++++++--
AMDiS/src/ProblemVec.cc | 42 ++++++++--------
AMDiS/src/ProblemVec.h | 2 +-
8 files changed, 98 insertions(+), 128 deletions(-)
diff --git a/AMDiS/src/BoundaryCondition.h b/AMDiS/src/BoundaryCondition.h
index 33d62b1a..94a8371a 100644
--- a/AMDiS/src/BoundaryCondition.h
+++ b/AMDiS/src/BoundaryCondition.h
@@ -46,21 +46,25 @@ namespace AMDiS {
rowFESpace(rowFESpace_),
colFESpace(colFESpace_)
{
- if(!colFESpace) colFESpace = rowFESpace;
+ if (!colFESpace)
+ colFESpace = rowFESpace;
}
/// Returns \ref boundaryType.
- inline BoundaryType getBoundaryType() {
+ inline BoundaryType getBoundaryType()
+ {
return boundaryType;
}
/// Returns \ref rowFESpace.
- inline const FiniteElemSpace *getRowFESpace() {
+ inline const FiniteElemSpace *getRowFESpace()
+ {
return rowFESpace;
}
/// Returns \ref rowFESpace.
- inline const FiniteElemSpace *getColFESpace() {
+ inline const FiniteElemSpace *getColFESpace()
+ {
return colFESpace;
}
@@ -106,13 +110,27 @@ namespace AMDiS {
}
/** \brief
- * Returns whether the condition must be treated as dirichlet condition
+ * Returns whether the condition must be treated as Dirichlet condition
* while assemblage.
*/
- virtual bool isDirichlet() {
+ virtual bool isDirichlet()
+ {
return false;
}
+ /** \brief
+ * In some situations it may be required to set Dirichlet boundary conditions,
+ * but not to apply them to the matrix. This is for example the case, if the
+ * boundary condition is set to a couple matrix. Then, the boundary conditions
+ * must be applied to the couple matrix, but they are set to all matrices in this
+ * row (to ensure that there are no other element entries in the Dirichlet boundary
+ * condition rows).
+ */
+ virtual bool applyBoundaryCondition()
+ {
+ return true;
+ }
+
protected:
/** \brief
* Speciefies for which parts of the boundary the condition holds.
diff --git a/AMDiS/src/BoundaryManager.h b/AMDiS/src/BoundaryManager.h
index 6aa72c70..13bfd66d 100644
--- a/AMDiS/src/BoundaryManager.h
+++ b/AMDiS/src/BoundaryManager.h
@@ -48,10 +48,10 @@ namespace AMDiS {
~BoundaryManager();
/// Adds a local boundary condition to the list of managed conditions.
- void addBoundaryCondition(BoundaryCondition *localBC) {
+ void addBoundaryCondition(BoundaryCondition *localBC)
+ {
BoundaryType type = localBC->getBoundaryType();
- TEST_EXIT(localBCs[type] == NULL)
- ("there is already a condition for this type\n");
+ TEST_EXIT(localBCs[type] == NULL)("there is already a condition for this type\n");
localBCs[type] = localBC;
}
@@ -83,15 +83,18 @@ namespace AMDiS {
DOFMatrix *matrix,
const DOFVectorBase<double> *dv);
- inline BoundaryCondition *getBoundaryCondition(BoundaryType type) {
+ inline BoundaryCondition *getBoundaryCondition(BoundaryType type)
+ {
return localBCs[type];
}
- const std::map<BoundaryType, BoundaryCondition*>& getBoundaryConditionMap() {
+ const std::map<BoundaryType, BoundaryCondition*>& getBoundaryConditionMap()
+ {
return localBCs;
}
- void setBoundaryConditionMap(const std::map<BoundaryType, BoundaryCondition*>& bcs) {
+ void setBoundaryConditionMap(const std::map<BoundaryType, BoundaryCondition*>& bcs)
+ {
localBCs = bcs;
}
diff --git a/AMDiS/src/DOFMatrix.cc b/AMDiS/src/DOFMatrix.cc
index c1f3116b..04716a56 100644
--- a/AMDiS/src/DOFMatrix.cc
+++ b/AMDiS/src/DOFMatrix.cc
@@ -100,9 +100,11 @@ namespace AMDiS {
typedef traits::range_generator<nz, cursor_type>::type icursor_type;
std::cout.precision(10);
- for (cursor_type cursor = begin<major>(matrix), cend = end<major>(matrix); cursor != cend; ++cursor)
- for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); icursor != icend; ++icursor)
- std::cout << row(*icursor) << " " << col(*icursor) << " " << value(*icursor) << "\n";
+ for (cursor_type cursor = begin<major>(matrix), cend = end<major>(matrix); cursor != cend; ++cursor) {
+ for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); icursor != icend; ++icursor)
+ std::cout << "(" << row(*icursor) << "," << col(*icursor) << "," << value(*icursor) << ") ";
+ std::cout << "\n";
+ }
}
bool DOFMatrix::symmetric()
@@ -202,21 +204,19 @@ namespace AMDiS {
// === Add element matrix to the global matrix using the indices mapping. ===
- DegreeOfFreedom row, col;
- double entry;
for (int i = 0; i < nRow; i++) { // for all rows of element matrix
- row = rowIndices[i];
+ DegreeOfFreedom row = rowIndices[i];
BoundaryCondition *condition =
bound ? boundaryManager->getBoundaryCondition(bound[i]) : NULL;
if (condition && condition->isDirichlet()) {
- if (!coupleMatrix)
+ if (condition->applyBoundaryCondition())
applyDBCs.insert(static_cast<int>(row));
} else
for (int j = 0; j < nCol; j++) { // for all columns
- col = colIndices[j];
- entry = elMat[i][j];
+ DegreeOfFreedom col = colIndices[j];
+ double entry = elMat[i][j];
if (add)
ins[row][col]+= sign * entry;
@@ -331,19 +331,17 @@ namespace AMDiS {
// call the operatos cleanup procedures
for (std::vector<Operator*>::iterator it = operators.begin();
it != operators.end();
- ++it) {
+ ++it)
(*it)->finishAssembling();
- }
}
// Should work as before
Flag DOFMatrix::getAssembleFlag()
{
Flag fillFlag(0);
- std::vector<Operator*>::iterator op;
- for(op = operators.begin(); op != operators.end(); ++op) {
+ for (std::vector<Operator*>::iterator op = operators.begin();
+ op != operators.end(); ++op)
fillFlag |= (*op)->getFillFlag();
- }
return fillFlag;
}
@@ -381,71 +379,6 @@ namespace AMDiS {
rows->clear();
}
- void DOFMatrix::createPictureFile(const char* filename, int dim)
- {
- TEST_EXIT(0)("Not yet re-implemented.");
-
-#if 0
- png_structp png_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING,
- NULL, NULL, NULL);
-
- if (!png_ptr)
- return;
-
- png_bytep rowPointers[dim];
- for (int i = 0; i < dim; i++) {
- rowPointers[i] = (png_byte*)png_malloc(png_ptr, dim);
-
- for (int j = 0; j < dim; j++)
- rowPointers[i][j] = 255;
- }
-
- double scalFactor = static_cast<double>(dim) / static_cast<double>(matrix.size());
-
- for (int i = 0; i < static_cast<int>(matrix.size()); i++) {
- int pi = static_cast<int>(static_cast<double>(i) * scalFactor);
-
- TEST_EXIT_DBG((pi >= 0) && (pi < dim))("PI");
-
- for (int j = 0; j < static_cast<int>(matrix[i].size()); j++) {
-
- int pj = static_cast<int>(static_cast<double>(matrix[i][j].col) * scalFactor);
-
- TEST_EXIT_DBG((pj >= 0) && (pj < dim))("PJ");
-
- rowPointers[pi][pj] = 0;
- }
- }
-
- FILE *fp = fopen(filename, "wb");
- TEST_EXIT(fp)("Cannot open file for writing matrix picture file!\n");
-
- png_infop info_ptr = png_create_info_struct(png_ptr);
-
- if (!info_ptr) {
- png_destroy_write_struct(&png_ptr, (png_infopp)NULL);
- return;
- }
-
- png_init_io(png_ptr, fp);
-
- png_set_IHDR(png_ptr, info_ptr, dim, dim, 8,
- PNG_COLOR_TYPE_GRAY,
- PNG_INTERLACE_NONE,
- PNG_COMPRESSION_TYPE_DEFAULT,
- PNG_FILTER_TYPE_DEFAULT);
-
- png_set_rows(png_ptr, info_ptr, rowPointers);
-
- png_write_png(png_ptr, info_ptr, PNG_TRANSFORM_IDENTITY, NULL);
-
- png_destroy_write_struct(&png_ptr, &info_ptr);
-
- fclose(fp);
-#endif
- }
-
-
int DOFMatrix::memsize()
{
return (num_rows(matrix) + matrix.nnz()) * sizeof(base_matrix_type::size_type)
diff --git a/AMDiS/src/DOFMatrix.h b/AMDiS/src/DOFMatrix.h
index a72f9ebd..5771b2c1 100644
--- a/AMDiS/src/DOFMatrix.h
+++ b/AMDiS/src/DOFMatrix.h
@@ -347,8 +347,6 @@ namespace AMDiS {
boundaryManager = bm;
}
- void createPictureFile(const char* filename, int dim);
-
private:
template <typename T>
void s_write(::std::ostream &out, const T& value)
diff --git a/AMDiS/src/DirichletBC.cc b/AMDiS/src/DirichletBC.cc
index 15652327..1c623dfd 100644
--- a/AMDiS/src/DirichletBC.cc
+++ b/AMDiS/src/DirichletBC.cc
@@ -9,10 +9,12 @@ namespace AMDiS {
DirichletBC::DirichletBC(BoundaryType type,
AbstractFunction<double, WorldVector<double> > *fct,
FiniteElemSpace *rowFESpace,
- FiniteElemSpace *colFESpace)
+ FiniteElemSpace *colFESpace,
+ bool apply)
: BoundaryCondition(type, rowFESpace, colFESpace),
f(fct),
- dofVec(NULL)
+ dofVec(NULL),
+ applyBC(apply)
{
worldCoords.resize(omp_get_overall_max_threads());
}
@@ -21,7 +23,8 @@ namespace AMDiS {
DOFVectorBase<double> *vec)
: BoundaryCondition(type, vec->getFESpace()),
f(NULL),
- dofVec(vec)
+ dofVec(vec),
+ applyBC(true)
{
worldCoords.resize(omp_get_overall_max_threads());
}
@@ -34,8 +37,7 @@ namespace AMDiS {
{
FUNCNAME("DirichletBC::fillBoundaryCondition()");
- TEST_EXIT_DBG(matrix->getRowFESpace() == rowFESpace)
- ("invalid row fe space\n");
+ TEST_EXIT_DBG(matrix->getRowFESpace() == rowFESpace)("invalid row fe space\n");
}
void DirichletBC::fillBoundaryCondition(DOFVectorBase<double>* vector,
@@ -46,8 +48,7 @@ namespace AMDiS {
{
FUNCNAME("DirichletBC::fillBoundaryCondition()");
- TEST_EXIT_DBG(vector->getFESpace() == rowFESpace)
- ("invalid row fe space\n");
+ TEST_EXIT_DBG(vector->getFESpace() == rowFESpace)("invalid row fe space\n");
const BasisFunction *basFcts = rowFESpace->getBasisFcts();
int myRank = omp_get_thread_num();
@@ -59,9 +60,8 @@ namespace AMDiS {
double fAtCoords = (*f)(worldCoords[myRank]);
(*vector)[dofIndices[i]] = fAtCoords;
}
- if (dofVec) {
+ if (dofVec)
(*vector)[dofIndices[i]] = (*dofVec)[dofIndices[i]];
- }
}
}
}
diff --git a/AMDiS/src/DirichletBC.h b/AMDiS/src/DirichletBC.h
index 082a8fc1..fc59f594 100644
--- a/AMDiS/src/DirichletBC.h
+++ b/AMDiS/src/DirichletBC.h
@@ -44,7 +44,8 @@ namespace AMDiS {
DirichletBC(BoundaryType type,
AbstractFunction<double, WorldVector<double> > *fct,
FiniteElemSpace *rowFESpace,
- FiniteElemSpace *colFESpace = NULL);
+ FiniteElemSpace *colFESpace = NULL,
+ bool apply = true);
/// Constructor.
DirichletBC(BoundaryType type,
@@ -72,15 +73,25 @@ namespace AMDiS {
return 0.0;
}
- bool isDirichlet() {
+ /// Because this is a Dirichlet boundary condition, always return true.
+ bool isDirichlet()
+ {
return true;
}
- inline AbstractFunction<double, WorldVector<double> > *getF() {
+ /// Returns \ref applyBC.
+ bool applyBoundaryCondition()
+ {
+ return applyBC;
+ }
+
+ inline AbstractFunction<double, WorldVector<double> > *getF()
+ {
return f;
}
- inline DOFVectorBase<double> *getDOFVector() {
+ inline DOFVectorBase<double> *getDOFVector()
+ {
return dofVec;
}
@@ -88,10 +99,17 @@ namespace AMDiS {
/// Function which is evaluated at world coords of Dirichlet dofs.
AbstractFunction<double, WorldVector<double> > *f;
+ ///
std::vector<WorldVector<double> > worldCoords;
/// DOFVector containing the boundary values
DOFVectorBase<double> *dofVec;
+
+ /** \brief
+ * Defines, if the boundary condition must be applied to the matrix. See
+ * comment of \ref BoundaryCondition::applyBoundaryCondition.
+ */
+ bool applyBC;
};
}
diff --git a/AMDiS/src/ProblemVec.cc b/AMDiS/src/ProblemVec.cc
index 1159d6cd..6b7bf79b 100644
--- a/AMDiS/src/ProblemVec.cc
+++ b/AMDiS/src/ProblemVec.cc
@@ -916,23 +916,27 @@ namespace AMDiS {
}
}
- void ProblemVec::addDirichletBC(BoundaryType type, int system,
+ void ProblemVec::addDirichletBC(BoundaryType type, int row, int col,
AbstractFunction<double, WorldVector<double> >* b)
{
FUNCNAME("ProblemVec::addDirichletBC()");
- DirichletBC *dirichlet = new DirichletBC(type,
- b,
- componentSpaces[system]);
- for (int i = 0; i < nComponents; i++) {
- if (systemMatrix && (*systemMatrix)[system][i]) {
- (*systemMatrix)[system][i]->getBoundaryManager()->addBoundaryCondition(dirichlet);
- }
- }
+ DirichletBC *dirichletApply =
+ new DirichletBC(type, b, componentSpaces[row], componentSpaces[col], true);
+ DirichletBC *dirichletNotApply =
+ new DirichletBC(type, b, componentSpaces[row], componentSpaces[col], false);
+
+ for (int i = 0; i < nComponents; i++)
+ if (systemMatrix && (*systemMatrix)[row][i])
+ if (i == col)
+ (*systemMatrix)[row][i]->getBoundaryManager()->addBoundaryCondition(dirichletApply);
+ else
+ (*systemMatrix)[row][i]->getBoundaryManager()->addBoundaryCondition(dirichletNotApply);
+
if (rhs)
- rhs->getDOFVector(system)->getBoundaryManager()->addBoundaryCondition(dirichlet);
+ rhs->getDOFVector(row)->getBoundaryManager()->addBoundaryCondition(dirichletApply);
if (solution)
- solution->getDOFVector(system)->getBoundaryManager()->addBoundaryCondition(dirichlet);
+ solution->getDOFVector(row)->getBoundaryManager()->addBoundaryCondition(dirichletApply);
}
void ProblemVec::addNeumannBC(BoundaryType type, int row, int col,
@@ -941,9 +945,8 @@ namespace AMDiS {
FUNCNAME("ProblemVec::addNeumannBC()");
NeumannBC *neumann =
- new NeumannBC(type, n,
- componentSpaces[row],
- componentSpaces[col]);
+ new NeumannBC(type, n, componentSpaces[row], componentSpaces[col]);
+
if (rhs)
rhs->getDOFVector(row)->getBoundaryManager()->addBoundaryCondition(neumann);
}
@@ -955,15 +958,12 @@ namespace AMDiS {
FUNCNAME("ProblemVec::addRobinBC()");
RobinBC *robin =
- new RobinBC(type, n, r,
- componentSpaces[row],
- componentSpaces[col]);
- if (rhs)
- rhs->getDOFVector(row)->getBoundaryManager()->addBoundaryCondition(robin);
+ new RobinBC(type, n, r, componentSpaces[row], componentSpaces[col]);
- if (systemMatrix && (*systemMatrix)[row][col]) {
+ if (systemMatrix && (*systemMatrix)[row][col])
(*systemMatrix)[row][col]->getBoundaryManager()->addBoundaryCondition(robin);
- }
+ if (rhs)
+ rhs->getDOFVector(row)->getBoundaryManager()->addBoundaryCondition(robin);
}
void ProblemVec::addPeriodicBC(BoundaryType type, int row, int col)
diff --git a/AMDiS/src/ProblemVec.h b/AMDiS/src/ProblemVec.h
index da66a0dd..3300fce3 100644
--- a/AMDiS/src/ProblemVec.h
+++ b/AMDiS/src/ProblemVec.h
@@ -208,7 +208,7 @@ namespace AMDiS {
bool fake = false);
/// Adds dirichlet boundary conditions.
- virtual void addDirichletBC(BoundaryType type, int system,
+ virtual void addDirichletBC(BoundaryType type, int row, int col,
AbstractFunction<double, WorldVector<double> > *b);
/// Adds neumann boundary conditions.
--
GitLab