diff --git a/AMDiS/bin/Makefile.am b/AMDiS/bin/Makefile.am index c650107263f7e31f7ab407b318baa8de2b9234cc..6b3303d4e07391d3893724ae7fe1ebbacae3f73e 100644 --- a/AMDiS/bin/Makefile.am +++ b/AMDiS/bin/Makefile.am @@ -19,7 +19,8 @@ if USE_PARALLEL_AMDIS $(PARALLEL_DIR)/ParallelDomainProblem.h $(PARALLEL_DIR)/ParallelDomainProblem.cc \ $(PARALLEL_DIR)/ParMetisPartitioner.h $(PARALLEL_DIR)/ParMetisPartitioner.cc \ $(PARALLEL_DIR)/PartitionElementData.h - PARALLEL_INCLUDES = -I$(MPI_DIR)/include -I$(PARMETIS_DIR) + PARALLEL_INCLUDES = -I/work/home7/witkowsk/local/include \ + -I$(MPI_DIR)/include -I$(PARMETIS_DIR) else PARALLEL_AMDIS_SOURCES = PARALLEL_INCLUDES = @@ -27,6 +28,10 @@ endif libamdis_la_CXXFLAGS = +if USE_PARALLEL_AMDIS + libamdis_la_CXXFLAGS += -DHAVE_PARALLEL_AMDIS=1 +endif + if ENABLE_UMFPACK libamdis_la_CXXFLAGS += -DHAVE_UMFPACK=1 AMDIS_INCLUDES += -I$(LIB_DIR)/UFconfig \ diff --git a/AMDiS/bin/Makefile.in b/AMDiS/bin/Makefile.in index 3c77cd69ebc7bf1843f5297a72da10bf7c27e7bf..623d04b026aed80e68662a65fcd298850219bd7b 100644 --- a/AMDiS/bin/Makefile.in +++ b/AMDiS/bin/Makefile.in @@ -36,17 +36,18 @@ PRE_UNINSTALL = : POST_UNINSTALL = : build_triplet = @build@ host_triplet = @host@ -@ENABLE_UMFPACK_TRUE@am__append_1 = -DHAVE_UMFPACK=1 -@ENABLE_UMFPACK_TRUE@am__append_2 = -I$(LIB_DIR)/UFconfig \ +@USE_PARALLEL_AMDIS_TRUE@am__append_1 = -DHAVE_PARALLEL_AMDIS=1 +@ENABLE_UMFPACK_TRUE@am__append_2 = -DHAVE_UMFPACK=1 +@ENABLE_UMFPACK_TRUE@am__append_3 = -I$(LIB_DIR)/UFconfig \ @ENABLE_UMFPACK_TRUE@ -I$(LIB_DIR)/AMD/Include \ @ENABLE_UMFPACK_TRUE@ -I$(LIB_DIR)/UMFPACK/Include -@ENABLE_MKL_TRUE@am__append_3 = -DHAVE_MKL=1 -I${MKL_INC} -@ENABLE_DUNE_TRUE@am__append_4 = -DHAVE_DUNE=1 -@ENABLE_DUNE_TRUE@am__append_5 = -I$(DUNE_DIR) -@ENABLE_BOOST_TRUE@am__append_6 = -DHAVE_BOOST=1 -@AMDIS_DEBUG_TRUE@am__append_7 = -g -O0 -Wall -DDEBUG=1 $(OPENMP_FLAG) -ftemplate-depth-30 $(INCLUDES) #-pedantic -@AMDIS_DEBUG_FALSE@am__append_8 = -O2 -Wall -DDEBUG=0 $(OPENMP_FLAG) -ftemplate-depth-30 $(INCLUDES) #-pedantic +@ENABLE_MKL_TRUE@am__append_4 = -DHAVE_MKL=1 -I${MKL_INC} +@ENABLE_DUNE_TRUE@am__append_5 = -DHAVE_DUNE=1 +@ENABLE_DUNE_TRUE@am__append_6 = -I$(DUNE_DIR) +@ENABLE_BOOST_TRUE@am__append_7 = -DHAVE_BOOST=1 +@AMDIS_DEBUG_TRUE@am__append_8 = -g -O0 -Wall -DDEBUG=1 $(OPENMP_FLAG) -ftemplate-depth-30 $(INCLUDES) #-pedantic +@AMDIS_DEBUG_FALSE@am__append_9 = -O2 -Wall -DDEBUG=0 $(OPENMP_FLAG) -ftemplate-depth-30 $(INCLUDES) #-pedantic subdir = bin DIST_COMMON = $(srcdir)/Makefile.am $(srcdir)/Makefile.in ACLOCAL_M4 = $(top_srcdir)/aclocal.m4 @@ -484,7 +485,7 @@ SOURCE_DIR = ../src LIB_DIR = ../lib PARALLEL_DIR = $(SOURCE_DIR) PARMETIS_DIR = ../lib/ParMetis-3.1 -AMDIS_INCLUDES = -I$(SOURCE_DIR) $(am__append_2) $(am__append_5) +AMDIS_INCLUDES = -I$(SOURCE_DIR) $(am__append_3) $(am__append_6) @USE_PARALLEL_AMDIS_FALSE@PARALLEL_AMDIS_SOURCES = @USE_PARALLEL_AMDIS_TRUE@PARALLEL_AMDIS_SOURCES = \ @USE_PARALLEL_AMDIS_TRUE@ $(PARALLEL_DIR)/ConditionalEstimator.h $(PARALLEL_DIR)/ConditionalEstimator.cc \ @@ -498,9 +499,12 @@ AMDIS_INCLUDES = -I$(SOURCE_DIR) $(am__append_2) $(am__append_5) @USE_PARALLEL_AMDIS_TRUE@ $(PARALLEL_DIR)/PartitionElementData.h @USE_PARALLEL_AMDIS_FALSE@PARALLEL_INCLUDES = -@USE_PARALLEL_AMDIS_TRUE@PARALLEL_INCLUDES = -I$(MPI_DIR)/include -I$(PARMETIS_DIR) -libamdis_la_CXXFLAGS = $(am__append_1) $(am__append_3) $(am__append_4) \ - $(am__append_6) $(am__append_7) $(am__append_8) +@USE_PARALLEL_AMDIS_TRUE@PARALLEL_INCLUDES = -I/work/home7/witkowsk/local/include \ +@USE_PARALLEL_AMDIS_TRUE@ -I$(MPI_DIR)/include -I$(PARMETIS_DIR) + +libamdis_la_CXXFLAGS = $(am__append_1) $(am__append_2) $(am__append_4) \ + $(am__append_5) $(am__append_7) $(am__append_8) \ + $(am__append_9) INCLUDES = $(AMDIS_INCLUDES) $(PARALLEL_INCLUDES) libamdis_la_SOURCES = \ $(PARALLEL_AMDIS_SOURCES) \ diff --git a/AMDiS/libtool b/AMDiS/libtool index c1ffc52dc2a18e9aea6a1c54bd62b799a57e937c..fcd0e2f7c8e2ddac64e11b952c5570c79e95b93b 100755 --- a/AMDiS/libtool +++ b/AMDiS/libtool @@ -44,7 +44,7 @@ available_tags=" CXX F77" # ### BEGIN LIBTOOL CONFIG -# Libtool was configured on host p2s148: +# Libtool was configured on host p2s256: # Shell to use when invoking shell scripts. SHELL="/bin/sh" @@ -6760,7 +6760,7 @@ build_old_libs=`case $build_libtool_libs in yes) $echo no;; *) $echo yes;; esac` # End: # ### BEGIN LIBTOOL TAG CONFIG: CXX -# Libtool was configured on host p2s148: +# Libtool was configured on host p2s256: # Shell to use when invoking shell scripts. SHELL="/bin/sh" @@ -7065,7 +7065,7 @@ include_expsyms="" # ### BEGIN LIBTOOL TAG CONFIG: F77 -# Libtool was configured on host p2s148: +# Libtool was configured on host p2s256: # Shell to use when invoking shell scripts. SHELL="/bin/sh" diff --git a/AMDiS/src/AdaptStationary.h b/AMDiS/src/AdaptStationary.h index 0c30c6dbed973d1a09f5df52ae40e157f19ed58c..8bf24778ebf2fd9590ebb905b49e122b6637c213 100644 --- a/AMDiS/src/AdaptStationary.h +++ b/AMDiS/src/AdaptStationary.h @@ -55,27 +55,19 @@ namespace AMDiS { public: MEMORY_MANAGED(AdaptStationary); - /** \brief - * Creates a AdaptStationary object with given name. - */ + /// Creates a AdaptStationary object with given name. AdaptStationary(const std::string &name, ProblemIterationInterface *prob, AdaptInfo *info); - /** \brief - * Destructor - */ + /// Destructor virtual ~AdaptStationary() {} - /** \brief - * Implementation of AdaptBase::adapt() - */ + /// Implementation of AdaptBase::adapt() virtual int adapt(); protected: - /** \brief - * Initialisation - */ + /// Initialisation void initialize(); }; diff --git a/AMDiS/src/DOFMatrix.cc b/AMDiS/src/DOFMatrix.cc index 8138d34ab9ae0a21b07743bcc7de5af8b7aec4de..98c48edc5d1325eb23f4c451bfa5297fd52211d5 100644 --- a/AMDiS/src/DOFMatrix.cc +++ b/AMDiS/src/DOFMatrix.cc @@ -47,6 +47,10 @@ namespace AMDiS { colFESpace->getBasisFcts()->getNumber()); applyDBCs.clear(); + +#ifdef HAVE_PARALLEL_AMDIS + applicationOrdering = NULL; +#endif } DOFMatrix::DOFMatrix(const DOFMatrix& rhs) diff --git a/AMDiS/src/DOFMatrix.h b/AMDiS/src/DOFMatrix.h index da73d1dec1678d3e23b546845e3343a3d91c9f69..49d5ea3c53b823b4de7d7cb95ca1a432240c3342 100644 --- a/AMDiS/src/DOFMatrix.h +++ b/AMDiS/src/DOFMatrix.h @@ -30,6 +30,7 @@ #include <set> #include <memory> #include <list> + #include "Global.h" #include "Flag.h" #include "RCNeighbourList.h" @@ -40,6 +41,10 @@ #include "Boundary.h" #include "Serializable.h" +#ifdef HAVE_PARALLEL_AMDIS +#include "petscao.h" +#endif + namespace AMDiS { // ============================================================================ @@ -420,23 +425,17 @@ namespace AMDiS { return rowFESpace; } - /** \brief - * Returns const \ref colFESpace - */ + /// Returns const \ref colFESpace const FiniteElemSpace* getColFESpace() const { return colFESpace; } - /** \brief - * Returns const \ref rowFESpace - */ + /// Returns const \ref rowFESpace const FiniteElemSpace* getFESpace() const { return rowFESpace; } - /** \brief - * Returns number of rows (\ref matrix.size()) - */ + /// Returns number of rows (\ref matrix.size()) inline int getSize() const { return matrix.size(); } @@ -455,58 +454,44 @@ namespace AMDiS { */ int getNumCols() const; - /** \brief - * Returns \ref name - */ + /// Returns \ref name inline const std::string& getName() const { return name; } - /** \brief - * Resizes \ref matrix to n rows - */ + /// Resizes \ref matrix to n rows inline void resize(int n) { TEST_EXIT_DBG(n >= 0)("Can't resize DOFMatrix to negative size\n"); matrix.resize(n); } - /** \brief - * Returns \ref matrix[i] which is the i-th row - */ + /// Returns \ref matrix[i] which is the i-th row inline const std::vector<MatEntry>& operator[](int i) const { TEST_EXIT_DBG((i >= 0) && (i < (static_cast<int>(matrix.size())))) ("Illegal matrix index %d.\n",i); return matrix[i]; } - /** \brief - * Returns \ref matrix[i] which is the i-th row - */ + /// Returns \ref matrix[i] which is the i-th row inline std::vector<MatEntry>& operator[](int i) { TEST_EXIT_DBG((i >= 0) && (i < (static_cast<int>(matrix.size())))) ("Illegal vector index %d.\n", i); return matrix[i]; } - /** \brief - * Access to \ref matrix[a][b].entry - */ + /// Access to \ref matrix[a][b].entry inline double physAcc(DegreeOfFreedom a, DegreeOfFreedom b) const { return matrix[a][b].entry; } - /** \brief - * Access to \ref matrix[a][b].col - */ + /// Access to \ref matrix[a][b].col inline DegreeOfFreedom physToLogIndex(DegreeOfFreedom a, DegreeOfFreedom b) const { return matrix[a][b].col; } - /** \brief - * Returns physical column index of logical index b in row a - */ + /// Returns physical column index of logical index b in row a inline DegreeOfFreedom logToPhysIndex(DegreeOfFreedom a, DegreeOfFreedom b) const { int size = static_cast<int>(matrix[a].size()); @@ -517,19 +502,13 @@ namespace AMDiS { return -1; } - /** \brief - * Returns value at logical indices a,b - */ + /// Returns value at logical indices a,b double logAcc(DegreeOfFreedom a, DegreeOfFreedom b) const; - /** \brief - * Changes col at logical indices a,b to c - */ + /// Changes col at logical indices a,b to c void changeColOfEntry(DegreeOfFreedom a, DegreeOfFreedom b, DegreeOfFreedom c); - /** \brief - * Changes col of \ref matrix[a][b] to c - */ + /// Changes col of \ref matrix[a][b] to c inline void changePhysColOfEntry(DegreeOfFreedom a, DegreeOfFreedom b, DegreeOfFreedom c) @@ -565,24 +544,16 @@ namespace AMDiS { void removeRowsWithDBC(std::set<int> *rows); - /** \brief - * Prints \ref matrix to stdout - */ + /// Prints \ref matrix to stdout void print() const; - /** \brief - * Prints a row of \ref matrix to stdout - */ + /// Prints a row of \ref matrix to stdout void printRow(int i) const; - /** \brief - * Removes all matrix entries - */ + /// Removes all matrix entries void clear(); - /** \brief - * Test whether \ref matrix is symmetric. Exits if not. - */ + /// Test whether \ref matrix is symmetric. Exits if not. void test(); bool symmetric(); @@ -612,9 +583,7 @@ namespace AMDiS { } void createPictureFile(const char* filename, int dim); - - // ===== Serializable implementation ===== - + void serialize(std::ostream &out) { unsigned int matrixSize = matrix.size(); @@ -640,13 +609,15 @@ namespace AMDiS { } } + /// int memsize(); - protected: - /** \brief - * Used by updateMatrix - */ - //static int assembleFct(ElInfo *el_info); +#ifdef HAVE_PARALLEL_AMDIS + /// Sets the petsc application ordering object to map dof indices. + void useApplicationOrdering(AO *ao) { + applicationOrdering = ao; + } +#endif protected: /** \brief @@ -661,19 +632,13 @@ namespace AMDiS { */ const FiniteElemSpace *colFESpace; - /** \brief - * Name of the DOFMatrix - */ + /// Name of the DOFMatrix std::string name; - /** \brief - * Sparse matrix stored in an STL vector of MatrixRow objects - */ + /// Sparse matrix stored in an STL vector of MatrixRow objects std::vector<MatrixRow> matrix; - /** \brief - * Used while mesh traversal - */ + /// Used while mesh traversal static DOFMatrix *traversePtr; /** \brief @@ -688,28 +653,26 @@ namespace AMDiS { */ std::vector<double*> operatorFactor; - /** \brief - * - */ + /// std::vector<double*> operatorEstFactor; - /** \brief - * - */ + /// BoundaryManager *boundaryManager; - /** \brief - * - */ + /// bool coupleMatrix; - /** \brief - * Temporary variable used in assemble() - */ + /// Temporary variable used in assemble() ElementMatrix *elementMatrix; + /// std::set<int> applyDBCs; +#ifdef HAVE_PARALLEL_AMDIS + /// Petsc application ordering to map dof indices. + AO *applicationOrdering; +#endif + friend class DOFAdmin; friend class DOFVector<double>; friend class DOFVector<unsigned char>; @@ -718,7 +681,8 @@ namespace AMDiS { }; - inline DegreeOfFreedom logToPhysIndex(DOFMatrix *m, DegreeOfFreedom a, DegreeOfFreedom b) + inline DegreeOfFreedom logToPhysIndex(DOFMatrix *m, + DegreeOfFreedom a, DegreeOfFreedom b) { return m->logToPhysIndex(a, b); } @@ -729,16 +693,11 @@ namespace AMDiS { double max(std::vector<MatEntry> *row); - /** \brief - * Addes two matrices. The result matrix is overwritten. - */ + /// Addes two matrices. The result matrix is overwritten. void addDOFMatrix(DOFMatrix *result, const DOFMatrix *a, const DOFMatrix *b); - /** \brief - * Addes matrix a to matrix result. - */ + /// Addes matrix a to matrix result. void addDOFMatrix(DOFMatrix *result, const DOFMatrix *a); - } #endif // AMDIS_DOFMATRIX_H diff --git a/AMDiS/src/DOFVector.h b/AMDiS/src/DOFVector.h index 5b1ae2cee1ec9b3dd024a49a99ed6f2e64139aef..ebd581b0edb1bebb80c51afb4853fae5e220f048 100644 --- a/AMDiS/src/DOFVector.h +++ b/AMDiS/src/DOFVector.h @@ -26,6 +26,10 @@ // ===== includes ============================================================ // =========================================================================== +#include <vector> +#include <memory> +#include <list> + #include "FixVec.h" #include "Global.h" #include "Flag.h" @@ -39,9 +43,10 @@ #include "Serializable.h" #include "DOFMatrix.h" #include "BasisFunction.h" -#include <vector> -#include <memory> -#include <list> + +#ifdef HAVE_PARALLEL_AMDIS +#include "petscao.h" +#endif namespace AMDiS { @@ -584,24 +589,16 @@ namespace AMDiS { */ T max() const; - /** \brief - * Returns absolute maximum of DOFVector - */ + /// Returns absolute maximum of DOFVector T absMax() const; - /** \brief - * Used by interpol while mesh traversal - */ + /// Used by interpol while mesh traversal static int interpolFct(ElInfo* elinfo); - /** \brief - * Prints \ref vec to stdout - */ + /// Prints \ref vec to stdout void print() const; - /** \brief - * - */ + /// int calcMemoryUsage() const; /** \brief @@ -612,8 +609,6 @@ namespace AMDiS { void interpol(DOFVector<T> *v, double factor); - // ===== Serializable implementation ===== - void serialize(std::ostream &out) { unsigned int size = vec.size(); out.write(reinterpret_cast<const char*>(&size), sizeof(unsigned int)); @@ -635,6 +630,13 @@ namespace AMDiS { getRecoveryGradient(DOFVector<WorldVector<T> >*) const; +#ifdef HAVE_PARALLEL_AMDIS + /// Sets the petsc application ordering object to map dof indices. + void useApplicationOrdering(AO *ao) { + applicationOrdering = ao; + } +#endif + protected: /** \brief * Used by Int while mesh traversal @@ -662,45 +664,34 @@ namespace AMDiS { static int DoubleWell_fct(ElInfo* elinfo); protected: - /** \brief - * Name of this DOFVector - */ + /// Name of this DOFVector std::string name; - /** \brief - * FiniteElemSpace of the vector - */ + /// FiniteElemSpace of the vector const FiniteElemSpace *feSpace; - /** \brief - * Data container - */ + /// Data container std::vector<T> vec; - /** \brief - * Specifies what operation should be performed after coarsening - */ + /// Specifies what operation should be performed after coarsening CoarsenOperation coarsenOperation; - /** \brief - * Used by \ref interpol - */ + /// Used by \ref interpol AbstractFunction<T, WorldVector<double> > *interFct; - /** \brief - * Used for mesh traversal - */ + /// Used for mesh traversal static DOFVector<T> *traverseVector; +#ifdef HAVE_PARALLEL_AMDIS + /// Petsc application ordering to map dof indices. + AO *applicationOrdering; +#endif + protected: - /** \brief - * Used while calculating vector norms - */ + /// Used while calculating vector norms static FastQuadrature *quad_fast; - /** \brief - * Stores the last calculated vector norm - */ + /// Stores the last calculated vector norm static double norm; }; @@ -744,9 +735,7 @@ namespace AMDiS { feSpace->getAdmin()->addDOFContainer(this); } - /** \brief - * Deregisters itself at DOFAdmin. - */ + /// Deregisters itself at DOFAdmin. ~DOFVectorDOF() { feSpace->getAdmin()->removeDOFContainer(this); } @@ -759,16 +748,12 @@ namespace AMDiS { return DOFVector<DegreeOfFreedom>::operator[](i); } - /** \brief - * Implements DOFIndexedBase::getSize() - */ + /// Implements DOFIndexedBase::getSize() int getSize() const { return DOFVector<DegreeOfFreedom>::getSize(); } - /** \brief - * Implements DOFIndexedBase::resize() - */ + /// Implements DOFIndexedBase::resize() void resize(int size) { DOFVector<DegreeOfFreedom>::resize(size); } diff --git a/AMDiS/src/DOFVector.hh b/AMDiS/src/DOFVector.hh index 8c29fab6fbbad835d47e18efa9c38584e6a7c847..d476dd76feb56f23e7aae43f31ae27e00a9edaee 100644 --- a/AMDiS/src/DOFVector.hh +++ b/AMDiS/src/DOFVector.hh @@ -66,6 +66,10 @@ namespace AMDiS { coarsenOperation(NO_OPERATION) { init(f, n); + +#ifdef HAVE_PARALLEL_AMDIS + applicationOrdering = NULL; +#endif } template<typename T> diff --git a/AMDiS/src/Element.cc b/AMDiS/src/Element.cc index 57d5e53ff299e75003f68726254714db7b604fd7..c8ac2cb33abafa45ed22fe7fef91cd467592bcbd 100644 --- a/AMDiS/src/Element.cc +++ b/AMDiS/src/Element.cc @@ -41,7 +41,6 @@ namespace AMDiS { newCoord = NULL; elementData = NULL; mark = 0; - hidden = false; if (mesh) { setDOFPtrs(); @@ -139,7 +138,6 @@ namespace AMDiS { el->mesh = mesh; el->index = index; el->mark = mark; - el->hidden = hidden; if (newCoord) { WorldVector<double> *nc = NEW WorldVector<double>(); *nc = *newCoord; @@ -422,9 +420,6 @@ namespace AMDiS { // write mark out.write(reinterpret_cast<const char*>(&mark), sizeof(signed char)); - // write hidden flag - out.write(reinterpret_cast<const char*>(&hidden), sizeof(hidden)); - // write newCoord if (newCoord) { out << "WorldVector\n"; @@ -506,9 +501,6 @@ namespace AMDiS { // read mark in.read(reinterpret_cast<char*>(&mark), sizeof(signed char)); - // read hidden - in.read(reinterpret_cast<char*>(&hidden), sizeof(bool)); - // read newCoord in >> typeName; in.get(); diff --git a/AMDiS/src/Element.h b/AMDiS/src/Element.h index c2bf1d84cf01631f7093a9b35f61bcf1c0828170..8e9700c0677881099d854949d523d23f33439df2 100644 --- a/AMDiS/src/Element.h +++ b/AMDiS/src/Element.h @@ -230,16 +230,6 @@ namespace AMDiS { return newCoord; } - /// Returns \ref hidden, i.e., if the element is visible for mesh traverses. - inline bool isHidden() { - return hidden; - } - - /// Sets \ref hidden, i.e., wether the element is hidden for mesh traverses or not. - inline void setHidden(bool b) { - hidden = b; - } - /** \} */ // ===== setting methods ====================================================== @@ -529,12 +519,6 @@ namespace AMDiS { /// Pointer to Element's leaf data ElementData* elementData; - /* \brief - * If true, this element (and therefor also all it's children) is hidden for - * all mesh traverses. If the value is set to false, it may be traversed. - */ - bool hidden; - /** \brief * This map is used for deletion of all DOFs of all elements of a mesh. Once * a DOF-vector (all DOFS at a node, edge, etc.) is deleted, its address is diff --git a/AMDiS/src/MacroElement.cc b/AMDiS/src/MacroElement.cc index 561189177cc2c0d1904b000cbabe37ea3546417e..2e1fe7883a38f909f33deb436400933cc9450471 100644 --- a/AMDiS/src/MacroElement.cc +++ b/AMDiS/src/MacroElement.cc @@ -26,9 +26,8 @@ namespace AMDiS { MacroElement::~MacroElement() { - if (element) { + if (element) DELETE element; - } } MacroElement& MacroElement::operator=(const MacroElement &el) diff --git a/AMDiS/src/MacroReader.cc b/AMDiS/src/MacroReader.cc index b01b5e4399d41e7b8e9cbc10ffba06b6fa3bcaa6..2409b0ecfd262a915b3fa1c970d61b52d056de91 100644 --- a/AMDiS/src/MacroReader.cc +++ b/AMDiS/src/MacroReader.cc @@ -1100,26 +1100,8 @@ namespace AMDiS { /*********************************************************************/ n_neigh = 1; - // static int faceOfEdge[6][2] = { {2, 3}, - // {1, 3}, - // {1, 2}, - // {0, 3}, - // {0, 2}, - // {0, 1} }; - - // bool periodicEdge; - // if(periodic[faceOfEdge[k][0]] || periodic[faceOfEdge[k][1]]) { - // periodicEdge = true; - // } else { - // periodicEdge = false; - // } - if (newEdge(mesh, (*(mel+i)), k, &n_neigh/*, periodicEdge*/)) { - // if(periodicEdge) { - // mesh->incrementNumberOfEdges(2); - // } else { mesh->incrementNumberOfEdges(1); - // } max_n_neigh = max(max_n_neigh, n_neigh); } } diff --git a/AMDiS/src/Mesh.cc b/AMDiS/src/Mesh.cc index 44f81033db6a367fa40615fc6a5752b0a52ae8a8..d3b380b750f0d57b412e132edf1619819c2c0710 100644 --- a/AMDiS/src/Mesh.cc +++ b/AMDiS/src/Mesh.cc @@ -262,9 +262,51 @@ namespace AMDiS { } } - void Mesh::addMacroElement(MacroElement* m) { - macroElements.push_back(m); - m->setIndex(macroElements.size()); + void Mesh::addMacroElement(MacroElement* me) + { + macroElements.push_back(me); + me->setIndex(macroElements.size()); + } + + void Mesh::removeMacroElement(MacroElement* me) + { + FUNCNAME("Mesh::removeMacroElement()"); + + TEST_EXIT(dim == 2)("Not yet implemented!\n"); + + bool found = false; + + // Remove the macro element from mesh's list of all macro elements. + for (std::deque<MacroElement*>::iterator it = macroElements.begin(); + it != macroElements.end(); + ++it) { + if (*it == me) { + macroElements.erase(it, it + 1); + found = true; + break; + } + } + + TEST_EXIT(found)("Cannot find MacroElement that should be removed!\n"); + + // Go through all neighbours of the macro element and remove this macro element + // to be neighbour of some other macro element. + for (int i = 0; i < dim; i++) { + if (me->getNeighbour(i)) { + for (int j = 0; j < dim; j++) { + if (me->getNeighbour(i)->getNeighbour(j) == me) { + me->getNeighbour(i)->setNeighbour(j, NULL); + } + } + } else { + // There is no neighbour at this edge, so we have to decrease the number + // of edges in the mesh. + nEdges--; + } + } + + nLeaves--; + nElements--; } int Mesh::traverse(int level, Flag flag, int (*el_fct)(ElInfo*)) @@ -496,7 +538,7 @@ namespace AMDiS { FUNCNAME("Mesh::freeDOF()"); TEST_EXIT_DBG(position >= CENTER && position <= FACE) - ("unknown position %d\n",position); + ("unknown position %d\n", position); int ndof = nDOF[position]; if (ndof) { @@ -512,7 +554,7 @@ namespace AMDiS { } TEST_EXIT_DBG(ndof <= MAX_DOF) - ("ndof too big: ndof=%d, MAX_DOF=%d\n",ndof,MAX_DOF); + ("ndof too big: ndof=%d, MAX_DOF=%d\n", ndof, MAX_DOF); for (int i = 0; i < static_cast<int>(admin.size()); i++) { DOFAdmin *localAdmin = admin[i]; @@ -571,8 +613,6 @@ namespace AMDiS { }; } - - bool Mesh::findElInfoAtPoint(const WorldVector<double>& xy, ElInfo *el_info, DimVec<double>& bary, @@ -937,10 +977,10 @@ namespace AMDiS { node.deserialize(in); // read admins - int i, size; + int size; in.read(reinterpret_cast<char*>(&size), sizeof(int)); admin.resize(size, NULL); - for (i = 0; i < size; i++) { + for (int i = 0; i < size; i++) { if (!admin[i]) { admin[i] = NEW DOFAdmin(this); } @@ -952,13 +992,13 @@ namespace AMDiS { std::vector< std::vector<int> > neighbourIndices(size); - for (i = 0; i < static_cast<int>(macroElements.size()); i++) { + for (int i = 0; i < static_cast<int>(macroElements.size()); i++) { if (macroElements[i]) { DELETE macroElements[i]; } } macroElements.resize(size); - for(i = 0; i < size; i++) { + for (int i = 0; i < size; i++) { macroElements[i] = NEW MacroElement(dim); macroElements[i]->writeNeighboursTo(&(neighbourIndices[i])); macroElements[i]->deserialize(in); @@ -971,9 +1011,9 @@ namespace AMDiS { in.read(reinterpret_cast<char*>(&initialized), sizeof(bool)); // set neighbour pointer in macro elements - int j, neighs = getGeo(NEIGH); - for(i = 0; i < static_cast<int>(macroElements.size()); i++) { - for(j = 0; j < neighs; j++) { + int neighs = getGeo(NEIGH); + for (int i = 0; i < static_cast<int>(macroElements.size()); i++) { + for (int j = 0; j < neighs; j++) { int index = neighbourIndices[i][j]; if(index != -1) { macroElements[i]->setNeighbour(j, macroElements[index]); @@ -986,7 +1026,7 @@ namespace AMDiS { // set mesh pointer in elements TraverseStack stack; ElInfo *elInfo = stack.traverseFirst(this, -1, CALL_EVERY_EL_PREORDER); - while(elInfo) { + while (elInfo) { elInfo->getElement()->setMesh(this); elInfo = stack.traverseNext(elInfo); } @@ -1057,8 +1097,7 @@ namespace AMDiS { void Mesh::clearMacroFileInfo() { - macroFileInfo->clear(getNumberOfEdges(), - getNumberOfVertices()); + macroFileInfo->clear(nEdges, nVertices); DELETE macroFileInfo; macroFileInfo = NULL; } diff --git a/AMDiS/src/Mesh.h b/AMDiS/src/Mesh.h index d107be226e2014ac1785367cd4a224c46011a286..d9d87abca22f904a2b9ee5f88887310d7afc9d8c 100644 --- a/AMDiS/src/Mesh.h +++ b/AMDiS/src/Mesh.h @@ -414,6 +414,13 @@ namespace AMDiS { /// Adds a macro element to the mesh void addMacroElement(MacroElement* me); + /* \brief + * Removes a macro element from the mesh. This works only for the case, that + * there are no global or local refinements, i.e., all macro elements have no + * children. + */ + void removeMacroElement(MacroElement* me); + /// Frees the array of DOF pointers (see \ref createDOFPtrs) void freeDOFPtrs(DegreeOfFreedom **ptrs); diff --git a/AMDiS/src/ParallelDomainProblem.cc b/AMDiS/src/ParallelDomainProblem.cc index 7ccc6de1833b6e58fb7fb7961c9229a8a441ebc2..ac9b5573d9353f4dd9465edb5896156a9438d413 100644 --- a/AMDiS/src/ParallelDomainProblem.cc +++ b/AMDiS/src/ParallelDomainProblem.cc @@ -12,10 +12,13 @@ namespace AMDiS { ParallelDomainProblemBase::ParallelDomainProblemBase(const std::string& name, - ProblemIterationInterface *iterationIF, - ProblemTimeInterface *timeIF, + ProblemIterationInterface *iIF, + ProblemTimeInterface *tIF, Mesh *m) - : mesh(m) + : iterationIF(iIF), + timeIF(tIF), + mesh(m), + nRankDOFs(0) { mpiRank = MPI::COMM_WORLD.Get_rank(); mpiSize = MPI::COMM_WORLD.Get_size(); @@ -23,6 +26,11 @@ namespace AMDiS { partitioner = new ParMetisPartitioner(mesh, &mpiComm); } + Flag ParallelDomainProblemBase::oneIteration(AdaptInfo *adaptInfo, Flag toDo) + { + return iterationIF->oneIteration(adaptInfo, toDo); + } + void ParallelDomainProblemBase::initParallelization(AdaptInfo *adaptInfo) { if (mpiSize <= 1) @@ -35,29 +43,99 @@ namespace AMDiS { // and now partition the mesh partitionMesh(adaptInfo); + std::vector<int> rankDofs; + std::map<int, std::set<int> > partitionDofs; TraverseStack stack; ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL); int nLeaves = 0; 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) { - element->setHidden(true); } else { - nLeaves++; } - + + // Determine to each dof the partition(s) it corresponds to. + for (int i = 0; i < 3; i++) + partitionDofs[element->getDOF(i, 0)].insert(partitionVec[element->getIndex()]); + elInfo = stack.traverseNext(elInfo); } - mesh->setNumberOfLeaves(nLeaves); + + for (std::map<int, 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 { + bool insert = true; + for (std::set<int>::iterator itpart2 = it->second.begin(); + itpart2 != it->second.end(); + ++itpart2) { + if (*itpart2 > mpiRank) { + insert = false; + break; + } + } + if (insert) { + rankDofs.push_back(it->first); + } + } + } + } + } + + bool removed = false; + do { + removed = false; + 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) { + mesh->removeMacroElement(*it); + removed = true; + break; + } + } + } while (removed); + + int *gOrder = (int*)(malloc(sizeof(int) * rankDofs.size())); + int *lOrder = (int*)(malloc(sizeof(int) * rankDofs.size())); + + for (std::vector<int>::iterator it = rankDofs.begin(); + it != rankDofs.end(); ++it) { + gOrder[nRankDOFs++] = *it; + } + + int rstart = 0; + MPI_Scan(&nRankDOFs, &rstart, 1, MPI_INT, MPI_SUM, PETSC_COMM_WORLD); + + for (int i = 0; i < nRankDOFs; i++) { + lOrder[i] = rstart - nRankDOFs + i; + } + + AOCreateBasic(PETSC_COMM_WORLD, nRankDOFs, gOrder, lOrder, &applicationOrdering); + + free(gOrder); + free(lOrder); } void ParallelDomainProblemBase::exitParallelization(AdaptInfo *adaptInfo) { + AODestroy(applicationOrdering); } double ParallelDomainProblemBase::setElemWeights(AdaptInfo *adaptInfo) diff --git a/AMDiS/src/ParallelDomainProblem.h b/AMDiS/src/ParallelDomainProblem.h index 40903a96e0bc3fc09079a485c2693f03e4ce20d7..77995799be3cb5ffeb6ba5f31af4c3eed3356440 100644 --- a/AMDiS/src/ParallelDomainProblem.h +++ b/AMDiS/src/ParallelDomainProblem.h @@ -28,6 +28,7 @@ #include "ProblemTimeInterface.h" #include "ProblemIterationInterface.h" #include "AdaptInfo.h" +#include "petscao.h" #include "mpi.h" namespace AMDiS { @@ -67,7 +68,15 @@ namespace AMDiS { virtual void transferInitialSolution(AdaptInfo *adaptInfo) {} - virtual Flag oneIteration(AdaptInfo *adaptInfo, Flag toDo = FULL_ITERATION) {} + virtual void beginIteration(AdaptInfo *adaptInfo) { + iterationIF->beginIteration(adaptInfo); + } + + virtual Flag oneIteration(AdaptInfo *adaptInfo, Flag toDo = FULL_ITERATION); + + virtual void endIteration(AdaptInfo *adaptInfo) { + iterationIF->endIteration(adaptInfo); + } virtual int getNumProblems() {} @@ -75,6 +84,16 @@ namespace AMDiS { return name; } + /// Returns pointer to \ref applicationOrdering. + AO* getApplicationOrdering() { + return &applicationOrdering; + } + + /// Returns \ref nRankDOFs, the number of DOFs in the rank mesh. + int getNumberRankDOFs() { + return nRankDOFs; + } + virtual ProblemStatBase *getProblem(int number = 0) {} virtual void serialize(std::ostream&) {} @@ -82,6 +101,12 @@ namespace AMDiS { virtual void deserialize(std::istream&) {} protected: + /// + ProblemIterationInterface *iterationIF; + + /// + ProblemTimeInterface *timeIF; + /// The rank of the current process. int mpiRank; @@ -121,6 +146,12 @@ namespace AMDiS { * of the parition it corresponds to is stored. */ std::map<int, int> oldPartitionVec; + + /// Petsc application ordering + AO applicationOrdering; + + /// Number of DOFs in the rank mesh. + int nRankDOFs; }; class ParallelDomainProblemScal : public ParallelDomainProblemBase diff --git a/AMDiS/src/ProblemIterationInterface.h b/AMDiS/src/ProblemIterationInterface.h index 742f9be822480e6a385396c30f6fbe58f7e543f6..f3a2c538739104583c1ec0c77638aef4f775b781 100644 --- a/AMDiS/src/ProblemIterationInterface.h +++ b/AMDiS/src/ProblemIterationInterface.h @@ -55,10 +55,8 @@ namespace AMDiS { public: virtual ~ProblemIterationInterface() {}; - /** \brief - * Called before each adaption loop iteration. - */ - virtual void beginIteration(AdaptInfo *adaptInfo) {}; + /// Called before each adaption loop iteration. + virtual void beginIteration(AdaptInfo *adaptInfo) {} /** \brief * Determines the execution order of the single adaption steps. If adapt is @@ -67,14 +65,10 @@ namespace AMDiS { */ virtual Flag oneIteration(AdaptInfo *adaptInfo, Flag toDo = FULL_ITERATION) = 0; - /** \brief - * Called after each adaption loop iteration. - */ - virtual void endIteration(AdaptInfo *adaptInfo) {}; + /// Called after each adaption loop iteration. + virtual void endIteration(AdaptInfo *adaptInfo) {} - /** \brief - * Returns number of managed problems - */ + /// Returns number of managed problems virtual int getNumProblems() = 0; /** \brief @@ -83,24 +77,18 @@ namespace AMDiS { */ virtual ProblemStatBase *getProblem(int number = 0) = 0; - /** \brief - * Returns the problem with the given name. - */ - virtual ProblemStatBase *getProblem(const std::string& name) { return NULL; }; + /// Returns the problem with the given name. + virtual ProblemStatBase *getProblem(const std::string& name) { + return NULL; + } - /** \brief - * Returns the name of the problem. - */ + /// Returns the name of the problem. virtual const std::string& getName() = 0; - /** \brief - * Function that serializes the problem plus information about the iteration. - */ + /// Function that serializes the problem plus information about the iteration. virtual void serialize(std::ostream &out) = 0; - /** \brief - * Function that deserializes the problem plus information about the iteration. - */ + /// Function that deserializes the problem plus information about the iteration. virtual void deserialize(std::istream &in) = 0; }; diff --git a/AMDiS/src/ProblemScal.cc b/AMDiS/src/ProblemScal.cc index b9fb370f3790c0b462026472e578c76f4abb2725..ba622f8d0a3d6af11f78ef142a727fa6c1bc438a 100644 --- a/AMDiS/src/ProblemScal.cc +++ b/AMDiS/src/ProblemScal.cc @@ -147,6 +147,17 @@ namespace AMDiS { rhs->getBoundaryManager()->addBoundaryCondition(periodic); } +#ifdef HAVE_PARALLEL_AMDIS + void ProblemScal::useApplicationOrdering(AO *ao, int nDofs) + { + applicationOrdering = ao; + nRankDOFs = nDofs; + systemMatrix->useApplicationOrdering(ao); + rhs->useApplicationOrdering(ao); + } +#endif + + void ProblemScal::createMesh() { TEST_EXIT(Parameters::initialized())("parameters not initialized\n"); @@ -169,16 +180,16 @@ namespace AMDiS { switch (dim) { case 1: - coarseningManager_ = NEW CoarseningManager1d(); - refinementManager_ = NEW RefinementManager1d(); + coarseningManager = NEW CoarseningManager1d(); + refinementManager = NEW RefinementManager1d(); break; case 2: - coarseningManager_ = NEW CoarseningManager2d(); - refinementManager_ = NEW RefinementManager2d(); + coarseningManager = NEW CoarseningManager2d(); + refinementManager = NEW RefinementManager2d(); break; case 3: - coarseningManager_ = NEW CoarseningManager3d(); - refinementManager_ = NEW RefinementManager3d(); + coarseningManager = NEW CoarseningManager3d(); + refinementManager = NEW RefinementManager3d(); break; default: ERROR_EXIT("invalid dim!\n"); @@ -189,8 +200,8 @@ namespace AMDiS { void ProblemScal::setMeshFromProblemVec(ProblemVec* pv, int i) { mesh = pv->getMesh(i); - coarseningManager_ = pv->getCoarseningManager(i); - refinementManager_ = pv->getRefinementManager(i); + coarseningManager = pv->getCoarseningManager(i); + refinementManager = pv->getRefinementManager(i); } @@ -207,14 +218,14 @@ namespace AMDiS { Flag ProblemScal::refineMesh(AdaptInfo *adaptInfo) { - return refinementManager_->refineMesh(mesh); + return refinementManager->refineMesh(mesh); } Flag ProblemScal::coarsenMesh(AdaptInfo *adaptInfo) { if (adaptInfo->isCoarseningAllowed(0)) { - return coarseningManager_->coarsenMesh(mesh); + return coarseningManager->coarsenMesh(mesh); } else { WARNING("coarsening not allowed\n"); return 0; @@ -273,8 +284,8 @@ namespace AMDiS { adoptFlag.isSet(INIT_FE_SPACE))) { TEST_EXIT(!mesh)("mesh already created\n"); mesh = adoptProblem->getMesh(); - refinementManager_ = adoptProblem->refinementManager_; - coarseningManager_ = adoptProblem->coarseningManager_; + refinementManager = adoptProblem->refinementManager; + coarseningManager = adoptProblem->coarseningManager; } } @@ -285,10 +296,12 @@ namespace AMDiS { if (feSpace_) { WARNING("feSpace already created\n"); } else { - if (initFlag.isSet(INIT_FE_SPACE) || (initFlag.isSet(INIT_SYSTEM) && !adoptFlag.isSet(INIT_FE_SPACE))) { + if (initFlag.isSet(INIT_FE_SPACE) || + (initFlag.isSet(INIT_SYSTEM) && !adoptFlag.isSet(INIT_FE_SPACE))) { createFESpace(); } - if (adoptProblem && (adoptFlag.isSet(INIT_FE_SPACE) || adoptFlag.isSet(INIT_SYSTEM))) { + if (adoptProblem && + (adoptFlag.isSet(INIT_FE_SPACE) || adoptFlag.isSet(INIT_SYSTEM))) { TEST_EXIT(!feSpace_)("feSpace already created"); feSpace_ = dynamic_cast<ProblemScal*>(adoptProblem)->getFESpace(); } @@ -412,7 +425,7 @@ namespace AMDiS { if (initFlag.isSet(INIT_GLOBAL_REFINES)) { int globalRefinements = 0; GET_PARAMETER(0, mesh->getName() + "->global refinements", "%d", &globalRefinements); - refinementManager_->globalRefine(mesh, globalRefinements); + refinementManager->globalRefine(mesh, globalRefinements); } } } @@ -576,7 +589,11 @@ namespace AMDiS { mesh->dofCompress(); MSG("%d DOFs for %s\n", +#ifdef HAVE_PARALLEL_AMDIS + nRankDOFs, +#else feSpace_->getAdmin()->getUsedSize(), +#endif feSpace_->getName().c_str()); Flag assembleFlag = @@ -603,20 +620,20 @@ namespace AMDiS { ElInfo *elInfo = stack.traverseFirst(mesh, -1, assembleFlag); while (elInfo) { - if (useGetBound) { - bound = GET_MEMORY(BoundaryType, feSpace_->getBasisFcts()->getNumber()); - feSpace_->getBasisFcts()->getBound(elInfo, bound); - } else { - bound = NULL; - } - - systemMatrix->assemble(1.0, elInfo, bound); - rhs->assemble(1.0, elInfo, bound); - - if (useGetBound) { - FREE_MEMORY(bound, BoundaryType, feSpace_->getBasisFcts()->getNumber()); - } - + if (useGetBound) { + bound = GET_MEMORY(BoundaryType, feSpace_->getBasisFcts()->getNumber()); + feSpace_->getBasisFcts()->getBound(elInfo, bound); + } else { + bound = NULL; + } + + systemMatrix->assemble(1.0, elInfo, bound); + rhs->assemble(1.0, elInfo, bound); + + if (useGetBound) { + FREE_MEMORY(bound, BoundaryType, feSpace_->getBasisFcts()->getNumber()); + } + if (systemMatrix->getBoundaryManager()) systemMatrix->getBoundaryManager()->fillBoundaryConditions(elInfo, systemMatrix); if (rhs->getBoundaryManager()) diff --git a/AMDiS/src/ProblemScal.h b/AMDiS/src/ProblemScal.h index 6cecdf75f6b48a8d9a6a827643901b861935c485..f574fdbf474b41397cb16523aded7342ba4004b8 100644 --- a/AMDiS/src/ProblemScal.h +++ b/AMDiS/src/ProblemScal.h @@ -22,6 +22,8 @@ #ifndef AMDIS_PROBLEMSCAL_H #define AMDIS_PROBLEMSCAL_H +#include <list> + #include "Global.h" #include "ProblemStatBase.h" #include "MemoryManager.h" @@ -30,7 +32,10 @@ #include "Boundary.h" #include "StandardProblemIteration.h" #include "ProblemVec.h" -#include <list> + +#ifdef HAVE_PARALLEL_AMDIS +#include "petscao.h" +#endif namespace AMDiS { @@ -58,7 +63,7 @@ namespace AMDiS { ProblemScal(const std::string &nameStr, ProblemIterationInterface *problemIteration = NULL) : StandardProblemIteration(this), - name(nameStr), + name(nameStr), feSpace_(NULL), mesh(NULL), marker(NULL), @@ -71,56 +76,42 @@ namespace AMDiS { leftPrecon(NULL), rightPrecon(NULL), useGetBound(true), - refinementManager_(NULL), - coarseningManager_(NULL), + refinementManager(NULL), + coarseningManager(NULL), +#ifdef HAVE_PARALLEL_AMDIS + applicationOrdering(NULL), + nRankDOFs(0), +#endif info_(10) {} - /** \brief - * Destructor - */ + /// Destructor virtual ~ProblemScal(); - /** \brief - * Initialisation of the problem. - */ + /// Initialisation of the problem. virtual void initialize(Flag initFlag, ProblemScal *adoptProblem = NULL, Flag adoptFlag = INIT_NOTHING); - /** \brief - * Used in \ref initialize(). - */ + /// Used in \ref initialize(). virtual void createMesh(); - /** \brief - * Used in \ref initialize(). - */ + /// Used in \ref initialize(). virtual void createFESpace(); - /** \brief - * Used in \ref initialize(). - */ + /// Used in \ref initialize(). virtual void createMatricesAndVectors(); - /** \brief - * Used in \ref initialize(). - */ + /// Used in \ref initialize(). virtual void createSolver(); - /** \brief - * Used in \ref initialize(). - */ + /// Used in \ref initialize(). virtual void createEstimator(); - /** \brief - * Used in \ref initialize(). - */ + /// Used in \ref initialize(). virtual void createMarker(); - /** \brief - * Used in \ref initialize(). - */ + /// Used in \ref initialize(). virtual void createFileWriter(); /** \brief @@ -173,16 +164,12 @@ namespace AMDiS { bool assembleMatrix = true, bool assembleVector = true); - /** \brief - * Returns number of managed problems - */ + /// Returns number of managed problems virtual int getNumProblems() { return 1; } - /** \brief - * Implementation of ProblemStatBase::getNumComponents() - */ + /// Implementation of ProblemStatBase::getNumComponents() virtual int getNumComponents() { return 1; } @@ -195,61 +182,50 @@ namespace AMDiS { return this; } - /** \brief - * Writes output files. - */ + /// Writes output files. void writeFiles(AdaptInfo *adaptInfo, bool force); - /** \brief - * Interpolates fct to \ref solution. - */ - void - interpolInitialSolution(AbstractFunction<double, WorldVector<double> > *fct); + /// Interpolates fct to \ref solution. + void interpolInitialSolution(AbstractFunction<double, WorldVector<double> > *fct); /// Adds an Operator to \ref systemMatrix. void addMatrixOperator(Operator *op, - double *factor = NULL, - double *estFactor = NULL); + double *factor = NULL, double *estFactor = NULL); /// Adds an Operator to \ref rhs. void addVectorOperator(Operator *op, - double *factor = NULL, - double *estFactor = NULL); + double *factor = NULL, double *estFactor = NULL); - /** \brief - * Adds dirichlet boundary conditions. - */ + /// Adds dirichlet boundary conditions. virtual void addDirichletBC(BoundaryType type, AbstractFunction<double, WorldVector<double> >* b); - /** \brief - * Adds dirichlet boundary conditions. - */ + /// Adds dirichlet boundary conditions. virtual void addDirichletBC(BoundaryType type, DOFVector<double> *vec); - /** \brief - * Adds robin boundary conditions. - */ + /// Adds robin boundary conditions. virtual void addRobinBC(BoundaryType type, AbstractFunction<double, WorldVector<double> > *n, AbstractFunction<double, WorldVector<double> > *r); - /** \brief - * Adds robin boundary conditions. - */ + /// Adds robin boundary conditions. virtual void addRobinBC(BoundaryType type, DOFVector<double> *n, DOFVector<double> *r); - /** \brief - * Adds Neumann boundary conditions. - */ + /// Adds Neumann boundary conditions. virtual void addNeumannBC(BoundaryType type, AbstractFunction<double, WorldVector<double> > *n); + /// Adds periodic boundary conditions. virtual void addPeriodicBC(BoundaryType type); + +#ifdef HAVE_PARALLEL_AMDIS + /// Sets the petsc application ordering object to map dof indices. + void useApplicationOrdering(AO *ao, int nDofs); +#endif // ===== getting-methods ====================================================== @@ -257,53 +233,39 @@ namespace AMDiS { * \{ */ - /** \brief - * Returns \ref solution. - */ + /// Returns \ref solution. inline DOFVector<double>* getSolution() { return solution; } - /** \brief - * Returns \ref rhs. - */ + /// Returns \ref rhs. inline DOFVector<double>* getRHS() { return rhs; } - /** \brief - * Returns \ref systemMatrix. - */ + /// Returns \ref systemMatrix. inline DOFMatrix *getSystemMatrix() { return systemMatrix; } - /** \brief - * Returns \ref mesh - */ + /// Returns \ref mesh inline Mesh* getMesh() { return mesh; } - /** \brief - * Returns \ref feSpace_. - */ + /// Returns \ref feSpace_. inline FiniteElemSpace* getFESpace() { return feSpace_; } - /** \brief - * Returns \ref estimator_. - */ + /// Returns \ref estimator_. inline Estimator* getEstimator() { return estimator; } - /** \brief - * Returns \ref refinementManager_. - */ + /// Returns \ref refinementManager. inline RefinementManager* getRefinementManager() { - return refinementManager_; + return refinementManager; } /// Returns \ref solver. @@ -338,7 +300,7 @@ namespace AMDiS { /// inline CoarseningManager *getCoarseningManager() { - return coarseningManager_; + return coarseningManager; } /** \} */ @@ -349,47 +311,35 @@ namespace AMDiS { * \{ */ - /** \brief - * Sets \ref mesh - */ + /// Sets \ref mesh void setMeshFromProblemVec(ProblemVec* pv, int i = 0); - /** \brief - * Sets \ref feSpace_. - */ + /// Sets \ref feSpace_. inline void setFESpace(FiniteElemSpace* fe) { feSpace_ = fe; } - /** \brief - * Sets \ref estimator_. - */ + /// Sets \ref estimator_. inline void setEstimator(Estimator* est) { estimator = est; } - /** \brief - * Sets \ref solver. - */ + /// Sets \ref solver. inline void setSolver(OEMSolver<DOFVector<double> >* sol) { solver = sol; } - /** \brief - * Sets \ref leftPrecon. - */ + /// Sets \ref leftPrecon. inline void setLeftPrecon(Preconditioner<DOFVector<double> > *p) { leftPrecon = p; } - /** \brief - * Sets \ref rightPrecon. - */ + /// Sets \ref rightPrecon. inline void setRightPrecon(Preconditioner<DOFVector<double> > *p) { rightPrecon = p; } - + /// Sets \ref marker. inline void setMarker(Marker *m) { marker = m; } @@ -398,16 +348,10 @@ namespace AMDiS { void writeResidualMesh(AdaptInfo *adaptInfo, const std::string name); - // ===== Serializable implementation ===== - - /** \brief - * serialization - */ + /// serialization virtual void serialize(std::ostream &out); - /** \brief - * deserialization - */ + /// deserialization virtual void deserialize(std::istream &in); std::vector<FileWriterInterface*> getFileWriterList() { @@ -462,17 +406,26 @@ namespace AMDiS { * If new refinement algorithms should be realized, one has to override * RefinementManager and give one instance of it to AdaptStationary. */ - RefinementManager *refinementManager_; + RefinementManager *refinementManager; /** \brief * All actions of mesh coarsening are performed by coarseningManager. * If new coarsening algorithms should be realized, one has to override * CoarseningManager and give one instance of it to AdaptStationary. */ - CoarseningManager *coarseningManager_; + CoarseningManager *coarseningManager; +#ifdef HAVE_PARALLEL_AMDIS + /// Petsc application ordering to map dof indices. + AO *applicationOrdering; + + /// Number of DOFs in the rank partition. + int nRankDOFs; +#endif + /// Info level. int info_; + }; } diff --git a/AMDiS/src/RefinementManager.h b/AMDiS/src/RefinementManager.h index 23a9e843db2aebc5a3fb68156fa7757719c81ba9..69fab8829ea750c708cf451657397bbe8dde6d10 100644 --- a/AMDiS/src/RefinementManager.h +++ b/AMDiS/src/RefinementManager.h @@ -22,19 +22,11 @@ #ifndef AMDIS_REFINEMENTMANAGER_H #define AMDIS_REFINEMENTMANAGER_H -// ============================================================================ -// ===== includes ============================================================= -// ============================================================================ - #include "Global.h" #include "Flag.h" namespace AMDiS { - // ============================================================================ - // ===== forward declarations ================================================= - // ============================================================================ - class Element; class Line; class Triangle; @@ -44,9 +36,6 @@ namespace AMDiS { class TraverseStack; class RCNeighbourList; - // ============================================================================ - // ===== class RefinementManager ============================================== - // ============================================================================ /** \ingroup Adaption * \brief @@ -57,18 +46,15 @@ namespace AMDiS { class RefinementManager { public: - /** \brief - * Constructs a RefinementManager which belongs to aMesh - */ + /// Constructs a RefinementManager which belongs to aMesh RefinementManager() - : mesh(NULL), newCoords(false), + : mesh(NULL), + newCoords(false), stack(NULL) - {}; + {} - /** \brief - * destructor - */ - virtual ~RefinementManager() {}; + /// Destructor + virtual ~RefinementManager() {} /** \brief * Generates new coordinates on curved boundaries. Can be overriden by @@ -77,7 +63,7 @@ namespace AMDiS { virtual void setNewCoords() { FUNCNAME("RefinementManager::setNewCoords"); ERROR_EXIT("called for base class!\n"); - }; + } /** \brief * Fulfills the refinement for all positive marked elements of the mesh. @@ -88,21 +74,17 @@ namespace AMDiS { */ virtual Flag refineMesh(Mesh *aMesh); - /** \brief - * All elements of the mesh will be refined. - */ + /// All elements of the mesh will be refined. Flag globalRefine(Mesh *aMesh, int mark); - /** \brief - * set \ref newCoords - */ + /// Set \ref newCoords inline void newCoord(bool nc) { newCoords = nc; - }; + } inline bool newCoord() { return newCoords; - }; + } /** \brief * Implements the refinement of el_info->el. Can be overriden by sub @@ -112,51 +94,39 @@ namespace AMDiS { FUNCNAME("RefinementManager::refineFunction()"); ERROR_EXIT("called for base class!\n"); return NULL; - }; + } inline void setMesh(Mesh *m) { mesh = m; - }; + } inline void setStack(TraverseStack *s) { stack = s; - }; + } inline TraverseStack *getStack() { return stack; - }; + } protected: - /** \brief - * Used by \ref globalRefine - */ + /// Used by \ref globalRefine static int globalRefineFct(ElInfo* elinfo); protected: - /** \brief - * The Mesh to be refined - */ + /// The Mesh to be refined Mesh *mesh; - /** \brief - * Number of new vertices on a boundary edge - */ + /// Number of new vertices on a boundary edge bool newCoords; - /** \brief - * Still more refinement to do? - */ + /// Still more refinement to do? static bool doMoreRecursiveRefine; - /** \brief - * number of DOFVectors which must be interpolated during refinement - */ + /// Number of DOFVectors which must be interpolated during refinement static int callRefineInterpol; - /** \brief - * Used for non recursive traversal - */ + /// Used for non recursive traversal TraverseStack* stack; /** \brief @@ -165,9 +135,7 @@ namespace AMDiS { */ static RefinementManager* traversePtr; - /** \brief - * Used by globalRefine to remember the given mark value - */ + /// Used by globalRefine to remember the given mark value static int globalMark; }; diff --git a/AMDiS/src/RefinementManager2d.h b/AMDiS/src/RefinementManager2d.h index fd36e4187df3f4879a323c246ca536f24ec94353..ab541dfcb658c1d14b6f287ce51f7cfdc940cbc0 100644 --- a/AMDiS/src/RefinementManager2d.h +++ b/AMDiS/src/RefinementManager2d.h @@ -26,10 +26,6 @@ namespace AMDiS { - // ============================================================================ - // ===== class RefinementManager2d ============================================ - // ============================================================================ - /** \ingroup Adaption * \brief * Implements a RefinementManager for 2-dimensional meshes. @@ -39,27 +35,19 @@ namespace AMDiS { public: MEMORY_MANAGED(RefinementManager2d); - /** \brief - * Calls base class constructor. - */ + /// Calls base class constructor. RefinementManager2d() : RefinementManager() - {}; + {} - /** \brief - * destructor - */ - virtual ~RefinementManager2d() {}; + /// Destructor + virtual ~RefinementManager2d() {} protected: - /** \brief - * Used by \ref setNewCoords - */ + /// Used by \ref setNewCoords static int newCoordsFct(ElInfo *el_info); - /** \brief - * Implements RefinementManager::setNewCoords - */ + /// Implements RefinementManager::setNewCoords void setNewCoords(); /** \brief @@ -73,20 +61,14 @@ namespace AMDiS { int getRefinePatch(ElInfo **el_info, DegreeOfFreedom *edge[2], int dir, RCNeighbourList* refineList, int *n_neigh); - /** \brief - * Refines all elements in the patch. - */ + /// Refines all elements in the patch. DegreeOfFreedom refinePatch(DegreeOfFreedom *edge[2], RCNeighbourList* refineList, int n_neigh, bool bound); - /** \brief - * Implements RefinementManager::refineFunction. - */ + /// Implements RefinementManager::refineFunction. ElInfo* refineFunction(ElInfo* el_info); - /** \brief - * Refines one Triangle. - */ + /// Refines one Triangle. void bisectTriangle(Triangle *el, DegreeOfFreedom* newDofs[3]); }; diff --git a/AMDiS/src/Traverse.cc b/AMDiS/src/Traverse.cc index 40b7e74e16304a79bc3f4826a369f076195c5212..b767e0aed517a06cc872de19ac7f8bb91ecbaabc 100644 --- a/AMDiS/src/Traverse.cc +++ b/AMDiS/src/Traverse.cc @@ -132,11 +132,10 @@ namespace AMDiS { if (stack_used == 0) { /* first call */ currentMacro = traverse_mesh->firstMacroElement(); - // Search for the first macro element that corresponds to the current thread - // and that is not hidden. So increment as long as one of this conditions does - // not hold and the current macro element is not the last one. - while (((*currentMacro)->getIndex() % maxThreads_ != myThreadId_ || - (*currentMacro)->getElement()->isHidden()) && + // Search for the first macro element that corresponds to the current thread. + // So increment as long as one of this condition does not hold and the + // current macro element is not the last one. + while ((*currentMacro)->getIndex() % maxThreads_ != myThreadId_ && currentMacro != traverse_mesh->endOfMacroElements()) { currentMacro++; } @@ -166,13 +165,12 @@ namespace AMDiS { /* goto next macro element */ if (stack_used < 1) { - // Again, search for the next macro element that corresponds to the current - // thread and that is not hidden. + // Again, search for the next macro element that corresponds to the + // current thread. do { currentMacro++; - } while (currentMacro != traverse_mesh->endOfMacroElements() && - ((*currentMacro)->getIndex() % maxThreads_ != myThreadId_ || - (*currentMacro)->getElement()->isHidden())); + } while (currentMacro != traverse_mesh->endOfMacroElements() && + (*currentMacro)->getIndex() % maxThreads_ != myThreadId_); if (currentMacro == traverse_mesh->endOfMacroElements()) { return NULL; @@ -198,49 +196,14 @@ namespace AMDiS { int visitedChildren = info_stack[stack_used]; if (visitedChildren == 0) { - // No children of this element have been visited. - - if (el->getFirstChild()->isHidden()) { - // The first child is hidden, so omit it. - info_stack[stack_used]++; - visitedChildren++; - - // And check the second child. - if (el->getSecondChild()->isHidden()) { - // Also the second child is hidden, so omit it too and start the - // traverse procedure again to go up into the stack. - - info_stack[stack_used]++; - - return traverseLeafElement(); - } else { - // Only the first child is hidden, so go down on the second one. - el = const_cast<Element*>(el->getSecondChild()); - } - } else { - // The first child is not hidden, so go down on it. - el = const_cast<Element*>(el->getFirstChild()); - } - + el = const_cast<Element*>(el->getFirstChild()); } else { - // The first child has been visited before, so check the second one. - if (el->getSecondChild()->isHidden()) { - // The second one is hidden, so omit it and start the traverse - // procedure again to go up into the stack; - info_stack[stack_used]++; - - return traverseLeafElement(); - } else { - // The secod child is not hidden, so go down on it. - el = const_cast<Element*>(el->getSecondChild()); - } + el = const_cast<Element*>(el->getSecondChild()); } - info_stack[stack_used]++; elinfo_stack[stack_used + 1]->fillElInfo(visitedChildren, elinfo_stack[stack_used]); - info_stack[stack_used + 1] = 0; stack_used++; TEST_EXIT_DBG(stack_used < stack_size)