diff --git a/AMDiS/src/DOFMatrix.cc b/AMDiS/src/DOFMatrix.cc
index 843b49f7f5ae1f2e9f0a9a58f47a611da6956475..08c978665826e923cff16310dae19101396d1fbe 100644
--- a/AMDiS/src/DOFMatrix.cc
+++ b/AMDiS/src/DOFMatrix.cc
@@ -14,7 +14,6 @@
#include "BoundaryCondition.h"
#include "BoundaryManager.h"
#include "Assembler.h"
-#include "Utilities.h"
namespace AMDiS {
diff --git a/AMDiS/src/ParallelDomainProblem.cc b/AMDiS/src/ParallelDomainProblem.cc
index b77bdd56a0e720335aaa642ddbe0982285e8b7bd..96bd9c810daed4d89576b4e1350ef67af59902a6 100644
--- a/AMDiS/src/ParallelDomainProblem.cc
+++ b/AMDiS/src/ParallelDomainProblem.cc
@@ -276,6 +276,10 @@ namespace AMDiS {
void ParallelDomainProblemBase::createInteriorBoundaryInfo(std::vector<const DegreeOfFreedom*>& rankDOFs,
std::map<const DegreeOfFreedom*, int>& boundaryDOFs)
{
+ FUNCNAME("ParallelDomainProblemBase::createInteriorBoundaryInfo()");
+
+ // === First, create all the information about the interior boundaries. ===
+
TraverseStack stack;
ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_NEIGH);
while (elInfo) {
@@ -323,20 +327,86 @@ namespace AMDiS {
/// === And add the part of the interior boundary. ===
AtomicBoundary& bound =
- interiorBoundary.getNewAtomicBoundary(ranksBoundary ? mpiRank :
- partitionVec[elInfo->getNeighbour(i)->getIndex()]);
+ (ranksBoundary ?
+ myIntBoundary.getNewAtomicBoundary(partitionVec[elInfo->getNeighbour(i)->getIndex()]) :
+ otherIntBoundary.getNewAtomicBoundary(partitionVec[elInfo->getNeighbour(i)->getIndex()]));
+
bound.rankObject.el = element;
bound.rankObject.subObjAtBoundary = EDGE;
bound.rankObject.ithObjAtBoundary = i;
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;
}
}
}
elInfo = stack.traverseNext(elInfo);
}
+
+
+ // === Once we have this information, we must care about their order. Eventually ===
+ // === all the boundaries have to be in the same order on both ranks (because ===
+ // === each boundary is shared by exactly two ranks). ===
+
+ std::vector<int*> sendBuffers;
+ std::vector<int*> recvBuffers;
+
+ for (std::map<int, std::vector<AtomicBoundary> >::iterator rankIt = myIntBoundary.boundary.begin();
+ rankIt != myIntBoundary.boundary.end();
+ ++rankIt) {
+ int* buffer = new int[rankIt->second.size()];
+ for (int i = 0; i < static_cast<int>(rankIt->second.size()); i++)
+ buffer[i] = (rankIt->second)[i].rankObject.el->getIndex();
+ sendBuffers.push_back(buffer);
+
+ mpiComm.Isend(buffer, rankIt->second.size(), MPI_INT, rankIt->first, 0);
+ }
+
+ for (std::map<int, std::vector<AtomicBoundary> >::iterator rankIt = otherIntBoundary.boundary.begin();
+ rankIt != otherIntBoundary.boundary.end();
+ ++rankIt) {
+ int *buffer = new int[rankIt->second.size()];
+ recvBuffers.push_back(buffer);
+
+ mpiComm.Irecv(buffer, rankIt->second.size(), MPI_INT, rankIt->first, 0);
+ }
+
+ mpiComm.Barrier();
+
+ int i = 0;
+ for (std::map<int, std::vector<AtomicBoundary> >::iterator rankIt = otherIntBoundary.boundary.begin();
+ rankIt != otherIntBoundary.boundary.end();
+ ++rankIt) {
+
+ // === We have received from rank "rankIt->first" the ordered list of element ===
+ // === indices. We now have to sort the corresponding list in this rank to ===
+ // === get the same order. ===
+
+ for (int j = 0; j < static_cast<int>(rankIt->second.size()); j++) {
+ // If the expected object is not at place, search for it.
+ if ((rankIt->second)[j].neighbourObject.el->getIndex() != recvBuffers[i][j]) {
+ int k = j + 1;
+ for (; k < static_cast<int>(rankIt->second.size()); k++)
+ if ((rankIt->second)[k].neighbourObject.el->getIndex() == recvBuffers[i][j])
+ 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");
+
+ // Swap the current with the found element.
+ AtomicBoundary tmpBound = (rankIt->second)[k];
+ (rankIt->second)[k] = (rankIt->second)[j];
+ (rankIt->second)[j] = tmpBound;
+ }
+ }
+
+ delete [] recvBuffers[i++];
+ }
+
+ for (int i = 0; i < static_cast<int>(sendBuffers.size()); i++)
+ delete [] sendBuffers[i];
}
@@ -527,8 +597,10 @@ namespace AMDiS {
void ParallelDomainProblemBase::updateLocalGlobalNumbering(int& nRankDOFs,
int& nOverallDOFs)
{
+ FUNCNAME("ParallelDomainProblemBase::updateLocalGlobalNumbering()");
+
std::set<const DegreeOfFreedom*> rankDOFs;
- std::set<const DegreeOfFreedom*> boundaryDOFs;
+ std::map<const DegreeOfFreedom*, int> boundaryDOFs;
/// === Get all DOFs in ranks partition. ===
@@ -545,20 +617,105 @@ namespace AMDiS {
}
- // === Traverse on interior boundaries and move all not ranked owned DOFs from
- // rankDOFs to boundaryDOFs === //
+ // === Traverse on interior boundaries and move all not ranked owned DOFs from ===
+ // === rankDOFs to boundaryDOFs ===
for (std::map<int, std::vector<AtomicBoundary> >::iterator it =
- interiorBoundary.boundary.begin();
- it != interiorBoundary.boundary.end();
+ myIntBoundary.boundary.begin();
+ it != myIntBoundary.boundary.end();
++it) {
- for (int i = 0; i < it->second.size(); i++) {
-
+ 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");
+ }
+
+ std::vector<const DegreeOfFreedom*> boundDOFs;
+ addAllDOFs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary,
+ boundDOFs);
}
}
}
+ void ParallelDomainProblemBase::addAllDOFs(Element *el, int ithEdge,
+ std::vector<const DegreeOfFreedom*>& dofs,
+ bool addVertices)
+ {
+ FUNCNAME("ParallelDomainProblemBase::addAllDOFs()");
+
+ 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");
+ }
+
+ dofs.push_back(boundDOF1);
+ }
+
+ switch (ithEdge) {
+ case 0:
+ if (el->getSecondChild() && el->getSecondChild()->getFirstChild()) {
+ addAllDOFs(el->getSecondChild()->getFirstChild(), 0, dofs, false);
+ dofs.push_back(el->getSecondChild()->getFirstChild()->getDOF(2));
+ addAllDOFs(el->getSecondChild()->getSecondChild(), 1, dofs, false);
+ }
+ break;
+ case 1:
+ if (el->getFirstChild() && el->getFirstChild()->getFirstChild()) {
+ addAllDOFs(el->getFirstChild()->getFirstChild(), 0, dofs, false);
+ dofs.push_back(el->getFirstChild()->getFirstChild()->getDOF(2));
+ addAllDOFs(el->getFirstChild()->getSecondChild(), 1, dofs, false);
+ }
+ break;
+ case 2:
+ if (el->getFirstChild()) {
+ addAllDOFs(el->getFirstChild(), 0, dofs, false);
+ dofs.push_back(el->getFirstChild()->getDOF(2));
+ addAllDOFs(el->getSecondChild(), 1, dofs, false);
+ }
+ break;
+ default:
+ ERROR_EXIT("Should never happen!\n");
+ }
+
+ if (addVertices)
+ dofs.push_back(boundDOF2);
+ }
+
+
void ParallelDomainProblemBase::createDOFMemberInfo(
std::map<const DegreeOfFreedom*, std::set<int> >& partitionDOFs,
std::vector<const DegreeOfFreedom*>& rankDOFs,
diff --git a/AMDiS/src/ParallelDomainProblem.h b/AMDiS/src/ParallelDomainProblem.h
index ee2b9c24723cc9a4e201c6094ab580a6ad6273bf..5fec46acccb9caddbbe2c5500933c767f2ce47e9 100644
--- a/AMDiS/src/ParallelDomainProblem.h
+++ b/AMDiS/src/ParallelDomainProblem.h
@@ -137,9 +137,15 @@ namespace AMDiS {
void updateLocalGlobalNumbering(int& nRankDOFs, int& nOverallDOFs);
+ void addAllDOFs(Element *el, int ithEdge,
+ std::vector<const DegreeOfFreedom*>& dofs,
+ bool addVertices = true);
+
/** \brief
* This function traverses the whole mesh, i.e. before it is really partitioned,
- * and collects information about which DOF corresponds to which rank.
+ * and collects information about which DOF corresponds to which rank. Can only
+ * 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.
@@ -215,10 +221,20 @@ namespace AMDiS {
int nRankDOFs;
/** \brief
- * Defines the interioir boundaries of the domain that result from partitioning
- * the whole mesh.
+ * Defines the interior boundaries of the domain that result from partitioning
+ * the whole mesh. Contains only the boundaries, which are owned by the rank, i.e.,
+ * the object gives for every neighbour rank i the boundaries this rank owns and
+ * shares with rank i.
+ */
+ InteriorBoundary myIntBoundary;
+
+ /** \brief
+ * Defines the interior boundaries of the domain that result from partitioning
+ * the whole mesh. Contains only the boundaries, which are not owned by the rank,
+ * i.e., the object gives for every neighbour rank i the boundaries that are
+ * owned by rank i and are shared with this rank.
*/
- InteriorBoundary interiorBoundary;
+ InteriorBoundary otherIntBoundary;
/** \brief
* This map contains for each rank the list of dofs the current rank must send
diff --git a/AMDiS/src/UmfPackSolver.h b/AMDiS/src/UmfPackSolver.h
index 62a88f09a82e0cc255bb24fbb28b4ee8b9603a40..6add42d9036f3f9c3c56fe44c8fe6dfbb3f6a3fe 100644
--- a/AMDiS/src/UmfPackSolver.h
+++ b/AMDiS/src/UmfPackSolver.h
@@ -51,7 +51,8 @@ namespace AMDiS {
/** \brief
* Returns a new UmfPackSolver object.
*/
- OEMSolver* create() {
+ OEMSolver* create()
+ {
return new UmfPackSolver(this->name);
}
};
@@ -67,7 +68,11 @@ namespace AMDiS {
}
/// Destructor
- ~UmfPackSolver() { if (solver) delete solver;}
+ ~UmfPackSolver()
+ {
+ if (solver)
+ delete solver;
+ }
/// Solves the system directly
int solveSystem(const DOFMatrix::base_matrix_type& A,
diff --git a/AMDiS/src/Utilities.h b/AMDiS/src/Utilities.h
deleted file mode 100644
index 143704ebc7859f6304a584125c78bd2b51e0cff9..0000000000000000000000000000000000000000
--- a/AMDiS/src/Utilities.h
+++ /dev/null
@@ -1,44 +0,0 @@
-// ============================================================================
-// == ==
-// == AMDiS - Adaptive multidimensional simulations ==
-// == ==
-// ============================================================================
-// == ==
-// == crystal growth group ==
-// == ==
-// == Stiftung caesar ==
-// == Ludwig-Erhard-Allee 2 ==
-// == 53175 Bonn ==
-// == germany ==
-// == ==
-// ============================================================================
-// == ==
-// == http://www.caesar.de/cg/AMDiS ==
-// == ==
-// ============================================================================
-
-/** \file Utilities.h */
-
-#ifndef AMDIS_UTILITIES_H
-#define AMDIS_UTILITIES_H
-
-#include <vector>
-#include <algorithm>
-
-
-
-template <typename T>
-std::vector<T> inline invert_reorder(const std::vector<T>& v)
-{
- T my_max= *std::max_element(v.begin(), v.end()) + 1;
- std::vector<T> u(my_max, T(-1));
- for (int i= 0; i < v.size(); i++)
- if (v[i] >= 0) {
- assert(u[v[i]] == -1); // check double insertion
- u[v[i]]= T(i);
- }
- assert(find(u.begin(), u.end(), T(-1)) == u.end()); // check if all entries are set
- return u;
-}
-
-#endif // AMDIS_UTILITIES_H
diff --git a/AMDiS/src/ValueWriter.h b/AMDiS/src/ValueWriter.h
index eb77eef87c56d794cebfc6d5db2ead8646a50206..6f6b9d66f2d6c115a71e09b3103c3309c6140834 100644
--- a/AMDiS/src/ValueWriter.h
+++ b/AMDiS/src/ValueWriter.h
@@ -29,18 +29,10 @@
#include "Flag.h"
#include "Mesh.h"
#include "DataCollector.h"
+#include "AMDiS_fwd.h"
namespace AMDiS {
- class DOFAdmin;
- class ElInfo;
-
- template<typename T> class DOFVector;
-
- // ============================================================================
- // ===== class ValueWriter ====================================================
- // ============================================================================
-
/** \ingroup Output
* \brief
* ValueWriter is a static class which writes the values of a DOFVector
@@ -53,9 +45,7 @@ namespace AMDiS {
class ValueWriter
{
public:
- /** \brief
- * Writes DOFVector values to values->feSpace->mesh.
- */
+ /// Writes DOFVector values to values->feSpace->mesh.
static void writeValues(DataCollector *dc,
const std::string filename,
double time = 0.0,
@@ -64,34 +54,22 @@ namespace AMDiS {
bool (*writeElem)(ElInfo*) = NULL);
protected:
- /** \brief
- * File to which the values should be written
- */
+ /// File to which the values should be written
static FILE *valueFile;
- /** \brief
- * Global interpolation point indexing
- */
+ /// Global interpolation point indexing
static DOFVector<int> *interpPointInd;
- /** \brief
- * list of coords for each dof
- */
+ /// list of coords for each dof
static DOFVector< std::list<WorldVector<double> > > *dofCoords;
- /** \brief
- * DOFAdmin of values
- */
+ /// DOFAdmin of values
static const DOFAdmin *admin;
- /** \brief
- * Pointer to have a global access to values
- */
+ /// Pointer to have a global access to values
static DOFVector<double> *valueVec;
- /** \brief
- * Total number of interpolation points.
- */
+ /// Total number of interpolation points.
static int ni;
};
diff --git a/AMDiS/src/VertexInfo.h b/AMDiS/src/VertexInfo.h
index 75f1ae756d57751e16b4947d9352a72fdb8c6c80..5337385235d74728012626ff8a89a81e44445fd9 100644
--- a/AMDiS/src/VertexInfo.h
+++ b/AMDiS/src/VertexInfo.h
@@ -26,35 +26,27 @@
namespace AMDiS {
- /** \brief
- * Stores coordinates and output index for one vertex.
- */
+ /// Stores coordinates and output index for one vertex.
class VertexInfo
{
public:
- /** \brief
- * Coordinates for this vertex.
- */
+ /// Coordinates for this vertex.
WorldVector<double> coords;
- /** \brief
- * Index for the output file.
- */
+ /// Index for the output file.
int outputIndex;
- /** \brief
- * Used to check, whether coords are already stored for a given dof.
- */
- bool operator==(const WorldVector<double>& c) {
+ /// Used to check, whether coords are already stored for a given dof.
+ bool operator==(const WorldVector<double>& c)
+ {
return (c == coords);
- };
+ }
- /** \brief
- * Used to check, whether coords are already stored for a given dof.
- */
- bool operator!=(const WorldVector<double>& c) {
+ /// Used to check, whether coords are already stored for a given dof.
+ bool operator!=(const WorldVector<double>& c)
+ {
return (c != coords);
- };
+ }
};
}
diff --git a/AMDiS/src/VertexVector.h b/AMDiS/src/VertexVector.h
index bd62d53ce6a5c69675d2fa981cb440f416adb5a5..6601f6920bfaab76faabe8698a9400395665886a 100644
--- a/AMDiS/src/VertexVector.h
+++ b/AMDiS/src/VertexVector.h
@@ -21,24 +21,27 @@ namespace AMDiS {
~VertexVector();
- const DOFAdmin *getAdmin() { return admin; };
+ const DOFAdmin *getAdmin()
+ {
+ return admin;
+ }
- void resize(int size) {
- int i, oldSize = static_cast<int>(vec.size());
+ void resize(int size)
+ {
+ int oldSize = static_cast<int>(vec.size());
vec.resize(size);
- for(i = oldSize; i < size; i++) {
+ for (int i = oldSize; i < size; i++)
vec[i] = i;
- }
}
void set(DegreeOfFreedom val);
- void compressDOFContainer(int size, std::vector<DegreeOfFreedom> &newDOF) {
+ void compressDOFContainer(int size, std::vector<DegreeOfFreedom> &newDOF)
+ {
DOFContainer::compressDOFContainer(size, newDOF);
int totalSize = getAdmin()->getSize();
- for(int i = size; i < totalSize; i++) {
+ for(int i = size; i < totalSize; i++)
vec[i] = i;
- }
}
protected: