Commit 4f99cef1 authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

Solved some bugs in parallel code.

parent 9420db7e
......@@ -11,6 +11,7 @@
#include "DOFVector.h"
#include "DOFIterator.h"
#include "Serializer.h"
#include "VertexVector.h"
namespace AMDiS {
......@@ -24,6 +25,7 @@ namespace AMDiS {
init();
}
DOFAdmin::DOFAdmin(Mesh* m, std::string aName)
: name(aName),
mesh(m),
......@@ -33,9 +35,11 @@ namespace AMDiS {
init();
}
DOFAdmin::~DOFAdmin()
{}
void DOFAdmin::init()
{
firstHole = 0;
......@@ -46,6 +50,7 @@ namespace AMDiS {
dofFree.clear();
}
DOFAdmin& DOFAdmin::operator=(const DOFAdmin& src)
{
if (this != &src) {
......@@ -67,10 +72,7 @@ namespace AMDiS {
return *this;
}
/****************************************************************************/
/* use a bit vector to indicate used/unused dofs */
/* storage needed: one bit per dof */
/****************************************************************************/
bool DOFAdmin::operator==(const DOFAdmin& ad) const
{
......@@ -90,9 +92,10 @@ namespace AMDiS {
ERROR_EXIT("TODO\n");
}
void DOFAdmin::freeDOFIndex(int dof)
void DOFAdmin::freeDofIndex(int dof)
{
FUNCNAME("DOFAdmin::freeDOFIndex()");
FUNCNAME("DOFAdmin::freeDofIndex()");
TEST_EXIT_DBG(usedCount > 0)("no dofs in use\n");
TEST_EXIT_DBG((dof >= 0) && (dof < size))("invalid dof index %d\n",dof);
......@@ -107,7 +110,7 @@ namespace AMDiS {
std::list<DOFContainer*>::iterator dcend = dofContainerList.end();
for (dc = dofContainerList.begin(); dc != dcend; ++dc)
(*dc)->freeDOFIndex(dof);
(*dc)->freeDofIndex(dof);
dofFree[dof] = true;
......@@ -118,6 +121,7 @@ namespace AMDiS {
holeCount++;
}
int DOFAdmin::getDOFIndex()
{
FUNCNAME("DOFAdmin::getDOFIndex()");
......@@ -157,6 +161,7 @@ namespace AMDiS {
return dof;
}
void DOFAdmin::enlargeDOFLists(int minsize)
{
FUNCNAME("DOFAdmin::enlargeDOFLists()");
......@@ -184,6 +189,7 @@ namespace AMDiS {
(*di)->resize(newval);
}
void DOFAdmin::addDOFIndexed(DOFIndexedBase* dofIndexed)
{
FUNCNAME("DOFAdmin::addDOFIndexed()");
......@@ -201,6 +207,7 @@ namespace AMDiS {
}
}
void DOFAdmin::removeDOFIndexed(DOFIndexedBase* dofIndexed)
{
FUNCNAME("DOFAdmin::removeDOFIndexed()");
......@@ -221,16 +228,20 @@ namespace AMDiS {
}
}
TEST_EXIT(removed)("DOFIndexed not in list\n");
TEST_EXIT(removed)("DOFIndexed not in list!\n");
}
void DOFAdmin::addDOFContainer(DOFContainer* cont)
{
FUNCNAME("DOFAdmin::addDOFContainer()");
TEST_EXIT_DBG(cont)("no container\n");
dofContainerList.push_back(cont);
}
void DOFAdmin::removeDOFContainer(DOFContainer* cont)
{
FUNCNAME("DOFAdmin::removeDOFContainer()");
......@@ -243,7 +254,8 @@ namespace AMDiS {
return;
}
}
ERROR("container not in list\n");
ERROR("Container not in list!\n");
}
......@@ -262,7 +274,7 @@ namespace AMDiS {
// mark used dofs
DOFIteratorBase it(this, USED_DOFS);
for (it.reset(); !it.end(); ++it)
new_dof[it.getDOFIndex()] = 1;
new_dof[it.getDOFIndex()] = 1;
int n = 0, last = 0;
for (int i = 0; i < size; i++) { /* create a MONOTONE compress */
......@@ -288,7 +300,7 @@ namespace AMDiS {
// get index of first changed dof
int first = last;
for (int i = 0; i<size; i++) {
for (int i = 0; i < size; i++) {
if (new_dof[i] < i && new_dof[i] >= 0) {
first = i;
break;
......@@ -303,8 +315,8 @@ namespace AMDiS {
std::list<DOFContainer*>::iterator dc;
std::list<DOFContainer*>::iterator endc = dofContainerList.end();
for (dc = dofContainerList.begin(); dc != endc; dc++)
(*dc)->compressDOFContainer(n, new_dof);
for (dc = dofContainerList.begin(); dc != endc; ++dc)
(*dc)->compressDofContainer(n, new_dof);
}
......@@ -317,6 +329,7 @@ namespace AMDiS {
nrDOF[i] = v;
}
void DOFAdmin::setNumberOfPreDOFs(int i, int v)
{
FUNCNAME("DOFAdmin::setNumberOfPreDOFs()");
......@@ -326,11 +339,13 @@ namespace AMDiS {
nr0DOF[i] = v;
}
int DOFAdmin::calcMemoryUsage()
{
return sizeof(DOFAdmin);
}
void DOFAdmin::serialize(std::ostream &out)
{
// write name
......@@ -352,7 +367,8 @@ namespace AMDiS {
nrDOF.serialize(out);
nr0DOF.serialize(out);
}
}
void DOFAdmin::deserialize(std::istream &in)
{
......
......@@ -95,6 +95,11 @@ namespace AMDiS {
/// Removes the given DOFContainer object from DOFAdmin.
void removeDOFContainer(DOFContainer* dofContainer);
std::list<DOFContainer*>& getDCL()
{
return dofContainerList;
}
/** \brief
* Removes all holes of unused DOF indices by compressing the used range of
* indices (it does not resize the vectors). While the global index of a DOF
......@@ -257,7 +262,7 @@ namespace AMDiS {
int getDOFIndex();
/// Frees index dof. Used by Mesh::getDOF()
void freeDOFIndex(int dof);
void freeDofIndex(int dof);
///
void serialize(std::ostream &out);
......
......@@ -44,22 +44,25 @@ namespace AMDiS {
*/
virtual DegreeOfFreedom& operator[](int i) = 0;
virtual void freeDOFIndex(DegreeOfFreedom dof) {}
virtual void freeDofIndex(DegreeOfFreedom dof) {}
/** \brief
* Used by DOFAdmin to actualize the DOF indices in this container after
* DOF compression.
*/
virtual void compressDOFContainer(int size, std::vector<DegreeOfFreedom> &newDOF)
virtual void compressDofContainer(int size, std::vector<DegreeOfFreedom> &newDOF)
{
FUNCNAME("DOFContainer::compressDofContainer()");
for (int i = 0; i < size; i++) {
int j = newDOF[operator[](i)];
if (j >= 0)
operator[](i) = j;
else
ERROR_EXIT("invalid dof in dof container\n");
ERROR_EXIT("Invalid DOF %d in DOF container! (%d %d)\n", j, i, operator[](i));
}
}
};
}
......
......@@ -447,7 +447,7 @@ namespace AMDiS {
{
FUNCNAME("DOFVector<T>::operator[]");
TEST_EXIT_DBG(i >= 0 && i < static_cast<int>(vec.size()))
("Illegal vector index %d.\n",i);
("Illegal vector index %d.\n", i);
return vec[i];
}
......@@ -457,7 +457,7 @@ namespace AMDiS {
FUNCNAME("DOFVector<T>::operator[]");
TEST_EXIT_DBG(i >= 0 && i < static_cast<int>(vec.size()))
("Illegal vector index %d.\n",i);
("Illegal vector index %d.\n", i);
return vec[i];
}
......@@ -639,7 +639,7 @@ namespace AMDiS {
/** \brief
* Implements DOFContainer::operator[]() by calling
* DOFVEctor<DegreeOfFreedom>::operator[]()
* DOFVector<DegreeOfFreedom>::operator[]()
*/
DegreeOfFreedom& operator[](DegreeOfFreedom i)
{
......
......@@ -191,7 +191,6 @@ namespace AMDiS {
}
return NULL;
}
......@@ -486,6 +485,32 @@ namespace AMDiS {
}
void printElementRefinementSequenz(Mesh *mesh, Element *el)
{
FUNCNAME("printElementRefinementSequenz()");
int elIndex = el->getIndex();
std::vector<int> refSeq;
Element *parent = getParentElement(mesh, el);
while (parent) {
if (parent->getChild(0) == el)
refSeq.push_back(0);
else
refSeq.push_back(1);
el = parent;
parent = getParentElement(mesh, el);
}
std::stringstream oss;
for (int i = refSeq.size() - 1; i >= 0; i--)
oss << refSeq[i];
MSG("Ref-Seq for element %d: %s\n", elIndex, oss.str().c_str());
}
int getLocalNeighbourIndex(Mesh *mesh, int elIndex, int neighIndex)
{
FUNCNAME("debug::getLocalNeighbourIndex()");
......
......@@ -137,6 +137,8 @@ namespace AMDiS {
void printElementHierarchie(Mesh *mesh, int elIndex);
void printElementRefinementSequenz(Mesh *mesh, Element *el);
int getLocalNeighbourIndex(Mesh *mesh, int elIndex, int neighIndex);
/** \brief
......
......@@ -590,7 +590,7 @@ namespace AMDiS {
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]);
localAdmin->freeDofIndex(dof[n0 + j]);
}
delete [] dof;
......
......@@ -51,7 +51,7 @@ namespace AMDiS {
// === Get the refinement patch. ===
getRefinePatch(&elInfo, edge, 0, refineList, &n_neigh);
// === Check for periodic boundary ===
......
......@@ -740,6 +740,7 @@ namespace AMDiS {
/* first, goto to leaf level, if necessary... */
if (!(el->isLeaf()) && neighbour < 2) {
if (stack_used >= static_cast<int>(elinfo_stack.size()) - 1)
enlargeTraverseStack();
......@@ -780,6 +781,7 @@ namespace AMDiS {
while (stack_used > 1) {
stack_used--;
int elIsIthChild = info_stack[stack_used];
if (traverse_fill_flag.isSet(Mesh::CALL_REVERSE_MODE) && elIsIthChild != 0)
elIsIthChild = (elIsIthChild == 1 ? 2 : 1);
......@@ -790,8 +792,8 @@ namespace AMDiS {
nb = coarse_nb[elIsIthChild][nb];
if (nb == -1)
break;
break;
TEST_EXIT_DBG(nb >= 0)("invalid coarse_nb %d\n",nb);
}
......@@ -818,7 +820,7 @@ namespace AMDiS {
if (traverse_mel == NULL)
return NULL;
nb = i;
stack_used = 1;
elinfo_stack[stack_used]->fillMacroInfo(const_cast<MacroElement*>(traverse_mel));
info_stack[stack_used] = 0;
......@@ -834,15 +836,27 @@ namespace AMDiS {
if (stack_used >= stack_size - 1)
enlargeTraverseStack();
int i = 2 - info_stack[stack_used];
info_stack[stack_used] = i + 1;
int fillIthChild = i;
if (traverse_fill_flag.isSet(Mesh::CALL_REVERSE_MODE))
TEST_EXIT_DBG(info_stack[stack_used] == 1 || info_stack[stack_used] == 2)
("Should not happen!\n");
int fillIthChild = -1;
if (info_stack[stack_used] == 1) {
info_stack[stack_used] = 2;
fillIthChild = 1;
nb = 0;
} else {
info_stack[stack_used] = 1;
fillIthChild = 0;
nb = 1;
}
if (traverse_fill_flag.isSet(Mesh::CALL_REVERSE_MODE)) {
fillIthChild = 1 - fillIthChild;
nb = 1 - nb;
}
elinfo_stack[stack_used + 1]->fillElInfo(fillIthChild, elinfo_stack[stack_used]);
stack_used++;
nb = 1 - i;
}
/****************************************************************************/
......@@ -861,12 +875,16 @@ namespace AMDiS {
if (stack_used >= stack_size - 1)
enlargeTraverseStack();
if (traverse_fill_flag.isSet(Mesh::CALL_REVERSE_MODE)) {
int t = 2 - nb;
info_stack[stack_used] = (t == 2 ? 1 : 2);
} else {
info_stack[stack_used] = 2 - nb;
}
int fillIthChild = 1 - nb;
if (traverse_fill_flag.isSet(Mesh::CALL_REVERSE_MODE))
fillIthChild = 1 - fillIthChild;
elinfo_stack[stack_used + 1]->fillElInfo(fillIthChild, elinfo_stack[stack_used]);
info_stack[stack_used] = 2 - nb;
stack_used++;
nb = 2;
}
......@@ -882,12 +900,16 @@ namespace AMDiS {
info_stack[stack_used] = i + 1;
int fillIthChild = i;
if (traverse_fill_flag.isSet(Mesh::CALL_REVERSE_MODE))
nb = i;
if (traverse_fill_flag.isSet(Mesh::CALL_REVERSE_MODE)) {
fillIthChild = 1 - fillIthChild;
nb = 1 - i;
}
elinfo_stack[stack_used + 1]->fillElInfo(fillIthChild, elinfo_stack[stack_used]);
stack_used++;
nb = i;
elinfo = elinfo_stack[stack_used];
el = dynamic_cast<Triangle*>(const_cast<Element*>(elinfo->getElement()));
......
......@@ -36,9 +36,9 @@ namespace AMDiS {
void set(DegreeOfFreedom val);
void compressDOFContainer(int size, std::vector<DegreeOfFreedom> &newDOF)
void compressDofContainer(int size, std::vector<DegreeOfFreedom> &newDOF)
{
DOFContainer::compressDOFContainer(size, newDOF);
DOFContainer::compressDofContainer(size, newDOF);
int totalSize = getAdmin()->getSize();
for (int i = size; i < totalSize; i++)
vec[i] = i;
......
......@@ -425,6 +425,8 @@ namespace AMDiS {
void MeshDistributor::removePeriodicBoundaryConditions()
{
FUNCNAME("MeshDistributor::removePeriodicBoundaryConditions()");
// Remove periodic boundaries in boundary manager on matrices and vectors.
for (unsigned int i = 0; i < probStat.size(); i++) {
int nComponents = probStat[i]->getNumComponents();
......@@ -451,6 +453,12 @@ namespace AMDiS {
elInfo->getElement()->deleteElementData(PERIODIC);
elInfo = stack.traverseNext(elInfo);
}
// Remove periodic vertex associations
for (std::map<BoundaryType, VertexVector*>::iterator it = mesh->getPeriodicAssociations().begin();
it != mesh->getPeriodicAssociations().end(); ++it)
const_cast<DOFAdmin&>(mesh->getDOFAdmin(0)).removeDOFContainer(dynamic_cast<DOFContainer*>(it->second));
mesh->getPeriodicAssociations().clear();
}
......@@ -572,7 +580,7 @@ namespace AMDiS {
bool meshFitTogether = true;
for (RankToBoundMap::iterator it = allBound.begin(); it != allBound.end(); ++it) {
MeshCodeVec &recvCodes = stdMpi.getRecvData()[it->first];
int i = 0;
......@@ -590,7 +598,7 @@ namespace AMDiS {
boundIt->rankObj.subObj,
boundIt->rankObj.ithObj,
boundIt->rankObj.elType,
boundIt->neighObj.reverseMode);
boundIt->rankObj.reverseMode);
if (b)
meshFitTogether = false;
......@@ -640,7 +648,8 @@ namespace AMDiS {
TEST_EXIT_DBG(elInfo->getElement() == el)("This should not happen!\n");
meshChanged = fitElementToMeshCode2(code, stack, subObj, ithObj, elType, reverseMode);
meshChanged =
fitElementToMeshCode2(code, stack, subObj, ithObj, elType, reverseMode);
return meshChanged;
}
......@@ -999,9 +1008,11 @@ namespace AMDiS {
std::map<int, Element*> elIndexMap;
std::map<int, int> elIndexTypeMap;
std::map<std::pair<DegreeOfFreedom, DegreeOfFreedom>, BoundaryType> periodicDofs;
std::map<DegreeOfFreedom, std::set<BoundaryType> > periodicDofAssoc;
std::map<std::pair<DofEdge, DofEdge>, BoundaryType> periodicEdges;
std::map<std::pair<DofFace, DofFace>, BoundaryType> periodicFaces;
std::map<std::pair<DegreeOfFreedom, DegreeOfFreedom>, BoundaryType> periodicDofs;
// === PHASE 1 ===
......@@ -1052,6 +1063,7 @@ namespace AMDiS {
TEST_EXIT_DBG(neighDof > -1)("Should not happen!\n");
periodicDofs[std::make_pair(dof, neighDof)] = boundaryType;
periodicDofAssoc[dof].insert(boundaryType);
}
}
}
......@@ -1076,6 +1088,41 @@ namespace AMDiS {
elObjects.createRankData();
// === Search for interectly connected vertices in periodic boundaries. ===
if (periodicDofs.size() > 0) {
TEST_EXIT(mesh->getDim() == 2)("Parallel boundaries not yet supported in 3D!\n");
std::vector<DegreeOfFreedom> twoPeriodicDof;
for (std::map<DegreeOfFreedom, std::set<BoundaryType> >::iterator it = periodicDofAssoc.begin();
it != periodicDofAssoc.end(); ++it) {
TEST_EXIT_DBG(it->second.size() < 3)("Should not happen!\n");
if (it->second.size() == 2)
twoPeriodicDof.push_back(it->first);
}
TEST_EXIT_DBG(twoPeriodicDof.size() == 0 || twoPeriodicDof.size() == 4)
("Should not happen (%d)!\n", twoPeriodicDof.size());
if (twoPeriodicDof.size() == 4) {
for (int i = 0; i < 4; i++) {
for (int j = 0; j < 4; j++) {
if (i == j)
continue;
std::pair<DegreeOfFreedom, DegreeOfFreedom> perDofs =
std::make_pair(twoPeriodicDof[i], twoPeriodicDof[j]);
if (periodicDofs.count(perDofs) == 0)
periodicDofs[perDofs] = 3;
}
}
}
}
// === PHASE 3 ===
for (int geoPos = 0; geoPos < mesh->getDim(); geoPos++) {
......@@ -1813,75 +1860,6 @@ namespace AMDiS {
}
}
if (dofFromRank.size() > 0) {
TEST_EXIT_DBG(mesh->getDim() == 2)
("Periodic boundary corner problem must be generalized to 3d!\n");
}
MPI::Request request[min(static_cast<int>(periodicBoundary.boundary.size() * 2), 4)];
int requestCounter = 0;
std::vector<int*> sendBuffers, recvBuffers;
for (std::map<DegreeOfFreedom, std::set<int> >::iterator it = dofFromRank.begin();
it != dofFromRank.end(); ++it) {
if (it->second.size() == 2) {
TEST_EXIT_DBG(periodicDofAssociations[it->first].size() == 2)
("DOF %d has only %d periodic associations!\n",
it->first, periodicDofAssociations[it->first].size());
int type0 = *(periodicDofAssociations[it->first].begin());
int type1 = *(++(periodicDofAssociations[it->first].begin()));
int *sendbuf = new int[2];
sendbuf[0] = periodicDof[type0][it->first];
sendbuf[1] = periodicDof[type1][it->first];
request[requestCounter++] =
mpiComm.Isend(sendbuf, 2, MPI_INT, *(it->second.begin()), 0);
request[requestCounter++] =
mpiComm.Isend(sendbuf, 2, MPI_INT, *(++(it->second.begin())), 0);
sendBuffers.push_back(sendbuf);
int *recvbuf1 = new int[2];
int *recvbuf2 = new int[2];
request[requestCounter++] =
mpiComm.Irecv(recvbuf1, 2, MPI_INT, *(it->second.begin()), 0);
request[requestCounter++] =
mpiComm.Irecv(recvbuf2, 2, MPI_INT, *(++(it->second.begin())), 0);
recvBuffers.push_back(recvbuf1);
recvBuffers.push_back(recvbuf2);
}
}
MPI::Request::Waitall(requestCounter, request);
int i = 0;