Commit 566d2c61 authored by Thomas Witkowski's avatar Thomas Witkowski

First update of parallel AMDiS with mixed finite elements. Not tested, and there ARE BUGS.

parent 15d34568
......@@ -150,7 +150,7 @@ namespace AMDiS {
return name;
}
/// Returns \ref nDof[i], i.e., the number of dofs for the position i.
/// Returns \ref nDof[i], i.e., the number of DOFs for the position i.
inline const int getNumberOfDofs(int i) const
{
return nDof[i];
......
......@@ -492,7 +492,10 @@ namespace AMDiS {
int nnzPerRow;
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
std::map<DegreeOfFreedom, bool> *rankDofs;
/// Stores for the DOFs of the row FE spaces whether they are owned by the
/// rank or not. This is used to ensure that Dirichlet BC is handled
/// correctly in parallel computations.
std::map<DegreeOfFreedom, bool> *rankDofs;
#endif
/// Inserter object: implemented as pointer, allocated and deallocated as needed
......
......@@ -25,7 +25,8 @@ namespace AMDiS {
namespace debug {
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
void writeLocalElementDofs(int rank, int elIdx, FiniteElemSpace *feSpace)
void writeLocalElementDofs(int rank, int elIdx,
const FiniteElemSpace *feSpace)
{
using boost::lexical_cast;
......@@ -37,7 +38,8 @@ namespace AMDiS {
}
void writeDofMesh(int rank, DegreeOfFreedom dof, FiniteElemSpace *feSpace)
void writeDofMesh(int rank, DegreeOfFreedom dof,
const FiniteElemSpace *feSpace)
{
using boost::lexical_cast;
......@@ -50,7 +52,7 @@ namespace AMDiS {
}
void writeMesh(FiniteElemSpace *feSpace, int rank, std::string filename)
void writeMesh(const FiniteElemSpace *feSpace, int rank, std::string filename)
{
using boost::lexical_cast;
......@@ -64,7 +66,7 @@ namespace AMDiS {
#endif
void writeDofIndexMesh(FiniteElemSpace *feSpace)
void writeDofIndexMesh(const FiniteElemSpace *feSpace)
{
DOFVector<double> tmp(feSpace, "tmp");
DOFIterator<double> it(&tmp, USED_DOFS);
......@@ -74,7 +76,7 @@ namespace AMDiS {
}
void colorEdgeInMesh(FiniteElemSpace *feSpace,
void colorEdgeInMesh(const FiniteElemSpace *feSpace,
Element *el,
int localEdgeNo,
std::string filename)
......@@ -190,7 +192,8 @@ namespace AMDiS {
}
Element* getDofIndexElement(FiniteElemSpace *feSpace, DegreeOfFreedom dof)
Element* getDofIndexElement(const FiniteElemSpace *feSpace,
DegreeOfFreedom dof)
{
const BasisFunction* basFcts = feSpace->getBasisFcts();
int nBasFcts = basFcts->getNumber();
......@@ -322,7 +325,7 @@ namespace AMDiS {
}
void printInfoByDof(FiniteElemSpace *feSpace, DegreeOfFreedom dof)
void printInfoByDof(const FiniteElemSpace *feSpace, DegreeOfFreedom dof)
{
FUNCNAME("debug::printInfoByDof()");
......@@ -387,7 +390,7 @@ namespace AMDiS {
}
void printAllDofCoords(FiniteElemSpace *feSpace)
void printAllDofCoords(const FiniteElemSpace *feSpace)
{
FUNCNAME("printAllDofCoords()");
......@@ -405,7 +408,8 @@ namespace AMDiS {
}
void getAllDofs(FiniteElemSpace *feSpace, std::set<const DegreeOfFreedom*>& dofs)
void getAllDofs(const FiniteElemSpace *feSpace,
std::set<const DegreeOfFreedom*>& dofs)
{
FUNCNAME("getAllDofs()");
......@@ -793,7 +797,7 @@ namespace AMDiS {
}
void testDofsByCoords(FiniteElemSpace *feSpace,
void testDofsByCoords(const FiniteElemSpace *feSpace,
DofContainer &dofs0, DofContainer &dofs1)
{
FUNCNAME("debug::testDofsByCoords()");
......
......@@ -41,9 +41,13 @@ namespace AMDiS {
typedef std::map<int, DofContainer> ElementIdxToDofs;
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
void writeLocalElementDofs(int rank, int elIdx, FiniteElemSpace *feSpace);
void writeLocalElementDofs(int rank,
int elIdx,
const FiniteElemSpace *feSpace);
void writeMesh(FiniteElemSpace *feSpace, int rank, std::string filename);
void writeMesh(const FiniteElemSpace *feSpace,
int rank,
std::string filename);
/** \brief
* Writes a vtu file with the mesh, where all DOFs are set to zero, and only
......@@ -55,7 +59,9 @@ namespace AMDiS {
* \param[in] dof Defines the DOF, which value is set to one in the mesh file.
* \param[in] feSpace The FE space to be used.
*/
void writeDofMesh(int rank, DegreeOfFreedom dof, FiniteElemSpace *feSpace);
void writeDofMesh(int rank,
DegreeOfFreedom dof,
const FiniteElemSpace *feSpace);
#endif
/** \brief
......@@ -64,9 +70,9 @@ namespace AMDiS {
*
* \param[in] feSpace The FE space to be used.
*/
void writeDofIndexMesh(FiniteElemSpace *feSpace);
void writeDofIndexMesh(const FiniteElemSpace *feSpace);
void colorEdgeInMesh(FiniteElemSpace *feSpace,
void colorEdgeInMesh(const FiniteElemSpace *feSpace,
Element *el,
int localEdgeNo,
std::string filename);
......@@ -92,7 +98,8 @@ namespace AMDiS {
Mesh *mesh,
int elIndex);
Element* getDofIndexElement(FiniteElemSpace *feSpace, DegreeOfFreedom dof);
Element* getDofIndexElement(const FiniteElemSpace *feSpace,
DegreeOfFreedom dof);
Element* getLevel0ParentElement(Mesh *mesh, Element *el);
......@@ -108,13 +115,14 @@ namespace AMDiS {
void printElementCoords(const FiniteElemSpace *feSpace, Element *el);
void printInfoByDof(FiniteElemSpace *feSpace, DegreeOfFreedom dof);
void printInfoByDof(const FiniteElemSpace *feSpace, DegreeOfFreedom dof);
void printMatValuesStatistics(Matrix<DOFMatrix*> *mat);
void printAllDofCoords(FiniteElemSpace *feSpace);
void printAllDofCoords(const FiniteElemSpace *feSpace);
void getAllDofs(FiniteElemSpace *feSpace, std::set<const DegreeOfFreedom*>& dofs);
void getAllDofs(const FiniteElemSpace *feSpace,
std::set<const DegreeOfFreedom*>& dofs);
/** \brief
* Creates a text file storing the value of a sparse matrix. Each line of the file
......@@ -200,7 +208,7 @@ namespace AMDiS {
const DegreeOfFreedom* dof3,
DofContainer &vec);
void testDofsByCoords(FiniteElemSpace *feSpace,
void testDofsByCoords(const FiniteElemSpace *feSpace,
DofContainer &dofs0, DofContainer &dofs1);
}
......
......@@ -625,7 +625,7 @@ namespace AMDiS {
}
void Element::getAllDofs(FiniteElemSpace* feSpace,
void Element::getAllDofs(const FiniteElemSpace* feSpace,
BoundaryObject bound,
DofContainer& dofs)
{
......
......@@ -421,7 +421,7 @@ namespace AMDiS {
* which all vertex dofs are assembled.
* \param[out] dofs List of dofs, where the result is stored.
*/
virtual void getNodeDofs(FiniteElemSpace* feSpace,
virtual void getNodeDofs(const FiniteElemSpace* feSpace,
BoundaryObject bound,
DofContainer& dofs) const = 0;
......@@ -436,11 +436,11 @@ namespace AMDiS {
* all non vertex dofs are assembled.
* \param[out] dofs All dofs are put to this dof list.
*/
virtual void getHigherOrderDofs(FiniteElemSpace* feSpace,
virtual void getHigherOrderDofs(const FiniteElemSpace* feSpace,
BoundaryObject bound,
DofContainer& dofs) const = 0;
void getAllDofs(FiniteElemSpace* feSpace,
void getAllDofs(const FiniteElemSpace* feSpace,
BoundaryObject bound,
DofContainer& dofs);
......
......@@ -36,14 +36,15 @@ namespace AMDiS {
// Get geo index of vertices in the given dimension.
posIndex = INDEX_OF_DIM(pos, dim);
// Get number of dofs per vertex (should be one in all cases).
// Get number of DOFs per vertex (should be one in all cases).
nDofs = admin->getNumberOfDofs(posIndex);
TEST_EXIT_DBG(nDofs != 0)("Mh, I've to think about this situation!\n");
// Calculate displacement. Is used if there is more than one dof admin on the mesh.
// Calculate displacement. Is used if there is more than one DOF admin
// on the mesh.
n0 = admin->getNumberOfPreDofs(posIndex);
// Get first dof index position for vertices.
// Get first DOF index position for vertices.
node0 = mesh->getNode(posIndex);
// Get number of vertices in this dimension.
nElements = Global::getGeo(posIndex, mesh->getDim());
......@@ -55,11 +56,11 @@ namespace AMDiS {
bool ElementDofIterator::next()
{
// First iterate over the dofs of one element (vertex, edge, face).
// First iterate over the DOFs of one element (vertex, edge, face).
dofPos++;
if (dofPos >= nDofs) {
// We are finished with all dofs of on element. Go to the next one.
// We are finished with all DOFs of on element. Go to the next one.
dofPos = 0;
elementPos++;
......@@ -72,33 +73,35 @@ namespace AMDiS {
return false;
// Increase position, i.e., go from vertices to edges to faces and search
// for the next position with dofs.
// for the next position with DOFs.
do {
pos++;
// Get geo index posistion.
posIndex = INDEX_OF_DIM(pos, dim);
// Get number of dofs in this position.
// Get number of DOFs in this position.
nDofs = admin->getNumberOfDofs(posIndex);
} while (nDofs == 0 && pos < dim);
if (nDofs > 0 && pos <= dim) {
// We have found on more position with dofs.
// We have found on more position with DOFs.
// Get number of elements in this position, i.e, the number of vertices,.
// edges and faces in the given dimension.
// Get number of elements in this position, i.e, the number of
// vertices, edges and faces in the given dimension.
nElements = Global::getGeo(posIndex, dim);
// Calculate displacement. Is used if there is more than one dof admin on the mesh.
// Calculate displacement. Is used if there is more than one DOF
// admin on the mesh.
n0 = admin->getNumberOfPreDofs(posIndex);
// Get first dof index position for the geo index position.
// Get first DOF index position for the geo index position.
node0 = mesh->getNode(posIndex);
if (inOrder)
orderPosition = basisFcts->orderOfPositionIndices(element, posIndex, 0);
orderPosition =
basisFcts->orderOfPositionIndices(element, posIndex, 0);
} else {
// That's all, we jave traversed all dofs of the mesh element.
// That's all, we jave traversed all DOFs of the mesh element.
return false;
}
......
......@@ -30,7 +30,7 @@
namespace AMDiS {
/** \brief
* This class implements an iterator to iterate over all dofs of one element
* This class implements an iterator to iterate over all DOFs of one element
* independet of dimension and the degree of basis functions.
*
* Should be used in the following way:
......@@ -57,12 +57,12 @@ namespace AMDiS {
/// Start a new traverse with the given element.
void reset(const Element* el);
/// Go to next dof. Returns false, if there is dof anymore.
/// Go to next DOF. Returns false, if there is no DOF anymore.
bool next();
bool nextStrict();
/// Returns the dof index of the current dof.
/// Returns the DOF index of the current DOF.
inline DegreeOfFreedom getDof()
{
if (inOrder)
......@@ -71,7 +71,7 @@ namespace AMDiS {
return dofs[node0 + elementPos][n0 + dofPos];
}
/// Returns a pointer to the current dof.
/// Returns a pointer to the current DOF.
inline const DegreeOfFreedom* getDofPtr()
{
if (inOrder)
......@@ -80,13 +80,22 @@ namespace AMDiS {
return &dofs[node0 + elementPos][n0 + dofPos];
}
/// Returns \ref pos, the current position (vertex, edge, face) of the traverse.
/// Returns a pointer to the starting position of the current DOF
/// array. Makes only sence, if \ref nextStrict() is used for traverse.
inline const DegreeOfFreedom* getStartDof()
{
return dofs[node0 + elementPos];
}
/// Returns \ref pos, the current position (vertex, edge, face) of
/// the traverse.
inline int getCurrentPos()
{
return pos;
}
/// Returns \ref elementPos, the number of vertex, edge or face that is traversed.
/// Returns \ref elementPos, the number of vertex, edge or face that
/// is traversed.
inline int getCurrentElementPos()
{
return elementPos;
......@@ -99,12 +108,12 @@ namespace AMDiS {
protected:
/// The dof admin for which dofs should be traversed.
/// The DOF admin for which DOFs should be traversed.
const DOFAdmin* admin;
const BasisFunction* basisFcts;
/// Pointer to the dofs that should be traversed.
/// Pointer to the DOFs that should be traversed.
const DegreeOfFreedom **dofs;
/// Mesh on which the element is defined.
......@@ -126,32 +135,28 @@ namespace AMDiS {
GeoIndex posIndex;
/** \brief
* Number of dofs at the current traverse position. Examples: independent of
* dimension and degree of basis functions there is only one dof per vertex.
* But in 2d and with 3rd degree lagrange basis functions there are two
* dofs per edge.
* Number of DOFs at the current traverse position. Examples: independent of
* dimension and degree of basis functions there is only one DOF per vertex.
* But in 2D and with 3rd degree lagrange basis functions there are two
* DOFs per edge.
*/
int nDofs;
/** \brief
* Displacement of dof indices. Used of more than one dof admin is defined
* on the mesh.
*/
/// Displacement of DOF indices. Used if more than one DOF admin is defined
/// on the mesh.
int n0;
/// Dof index of the first dof at this geo index position.
/// DOF index of the first DOF at this geo index position.
int node0;
/** \brief
* Number of elements in the current geo position. Examples: 3 vertices in 2d,
* 1 face in 2d, 4 faces in 3d, etc.
*/
/// Number of elements in the current geo position. Examples: 3 vertices in
/// 2d, 1 face in 2d, 4 faces in 3d, etc.
int nElements;
/// Current element, i.e., ith vertex, edge or face, that is traversed.
int elementPos;
/// Currrent dof that is traversed on the current element;
/// Currrent DOF that is traversed on the current element;
int dofPos;
};
}
......
......@@ -173,13 +173,13 @@ namespace AMDiS {
return "Line";
}
void getNodeDofs(FiniteElemSpace*, BoundaryObject, DofContainer&) const
void getNodeDofs(const FiniteElemSpace*, BoundaryObject, DofContainer&) const
{
FUNCNAME("Line::getNodeDofs()");
ERROR_EXIT("Not yet implemented!\n");
}
void getHigherOrderDofs(FiniteElemSpace*, BoundaryObject, DofContainer&) const
void getHigherOrderDofs(const FiniteElemSpace*, BoundaryObject, DofContainer&) const
{
FUNCNAME("Line::getHigherOrderDofs()");
ERROR_EXIT("Not yet implemented!\n");
......
......@@ -271,21 +271,32 @@ namespace AMDiS {
void Mesh::removeMacroElements(std::set<MacroElement*>& macros,
const FiniteElemSpace *feSpace)
vector<const FiniteElemSpace*>& feSpaces)
{
FUNCNAME("Mesh::removeMacroElement()");
typedef map<const DegreeOfFreedom*, std::set<MacroElement*> > DofElMap;
typedef map<const DegreeOfFreedom*, GeoIndex> DofPosMap;
TEST_EXIT(admin.size() == 1)("Not yet implemented for multiple admins!\n");
TEST_EXIT(admin[0])("There is something wrong!\n");
TEST_EXIT(feSpaces.size() > 0)("Should not happen!\n");
// === Search for the FE space with the highest degree of polynomials. ===
// === Using this FE space ensures that deleting DOFs defined on it, ===
// === also DOFs of lower order FE spaces will be deleted correctly. ===
const FiniteElemSpace *feSpace = feSpaces[0];
for (unsigned int i = 1; i < feSpaces.size(); i++)
if (feSpaces[i]->getBasisFcts()->getDegree() >
feSpace->getBasisFcts()->getDegree())
feSpace = feSpaces[i];
// === Determine to all DOFs in mesh the macro elements where the DOF ===
// === is part of. ===
// Map that stores for each dof pointer (which may have a list of dofs)
// all macro element indices that own this dof.
// Map that stores for each DOF pointer (which may have a list of DOFs)
// all macro element indices that own this DOF.
DofElMap dofsOwner;
DofPosMap dofsPosIndex;
......@@ -295,8 +306,8 @@ namespace AMDiS {
while (elInfo) {
elDofIter.reset(elInfo->getElement());
do {
dofsOwner[elDofIter.getDofPtr()].insert(elInfo->getMacroElement());
dofsPosIndex[elDofIter.getDofPtr()] = elDofIter.getPosIndex();
dofsOwner[elDofIter.getStartDof()].insert(elInfo->getMacroElement());
dofsPosIndex[elDofIter.getStartDof()] = elDofIter.getPosIndex();
} while (elDofIter.nextStrict());
elInfo = stack.traverseNext(elInfo);
......@@ -615,15 +626,14 @@ namespace AMDiS {
("ndof too big: ndof = %d, MAX_DOF = %d\n", ndof, MAX_DOF);
for (unsigned int i = 0; i < admin.size(); i++) {
DOFAdmin *localAdmin = admin[i];
int n = localAdmin->getNumberOfDofs(position);
int n0 = localAdmin->getNumberOfPreDofs(position);
int n = admin[i]->getNumberOfDofs(position);
int n0 = admin[i]->getNumberOfPreDofs(position);
TEST_EXIT_DBG(n + n0 <= ndof)
("n = %d, n0 = %d too large: ndof = %d\n", n, n0, ndof);
for (int j = 0; j < n; j++)
localAdmin->freeDofIndex(dof[n0 + j]);
admin[i]->freeDofIndex(dof[n0 + j]);
}
delete [] dof;
......@@ -974,7 +984,7 @@ namespace AMDiS {
}
void Mesh::getAllDofs(FiniteElemSpace *feSpace,
void Mesh::getAllDofs(const FiniteElemSpace *feSpace,
std::set<const DegreeOfFreedom*>& allDofs)
{
FUNCNAME("Mesh::getAllDofs()");
......@@ -1381,7 +1391,7 @@ namespace AMDiS {
FiniteElemSpace *feSpace =
FiniteElemSpace::provideFeSpace(localAdmin, basFcts, &testMesh, "tmp");
DataCollector dc(feSpace);
DataCollector<> dc(feSpace);
MacroWriter::writeMacro(&dc, newMacroFilename.str().c_str());
if (periodicFilename != "")
......
......@@ -409,12 +409,12 @@ namespace AMDiS {
void addMacroElement(MacroElement* me);
/* \brief
* Removes a set of macro elements from the mesh. This works only for the case,
* that there are no global or local refinements, i.e., all macro elements have
* no children.
* Removes a set of macro elements from the mesh. This works only for the
* case, that there are no global or local refinements, i.e., all macro
* elements have no children.
*/
void removeMacroElements(std::set<MacroElement*>& macros,
const FiniteElemSpace* feSpace);
vector<const FiniteElemSpace*>& feSpaces);
/// Frees the array of DOF pointers (see \ref createDofPtrs)
void freeDofPtrs(DegreeOfFreedom **ptrs);
......@@ -515,7 +515,7 @@ namespace AMDiS {
* @param[in] feSpace The FE space to be used for collecting DOFs.
* @param[out] allDofs The set which is filled with all DOFs.
*/
void getAllDofs(FiniteElemSpace *feSpace,
void getAllDofs(const FiniteElemSpace *feSpace,
std::set<const DegreeOfFreedom*>& allDofs);
/// Returns FILL_ANY_?D
......
......@@ -350,7 +350,7 @@ namespace AMDiS {
return rhs;
}
inline const DOFVector<double>* getRhsVector(int i = 0)
inline DOFVector<double>* getRhsVector(int i = 0)
{
return rhs->getDOFVector(i);
}
......
......@@ -434,7 +434,8 @@ namespace AMDiS {
FUNCNAME("MatrixGradient_SOT::weakEval()");
TEST_EXIT_DBG(f)("No function f!\n");
TEST_EXIT_DBG(gradAtQPs)("Operator was not initialized correctly!\n");
TEST_EXIT_DBG(num_rows(gradAtQPs))
("Operator was not initialized correctly!\n");
int nPoints = grdUhAtQP.size();
for (int iq = 0; iq < nPoints; iq++) {
......
......@@ -200,7 +200,7 @@ namespace AMDiS {
}
void Tetrahedron::getNodeDofs(FiniteElemSpace* feSpace,
void Tetrahedron::getNodeDofs(const FiniteElemSpace* feSpace,
BoundaryObject bound,
DofContainer& dofs) const
{
......@@ -222,7 +222,7 @@ namespace AMDiS {
}
void Tetrahedron::getNodeDofsAtFace(FiniteElemSpace* feSpace,
void Tetrahedron::getNodeDofsAtFace(const FiniteElemSpace* feSpace,
BoundaryObject bound,
DofContainer& dofs) const
{
......@@ -285,7 +285,7 @@ namespace AMDiS {
}
void Tetrahedron::getNodeDofsAtEdge(FiniteElemSpace* feSpace,
void Tetrahedron::getNodeDofsAtEdge(const FiniteElemSpace* feSpace,
BoundaryObject bound,
DofContainer& dofs) const
{
......@@ -359,7 +359,7 @@ namespace AMDiS {
}
void Tetrahedron::getHigherOrderDofs(FiniteElemSpace* feSpace,
void Tetrahedron::getHigherOrderDofs(const FiniteElemSpace* feSpace,
BoundaryObject bound,
DofContainer& dofs) const
{
......
......@@ -136,19 +136,19 @@ namespace AMDiS {
return "Tetrahedron";
}
void getNodeDofs(FiniteElemSpace* feSpace,
void getNodeDofs(const FiniteElemSpace* feSpace,
BoundaryObject bound,
DofContainer& dofs) const;
void getNodeDofsAtFace(FiniteElemSpace* feSpace,
void getNodeDofsAtFace(const FiniteElemSpace* feSpace,
BoundaryObject bound,
DofContainer& dofs) const;
void getNodeDofsAtEdge(FiniteElemSpace* feSpace,
void getNodeDofsAtEdge(const FiniteElemSpace* feSpace,
BoundaryObject bound,
DofContainer& dofs) const;
void getHigherOrderDofs(FiniteElemSpace* feSpace,
void getHigherOrderDofs(const FiniteElemSpace* feSpace,
BoundaryObject bound,
DofContainer& dofs) const;
......
......@@ -87,7 +87,7 @@ namespace AMDiS {
}
void Triangle::getNodeDofs(FiniteElemSpace* feSpace,
void Triangle::getNodeDofs(const FiniteElemSpace* feSpace,
BoundaryObject bound,
DofContainer& dofs) const
{
......@@ -111,10 +111,10 @@ namespace AMDiS {
child[1]->getSecondChild()->getNodeDofs(feSpace, nextBound, dofs);
dofs.push_back(child[1]->getFirstChild()->getDof(2));
nextBound.ithObj = 0;
child[1]->getFirstChild()->getNodeDofs(feSpace, nextBound, dofs);
child[1]->getFirstChild()->getNodeDofs(feSpace, nextBound, dofs);
} else {
nextBound.ithObj = 0;
child[1]->getFirstChild()->getNodeDofs(feSpace, nextBound, dofs);
child[1]->getFirstChild()->getNodeDofs(feSpace, nextBound, dofs);
dofs.push_back(child[1]->getFirstChild()->getDof(2));
nextBound.ithObj = 1;
child[1]->getSecondChild()->getNodeDofs(feSpace, nextBound, dofs);
......@@ -164,7 +164,7 @@ namespace AMDiS {
}
void Triangle::getHigherOrderDofs(FiniteElemSpace* feSpace,
void Triangle::getHigherOrderDofs(const FiniteElemSpace* feSpace,