Commit c520ca87 authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

Work on pdd and some other small code refactorings.

parent 073ad4f8
......@@ -44,7 +44,7 @@ available_tags=" CXX F77"
# ### BEGIN LIBTOOL CONFIG
# Libtool was configured on host p1d062:
# Libtool was configured on host p2q004:
# Shell to use when invoking shell scripts.
SHELL="/bin/sh"
......@@ -6760,7 +6760,7 @@ build_old_libs=`case $build_libtool_libs in yes) $echo no;; *) $echo yes;; esac`
# End:
# ### BEGIN LIBTOOL TAG CONFIG: CXX
# Libtool was configured on host p1d062:
# Libtool was configured on host p2q004:
# Shell to use when invoking shell scripts.
SHELL="/bin/sh"
......@@ -7065,7 +7065,7 @@ include_expsyms=""
# ### BEGIN LIBTOOL TAG CONFIG: F77
# Libtool was configured on host p1d062:
# Libtool was configured on host p2q004:
# Shell to use when invoking shell scripts.
SHELL="/bin/sh"
......
......@@ -84,7 +84,6 @@
#include "SystemVector.h"
#include "TecPlotWriter.h"
#include "Tetrahedron.h"
#include "TimedObject.h"
#include "Traverse.h"
#include "Triangle.h"
#include "ValueWriter.h"
......
......@@ -213,7 +213,7 @@ namespace AMDiS {
/// Sets \ref firstHole
inline void setFirstHole(int i)
{
TEST_EXIT_DBG(!dofFree[i])("There is no hole!\n");
TEST_EXIT_DBG(dofFree[i])("There is no hole!\n");
firstHole = i;
}
......
......@@ -216,7 +216,8 @@ namespace AMDiS {
applyDBCs.insert(static_cast<int>(row));
#endif
}
} else
} else {
for (int j = 0; j < nCol; j++) { // for all columns
DegreeOfFreedom col = colIndices[j];
double entry = elMat[i][j];
......@@ -224,8 +225,9 @@ namespace AMDiS {
if (add)
ins[row][col] += sign * entry;
else
ins[row][col] = sign * entry;
ins[row][col] = sign * entry;
}
}
}
}
......@@ -250,18 +252,6 @@ namespace AMDiS {
(*it)->getElementMatrix(elInfo, elementMatrix, *factorIt ? **factorIt : 1.0);
addElementMatrix(factor, elementMatrix, bound, elInfo, NULL);
// if (MPI::COMM_WORLD.Get_rank() == 0 && elInfo->getElement()->getIndex() == 53) {
// std::cout << elementMatrix << std::endl;
// rowFESpace->getBasisFcts()->getLocalIndicesVec(elInfo->getElement(),
// rowFESpace->getAdmin(),
// &rowIndices);
// rowIndices.print();
// print();
// }
}
void DOFMatrix::assemble(double factor, ElInfo *elInfo, const BoundaryType *bound,
......
......@@ -52,13 +52,9 @@ namespace AMDiS {
// Set of all DOFs of the rank.
std::vector<const DegreeOfFreedom*> rankDOFs;
// Set of all interior boundary DOFs in ranks partition which are owned by
// another rank.
std::map<const DegreeOfFreedom*, int> boundaryDOFs;
// Number of DOFs in ranks partition that are owned by the rank.
int nRankDOFs = 0;
// Number of DOFs in ranks partition that are at an interior boundary and are
// owned by other ranks.
// Number of all DOFs in the macro mesh.
int nOverallDOFs = 0;
createLocalGlobalNumbering(rankDOFs, boundaryDOFs, nRankDOFs, nOverallDOFs);
......@@ -74,7 +70,7 @@ namespace AMDiS {
removeMacroElements();
/// === Reset all DOFAdmins of the mesh. ===
// === Reset all DOFAdmins of the mesh. ===
int nAdmins = mesh->getNumberOfDOFAdmin();
for (int i = 0; i < nAdmins; i++) {
......@@ -85,21 +81,20 @@ namespace AMDiS {
for (int j = 0; j < static_cast<int>(mapLocalGlobalDOFs.size()); j++)
admin.setDOFFree(j, false);
admin.setUsedSize(mapLocalGlobalDOFs.size() - 1);
admin.setUsedSize(mapLocalGlobalDOFs.size());
admin.setUsedCount(mapLocalGlobalDOFs.size());
admin.setFirstHole(mapLocalGlobalDOFs.size());
}
// === Global refinements. ===
/// === Global refinements. ===
refinementManager->globalRefine(mesh, 1);
// refinementManager->globalRefine(mesh, 1);
updateLocalGlobalNumbering(nRankDOFs, nOverallDOFs);
// updateLocalGlobalNumbering(nRankDOFs, nOverallDOFs);
exit(0);
// exit(0);
/// === Create petsc matrix. ===
// === Create petsc matrix. ===
int ierr;
ierr = MatCreate(PETSC_COMM_WORLD, &petscMatrix);
......@@ -147,13 +142,14 @@ namespace AMDiS {
MatAssemblyBegin(petscMatrix, MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(petscMatrix, MAT_FINAL_ASSEMBLY);
DOFVector<double>::Iterator dofIt(vec, USED_DOFS);
for (dofIt.reset(); !dofIt.end(); ++dofIt) {
int index = mapLocalGlobalDOFs[dofIt.getDOFIndex()];
double value = *dofIt;
VecSetValues(petscRhsVec, 1, &index, &value, ADD_VALUES);
}
}
}
......@@ -185,35 +181,40 @@ namespace AMDiS {
std::vector<double*> recvBuffers(recvDofs.size());
int i = 0;
for (std::map<int, std::vector<DegreeOfFreedom> >::iterator sendIt = sendDofs.begin();
for (std::map<int, std::vector<const DegreeOfFreedom*> >::iterator
sendIt = sendDofs.begin();
sendIt != sendDofs.end();
++sendIt, i++) {
sendBuffers[i] = new double[sendIt->second.size()];
int nSendDOFs = sendIt->second.size();
sendBuffers[i] = new double[nSendDOFs];
for (int j = 0; j < static_cast<int>(sendIt->second.size()); j++)
sendBuffers[i][j] = (*vec)[(sendIt->second)[j]];
for (int j = 0; j < nSendDOFs; j++)
sendBuffers[i][j] = (*vec)[(sendIt->second)[j][0]];
mpiComm.Isend(sendBuffers[i], sendIt->second.size(), MPI_DOUBLE, sendIt->first, 0);
mpiComm.Isend(sendBuffers[i], nSendDOFs, MPI_DOUBLE, sendIt->first, 0);
}
i = 0;
for (std::map<int, std::vector<DegreeOfFreedom> >::iterator recvIt = recvDofs.begin();
for (std::map<int, std::vector<const DegreeOfFreedom*> >::iterator
recvIt = recvDofs.begin();
recvIt != recvDofs.end();
++recvIt, i++) {
recvBuffers[i] = new double[recvIt->second.size()];
int nRecvDOFs = recvIt->second.size();
recvBuffers[i] = new double[nRecvDOFs];
mpiComm.Irecv(recvBuffers[i], recvIt->second.size(), MPI_DOUBLE, recvIt->first, 0);
mpiComm.Irecv(recvBuffers[i], nRecvDOFs, MPI_DOUBLE, recvIt->first, 0);
}
mpiComm.Barrier();
i = 0;
for (std::map<int, std::vector<DegreeOfFreedom> >::iterator recvIt = recvDofs.begin();
for (std::map<int, std::vector<const DegreeOfFreedom*> >::iterator
recvIt = recvDofs.begin();
recvIt != recvDofs.end();
++recvIt, i++) {
for (int j = 0; j < static_cast<int>(recvIt->second.size()); j++)
(*vec)[(recvIt->second)[j]] = recvBuffers[i][j];
(*vec)[(recvIt->second)[j][0]] = recvBuffers[i][j];
delete [] recvBuffers[i];
}
......@@ -324,7 +325,7 @@ namespace AMDiS {
bool isRankDOF2 = (find(rankDOFs.begin(), rankDOFs.end(), boundDOF2) != rankDOFs.end());
bool ranksBoundary = isRankDOF1 || isRankDOF2;
/// === And add the part of the interior boundary. ===
// === And add the part of the interior boundary. ===
AtomicBoundary& bound =
(ranksBoundary ?
......@@ -337,7 +338,6 @@ namespace AMDiS {
bound.neighbourObject.el = elInfo->getNeighbour(i);
bound.neighbourObject.subObjAtBoundary = EDGE;
bound.neighbourObject.ithObjAtBoundary = -1;
std::cout << "ADD IN " << mpiRank << ": " << element->getIndex() << " " << elInfo->getNeighbour(i)->getIndex() << std::endl;
}
}
}
......@@ -393,7 +393,8 @@ namespace AMDiS {
break;
// The element must always be found, because the list is just in another order.
TEST_EXIT(k < rankIt->second.size())("Should never happen!\n");
TEST_EXIT(k < static_cast<int>(rankIt->second.size()))
("Should never happen!\n");
// Swap the current with the found element.
AtomicBoundary tmpBound = (rankIt->second)[k];
......@@ -427,12 +428,15 @@ namespace AMDiS {
}
void ParallelDomainProblemBase::createLocalGlobalNumbering(std::vector<const DegreeOfFreedom*>& rankDOFs,
std::map<const DegreeOfFreedom*, int>& boundaryDOFs,
int& nRankDOFs,
int& nOverallDOFs)
void ParallelDomainProblemBase::createLocalGlobalNumbering(
std::vector<const DegreeOfFreedom*>& rankDOFs,
std::map<const DegreeOfFreedom*, int>& boundaryDOFs,
int& nRankDOFs,
int& nOverallDOFs)
{
/// === Get rank information about DOFs. ===
FUNCNAME("ParallelDomainProblemBase::createLocalGlobalNumbering()");
// === Get rank information about DOFs. ===
// Stores to each DOF pointer the set of ranks the DOF is part of.
std::map<const DegreeOfFreedom*, std::set<int> > partitionDOFs;
......@@ -442,17 +446,25 @@ namespace AMDiS {
nRankDOFs = rankDOFs.size();
nOverallDOFs = partitionDOFs.size();
// === Get starting position for global rank dof ordering ====
// === Get starting position for global rank dof ordering. ====
int rstart = 0;
MPI_Scan(&nRankDOFs, &rstart, 1, MPI_INT, MPI_SUM, PETSC_COMM_WORLD);
rstart -= nRankDOFs;
/// === Create information which dof indices must be send and which must be received. ===
// === Create information which dof indices must be send and which must ===
// === be received. ===
// Stores to each rank a map from DOF pointers (which are all owned by the rank
// and lie on an interior boundary) to new global DOF indices.
std::map<int, std::map<const DegreeOfFreedom*, DegreeOfFreedom> > sendNewDofs;
std::map<int, std::map<DegreeOfFreedom, DegreeOfFreedom> > sendNewDofs;
std::map<int, std::vector<DegreeOfFreedom> > recvNewDofs;
// Maps to each rank the number of new DOF indices this rank will receive from.
// All these DOFs are at an interior boundary on this rank, but are owned by
// another rank.
std::map<int, int> recvNewDofs;
for (std::map<const DegreeOfFreedom*, int>::iterator it = boundaryDOFs.begin();
it != boundaryDOFs.end();
......@@ -461,11 +473,9 @@ namespace AMDiS {
if (it->second == mpiRank) {
// If the boundary dof is a rank dof, it must be send to other ranks.
// old global index
int oldDofIndex = (it->first)[0];
// search for new dof index in ranks partition for this boundary dof
int newDofIndex = 0;
for (int i = 0; i < static_cast<int>(rankDOFs.size()); i++) {
DegreeOfFreedom newDofIndex = 0;
for (int i = 0; i < nRankDOFs; i++) {
if (rankDOFs[i] == it->first) {
newDofIndex = rstart + i;
break;
......@@ -477,64 +487,79 @@ namespace AMDiS {
itRanks != partitionDOFs[it->first].end();
++itRanks) {
if (*itRanks != mpiRank)
sendNewDofs[*itRanks][oldDofIndex] = newDofIndex;
sendNewDofs[*itRanks][it->first] = newDofIndex;
}
} else {
// If the boundary dof is not a rank dof, its new dof index, and later
// also the dof values, must be received from another rank.
recvNewDofs[it->second].push_back((it->first)[0]);
if (recvNewDofs.find(it->second) == recvNewDofs.end())
recvNewDofs[it->second] = 1;
else
recvNewDofs[it->second]++;
}
}
/// === Send and receive the dof indices at boundary. ===
// === Send and receive the dof indices at boundary. ===
std::vector<int*> sendBuffers(sendNewDofs.size());
std::vector<int*> recvBuffers(recvNewDofs.size());
sendDofs.clear();
recvDofs.clear();
int i = 0;
for (std::map<int, std::map<DegreeOfFreedom, DegreeOfFreedom> >::iterator sendIt = sendNewDofs.begin();
for (std::map<int, std::map<const DegreeOfFreedom*, DegreeOfFreedom> >::iterator
sendIt = sendNewDofs.begin();
sendIt != sendNewDofs.end();
++sendIt, i++) {
sendBuffers[i] = new int[sendIt->second.size() * 2];
int nSendDOFs = sendIt->second.size() * 2;
sendBuffers[i] = new int[nSendDOFs];
int c = 0;
for (std::map<DegreeOfFreedom, DegreeOfFreedom>::iterator dofIt = sendIt->second.begin();
for (std::map<const DegreeOfFreedom*, DegreeOfFreedom>::iterator
dofIt = sendIt->second.begin();
dofIt != sendIt->second.end();
++dofIt, c += 2) {
sendBuffers[i][c] = dofIt->first;
sendBuffers[i][c + 1] = dofIt->second;
++dofIt) {
sendBuffers[i][c++] = (dofIt->first)[0];
sendBuffers[i][c++] = dofIt->second;
sendDofs[sendIt->first].push_back(dofIt->second);
sendDofs[sendIt->first].push_back(dofIt->first);
}
mpiComm.Isend(sendBuffers[i], sendIt->second.size() * 2, MPI_INT, sendIt->first, 0);
mpiComm.Isend(sendBuffers[i], nSendDOFs, MPI_INT, sendIt->first, 0);
}
i = 0;
for (std::map<int, std::vector<DegreeOfFreedom> >::iterator recvIt = recvNewDofs.begin();
for (std::map<int, int>::iterator recvIt = recvNewDofs.begin();
recvIt != recvNewDofs.end();
++recvIt, i++) {
recvBuffers[i] = new int[recvIt->second.size() * 2];
int nRecvDOFs = recvIt->second * 2;
recvBuffers[i] = new int[nRecvDOFs];
mpiComm.Irecv(recvBuffers[i], recvIt->second.size() * 2, MPI_INT, recvIt->first, 0);
mpiComm.Irecv(recvBuffers[i], nRecvDOFs, MPI_INT, recvIt->first, 0);
}
mpiComm.Barrier();
/// === Delete send buffers. ===
// === Delete send buffers. ===
i = 0;
for (std::map<int, std::map<DegreeOfFreedom, DegreeOfFreedom> >::iterator sendIt = sendNewDofs.begin();
for (std::map<int, std::map<const DegreeOfFreedom*, DegreeOfFreedom> >::iterator
sendIt = sendNewDofs.begin();
sendIt != sendNewDofs.end();
++sendIt, i++)
delete [] sendBuffers[i];
/// === Change dof indices for rank partition. ===
// === Change dof indices for rank partition. ===
mapLocalGlobalDOFs.clear();
mapGlobalLocalDOFs.clear();
isRankDOF.clear();
for (int i = 0; i < static_cast<int>(rankDOFs.size()); i++) {
for (int i = 0; i < nRankDOFs; i++) {
const_cast<DegreeOfFreedom*>(rankDOFs[i])[0] = i;
mapLocalGlobalDOFs[i] = rstart + i;
mapGlobalLocalDOFs[rstart + i] = i;
......@@ -542,55 +567,67 @@ namespace AMDiS {
}
/// === Change dof indices at boundary from other ranks. ===
// === Change dof indices at boundary from other ranks. ===
i = 0;
for (std::map<int, std::vector<DegreeOfFreedom> >::iterator recvIt = recvNewDofs.begin();
for (std::map<int, int>::iterator recvIt = recvNewDofs.begin();
recvIt != recvNewDofs.end();
++recvIt, i++) {
for (int j = 0; j < static_cast<int>(recvIt->second.size()); j++) {
for (int j = 0; j < recvIt->second; j++) {
int oldDof = recvBuffers[i][j * 2];
int newDof = recvBuffers[i][j * 2 + 1];
int newLocalDof = mapLocalGlobalDOFs.size();
DegreeOfFreedom oldDof = recvBuffers[i][j * 2];
DegreeOfFreedom newGlobalDof = recvBuffers[i][j * 2 + 1];
DegreeOfFreedom newLocalDof = mapLocalGlobalDOFs.size();
recvDofs[recvIt->first].push_back(newDof);
bool found = false;
for (std::map<const DegreeOfFreedom*, int>::iterator dofIt = boundaryDOFs.begin();
dofIt != boundaryDOFs.end();
++dofIt) {
if ((dofIt->first)[0] == oldDof) {
recvDofs[recvIt->first].push_back(dofIt->first);
const_cast<DegreeOfFreedom*>(dofIt->first)[0] = newLocalDof;
mapLocalGlobalDOFs[newLocalDof] = newDof;
mapGlobalLocalDOFs[newDof] = newLocalDof;
mapLocalGlobalDOFs[newLocalDof] = newGlobalDof;
mapGlobalLocalDOFs[newGlobalDof] = newLocalDof;
isRankDOF[newLocalDof] = false;
found = true;
break;
}
}
TEST_EXIT(found)("Should not happen!\n");
}
delete [] recvBuffers[i];
}
#if 1
// === Create local information from sendDofs and recvDofs
/// === Create local information from sendDofs and recvDofs
for (std::map<int, std::vector<DegreeOfFreedom> >::iterator it = sendDofs.begin();
for (std::map<int, std::vector<const DegreeOfFreedom*> >::iterator
it = sendDofs.begin();
it != sendDofs.end();
++it)
for (std::vector<DegreeOfFreedom>::iterator dofIt = it->second.begin();
for (std::vector<const DegreeOfFreedom*>::iterator dofIt = it->second.begin();
dofIt != it->second.end();
dofIt++)
*dofIt = mapGlobalLocalDOFs[*dofIt];
dofIt++) {
// std::cout << "AND SET S " << (*dofIt)[0] << " = " << mapGlobalLocalDOFs[(*dofIt)[0]] << std::endl;
// const_cast<DegreeOfFreedom*>(*dofIt)[0] = mapGlobalLocalDOFs[(*dofIt)[0]];
}
for (std::map<int, std::vector<DegreeOfFreedom> >::iterator it = recvDofs.begin();
for (std::map<int, std::vector<const DegreeOfFreedom*> >::iterator
it = recvDofs.begin();
it != recvDofs.end();
++it)
for (std::vector<DegreeOfFreedom>::iterator dofIt = it->second.begin();
for (std::vector<const DegreeOfFreedom*>::iterator dofIt = it->second.begin();
dofIt != it->second.end();
dofIt++)
*dofIt = mapGlobalLocalDOFs[*dofIt];
dofIt++) {
// std::cout << "AND SET R " << (*dofIt)[0] << " = " << mapGlobalLocalDOFs[(*dofIt)[0]] << std::endl;
// const_cast<DegreeOfFreedom*>(*dofIt)[0] = mapGlobalLocalDOFs[(*dofIt)[0]];
}
#endif
}
......@@ -600,10 +637,9 @@ namespace AMDiS {
FUNCNAME("ParallelDomainProblemBase::updateLocalGlobalNumbering()");
std::set<const DegreeOfFreedom*> rankDOFs;
std::map<const DegreeOfFreedom*, int> boundaryDOFs;
std::map<const DegreeOfFreedom*, int> newBoundaryDOFs;
/// === Get all DOFs in ranks partition. ===
// === Get all DOFs in ranks partition. ===
TraverseStack stack;
ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
......@@ -620,6 +656,9 @@ namespace AMDiS {
// === Traverse on interior boundaries and move all not ranked owned DOFs from ===
// === rankDOFs to boundaryDOFs ===
std::map<int, std::vector<const DegreeOfFreedom*> > sendNewDOFs;
std::map<int, std::vector<const DegreeOfFreedom*> > recvNewDOFs;
for (std::map<int, std::vector<AtomicBoundary> >::iterator it =
myIntBoundary.boundary.begin();
it != myIntBoundary.boundary.end();
......@@ -647,72 +686,147 @@ namespace AMDiS {
ERROR_EXIT("Should never happen!\n");
}
TEST_EXIT(boundaryDOFs.find(dof1) != boundaryDOFs.end())
("Should never happen!\n");
TEST_EXIT(boundaryDOFs.find(dof2) != boundaryDOFs.end())
("Should never happen!\n");
newBoundaryDOFs[dof1] = boundaryDOFs[dof1];
newBoundaryDOFs[dof2] = boundaryDOFs[dof2];
std::vector<const DegreeOfFreedom*> boundDOFs;
addAllDOFs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary,
boundDOFs);
for (int i = 0; i < static_cast<int>(boundDOFs.size()); i++) {
newBoundaryDOFs[boundDOFs[i]] = mpiRank;
sendNewDOFs[it->first].push_back(boundDOFs[i]);
}
}
}
}
void ParallelDomainProblemBase::addAllDOFs(Element *el, int ithEdge,
std::vector<const DegreeOfFreedom*>& dofs,
bool addVertices)
{
FUNCNAME("ParallelDomainProblemBase::addAllDOFs()");
for (std::map<int, std::vector<AtomicBoundary> >::iterator it =
otherIntBoundary.boundary.begin();
it != otherIntBoundary.boundary.end();
++it) {
for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
boundIt != it->second.end();
++boundIt) {
const DegreeOfFreedom *dof1 = NULL;
const DegreeOfFreedom *dof2 = NULL;
switch (boundIt->rankObject.ithObjAtBoundary) {
case 0:
dof1 = boundIt->rankObject.el->getDOF(1);
dof2 = boundIt->rankObject.el->getDOF(2);
break;
case 1:
dof1 = boundIt->rankObject.el->getDOF(0);
dof2 = boundIt->rankObject.el->getDOF(2);
break;
case 2:
dof1 = boundIt->rankObject.el->getDOF(0);
dof2 = boundIt->rankObject.el->getDOF(1);
break;
default:
ERROR_EXIT("Should never happen!\n");
}
TEST_EXIT(boundaryDOFs.find(dof1) != boundaryDOFs.end())
("Should never happen!\n");
TEST_EXIT(boundaryDOFs.find(dof2) != boundaryDOFs.end())
("Should never happen!\n");
rankDOFs.erase(dof1);
rankDOFs.erase(dof2);
newBoundaryDOFs[dof1] = boundaryDOFs[dof1];
newBoundaryDOFs[dof2] = boundaryDOFs[dof2];
const DegreeOfFreedom* boundDOF1 = NULL;
const DegreeOfFreedom* boundDOF2 = NULL;
if (addVertices) {
switch (ithEdge) {
case 0:
boundDOF1 = el->getDOF(1);
boundDOF2 = el->getDOF(2);
break;
case 1:
boundDOF1 = el->getDOF(0);
boundDOF2 = el->getDOF(2);
break;
case 2:
boundDOF1 = el->getDOF(0);
boundDOF2 = el->getDOF(1);
break;
default:
ERROR_EXIT("Should never happen!\n");
std::vector<const DegreeOfFreedom*> boundDOFs;
addAllDOFs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary,
boundDOFs);
for (int i = 0; i < static_cast<int>(boundDOFs.size()); i++) {
rankDOFs.erase(boundDOFs[i]);
newBoundaryDOFs[boundDOFs[i]] = it->first;
recvNewDOFs[it->first].push_back(boundDOFs[i]);
}
}
dofs.push_back(boundDOF1);