Commit 5cede7ac authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

Dirichlet boundary conditions fixed.

parent 3ae1479f
...@@ -46,21 +46,25 @@ namespace AMDiS { ...@@ -46,21 +46,25 @@ namespace AMDiS {
rowFESpace(rowFESpace_), rowFESpace(rowFESpace_),
colFESpace(colFESpace_) colFESpace(colFESpace_)
{ {
if(!colFESpace) colFESpace = rowFESpace; if (!colFESpace)
colFESpace = rowFESpace;
} }
/// Returns \ref boundaryType. /// Returns \ref boundaryType.
inline BoundaryType getBoundaryType() { inline BoundaryType getBoundaryType()
{
return boundaryType; return boundaryType;
} }
/// Returns \ref rowFESpace. /// Returns \ref rowFESpace.
inline const FiniteElemSpace *getRowFESpace() { inline const FiniteElemSpace *getRowFESpace()
{
return rowFESpace; return rowFESpace;
} }
/// Returns \ref rowFESpace. /// Returns \ref rowFESpace.
inline const FiniteElemSpace *getColFESpace() { inline const FiniteElemSpace *getColFESpace()
{
return colFESpace; return colFESpace;
} }
...@@ -106,13 +110,27 @@ namespace AMDiS { ...@@ -106,13 +110,27 @@ namespace AMDiS {
} }
/** \brief /** \brief
* Returns whether the condition must be treated as dirichlet condition * Returns whether the condition must be treated as Dirichlet condition
* while assemblage. * while assemblage.
*/ */
virtual bool isDirichlet() { virtual bool isDirichlet()
{
return false; 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: protected:
/** \brief /** \brief
* Speciefies for which parts of the boundary the condition holds. * Speciefies for which parts of the boundary the condition holds.
......
...@@ -48,10 +48,10 @@ namespace AMDiS { ...@@ -48,10 +48,10 @@ namespace AMDiS {
~BoundaryManager(); ~BoundaryManager();
/// Adds a local boundary condition to the list of managed conditions. /// Adds a local boundary condition to the list of managed conditions.
void addBoundaryCondition(BoundaryCondition *localBC) { void addBoundaryCondition(BoundaryCondition *localBC)
{
BoundaryType type = localBC->getBoundaryType(); BoundaryType type = localBC->getBoundaryType();
TEST_EXIT(localBCs[type] == NULL) TEST_EXIT(localBCs[type] == NULL)("there is already a condition for this type\n");
("there is already a condition for this type\n");
localBCs[type] = localBC; localBCs[type] = localBC;
} }
...@@ -83,15 +83,18 @@ namespace AMDiS { ...@@ -83,15 +83,18 @@ namespace AMDiS {
DOFMatrix *matrix, DOFMatrix *matrix,
const DOFVectorBase<double> *dv); const DOFVectorBase<double> *dv);
inline BoundaryCondition *getBoundaryCondition(BoundaryType type) { inline BoundaryCondition *getBoundaryCondition(BoundaryType type)
{
return localBCs[type]; return localBCs[type];
} }
const std::map<BoundaryType, BoundaryCondition*>& getBoundaryConditionMap() { const std::map<BoundaryType, BoundaryCondition*>& getBoundaryConditionMap()
{
return localBCs; return localBCs;
} }
void setBoundaryConditionMap(const std::map<BoundaryType, BoundaryCondition*>& bcs) { void setBoundaryConditionMap(const std::map<BoundaryType, BoundaryCondition*>& bcs)
{
localBCs = bcs; localBCs = bcs;
} }
......
...@@ -100,9 +100,11 @@ namespace AMDiS { ...@@ -100,9 +100,11 @@ namespace AMDiS {
typedef traits::range_generator<nz, cursor_type>::type icursor_type; typedef traits::range_generator<nz, cursor_type>::type icursor_type;
std::cout.precision(10); std::cout.precision(10);
for (cursor_type cursor = begin<major>(matrix), cend = end<major>(matrix); cursor != cend; ++cursor) 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) for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); icursor != icend; ++icursor)
std::cout << row(*icursor) << " " << col(*icursor) << " " << value(*icursor) << "\n"; std::cout << "(" << row(*icursor) << "," << col(*icursor) << "," << value(*icursor) << ") ";
std::cout << "\n";
}
} }
bool DOFMatrix::symmetric() bool DOFMatrix::symmetric()
...@@ -202,21 +204,19 @@ namespace AMDiS { ...@@ -202,21 +204,19 @@ namespace AMDiS {
// === Add element matrix to the global matrix using the indices mapping. === // === 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 for (int i = 0; i < nRow; i++) { // for all rows of element matrix
row = rowIndices[i]; DegreeOfFreedom row = rowIndices[i];
BoundaryCondition *condition = BoundaryCondition *condition =
bound ? boundaryManager->getBoundaryCondition(bound[i]) : NULL; bound ? boundaryManager->getBoundaryCondition(bound[i]) : NULL;
if (condition && condition->isDirichlet()) { if (condition && condition->isDirichlet()) {
if (!coupleMatrix) if (condition->applyBoundaryCondition())
applyDBCs.insert(static_cast<int>(row)); applyDBCs.insert(static_cast<int>(row));
} else } else
for (int j = 0; j < nCol; j++) { // for all columns for (int j = 0; j < nCol; j++) { // for all columns
col = colIndices[j]; DegreeOfFreedom col = colIndices[j];
entry = elMat[i][j]; double entry = elMat[i][j];
if (add) if (add)
ins[row][col]+= sign * entry; ins[row][col]+= sign * entry;
...@@ -331,19 +331,17 @@ namespace AMDiS { ...@@ -331,19 +331,17 @@ namespace AMDiS {
// call the operatos cleanup procedures // call the operatos cleanup procedures
for (std::vector<Operator*>::iterator it = operators.begin(); for (std::vector<Operator*>::iterator it = operators.begin();
it != operators.end(); it != operators.end();
++it) { ++it)
(*it)->finishAssembling(); (*it)->finishAssembling();
}
} }
// Should work as before // Should work as before
Flag DOFMatrix::getAssembleFlag() Flag DOFMatrix::getAssembleFlag()
{ {
Flag fillFlag(0); Flag fillFlag(0);
std::vector<Operator*>::iterator op; for (std::vector<Operator*>::iterator op = operators.begin();
for(op = operators.begin(); op != operators.end(); ++op) { op != operators.end(); ++op)
fillFlag |= (*op)->getFillFlag(); fillFlag |= (*op)->getFillFlag();
}
return fillFlag; return fillFlag;
} }
...@@ -381,71 +379,6 @@ namespace AMDiS { ...@@ -381,71 +379,6 @@ namespace AMDiS {
rows->clear(); 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() int DOFMatrix::memsize()
{ {
return (num_rows(matrix) + matrix.nnz()) * sizeof(base_matrix_type::size_type) return (num_rows(matrix) + matrix.nnz()) * sizeof(base_matrix_type::size_type)
......
...@@ -347,8 +347,6 @@ namespace AMDiS { ...@@ -347,8 +347,6 @@ namespace AMDiS {
boundaryManager = bm; boundaryManager = bm;
} }
void createPictureFile(const char* filename, int dim);
private: private:
template <typename T> template <typename T>
void s_write(::std::ostream &out, const T& value) void s_write(::std::ostream &out, const T& value)
......
...@@ -9,10 +9,12 @@ namespace AMDiS { ...@@ -9,10 +9,12 @@ namespace AMDiS {
DirichletBC::DirichletBC(BoundaryType type, DirichletBC::DirichletBC(BoundaryType type,
AbstractFunction<double, WorldVector<double> > *fct, AbstractFunction<double, WorldVector<double> > *fct,
FiniteElemSpace *rowFESpace, FiniteElemSpace *rowFESpace,
FiniteElemSpace *colFESpace) FiniteElemSpace *colFESpace,
bool apply)
: BoundaryCondition(type, rowFESpace, colFESpace), : BoundaryCondition(type, rowFESpace, colFESpace),
f(fct), f(fct),
dofVec(NULL) dofVec(NULL),
applyBC(apply)
{ {
worldCoords.resize(omp_get_overall_max_threads()); worldCoords.resize(omp_get_overall_max_threads());
} }
...@@ -21,7 +23,8 @@ namespace AMDiS { ...@@ -21,7 +23,8 @@ namespace AMDiS {
DOFVectorBase<double> *vec) DOFVectorBase<double> *vec)
: BoundaryCondition(type, vec->getFESpace()), : BoundaryCondition(type, vec->getFESpace()),
f(NULL), f(NULL),
dofVec(vec) dofVec(vec),
applyBC(true)
{ {
worldCoords.resize(omp_get_overall_max_threads()); worldCoords.resize(omp_get_overall_max_threads());
} }
...@@ -34,8 +37,7 @@ namespace AMDiS { ...@@ -34,8 +37,7 @@ namespace AMDiS {
{ {
FUNCNAME("DirichletBC::fillBoundaryCondition()"); FUNCNAME("DirichletBC::fillBoundaryCondition()");
TEST_EXIT_DBG(matrix->getRowFESpace() == rowFESpace) TEST_EXIT_DBG(matrix->getRowFESpace() == rowFESpace)("invalid row fe space\n");
("invalid row fe space\n");
} }
void DirichletBC::fillBoundaryCondition(DOFVectorBase<double>* vector, void DirichletBC::fillBoundaryCondition(DOFVectorBase<double>* vector,
...@@ -46,8 +48,7 @@ namespace AMDiS { ...@@ -46,8 +48,7 @@ namespace AMDiS {
{ {
FUNCNAME("DirichletBC::fillBoundaryCondition()"); FUNCNAME("DirichletBC::fillBoundaryCondition()");
TEST_EXIT_DBG(vector->getFESpace() == rowFESpace) TEST_EXIT_DBG(vector->getFESpace() == rowFESpace)("invalid row fe space\n");
("invalid row fe space\n");
const BasisFunction *basFcts = rowFESpace->getBasisFcts(); const BasisFunction *basFcts = rowFESpace->getBasisFcts();
int myRank = omp_get_thread_num(); int myRank = omp_get_thread_num();
...@@ -59,9 +60,8 @@ namespace AMDiS { ...@@ -59,9 +60,8 @@ namespace AMDiS {
double fAtCoords = (*f)(worldCoords[myRank]); double fAtCoords = (*f)(worldCoords[myRank]);
(*vector)[dofIndices[i]] = fAtCoords; (*vector)[dofIndices[i]] = fAtCoords;
} }
if (dofVec) { if (dofVec)
(*vector)[dofIndices[i]] = (*dofVec)[dofIndices[i]]; (*vector)[dofIndices[i]] = (*dofVec)[dofIndices[i]];
}
} }
} }
} }
......
...@@ -44,7 +44,8 @@ namespace AMDiS { ...@@ -44,7 +44,8 @@ namespace AMDiS {
DirichletBC(BoundaryType type, DirichletBC(BoundaryType type,
AbstractFunction<double, WorldVector<double> > *fct, AbstractFunction<double, WorldVector<double> > *fct,
FiniteElemSpace *rowFESpace, FiniteElemSpace *rowFESpace,
FiniteElemSpace *colFESpace = NULL); FiniteElemSpace *colFESpace = NULL,
bool apply = true);
/// Constructor. /// Constructor.
DirichletBC(BoundaryType type, DirichletBC(BoundaryType type,
...@@ -72,15 +73,25 @@ namespace AMDiS { ...@@ -72,15 +73,25 @@ namespace AMDiS {
return 0.0; return 0.0;
} }
bool isDirichlet() { /// Because this is a Dirichlet boundary condition, always return true.
bool isDirichlet()
{
return true; return true;
} }
inline AbstractFunction<double, WorldVector<double> > *getF() { /// Returns \ref applyBC.
bool applyBoundaryCondition()
{
return applyBC;
}
inline AbstractFunction<double, WorldVector<double> > *getF()
{
return f; return f;
} }
inline DOFVectorBase<double> *getDOFVector() { inline DOFVectorBase<double> *getDOFVector()
{
return dofVec; return dofVec;
} }
...@@ -88,10 +99,17 @@ namespace AMDiS { ...@@ -88,10 +99,17 @@ namespace AMDiS {
/// Function which is evaluated at world coords of Dirichlet dofs. /// Function which is evaluated at world coords of Dirichlet dofs.
AbstractFunction<double, WorldVector<double> > *f; AbstractFunction<double, WorldVector<double> > *f;
///
std::vector<WorldVector<double> > worldCoords; std::vector<WorldVector<double> > worldCoords;
/// DOFVector containing the boundary values /// DOFVector containing the boundary values
DOFVectorBase<double> *dofVec; DOFVectorBase<double> *dofVec;
/** \brief
* Defines, if the boundary condition must be applied to the matrix. See
* comment of \ref BoundaryCondition::applyBoundaryCondition.
*/
bool applyBC;
}; };
} }
......
...@@ -916,23 +916,27 @@ namespace AMDiS { ...@@ -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) AbstractFunction<double, WorldVector<double> >* b)
{ {
FUNCNAME("ProblemVec::addDirichletBC()"); FUNCNAME("ProblemVec::addDirichletBC()");
DirichletBC *dirichlet = new DirichletBC(type, DirichletBC *dirichletApply =
b, new DirichletBC(type, b, componentSpaces[row], componentSpaces[col], true);
componentSpaces[system]); DirichletBC *dirichletNotApply =
for (int i = 0; i < nComponents; i++) { new DirichletBC(type, b, componentSpaces[row], componentSpaces[col], false);
if (systemMatrix && (*systemMatrix)[system][i]) {
(*systemMatrix)[system][i]->getBoundaryManager()->addBoundaryCondition(dirichlet); 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) if (rhs)
rhs->getDOFVector(system)->getBoundaryManager()->addBoundaryCondition(dirichlet); rhs->getDOFVector(row)->getBoundaryManager()->addBoundaryCondition(dirichletApply);
if (solution) if (solution)
solution->getDOFVector(system)->getBoundaryManager()->addBoundaryCondition(dirichlet); solution->getDOFVector(row)->getBoundaryManager()->addBoundaryCondition(dirichletApply);
} }
void ProblemVec::addNeumannBC(BoundaryType type, int row, int col, void ProblemVec::addNeumannBC(BoundaryType type, int row, int col,
...@@ -941,9 +945,8 @@ namespace AMDiS { ...@@ -941,9 +945,8 @@ namespace AMDiS {
FUNCNAME("ProblemVec::addNeumannBC()"); FUNCNAME("ProblemVec::addNeumannBC()");
NeumannBC *neumann = NeumannBC *neumann =
new NeumannBC(type, n, new NeumannBC(type, n, componentSpaces[row], componentSpaces[col]);
componentSpaces[row],
componentSpaces[col]);
if (rhs) if (rhs)
rhs->getDOFVector(row)->getBoundaryManager()->addBoundaryCondition(neumann); rhs->getDOFVector(row)->getBoundaryManager()->addBoundaryCondition(neumann);
} }
...@@ -955,15 +958,12 @@ namespace AMDiS { ...@@ -955,15 +958,12 @@ namespace AMDiS {
FUNCNAME("ProblemVec::addRobinBC()"); FUNCNAME("ProblemVec::addRobinBC()");
RobinBC *robin = RobinBC *robin =
new RobinBC(type, n, r, new RobinBC(type, n, r, componentSpaces[row], componentSpaces[col]);
componentSpaces[row],
componentSpaces[col]);
if (rhs)
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);
} if (rhs)
rhs->getDOFVector(row)->getBoundaryManager()->addBoundaryCondition(robin);
} }
void ProblemVec::addPeriodicBC(BoundaryType type, int row, int col) void ProblemVec::addPeriodicBC(BoundaryType type, int row, int col)
......
...@@ -208,7 +208,7 @@ namespace AMDiS { ...@@ -208,7 +208,7 @@ namespace AMDiS {
bool fake = false); bool fake = false);
/// Adds dirichlet boundary conditions. /// 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); AbstractFunction<double, WorldVector<double> > *b);
/// Adds neumann boundary conditions. /// Adds neumann boundary conditions.
......
Supports Markdown
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