Commit 531be06c authored by Thomas Witkowski's avatar Thomas Witkowski

Some large changes in code, seems to work....

parent fe752c0c
......@@ -98,6 +98,7 @@ namespace AMDiS {
class VertexVector;
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
class ElementObjectDatabase;
class FeSpaceDofMap;
class MeshLevelData;
#endif
......
......@@ -10,15 +10,45 @@
// See also license.opensource.txt in the distribution.
#include "DofComm.h"
#include "parallel/DofComm.h"
#include "parallel/InteriorBoundary.h"
#include "FiniteElemSpace.h"
namespace AMDiS {
using namespace std;
void DofComm::removeEmpty()
int DofComm::getNumberDofs(int level, const FiniteElemSpace *feSpace)
{
FUNCNAME("DofComm::removeEmpty()");
FUNCNAME("DofComm::getNumberDofs()");
TEST_EXIT_DBG(level < data.size())("Should not happen!\n");
DofContainerSet dofs;
for (DataIter rankIt = data[level].begin();
rankIt != data[level].end(); ++rankIt)
for (FeMapIter feIt = rankIt->second.begin();
feIt != rankIt->second.end(); ++feIt)
if (feIt->first == feSpace)
dofs.insert(feIt->second.begin(), feIt->second.end());
return static_cast<int>(dofs.size());
}
void DofComm::create(RankToBoundMap &boundary)
{
// === Fill data. ===
for (unsigned int i = 0; i < feSpaces.size(); i++)
for (int level = 0; level < nLevel; level++)
for (InteriorBoundary::iterator it(boundary, level); !it.end(); ++it)
it->rankObj.el->getAllDofs(feSpaces[i], it->rankObj,
data[level][it.getRank()][feSpaces[i]]);
// === Remove empty data containers. ===
for (unsigned int i = 0; i < data.size(); i++) {
DataIter dit = data[i].begin();
......@@ -42,25 +72,6 @@ namespace AMDiS {
}
int DofComm::getNumberDofs(int level, const FiniteElemSpace *feSpace)
{
FUNCNAME("DofComm::getNumberDofs()");
TEST_EXIT_DBG(level < data.size())("Should not happen!\n");
DofContainerSet dofs;
for (DataIter rankIt = data[level].begin();
rankIt != data[level].end(); ++rankIt)
for (FeMapIter feIt = rankIt->second.begin();
feIt != rankIt->second.end(); ++feIt)
if (feIt->first == feSpace)
dofs.insert(feIt->second.begin(), feIt->second.end());
return static_cast<int>(dofs.size());
}
bool DofComm::Iterator::setNextFeMap()
{
FUNCNAME("DofComm::Iterator::setNextFeMap()");
......
......@@ -24,6 +24,7 @@
#define AMDIS_DOF_COMM_H
#include <map>
#include "parallel/ParallelTypes.h"
#include "FiniteElemSpace.h"
#include "Global.h"
......@@ -52,12 +53,13 @@ namespace AMDiS {
return data[level][rank][feSpace];
}
void removeEmpty();
void init(int nLevels = 1)
void init(int n, vector<const FiniteElemSpace*> &fe)
{
nLevel = n;
feSpaces = fe;
data.clear();
data.resize(nLevels);
data.resize(nLevel);
}
DataType& getData(int level = 0)
......@@ -67,9 +69,15 @@ namespace AMDiS {
int getNumberDofs(int level, const FiniteElemSpace *feSpace);
void create(RankToBoundMap &boundary);
protected:
LevelDataType data;
int nLevel;
vector<const FiniteElemSpace*> feSpaces;
friend class Iterator;
public:
......
......@@ -15,6 +15,49 @@
namespace AMDiS {
void ElementObjectDatabase::create()
{
// === Fills macro element data structures. ===
TraverseStack stack;
ElInfo *elInfo = stack.traverseFirst(mesh, -1,
Mesh::CALL_LEAF_EL |
Mesh::FILL_NEIGH |
Mesh::FILL_BOUND);
while (elInfo) {
TEST_EXIT_DBG(elInfo->getLevel() == 0)("Should not happen!\n");
Element *el = elInfo->getElement();
macroElIndexMap.insert(make_pair(el->getIndex(), el));
macroElIndexTypeMap.insert(make_pair(el->getIndex(), elInfo->getType()));
// Add all sub object of the element to the variable elObjDb.
addElement(elInfo);
elInfo = stack.traverseNext(elInfo);
}
// Create periodic data, if there are periodic boundary conditions.
createPeriodicData();
// Create data about the reverse modes of neighbouring elements.
createReverseModeData();
}
void ElementObjectDatabase::createMacroElementInfo(vector<MacroElement*> &mel)
{
macroElIndexMap.clear();
macroElIndexTypeMap.clear();
for (vector<MacroElement*>::iterator it = mel.begin();
it != mel.end(); ++it) {
macroElIndexMap.insert(make_pair((*it)->getIndex(), (*it)->getElement()));
macroElIndexTypeMap.insert(make_pair((*it)->getIndex(), (*it)->getElType()));
}
}
void ElementObjectDatabase::addElement(ElInfo *elInfo)
{
FUNCNAME("ElementObjectDatabase::addElement()");
......@@ -128,7 +171,7 @@ namespace AMDiS {
}
void ElementObjectDatabase::createPeriodicData(const FiniteElemSpace *feSpace)
void ElementObjectDatabase::createPeriodicData()
{
FUNCNAME("ElementObjectDatabase::createPeriodicData()");
......@@ -223,9 +266,9 @@ namespace AMDiS {
("Should not happen!\n");
periodicVertices[make_pair(dof0, dof3)] =
provideConnectedPeriodicBoundary(feSpace->getAdmin(), type0, type1);
provideConnectedPeriodicBoundary(type0, type1);
periodicVertices[make_pair(dof3, dof0)] =
provideConnectedPeriodicBoundary(feSpace->getAdmin(), type0, type1);
provideConnectedPeriodicBoundary(type0, type1);
for (unsigned int j = i + 1; j < multPeriodicDof2.size(); j++)
if (multPeriodicDof2[j] == dof3)
......@@ -244,7 +287,7 @@ namespace AMDiS {
make_pair(multPeriodicDof3[j], multPeriodicDof3[i]);
if (periodicVertices.count(perDofs0) == 0) {
BoundaryType b = getNewBoundaryType(feSpace->getAdmin());
BoundaryType b = getNewBoundaryType();
periodicVertices[perDofs0] = b;
periodicVertices[perDofs1] = b;
}
......@@ -275,8 +318,7 @@ namespace AMDiS {
BoundaryType type0 = periodicEdges[make_pair(edge0, edge1)];
BoundaryType type1 = periodicEdges[make_pair(edge0, edge2)];
BoundaryType type2 =
provideConnectedPeriodicBoundary(feSpace->getAdmin(), type0, type1);
BoundaryType type2 = provideConnectedPeriodicBoundary(type0, type1);
periodicEdges[perEdge0] = type2;
periodicEdges[perEdge1] = type2;
......@@ -321,7 +363,7 @@ namespace AMDiS {
}
BoundaryType ElementObjectDatabase::getNewBoundaryType(DOFAdmin *admin)
BoundaryType ElementObjectDatabase::getNewBoundaryType()
{
FUNCNAME("ElementObjectDatabase::getNewBoundaryType()");
......@@ -334,15 +376,14 @@ namespace AMDiS {
newPeriodicBoundaryType--;
mesh->getPeriodicAssociations()[newPeriodicBoundaryType] =
new VertexVector(admin, "");
new VertexVector(feSpace->getAdmin(), "");
return newPeriodicBoundaryType;
}
BoundaryType
ElementObjectDatabase::provideConnectedPeriodicBoundary(DOFAdmin *admin,
BoundaryType b0,
ElementObjectDatabase::provideConnectedPeriodicBoundary(BoundaryType b0,
BoundaryType b1)
{
FUNCNAME("ElementObjectDatabase::provideConnectedPeriodicBoundary()");
......@@ -351,13 +392,13 @@ namespace AMDiS {
(b0 <= b1 ? make_pair(b0, b1) : make_pair(b1, b0));
if (bConnMap.count(bConn) == 0) {
BoundaryType newPeriodicBoundaryType = getNewBoundaryType(admin);
BoundaryType newPeriodicBoundaryType = getNewBoundaryType();
VertexVector &vecB0 = mesh->getPeriodicAssociations(b0);
VertexVector &vecB1 = mesh->getPeriodicAssociations(b1);
VertexVector &vecC = mesh->getPeriodicAssociations(newPeriodicBoundaryType);
DOFIteratorBase it(const_cast<DOFAdmin*>(admin), USED_DOFS);
DOFIteratorBase it(const_cast<DOFAdmin*>(feSpace->getAdmin()), USED_DOFS);
for (it.reset(); !it.end(); ++it) {
if (!it.isDofFree()) {
......@@ -478,9 +519,7 @@ namespace AMDiS {
}
void ElementObjectDatabase::createReverseModeData(const FiniteElemSpace* feSpace,
map<int, Element*> &elIndexMap,
map<int, int> &elIndexTypeMap)
void ElementObjectDatabase::createReverseModeData()
{
FUNCNAME("ElementObjectDatabase::createReverseModeData()");
......@@ -498,13 +537,13 @@ namespace AMDiS {
vector<ElementObjectData>& els = edgeIt->second;
for (unsigned int i = 0; i < els.size(); i++) {
BoundaryObject obj0(elIndexMap[els[i].elIndex],
elIndexTypeMap[els[i].elIndex],
BoundaryObject obj0(macroElIndexMap[els[i].elIndex],
macroElIndexTypeMap[els[i].elIndex],
EDGE, els[i].ithObject);
for (unsigned int j = i + 1; j < els.size(); j++) {
BoundaryObject obj1(elIndexMap[els[j].elIndex],
elIndexTypeMap[els[j].elIndex],
BoundaryObject obj1(macroElIndexMap[els[j].elIndex],
macroElIndexTypeMap[els[j].elIndex],
EDGE, els[j].ithObject);
bool reverseMode =
......@@ -522,13 +561,13 @@ namespace AMDiS {
vector<ElementObjectData>& els = faceIt->second;
for (unsigned int i = 0; i < els.size(); i++) {
BoundaryObject obj0(elIndexMap[els[i].elIndex],
elIndexTypeMap[els[i].elIndex],
BoundaryObject obj0(macroElIndexMap[els[i].elIndex],
macroElIndexTypeMap[els[i].elIndex],
FACE, els[i].ithObject);
for (unsigned int j = i + 1; j < els.size(); j++) {
BoundaryObject obj1(elIndexMap[els[j].elIndex],
elIndexTypeMap[els[j].elIndex],
BoundaryObject obj1(macroElIndexMap[els[j].elIndex],
macroElIndexTypeMap[els[j].elIndex],
FACE, els[j].ithObject);
bool reverseMode =
......@@ -549,13 +588,13 @@ namespace AMDiS {
vector<ElementObjectData> &edges1 = edgeElements[edgeIt->first.second];
for (unsigned int i = 0; i < edges0.size(); i++) {
BoundaryObject obj0(elIndexMap[edges0[i].elIndex],
elIndexTypeMap[edges0[i].elIndex],
BoundaryObject obj0(macroElIndexMap[edges0[i].elIndex],
macroElIndexTypeMap[edges0[i].elIndex],
EDGE, edges0[i].ithObject);
for (unsigned int j = 0; j < edges1.size(); j++) {
BoundaryObject obj1(elIndexMap[edges1[j].elIndex],
elIndexTypeMap[edges1[j].elIndex],
BoundaryObject obj1(macroElIndexMap[edges1[j].elIndex],
macroElIndexTypeMap[edges1[j].elIndex],
EDGE, edges1[j].ithObject);
bool reverseMode =
......@@ -575,11 +614,11 @@ namespace AMDiS {
TEST_EXIT_DBG(faces0.size() == faces1.size() == 1)("Should not happen!\n");
BoundaryObject obj0(elIndexMap[faces0[0].elIndex],
elIndexTypeMap[faces0[0].elIndex],
BoundaryObject obj0(macroElIndexMap[faces0[0].elIndex],
macroElIndexTypeMap[faces0[0].elIndex],
FACE, faces0[0].ithObject);
BoundaryObject obj1(elIndexMap[faces1[0].elIndex],
elIndexTypeMap[faces1[0].elIndex],
BoundaryObject obj1(macroElIndexMap[faces1[0].elIndex],
macroElIndexTypeMap[faces1[0].elIndex],
FACE, faces1[0].ithObject);
bool reverseMode =
......
......@@ -103,39 +103,25 @@ namespace AMDiS {
class ElementObjectDatabase {
public:
ElementObjectDatabase()
: mesh(NULL),
: feSpace(NULL),
mesh(NULL),
iterGeoPos(CENTER)
{}
/// Set the mesh that should be used for the database.
void setMesh(Mesh *m)
void setFeSpace(const FiniteElemSpace *fe)
{
mesh = m;
feSpace = fe;
mesh = feSpace->getMesh();
}
Mesh* getMesh()
{
return mesh;
}
void create();
/** \brief
* Adds an element to the object database. If the element is part of a
* periodic boundary, all information about subobjects of the element on
* this boundary are collected.
*
* \param[in] elInfo ElInfo object of the element.
*/
void addElement(ElInfo *elInfo);
/** \brief
* 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
* search for indirectly connected vertices in periodic boundaries. This is
* only the case, if there are more than one boundary conditions. Then, e.g.,
* in 2D, all edges of a square are iterectly connected. In 3D, if the macro
* mesh is a box, all eight vertex nodes and always four of the 12 edges are
* indirectly connected.
*/
void createPeriodicData(const FiniteElemSpace *feSpace);
void createMacroElementInfo(vector<MacroElement*> &mel);
/** \brief
* Create for a filled object database the membership information for all
......@@ -149,22 +135,6 @@ namespace AMDiS {
MeshLevelData& levelData);
/** \brief
* Creates on all boundaries the reverse mode flag.
*
* \param[in] feSpace An arbitrary FE space defined on the mesh.
* Is used to get the orientation of the DOFs on
* elements.
* \param[in] elIndexMap Maps an element index to the pointer to the
* element.
* \param[in] elIndexTypeMap Maps an element index to its type id (only
* relevant in 3D).
*/
void createReverseModeData(const FiniteElemSpace* feSpace,
map<int, Element*> &elIndexMap,
map<int, int> &elIndexTypeMap);
/** \brief
* Iterates over all elements for one geometrical index, i.e., over all
* vertices, edges or faces in the mesh. The function returns true, if the
......@@ -480,6 +450,16 @@ namespace AMDiS {
return (t >= smallestPeriodicBcType);
}
inline Element* getElementPtr(int index)
{
return macroElIndexMap[index];
}
inline int getElementType(int index)
{
return macroElIndexTypeMap[index];
}
/// Write the element database to disk.
void serialize(ostream &out);
......@@ -487,6 +467,15 @@ namespace AMDiS {
void deserialize(istream &in);
protected:
/** \brief
* Adds an element to the object database. If the element is part of a
* periodic boundary, all information about subobjects of the element on
* this boundary are collected.
*
* \param[in] elInfo ElInfo object of the element.
*/
void addElement(ElInfo *elInfo);
/// Adds the i-th DOF vertex of an element to the object database.
void addVertex(Element *el, int ith)
{
......@@ -520,10 +509,24 @@ namespace AMDiS {
faceLocalMap[elObj] = face;
}
BoundaryType getNewBoundaryType(DOFAdmin *admin);
/** \brief
* 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
* search for indirectly connected vertices in periodic boundaries. This is
* only the case, if there are more than one boundary conditions. Then, e.g.,
* in 2D, all edges of a square are iterectly connected. In 3D, if the macro
* mesh is a box, all eight vertex nodes and always four of the 12 edges are
* indirectly connected.
*/
void createPeriodicData();
BoundaryType provideConnectedPeriodicBoundary(DOFAdmin *admin,
BoundaryType b0,
/// Creates on all boundaries the reverse mode flag.
void createReverseModeData();
BoundaryType getNewBoundaryType();
BoundaryType provideConnectedPeriodicBoundary(BoundaryType b0,
BoundaryType b1);
/// Some auxiliary function to write the element object database to disk.
......@@ -539,10 +542,12 @@ namespace AMDiS {
void deserialize(istream &in, map<int, ElementObjectData>& data);
private:
const FiniteElemSpace* feSpace;
/// The mesh that is used to store all its element information in
/// the database.
Mesh *mesh;
/// Maps to each vertex DOF all element objects that represent this vertex.
map<DegreeOfFreedom, vector<ElementObjectData> > vertexElements;
......@@ -630,6 +635,12 @@ namespace AMDiS {
map<pair<ElementObjectData, ElementObjectData>, bool> edgeReverseMode;
map<pair<ElementObjectData, ElementObjectData>, bool> faceReverseMode;
/// Maps to each macro element index a pointer to the corresponding element.
map<int, Element*> macroElIndexMap;
/// Maps to each macro element index the type of this element.
map<int, int> macroElIndexTypeMap;
};
}
......
......@@ -11,6 +11,7 @@
#include "parallel/InteriorBoundary.h"
#include "parallel/ElementObjectDatabase.h"
#include "FiniteElemSpace.h"
#include "BasisFunction.h"
#include "Serializer.h"
......@@ -18,50 +19,352 @@
namespace AMDiS {
AtomicBoundary& InteriorBoundary::getNewAtomic(int rank)
AtomicBoundary& InteriorBoundary::getNewAtomicOwn(int rank)
{
int size = boundary[rank].size();
boundary[rank].resize(size + 1);
return boundary[rank][size];
int size = own[rank].size();
own[rank].resize(size + 1);
return own[rank][size];
}
AtomicBoundary& InteriorBoundary::getNewAtomicOther(int rank)
{
int size = other[rank].size();
other[rank].resize(size + 1);
return other[rank][size];
}
bool InteriorBoundary::operator==(const InteriorBoundary& other) const
AtomicBoundary& InteriorBoundary::getNewAtomicPer(int rank)
{
InteriorBoundary& other2 = const_cast<InteriorBoundary&>(other);
int size = periodic[rank].size();
periodic[rank].resize(size + 1);
return periodic[rank][size];
}
void InteriorBoundary::create(MPI::Intracomm &mpiComm,
ElementObjectDatabase &elObjDb)
{
FUNCNAME("InteriorBoundary::clear()");
own.clear();
other.clear();
periodic.clear();
Mesh *mesh = elObjDb.getMesh();
TEST_EXIT_DBG(mesh)("Should not happen!\n");
int mpiRank = mpiComm.Get_rank();
if (boundary.size() != other2.boundary.size())
return false;
// === Create interior boundary data structure. ===
for (int geoPos = 0; geoPos < mesh->getDim(); geoPos++) {
GeoIndex geoIndex = INDEX_OF_DIM(geoPos, mesh->getDim());
for (unsigned int level = 0; level < boundary.size(); level++) {
for (RankToBoundMap::const_iterator it = boundary.begin();
it != boundary.end(); ++it) {
if (other2.boundary.count(it->first) == 0)
return false;
while (elObjDb.iterate(geoIndex)) {
map<int, ElementObjectData>& objData = elObjDb.getIterateData();
if (!(objData.count(mpiRank) && objData.size() > 1))
continue;
int owner = elObjDb.getIterateOwner();
ElementObjectData& rankBoundEl = objData[mpiRank];
AtomicBoundary bound;
bound.maxLevel = elObjDb.getIterateMaxLevel();
bound.rankObj.el = elObjDb.getElementPtr(rankBoundEl.elIndex);
bound.rankObj.elIndex = rankBoundEl.elIndex;
bound.rankObj.elType = elObjDb.getElementType(rankBoundEl.elIndex);
bound.rankObj.subObj = geoIndex;
bound.rankObj.ithObj = rankBoundEl.ithObject;
if (geoIndex == FACE) {
for (int edgeNo = 0; edgeNo < 3; edgeNo++) {
int edgeOfFace =
bound.rankObj.el->getEdgeOfFace(bound.rankObj.ithObj, edgeNo);
bound.rankObj.excludedSubstructures.push_back(make_pair(EDGE, edgeOfFace));
}
}
if (other2.boundary[it->first].size() != it->second.size())
return false;
for (unsigned int i = 0; i < it->second.size(); i++) {
std::vector<AtomicBoundary>::iterator bIt =
find(other2.boundary[it->first].begin(),
other2.boundary[it->first].end(), it->second[i]);
if (owner == mpiRank) {
for (map<int, ElementObjectData>::iterator it2 = objData.begin();
it2 != objData.end(); ++it2) {
if (it2->first == mpiRank)
continue;
bound.neighObj.el = elObjDb.getElementPtr(it2->second.elIndex);
bound.neighObj.elIndex = it2->second.elIndex;
bound.neighObj.elType = elObjDb.getElementType(it2->second.elIndex);
bound.neighObj.subObj = geoIndex;
bound.neighObj.ithObj = it2->second.ithObject;
bound.type = INTERIOR;
AtomicBoundary& b = getNewAtomicOwn(it2->first);
b = bound;
if (geoIndex == EDGE)
b.neighObj.reverseMode =
elObjDb.getEdgeReverseMode(rankBoundEl, it2->second);
if (geoIndex == FACE)
b.neighObj.reverseMode =
elObjDb.getFaceReverseMode(rankBoundEl, it2->second);
}
} else {
TEST_EXIT_DBG(objData.count(owner) == 1)
("Should not happen!\n");
ElementObjectData& ownerBoundEl = objData[owner];
if (bIt == other2.boundary[it->first].end())
return false;
}
bound.neighObj.el = elObjDb.getElementPtr(ownerBoundEl.elIndex);
bound.neighObj.elIndex = ownerBoundEl.elIndex;
bound.neighObj.elType = -1;
bound.neighObj.subOb