Commit 7856d906 authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

Periodic boundary conditions with mixed finite elements in parallel should work now.

parent 6bfde321
...@@ -57,7 +57,7 @@ namespace AMDiS { ...@@ -57,7 +57,7 @@ namespace AMDiS {
FUNCNAME("DOFMatrix::DOFMatrix()"); FUNCNAME("DOFMatrix::DOFMatrix()");
TEST_EXIT(rowFeSpace)("No fe space for row!\n"); TEST_EXIT(rowFeSpace)("No fe space for row!\n");
if (!colFeSpace) if (!colFeSpace)
colFeSpace = rowFeSpace; colFeSpace = rowFeSpace;
...@@ -211,16 +211,18 @@ namespace AMDiS { ...@@ -211,16 +211,18 @@ namespace AMDiS {
using namespace mtl; using namespace mtl;
#if 0 #if 0
std::cout << "----- PRINT MAT " << rowElInfo->getElement()->getIndex() << "--------" << std::endl; if (MPI::COMM_WORLD.Get_rank() == 0) {
std::cout << elMat << std::endl; std::cout << "----- PRINT MAT " << rowElInfo->getElement()->getIndex() << "--------" << std::endl;
std::cout << "rows: "; std::cout << elMat << std::endl;
for (int i = 0; i < rowIndices.size(); i++) std::cout << "rows: ";
std::cout << rowIndices[i] << " "; for (int i = 0; i < rowIndices.size(); i++)
std::cout << std::endl; std::cout << rowIndices[i] << " ";
std::cout << "cols: "; std::cout << std::endl;
for (int i = 0; i < colIndices.size(); i++) std::cout << "cols: ";
std::cout << colIndices[i] << " "; for (int i = 0; i < colIndices.size(); i++)
std::cout << std::endl; std::cout << colIndices[i] << " ";
std::cout << std::endl;
}
#endif #endif
for (int i = 0; i < nRow; i++) { for (int i = 0; i < nRow; i++) {
...@@ -229,7 +231,7 @@ namespace AMDiS { ...@@ -229,7 +231,7 @@ namespace AMDiS {
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 (condition->applyBoundaryCondition()) { if (condition->applyBoundaryCondition()) {
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
if ((*rankDofs)[rowIndices[i]]) if ((*rankDofs)[rowIndices[i]])
...@@ -241,8 +243,9 @@ namespace AMDiS { ...@@ -241,8 +243,9 @@ namespace AMDiS {
} else { } else {
for (int j = 0; j < nCol; j++) { for (int j = 0; j < nCol; j++) {
DegreeOfFreedom col = colIndices[j]; DegreeOfFreedom col = colIndices[j];
if (fabs(elMat[i][j]) > 1e-10) if (fabs(elMat[i][j]) > 1e-10) {
ins[row][col] += elMat[i][j]; ins[row][col] += elMat[i][j];
}
} }
} }
} }
...@@ -259,7 +262,9 @@ namespace AMDiS { ...@@ -259,7 +262,9 @@ namespace AMDiS {
{} {}
void DOFMatrix::assemble(double factor, ElInfo *elInfo, const BoundaryType *bound) void DOFMatrix::assemble(double factor,
ElInfo *elInfo,
const BoundaryType *bound)
{ {
FUNCNAME("DOFMatrix::assemble()"); FUNCNAME("DOFMatrix::assemble()");
...@@ -278,7 +283,9 @@ namespace AMDiS { ...@@ -278,7 +283,9 @@ namespace AMDiS {
} }
void DOFMatrix::assemble(double factor, ElInfo *elInfo, const BoundaryType *bound, void DOFMatrix::assemble(double factor,
ElInfo *elInfo,
const BoundaryType *bound,
Operator *op) Operator *op)
{ {
FUNCNAME("DOFMatrix::assemble()"); FUNCNAME("DOFMatrix::assemble()");
......
...@@ -1356,6 +1356,10 @@ namespace AMDiS { ...@@ -1356,6 +1356,10 @@ namespace AMDiS {
{ {
FUNCNAME("ProblemStat::addMatrixOperator()"); FUNCNAME("ProblemStat::addMatrixOperator()");
TEST_EXIT(i < nComponents && j < nComponents)
("Cannot add matrix operator at position %d/%d. The stationary problem has only %d components!\n",
i, j, nComponents);
TEST_EXIT(!boundaryConditionSet) TEST_EXIT(!boundaryConditionSet)
("Do not add operators after boundary conditions were set!\n"); ("Do not add operators after boundary conditions were set!\n");
...@@ -1402,6 +1406,10 @@ namespace AMDiS { ...@@ -1402,6 +1406,10 @@ namespace AMDiS {
{ {
FUNCNAME("ProblemStat::addVectorOperator()"); FUNCNAME("ProblemStat::addVectorOperator()");
TEST_EXIT(i < nComponents)
("Cannot add vector operator at position %d. The stationary problem has only %d components!\n",
i, nComponents);
TEST_EXIT(!boundaryConditionSet) TEST_EXIT(!boundaryConditionSet)
("Do not add operators after boundary conditions were set!\n"); ("Do not add operators after boundary conditions were set!\n");
......
...@@ -64,17 +64,13 @@ namespace AMDiS { ...@@ -64,17 +64,13 @@ namespace AMDiS {
/// Destructor /// Destructor
virtual ~SubAssembler() {} virtual ~SubAssembler() {}
/** \brief /// Calculates the element matrix for elInfo and adds it to mat. Memory for
* Calculates the element matrix for elInfo and adds it to mat. Memory /// mat must be provided by the caller.
* for mat must be provided by the caller.
*/
virtual void calculateElementMatrix(const ElInfo *elInfo, virtual void calculateElementMatrix(const ElInfo *elInfo,
ElementMatrix& mat) = 0; ElementMatrix& mat) = 0;
/** \brief /// Calculates the element vector for elInfo and adds it to vec. Memory for
* Calculates the element vector for elInfo and adds it to vec. Memory /// vec must be provided by the caller.
* for vec must be provided by the caller.
*/
virtual void calculateElementVector(const ElInfo *elInfo, virtual void calculateElementVector(const ElInfo *elInfo,
ElementVector& vec) = 0; ElementVector& vec) = 0;
...@@ -104,10 +100,8 @@ namespace AMDiS { ...@@ -104,10 +100,8 @@ namespace AMDiS {
WorldVector<double>* getCoordsAtQPs(const ElInfo* elInfo, WorldVector<double>* getCoordsAtQPs(const ElInfo* elInfo,
Quadrature *quad = NULL); Quadrature *quad = NULL);
/** \brief /// DOFVector dv evaluated at quadrature points.
* DOFVector dv evaluated at quadrature points. /// Used by \ref OperatorTerm::initElement().
* Used by \ref OperatorTerm::initElement().
*/
template<typename T> template<typename T>
void getVectorAtQPs(DOFVectorBase<T>* dv, void getVectorAtQPs(DOFVectorBase<T>* dv,
const ElInfo* elInfo, const ElInfo* elInfo,
...@@ -122,10 +116,8 @@ namespace AMDiS { ...@@ -122,10 +116,8 @@ namespace AMDiS {
Quadrature *quad, Quadrature *quad,
mtl::dense_vector<T>& vecAtQPs); mtl::dense_vector<T>& vecAtQPs);
/** \brief /// Gradients of DOFVector dv evaluated at quadrature points.
* Gradients of DOFVector dv evaluated at quadrature points. /// Used by \ref OperatorTerm::initElement().
* Used by \ref OperatorTerm::initElement().
*/
template<typename T> template<typename T>
void getGradientsAtQPs( DOFVectorBase<T>* dv, void getGradientsAtQPs( DOFVectorBase<T>* dv,
const ElInfo* elInfo, const ElInfo* elInfo,
...@@ -139,11 +131,10 @@ namespace AMDiS { ...@@ -139,11 +131,10 @@ namespace AMDiS {
Quadrature *quad, Quadrature *quad,
mtl::dense_vector<typename GradientType<T>::type>& grdAtQPs); mtl::dense_vector<typename GradientType<T>::type>& grdAtQPs);
/** \brief /// The comp'th component of the derivative of DOFVector dv evaluated at
* The comp'th component of the derivative of DOFVector dv evaluated at quadrature points. /// quadrature points. Used by \ref OperatorTerm::initElement().
* Used by \ref OperatorTerm::initElement(). /// Attention: not caching at the moment! Using cache if gradients for read
* Attention: not caching at the moment! Using cache if gradients for read but not for write /// but not for write.
*/
template<typename T> template<typename T>
void getDerivativeAtQPs(DOFVectorBase<T>* dv, void getDerivativeAtQPs(DOFVectorBase<T>* dv,
const ElInfo* elInfo, const ElInfo* elInfo,
...@@ -159,11 +150,9 @@ namespace AMDiS { ...@@ -159,11 +150,9 @@ namespace AMDiS {
int comp, int comp,
mtl::dense_vector<T>& grdAtQPs); mtl::dense_vector<T>& grdAtQPs);
/** \brief /// Called once for each ElInfo when \ref calculateElementMatrix() or
* Called once for each ElInfo when \ref calculateElementMatrix() or /// \ref calculateElementVector() is called for the first time for this
* \ref calculateElementVector() is called for the first time for this /// Element.
* Element.
*/
virtual void initElement(const ElInfo *smallElInfo, virtual void initElement(const ElInfo *smallElInfo,
const ElInfo *largeElInfo = NULL, const ElInfo *largeElInfo = NULL,
Quadrature *quad = NULL); Quadrature *quad = NULL);
...@@ -202,16 +191,12 @@ namespace AMDiS { ...@@ -202,16 +191,12 @@ namespace AMDiS {
/// Column FiniteElemSpace. /// Column FiniteElemSpace.
const FiniteElemSpace *colFeSpace; const FiniteElemSpace *colFeSpace;
/** \brief /// Number of rows of the element matrix and length of the element
* Number of rows of the element matrix and length of the element /// vector. Is equal to the number of row basis functions
* vector. Is equal to the number of row basis functions
*/
int nRow; int nRow;
/** \brief /// Number of columns of the element matrix. Is equal to the number
* Number of columns of the element matrix. Is equal to the number /// of column basis functions
* of column basis functions
*/
int nCol; int nCol;
/// Used for \ref getVectorAtQPs() and \ref getGradientsAtQPs(). /// Used for \ref getVectorAtQPs() and \ref getGradientsAtQPs().
...@@ -230,16 +215,15 @@ namespace AMDiS { ...@@ -230,16 +215,15 @@ namespace AMDiS {
std::map<const void*, ValuesAtQPs* > cachedValuesAtQPs; std::map<const void*, ValuesAtQPs* > cachedValuesAtQPs;
std::map<const void*, ValuesAtQPs* > cachedGradientsAtQPs; std::map<const void*, ValuesAtQPs* > cachedGradientsAtQPs;
/** \brief /// Set and updated by \ref initElement() for each ElInfo.
* Set and updated by \ref initElement() for each ElInfo. /// coordsAtQPs[i] points to the coordinates of the i-th quadrature point.
* coordsAtQPs[i] points to the coordinates of the i-th quadrature point.
*/
WorldVector<double> *coordsAtQPs; WorldVector<double> *coordsAtQPs;
/// Used for \ref getCoordsAtQPs(). /// Used for \ref getCoordsAtQPs().
bool coordsValid; bool coordsValid;
/// Used for \ref getCoordsAtQP(). Stores the number of allocated WorldVectors. /// Used for \ref getCoordsAtQP(). Stores the number of allocated
/// WorldVectors.
int coordsNumAllocated; int coordsNumAllocated;
/// Quadrature object to be used for assembling. /// Quadrature object to be used for assembling.
......
...@@ -376,6 +376,48 @@ namespace AMDiS { ...@@ -376,6 +376,48 @@ namespace AMDiS {
return; return;
break; break;
case EDGE: case EDGE:
{
// === Create boundary information objects for children elements. ===
BoundaryObject nextBound0 = bound;
prepareNextBound(nextBound0, 0);
BoundaryObject nextBound1 = bound;
prepareNextBound(nextBound1, 1);
// === Check for boundary on children elements. ===
if ((nextBound0.ithObj >= 0 || nextBound1.ithObj >= 0) && child[0]) {
// So, the edge is contained in at least on of the children and the
// element is also refined. Then we have go down further in refinement
// hierarchie.
if (bound.reverseMode) {
if (nextBound1.ithObj >= 0)
child[1]->getHigherOrderDofs(feSpace, nextBound1, dofs);
if (nextBound0.ithObj >= 0)
child[0]->getHigherOrderDofs(feSpace, nextBound0, dofs);
} else {
if (nextBound0.ithObj >= 0)
child[0]->getHigherOrderDofs(feSpace, nextBound0, dofs);
if (nextBound1.ithObj >= 0)
child[1]->getHigherOrderDofs(feSpace, nextBound1, dofs);
}
} else {
// Either the edge is not contained in further refined children, or
// the element is not refined further on this edge. Then we can get
// all the DOFs on this edge.
ElementDofIterator elDofIter(feSpace, true);
elDofIter.reset(this);
do {
if (elDofIter.getCurrentPos() == 1 &&
elDofIter.getCurrentElementPos() == bound.ithObj)
dofs.push_back(elDofIter.getDofPtr());
} while(elDofIter.next());
}
}
break; break;
case FACE: case FACE:
{ {
...@@ -403,10 +445,8 @@ namespace AMDiS { ...@@ -403,10 +445,8 @@ namespace AMDiS {
elDofIter.reset(this); elDofIter.reset(this);
do { do {
if (elDofIter.getCurrentPos() == 2 && if (elDofIter.getCurrentPos() == 2 &&
elDofIter.getCurrentElementPos() == bound.ithObj) { elDofIter.getCurrentElementPos() == bound.ithObj)
ERROR_EXIT("Check this, if it will really work!\n");
dofs.push_back(elDofIter.getDofPtr()); dofs.push_back(elDofIter.getDofPtr());
}
} while(elDofIter.next()); } while(elDofIter.next());
} }
} }
......
...@@ -324,10 +324,10 @@ namespace AMDiS { ...@@ -324,10 +324,10 @@ namespace AMDiS {
static const int sideOfChild[3][2][4]; static const int sideOfChild[3][2][4];
/** \brief /** \brief
* edgeOfChild[elType][i][j] is the local edge number of the j-th edge within the * edgeOfChild[elType][i][j] is the local edge number of the j-th edge within
* i-th children of an element of elType. If the value is -1, the edge is not * the i-th children of an element of elType. If the value is -1, the edge is
* included in the element's child. Note that the 0 edge is included in both * not included in the element's child. Note that the 0 edge is included in
* children only by its half. * both children only by its half.
*/ */
static const int edgeOfChild[3][2][6]; static const int edgeOfChild[3][2][6];
......
...@@ -169,7 +169,7 @@ namespace AMDiS { ...@@ -169,7 +169,7 @@ namespace AMDiS {
break; break;
default: default:
ERROR_EXIT("Should never happen!\n"); ERROR_EXIT("Should never happen!\n");
} }
} }
......
...@@ -77,9 +77,9 @@ namespace AMDiS { ...@@ -77,9 +77,9 @@ namespace AMDiS {
newAssembler = new StandardZOA(op, assembler, quad); newAssembler = new StandardZOA(op, assembler, quad);
} else { } else {
if (pwConst) if (pwConst)
newAssembler = new PrecalcZOA(op, assembler, quad); newAssembler = new PrecalcZOA(op, assembler, quad);
else else
newAssembler = new FastQuadZOA(op, assembler, quad); newAssembler = new FastQuadZOA(op, assembler, quad);
} }
subAssemblers->push_back(newAssembler); subAssemblers->push_back(newAssembler);
......
...@@ -85,9 +85,12 @@ namespace AMDiS { ...@@ -85,9 +85,12 @@ namespace AMDiS {
periodicFaces[make_pair(face0, face1)] = elInfo->getBoundary(i); periodicFaces[make_pair(face0, face1)] = elInfo->getBoundary(i);
/// Add all three vertices of the face to be periodic. /// Add all three vertices of the face to be periodic.
periodicVertices[make_pair(face0.get<0>(), face1.get<0>())] = boundaryType; periodicVertices[make_pair(face0.get<0>(), face1.get<0>())] =
periodicVertices[make_pair(face0.get<1>(), face1.get<1>())] = boundaryType; boundaryType;
periodicVertices[make_pair(face0.get<2>(), face1.get<2>())] = boundaryType; periodicVertices[make_pair(face0.get<1>(), face1.get<1>())] =
boundaryType;
periodicVertices[make_pair(face0.get<2>(), face1.get<2>())] =
boundaryType;
periodicDofAssoc[face0.get<0>()].insert(boundaryType); periodicDofAssoc[face0.get<0>()].insert(boundaryType);
periodicDofAssoc[face0.get<1>()].insert(boundaryType); periodicDofAssoc[face0.get<1>()].insert(boundaryType);
...@@ -131,8 +134,8 @@ namespace AMDiS { ...@@ -131,8 +134,8 @@ namespace AMDiS {
TEST_EXIT_DBG(mesh)("Mesh not set!\n"); TEST_EXIT_DBG(mesh)("Mesh not set!\n");
// === Return, if there are no periodic vertices, i.e., there are no no === // === Return, if there are no periodic vertices, i.e., there are no ===
// === periodic boundaries in the mesh. === // === periodic boundaries in the mesh. ===
if (periodicVertices.size() == 0) if (periodicVertices.size() == 0)
return; return;
...@@ -141,9 +144,9 @@ namespace AMDiS { ...@@ -141,9 +144,9 @@ namespace AMDiS {
// === Get all vertex DOFs that have multiple periodic associations. === // === Get all vertex DOFs that have multiple periodic associations. ===
// We group all vertices together, that have either two or three periodic // We group all vertices together, that have either two or three periodic
// associations. For rectangular domains in 2D, the four corner vertices have all // associations. For rectangular domains in 2D, the four corner vertices have
// two periodic associations. For box domains in 3D, the eight corner vertices // all two periodic associations. For box domains in 3D, the eight corner
// have all three periodic associations. // vertices have all three periodic associations.
vector<DegreeOfFreedom> multPeriodicDof2, multPeriodicDof3; vector<DegreeOfFreedom> multPeriodicDof2, multPeriodicDof3;
for (map<DegreeOfFreedom, std::set<BoundaryType> >::iterator it = periodicDofAssoc.begin(); for (map<DegreeOfFreedom, std::set<BoundaryType> >::iterator it = periodicDofAssoc.begin();
......
...@@ -126,11 +126,11 @@ namespace AMDiS { ...@@ -126,11 +126,11 @@ namespace AMDiS {
/** \brief /** \brief
* Creates final data of the periodic boundaries. Must be called after all * Creates final data of the periodic boundaries. Must be called after all
* elements of the mesh are added to the object database. Then this functions * elements of the mesh are added to the object database. Then this functions
* search for interectly connected vertices in periodic boundaries. This is only * search for indirectly connected vertices in periodic boundaries. This is
* the case, if there are more than one boundary conditions. Then, e.g., in 2D, * only the case, if there are more than one boundary conditions. Then, e.g.,
* all edges of a square are iterectly connected. In 3D, if the macro mesh is a * in 2D, all edges of a square are iterectly connected. In 3D, if the macro
* box, all eight vertex nodes and always four of the 12 edges are iterectly * mesh is a box, all eight vertex nodes and always four of the 12 edges are
* connected. * indirectly connected.
*/ */
void createPeriodicData(const FiniteElemSpace *feSpace); void createPeriodicData(const FiniteElemSpace *feSpace);
......
...@@ -289,17 +289,17 @@ namespace AMDiS { ...@@ -289,17 +289,17 @@ namespace AMDiS {
// We have to remove the VertexVectors, which contain periodic assoiciations, // We have to remove the VertexVectors, which contain periodic assoiciations,
// because they are not valid anymore after some macro elements have been removed // because they are not valid anymore after some macro elements have been
// and the corresponding DOFs were deleted. // removed and the corresponding DOFs were deleted.
for (map<BoundaryType, VertexVector*>::iterator it = mesh->getPeriodicAssociations().begin(); for (map<BoundaryType, VertexVector*>::iterator it = mesh->getPeriodicAssociations().begin();
it != mesh->getPeriodicAssociations().end(); ++it) it != mesh->getPeriodicAssociations().end(); ++it)
const_cast<DOFAdmin&>(mesh->getDofAdmin(0)).removeDOFContainer(dynamic_cast<DOFContainer*>(it->second)); const_cast<DOFAdmin&>(mesh->getDofAdmin(0)).removeDOFContainer(dynamic_cast<DOFContainer*>(it->second));
updateLocalGlobalNumbering(); updateLocalGlobalNumbering();
// === In 3D we have to make some test, if the resulting mesh is valid. If === // === In 3D we have to make some test, if the resulting mesh is valid. ===
// === it is not valid, there is no possiblity yet to fix this problem, just === // === If it is not valid, there is no possiblity yet to fix this ===
// === exit with an error message. === // === problem, just exit with an error message. ===
check3dValidMesh(); check3dValidMesh();
...@@ -320,7 +320,7 @@ namespace AMDiS { ...@@ -320,7 +320,7 @@ namespace AMDiS {
// === Create periodic DOF mapping, if there are periodic boundaries. === // === Create periodic DOF mapping, if there are periodic boundaries. ===
createPeriodicMap(feSpaces[0]); createPeriodicMap();
#if (DEBUG != 0) #if (DEBUG != 0)
ParallelDebug::testPeriodicBoundary(*this); ParallelDebug::testPeriodicBoundary(*this);
...@@ -335,25 +335,20 @@ namespace AMDiS { ...@@ -335,25 +335,20 @@ namespace AMDiS {
refineManager->globalRefine(mesh, globalRefinement); refineManager->globalRefine(mesh, globalRefinement);
updateLocalGlobalNumbering(); updateLocalGlobalNumbering();
// === Update periodic mapping, if there are periodic boundaries. === // === Update periodic mapping, if there are periodic boundaries. ===
createPeriodicMap(feSpaces[0]); createPeriodicMap();
#if (DEBUG != 0) #if (DEBUG != 0)
ParallelDebug::testPeriodicBoundary(*this); ParallelDebug::testPeriodicBoundary(*this);
#endif #endif
} }
// Set DOF rank information to all matrices and vectors.
/// === Set DOF rank information to all matrices and vectors. ===
setRankDofs(); setRankDofs();
// Remove periodic boundary conditions in sequential problem definition.
// === Remove periodic boundary conditions in sequential problem definition. ===
removePeriodicBoundaryConditions(); removePeriodicBoundaryConditions();
initialized = true; initialized = true;
...@@ -364,6 +359,9 @@ namespace AMDiS { ...@@ -364,6 +359,9 @@ namespace AMDiS {
{ {
FUNCNAME("MeshDistributor::addProblemStat()"); FUNCNAME("MeshDistributor::addProblemStat()");
TEST_EXIT_DBG(probStat->getFeSpaces().size())
("No FE spaces in stationary problem!\n");
// === Add FE spaces from stationary problem to mesh distributor. === // === Add FE spaces from stationary problem to mesh distributor. ===
for (unsigned int i = 0; i < probStat->getFeSpaces().size(); i++) { for (unsigned int i = 0; i < probStat->getFeSpaces().size(); i++) {
...@@ -845,21 +843,20 @@ namespace AMDiS { ...@@ -845,21 +843,20 @@ namespace AMDiS {
#endif #endif
// === Because the mesh has been changed, update the DOF numbering and mappings. === // Because the mesh has been changed, update the DOF numbering and mappings.
updateLocalGlobalNumbering(); updateLocalGlobalNumbering();
// Update periodic mapping, if there are periodic boundaries.
createPeriodicMap();
// === Update periodic mapping, if there are periodic boundaries. ===
createPeriodicMap(feSpaces[0]);