diff --git a/AMDiS/src/CoarseningManager.cc b/AMDiS/src/CoarseningManager.cc index 69fbe11f3b323576810944799dada002b58eb573..ed5dd54e99a1c2d04875ae0d9fc8714260843fce 100644 --- a/AMDiS/src/CoarseningManager.cc +++ b/AMDiS/src/CoarseningManager.cc @@ -89,12 +89,8 @@ namespace AMDiS { Flag CoarseningManager::coarsenMesh(Mesh *aMesh) { - int n_elements; - ElInfo *el_info; - mesh = aMesh; - - n_elements = mesh->getNumberOfLeaves(); + int n_elements = mesh->getNumberOfLeaves(); spreadCoarsenMark(); @@ -102,14 +98,10 @@ namespace AMDiS { do { doMore = false; - el_info = stack->traverseFirst(mesh, -1, - Mesh::CALL_EVERY_EL_POSTORDER | Mesh::FILL_NEIGH); + ElInfo* el_info = + stack->traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_POSTORDER | Mesh::FILL_NEIGH); while (el_info) { - int idx = el_info->getElement()->getIndex(); - - // if (idx != 2288 && idx != 2283) - coarsenFunction(el_info); - + coarsenFunction(el_info); el_info = stack->traverseNext(el_info); } } while (doMore); diff --git a/AMDiS/src/DOFMatrix.cc b/AMDiS/src/DOFMatrix.cc index 8b71e9a9e13ee1f3f3ae2af9c5cf830fa7051aa1..843b49f7f5ae1f2e9f0a9a58f47a611da6956475 100644 --- a/AMDiS/src/DOFMatrix.cc +++ b/AMDiS/src/DOFMatrix.cc @@ -57,10 +57,6 @@ namespace AMDiS { colIndices.resize(nCol); applyDBCs.clear(); - -#ifdef HAVE_PARALLEL_DOMAIN_AMDIS - applicationOrdering = NULL; -#endif } DOFMatrix::DOFMatrix(const DOFMatrix& rhs) @@ -213,17 +209,23 @@ namespace AMDiS { bound ? boundaryManager->getBoundaryCondition(bound[i]) : NULL; if (condition && condition->isDirichlet()) { - if (condition->applyBoundaryCondition()) + if (condition->applyBoundaryCondition()) { +#ifdef HAVE_PARALLEL_DOMAIN_AMDIS + if (isRankDOF[row]) + applyDBCs.insert(static_cast<int>(row)); +#else applyDBCs.insert(static_cast<int>(row)); +#endif + } } else for (int j = 0; j < nCol; j++) { // for all columns DegreeOfFreedom col = colIndices[j]; double entry = elMat[i][j]; if (add) - ins[row][col]+= sign * entry; + ins[row][col] += sign * entry; else - ins[row][col]= sign * entry; + ins[row][col] = sign * entry; } } } @@ -249,6 +251,18 @@ 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, diff --git a/AMDiS/src/DOFMatrix.h b/AMDiS/src/DOFMatrix.h index cc680eab408d8fd3e4c824d248b6ce98d67bb7b2..a2a114281ff7034f008975b958b75700cdb464c3 100644 --- a/AMDiS/src/DOFMatrix.h +++ b/AMDiS/src/DOFMatrix.h @@ -37,10 +37,6 @@ #include "Boundary.h" #include "Serializable.h" -#ifdef HAVE_PARALLEL_DOMAIN_AMDIS -#include "petscao.h" -#endif - namespace AMDiS { /** \ingroup DOFAdministration @@ -447,10 +443,9 @@ namespace AMDiS { int memsize(); #ifdef HAVE_PARALLEL_DOMAIN_AMDIS - /// Sets the petsc application ordering object to map dof indices. - void useApplicationOrdering(AO *ao) + void setIsRankDOF(std::map<DegreeOfFreedom, bool>& dofmap) { - applicationOrdering = ao; + isRankDOF = dofmap; } #endif @@ -525,8 +520,7 @@ namespace AMDiS { std::set<int> applyDBCs; #ifdef HAVE_PARALLEL_DOMAIN_AMDIS - /// Petsc application ordering to map dof indices. - AO *applicationOrdering; + std::map<DegreeOfFreedom, bool> isRankDOF; #endif /// Inserter object: implemented as pointer, allocated and deallocated as needed diff --git a/AMDiS/src/FileWriter.h b/AMDiS/src/FileWriter.h index 228bf1715f942eb44992e93e4374e773231b63b4..17c0fd4f39ef709305329c0d44f4606175ed5f53 100644 --- a/AMDiS/src/FileWriter.h +++ b/AMDiS/src/FileWriter.h @@ -82,11 +82,13 @@ namespace AMDiS { writeElement = writeElem; } - std::string getFilename() { + std::string getFilename() + { return filename; } - void setFilename(std::string n) { + void setFilename(std::string n) + { filename = n; } diff --git a/AMDiS/src/FirstOrderAssembler.h b/AMDiS/src/FirstOrderAssembler.h index 723ee31e6816af5478d8065ee67296209567499e..d6e06b6d9336ffcecef1978a01e7dea9c4b9db7d 100644 --- a/AMDiS/src/FirstOrderAssembler.h +++ b/AMDiS/src/FirstOrderAssembler.h @@ -142,7 +142,8 @@ namespace AMDiS { void calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat); /// Implements SubAssembler::calculateElementVector(). - void calculateElementVector(const ElInfo*, ElementVector&) { + void calculateElementVector(const ElInfo*, ElementVector&) + { ERROR_EXIT("should not be called\n"); } }; @@ -193,7 +194,8 @@ namespace AMDiS { void calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat); /// Implements SubAssembler::calculateElementVector(). - void calculateElementVector(const ElInfo*, ElementVector&) { + void calculateElementVector(const ElInfo*, ElementVector&) + { ERROR_EXIT("should not be called\n"); } diff --git a/AMDiS/src/FixVec.h b/AMDiS/src/FixVec.h index cffc0ce9ac14c8a4040a5fe9c03a3af37ddc29d1..62dc71023a46fba84092dfbe87673a1f5eb27057 100644 --- a/AMDiS/src/FixVec.h +++ b/AMDiS/src/FixVec.h @@ -89,28 +89,31 @@ namespace AMDiS { } /// Initialisation for dim. - inline void init(int dim) { + inline void init(int dim) + { this->resize(calcSize(dim)); } /// Initialisation for size - inline void initSize(int size) { + inline void initSize(int size) + { this->resize(size); } /// Returns the \ref size_ of the FixVec. - inline int size() const { + inline int size() const + { return this->getSize(); } protected: /// Determines needed vector size. - static int calcSize(int dim) { - if (dim < 0) { + static int calcSize(int dim) + { + if (dim < 0) return Global::getGeo(WORLD); - } else { + else return Global::getGeo(d, dim); - } } public: diff --git a/AMDiS/src/FixVecConvert.h b/AMDiS/src/FixVecConvert.h index 358de8979e1b0427d9a85e82d5e115285a7105e6..470d5de7989adfb3ce0bde78e1704db5f0570ff4 100644 --- a/AMDiS/src/FixVecConvert.h +++ b/AMDiS/src/FixVecConvert.h @@ -30,7 +30,8 @@ namespace AMDiS { class VecConv { public: - static FixVec<T,d1>& convertVec(FixVec<T,d2>& rhs, Mesh* mesh) { + static FixVec<T,d1>& convertVec(FixVec<T,d2>& rhs, Mesh* mesh) + { TEST_EXIT(mesh->getGeo(d1)==mesh->getGeo(d2))("Incompatible dimensions.\n"); return reinterpret_cast<FixVec<T,d1>&>(rhs); } diff --git a/AMDiS/src/Flag.h b/AMDiS/src/Flag.h index 220625a76213c1359ef37a04c437f29a4423b224..15524ffc899323fdb97aad3a2c64d9842c29fabd 100644 --- a/AMDiS/src/Flag.h +++ b/AMDiS/src/Flag.h @@ -48,129 +48,150 @@ namespace AMDiS { inline ~Flag() {} /// Compares two Flags - inline bool operator==(const Flag& f) const { + inline bool operator==(const Flag& f) const + { return (flags == f.flags); } /// Compares two Flags - bool operator!=(const Flag& f) const { + bool operator!=(const Flag& f) const + { return !(f == *this); } /// Assignment operator - inline Flag& operator=(const Flag& f) { + inline Flag& operator=(const Flag& f) + { if (this != &f) flags = f.flags; return *this; } /// Typecast - inline operator bool() const { + inline operator bool() const + { return isAnySet(); } /// Set \ref flags - inline void setFlags(const unsigned long f) { + inline void setFlags(const unsigned long f) + { flags = f; } /// Set \ref flags - inline void setFlags(const Flag& f) { + inline void setFlags(const Flag& f) + { flags = f.flags; } /// Sets \ref flags to \ref flags | f - inline void setFlag(const unsigned long f) { + inline void setFlag(const unsigned long f) + { flags |= f; } /// Sets \ref flags to \ref flags | f.flags - inline void setFlag(const Flag& f) { + inline void setFlag(const Flag& f) + { flags |= f.flags; } /// Sets \ref flags to \ref flags & ~f - inline void unsetFlag(const unsigned long f) { + inline void unsetFlag(const unsigned long f) + { flags &= ~f; } /// Sets \ref flags to \ref flags & ~f.flags - inline void unsetFlag(const Flag& f) { + inline void unsetFlag(const Flag& f) + { flags &= ~f.flags; } - inline const unsigned long getFlags() const { + inline const unsigned long getFlags() const + { return flags; } /// Returns \ref flags | f.flags - inline Flag operator+(const Flag& f) const { + inline Flag operator+(const Flag& f) const + { Flag r(flags); r.setFlag(f); return r; } /// Returns \ref flags & ~f.flags - inline Flag operator-(const Flag& f) const { + inline Flag operator-(const Flag& f) const + { Flag r(flags); r.unsetFlag(f); return r; } /// Returns \ref flags | f.flags - inline Flag operator|(const Flag& f) const { + inline Flag operator|(const Flag& f) const + { Flag r(flags); r.setFlag(f); return r; } /// Returns \ref flags & f.flags - inline Flag operator&(const Flag& f) const { + inline Flag operator&(const Flag& f) const + { Flag r(flags); r.flags &= f.flags; return r; } /// Sets \ref flags to \ref flags &= f.flags - inline Flag operator&=(const Flag& f) { + inline Flag operator&=(const Flag& f) + { flags &= f.flags; return *this; } /// Returns \ref flags ^ f.flags - inline Flag operator^(const Flag& f) const { + inline Flag operator^(const Flag& f) const + { Flag r(flags); r.flags ^=f.flags; return r; } /// Sets \ref flags to \ref flags & f.flags - inline Flag& operator|=(const Flag& f) { - if (this != &f) { + inline Flag& operator|=(const Flag& f) + { + if (this != &f) flags |= f.flags; - }; return *this; } /// Returns ~\ref flags - inline Flag operator~() const { + inline Flag operator~() const + { Flag r; r.flags = ~flags; return r; } /// Checks whether all set bits of f.flags are set in \ref flags too. - inline bool isSet(const Flag& f) const { + inline bool isSet(const Flag& f) const + { return ((flags&f.flags) == f.flags); } /// Returns !\ref isSet(f) - inline bool isUnset(const Flag& f) const { + inline bool isUnset(const Flag& f) const + { return !isSet(f); } /// Returns true if \ref flags != 0 - inline bool isAnySet() const { + inline bool isAnySet() const + { return (flags != 0); } diff --git a/AMDiS/src/GNUPlotWriter.h b/AMDiS/src/GNUPlotWriter.h index b07c9104111fe7da0d6ed49337b3e8686a27d983..bd4b9cc6434caa5643ffbd22f1c438b165ccc1a7 100644 --- a/AMDiS/src/GNUPlotWriter.h +++ b/AMDiS/src/GNUPlotWriter.h @@ -25,49 +25,33 @@ #include <vector> #include <string> #include "FileWriter.h" +#include "AMDiS_fwd.h" namespace AMDiS { - template<typename T> class DOFVector; - class FiniteElemSpace; - - /** \brief - * - */ + /// class GNUPlotWriter : public FileWriterInterface { public: - /** \brief - * Constructor - */ + /// GNUPlotWriter(std::string filename, FiniteElemSpace *feSpace, std::vector<DOFVector<double>*> &dofVectors); - /** \brief - * Destructor - */ - virtual ~GNUPlotWriter() {}; + /// + virtual ~GNUPlotWriter() {} - /** \brief - * Implementation of FileWriter::writeFiles() - */ + /// virtual void writeFiles(AdaptInfo *adaptInfo, bool force); protected: - /** \brief - * Contains the mesh - */ + /// Contains the mesh FiniteElemSpace *feSpace_; - /** \brief - * vector of dof vectors to write - */ + /// vector of dof vectors to write std::vector<DOFVector<double>*> dofVectors_; - /** \brief - * file name - */ + /// file name std::string filename_; }; diff --git a/AMDiS/src/Global.h b/AMDiS/src/Global.h index b6778ff3524d8535ffb2c8997977274bfa616743..6256f13a595cb04371daa8efe0e02b365961c68a 100644 --- a/AMDiS/src/Global.h +++ b/AMDiS/src/Global.h @@ -187,25 +187,43 @@ namespace AMDiS { int ntime); /// Sets \ref msgInfo - static void setMsgInfo(int info) { msgInfo = info; } + static void setMsgInfo(int info) + { + msgInfo = info; + } /// Returns \ref msgInfo - static int getMsgInfo() { return msgInfo; } + static int getMsgInfo() + { + return msgInfo; + } /// Sets \ref msgWait - static void setMsgWait(bool wait) { msgWait = wait; } + static void setMsgWait(bool wait) + { + msgWait = wait; + } /// Returns \ref msgWait - static bool getMsgWait() { return msgWait; } + static bool getMsgWait() + { + return msgWait; + } /// Waits for enter if w is true static void wait(bool w); /// Returns \ref out - static std::ostream *getOutStream() { return out; } + static std::ostream *getOutStream() + { + return out; + } /// Returns \ref error - static std::ostream *getErrorStream() { return error; } + static std::ostream *getErrorStream() + { + return error; + } protected: /// Message stram @@ -341,14 +359,16 @@ namespace AMDiS { * returns a pointer to \ref referenceElement [dim]. With this pointer you * can get information about the element via Element's getGeo method. */ - static const Element *getReferenceElement(int dim) { + static const Element *getReferenceElement(int dim) + { FUNCNAME("Global::getReferenceElement()"); TEST_EXIT((dim > 0) && (dim < 4))("invalid dim: %d\n", dim); return referenceElement[dim]; } /// returns geometrical information. Currently this is only dimOfWorld. - static inline int getGeo(GeoIndex p) { + static inline int getGeo(GeoIndex p) + { if (WORLD == p) return dimOfWorld; @@ -360,7 +380,8 @@ namespace AMDiS { * returns geometrical information about elements of the dimension dim. * getGeo(VERTEX, 3) returns 4 because a Tetrahedron has 4 vertices. */ - static inline int getGeo(GeoIndex p, int dim) { + static inline int getGeo(GeoIndex p, int dim) + { TEST_EXIT_DBG((p >= MINPART) && (p <= MAXPART)) ("Calling for invalid geometry value %d\n",p); TEST_EXIT_DBG((dim >= 0) && (dim < 4)) diff --git a/AMDiS/src/GridWriter.h b/AMDiS/src/GridWriter.h index 1d5f8e94597f783d9a61ddf420a5b6582d6ea0ec..0cce47b9d6c56746f5f94f0f0fbfb5d5bb34b691 100644 --- a/AMDiS/src/GridWriter.h +++ b/AMDiS/src/GridWriter.h @@ -23,16 +23,10 @@ #define AMDIS_GRID_WRITER_H #include "Global.h" +#include "AMDiS_fwd.h" namespace AMDiS { - template<typename T> class WorldVector; - template<typename T> class DOFVector; - - // ============================================================================ - // ===== class GridWriter ===================================================== - // ============================================================================ - /** \ingroup Output * \brief * Produces a grid-based output of a dof-vector @@ -62,9 +56,7 @@ namespace AMDiS { int *numPoints, double *dist, DOFVector<T> *vec, - char* filename); - - private: + char* filename); }; } diff --git a/AMDiS/src/Lagrange.h b/AMDiS/src/Lagrange.h index 97621586cc266ce55d8a08960187de1531a8fe4c..b0f23c249ac65e2a81444f4254240377880cb4b1 100644 --- a/AMDiS/src/Lagrange.h +++ b/AMDiS/src/Lagrange.h @@ -61,38 +61,27 @@ namespace AMDiS { */ static Lagrange* getLagrange(int dim, int degree); - /** \brief - * Implements BasisFunction::interpol - */ + /// Implements BasisFunction::interpol const double *interpol(const ElInfo *, int, const int *, AbstractFunction<double, WorldVector<double> >*, double *); - - /** \brief - * Implements BasisFunction::interpol - */ + /// Implements BasisFunction::interpol const WorldVector<double> *interpol(const ElInfo *, int, const int *b_no, AbstractFunction<WorldVector<double>, WorldVector<double> >*, WorldVector<double> *); - /** \brief - * Returns the barycentric coordinates of the i-th basis function. - */ + /// Returns the barycentric coordinates of the i-th basis function. DimVec<double> *getCoords(int i) const; - /** \brief - * Implements BasisFunction::getDOFIndices - */ + /// Implements BasisFunction::getDOFIndices const DegreeOfFreedom* getDOFIndices(const Element*, const DOFAdmin&, DegreeOfFreedom*) const; - /** \brief - * Implements BasisFunction::getBound - */ + /// Implements BasisFunction::getBound void getBound(const ElInfo*, BoundaryType *) const; /** \brief @@ -158,44 +147,29 @@ namespace AMDiS { */ void setBary(); - /** \brief - * Recursive calculation of coordinates. Used by \ref setBary - */ + /// Recursive calculation of coordinates. Used by \ref setBary void createCoords(int* coordInd, int numCoords, int dimIndex, int rest, DimVec<double>* vec = NULL); - /** \brief - * Used by \ref setBary - */ + /// Used by \ref setBary int** getIndexPermutations(int numIndices) const; - /** \brief - * Implements BasisFunction::setNDOF - */ + /// Implements BasisFunction::setNDOF void setNDOF(); - /** \brief - * Sets used function pointers - */ + /// Sets used function pointers void setFunctionPointer(); - /** \brief - * Used by \ref getDOFIndices and \ref getVec - */ + /// Used by \ref getDOFIndices and \ref getVec int* orderOfPositionIndices(const Element* el, GeoIndex position, int positionIndex) const; - /** \brief - * Calculates the number of DOFs needed for Lagrange of the given dim and - * degree. - */ + /// Calculates the number of DOFs needed for Lagrange of the given dim and degree. static int getNumberOfDOFs(int dim, int degree); private: - /** \brief - * barycentric coordinates of the locations of all basis functions - */ + /// barycentric coordinates of the locations of all basis functions std::vector<DimVec<double>* > *bary; /** \name static dim-degree-arrays @@ -217,11 +191,8 @@ namespace AMDiS { protected: - /** \brief - * Pointer to the used refineInter function - */ - void (*refineInter_fct)(DOFIndexed<double> *, RCNeighbourList*, int, - BasisFunction*); + /// Pointer to the used refineInter function + void (*refineInter_fct)(DOFIndexed<double> *, RCNeighbourList*, int, BasisFunction*); /** \name refineInter functions * \{ @@ -250,11 +221,8 @@ namespace AMDiS { BasisFunction*); /** \} */ - /** \brief - * Pointer to the used coarseRestr function - */ - void (*coarseRestr_fct)(DOFIndexed<double> *, RCNeighbourList*, int, - BasisFunction*); + /// Pointer to the used coarseRestr function + void (*coarseRestr_fct)(DOFIndexed<double> *, RCNeighbourList*, int, BasisFunction*); /** \name coarseRestr functions * \{ @@ -283,11 +251,8 @@ namespace AMDiS { BasisFunction*); /** \} */ - /** \brief - * Pointer to the used coarseInter function - */ - void (*coarseInter_fct)(DOFIndexed<double> *, RCNeighbourList*, int, - BasisFunction*); + /// Pointer to the used coarseInter function + void (*coarseInter_fct)(DOFIndexed<double> *, RCNeighbourList*, int, BasisFunction*); /** \name coarseInter functions * \{ @@ -317,12 +282,7 @@ namespace AMDiS { /** \} */ - - // ===== Phi ========================================================= - - /** \brief - * AbstractFunction which implements lagrange basis functions - */ + /// AbstractFunction which implements lagrange basis functions class Phi : public BasFctType { public: @@ -333,26 +293,19 @@ namespace AMDiS { */ Phi(Lagrange* owner, GeoIndex position, int positionIndex, int nodeIndex); - /** \brief - * Destructor - */ + /// Destructor virtual ~Phi(); private: - /** \brief - * vertices needed for evaluation of this function - */ + /// vertices needed for evaluation of this function int* vertices; - /** \brief - * Pointer to the evaluating function - */ + /// Pointer to the evaluating function double (*func)(const DimVec<double>& lambda, int* vert); - /** \brief - * Returns \ref func(lambda, vertices) - */ - inline double operator()(const DimVec<double>& lambda) const { + /// Returns \ref func(lambda, vertices) + inline double operator()(const DimVec<double>& lambda) const + { return func(lambda, vertices); } @@ -362,72 +315,84 @@ namespace AMDiS { // ====== Lagrange, degree = 0 ===================================== // center - inline static double phi0c(const DimVec<double>&, int*) { + inline static double phi0c(const DimVec<double>&, int*) + { return 1.0; } // ====== Lagrange, degree = 1 ===================================== // vertex - inline static double phi1v(const DimVec<double>& lambda, int* vertices) { + inline static double phi1v(const DimVec<double>& lambda, int* vertices) + { return lambda[vertices[0]]; } // ====== Lagrange, degree = 2 ===================================== // vertex - inline static double phi2v(const DimVec<double>& lambda, int* vertices) { + inline static double phi2v(const DimVec<double>& lambda, int* vertices) + { return lambda[vertices[0]] * (2.0 * lambda[vertices[0]] - 1.0); } // edge - inline static double phi2e(const DimVec<double>& lambda, int* vertices) { + inline static double phi2e(const DimVec<double>& lambda, int* vertices) + { return (4.0 * lambda[vertices[0]] * lambda[vertices[1]]); } // ====== Lagrange, degree = 3 ===================================== // vertex - inline static double phi3v(const DimVec<double>& lambda, int* vertices) { + inline static double phi3v(const DimVec<double>& lambda, int* vertices) + { return (4.5 * (lambda[vertices[0]] - 1.0) * lambda[vertices[0]] + 1.0) * lambda[vertices[0]]; } // edge - inline static double phi3e(const DimVec<double>& lambda, int* vertices) { + inline static double phi3e(const DimVec<double>& lambda, int* vertices) + { return (13.5 * lambda[vertices[0]] - 4.5) * lambda[vertices[0]] * lambda[vertices[1]]; } // face - inline static double phi3f(const DimVec<double>& lambda, int* vertices) { + inline static double phi3f(const DimVec<double>& lambda, int* vertices) + { return 27.0 * lambda[vertices[0]] * lambda[vertices[1]] * lambda[vertices[2]]; } // ====== Lagrange, degree = 4 ====================================== // vertex - inline static double phi4v(const DimVec<double>& lambda, int* vertices) { + inline static double phi4v(const DimVec<double>& lambda, int* vertices) + { return (((32.0 * lambda[vertices[0]] - 48.0) * lambda[vertices[0]] + 22.0) * lambda[vertices[0]] - 3.0) * lambda[vertices[0]] / 3.0; } // edge - inline static double phi4e0(const DimVec<double>& lambda, int* vertices) { + inline static double phi4e0(const DimVec<double>& lambda, int* vertices) + { return ((128.0 * lambda[vertices[0]] - 96.0) * lambda[vertices[0]] + 16.0) * lambda[vertices[0]] * lambda[vertices[1]] / 3.0; } - inline static double phi4e1(const DimVec<double>& lambda, int* vertices) { + inline static double phi4e1(const DimVec<double>& lambda, int* vertices) + { return (4.0 * lambda[vertices[0]] - 1.0) * lambda[vertices[0]] * (4.0 * lambda[vertices[1]] - 1.0) * lambda[vertices[1]] * 4.0; } // face - inline static double phi4f(const DimVec<double>& lambda, int* vertices) { + inline static double phi4f(const DimVec<double>& lambda, int* vertices) + { return (4.0 * lambda[vertices[0]] - 1.0) * lambda[vertices[0]] * lambda[vertices[1]] * lambda[vertices[2]] * 32.0; } // center - inline static double phi4c(const DimVec<double>& lambda, int* vertices) { + inline static double phi4c(const DimVec<double>& lambda, int* vertices) + { return 256.0 * lambda[vertices[0]] * lambda[vertices[1]] * lambda[vertices[2]] * lambda[vertices[3]]; } @@ -438,8 +403,6 @@ namespace AMDiS { - // ===== GrdPhi ====================================================== - /** \brief * AbstractFunction which implements gradients of lagrange basis functions. * See \ref Phi @@ -452,11 +415,12 @@ namespace AMDiS { virtual ~GrdPhi(); private: int* vertices; - void (*func)(const DimVec<double>& lambda, - int* vertices_, + + void (*func)(const DimVec<double>& lambda, int* vertices_, DimVec<double>& result); - inline void operator()(const DimVec<double>& lambda, DimVec<double>& result) const { + inline void operator()(const DimVec<double>& lambda, DimVec<double>& result) const + { func(lambda, vertices, result); } @@ -598,9 +562,6 @@ namespace AMDiS { - - // ===== D2Phi ======================================================= - /** \brief * AbstractFunction which implements second derivatives of Lagrange basis * functions. See \ref Phi diff --git a/AMDiS/src/LeafData.h b/AMDiS/src/LeafData.h index 06ba8441b383b653171718e51981a9fcdfe07861..1d3161d760fa760d3c57860c28226f19dbbd64d1 100644 --- a/AMDiS/src/LeafData.h +++ b/AMDiS/src/LeafData.h @@ -44,25 +44,24 @@ namespace AMDiS { public LeafDataEstimatableInterface { public: - inline bool isOfType(int type) const { - if (type == ESTIMATABLE) { + inline bool isOfType(int type) const + { + if (type == ESTIMATABLE) return true; - } - + return false; } class Creator : public CreatorInterface<ElementData> { public: - ElementData* create() { + ElementData* create() + { return new LeafDataEstimatable; } }; - /** \brief - * constructor - */ + /// constructor LeafDataEstimatable(ElementData *decorated = NULL) : ElementData(decorated), errorEstimate(0.0) @@ -70,39 +69,33 @@ namespace AMDiS { /** \brief * Refinement of parent to child1 and child2. - * @return true: must this ElementData, else not allowed to delete it + * @return true: must this ElementData, else not allowed to delete it */ bool refineElementData(Element* parent, Element* child1, Element* child2, int elType); - /** \brief - * - */ void coarsenElementData(Element* parent, Element* thisChild, Element* otherChild, int elTypeParent); - /** \brief - * Sets \ref errorEstimate - */ - inline void setErrorEstimate(int, double est) { + /// Sets \ref errorEstimate + inline void setErrorEstimate(int, double est) + { errorEstimate = est; } - /** \brief - * Returns \ref errorEstimate - */ - inline double getErrorEstimate(int) { + /// Returns \ref errorEstimate + inline double getErrorEstimate(int) + { return errorEstimate; } - /** \brief - * Implements ElementData::clone(). - */ - virtual ElementData *clone() const { + /// Implements ElementData::clone(). + virtual ElementData *clone() const + { // create new estimatable leaf data LeafDataEstimatable *newObj = new LeafDataEstimatable(NULL); @@ -115,14 +108,14 @@ namespace AMDiS { return newObj; } - /** \brief - * Returns the name of element data type. - */ - inline std::string getTypeName() const { + /// Returns the name of element data type. + inline std::string getTypeName() const + { return "LeafDataEstimatable"; } - inline const int getTypeID() const { + inline const int getTypeID() const + { return ESTIMATABLE; } @@ -149,59 +142,50 @@ namespace AMDiS { class Creator : public CreatorInterface<ElementData> { public: - ElementData* create() { + ElementData* create() + { return new LeafDataEstimatableVec; } }; - inline bool isOfType(int type) const { + inline bool isOfType(int type) const + { if (type == ESTIMATABLE) return true; return false; } - /** \brief - * constructor - */ + /// constructor LeafDataEstimatableVec(ElementData *decorated = NULL) : ElementData(decorated) - { - } + {} - /** \brief - * Refinement of parent to child1 and child2. - */ + /// Refinement of parent to child1 and child2. bool refineElementData(Element* parent, Element* child1, Element* child2, int elType); - /** \brief - * - */ void coarsenElementData(Element* parent, Element* thisChild, Element* otherChild, int elTypeParent); - /** \brief - * Sets \ref errorEstimate - */ - inline void setErrorEstimate(int index, double est) { + /// Sets \ref errorEstimate + inline void setErrorEstimate(int index, double est) + { errorEstimate[index] = est; } - /** \brief - * Returns \ref errorEstimate - */ - inline double getErrorEstimate(int index) { + /// Returns \ref errorEstimate + inline double getErrorEstimate(int index) + { return errorEstimate[index]; } - /** \brief - * Implements ElementData::clone(). - */ - virtual ElementData *clone() const { + /// Implements ElementData::clone(). + virtual ElementData *clone() const + { // create new estimatable leaf data LeafDataEstimatableVec *newObj = new LeafDataEstimatableVec(NULL); @@ -221,7 +205,7 @@ namespace AMDiS { unsigned int size = errorEstimate.size(); out.write(reinterpret_cast<const char*>(&size), sizeof(unsigned int)); std::map<int, double>::iterator it; - for(it = errorEstimate.begin(); it != errorEstimate.end(); ++it) { + for (it = errorEstimate.begin(); it != errorEstimate.end(); ++it) { out.write(reinterpret_cast<const char*>(&(it->first)), sizeof(int)); out.write(reinterpret_cast<const char*>(&(it->second)), sizeof(double)); } @@ -232,7 +216,7 @@ namespace AMDiS { ElementData::deserialize(in); unsigned size; in.read(reinterpret_cast<char*>(&size), sizeof(unsigned int)); - for(unsigned int i = 0; i < size; i++) { + for (unsigned int i = 0; i < size; i++) { int index; double estimate; in.read(reinterpret_cast<char*>(&index), sizeof(int)); @@ -241,11 +225,13 @@ namespace AMDiS { } } - std::string getTypeName() const { + std::string getTypeName() const + { return "LeafDataEstimatableVec"; } - inline const int getTypeID() const { + inline const int getTypeID() const + { return ESTIMATABLE; } @@ -258,14 +244,10 @@ namespace AMDiS { public: virtual ~LeafDataCoarsenableInterface() {} - /** \brief - * Sets \ref coarseningError - */ + /// Sets \ref coarseningError virtual void setCoarseningErrorEstimate(int index, double est) = 0; - /** \brief - * Returns \ref coarseningError - */ + /// Returns \ref coarseningError virtual double getCoarseningErrorEstimate(int index) = 0; }; @@ -282,15 +264,14 @@ namespace AMDiS { } }; - inline bool isOfType(int type) const { + inline bool isOfType(int type) const + { if(type == COARSENABLE) return true; return false; } - /** \brief - * constructor - */ + /// constructor LeafDataCoarsenable(ElementData *decorated = NULL) : ElementData(decorated), coarseningError(0.00) @@ -299,25 +280,19 @@ namespace AMDiS { ~LeafDataCoarsenable() {} - /** \brief - * Refinement of parent to child1 and child2. - */ + /// Refinement of parent to child1 and child2. bool refineElementData(Element* parent, Element* child1, Element* child2, int elType); - /** \brief - * Refinement of parent to child1 and child2. - */ + /// Refinement of parent to child1 and child2. void coarsenElementData(Element* parent, Element* thisChild, Element* otherChild, int elTypeParent); - /** \brief - * Implements ElementData::clone(). - */ + /// Implements ElementData::clone(). inline ElementData *clone() const { // create new estimatable leaf data LeafDataCoarsenable *newObj = new LeafDataCoarsenable(NULL); @@ -329,22 +304,18 @@ namespace AMDiS { return newObj; } - /** \brief - * Sets \ref coarseningError - */ - virtual void setCoarseningErrorEstimate(int , double est) { + /// Sets \ref coarseningError + virtual void setCoarseningErrorEstimate(int , double est) + { coarseningError = est; } - /** \brief - * Returns \ref coarseningError - */ - virtual double getCoarseningErrorEstimate(int) { + /// Returns \ref coarseningError + virtual double getCoarseningErrorEstimate(int) + { return coarseningError; } - // ===== Serializable implementation ===== - void serialize(std::ostream& out) { ElementData::serialize(out); @@ -366,7 +337,7 @@ namespace AMDiS { } private: - double coarseningError; /**< \brief coarsening error */ + double coarseningError; }; @@ -378,28 +349,27 @@ namespace AMDiS { class Creator : public CreatorInterface<ElementData> { public: - ElementData* create() { + ElementData* create() + { return new LeafDataCoarsenableVec; - }; + } }; - inline bool isOfType(int type) const { - if(type == COARSENABLE) + inline bool isOfType(int type) const + { + if (type == COARSENABLE) return true; return false; - }; + } - /** \brief - * constructor - */ + /// constructor LeafDataCoarsenableVec(ElementData *decorated = NULL) : ElementData(decorated) - {}; + {} - /** \brief - * Implements ElementData::clone(). - */ - inline ElementData *clone() const { + /// Implements ElementData::clone(). + inline ElementData *clone() const + { // create new estimatable leaf data LeafDataCoarsenableVec *newObj = new LeafDataCoarsenableVec(NULL); @@ -411,39 +381,31 @@ namespace AMDiS { // return the clone return newObj; - }; + } - /** \brief - * Refinement of parent to child1 and child2. - */ + /// Refinement of parent to child1 and child2. bool refineElementData(Element* parent, Element* child1, Element* child2, int elType); - /** \brief - * Refinement of parent to child1 and child2. - */ + /// Refinement of parent to child1 and child2. void coarsenElementData(Element* parent, Element* thisChild, Element* otherChild, int elTypeParent); - /** \brief - * Sets \ref coarseningError - */ - virtual void setCoarseningErrorEstimate(int index, double est) { + /// Sets \ref coarseningError + virtual void setCoarseningErrorEstimate(int index, double est) + { coarseningError[index] = est; - }; + } - /** \brief - * Returns \ref coarseningError - */ - virtual double getCoarseningErrorEstimate(int index) { + /// Returns \ref coarseningError + virtual double getCoarseningErrorEstimate(int index) + { return coarseningError[index]; - }; - - // ===== Serializable implementation ===== + } void serialize(std::ostream& out) { @@ -451,35 +413,38 @@ namespace AMDiS { unsigned int size = coarseningError.size(); out.write(reinterpret_cast<const char*>(&size), sizeof(unsigned int)); std::map<int, double>::iterator it; - for(it = coarseningError.begin(); it != coarseningError.end(); ++it) { + for (it = coarseningError.begin(); it != coarseningError.end(); ++it) { out.write(reinterpret_cast<const char*>(&(it->first)), sizeof(int)); out.write(reinterpret_cast<const char*>(&(it->second)), sizeof(double)); } - }; + } void deserialize(std::istream& in) { ElementData::deserialize(in); - unsigned int i, size; + unsigned int size; in.read(reinterpret_cast<char*>(&size), sizeof(unsigned int)); - for(i = 0; i < size; i++) { + for (unsigned int i = 0; i < size; i++) { int index; double estimate; in.read(reinterpret_cast<char*>(&index), sizeof(int)); in.read(reinterpret_cast<char*>(&estimate), sizeof(double)); coarseningError[index] = estimate; } - }; + } std::string getTypeName() const { return "LeafDataCoarsenableVec"; - }; + } - inline const int getTypeID() const { return COARSENABLE; }; + inline const int getTypeID() const + { + return COARSENABLE; + } private: - std::map<int, double> coarseningError; /**< \brief coarsening error */ + std::map<int, double> coarseningError; }; @@ -490,13 +455,15 @@ namespace AMDiS { class Creator : public CreatorInterface<ElementData> { public: - ElementData* create() { + ElementData* create() + { return new LeafDataPeriodic; - }; + } }; - inline bool isOfType(int type) const { - if(type == PERIODIC) + inline bool isOfType(int type) const + { + if (type == PERIODIC) return true; return false; }; @@ -506,36 +473,38 @@ namespace AMDiS { { public: PeriodicInfo() - : periodicCoords(NULL) - {}; + : periodicCoords(NULL) + {} PeriodicInfo(int mode, BoundaryType t, int side, const DimVec<WorldVector<double> > *coords); - virtual ~PeriodicInfo() { - if(periodicCoords) delete periodicCoords; - }; + virtual ~PeriodicInfo() + { + if (periodicCoords) + delete periodicCoords; + } PeriodicInfo(const PeriodicInfo &rhs); - void serialize(std::ostream &out) { + void serialize(std::ostream &out) + { out.write(reinterpret_cast<const char*>(&periodicMode), sizeof(int)); out.write(reinterpret_cast<const char*>(&type), sizeof(BoundaryType)); out.write(reinterpret_cast<const char*>(&elementSide), sizeof(int)); - int i, size; - if(periodicCoords) { - size = periodicCoords->getSize(); + if (periodicCoords) { + int size = periodicCoords->getSize(); out.write(reinterpret_cast<const char*>(&size), sizeof(int)); - for(i = 0; i < size; i++) { + for (int i = 0; i < size; i++) { (*periodicCoords)[i].serialize(out); } } else { - size = 0; + int size = 0; out.write(reinterpret_cast<const char*>(&size), sizeof(int)); } - }; + } void deserialize(std::istream &in) { @@ -543,47 +512,44 @@ namespace AMDiS { in.read(reinterpret_cast<char*>(&type), sizeof(BoundaryType)); in.read(reinterpret_cast<char*>(&elementSide), sizeof(int)); - int i, size; + int size; in.read(reinterpret_cast<char*>(&size), sizeof(int)); - if(periodicCoords) delete periodicCoords; - if(size == 0) { + if (periodicCoords) + delete periodicCoords; + if (size == 0) { periodicCoords = NULL; } else { periodicCoords = new DimVec<WorldVector<double> >(size-1, NO_INIT); - for(i = 0; i < size; i++) { + for (int i = 0; i < size; i++) (*periodicCoords)[i].deserialize(in); - } } - }; + } int periodicMode; + BoundaryType type; + int elementSide; + DimVec<WorldVector<double> > *periodicCoords; }; public: - /** \brief - * constructor - */ + /// constructor LeafDataPeriodic(ElementData *decorated = NULL) : ElementData(decorated) - //newDOF(-1) - {}; + {} - /** \brief - * Destructor. - */ - ~LeafDataPeriodic() {}; + /// Destructor. + ~LeafDataPeriodic() {} - /** \brief - * Implements LeafData::clone(). - */ - inline ElementData *clone() const { + /// Implements LeafData::clone(). + inline ElementData *clone() const + { LeafDataPeriodic *newObj = new LeafDataPeriodic; newObj->decorated_ = ElementData::clone(); return newObj; - }; + } inline void addPeriodicInfo(int mode, BoundaryType type, @@ -592,55 +558,50 @@ namespace AMDiS { { PeriodicInfo periodicInfo(mode, type, side, coords); periodicInfoList.push_back(periodicInfo); - }; + } - inline std::list<PeriodicInfo>& getInfoList() { + inline std::list<PeriodicInfo>& getInfoList() + { return periodicInfoList; - }; - - // ===== Serializable implementation ===== + } void serialize(std::ostream& out) { ElementData::serialize(out); - //out.write(reinterpret_cast<const char*>(&newDOF), sizeof(DegreeOfFreedom)); unsigned int size = periodicInfoList.size(); out.write(reinterpret_cast<const char*>(&size), sizeof(unsigned int)); std::list<PeriodicInfo>::iterator it; - for(it = periodicInfoList.begin(); it != periodicInfoList.end(); ++it) { + for (it = periodicInfoList.begin(); it != periodicInfoList.end(); ++it) it->serialize(out); - } - }; + } void deserialize(std::istream& in) { ElementData::deserialize(in); - //in.read(reinterpret_cast<char*>(&newDOF), sizeof(DegreeOfFreedom)); unsigned int size; in.read(reinterpret_cast<char*>(&size), sizeof(unsigned int)); periodicInfoList.resize(size); std::list<PeriodicInfo>::iterator it; - for(it = periodicInfoList.begin(); it != periodicInfoList.end(); ++it) { + for (it = periodicInfoList.begin(); it != periodicInfoList.end(); ++it) it->deserialize(in); - } - }; + } std::string getTypeName() const { return "LeafDataPeriodic"; - }; + } - inline const int getTypeID() const { return PERIODIC; }; + inline const int getTypeID() const + { + return PERIODIC; + } bool refineElementData(Element* parent, Element* child1, Element* child2, int elType); - // inline void setNewDOF(DegreeOfFreedom nd) { newDOF = nd; }; - private: - //DegreeOfFreedom newDOF; std::list<PeriodicInfo> periodicInfoList; friend class LeafDataPeriodicRefinable; diff --git a/AMDiS/src/Line.h b/AMDiS/src/Line.h index b5088aacced09f6e59aa7ff9787073d689f8bb83..ed09dd452a74fb413d3e7bb59f5359629edab4cc 100644 --- a/AMDiS/src/Line.h +++ b/AMDiS/src/Line.h @@ -37,37 +37,32 @@ namespace AMDiS { class Line : public Element { public: - /** \brief - * calls base class contructor. - */ + /// calls base class contructor. Line(Mesh* aMesh) : Element(aMesh) {} - /** \brief - * implements Element::getVertexOfEdge - */ - inline int getVertexOfEdge(int i, int j) const { + /// implements Element::getVertexOfEdge + inline int getVertexOfEdge(int i, int j) const + { return vertexOfEdge[i][j]; } - /** \brief - * implements Element::getVertexOfPosition - */ + /// implements Element::getVertexOfPosition virtual int getVertexOfPosition(GeoIndex position, int positionIndex, int vertexIndex) const; - virtual int getPositionOfVertex(int side, int vertex) const { + virtual int getPositionOfVertex(int side, int vertex) const + { static int positionOfVertex[2][2] = {{0,-1},{-1,0}}; return positionOfVertex[side][vertex]; } - /** \brief - * implements Element::getGeo - */ - inline int getGeo(GeoIndex i) const { - switch(i) { + /// implements Element::getGeo + inline int getGeo(GeoIndex i) const + { + switch (i) { case VERTEX: case PARTS: case NEIGH: return 2; break; @@ -86,39 +81,35 @@ namespace AMDiS { } } - inline int getEdgeOfFace(int /* face */, int /*edge*/ ) const { + inline int getEdgeOfFace(int /* face */, int /*edge*/ ) const + { ERROR_EXIT("called for a line\n"); return 0; } - /** \brief - * implements Element::sortFaceIndices - */ + /// implements Element::sortFaceIndices const FixVec<int,WORLD>& sortFaceIndices(int face, FixVec<int,WORLD> *vec) const; - /** \brief - * implements Element::clone - */ - inline Element *clone() { + /// implements Element::clone + inline Element *clone() + { return new Line(mesh); } - /** \brief - * implements Element::getSideOfChild() - */ - virtual int getSideOfChild(int child, int side, int ) const { + /// implements Element::getSideOfChild() + virtual int getSideOfChild(int child, int side, int ) const + { FUNCNAME("Line::getSideOfChild()"); TEST_EXIT_DBG(child==0 || child==1)("child must be in (0,1)\n"); TEST_EXIT_DBG(side >= 0 && side <= 1)("side must be between 0 and 1\n"); return sideOfChild[child][side]; } - /** \brief - * implements Element::getVertexOfParent() - */ - virtual int getVertexOfParent(int child, int side, int) const { + /// implements Element::getVertexOfParent() + virtual int getVertexOfParent(int child, int side, int) const + { FUNCNAME("Line::getVertexOfParent()"); TEST_EXIT_DBG(child==0 || child==1)("child must be in (0,1)\n"); TEST_EXIT_DBG(side >= 0 && side <= 2)("side must be between 0 and 1\n"); @@ -126,40 +117,33 @@ namespace AMDiS { } - /** \brief - * implements Element::hasSide - */ - inline bool hasSide(Element* /*sideElem*/) const { + /// implements Element::hasSide + inline bool hasSide(Element* /*sideElem*/) const + { ERROR_EXIT("a Line has no side elements!\n"); return false; } - /** \brief - * implements Element::isLine. Returns true because this element is a Line - */ - inline bool isLine() const { + /// implements Element::isLine. Returns true because this element is a Line + inline bool isLine() const + { return true; } - /** \brief - * implements Element::isTriangle. Returns false because this element is a - * Line - */ - inline bool isTriangle() const { + /// implements Element::isTriangle. Returns false because this element is a Line + inline bool isTriangle() const + { return false; } - /** \brief - * implements Element::isTetrahedron. Returns false because this element is - * a Line - */ - inline bool isTetrahedron() const { + /// implements Element::isTetrahedron. Returns false because this element is a Line + inline bool isTetrahedron() const + { return false; } - - // ===== Serializable implementation ===== - std::string getTypeName() const { + std::string getTypeName() const + { return "Line"; } diff --git a/AMDiS/src/ParallelDomainProblem.cc b/AMDiS/src/ParallelDomainProblem.cc index 5dd470c79df08f49d58a87348979ccdb90d44506..0a290ef59073ffbdeb1035a25626aaec34298cd9 100644 --- a/AMDiS/src/ParallelDomainProblem.cc +++ b/AMDiS/src/ParallelDomainProblem.cc @@ -46,128 +46,37 @@ namespace AMDiS { partitionMesh(adaptInfo); - /// === Determine to each dof the set of partitions the dof belongs to. === + /// === Get rank information about DOFs. === - std::map<const DegreeOfFreedom*, std::set<int> > partitionDofs; - TraverseStack stack; - ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL); - while (elInfo) { - Element *element = elInfo->getElement(); + // Stores to each DOF pointer the set of ranks the DOF is part of. + std::map<const DegreeOfFreedom*, std::set<int> > partitionDOFs; + // 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; - // Determine to each dof the partition(s) it corresponds to. - for (int i = 0; i < 3; i++) - partitionDofs[element->getDOF(i)].insert(partitionVec[element->getIndex()]); - - elInfo = stack.traverseNext(elInfo); - } - - /// === Determine the set of ranks dofs and the dofs ownership at the boundary. === - - std::vector<const DegreeOfFreedom*> rankDofs; - for (std::map<const DegreeOfFreedom*, std::set<int> >::iterator it = partitionDofs.begin(); - it != partitionDofs.end(); - ++it) { - for (std::set<int>::iterator itpart1 = it->second.begin(); - itpart1 != it->second.end(); - ++itpart1) { - if (*itpart1 == mpiRank) { - if (it->second.size() == 1) { - rankDofs.push_back(it->first); - } else { - // This dof is at the ranks boundary. It is owned by the rank only if - // the rank number is the highest of all ranks containing this dof. - - bool insert = true; - int highestRank = mpiRank; - for (std::set<int>::iterator itpart2 = it->second.begin(); - itpart2 != it->second.end(); - ++itpart2) { - if (*itpart2 > mpiRank) - insert = false; - - if (*itpart2 > highestRank) - highestRank = *itpart2; - } - - if (insert) - rankDofs.push_back(it->first); - - boundaryDofs[it->first] = highestRank; - } - } - } - } + createDOFMemberInformation(partionDOFs, rankDOFs, boundaryDOFs); // === Create interior boundary information === - elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_NEIGH); - while (elInfo) { - Element *element = elInfo->getElement(); - - // Hidde elements which are not part of ranks partition. - PartitionElementData *partitionData = - dynamic_cast<PartitionElementData*>(element->getElementData(PARTITION_ED)); - if (partitionData->getPartitionStatus() == IN) { - for (int i = 0; i < 3; i++) { - if (!elInfo->getNeighbour(i)) - continue; - - PartitionElementData *neighbourPartitionData = - dynamic_cast<PartitionElementData*>(elInfo->getNeighbour(i)->getElementData(PARTITION_ED)); - if (neighbourPartitionData->getPartitionStatus() == OUT) { - AtomicBoundary& bound = interiorBoundary. - 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; - } - } - } - - elInfo = stack.traverseNext(elInfo); - } + createInteriorBoundaryInfo(); // === Remove all macro elements that are not part of the rank partition. === - std::vector<MacroElement*> macrosToRemove; - for (std::deque<MacroElement*>::iterator it = mesh->firstMacroElement(); - it != mesh->endOfMacroElements(); - ++it) { - PartitionElementData *partitionData = - dynamic_cast<PartitionElementData*> - ((*it)->getElement()->getElementData(PARTITION_ED)); - if (partitionData->getPartitionStatus() != IN) - macrosToRemove.push_back(*it); - } - - mesh->removeMacroElements(macrosToRemove); - - // === Create local and global dofs ordering. === + removeMacroElements(); - int *gOrder = (int*)(malloc(sizeof(int) * rankDofs.size())); - int *lOrder = (int*)(malloc(sizeof(int) * rankDofs.size())); - for (std::vector<const DegreeOfFreedom*>::iterator it = rankDofs.begin(); - it != rankDofs.end(); ++it) - gOrder[nRankDOFs++] = (*it)[0]; + // === Get starting position for global rank dof ordering ==== int rstart = 0; + int nRankDOFs = rankDofs.size(); MPI_Scan(&nRankDOFs, &rstart, 1, MPI_INT, MPI_SUM, PETSC_COMM_WORLD); rstart -= nRankDOFs; - - for (int i = 0; i < nRankDOFs; i++) - lOrder[i] = rstart + i; - - AOCreateBasic(PETSC_COMM_WORLD, nRankDOFs, gOrder, lOrder, &applicationOrdering); - - free(gOrder); - free(lOrder); + /// === Create information which dof indices must be send and which must be received. === std::map<int, std::map<DegreeOfFreedom, DegreeOfFreedom> > sendNewDofs; @@ -195,10 +104,8 @@ namespace AMDiS { for (std::set<int>::iterator itRanks = partitionDofs[it->first].begin(); itRanks != partitionDofs[it->first].end(); ++itRanks) { - if (*itRanks != mpiRank) { + if (*itRanks != mpiRank) sendNewDofs[*itRanks][oldDofIndex] = newDofIndex; - sendDofs[*itRanks].push_back(newDofIndex); - } } } else { // If the boundary dof is not a rank dof, its new dof index, and later @@ -207,6 +114,7 @@ namespace AMDiS { } } + /// === Send and receive the dof indices at boundary. === std::vector<int*> sendBuffers(sendNewDofs.size()); @@ -223,6 +131,8 @@ namespace AMDiS { ++dofIt, c += 2) { sendBuffers[i][c] = dofIt->first; sendBuffers[i][c + 1] = dofIt->second; + + sendDofs[sendIt->first].push_back(dofIt->second); } mpiComm.Isend(sendBuffers[i], sendIt->second.size() * 2, MPI_INT, sendIt->first, 0); @@ -233,21 +143,32 @@ namespace AMDiS { recvIt != recvNewDofs.end(); ++recvIt, i++) { recvBuffers[i] = new int[recvIt->second.size() * 2]; - + mpiComm.Irecv(recvBuffers[i], recvIt->second.size() * 2, MPI_INT, recvIt->first, 0); } mpiComm.Barrier(); + + /// === Delete send buffers. === - /// === Reset all DOFAdmins of the mesh. === + i = 0; + for (std::map<int, std::map<DegreeOfFreedom, DegreeOfFreedom> >::iterator sendIt = sendNewDofs.begin(); + sendIt != sendNewDofs.end(); + ++sendIt, i++) + delete [] sendBuffers[i]; + + + /// === Change dof indices for rank partition. === + + for (int i = 0; i < static_cast<int>(rankDofs.size()); i++) { + const_cast<DegreeOfFreedom*>(rankDofs[i])[0] = i; + mapLocalGlobalDOFs[i] = rstart + i; + mapGlobalLocalDOFs[rstart + i] = i; + isRankDOF[i] = true; + } - int nAdmins = mesh->getNumberOfDOFAdmin(); - for (int i = 0; i < nAdmins; i++) - for (int j = 0; j < mesh->getDOFAdmin(i).getSize(); j++) - const_cast<DOFAdmin&>(mesh->getDOFAdmin(i)).setDOFFree(j, true); - /// === Change dof indices at boundary from other ranks. === @@ -257,19 +178,21 @@ namespace AMDiS { ++recvIt, i++) { for (int j = 0; j < static_cast<int>(recvIt->second.size()); j++) { - for (std::map<const DegreeOfFreedom*, int>::iterator dofIt = boundaryDofs.begin(); - dofIt != boundaryDofs.end(); - ++dofIt) { - if ((dofIt->first)[0] == recvBuffers[i][j * 2]) { - int newDof = recvBuffers[i][j * 2 + 1]; - recvDofs[recvIt->first].push_back(newDof); + int oldDof = recvBuffers[i][j * 2]; + int newDof = recvBuffers[i][j * 2 + 1]; + int newLocalDof = mapLocalGlobalDOFs.size(); - const_cast<DegreeOfFreedom*>(dofIt->first)[0] = newDof; - - for (int k = 0; k < nAdmins; k++) - const_cast<DOFAdmin&>(mesh->getDOFAdmin(k)).setDOFFree(newDof, false); + recvDofs[recvIt->first].push_back(newDof); + for (std::map<const DegreeOfFreedom*, int>::iterator dofIt = boundaryDofs.begin(); + dofIt != boundaryDofs.end(); + ++dofIt) { + if ((dofIt->first)[0] == oldDof) { + const_cast<DegreeOfFreedom*>(dofIt->first)[0] = newLocalDof; + mapLocalGlobalDOFs[newLocalDof] = newDof; + mapGlobalLocalDOFs[newDof] = newLocalDof; + isRankDOF[newLocalDof] = false; break; } } @@ -278,24 +201,39 @@ namespace AMDiS { delete [] recvBuffers[i]; } - i = 0; - for (std::map<int, std::map<DegreeOfFreedom, DegreeOfFreedom> >::iterator sendIt = sendNewDofs.begin(); - sendIt != sendNewDofs.end(); - ++sendIt, i++) { - delete [] sendBuffers[i]; + + /// === Reset all DOFAdmins of the mesh. === + + int nAdmins = mesh->getNumberOfDOFAdmin(); + for (int i = 0; i < nAdmins; i++) { + for (int j = 0; j < mesh->getDOFAdmin(i).getSize(); j++) + const_cast<DOFAdmin&>(mesh->getDOFAdmin(i)).setDOFFree(j, true); + for (int j = 0; j < mapLocalGlobalDOFs.size(); j++) + const_cast<DOFAdmin&>(mesh->getDOFAdmin(i)).setDOFFree(j, false); } - /// === Change dof indices for rank partition. === - for (int i = 0; i < static_cast<int>(rankDofs.size()); i++) { - const_cast<DegreeOfFreedom*>(rankDofs[i])[0] = rstart + i; + /// === Create local information from sendDofs and recvDofs - for (int k = 0; k < nAdmins; k++) - const_cast<DOFAdmin&>(mesh->getDOFAdmin(k)).setDOFFree(rstart + i, false); + for (std::map<int, std::vector<DegreeOfFreedom> >::iterator it = sendDofs.begin(); + it != sendDofs.end(); + ++it) + for (std::vector<DegreeOfFreedom>::iterator dofIt = it->second.begin(); + dofIt != it->second.end(); + dofIt++) + *dofIt = mapGlobalLocalDOFs[*dofIt]; - } + for (std::map<int, std::vector<DegreeOfFreedom> >::iterator it = recvDofs.begin(); + it != recvDofs.end(); + ++it) + for (std::vector<DegreeOfFreedom>::iterator dofIt = it->second.begin(); + dofIt != it->second.end(); + dofIt++) + *dofIt = mapGlobalLocalDOFs[*dofIt]; + /// === Create petsc matrix. === + int ierr; ierr = MatCreate(PETSC_COMM_WORLD, &petscMatrix); ierr = MatSetSizes(petscMatrix, rankDofs.size(), rankDofs.size(), @@ -312,9 +250,7 @@ namespace AMDiS { } void ParallelDomainProblemBase::exitParallelization(AdaptInfo *adaptInfo) - { - AODestroy(applicationOrdering); - } + {} void ParallelDomainProblemBase::fillPetscMatrix(DOFMatrix *mat, DOFVector<double> *vec) @@ -330,13 +266,13 @@ namespace AMDiS { typedef traits::range_generator<major, Matrix>::type cursor_type; typedef traits::range_generator<nz, cursor_type>::type icursor_type; - std::cout.precision(10); for (cursor_type cursor = begin<major>(mat->getBaseMatrix()), cend = end<major>(mat->getBaseMatrix()); cursor != cend; ++cursor) for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); icursor != icend; ++icursor) if (value(*icursor) != 0.0) { - int r = row(*icursor); - int c = col(*icursor); + int r = mapLocalGlobalDOFs[row(*icursor)]; + int c = mapLocalGlobalDOFs[col(*icursor)]; double v = value(*icursor); + MatSetValues(petscMatrix, 1, &r, 1, &c, &v, ADD_VALUES); } @@ -346,7 +282,7 @@ namespace AMDiS { DOFVector<double>::Iterator dofIt(vec, USED_DOFS); for (dofIt.reset(); !dofIt.end(); ++dofIt) { - int index = dofIt.getDOFIndex(); + int index = mapLocalGlobalDOFs[dofIt.getDOFIndex()]; double value = *dofIt; VecSetValues(petscRhsVec, 1, &index, &value, ADD_VALUES); @@ -396,7 +332,7 @@ namespace AMDiS { for (std::map<int, std::vector<DegreeOfFreedom> >::iterator recvIt = recvDofs.begin(); recvIt != recvDofs.end(); ++recvIt, i++) { - recvBuffers[i] = new double[recvIt->second.size()]; + recvBuffers[i] = new double[recvIt->second.size()]; mpiComm.Irecv(recvBuffers[i], recvIt->second.size(), MPI_DOUBLE, recvIt->first, 0); } @@ -408,27 +344,14 @@ namespace AMDiS { for (std::map<int, std::vector<DegreeOfFreedom> >::iterator recvIt = recvDofs.begin(); recvIt != recvDofs.end(); ++recvIt, i++) { - for (int j = 0; j < recvIt->second.size(); j++) { - if (MPI::COMM_WORLD.Get_rank() == 0) - std::cout << (recvIt->second)[j] << " = " << recvBuffers[i][j] << std::endl; - + for (int j = 0; j < recvIt->second.size(); j++) (*vec)[(recvIt->second)[j]] = recvBuffers[i][j]; - } delete [] recvBuffers[i]; } for (int i = 0; i < sendBuffers.size(); i++) delete [] sendBuffers[i]; - - std::cout << "TEST" << std::endl; - - if (MPI::COMM_WORLD.Get_rank() == 0) { - DOFVector<double>::Iterator it(vec, USED_DOFS); - for (it.reset(); !it.end(); ++it) { - std::cout << it.getDOFIndex() << " = " << *it << std::endl; - } - } } double ParallelDomainProblemBase::setElemWeights(AdaptInfo *adaptInfo) @@ -479,6 +402,119 @@ namespace AMDiS { partitioner->fillCoarsePartitionVec(&partitionVec); } + void ParallelDomainProblemBase::createDOFMemberInfo( + std::map<const DegreeOfFreedom*, std::set<int>& partitionDOFs, + std::vector<const DegreeOfFreedom*>& rankDOFs, + std::map<const DegreeOfFreedom*, int>& boundaryDOFs) + { + /// === Determine to each dof the set of partitions the dof belongs to. === + + TraverseStack stack; + ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL); + while (elInfo) { + Element *element = elInfo->getElement(); + + // Determine to each dof the partition(s) it corresponds to. + for (int i = 0; i < 3; i++) + partitionDOFs[element->getDOF(i)].insert(partitionVec[element->getIndex()]); + + elInfo = stack.traverseNext(elInfo); + } + + /// === Determine the set of ranks dofs and the dofs ownership at the boundary. === + + for (std::map<const DegreeOfFreedom*, std::set<int> >::iterator it = partitionDOFs.begin(); + it != partitionDOFs.end(); + ++it) { + for (std::set<int>::iterator itpart1 = it->second.begin(); + itpart1 != it->second.end(); + ++itpart1) { + if (*itpart1 == mpiRank) { + if (it->second.size() == 1) { + rankDOFs.push_back(it->first); + } else { + // This dof is at the ranks boundary. It is owned by the rank only if + // the rank number is the highest of all ranks containing this dof. + + bool insert = true; + int highestRank = mpiRank; + for (std::set<int>::iterator itpart2 = it->second.begin(); + itpart2 != it->second.end(); + ++itpart2) { + if (*itpart2 > mpiRank) + insert = false; + + if (*itpart2 > highestRank) + highestRank = *itpart2; + } + + if (insert) + rankDOFs.push_back(it->first); + + boundaryDOFs[it->first] = highestRank; + } + + break; + } + } + } + } + + void ParallelDomainProblemBase::createInteriorBoundaryInfo() + { + TraverseStack stack; + ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_NEIGH); + while (elInfo) { + Element *element = elInfo->getElement(); + + // Hidde elements which are not part of ranks partition. + PartitionElementData *partitionData = + dynamic_cast<PartitionElementData*>(element->getElementData(PARTITION_ED)); + if (partitionData->getPartitionStatus() == IN) { + for (int i = 0; i < 3; i++) { + if (!elInfo->getNeighbour(i)) + continue; + + PartitionElementData *neighbourPartitionData = + dynamic_cast<PartitionElementData*>(elInfo->getNeighbour(i)->getElementData(PARTITION_ED)); + if (neighbourPartitionData->getPartitionStatus() == OUT) { + AtomicBoundary& bound = interiorBoundary. + 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; + } + } + } + + elInfo = stack.traverseNext(elInfo); + } + } + + void ParallelDomainProblemBase::removeMacroElements() + { + std::vector<MacroElement*> macrosToRemove; + for (std::deque<MacroElement*>::iterator it = mesh->firstMacroElement(); + it != mesh->endOfMacroElements(); + ++it) { + PartitionElementData *partitionData = + dynamic_cast<PartitionElementData*> + ((*it)->getElement()->getElementData(PARTITION_ED)); + if (partitionData->getPartitionStatus() != IN) + macrosToRemove.push_back(*it); + } + + mesh->removeMacroElements(macrosToRemove); + } + + void ParallelDomainProblemBase::createLocalGlobalOrder() + { + } + + ParallelDomainProblemScal::ParallelDomainProblemScal(const std::string& name, ProblemScal *problem, ProblemInstatScal *problemInstat) @@ -487,6 +523,13 @@ namespace AMDiS { { } + void ParallelDomainProblemScal::initParallelization(AdaptInfo *adaptInfo) + { + ParallelDomainProblemBase::initParallelization(adaptInfo); + + probScal->getSystemMatrix()->setIsRankDOF(isRankDOF); + } + Flag ParallelDomainProblemScal::oneIteration(AdaptInfo *adaptInfo, Flag toDo) { // return iterationIF->oneIteration(adaptInfo, toDo); diff --git a/AMDiS/src/ParallelDomainProblem.h b/AMDiS/src/ParallelDomainProblem.h index 9704091d7a1cf610302f925c94beb2a6938b7559..9e9fb32e5f6a67a5042ddce08dec691a9da4f3b0 100644 --- a/AMDiS/src/ParallelDomainProblem.h +++ b/AMDiS/src/ParallelDomainProblem.h @@ -55,7 +55,7 @@ namespace AMDiS { virtual ~ParallelDomainProblemBase() {} - void initParallelization(AdaptInfo *adaptInfo); + virtual void initParallelization(AdaptInfo *adaptInfo); void exitParallelization(AdaptInfo *adaptInfo); @@ -118,7 +118,35 @@ namespace AMDiS { virtual void serialize(std::ostream&) {} virtual void deserialize(std::istream&) {} + + + protected: + /** \brief + * This function traverses the whole mesh, i.e. before it is really partitioned, + * and collects information about which DOF corresponds to which rank. + * + * @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. + */ + void createDOFMemberInfo(std::map<const DegreeOfFreedom*, std::set<int>& partitionDOFs, + std::vector<const DegreeOfFreedom*>& rankDOFs, + std::map<const DegreeOfFreedom*, int>& boundaryDOFs); + + /** \brief + * Determine the interior boundaries, i.e. boundaries between ranks, and store + * all information about them in \ref interiorBoundary. + */ + void createInteriorBoundaryInfo(); + + /// Removes all macro elements from the mesh that are not part of ranks partition. + void removeMacroElements(); + + void createLocalGlobalOrder(); + protected: /// ProblemIterationInterface *iterationIF; @@ -169,9 +197,6 @@ namespace AMDiS { */ std::map<int, int> oldPartitionVec; - /// Petsc application ordering - AO applicationOrdering; - Mat petscMatrix; Vec petscRhsVec; @@ -187,8 +212,6 @@ namespace AMDiS { */ InteriorBoundary interiorBoundary; - std::map<const DegreeOfFreedom*, int> boundaryDofs; - /** \brief * This map contains for each rank the list of dofs the current rank must send * to exchange solution dofs at the interior boundaries. @@ -200,6 +223,19 @@ namespace AMDiS { * must receive solution values of dofs at the interior boundaries. */ std::map<int, std::vector<DegreeOfFreedom> > recvDofs; + + /// Maps local to global dof indices. + std::map<DegreeOfFreedom, DegreeOfFreedom> mapLocalGlobalDOFs; + + /// Maps global to local dof indices. + std::map<DegreeOfFreedom, DegreeOfFreedom> mapGlobalLocalDOFs; + + /** \brief + * Maps all DOFs in ranks partition to a bool value. If it is true, the DOF is + * owned by the rank. Otherwise, its an interior boundary DOF that is owned by + * another rank. + */ + std::map<DegreeOfFreedom, bool> isRankDOF; }; class ParallelDomainProblemScal : public ParallelDomainProblemBase @@ -209,6 +245,8 @@ namespace AMDiS { ProblemScal *problem, ProblemInstatScal *problemInstat); + void initParallelization(AdaptInfo *adaptInfo); + virtual Flag oneIteration(AdaptInfo *adaptInfo, Flag toDo = FULL_ITERATION); protected: diff --git a/AMDiS/src/ProblemScal.cc b/AMDiS/src/ProblemScal.cc index daa80b7d4520d116aa8778515148fbe3d77cad6e..d54b39111b470455c88b90de2dbb65cdf4d3fe62 100644 --- a/AMDiS/src/ProblemScal.cc +++ b/AMDiS/src/ProblemScal.cc @@ -145,7 +145,6 @@ namespace AMDiS { void ProblemScal::useApplicationOrdering(AO *ao) { applicationOrdering = ao; - systemMatrix->useApplicationOrdering(ao); rhs->useApplicationOrdering(ao); } #endif diff --git a/AMDiS/src/VtkWriter.hh b/AMDiS/src/VtkWriter.hh index 856285ac630bdfd714287c8e0e6f31b563b216a3..8039cea447ecfd9755fcd1bf93c5ebfe216474e5 100644 --- a/AMDiS/src/VtkWriter.hh +++ b/AMDiS/src/VtkWriter.hh @@ -1,6 +1,7 @@ #include <list> #include <vector> #include "DOFVector.h" +#include "mpi.h" namespace AMDiS { @@ -139,16 +140,15 @@ namespace AMDiS { DOFVector< std::list<WorldVector<double> > >::Iterator coordIt(dofCoords, USED_DOFS); file << std::fixed; - + // Write the values for all vertex DOFs. for (intPointIt.reset(), valueIt.reset(), coordIt.reset(); !intPointIt.end(); ++intPointIt, ++valueIt, ++coordIt) { - if (*intPointIt == -2) { - for (int i = 0; i < static_cast<int>(coordIt->size()); i++) - file << " " << (fabs(*valueIt) < 1e-40 ? 0.0 : *valueIt) << "\n"; - } + if (*intPointIt == -2) + for (int i = 0; i < static_cast<int>(coordIt->size()); i++) + file << " " << (fabs(*valueIt) < 1e-40 ? 0.0 : *valueIt) << "\n"; } // For the second dim case, write also the values of the interpolation points.