Commit 10141fab authored by Thomas Witkowski's avatar Thomas Witkowski

Multiple FE spaces in paralle computation. Work now for 2D with no periodic bc.

parent 5ba917be
......@@ -1260,12 +1260,6 @@ namespace AMDiS {
{
FUNCNAME("Mesh::checkParallelMacroFile()");
// In parallel computations, first the mesh must be checked if it is refined
// in an apropriate way.
TEST_EXIT(admin.size() == 1)("Not yet implemented!\n");
// === Create a temporary mesh and load the macro file to it. ===
Mesh testMesh(name, dim);
......
......@@ -208,7 +208,10 @@ namespace AMDiS {
switch (bound.subObj) {
case VERTEX:
dofs.push_back(dof[bound.ithObj]);
{
int n0 = feSpace->getAdmin()->getNumberOfPreDofs(VERTEX);
dofs.push_back(&(dof[bound.ithObj][n0]));
}
break;
case EDGE:
getNodeDofsAtEdge(feSpace, bound, dofs);
......@@ -250,6 +253,7 @@ namespace AMDiS {
case 2:
case 3:
{
int n0 = feSpace->getAdmin()->getNumberOfPreDofs(VERTEX);
BoundaryObject nextBound0 = bound, nextBound1 = bound;
prepareNextBound(nextBound0, 0);
prepareNextBound(nextBound1, 1);
......@@ -266,14 +270,14 @@ namespace AMDiS {
child[1]->getNodeDofs(feSpace, nextBound1, dofs);
if (addDof)
dofs.push_back(child[0]->getDof(3));
dofs.push_back(&(child[0]->getDof()[3][n0]));
child[0]->getNodeDofs(feSpace, nextBound0, dofs);
} else {
child[0]->getNodeDofs(feSpace, nextBound0, dofs);
if (addDof)
dofs.push_back(child[0]->getDof(3));
dofs.push_back(&(child[0]->getDof()[3][n0]));
child[1]->getNodeDofs(feSpace, nextBound1, dofs);
}
......@@ -294,6 +298,7 @@ namespace AMDiS {
if (!child[0])
return;
int n0 = feSpace->getAdmin()->getNumberOfPreDofs(VERTEX);
BoundaryObject nextBound0 = bound, nextBound1 = bound;
prepareNextBound(nextBound0, 0);
prepareNextBound(nextBound1, 1);
......@@ -305,17 +310,18 @@ namespace AMDiS {
if (bound.reverseMode) {
child[1]->getNodeDofs(feSpace, nextBound0, dofs);
dofs.push_back(child[0]->getDof(3));
dofs.push_back(&(child[0]->getDof()[3][n0]));
child[0]->getNodeDofs(feSpace, nextBound1, dofs);
} else {
child[0]->getNodeDofs(feSpace, nextBound0, dofs);
dofs.push_back(child[0]->getDof(3));
dofs.push_back(&(child[0]->getDof()[3][n0]));
child[1]->getNodeDofs(feSpace, nextBound1, dofs);
}
break;
case 5:
TEST_EXIT_DBG(nextBound0.ithObj == nextBound1.ithObj)("Should not happen!\n");
TEST_EXIT_DBG(nextBound0.ithObj == nextBound1.ithObj)
("Should not happen!\n");
if (nextBound0.ithObj != -1)
child[0]->getNodeDofs(feSpace, nextBound0, dofs);
......
......@@ -93,8 +93,11 @@ namespace AMDiS {
{
FUNCNAME("Triangle::getNodeDofs()");
// Get displacement if more than one FE space is defined on mesh.
int n0 = feSpace->getAdmin()->getNumberOfPreDofs(VERTEX);
if (bound.subObj == VERTEX) {
dofs.push_back(dof[bound.ithObj]);
dofs.push_back(&(dof[bound.ithObj][n0]));
return;
}
......@@ -106,16 +109,18 @@ namespace AMDiS {
case 0:
{
if (child[1] && child[1]->getFirstChild()) {
const DegreeOfFreedom** elDofs = child[1]->getFirstChild()->getDof();
if (bound.reverseMode) {
nextBound.ithObj = 1;
child[1]->getSecondChild()->getNodeDofs(feSpace, nextBound, dofs);
dofs.push_back(child[1]->getFirstChild()->getDof(2));
dofs.push_back(&(elDofs[2][n0]));
nextBound.ithObj = 0;
child[1]->getFirstChild()->getNodeDofs(feSpace, nextBound, dofs);
} else {
nextBound.ithObj = 0;
child[1]->getFirstChild()->getNodeDofs(feSpace, nextBound, dofs);
dofs.push_back(child[1]->getFirstChild()->getDof(2));
dofs.push_back(&(elDofs[2][n0]));
nextBound.ithObj = 1;
child[1]->getSecondChild()->getNodeDofs(feSpace, nextBound, dofs);
}
......@@ -125,16 +130,18 @@ namespace AMDiS {
case 1:
{
if (child[0] && child[0]->getFirstChild()) {
const DegreeOfFreedom** elDofs = child[0]->getFirstChild()->getDof();
if (bound.reverseMode) {
nextBound.ithObj = 1;
child[0]->getSecondChild()->getNodeDofs(feSpace, nextBound, dofs);
dofs.push_back(child[0]->getFirstChild()->getDof(2));
dofs.push_back(&(elDofs[2][n0]));
nextBound.ithObj = 0;
child[0]->getFirstChild()->getNodeDofs(feSpace, nextBound, dofs);
} else {
nextBound.ithObj = 0;
child[0]->getFirstChild()->getNodeDofs(feSpace, nextBound, dofs);
dofs.push_back(child[0]->getFirstChild()->getDof(2));
dofs.push_back(&(elDofs[2][n0]));
nextBound.ithObj = 1;
child[0]->getSecondChild()->getNodeDofs(feSpace, nextBound, dofs);
}
......@@ -143,16 +150,18 @@ namespace AMDiS {
break;
case 2:
if (child[0]) {
const DegreeOfFreedom** elDofs = child[0]->getDof();
if (bound.reverseMode) {
nextBound.ithObj = 1;
child[1]->getNodeDofs(feSpace, nextBound, dofs);
dofs.push_back(child[0]->getDof(2));
dofs.push_back(&(elDofs[2][n0]));
nextBound.ithObj = 0;
child[0]->getNodeDofs(feSpace, nextBound, dofs);
} else {
nextBound.ithObj = 0;
child[0]->getNodeDofs(feSpace, nextBound, dofs);
dofs.push_back(child[0]->getDof(2));
dofs.push_back(&(elDofs[2][n0]));
nextBound.ithObj = 1;
child[1]->getNodeDofs(feSpace, nextBound, dofs);
}
......
......@@ -247,6 +247,10 @@ namespace AMDiS {
#if HAVE_PARALLEL_DOMAIN_AMDIS
paraFilename += timeStr;
postfix += timeStr + paraviewFileExt;
#endif
} else {
#if HAVE_PARALLEL_DOMAIN_AMDIS
postfix += paraviewFileExt;
#endif
}
......
......@@ -25,13 +25,13 @@ namespace AMDiS {
int nElements = (*dataCollector)[0]->getNumberElements();
int vertices = (*dataCollector)[0]->getMesh()->getGeo(VERTEX);
if ((dim == 2) && (degree == 2)) {
if (dim == 2 && degree == 2) {
nVertices += (*dataCollector)[0]->getNumberInterpPoints();
nElements *= 4;
} else if ((dim == 2) && (degree == 3)) {
} else if (dim == 2 && degree == 3) {
nVertices += (*dataCollector)[0]->getNumberInterpPoints();
nElements *= 9;
} else if ((dim == 2) && (degree == 4)) {
} else if (dim == 2 && degree == 4) {
nVertices += (*dataCollector)[0]->getNumberInterpPoints();
nElements *= 16;
}
......
......@@ -9,6 +9,8 @@
//
// See also license.opensource.txt in the distribution.
#ifdef HAVE_BDDC_ML
extern "C" {
#include <bddcml_interface_c.h>
}
......@@ -18,8 +20,6 @@ extern "C" {
namespace AMDiS {
#ifdef HAVE_BDDC_ML
void BddcMlSolver::fillPetscMatrix(Matrix<DOFMatrix*> *m)
{
FUNCNAME("BddcMlSolver::fillPetscMatrix()");
......@@ -339,6 +339,7 @@ namespace AMDiS {
}
}
#endif
}
......@@ -97,6 +97,19 @@ namespace AMDiS {
++feMapIter;
}
inline void next()
{
++feMapIter;
if (feMapIter == dataIter->second.end()) {
++dataIter;
setNextFeMap();
} else {
dofIter = feMapIter->second.begin();
dofCounter = 0;
}
}
inline void beginDofIter(const FiniteElemSpace *fe = NULL)
{
FUNCNAME("DofComm::Iterator::beginDofIter()");
......
......@@ -438,6 +438,12 @@ namespace AMDiS {
return faceReverseMode[make_pair(obj0, obj1)];
}
/// Returns true if there is periodic data.
bool hasPeriodicData()
{
return (periodicVertices.size() != 0);
}
/// Write the element database to disk.
void serialize(ostream &out);
......
......@@ -385,11 +385,6 @@ namespace AMDiS {
mesh = feSpaces[0]->getMesh();
info = probStat->getInfo();
TEST_EXIT(mesh->getNumberOfDOFAdmin() == 1)
("Only meshes with one DOFAdmin are supported!\n");
TEST_EXIT(mesh->getDofAdmin(0).getNumberOfPreDofs(VERTEX) == 0)
("Wrong pre dof number for DOFAdmin!\n");
switch (mesh->getDim()) {
case 2:
refineManager = new RefinementManager2d();
......@@ -1448,8 +1443,10 @@ namespace AMDiS {
// === Create periodic data, if there are periodic boundary conditions. ===
TEST_EXIT(feSpaces.size() == 1)
("Sebastian: Na, dass funktioniert auch noch nicht mit mehreren FE spaces. Du weisst schon, wen du jetzt mobben kannst :)!\n");
if (elObjects.hasPeriodicData()) {
TEST_EXIT(feSpaces.size() == 1)
("Sebastian: Na, dass funktioniert auch noch nicht mit mehreren FE spaces. Du weisst schon, wen du jetzt mobben kannst :)!\n");
}
elObjects.createPeriodicData(feSpaces[0]);
// === Create data about the reverse modes of neighbouring elements. ===
......@@ -1830,16 +1827,18 @@ namespace AMDiS {
}
}
} else {
for (InteriorBoundary::iterator it(myIntBoundary); !it.end(); ++it)
for (InteriorBoundary::iterator it(myIntBoundary); !it.end(); ++it) {
it->rankObj.el->getAllDofs(feSpace, it->rankObj,
sendDofs.getDofCont(it.getRank(), feSpace));
}
for (InteriorBoundary::iterator it(otherIntBoundary); !it.end(); ++it)
for (InteriorBoundary::iterator it(otherIntBoundary); !it.end(); ++it) {
it->rankObj.el->getAllDofs(feSpace, it->rankObj,
recvDofs.getDofCont(it.getRank(), feSpace));
}
}
// === Delete all empty DOF send and recv positions ===
sendDofs.removeEmpty();
......@@ -1920,7 +1919,6 @@ namespace AMDiS {
sort(rankDofs.begin(), rankDofs.end(), cmpDofsByValue);
int nRankAllDofs = rankDofs.size();
// === Traverse interior boundaries and get all DOFs on them. ===
createBoundaryDofs(feSpace);
......
......@@ -57,13 +57,9 @@ namespace AMDiS {
// === Create PETSc matrix with the computed nnz data structure. ===
// MatCreateMPIAIJ(PETSC_COMM_WORLD, nRankRows, nRankRows,
// nOverallRows, nOverallRows,
// 0, d_nnz, 0, o_nnz, &petscMatrix);
MatCreateMPIAIJ(PETSC_COMM_WORLD, nRankRows, nRankRows,
nOverallRows, nOverallRows,
30, PETSC_NULL, 30, PETSC_NULL, &petscMatrix);
MatCreateMPIAIJ(PETSC_COMM_WORLD, nRankRows, nRankRows,
nOverallRows, nOverallRows,
0, d_nnz, 0, o_nnz, &petscMatrix);
#if (DEBUG != 0)
MSG("Fill petsc matrix 1 needed %.5f seconds\n", MPI::Wtime() - wtime);
......@@ -138,7 +134,6 @@ namespace AMDiS {
{
FUNCNAME("PetscSolverGlobalMatrix::solvePetscMatrix()");
const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);
int nComponents = vec.getSize();
// === Set old solution to be initiual guess for PETSc solver. ===
......@@ -162,7 +157,8 @@ namespace AMDiS {
int c = 0;
for (int i = 0; i < nComponents; i++) {
DOFVector<double> &dv = *(vec.getDOFVector(i));
int nRankDofs = meshDistributor->getNumberRankDofs(dv.getFeSpace());
const FiniteElemSpace *feSpace = dv.getFeSpace();
int nRankDofs = meshDistributor->getNumberRankDofs(feSpace);
for (int j = 0; j < nRankDofs; j++)
dv[meshDistributor->mapLocalToDofIndex(feSpace, j)] = vecPointer[c++];
......@@ -474,9 +470,9 @@ namespace AMDiS {
TEST_EXIT_DBG(!d_nnz)("There is something wrong!\n");
TEST_EXIT_DBG(!o_nnz)("There is something wrong!\n");
const FiniteElemSpace* feSpace = meshDistributor->getFeSpace(0);
int nComponents = mat->getNumRows();
int nRankRows = meshDistributor->getNumberRankDofs(feSpace) * nComponents;
vector<const FiniteElemSpace*> feSpaces = getFeSpaces(mat);
int nRankRows = meshDistributor->getNumberRankDofs(feSpaces);
int rankStartIndex = meshDistributor->getStartDofs(feSpaces);
d_nnz = new int[nRankRows];
o_nnz = new int[nRankRows];
for (int i = 0; i < nRankRows; i++) {
......@@ -490,38 +486,50 @@ namespace AMDiS {
typedef vector<pair<int, int> > MatrixNnzEntry;
typedef map<int, DofContainer> RankToDofContainer;
// Stores to each rank a list of nnz entries (i.e. pairs of row and column index)
// that this rank will send to. These nnz entries will be assembled on this rank,
// but because the row DOFs are not DOFs of this rank they will be send to the
// owner of the row DOFs.
// Stores to each rank a list of nnz entries (i.e. pairs of row and column
// index) that this rank will send to. These nnz entries will be assembled
// on this rank, but because the row DOFs are not DOFs of this rank they
// will be send to the owner of the row DOFs.
map<int, MatrixNnzEntry> sendMatrixEntry;
// Maps to each DOF that must be send to another rank the rank number of the
// receiving rank.
map<DegreeOfFreedom, int> sendDofToRank;
map<pair<DegreeOfFreedom, int>, int> sendDofToRank;
// First, create for all ranks we send data to MatrixNnzEntry object
// with 0 entries.
for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace);
!it.end(); it.nextRank()) {
sendMatrixEntry[it.getRank()].resize(0);
for (; !it.endDofIter(); it.nextDof())
sendDofToRank[it.getDofIndex()] = it.getRank();
// First, create for all ranks, to which we send data to, MatrixNnzEntry
// object with 0 entries.
for (unsigned int i = 0; i < feSpaces.size(); i++) {
for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpaces[i]);
!it.end(); it.nextRank()) {
sendMatrixEntry[it.getRank()].resize(0);
for (; !it.endDofIter(); it.nextDof())
sendDofToRank[make_pair(it.getDofIndex(), i)] = it.getRank();
}
}
// Create list of ranks from which we receive data from.
std::set<int> recvFromRank;
for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpace);
!it.end(); it.nextRank())
recvFromRank.insert(it.getRank());
for (unsigned int i = 0; i < feSpaces.size(); i++)
for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpaces[i]);
!it.end(); it.nextRank())
recvFromRank.insert(it.getRank());
// === Traverse matrices to create nnz data. ===
int nComponents = mat->getNumRows();
for (int i = 0; i < nComponents; i++) {
for (int j = 0; j < nComponents; j++) {
if (!(*mat)[i][j])
continue;
TEST_EXIT_DBG((*mat)[i][j]->getRowFeSpace() == feSpaces[i])
("Should not happen!\n");
TEST_EXIT_DBG((*mat)[i][j]->getColFeSpace() == feSpaces[j])
("Should not happen!\n");
Matrix bmat = (*mat)[i][j]->getBaseMatrix();
traits::col<Matrix>::type col(bmat);
......@@ -532,42 +540,25 @@ namespace AMDiS {
for (cursor_type cursor = begin<row>(bmat),
cend = end<row>(bmat); cursor != cend; ++cursor) {
int globalRowDof = meshDistributor->mapLocalToGlobal(feSpace, *cursor);
vector<int> rows;
rows.push_back(globalRowDof);
vector<int> rowRank;
if (meshDistributor->getIsRankDof(feSpace, *cursor)) {
rowRank.push_back(meshDistributor->getMpiRank());
} else {
// Find out who is the member of this DOF.
TEST_EXIT_DBG(sendDofToRank.count(*cursor))("DOF %d has no receiver!\n", *cursor);
rowRank.push_back(sendDofToRank[*cursor]);
}
// Map the local row number to the global DOF index and create from it
// the global PETSc row index of this DOF.
int petscRowIdx = globalRowDof * nComponents + i;
// The corresponding global matrix row index of the current row DOF.
int petscRowIdx = dofToGlobalIndex[make_pair(*cursor, i)];
if (meshDistributor->getIsRankDof(feSpace, *cursor)) {
if (meshDistributor->getIsRankDof(feSpaces[i], *cursor)) {
// === The current row DOF is a rank dof, so create the corresponding ===
// === nnz values directly on rank's nnz data. ===
// === The current row DOF is a rank DOF, so create the ===
// === corresponding nnz values directly on rank's nnz data. ===
// This is the local row index of the local PETSc matrix.
int localPetscRowIdx =
petscRowIdx - meshDistributor->getStartDofs(feSpace) * nComponents;
int localPetscRowIdx = petscRowIdx - rankStartIndex;
TEST_EXIT_DBG(localPetscRowIdx >= 0 && localPetscRowIdx < nRankRows)
("Should not happen! \n Debug info: localRowIdx = %d globalRowIndx = %d petscRowIdx = %d localPetscRowIdx = %d rStart = %d nCompontens = %d nRankRows = %d\n",
*cursor,
meshDistributor->mapLocalToGlobal(feSpace, *cursor),
meshDistributor->mapLocalToGlobal(feSpaces[i], *cursor),
petscRowIdx,
localPetscRowIdx,
meshDistributor->getStartDofs(feSpace),
rankStartIndex,
nComponents,
nRankRows);
......@@ -575,33 +566,33 @@ namespace AMDiS {
// Traverse all non zero entries in this row.
for (icursor_type icursor = begin<nz>(cursor),
icend = end<nz>(cursor); icursor != icend; ++icursor) {
int petscColIdx =
meshDistributor->mapLocalToGlobal(feSpace, col(*icursor)) * nComponents + j;
int petscColIdx = dofToGlobalIndex[make_pair(col(*icursor), j)];
if (value(*icursor) != 0.0 || petscRowIdx == petscColIdx) {
// The row DOF is a rank DOF, if also the column is a rank DOF,
// increment the d_nnz values for this row, otherwise the o_nnz value.
if (petscColIdx >= meshDistributor->getStartDofs(feSpace) * nComponents &&
petscColIdx < meshDistributor->getStartDofs(feSpace) * nComponents + nRankRows)
// increment the d_nnz values for this row, otherwise the
// o_nnz value.
if (petscColIdx >= rankStartIndex &&
petscColIdx < rankStartIndex + nRankRows)
d_nnz[localPetscRowIdx]++;
else
o_nnz[localPetscRowIdx]++;
}
}
} else {
// === The current row DOF is not a rank dof, i.e., it will be created ===
// === on this rank, but after this it will be send to another rank ===
// === matrix. So we need to send also the corresponding nnz structure ===
// === of this row to the corresponding rank. ===
// === The current row DOF is not a rank DOF, i.e., its values ===
// === are also created on this rank, but afterthere they will ===
// === be send to another rank. So we need to send also the ===
// === corresponding nnz structure of this row to the corres- ===
// === ponding rank. ===
// Send all non zero entries to the member of the row DOF.
int sendToRank = sendDofToRank[*cursor];
int sendToRank = sendDofToRank[make_pair(*cursor, i)];
for (icursor_type icursor = begin<nz>(cursor),
icend = end<nz>(cursor); icursor != icend; ++icursor) {
if (value(*icursor) != 0.0) {
int petscColIdx =
meshDistributor->mapLocalToGlobal(feSpace, col(*icursor)) * nComponents + j;
int petscColIdx = dofToGlobalIndex[make_pair(col(*icursor), j)];
sendMatrixEntry[sendToRank].
push_back(make_pair(petscRowIdx, petscColIdx));
......@@ -623,8 +614,8 @@ namespace AMDiS {
stdMpi.startCommunication();
// === Evaluate the nnz structure this rank got from other ranks and add it to ===
// === the PETSc nnz data structure. ===
// === Evaluate the nnz structure this rank got from other ranks and add ===
// === it to the PETSc nnz data structure. ===
for (map<int, MatrixNnzEntry>::iterator it = stdMpi.getRecvData().begin();
it != stdMpi.getRecvData().end(); ++it) {
......@@ -633,14 +624,13 @@ namespace AMDiS {
int r = it->second[i].first;
int c = it->second[i].second;
int localRowIdx = r - meshDistributor->getStartDofs(feSpace) * nComponents;
int localRowIdx = r - rankStartIndex;
TEST_EXIT_DBG(localRowIdx >= 0 && localRowIdx < nRankRows)
("Got row index %d/%d (nRankRows = %d) from rank %d. Should not happen!\n",
r, localRowIdx, nRankRows, it->first);
if (c < meshDistributor->getStartDofs(feSpace) * nComponents ||
c >= meshDistributor->getStartDofs(feSpace) * nComponents + nRankRows)
if (c < rankStartIndex || c >= rankStartIndex + nRankRows)
o_nnz[localRowIdx]++;
else
d_nnz[localRowIdx]++;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment