diff --git a/AMDiS/bin/Makefile.am b/AMDiS/bin/Makefile.am
index 222d304a83038c48752d01463ae789fb1429c6f4..17f54f5646c0adb4b31c70af0d08d656b48f0d20 100644
--- a/AMDiS/bin/Makefile.am
+++ b/AMDiS/bin/Makefile.am
@@ -216,6 +216,7 @@ $(SOURCE_DIR)/SubElInfo.h $(SOURCE_DIR)/SubElInfo.cc \
$(SOURCE_DIR)/SolutionDataStorage.h $(SOURCE_DIR)/SolutionDataStorage.hh \
$(SOURCE_DIR)/InteriorBoundary.h $(SOURCE_DIR)/InteriorBoundary.cc \
$(SOURCE_DIR)/ElementDofIteartor.h $(SOURCE_DIR)/ElementDofIterator.cc \
+$(SOURCE_DIR)/Debug.h $(SOURCE_DIR)/Debug.cc \
$(SOURCE_DIR)/parareal/ProblemBase.h \
$(SOURCE_DIR)/parareal/AdaptParaReal.h $(SOURCE_DIR)/parareal/AdaptParaReal.cc
diff --git a/AMDiS/bin/Makefile.in b/AMDiS/bin/Makefile.in
index 55098118ddead6f3567eb660676f4374bf2d78d4..50edae05072cf0766b4ff36f7c8aa4ab4eef89cd 100644
--- a/AMDiS/bin/Makefile.in
+++ b/AMDiS/bin/Makefile.in
@@ -225,8 +225,8 @@ am__libamdis_la_SOURCES_DIST = $(PARALLEL_DIR)/ParallelDomainBase.h \
$(SOURCE_DIR)/InteriorBoundary.h \
$(SOURCE_DIR)/InteriorBoundary.cc \
$(SOURCE_DIR)/ElementDofIteartor.h \
- $(SOURCE_DIR)/ElementDofIterator.cc \
- $(SOURCE_DIR)/parareal/ProblemBase.h \
+ $(SOURCE_DIR)/ElementDofIterator.cc $(SOURCE_DIR)/Debug.h \
+ $(SOURCE_DIR)/Debug.cc $(SOURCE_DIR)/parareal/ProblemBase.h \
$(SOURCE_DIR)/parareal/AdaptParaReal.h \
$(SOURCE_DIR)/parareal/AdaptParaReal.cc
@USE_PARALLEL_DOMAIN_AMDIS_TRUE@am__objects_1 = libamdis_la-ParallelDomainBase.lo \
@@ -290,7 +290,8 @@ am_libamdis_la_OBJECTS = $(am__objects_2) libamdis_la-DOFIndexed.lo \
libamdis_la-PngWriter.lo libamdis_la-PovrayWriter.lo \
libamdis_la-DataCollector.lo libamdis_la-ScalableQuadrature.lo \
libamdis_la-SubElInfo.lo libamdis_la-InteriorBoundary.lo \
- libamdis_la-ElementDofIterator.lo libamdis_la-AdaptParaReal.lo
+ libamdis_la-ElementDofIterator.lo libamdis_la-Debug.lo \
+ libamdis_la-AdaptParaReal.lo
libamdis_la_OBJECTS = $(am_libamdis_la_OBJECTS)
libcompositeFEM_la_LIBADD =
am_libcompositeFEM_la_OBJECTS = libcompositeFEM_la-CFE_Integration.lo \
@@ -629,6 +630,7 @@ $(SOURCE_DIR)/SubElInfo.h $(SOURCE_DIR)/SubElInfo.cc \
$(SOURCE_DIR)/SolutionDataStorage.h $(SOURCE_DIR)/SolutionDataStorage.hh \
$(SOURCE_DIR)/InteriorBoundary.h $(SOURCE_DIR)/InteriorBoundary.cc \
$(SOURCE_DIR)/ElementDofIteartor.h $(SOURCE_DIR)/ElementDofIterator.cc \
+$(SOURCE_DIR)/Debug.h $(SOURCE_DIR)/Debug.cc \
$(SOURCE_DIR)/parareal/ProblemBase.h \
$(SOURCE_DIR)/parareal/AdaptParaReal.h $(SOURCE_DIR)/parareal/AdaptParaReal.cc
@@ -747,6 +749,7 @@ distclean-compile:
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-DOFMatrix.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-DOFVector.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-DataCollector.Plo@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-Debug.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-DirichletBC.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-DualTraverse.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-ElInfo.Plo@am__quote@
@@ -1532,6 +1535,13 @@ libamdis_la-ElementDofIterator.lo: $(SOURCE_DIR)/ElementDofIterator.cc
@AMDEP_TRUE@@am__fastdepCXX_FALSE@ DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
@am__fastdepCXX_FALSE@ $(LIBTOOL) --tag=CXX --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libamdis_la_CXXFLAGS) $(CXXFLAGS) -c -o libamdis_la-ElementDofIterator.lo `test -f '$(SOURCE_DIR)/ElementDofIterator.cc' || echo '$(srcdir)/'`$(SOURCE_DIR)/ElementDofIterator.cc
+libamdis_la-Debug.lo: $(SOURCE_DIR)/Debug.cc
+@am__fastdepCXX_TRUE@ if $(LIBTOOL) --tag=CXX --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libamdis_la_CXXFLAGS) $(CXXFLAGS) -MT libamdis_la-Debug.lo -MD -MP -MF "$(DEPDIR)/libamdis_la-Debug.Tpo" -c -o libamdis_la-Debug.lo `test -f '$(SOURCE_DIR)/Debug.cc' || echo '$(srcdir)/'`$(SOURCE_DIR)/Debug.cc; \
+@am__fastdepCXX_TRUE@ then mv -f "$(DEPDIR)/libamdis_la-Debug.Tpo" "$(DEPDIR)/libamdis_la-Debug.Plo"; else rm -f "$(DEPDIR)/libamdis_la-Debug.Tpo"; exit 1; fi
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@ source='$(SOURCE_DIR)/Debug.cc' object='libamdis_la-Debug.lo' libtool=yes @AMDEPBACKSLASH@
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@ DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
+@am__fastdepCXX_FALSE@ $(LIBTOOL) --tag=CXX --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libamdis_la_CXXFLAGS) $(CXXFLAGS) -c -o libamdis_la-Debug.lo `test -f '$(SOURCE_DIR)/Debug.cc' || echo '$(srcdir)/'`$(SOURCE_DIR)/Debug.cc
+
libamdis_la-AdaptParaReal.lo: $(SOURCE_DIR)/parareal/AdaptParaReal.cc
@am__fastdepCXX_TRUE@ if $(LIBTOOL) --tag=CXX --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libamdis_la_CXXFLAGS) $(CXXFLAGS) -MT libamdis_la-AdaptParaReal.lo -MD -MP -MF "$(DEPDIR)/libamdis_la-AdaptParaReal.Tpo" -c -o libamdis_la-AdaptParaReal.lo `test -f '$(SOURCE_DIR)/parareal/AdaptParaReal.cc' || echo '$(srcdir)/'`$(SOURCE_DIR)/parareal/AdaptParaReal.cc; \
@am__fastdepCXX_TRUE@ then mv -f "$(DEPDIR)/libamdis_la-AdaptParaReal.Tpo" "$(DEPDIR)/libamdis_la-AdaptParaReal.Plo"; else rm -f "$(DEPDIR)/libamdis_la-AdaptParaReal.Tpo"; exit 1; fi
diff --git a/AMDiS/src/Debug.cc b/AMDiS/src/Debug.cc
new file mode 100644
index 0000000000000000000000000000000000000000..9b48fca8992114ae871bee4c6f3b21c00db5bbd0
--- /dev/null
+++ b/AMDiS/src/Debug.cc
@@ -0,0 +1,139 @@
+#include <boost/lexical_cast.hpp>
+#include "Debug.h"
+#include "DOFVector.h"
+#include "MacroElement.h"
+#include "VtkWriter.h"
+
+namespace AMDiS {
+
+ namespace debug {
+
+#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
+ void writeLocalElementDofs(int rank, int elIdx, FiniteElemSpace *feSpace)
+ {
+ using boost::lexical_cast;
+
+ if (MPI::COMM_WORLD.Get_rank() == rank) {
+ DOFVector<double> tmp(feSpace, "tmp");
+ colorDofVectorByLocalElementDofs(tmp, feSpace->getMesh(), elIdx);
+ VtkWriter::writeFile(tmp, "tmp" + lexical_cast<std::string>(elIdx) + ".vtu");
+ }
+ }
+
+ void writeDofMesh(int rank, DegreeOfFreedom dof, FiniteElemSpace *feSpace)
+ {
+ using boost::lexical_cast;
+
+ if (MPI::COMM_WORLD.Get_rank() == rank) {
+ DOFVector<double> tmp(feSpace, "tmp");
+ tmp.set(0.0);
+ tmp[dof] = 1.0;
+ VtkWriter::writeFile(tmp, "dofmesh" + lexical_cast<std::string>(rank) + ".vtu");
+ }
+ }
+
+ void writeMesh(FiniteElemSpace *feSpace, int rank, std::string filename)
+ {
+ using boost::lexical_cast;
+
+ int myRank = MPI::COMM_WORLD.Get_rank();
+ if (rank == -1 || myRank == rank) {
+ DOFVector<double> tmp(feSpace, "tmp");
+ VtkWriter::writeFile(tmp, filename + lexical_cast<std::string>(myRank) + ".vtu");
+ }
+ }
+#endif
+
+ void colorDofVectorByLocalElementDofs(DOFVector<double>& vec, Element *el)
+ {
+ // === Get local indices of the given element. ===
+
+ const BasisFunction *basisFcts = vec.getFESpace()->getBasisFcts();
+ int nBasisFcts = basisFcts->getNumber();
+ std::vector<DegreeOfFreedom> localDofs(nBasisFcts);
+ basisFcts->getLocalIndices(el, vec.getFESpace()->getAdmin(), localDofs);
+
+ // === Set the values of the dof vector. ===
+
+ vec.set(0.0);
+ for (int i = 0; i < nBasisFcts; i++)
+ vec[localDofs[i]] = static_cast<double>(i);
+ }
+
+ bool colorDofVectorByLocalElementDofs(DOFVector<double>& vec, Mesh *mesh,
+ int elIndex)
+ {
+ FUNCNAME("colorDofVectorByLocalElementDofs()");
+
+ TraverseStack stack;
+ ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
+ while (elInfo) {
+ if (elInfo->getElement()->getIndex() == elIndex) {
+ colorDofVectorByLocalElementDofs(vec, elInfo->getElement());
+ return true;
+ }
+ elInfo = stack.traverseNext(elInfo);
+ }
+
+ return false;
+ }
+
+ Element* getDofIndexElement(FiniteElemSpace *feSpace, DegreeOfFreedom dof)
+ {
+ const BasisFunction* basFcts = feSpace->getBasisFcts();
+ int nBasFcts = basFcts->getNumber();
+ std::vector<DegreeOfFreedom> dofVec(nBasFcts);
+
+ TraverseStack stack;
+ ElInfo *elInfo = stack.traverseFirst(feSpace->getMesh(), -1,
+ Mesh::CALL_EVERY_EL_PREORDER);
+ while (elInfo) {
+ basFcts->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(), dofVec);
+ for (int i = 0; i < nBasFcts; i++)
+ if (dofVec[i] == dof)
+ return elInfo->getElement();
+
+ elInfo = stack.traverseNext(elInfo);
+ }
+
+ return NULL;
+ }
+
+ Element* getLevel0ParentElement(Mesh *mesh, Element *el)
+ {
+ TraverseStack stack;
+ ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
+ while (elInfo) {
+ if (elInfo->getElement() == el)
+ return elInfo->getMacroElement()->getElement();
+
+ elInfo = stack.traverseNext(elInfo);
+ }
+
+ return NULL;
+ }
+
+ void printInfoByDof(FiniteElemSpace *feSpace, DegreeOfFreedom dof)
+ {
+ Element *el = getDofIndexElement(feSpace, dof);
+ Element *parEl = getLevel0ParentElement(feSpace->getMesh(), el);
+
+ std::cout << "DOF-INFO: dof = " << dof
+ << " elidx = " << el->getIndex()
+ << " pelidx = " << parEl->getIndex() << std::endl;
+
+ TraverseStack stack;
+ ElInfo *elInfo = stack.traverseFirst(feSpace->getMesh(), -1,
+ Mesh::CALL_EVERY_EL_PREORDER);
+ while (elInfo) {
+ if (elInfo->getElement()->getIndex() == parEl->getIndex())
+ std::cout << "EL INFO TO " << parEl->getIndex() << ": "
+ << elInfo->getType() << std::endl;
+
+ elInfo = stack.traverseNext(elInfo);
+ }
+ }
+
+ } // namespace debug
+
+} // namespace AMDiS
diff --git a/AMDiS/src/Debug.h b/AMDiS/src/Debug.h
new file mode 100644
index 0000000000000000000000000000000000000000..750b8f5ce44ce4fa3586a74126de3a8621cbc2a1
--- /dev/null
+++ b/AMDiS/src/Debug.h
@@ -0,0 +1,53 @@
+// ============================================================================
+// == ==
+// == AMDiS - Adaptive multidimensional simulations ==
+// == ==
+// ============================================================================
+// == ==
+// == TU Dresden ==
+// == ==
+// == Institut f�r Wissenschaftliches Rechnen ==
+// == Zellescher Weg 12-14 ==
+// == 01069 Dresden ==
+// == germany ==
+// == ==
+// ============================================================================
+// == ==
+// == https://gforge.zih.tu-dresden.de/projects/amdis/ ==
+// == ==
+// ============================================================================
+
+/** \file Debug.h */
+
+#ifndef AMDIS_DEBUG_H
+#define AMDIS_DEBUG_H
+
+#include "AMDiS_fwd.h"
+#include "Global.h"
+
+namespace AMDiS {
+
+ namespace debug {
+#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
+ void writeLocalElementDofs(int rank, int elIdx, FiniteElemSpace *feSpace);
+
+ void writeMesh(FiniteElemSpace *feSpace, int rank, std::string filename);
+
+ void writeDofMesh(int rank, DegreeOfFreedom dof, FiniteElemSpace *feSpace);
+#endif
+
+ void colorDofVectorByLocalElementDofs(DOFVector<double>& vec, Element *el);
+
+ bool colorDofVectorByLocalElementDofs(DOFVector<double>& vec,
+ Mesh *mesh,
+ int elIndex);
+
+ Element* getDofIndexElement(FiniteElemSpace *feSpace, DegreeOfFreedom dof);
+
+ Element* getLevel0ParentElement(Mesh *mesh, Element *el);
+
+ void printInfoByDof(FiniteElemSpace *feSpace, DegreeOfFreedom dof);
+ }
+}
+
+#endif
diff --git a/AMDiS/src/ParallelDomainBase.cc b/AMDiS/src/ParallelDomainBase.cc
index 1d7504653fe7289dd1e3cc74c8b9118583f22c20..083802e1ff2d3d0fb9e1e255e5c9e64914158607 100644
--- a/AMDiS/src/ParallelDomainBase.cc
+++ b/AMDiS/src/ParallelDomainBase.cc
@@ -1,5 +1,4 @@
#include <algorithm>
-#include <boost/lexical_cast.hpp>
#include "ParallelDomainBase.h"
#include "ParMetisPartitioner.h"
@@ -20,6 +19,7 @@
#include "VertexVector.h"
#include "StdMpi.h"
#include "MeshStructure.h"
+#include "Debug.h"
#include "petscksp.h"
@@ -119,6 +119,10 @@ namespace AMDiS {
removeMacroElements();
+ // === Reset all DOFAdmins of the mesh. ===
+
+ updateDofAdmins();
+
// === If in debug mode, make some tests. ===
#if (DEBUG != 0)
@@ -127,11 +131,9 @@ namespace AMDiS {
dbgTestInteriorBoundary();
dbgTestCommonDofs(true);
MSG("Debug mode tests finished!\n");
-#endif
- // === Reset all DOFAdmins of the mesh. ===
-
- updateDofAdmins();
+ debug::writeMesh(feSpace, -1, "macromesh");
+#endif
// === Create periodic dof mapping, if there are periodic boundaries. ===
@@ -167,17 +169,17 @@ namespace AMDiS {
// There must be always more allocated DOFs than used DOFs in DOFAdmin. Otherwise,
// it is not possible to define a first DOF hole.
- if (static_cast<unsigned int>(admin.getSize()) == mapLocalGlobalDOFs.size())
+ if (static_cast<unsigned int>(admin.getSize()) == mapLocalGlobalDofs.size())
admin.enlargeDOFLists();
for (int j = 0; j < admin.getSize(); j++)
admin.setDOFFree(j, true);
- for (unsigned int j = 0; j < mapLocalGlobalDOFs.size(); j++)
+ for (unsigned int j = 0; j < mapLocalGlobalDofs.size(); j++)
admin.setDOFFree(j, false);
- admin.setUsedSize(mapLocalGlobalDOFs.size());
- admin.setUsedCount(mapLocalGlobalDOFs.size());
- admin.setFirstHole(mapLocalGlobalDOFs.size());
+ admin.setUsedSize(mapLocalGlobalDofs.size());
+ admin.setUsedCount(mapLocalGlobalDofs.size());
+ admin.setFirstHole(mapLocalGlobalDofs.size());
}
}
@@ -236,7 +238,7 @@ namespace AMDiS {
values.clear();
// Global index of the current row dof.
- int globalRowDof = mapLocalGlobalDOFs[*cursor];
+ int globalRowDof = mapLocalGlobalDofs[*cursor];
// Test if the current row dof is a periodic dof.
bool periodicRow = (periodicDof.count(globalRowDof) > 0);
@@ -251,7 +253,7 @@ namespace AMDiS {
// Set only non null values.
if (value(*icursor) != 0.0) {
// Global index of the current column index.
- int globalColDof = mapLocalGlobalDOFs[col(*icursor)];
+ int globalColDof = mapLocalGlobalDofs[col(*icursor)];
// Calculate the exact position of the column index in the petsc matrix.
int colIndex = globalColDof * dispMult + dispAddCol;
@@ -337,7 +339,7 @@ namespace AMDiS {
DOFVector<double>::Iterator dofIt(vec, USED_DOFS);
for (dofIt.reset(); !dofIt.end(); ++dofIt) {
// Calculate global row index of the dof.
- int globalRow = mapLocalGlobalDOFs[dofIt.getDOFIndex()];
+ int globalRow = mapLocalGlobalDofs[dofIt.getDOFIndex()];
// Calculate petsc index of the row dof.
int index = globalRow * dispMult + dispAdd;
@@ -368,6 +370,8 @@ namespace AMDiS {
clock_t first = clock();
+ // === Create PETSc vector (rhs, solution and a temporary vector). ===
+
VecCreate(PETSC_COMM_WORLD, &petscRhsVec);
VecSetSizes(petscRhsVec, nRankRows, nOverallRows);
VecSetType(petscRhsVec, VECMPI);
@@ -381,19 +385,25 @@ namespace AMDiS {
VecSetType(petscTmpVec, VECMPI);
using mtl::tag::row; using mtl::tag::nz; using mtl::begin; using mtl::end;
- namespace traits= mtl::traits;
+ namespace traits = mtl::traits;
typedef DOFMatrix::base_matrix_type Matrix;
+ typedef std::vector<std::pair<int, int> > MatrixNnzEntry;
int d_nnz[nRankRows];
int o_nnz[nRankRows];
-
- std::map<int, std::vector< std::pair<int, int> > > sendMatrixEntry;
-
for (int i = 0; i < nRankRows; i++) {
d_nnz[i] = 0;
o_nnz[i] = 0;
}
+ MSG("Is rank DOF: %d\n", isRankDof[6717 - rstart * nComponents]);
+
+ // Stores to each rank a list of nnz entries (i.e. pairs of row and column index)
+ // that this rank will send to. This 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.
+ std::map<int, MatrixNnzEntry> sendMatrixEntry;
+
for (int i = 0; i < nComponents; i++) {
for (int j = 0; j < nComponents; j++) {
if ((*mat)[i][j]) {
@@ -408,41 +418,63 @@ namespace AMDiS {
for (cursor_type cursor = begin<row>(bmat),
cend = end<row>(bmat); cursor != cend; ++cursor) {
- int r = mapLocalGlobalDOFs[*cursor] * nComponents + i;
+ // Map the local row number to the global DOF index and create from it
+ // the global PETSc row index of this DOF.
+ int petscRowIdx = mapLocalGlobalDofs[*cursor] * nComponents + i;
if (isRankDof[*cursor]) {
- r -= rstart * nComponents;
+
+ // === 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 - rstart * nComponents;
#if (DEBUG != 0)
- if (r < 0 || r >= nRankRows) {
+ if (localPetscRowIdx < 0 || localPetscRowIdx >= nRankRows) {
std::cout << "ERROR in rank: " << mpiRank << std::endl;
- std::cout << " Wrong r = " << r << " " << *cursor << " "
- << mapLocalGlobalDOFs[*cursor] << " "
+ std::cout << " Wrong r = " << localPetscRowIdx << " " << *cursor
+ << " " << mapLocalGlobalDofs[*cursor] << " "
<< nRankRows << std::endl;
ERROR_EXIT("Should not happen!\n");
}
#endif
+ // Traverse all non zero entries in this row.
for (icursor_type icursor = begin<nz>(cursor),
icend = end<nz>(cursor); icursor != icend; ++icursor) {
if (value(*icursor) != 0.0) {
- int c = mapLocalGlobalDOFs[col(*icursor)] * nComponents + j;
+ int petscColIdx = mapLocalGlobalDofs[col(*icursor)] * nComponents + j;
- if (c >= rstart * nComponents &&
- c < rstart * nComponents + nRankRows)
- d_nnz[r]++;
+ // 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 >= rstart * nComponents &&
+ petscColIdx < rstart * nComponents + nRankRows)
+ d_nnz[localPetscRowIdx]++;
else
- o_nnz[r]++;
+ o_nnz[localPetscRowIdx]++;
}
}
} else {
- int sendToRank = -1;
+
+ // === 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. ===
+ // Find out who is the member of this DOF.
+ int sendToRank = -1;
for (RankToDofContainer::iterator it = recvDofs.begin();
it != recvDofs.end(); ++it) {
for (DofContainer::iterator dofIt = it->second.begin();
dofIt != it->second.end(); ++dofIt) {
if (**dofIt == *cursor) {
+
+ if (petscRowIdx == 6717) {
+ debug::writeDofMesh(mpiRank, *cursor, feSpace);
+ MSG("SEND DOF TO: %d/%d\n", it->first, *cursor);
+ }
+
sendToRank = it->first;
break;
}
@@ -454,117 +486,56 @@ namespace AMDiS {
TEST_EXIT_DBG(sendToRank != -1)("Should not happen!\n");
+ // Send all non zero entries to the member of the row DOF.
for (icursor_type icursor = begin<nz>(cursor),
icend = end<nz>(cursor); icursor != icend; ++icursor) {
if (value(*icursor) != 0.0) {
- int c = mapLocalGlobalDOFs[col(*icursor)] * nComponents + j;
+ int petscColIdx = mapLocalGlobalDofs[col(*icursor)] * nComponents + j;
- sendMatrixEntry[sendToRank].push_back(std::make_pair(r, c));
+ sendMatrixEntry[sendToRank].
+ push_back(std::make_pair(petscRowIdx, petscColIdx));
}
}
- }
- }
- }
- }
- }
-
- MPI::Request request[sendDofs.size() + recvDofs.size()];
- int requestCounter = 0;
-
- std::vector<int*> sendBuffers;
- sendBuffers.reserve(recvDofs.size());
-
- for (RankToDofContainer::iterator it = recvDofs.begin();
- it != recvDofs.end(); ++it) {
- int nSend = sendMatrixEntry[it->first].size();
-
- request[requestCounter++] = mpiComm.Isend(&nSend, 1, MPI_INT, it->first, 0);
-
- if (nSend > 0) {
- int *sendbuf = new int[nSend * 2];
- for (int i = 0; i < nSend; i++) {
- sendbuf[i * 2] = sendMatrixEntry[it->first][i].first;
- sendbuf[i * 2 + 1] = sendMatrixEntry[it->first][i].second;
- }
- sendBuffers.push_back(sendbuf);
- } else {
- sendBuffers.push_back(NULL);
- }
- }
-
- std::vector<int> recvSize;
- recvSize.reserve(sendDofs.size());
-
- int i = 0;
- for (RankToDofContainer::iterator it = sendDofs.begin();
- it != sendDofs.end(); ++it)
- request[requestCounter++] =
- mpiComm.Irecv(&(recvSize[i++]), 1, MPI_INT, it->first, 0);
-
- MPI::Request::Waitall(requestCounter, request);
-
- requestCounter = 0;
-
- i = 0;
- for (RankToDofContainer::iterator it = recvDofs.begin();
- it != recvDofs.end(); ++it) {
- int nSend = sendMatrixEntry[it->first].size();
-
- if (nSend > 0)
- request[requestCounter++] =
- mpiComm.Isend(sendBuffers[i], nSend * 2, MPI_INT, it->first, 0);
- i++;
+ } // if (isRankDof[*cursor]) ... else ...
+ } // for each row in mat[i][j]
+ } // if mat[i][j]
+ }
}
- std::vector<int*> recvBuffers;
- recvBuffers.reserve(sendDofs.size());
-
- i = 0;
- for (RankToDofContainer::iterator it = sendDofs.begin();
- it != sendDofs.end(); ++it) {
- if (recvSize[i] > 0) {
- int *recvbuf = new int[recvSize[i] * 2];
- recvBuffers.push_back(recvbuf);
+ // === Send and recv the nnz row structure to/from other ranks. ===
- request[requestCounter++] =
- mpiComm.Irecv(recvbuf, recvSize[i] * 2, MPI_INT, it->first, 0);
- } else {
- recvBuffers.push_back(NULL);
- }
-
- i++;
- }
-
- MPI::Request::Waitall(requestCounter, request);
+ StdMpi<MatrixNnzEntry, MatrixNnzEntry> stdMpi(mpiComm);
+ stdMpi.prepareCommunication(true);
+ stdMpi.send(sendMatrixEntry);
+ stdMpi.recv(sendDofs);
+ stdMpi.startCommunication<int>(MPI_INT);
- for (int j = 0; j < static_cast<int>(sendBuffers.size()); j++)
- if (sendBuffers[j])
- delete [] sendBuffers[j];
+ // === Evaluate the nnz structure this rank got from other ranks and add it to ===
+ // === the PETSc nnz data structure. ===
- i = 0;
- for (RankToDofContainer::iterator it = sendDofs.begin();
- it != sendDofs.end(); ++it) {
- if (recvSize[i] > 0) {
- for (int j = 0; j < recvSize[i]; j++) {
- int r = recvBuffers[i][j * 2];
- int c = recvBuffers[i][j * 2 + 1];
+ for (std::map<int, MatrixNnzEntry>::iterator it = stdMpi.getRecvData().begin();
+ it != stdMpi.getRecvData().end(); ++it) {
+ if (it->second.size() > 0) {
+ for (unsigned int i = 0; i < it->second.size(); i++) {
+ int r = it->second[i].first;
+ int c = it->second[i].second;
- r -= rstart * nComponents;
+ int localRowIdx = r - rstart * nComponents;
- TEST_EXIT_DBG(r >= 0 && r < nRankRows)("Should not happen!\n");
+ 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 < rstart * nComponents || c >= rstart * nComponents + nRankRows)
- o_nnz[r]++;
+ o_nnz[localRowIdx]++;
else
- d_nnz[r]++;
+ d_nnz[localRowIdx]++;
}
-
- delete [] recvBuffers[i];
}
+ }
- i++;
- }
+ // === Create PETSc matrix with the computed nnz data structure. ===
MatCreateMPIAIJ(PETSC_COMM_WORLD, nRankRows, nRankRows, nOverallRows, nOverallRows,
0, d_nnz, 0, o_nnz, &petscMatrix);
@@ -576,9 +547,7 @@ namespace AMDiS {
TEST_EXIT(b == rstart * nComponents + nRankRows)("Wrong matrix ownership range!\n");
#endif
- using mtl::tag::major; using mtl::tag::nz; using mtl::begin; using mtl::end;
- namespace traits= mtl::traits;
- typedef DOFMatrix::base_matrix_type Matrix;
+ // === Transfer values from DOF matrices to the PETSc matrix. ===
for (int i = 0; i < nComponents; i++)
for (int j = 0; j < nComponents; j++)
@@ -588,6 +557,8 @@ namespace AMDiS {
MatAssemblyBegin(petscMatrix, MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(petscMatrix, MAT_FINAL_ASSEMBLY);
+ // === Transfer values from DOF vector to the PETSc vector. ===
+
for (int i = 0; i < nComponents; i++)
setDofVector(petscRhsVec, vec->getDOFVector(i), nComponents, i);
@@ -789,7 +760,7 @@ namespace AMDiS {
#if (DEBUG != 0)
- writeMesh(feSpace);
+ debug::writeMesh(feSpace, -1, "mesh");
#endif
INFO(info, 8)("Parallel mesh adaption needed %.5f seconds\n",
@@ -798,6 +769,8 @@ namespace AMDiS {
// === Because the mesh has been changed, update the DOF numbering and mappings. ===
updateLocalGlobalNumbering();
+
+ exit(0);
}
@@ -1491,6 +1464,9 @@ namespace AMDiS {
for (DofContainer::iterator dofIt = it->second.begin();
dofIt != it->second.end(); ++dofIt) {
if (oldVertexDof[*dofIt]) {
+ if (mpiRank == 0 && **dofIt == 26) {
+ MSG("ADD RECV-1 26 TO RANK: %d\n", it->first);
+ }
recvDofs[it->first].push_back(*dofIt);
DofContainer::iterator eraseIt =
@@ -1548,9 +1524,13 @@ namespace AMDiS {
for (int i = 0; i < static_cast<int>(dofs.size()); i++) {
DofContainer::iterator eraseIt =
find(rankDofs.begin(), rankDofs.end(), dofs[i]);
- if (eraseIt != rankDofs.end())
- rankDofs.erase(eraseIt);
-
+ if (eraseIt != rankDofs.end())
+ rankDofs.erase(eraseIt);
+
+ if (mpiRank == 0 && *(dofs[i]) == 26) {
+ MSG("ADD RECV-2 26 TO RANK: %d\n", it->first);
+ }
+
dofsToRecv.push_back(dofs[i]);
}
}
@@ -1673,7 +1653,7 @@ namespace AMDiS {
DofIndexMap &rankOwnedDofsNewLocalIndex,
DofIndexMap &rankDofsNewGlobalIndex)
{
- mapLocalGlobalDOFs.clear();
+ mapLocalGlobalDofs.clear();
mapLocalToDofIndex.clear();
for (DofIndexMap::iterator dofIt = rankDofsNewLocalIndex.begin();
@@ -1682,7 +1662,7 @@ namespace AMDiS {
DegreeOfFreedom globalDof = rankDofsNewGlobalIndex[dofIt->first];
*const_cast<DegreeOfFreedom*>(dofIt->first) = localDof;
- mapLocalGlobalDOFs[localDof] = globalDof;
+ mapLocalGlobalDofs[localDof] = globalDof;
}
for (DofIndexMap::iterator dofIt = rankOwnedDofsNewLocalIndex.begin();
@@ -1821,7 +1801,7 @@ namespace AMDiS {
int *sendbuf = new int[dofs.size()];
for (int i = 0; i < static_cast<int>(dofs.size()); i++)
- sendbuf[i] = mapLocalGlobalDOFs[*(dofs[i])];
+ sendbuf[i] = mapLocalGlobalDofs[*(dofs[i])];
request[requestCounter++] =
mpiComm.Isend(sendbuf, dofs.size(), MPI_INT, it->first, 0);
@@ -1871,7 +1851,7 @@ namespace AMDiS {
// Added the received dofs to the mapping.
for (int j = 0; j < static_cast<int>(dofs.size()); j++) {
- int globalDofIndex = mapLocalGlobalDOFs[*(dofs[j])];
+ int globalDofIndex = mapLocalGlobalDofs[*(dofs[j])];
periodicDof[globalDofIndex].insert(recvBuffers[i][j]);
dofFromRank[globalDofIndex].insert(it->first);
}
@@ -2006,7 +1986,7 @@ namespace AMDiS {
serialize(out, sendDofs);
serialize(out, recvDofs);
- SerUtil::serialize(out, mapLocalGlobalDOFs);
+ SerUtil::serialize(out, mapLocalGlobalDofs);
SerUtil::serialize(out, mapLocalToDofIndex);
SerUtil::serialize(out, isRankDof);
@@ -2051,7 +2031,7 @@ namespace AMDiS {
deserialize(in, sendDofs, dofMap);
deserialize(in, recvDofs, dofMap);
- SerUtil::deserialize(in, mapLocalGlobalDOFs);
+ SerUtil::deserialize(in, mapLocalGlobalDofs);
SerUtil::deserialize(in, mapLocalToDofIndex);
SerUtil::deserialize(in, isRankDof);
@@ -2360,8 +2340,8 @@ namespace AMDiS {
if (rank == -1 || mpiRank == rank) {
std::cout << "====== DOF MAP LOCAL -> GLOBAL ====== " << std::endl;
- for (DofMapping::iterator it = mapLocalGlobalDOFs.begin();
- it != mapLocalGlobalDOFs.end(); it++) {
+ for (DofMapping::iterator it = mapLocalGlobalDofs.begin();
+ it != mapLocalGlobalDofs.end(); it++) {
DegreeOfFreedom localdof = -1;
if (mapLocalToDofIndex.count(it->first) > 0)
localdof = mapLocalToDofIndex[it->first];
@@ -2408,8 +2388,8 @@ namespace AMDiS {
std::cout << std::endl;
DegreeOfFreedom localdof = -1;
- for (DofMapping::iterator dofIt = mapLocalGlobalDOFs.begin();
- dofIt != mapLocalGlobalDOFs.end(); ++dofIt)
+ for (DofMapping::iterator dofIt = mapLocalGlobalDofs.begin();
+ dofIt != mapLocalGlobalDofs.end(); ++dofIt)
if (dofIt->second == it->first)
localdof = dofIt->first;
@@ -2468,114 +2448,4 @@ namespace AMDiS {
ElementFileWriter::writeFile(vec, feSpace, filename);
}
- void writeLocalElementDofs(int rank, int elIdx, FiniteElemSpace *feSpace)
- {
- using boost::lexical_cast;
-
- if (MPI::COMM_WORLD.Get_rank() == rank) {
- DOFVector<double> tmp(feSpace, "tmp");
- colorDofVectorByLocalElementDofs(tmp, feSpace->getMesh(), elIdx);
- VtkWriter::writeFile(tmp, "tmp" + lexical_cast<std::string>(elIdx) + ".vtu");
- }
- }
-
- void writeMesh(FiniteElemSpace *feSpace, int rank)
- {
- using boost::lexical_cast;
-
- int myRank = MPI::COMM_WORLD.Get_rank();
- if (rank == -1 || myRank == rank) {
- DOFVector<double> tmp(feSpace, "tmp");
- VtkWriter::writeFile(tmp, "mesh" + lexical_cast<std::string>(myRank) + ".vtu");
- }
- }
-
- void colorDofVectorByLocalElementDofs(DOFVector<double>& vec, Element *el)
- {
- // === Get local indices of the given element. ===
-
- const BasisFunction *basisFcts = vec.getFESpace()->getBasisFcts();
- int nBasisFcts = basisFcts->getNumber();
- std::vector<DegreeOfFreedom> localDofs(nBasisFcts);
- basisFcts->getLocalIndices(el, vec.getFESpace()->getAdmin(), localDofs);
-
- // === Set the values of the dof vector. ===
-
- vec.set(0.0);
- for (int i = 0; i < nBasisFcts; i++)
- vec[localDofs[i]] = static_cast<double>(i);
- }
-
- bool colorDofVectorByLocalElementDofs(DOFVector<double>& vec, Mesh *mesh, int elIndex)
- {
- FUNCNAME("colorDofVectorByLocalElementDofs()");
-
- TraverseStack stack;
- ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
- while (elInfo) {
- if (elInfo->getElement()->getIndex() == elIndex) {
- colorDofVectorByLocalElementDofs(vec, elInfo->getElement());
- return true;
- }
- elInfo = stack.traverseNext(elInfo);
- }
-
- return false;
- }
-
- Element* getDofIndexElement(FiniteElemSpace *feSpace, DegreeOfFreedom dof)
- {
- const BasisFunction* basFcts = feSpace->getBasisFcts();
- int nBasFcts = basFcts->getNumber();
- std::vector<DegreeOfFreedom> dofVec(nBasFcts);
-
- TraverseStack stack;
- ElInfo *elInfo = stack.traverseFirst(feSpace->getMesh(), -1,
- Mesh::CALL_EVERY_EL_PREORDER);
- while (elInfo) {
- basFcts->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(), dofVec);
- for (int i = 0; i < nBasFcts; i++)
- if (dofVec[i] == dof)
- return elInfo->getElement();
-
- elInfo = stack.traverseNext(elInfo);
- }
-
- return NULL;
- }
-
- Element* getLevel0ParentElement(Mesh *mesh, Element *el)
- {
- TraverseStack stack;
- ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
- while (elInfo) {
- if (elInfo->getElement() == el)
- return elInfo->getMacroElement()->getElement();
-
- elInfo = stack.traverseNext(elInfo);
- }
-
- return NULL;
- }
-
- void printInfoByDof(FiniteElemSpace *feSpace, DegreeOfFreedom dof)
- {
- Element *el = getDofIndexElement(feSpace, dof);
- Element *parEl = getLevel0ParentElement(feSpace->getMesh(), el);
-
- std::cout << "DOF-INFO: dof = " << dof
- << " elidx = " << el->getIndex()
- << " pelidx = " << parEl->getIndex() << std::endl;
-
- TraverseStack stack;
- ElInfo *elInfo = stack.traverseFirst(feSpace->getMesh(), -1,
- Mesh::CALL_EVERY_EL_PREORDER);
- while (elInfo) {
- if (elInfo->getElement()->getIndex() == parEl->getIndex())
- std::cout << "EL INFO TO " << parEl->getIndex() << ": "
- << elInfo->getType() << std::endl;
-
- elInfo = stack.traverseNext(elInfo);
- }
- }
}
diff --git a/AMDiS/src/ParallelDomainBase.h b/AMDiS/src/ParallelDomainBase.h
index 9a45327b610723bfd5f030c3fe409d5711480680..cdc9dceff388f30e072bb1edb74ea5801897dda5 100644
--- a/AMDiS/src/ParallelDomainBase.h
+++ b/AMDiS/src/ParallelDomainBase.h
@@ -178,8 +178,6 @@ namespace AMDiS {
return nRankDofs;
}
- void fillPetscMatrix(Matrix<DOFMatrix*> *mat, SystemVector *vec);
-
void solvePetscMatrix(SystemVector &vec);
virtual ProblemStatBase *getProblem(int number = 0)
@@ -218,7 +216,7 @@ namespace AMDiS {
/** \brief
* This function create new mappings from local to global indices,
- * \ref mapLocalGlobalDOFs, and from local to dof indices, \ref mapLocalToDofIndex.
+ * \ref mapLocalGlobalDofs, and from local to dof indices, \ref mapLocalToDofIndex.
* Furthermore, using the first argument the dof indices in ranks partition are
* changed.
*
@@ -239,11 +237,11 @@ namespace AMDiS {
* be used, if \ref partitionVec is set correctly. This is only the case, when
* the macro mesh is partitioned.
*
- * \param[out] partionDOFs Stores to each DOF pointer the set of ranks the DOF is
- * part of.
- * \param[out] rankDOFs Stores all rank DOFs.
- * \param[out] boundaryDOFs Stores all DOFs in ranks partition that are on an
- * interior boundary but correspond to another rank.
+ * \param[out] partitionDOFs Stores to each DOF pointer the set of ranks the DOF
+ * is part of.
+ * \param[out] rankDOFs Stores all rank DOFs.
+ * \param[out] boundaryDOFs Stores all DOFs in ranks partition that are on an
+ * interior boundary but correspond to another rank.
*/
void createDofMemberInfo(DofToPartitions& partitionDofs,
DofContainer& rankOwnedDofs,
@@ -251,6 +249,16 @@ namespace AMDiS {
DofToRank& boundaryDofs,
DofToBool& vertexDof);
+ /** \brief
+ * Create a PETSc matrix and PETSc vectors. The given DOF matrices are used to
+ * create the nnz structure of the PETSc matrix and the values are transfered to it.
+ * The given DOF vectors are used to the the values of the PETSc rhs vector.
+ *
+ * \param[in] mat
+ * \param[in] vec
+ */
+ void fillPetscMatrix(Matrix<DOFMatrix*> *mat, SystemVector *vec);
+
/// Takes a dof matrix and sends the values to the global petsc matrix.
void setDofMatrix(DOFMatrix* mat, int dispMult = 1,
int dispAddRow = 0, int dispAddCol = 0);
@@ -301,7 +309,7 @@ namespace AMDiS {
/** \brief
* This function is used for debugging only. It prints all information from
- * the local to global dof mapping, see \ref mapLocalGlobalDOFs.
+ * the local to global dof mapping, see \ref mapLocalGlobalDofs.
*
* \param rank If specified, only the information from the given rank is printed.
*/
@@ -537,13 +545,14 @@ namespace AMDiS {
RankToDofContainer sendDofs;
/** \brief
- * This map contains for each rank the list of dofs from which the current rank
- * must receive solution values of dofs at the interior boundaries.
+ * This map contains for each rank the list of DOFs from which the current rank
+ * will receive DOF values (i.e., this are all DOFs at an interior boundary). The
+ * DOF indices are given in rank's local numbering.
*/
RankToDofContainer recvDofs;
/// Maps local to global dof indices.
- DofMapping mapLocalGlobalDOFs;
+ DofMapping mapLocalGlobalDofs;
/// Maps local dof indices to real dof indices.
DofMapping mapLocalToDofIndex;
@@ -597,20 +606,6 @@ namespace AMDiS {
*/
long lastMeshChangeIndex;
};
-
- void writeLocalElementDofs(int rank, int elIdx, FiniteElemSpace *feSpace);
-
- void writeMesh(FiniteElemSpace *feSpace, int rank = -1);
-
- void colorDofVectorByLocalElementDofs(DOFVector<double>& vec, Element *el);
-
- bool colorDofVectorByLocalElementDofs(DOFVector<double>& vec, Mesh *mesh, int elIndex);
-
- Element* getDofIndexElement(FiniteElemSpace *feSpace, DegreeOfFreedom dof);
-
- Element* getLevel0ParentElement(Mesh *mesh, Element *el);
-
- void printInfoByDof(FiniteElemSpace *feSpace, DegreeOfFreedom dof);
}
#endif // AMDIS_PARALLELDOMAINBASE_H
diff --git a/AMDiS/src/StdMpi.h b/AMDiS/src/StdMpi.h
index c44fb5d00adcb1b99a9e3c0cee4d3315998ac0ce..fe17c99d36887ee3bdd6d5831c99b3a8152c7dd5 100644
--- a/AMDiS/src/StdMpi.h
+++ b/AMDiS/src/StdMpi.h
@@ -47,6 +47,11 @@ namespace AMDiS {
return data.size() * Global::getGeo(WORLD);
}
+ int intSizeOf(std::vector<std::pair<int, int> > &data)
+ {
+ return data.size() * 2;
+ }
+
void makeBuf(std::map<const DegreeOfFreedom*, DegreeOfFreedom> &data, int* buf)
{
int i = 0;
@@ -63,14 +68,13 @@ namespace AMDiS {
if (bufSize == 0)
return;
+ TEST_EXIT(bufSize % 2 == 0)("This should not happen!\n");
+
data.clear();
data.reserve(bufSize / 2);
- int i = 0;
- do {
- data.push_back(std::make_pair(buf[i], buf[i + 1]));
- i += 2;
- } while (i < bufSize);
+ for (int i = 0; i < (bufSize / 2); i++)
+ data.push_back(std::make_pair(buf[i * 2], buf[i * 2 + 1]));
}
void makeBuf(std::vector<MeshStructure> &data, unsigned long int *buf)
@@ -126,6 +130,14 @@ namespace AMDiS {
data[i][j] = buf[pos++];
}
+ void makeBuf(std::vector<std::pair<int, int> > &data, int *buf)
+ {
+ for (unsigned int i = 0; i < data.size(); i++) {
+ buf[i * 2] = data[i].first;
+ buf[i * 2 + 1] = data[i].second;
+ }
+ }
+
template<typename SendT, typename RecvT>
class StdMpi
{