diff --git a/AMDiS/src/DOFVector.hh b/AMDiS/src/DOFVector.hh index 141a3dc6ea91ff146138f0dad101f7acffa5d793..21c9357f56dab7f556ac8d30caec6360b6dd19e1 100644 --- a/AMDiS/src/DOFVector.hh +++ b/AMDiS/src/DOFVector.hh @@ -378,9 +378,9 @@ namespace AMDiS { template void DOFVector::interpol(AbstractFunction > *fct) { - FUNCNAME("interpol"); + FUNCNAME("DOFVector::interpol()"); - TEST_EXIT_DBG(fct)("no function to interpolate\n"); + TEST_EXIT_DBG(fct)("No function to interpolate!\n"); interFct = fct; diff --git a/AMDiS/src/Operator.h b/AMDiS/src/Operator.h index 00eb5f39e9bf5aefe9dd08e621fbc166e6c61e30..821767500eaa5d625986a19af3db19affdf7401a 100644 --- a/AMDiS/src/Operator.h +++ b/AMDiS/src/Operator.h @@ -570,9 +570,7 @@ namespace AMDiS { lalt_kl(Lambda, xi, xj, *(LALt[iq]), (*factor)); } - /** \brief - * Implenetation of SecondOrderTerm::eval(). - */ + /// Implenetation of SecondOrderTerm::eval(). void eval(int nPoints, const double *, const WorldVector *, @@ -585,9 +583,7 @@ namespace AMDiS { result[iq] += (*factor) * D2UhAtQP[iq][xi][xj] * fac; } - /** \brief - * Implenetation of SecondOrderTerm::weakEval(). - */ + /// Implenetation of SecondOrderTerm::weakEval(). void weakEval(int nPoints, const WorldVector *grdUhAtQP, WorldVector *result) const @@ -598,14 +594,10 @@ namespace AMDiS { } private: - /** \brief - * Directions for the derivatives. - */ + /// Directions for the derivatives. int xi, xj; - /** \brief - * Pointer to the factor. - */ + /// Pointer to the factor. double *factor; }; @@ -722,20 +714,14 @@ namespace AMDiS { setSymmetric(true); } - /** \brief - * Implementation of \ref OperatorTerm::initElement(). - */ + /// Implementation of \ref OperatorTerm::initElement(). void initElement(const ElInfo* elInfo, SubAssembler* subAssembler, Quadrature *quad = NULL); - /** \brief - * Implements SecondOrderTerm::getLALt(). - */ + /// Implements SecondOrderTerm::getLALt(). void getLALt(const ElInfo *elInfo, int nPoints, DimMat **LALt) const; - /** \brief - * Implenetation of SecondOrderTerm::eval(). - */ + /// Implenetation of SecondOrderTerm::eval(). void eval(int nPoints, const double *uhAtQP, const WorldVector *grdUhAtQP, @@ -743,23 +729,17 @@ namespace AMDiS { double *result, double factor) const; - /** \brief - * Implenetation of SecondOrderTerm::weakEval(). - */ + /// Implenetation of SecondOrderTerm::weakEval(). void weakEval(int nPoints, const WorldVector *grdUhAtQP, WorldVector *result) const; protected: - /** \brief - * Stores coordinates at quadrature points. Set in \ref initElement(). - */ + /// Stores coordinates at quadrature points. Set in \ref initElement(). WorldVector* coordsAtQPs; - /** \brief - * Function evaluated at quadrature points. - */ + /// Function evaluated at quadrature points. AbstractFunction > *g; }; @@ -2522,20 +2502,14 @@ namespace AMDiS { VecGradCoordsAtQP_ZOT(DOFVectorBase *dv, TertiaryAbstractFunction, WorldVector > *f); - /** \brief - * Implementation of \ref OperatorTerm::initElement(). - */ + /// Implementation of \ref OperatorTerm::initElement(). void initElement(const ElInfo* elInfo, SubAssembler* subAssembler, Quadrature *quad = NULL); - /** \brief - * Implements ZeroOrderTerm::getC(). - */ + /// Implements ZeroOrderTerm::getC(). void getC(const ElInfo *, int nPoints, std::vector &C) const; - /** \brief - * Implements ZeroOrderTerm::eval(). - */ + /// Implements ZeroOrderTerm::eval(). void eval(int nPoints, const double *uhAtQP, const WorldVector *grdUhAtQP, @@ -2544,26 +2518,18 @@ namespace AMDiS { double fac) const; protected: - /** \brief - * DOFVector to be evaluated at quadrature points. - */ + /// DOFVector to be evaluated at quadrature points. DOFVectorBase* vec; - /** \brief - * Vector v at quadrature points. - */ + /// Vector v at quadrature points. double *vecAtQPs; - /** \brief - * Gradient at quadrature points. - */ + /// Gradient at quadrature points. WorldVector* gradAtQPs; WorldVector* coordsAtQPs; - /** \brief - * Function for c. - */ + /// Function for c. TertiaryAbstractFunction, WorldVector > *f; }; @@ -2579,20 +2545,14 @@ namespace AMDiS { VecAndCoordsAtQP_ZOT(DOFVectorBase *dv, BinaryAbstractFunction > *f); - /** \brief - * Implementation of \ref OperatorTerm::initElement(). - */ + /// Implementation of \ref OperatorTerm::initElement(). void initElement(const ElInfo* elInfo, SubAssembler* subAssembler, Quadrature *quad = NULL); - /** \brief - * Implements ZeroOrderTerm::getC(). - */ + /// Implements ZeroOrderTerm::getC(). void getC(const ElInfo *, int nPoints, std::vector &C) const; - /** \brief - * Implements ZeroOrderTerm::eval(). - */ + /// Implements ZeroOrderTerm::eval(). void eval(int nPoints, const double *uhAtQP, const WorldVector *grdUhAtQP, @@ -2601,24 +2561,16 @@ namespace AMDiS { double fac) const; protected: - /** \brief - * DOFVector to be evaluated at quadrature points. - */ + /// DOFVector to be evaluated at quadrature points. DOFVectorBase* vec; - /** \brief - * Vector v at quadrature points. - */ + /// Vector v at quadrature points. double *vecAtQPs; - /** \brief - * Gradient at quadrature points. - */ + /// Gradient at quadrature points. WorldVector* coordsAtQPs; - /** \brief - * Function for c. - */ + /// Function for c. BinaryAbstractFunction > *f; }; @@ -2635,20 +2587,14 @@ namespace AMDiS { Vec2AndGradAtQP_ZOT(DOFVectorBase *dv1, DOFVectorBase *dv2, TertiaryAbstractFunction, double > *af); - /** \brief - * Implementation of \ref OperatorTerm::initElement(). - */ + /// Implementation of \ref OperatorTerm::initElement(). void initElement(const ElInfo* elInfo, SubAssembler* subAssembler, Quadrature *quad = NULL); - /** \brief - * Implements ZeroOrderTerm::getC(). - */ + /// Implements ZeroOrderTerm::getC(). void getC(const ElInfo *, int nPoints, std::vector &C) const; - /** \brief - * Implements ZeroOrderTerm::eval(). - */ + /// Implements ZeroOrderTerm::eval(). void eval(int nPoints, const double *uhAtQP, const WorldVector *grdUhAtQP, @@ -2657,26 +2603,18 @@ namespace AMDiS { double fac) const; protected: - /** \brief - * DOFVector to be evaluated at quadrature points. - */ + /// DOFVector to be evaluated at quadrature points. DOFVectorBase* vec1; DOFVectorBase* vec2; - /** \brief - * Vector v at quadrature points. - */ + /// Vector v at quadrature points. double *vecAtQPs1; double *vecAtQPs2; - /** \brief - * Gradient at quadrature points. - */ + /// Gradient at quadrature points. WorldVector *gradAtQPs; - /** \brief - * Function for c. - */ + /// Function for c. TertiaryAbstractFunction, double > *f; }; @@ -2694,20 +2632,14 @@ namespace AMDiS { FctGradient_ZOT(DOFVectorBase *dv, AbstractFunction > *f); - /** \brief - * Implementation of \ref OperatorTerm::initElement(). - */ + /// Implementation of \ref OperatorTerm::initElement(). void initElement(const ElInfo* elInfo, SubAssembler* subAssembler, Quadrature *quad = NULL); - /** \brief - * Implements SecondOrderTerm::getC(). - */ + /// Implements SecondOrderTerm::getC(). void getC(const ElInfo *elInfo, int nPoints, std::vector &C) const; - /** \brief - * Implements ZeroOrderTerm::eval(). - */ + /// Implements ZeroOrderTerm::eval(). void eval(int nPoints, const double *uhAtQP, const WorldVector *grdUhAtQP, @@ -2718,9 +2650,7 @@ namespace AMDiS { protected: DOFVectorBase* vec; - /** \brief - * Function wich maps \ref gradAtQPs to a double. - */ + /// Function wich maps \ref gradAtQPs to a double. AbstractFunction > *f; /** \brief @@ -2743,20 +2673,14 @@ namespace AMDiS { VecAndGradAtQP_ZOT(DOFVectorBase *dv, BinaryAbstractFunction > *f); - /** \brief - * Implementation of \ref OperatorTerm::initElement(). - */ + /// Implementation of \ref OperatorTerm::initElement(). void initElement(const ElInfo* elInfo, SubAssembler* subAssembler, Quadrature *quad = NULL); - /** \brief - * Implements ZeroOrderTerm::getC(). - */ + /// Implements ZeroOrderTerm::getC(). void getC(const ElInfo *, int nPoints, std::vector &C) const; - /** \brief - * Implements ZeroOrderTerm::eval(). - */ + /// Implements ZeroOrderTerm::eval(). void eval(int nPoints, const double *uhAtQP, const WorldVector *grdUhAtQP, @@ -2765,24 +2689,16 @@ namespace AMDiS { double fac) const; protected: - /** \brief - * DOFVector to be evaluated at quadrature points. - */ + /// DOFVector to be evaluated at quadrature points. DOFVectorBase* vec; - /** \brief - * Vector v at quadrature points. - */ + /// Vector v at quadrature points. double *vecAtQPs; - /** \brief - * Gradient at quadrature points. - */ + /// Gradient at quadrature points. WorldVector *gradAtQPs; - /** \brief - * Function for c. - */ + /// Function for c. BinaryAbstractFunction > *f; }; @@ -2798,21 +2714,15 @@ namespace AMDiS { VecAndGradAtQP_SOT(DOFVectorBase *dv, BinaryAbstractFunction > *af); - /** \brief - * Implementation of \ref OperatorTerm::initElement(). - */ + /// Implementation of \ref OperatorTerm::initElement(). void initElement(const ElInfo* elInfo, SubAssembler* subAssembler, Quadrature *quad = NULL); - /** \brief - * Implements SecondOrderTerm::getLALt(). - */ + /// Implements SecondOrderTerm::getLALt(). inline void getLALt(const ElInfo *elInfo, int nPoints, DimMat **LALt) const; - /** \brief - * Implements SecondOrderTerm::eval(). - */ + /// Implements SecondOrderTerm::eval(). void eval(int nPoints, const double *uhAtQP, const WorldVector *grdUhAtQP, @@ -2820,32 +2730,22 @@ namespace AMDiS { double *result, double factor) const; - /** \brief - * Implements SecondOrderTerm::weakEval(). - */ + /// Implements SecondOrderTerm::weakEval(). void weakEval(int nPoints, const WorldVector *grdUhAtQP, WorldVector *result) const; protected: - /** \brief - * DOFVector to be evaluated at quadrature points. - */ + /// DOFVector to be evaluated at quadrature points. DOFVectorBase* vec; - /** \brief - * Vector v at quadrature points. - */ + /// Vector v at quadrature points. double *vecAtQPs; - /** \brief - * Gradient at quadrature points. - */ + /// Gradient at quadrature points. WorldVector *gradAtQPs; - /** \brief - * Function for c. - */ + /// Function for c. BinaryAbstractFunction > *f; }; @@ -2916,21 +2816,15 @@ namespace AMDiS { : ZeroOrderTerm(af->getDegree()), g(af) {} - /** \brief - * Implementation of \ref OperatorTerm::initElement(). - */ + /// Implementation of \ref OperatorTerm::initElement(). void initElement(const ElInfo* elInfo, SubAssembler* subAssembler, Quadrature *quad = NULL); - /** \brief - * Implements ZeroOrderTerm::getC(). - */ + /// Implements ZeroOrderTerm::getC(). void getC(const ElInfo *, int nPoints, std::vector &C) const; - /** \brief - * Implements ZeroOrderTerm::eval(). - */ + /// Implements ZeroOrderTerm::eval(). void eval(int nPoints, const double *uhAtQP, const WorldVector *grdUhAtQP, @@ -2939,14 +2833,10 @@ namespace AMDiS { double fac) const; protected: - /** \brief - * Stores coordinates at quadrature points. Set in \ref initElement(). - */ + /// Stores coordinates at quadrature points. Set in \ref initElement(). WorldVector* coordsAtQPs; - /** \brief - * Function for c. - */ + /// Function for c. AbstractFunction > *g; }; @@ -2969,20 +2859,14 @@ namespace AMDiS { setSymmetric(xi == xj); } - /** \brief - * Implementation of \ref OperatorTerm::initElement(). - */ + /// Implementation of \ref OperatorTerm::initElement(). void initElement(const ElInfo* elInfo, SubAssembler* subAssembler, Quadrature *quad = NULL); - /** \brief - * Implements SecondOrderTerm::getLALt(). - */ + /// Implements SecondOrderTerm::getLALt(). inline void getLALt(const ElInfo *elInfo, int nPoints, DimMat **LALt) const; - /** \brief - * Implements SecondOrderTerm::eval(). - */ + /// Implements SecondOrderTerm::eval(). void eval(int nPoints, const double *uhAtQP, const WorldVector *grdUhAtQP, @@ -2990,27 +2874,19 @@ namespace AMDiS { double *result, double factor) const; - /** \brief - * Implements SecondOrderTerm::weakEval(). - */ + /// Implements SecondOrderTerm::weakEval(). void weakEval(int nPoints, const WorldVector *grdUhAtQP, WorldVector *result) const; private: - /** \brief - * Stores coordinates at quadrature points. Set in \ref initElement(). - */ + /// Stores coordinates at quadrature points. Set in \ref initElement(). WorldVector* coordsAtQPs; - /** \brief - * Function evaluated at quadrature points. - */ + /// Function evaluated at quadrature points. AbstractFunction > *g; - /** \brief - * Directions for the derivatives. - */ + /// Directions for the derivatives. int xi, xj; }; @@ -3030,20 +2906,14 @@ namespace AMDiS { VecAtQP_IJ_SOT(DOFVectorBase *dv, AbstractFunction *af, int x_i, int x_j); - /** \brief - * Implementation of \ref OperatorTerm::initElement(). - */ + /// Implementation of \ref OperatorTerm::initElement(). void initElement(const ElInfo* elInfo, SubAssembler* subAssembler, Quadrature *quad = NULL); - /** \brief - * Implements SecondOrderTerm::getLALt(). - */ + /// Implements SecondOrderTerm::getLALt(). void getLALt(const ElInfo *elInfo, int nPoints, DimMat **LALt) const; - /** \brief - * Implements SecondOrderTerm::eval(). - */ + /// Implements SecondOrderTerm::eval(). void eval(int nPoints, const double *uhAtQP, const WorldVector *grdUhAtQP, @@ -3051,34 +2921,23 @@ namespace AMDiS { double *result, double factor) const; - /** \brief - * Implements SecondOrderTerm::weakEval(). - */ + /// Implements SecondOrderTerm::weakEval(). void weakEval(int nPoints, const WorldVector *grdUhAtQP, WorldVector *result) const; protected: - /** \brief - * DOFVector to be evaluated at quadrature points. - */ + /// DOFVector to be evaluated at quadrature points. DOFVectorBase* vec; - /** \brief - * Pointer to an array containing the DOFVector evaluated at quadrature - * points. - */ + /// Pointer to an array containing the DOFVector evaluated at quadrature points. double* vecAtQPs; - /** \brief - * Function evaluated at \ref vecAtQPs. - */ + /// Function evaluated at \ref vecAtQPs. AbstractFunction *f; private: - /** \brief - * Directions for the derivatives. - */ + /// Directions for the derivatives. int xi, xj; }; @@ -3090,20 +2949,14 @@ namespace AMDiS { DOFVectorBase *dGrd, BinaryAbstractFunction > *f); - /** \brief - * Implementation of \ref OperatorTerm::initElement(). - */ + /// Implementation of \ref OperatorTerm::initElement(). void initElement(const ElInfo* elInfo, SubAssembler* subAssembler, Quadrature *quad = NULL); - /** \brief - * Implements ZeroOrderTerm::getC(). - */ + /// Implements ZeroOrderTerm::getC(). void getC(const ElInfo *, int nPoints, std::vector &C) const; - /** \brief - * Implements ZeroOrderTerm::eval(). - */ + /// Implements ZeroOrderTerm::eval(). void eval(int nPoints, const double *uhAtQP, const WorldVector *grdUhAtQP, @@ -3112,29 +2965,19 @@ namespace AMDiS { double fac) const; protected: - /** \brief - * DOFVector to be evaluated at quadrature points. - */ + /// DOFVector to be evaluated at quadrature points. DOFVectorBase* vec; - /** \brief - * Vector v at quadrature points. - */ + /// Vector v at quadrature points. double *vecAtQPs; - /** \brief - * First DOFVector whose gradient is evaluated at quadrature points. - */ + /// First DOFVector whose gradient is evaluated at quadrature points. DOFVectorBase* vecGrd; - /** \brief - * Gradient of first vector at quadrature points. - */ + /// Gradient of first vector at quadrature points. WorldVector *gradAtQPs; - /** \brief - * Function for c. - */ + /// Function for c. BinaryAbstractFunction > *f; }; @@ -3147,20 +2990,14 @@ namespace AMDiS { DOFVectorBase *dGrd2, TertiaryAbstractFunction, WorldVector > *af); - /** \brief - * Implementation of \ref OperatorTerm::initElement(). - */ + /// Implementation of \ref OperatorTerm::initElement(). void initElement(const ElInfo* elInfo, SubAssembler* subAssembler, Quadrature *quad = NULL); - /** \brief - * Implements ZeroOrderTerm::getC(). - */ + /// Implements ZeroOrderTerm::getC(). void getC(const ElInfo *, int nPoints, std::vector &C) const; - /** \brief - * Implements ZeroOrderTerm::eval(). - */ + /// Implements ZeroOrderTerm::eval(). void eval(int nPoints, const double *uhAtQP, const WorldVector *grdUhAtQP, @@ -3169,39 +3006,25 @@ namespace AMDiS { double fac) const; protected: - /** \brief - * DOFVector to be evaluated at quadrature points. - */ + /// DOFVector to be evaluated at quadrature points. DOFVectorBase* vec; - /** \brief - * Vector v at quadrature points. - */ + /// Vector v at quadrature points. double *vecAtQPs; - /** \brief - * First DOFVector whose gradient is evaluated at quadrature points. - */ + /// First DOFVector whose gradient is evaluated at quadrature points. DOFVectorBase* vecGrd1; - /** \brief - * Gradient of first vector at quadrature points. - */ + /// Gradient of first vector at quadrature points. WorldVector *grad1AtQPs; - /** \brief - * Second DOFVector whose gradient is evaluated at quadrature points. - */ + /// Second DOFVector whose gradient is evaluated at quadrature points. DOFVectorBase* vecGrd2; - /** \brief - * Gradient of second vector at quadrature points. - */ + /// Gradient of second vector at quadrature points. WorldVector *grad2AtQPs; - /** \brief - * Function for c. - */ + /// Function for c. TertiaryAbstractFunction, WorldVector > *f; }; @@ -3213,20 +3036,14 @@ namespace AMDiS { VecOfDOFVecsAtQP_ZOT(const std::vector*>& dv, AbstractFunction > *f); - /** \brief - * Implementation of \ref OperatorTerm::initElement(). - */ + /// Implementation of \ref OperatorTerm::initElement(). void initElement(const ElInfo* elInfo, SubAssembler* subAssembler, Quadrature *quad = NULL); - /** \brief - * Implements ZeroOrderTerm::getC(). - */ + /// Implements ZeroOrderTerm::getC(). void getC(const ElInfo *, int nPoints, std::vector &C) const; - /** \brief - * Implements ZeroOrderTerm::eval(). - */ + /// Implements ZeroOrderTerm::eval(). void eval(int nPoints, const double *uhAtQP, const WorldVector *grdUhAtQP, @@ -3235,19 +3052,13 @@ namespace AMDiS { double fac) const; protected: - /** \brief - * Vector of DOFVectors to be evaluated at quadrature points. - */ + /// Vector of DOFVectors to be evaluated at quadrature points. std::vector*> vecs; - /** \brief - * Vectors at quadrature points. - */ + /// Vectors at quadrature points. std::vector vecsAtQPs; - /** \brief - * Function for c. - */ + /// Function for c. AbstractFunction > *f; }; @@ -3258,20 +3069,14 @@ namespace AMDiS { VecOfGradientsAtQP_ZOT(const std::vector*>& dv, AbstractFunction*> > *af); - /** \brief - * Implementation of \ref OperatorTerm::initElement(). - */ + /// Implementation of \ref OperatorTerm::initElement(). void initElement(const ElInfo* elInfo, SubAssembler* subAssembler, Quadrature *quad = NULL); - /** \brief - * Implements ZeroOrderTerm::getC(). - */ + /// Implements ZeroOrderTerm::getC(). void getC(const ElInfo *, int nPoints, std::vector &C) const; - /** \brief - * Implements ZeroOrderTerm::eval(). - */ + /// Implements ZeroOrderTerm::eval(). void eval(int nPoints, const double *uhAtQP, const WorldVector *grdUhAtQP, @@ -3280,19 +3085,13 @@ namespace AMDiS { double fac) const; protected: - /** \brief - * Vector of DOFVectors to be evaluated at quadrature points. - */ + /// Vector of DOFVectors to be evaluated at quadrature points. std::vector*> vecs; - /** \brief - * Vectors at quadrature points. - */ + /// Vectors at quadrature points. std::vector*> gradsAtQPs; - /** \brief - * Function for c. - */ + /// Function for c. AbstractFunction*> > *f; }; @@ -3306,20 +3105,14 @@ namespace AMDiS { DOFVectorBase *vec1 = NULL, DOFVectorBase *vec2 = NULL); - /** \brief - * Implementation of \ref OperatorTerm::initElement(). - */ + /// Implementation of \ref OperatorTerm::initElement(). void initElement(const ElInfo* elInfo, SubAssembler* subAssembler, Quadrature *quad = NULL); - /** \brief - * Implements ZeroOrderTerm::getC(). - */ + /// Implements ZeroOrderTerm::getC(). void getC(const ElInfo *, int nPoints, std::vector &C) const; - /** \brief - * Implements ZeroOrderTerm::eval(). - */ + /// Implements ZeroOrderTerm::eval(). void eval(int nPoints, const double *uhAtQP, const WorldVector *grdUhAtQP, @@ -3328,14 +3121,10 @@ namespace AMDiS { double fac) const; protected: - /** \brief - * Vector of DOFVectors to be evaluated at quadrature points. - */ + /// Vector of DOFVectors to be evaluated at quadrature points. std::vector*> vecs; - /** \brief - * Vectors at quadrature points. - */ + /// Vectors at quadrature points. std::vector*> gradsAtQPs; }; @@ -3349,20 +3138,14 @@ namespace AMDiS { const std::vector*>& dv, BinaryAbstractFunction*> > *af); - /** \brief - * Implementation of \ref OperatorTerm::initElement(). - */ + /// Implementation of \ref OperatorTerm::initElement(). void initElement(const ElInfo* elInfo, SubAssembler* subAssembler, Quadrature *quad = NULL); - /** \brief - * Implements ZeroOrderTerm::getC(). - */ + /// Implements ZeroOrderTerm::getC(). void getC(const ElInfo *, int nPoints, std::vector &C) const; - /** \brief - * Implements ZeroOrderTerm::eval(). - */ + /// Implements ZeroOrderTerm::eval(). void eval(int nPoints, const double *uhAtQP, const WorldVector *grdUhAtQP, @@ -3507,20 +3290,14 @@ namespace AMDiS { std::vector, std::vector > > *f); - /** \brief - * Implementation of \ref OperatorTerm::initElement(). - */ + /// Implementation of \ref OperatorTerm::initElement(). void initElement(const ElInfo* elInfo, SubAssembler* subAssembler, Quadrature *quad = NULL); - /** \brief - * Implements ZeroOrderTerm::getC(). - */ + /// Implements ZeroOrderTerm::getC(). void getC(const ElInfo *, int nPoints, std::vector &C) const; - /** \brief - * Implements ZeroOrderTerm::eval(). - */ + /// Implements ZeroOrderTerm::eval(). void eval(int nPoints, const double *uhAtQP, const WorldVector *grdUhAtQP, diff --git a/AMDiS/src/ParallelDomainBase.cc b/AMDiS/src/ParallelDomainBase.cc index 578c9ebee8299e8cc31fb28e5cf1eef14d36cd7d..23aaf48b3a958099dc7156a9817d18e3b335bd70 100644 --- a/AMDiS/src/ParallelDomainBase.cc +++ b/AMDiS/src/ParallelDomainBase.cc @@ -503,7 +503,7 @@ namespace AMDiS { // === Send and recv the nnz row structure to/from other ranks. === - StdMpi stdMpi(mpiComm, true); + StdMpi stdMpi(mpiComm, true); stdMpi.send(sendMatrixEntry); stdMpi.recv(sendDofs); stdMpi.startCommunication(MPI_INT); @@ -611,7 +611,7 @@ namespace AMDiS { // === Synchronize dofs at common dofs, i.e., dofs that correspond to more === // === than one partition. === - synchVectors(vec); + synchVector(vec); // === Print information about solution process. === @@ -636,84 +636,73 @@ namespace AMDiS { KSPDestroy(solver); } - void ParallelDomainBase::synchVectors(SystemVector &vec) + + void ParallelDomainBase::synchVector(DOFVector &vec) { -#if 0 - StdMpi, std::vector > stdMpi(mpiComm); + StdMpi > stdMpi(mpiComm); for (RankToDofContainer::iterator sendIt = sendDofs.begin(); - sendIt != sendDofs.end(); ++sendIt, i++) { + sendIt != sendDofs.end(); ++sendIt) { std::vector dofs; - int nSendDOFs = sendIt->second.size(); - dofs.reserve(nComponents * sendIt->second.size()); + int nSendDofs = sendIt->second.size(); + dofs.reserve(nSendDofs); - for (int j = 0; j < nComponents; j++) { - DOFVector *dofvec = vec.getDOFVector(j); - for (int k = 0; k < nSendDOFs; k++) - dofs.push_back((*dofvec)[*((sendIt->second)[k])]); - } + for (int i = 0; i < nSendDofs; i++) + dofs.push_back(vec[*((sendIt->second)[i])]); stdMpi.send(sendIt->first, dofs); } - stdMpi.startCommunication(); -#endif + for (RankToDofContainer::iterator recvIt = recvDofs.begin(); + recvIt != recvDofs.end(); ++recvIt) + stdMpi.recv(recvIt->first, recvIt->second.size()); -#if 1 + stdMpi.startCommunication(MPI_DOUBLE); - std::vector sendBuffers(sendDofs.size()); - std::vector recvBuffers(recvDofs.size()); + for (RankToDofContainer::iterator recvIt = recvDofs.begin(); + recvIt != recvDofs.end(); ++recvIt) + for (unsigned int i = 0; i < recvIt->second.size(); i++) + vec[*(recvIt->second)[i]] = stdMpi.getRecvData(recvIt->first)[i]; + } - MPI::Request request[sendDofs.size() + recvDofs.size()]; - int requestCounter = 0; - - int i = 0; - for (RankToDofContainer::iterator sendIt = sendDofs.begin(); - sendIt != sendDofs.end(); ++sendIt, i++) { - int nSendDOFs = sendIt->second.size(); - sendBuffers[i] = new double[nSendDOFs * nComponents]; - int counter = 0; - for (int j = 0; j < nComponents; j++) { - DOFVector *dofvec = vec.getDOFVector(j); - for (int k = 0; k < nSendDOFs; k++) - sendBuffers[i][counter++] = (*dofvec)[*((sendIt->second)[k])]; + void ParallelDomainBase::synchVector(SystemVector &vec) + { + StdMpi > stdMpi(mpiComm); + + for (RankToDofContainer::iterator sendIt = sendDofs.begin(); + sendIt != sendDofs.end(); ++sendIt) { + std::vector dofs; + int nSendDofs = sendIt->second.size(); + dofs.reserve(nComponents * nSendDofs); + + for (int i = 0; i < nComponents; i++) { + DOFVector *dofvec = vec.getDOFVector(i); + for (int j = 0; j < nSendDofs; j++) + dofs.push_back((*dofvec)[*((sendIt->second)[j])]); } - request[requestCounter++] = mpiComm.Isend(sendBuffers[i], nSendDOFs * nComponents, - MPI_DOUBLE, sendIt->first, 0); + stdMpi.send(sendIt->first, dofs); } - i = 0; for (RankToDofContainer::iterator recvIt = recvDofs.begin(); - recvIt != recvDofs.end(); ++recvIt, i++) { - int nRecvDOFs = recvIt->second.size() * nComponents; - recvBuffers[i] = new double[nRecvDOFs]; + recvIt != recvDofs.end(); ++recvIt) + stdMpi.recv(recvIt->first, recvIt->second.size() * nComponents); - request[requestCounter++] = - mpiComm.Irecv(recvBuffers[i], nRecvDOFs, MPI_DOUBLE, recvIt->first, 0); - } - - MPI::Request::Waitall(requestCounter, request); + stdMpi.startCommunication(MPI_DOUBLE); - i = 0; for (RankToDofContainer::iterator recvIt = recvDofs.begin(); - recvIt != recvDofs.end(); ++recvIt, i++) { - int nRecvDOFs = recvIt->second.size(); + recvIt != recvDofs.end(); ++recvIt) { + int nRecvDofs = recvIt->second.size(); int counter = 0; - for (int j = 0; j < nComponents; j++) { - DOFVector *dofvec = vec.getDOFVector(j); - for (int k = 0; k < nRecvDOFs; k++) - (*dofvec)[*(recvIt->second)[k]] = recvBuffers[i][counter++]; + for (int i = 0; i < nComponents; i++) { + DOFVector *dofvec = vec.getDOFVector(i); + for (int j = 0; j < nRecvDofs; j++) + (*dofvec)[*(recvIt->second)[j]] = + stdMpi.getRecvData(recvIt->first)[counter++]; } - - delete [] recvBuffers[i]; } - - for (int i = 0; i < static_cast(sendBuffers.size()); i++) - delete [] sendBuffers[i]; -#endif } @@ -785,7 +774,7 @@ namespace AMDiS { } } - StdMpi stdMpi(mpiComm, true); + StdMpi stdMpi(mpiComm, true); stdMpi.send(sendCodes); stdMpi.recv(allBound); stdMpi.startCommunication(MPI_UNSIGNED_LONG); @@ -936,6 +925,117 @@ namespace AMDiS { { FUNCNAME("ParallelDomainBase::createInteriorBoundaryInfo()"); + // === First, create the interior boundaries based on macro element's === + // === neighbour informations. === + + createMacroElementInteriorBoundaryInfo(); + + // === Second, search the whole mesh for interior boundaries that consists of === + // === substructures of the macro elements. === + + createSubstructureInteriorBoundaryInfo(); + + + // === Once we have this information, we must care about the order of the atomic === + // === bounds in the three boundary handling object. Eventually all the bound- === + // === aries have to be in the same order on both ranks that share the bounday. === + + StdMpi > stdMpi(mpiComm); + stdMpi.send(myIntBoundary.boundary); + stdMpi.recv(otherIntBoundary.boundary); + stdMpi.startCommunication(MPI_INT); + + // === The information about all neighbouring boundaries has been received. So === + // === the rank tests if its own atomic boundaries are in the same order. If === + // === not, the atomic boundaries are swaped to the correct order. === + + for (RankToBoundMap::iterator rankIt = otherIntBoundary.boundary.begin(); + rankIt != otherIntBoundary.boundary.end(); ++rankIt) { + + // === We have received from rank "rankIt->first" the ordered list of element === + // === indices. We now have to sort the corresponding list in this rank to === + // === get the same order. === + + for (int j = 0; j < static_cast(rankIt->second.size()); j++) { + // If the expected object is not at place, search for it. + + BoundaryObject &recvedBound = stdMpi.getRecvData()[rankIt->first][j].rankObj; + + if ((rankIt->second)[j].neighObj != recvedBound) { + int k = j + 1; + for (; k < static_cast(rankIt->second.size()); k++) + if ((rankIt->second)[k].neighObj == recvedBound) + break; + + // The element must always be found, because the list is just in another order. + TEST_EXIT_DBG(k < static_cast(rankIt->second.size())) + ("Should never happen!\n"); + + // Swap the current with the found element. + AtomicBoundary tmpBound = (rankIt->second)[k]; + (rankIt->second)[k] = (rankIt->second)[j]; + (rankIt->second)[j] = tmpBound; + } + } + } + + // === Do the same for the periodic boundaries. === + + if (periodicBoundary.boundary.size() > 0) { + stdMpi.clear(); + + InteriorBoundary sendBounds, recvBounds; + for (RankToBoundMap::iterator rankIt = periodicBoundary.boundary.begin(); + rankIt != periodicBoundary.boundary.end(); ++rankIt) { + + TEST_EXIT_DBG(rankIt->first != mpiRank) + ("It is no possible to have an interior boundary within a rank partition!\n"); + + if (rankIt->first < mpiRank) + sendBounds.boundary[rankIt->first] = rankIt->second; + else + recvBounds.boundary[rankIt->first] = rankIt->second; + } + + stdMpi.send(sendBounds.boundary); + stdMpi.recv(recvBounds.boundary); + stdMpi.startCommunication(MPI_INT); + + for (RankToBoundMap::iterator rankIt = periodicBoundary.boundary.begin(); + rankIt != periodicBoundary.boundary.end(); ++rankIt) { + if (rankIt->first <= mpiRank) + continue; + + for (int j = 0; j < static_cast(rankIt->second.size()); j++) { + + BoundaryObject &recvedBound = stdMpi.getRecvData()[rankIt->first][j].rankObj; + + if (periodicBoundary.boundary[rankIt->first][j].neighObj != recvedBound) { + int k = j + 1; + for (; k < static_cast(rankIt->second.size()); k++) + if (periodicBoundary.boundary[rankIt->first][k].neighObj == recvedBound) + break; + + // The element must always be found, because the list is just in + // another order. + TEST_EXIT_DBG(k < static_cast(rankIt->second.size())) + ("Should never happen!\n"); + + // Swap the current with the found element. + AtomicBoundary tmpBound = (rankIt->second)[k]; + (rankIt->second)[k] = (rankIt->second)[j]; + (rankIt->second)[j] = tmpBound; + } + } + } + } // periodicBoundary.boundary.size() > 0 + } + + + void ParallelDomainBase::createMacroElementInteriorBoundaryInfo() + { + FUNCNAME("ParallelDomainBase::createMacroElementInteriorBoundaryInfo()"); + int nNeighbours = mesh->getGeo(NEIGH); int dim = mesh->getDim(); GeoIndex subObj = CENTER; @@ -1032,27 +1132,38 @@ namespace AMDiS { elInfo = stack.traverseNext(elInfo); } + } - // === Seach for all edges which are part of an interior boundary, but are not === - // === a part of some elements at the interior bouundary. === + void ParallelDomainBase::createSubstructureInteriorBoundaryInfo() + { + FUNCNAME("ParallelDomainBase::createSubstructureInteriorBoundaryInfo()"); + + // === Seach for all vertices/edges, which are part of an interior boundary, === + // === but are not a part of the interior boundaries that are created based on === + // === the information of macro elements. === + int dim = mesh->getDim(); const BasisFunction *basFcts = feSpace->getBasisFcts(); - int nBasFcts = basFcts->getNumber(); - std::vector localIndices(nBasFcts, 0); + std::vector localIndices(basFcts->getNumber(), 0); + + // Maps each DOF in the whole domain to the rank that ownes this DOF. + std::map dofOwner; + // Maps each DOF in ranks partition of the domain to the element object that + // contains this DOF. + std::map rankDofs; // Maps each edge in the whole domain to the rank that ownes this edge. std::map edgeOwner; - // Maps each edge in ranks part of the domain to the element object that contains - // this edge. + // Maps each edge in ranks partition of the domain to the element object that + // contains this edge. std::map rankEdges; - std::map dofOwner; - std::map rankDofs; - // === Traverse whole mesh and fill the maps edgeOwner and rankEdges. === + // === Traverse whole mesh and fill the maps defined above. === - elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL); + TraverseStack stack; + ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL); while (elInfo) { Element *el = elInfo->getElement(); basFcts->getLocalIndices(el, feSpace->getAdmin(), localIndices); @@ -1060,66 +1171,76 @@ namespace AMDiS { PartitionElementData *partitionData = dynamic_cast (el->getElementData(PARTITION_ED)); - // Traverse all edges of current element. - for (int i = 0; i < 6; i++) { - DegreeOfFreedom dof0 = localIndices[el->getVertexOfEdge(i, 0)]; - DegreeOfFreedom dof1 = localIndices[el->getVertexOfEdge(i, 1)]; - GlobalEdge edge = std::make_pair(min(dof0, dof1), max(dof0, dof1)); - - // Update the owner of the current edge. - edgeOwner[edge] = max(edgeOwner[edge], partitionVec[el->getIndex()]); - dofOwner[dof0] = max(dofOwner[dof0], partitionVec[el->getIndex()]); - dofOwner[dof1] = max(dofOwner[dof1], partitionVec[el->getIndex()]); - - // If the edge is part of an element that is part of rank's domain, add it - // to the set of all rank's edges. - if (partitionData && partitionData->getPartitionStatus() == IN) { - BoundaryObject b(el, elInfo->getType(), EDGE, i); - rankEdges[edge] = b; - rankDofs[dof0] = b; - rankDofs[dof1] = b; + // In 3d, traverse all edges of current element. + if (dim == 3) { + for (int i = 0; i < 6; i++) { + DegreeOfFreedom dof0 = localIndices[el->getVertexOfEdge(i, 0)]; + DegreeOfFreedom dof1 = localIndices[el->getVertexOfEdge(i, 1)]; + GlobalEdge edge = std::make_pair(min(dof0, dof1), max(dof0, dof1)); + + // Update the owner of the current edge. + edgeOwner[edge] = max(edgeOwner[edge], partitionVec[el->getIndex()]); + + // If the edge is part of an element that is part of rank's domain, add it + // to the set of all rank's edges. + if (partitionData && partitionData->getPartitionStatus() == IN) + rankEdges[edge] = BoundaryObject(el, elInfo->getType(), EDGE, i); } + } + + + // In 2d and 3d, traverse all vertices of current element. + for (int i = 0; i < 4; i++) { + DegreeOfFreedom dof0 = localIndices[i]; + dofOwner[dof0] = max(dofOwner[dof0], partitionVec[el->getIndex()]); + + if (partitionData && partitionData->getPartitionStatus() == IN) + rankDofs[dof0] = BoundaryObject(el, elInfo->getType(), VERTEX, i); } elInfo = stack.traverseNext(elInfo); } - // === Create a set of all edges at rank's interior boundaries. === + + // === Create a set of all edges and vertices at rank's interior boundaries. === // Stores all edges at rank's interior boundaries. std::set rankBoundaryEdges; - + // Stores all vertices at rank's interior boundaries. std::set rankBoundaryDofs; // First, traverse the rank owned elements af the interior boundaries. for (RankToBoundMap::iterator rankIt = myIntBoundary.boundary.begin(); rankIt != myIntBoundary.boundary.end(); ++rankIt) { - for (unsigned int i = 0; i < rankIt->second.size(); i++) { - Element *el = rankIt->second[i].rankObj.el; + for (std::vector::iterator boundIt = rankIt->second.begin(); + boundIt != rankIt->second.end(); ++boundIt) { + Element *el = boundIt->rankObj.el; basFcts->getLocalIndices(el, feSpace->getAdmin(), localIndices); - for (int j = 0; j < 3; j++) { - int edgeNo = el->getEdgeOfFace(rankIt->second[i].rankObj.ithObj, j); - DegreeOfFreedom dof0 = localIndices[el->getVertexOfEdge(edgeNo, 0)]; - DegreeOfFreedom dof1 = localIndices[el->getVertexOfEdge(edgeNo, 1)]; - GlobalEdge edge = std::make_pair(min(dof0, dof1), max(dof0, dof1)); - - // If the edge at rank's interior boundarie has a higher owner rank, than - // we have to remove this edge from the corresponding boundary element. - // Otherwise, it is part of the interior boundary and we add it to the set - // rankBoundaryEdges. - if (edgeOwner[edge] > mpiRank) - rankIt->second[i].rankObj.notIncludedSubStructures.push_back(std::pair(EDGE, edgeNo)); - else - rankBoundaryEdges.insert(edge); + if (dim == 3) { + for (int j = 0; j < 3; j++) { + int edgeNo = el->getEdgeOfFace(boundIt->rankObj.ithObj, j); + DegreeOfFreedom dof0 = localIndices[el->getVertexOfEdge(edgeNo, 0)]; + DegreeOfFreedom dof1 = localIndices[el->getVertexOfEdge(edgeNo, 1)]; + GlobalEdge edge = std::make_pair(min(dof0, dof1), max(dof0, dof1)); + + // If the edge at rank's interior boundarie has a higher owner rank, than + // we have to remove this edge from the corresponding boundary element. + // Otherwise, it is part of the interior boundary and we add it to the set + // rankBoundaryEdges. + if (edgeOwner[edge] > mpiRank) + boundIt->rankObj.notIncludedSubStructures.push_back(std::pair(EDGE, edgeNo)); + else + rankBoundaryEdges.insert(edge); + } } for (int j = 0; j < 3; j++) { - int dofNo = el->getVertexOfPosition(FACE, rankIt->second[i].rankObj.ithObj, j); + int dofNo = el->getVertexOfPosition(FACE, boundIt->rankObj.ithObj, j); DegreeOfFreedom dof = localIndices[dofNo]; if (dofOwner[dof] > mpiRank) - rankIt->second[i].rankObj.notIncludedSubStructures.push_back(std::pair(VERTEX, dofNo)); + boundIt->rankObj.notIncludedSubStructures.push_back(std::pair(VERTEX, dofNo)); else rankBoundaryDofs.insert(dof); } @@ -1129,29 +1250,30 @@ namespace AMDiS { // Now, do the same with all other elements at the interio boundaries. for (RankToBoundMap::iterator rankIt = otherIntBoundary.boundary.begin(); rankIt != otherIntBoundary.boundary.end(); ++rankIt) { - for (unsigned int i = 0; i < rankIt->second.size(); i++) { - Element *el = rankIt->second[i].rankObj.el; + for (std::vector::iterator boundIt = rankIt->second.begin(); + boundIt != rankIt->second.end(); ++boundIt) { + Element *el = boundIt->rankObj.el; basFcts->getLocalIndices(el, feSpace->getAdmin(), localIndices); for (int j = 0; j < 3; j++) { - int edgeNo = el->getEdgeOfFace(rankIt->second[i].rankObj.ithObj, j); + int edgeNo = el->getEdgeOfFace(boundIt->rankObj.ithObj, j); DegreeOfFreedom dof0 = localIndices[el->getVertexOfEdge(edgeNo, 0)]; DegreeOfFreedom dof1 = localIndices[el->getVertexOfEdge(edgeNo, 1)]; GlobalEdge edge = std::make_pair(min(dof0, dof1), max(dof0, dof1)); if (edgeOwner[edge] > rankIt->first) - rankIt->second[i].rankObj.notIncludedSubStructures.push_back(std::pair(EDGE, edgeNo)); + boundIt->rankObj.notIncludedSubStructures.push_back(std::pair(EDGE, edgeNo)); else rankBoundaryEdges.insert(edge); } for (int j = 0; j < 3; j++) { - int dofNo = el->getVertexOfPosition(FACE, rankIt->second[i].rankObj.ithObj, j); + int dofNo = el->getVertexOfPosition(FACE, boundIt->rankObj.ithObj, j); DegreeOfFreedom dof = localIndices[dofNo]; if (dofOwner[dof] > rankIt->first) - rankIt->second[i].rankObj.notIncludedSubStructures.push_back(std::pair(VERTEX, dofNo)); + boundIt->rankObj.notIncludedSubStructures.push_back(std::pair(VERTEX, dofNo)); else rankBoundaryDofs.insert(dof); } @@ -1168,10 +1290,10 @@ namespace AMDiS { // interior boundary with. std::map > recvEdges; // Stores to the map above the corresponding element objects that include the edge. - std::map > recvEdgesObj; + std::map > recvObjs; for (int i = 0; i < mpiSize; i++) { recvEdges[i].clear(); - recvEdgesObj[i].clear(); + recvObjs[i].clear(); } // Traverse all edges in rank's domain. @@ -1184,11 +1306,11 @@ namespace AMDiS { if (edgeOwner[it->first] > mpiRank && rankBoundaryEdges.count(it->first) == 0) { std::vector &ownerEdges = recvEdges[edgeOwner[it->first]]; - + // Check, we have not add this edge before. if (find(ownerEdges.begin(), ownerEdges.end(), it->first) == ownerEdges.end()) { ownerEdges.push_back(it->first); - recvEdgesObj[edgeOwner[it->first]].push_back(it->second); + recvObjs[edgeOwner[it->first]].push_back(it->second); rankBoundaryDofs.insert(it->first.first); rankBoundaryDofs.insert(it->first.second); @@ -1200,15 +1322,15 @@ namespace AMDiS { // === Send all edge interior boundary infos to the owner of the new edge === // === interior boundaries. === - StdMpi, std::vector > stdMpiEdge(mpiComm, true); + StdMpi > stdMpiEdge(mpiComm, true); stdMpiEdge.send(recvEdges); stdMpiEdge.recvFromAll(); stdMpiEdge.startCommunication(MPI_INT); - StdMpi, std::vector > stdMpiEdgeObj(mpiComm, true); - stdMpiEdgeObj.send(recvEdgesObj); - stdMpiEdgeObj.recvFromAll(); - stdMpiEdgeObj.startCommunication(MPI_INT); + StdMpi > stdMpiObjs(mpiComm, true); + stdMpiObjs.send(recvObjs); + stdMpiObjs.recvFromAll(); + stdMpiObjs.startCommunication(MPI_INT); // The owner of an edge interior boundary must send back information about the @@ -1224,12 +1346,12 @@ namespace AMDiS { for (std::map >::iterator it = stdMpiEdge.getRecvData().begin(); it != stdMpiEdge.getRecvData().end(); ++it) { for (unsigned int i = 0; i < it->second.size(); i++) { - TEST_EXIT_DBG(stdMpiEdgeObj.getRecvData(it->first).size() > i) + TEST_EXIT_DBG(stdMpiObjs.getRecvData(it->first).size() > i) ("Should not happen!\n"); AtomicBoundary& b = myIntBoundary.getNewAtomic(it->first); b.rankObj = rankEdges[it->second[i]]; - b.neighObj = stdMpiEdgeObj.getRecvData(it->first)[i]; + b.neighObj = stdMpiObjs.getRecvData(it->first)[i]; sendObjects[it->first].push_back(b.rankObj); @@ -1241,132 +1363,99 @@ namespace AMDiS { // === Send the information about the elements on the owners side of the new === // === edge interior boundaries. === - stdMpiEdgeObj.clear(); - stdMpiEdgeObj.send(sendObjects); - stdMpiEdgeObj.recvFromAll(); - stdMpiEdgeObj.startCommunication(MPI_INT); + stdMpiObjs.clear(); + stdMpiObjs.send(sendObjects); + stdMpiObjs.recvFromAll(); + stdMpiObjs.startCommunication(MPI_INT); // === To the last, all non-owners of a edge interior boundary will now create === // === the necessary boundary objects. === - for (std::map >::iterator it = recvEdgesObj.begin(); - it != recvEdgesObj.end(); ++it) { + for (std::map >::iterator it = recvObjs.begin(); + it != recvObjs.end(); ++it) { if (it->second.size() > 0) { - TEST_EXIT_DBG(it->second.size() == stdMpiEdgeObj.getRecvData(it->first).size()) + TEST_EXIT_DBG(it->second.size() == stdMpiObjs.getRecvData(it->first).size()) ("Should not happen!\n"); for (unsigned int i = 0; i < it->second.size(); i++) { AtomicBoundary& b = otherIntBoundary.getNewAtomic(it->first); b.rankObj = it->second[i]; - b.neighObj = stdMpiEdgeObj.getRecvData(it->first)[i]; + b.neighObj = stdMpiObjs.getRecvData(it->first)[i]; } } } - for (std::map::iterator it = rankDofs.begin(); + // === Create the new interior boundaries consisting only of vertices. As done === + // === above for edge interior boundaries, also the vertex interior boundaries === + // === are created on that ranks, which do not own this part of the boundary. === + + // Maps that stores to each rank number the DOF this rank has an interior + // boundary with. + std::map > recvGlobalDofs; + for (int i = 0; i < mpiSize; i++) { + recvGlobalDofs[i].clear(); + recvObjs[i].clear(); + } + + for (std::map::iterator it = rankDofs.begin(); it != rankDofs.end(); ++it) { if (dofOwner[it->first] > mpiRank && rankBoundaryDofs.count(it->first) == 0) { - MSG("HAB DICH: %d\n", it->first); + std::vector &ownerDofs = recvGlobalDofs[dofOwner[it->first]]; + + if (find(ownerDofs.begin(), ownerDofs.end(), it->first) == ownerDofs.end()) { + ownerDofs.push_back(it->first); + recvObjs[dofOwner[it->first]].push_back(it->second); + } } } - - // === Once we have this information, we must care about the order of the atomic === - // === bounds in the three boundary handling object. Eventually all the bound- === - // === aries have to be in the same order on both ranks that share the bounday. === + StdMpi > stdMpiDofs(mpiComm, true); + stdMpiDofs.send(recvGlobalDofs); + stdMpiDofs.recvFromAll(); + stdMpiDofs.startCommunication(MPI_INT); - StdMpi, std::vector > stdMpi(mpiComm); - stdMpi.send(myIntBoundary.boundary); - stdMpi.recv(otherIntBoundary.boundary); - stdMpi.startCommunication(MPI_INT); - - // === The information about all neighbouring boundaries has been received. So === - // === the rank tests if its own atomic boundaries are in the same order. If === - // === not, the atomic boundaries are swaped to the correct order. === + stdMpiObjs.clear(); + stdMpiObjs.send(recvObjs); + stdMpiObjs.recvFromAll(); + stdMpiObjs.startCommunication(MPI_INT); - for (RankToBoundMap::iterator rankIt = otherIntBoundary.boundary.begin(); - rankIt != otherIntBoundary.boundary.end(); ++rankIt) { - - // === We have received from rank "rankIt->first" the ordered list of element === - // === indices. We now have to sort the corresponding list in this rank to === - // === get the same order. === - - for (int j = 0; j < static_cast(rankIt->second.size()); j++) { - // If the expected object is not at place, search for it. - - BoundaryObject &recvedBound = stdMpi.getRecvData()[rankIt->first][j].rankObj; - - if ((rankIt->second)[j].neighObj != recvedBound) { - int k = j + 1; - for (; k < static_cast(rankIt->second.size()); k++) - if ((rankIt->second)[k].neighObj == recvedBound) - break; - - // The element must always be found, because the list is just in another order. - TEST_EXIT_DBG(k < static_cast(rankIt->second.size())) - ("Should never happen!\n"); + for (int i = 0; i < mpiSize; i++) + sendObjects[i].clear(); - // Swap the current with the found element. - AtomicBoundary tmpBound = (rankIt->second)[k]; - (rankIt->second)[k] = (rankIt->second)[j]; - (rankIt->second)[j] = tmpBound; - } + for (std::map >::iterator it = stdMpiDofs.getRecvData().begin(); + it != stdMpiDofs.getRecvData().end(); ++it) { + for (unsigned int i = 0; i < it->second.size(); i++) { + TEST_EXIT_DBG(stdMpiObjs.getRecvData(it->first).size() > i) + ("Should not happen!\n"); + + AtomicBoundary& b = myIntBoundary.getNewAtomic(it->first); + b.rankObj = rankDofs[it->second[i]]; + b.neighObj = stdMpiObjs.getRecvData(it->first)[i]; + + sendObjects[it->first].push_back(b.rankObj); } } - // === Do the same for the periodic boundaries. === + stdMpiObjs.clear(); + stdMpiObjs.send(sendObjects); + stdMpiObjs.recvFromAll(); + stdMpiObjs.startCommunication(MPI_INT); - if (periodicBoundary.boundary.size() > 0) { - stdMpi.clear(); - InteriorBoundary sendBounds, recvBounds; - for (RankToBoundMap::iterator rankIt = periodicBoundary.boundary.begin(); - rankIt != periodicBoundary.boundary.end(); ++rankIt) { - - TEST_EXIT_DBG(rankIt->first != mpiRank) - ("It is no possible to have an interior boundary within a rank partition!\n"); - - if (rankIt->first < mpiRank) - sendBounds.boundary[rankIt->first] = rankIt->second; - else - recvBounds.boundary[rankIt->first] = rankIt->second; - } - - stdMpi.send(sendBounds.boundary); - stdMpi.recv(recvBounds.boundary); - stdMpi.startCommunication(MPI_INT); + for (std::map >::iterator it = recvObjs.begin(); + it != recvObjs.end(); ++it) { + if (it->second.size() > 0) { + TEST_EXIT_DBG(it->second.size() == stdMpiObjs.getRecvData(it->first).size()) + ("Should not happen!\n"); - for (RankToBoundMap::iterator rankIt = periodicBoundary.boundary.begin(); - rankIt != periodicBoundary.boundary.end(); ++rankIt) { - if (rankIt->first <= mpiRank) - continue; - - for (int j = 0; j < static_cast(rankIt->second.size()); j++) { - - BoundaryObject &recvedBound = stdMpi.getRecvData()[rankIt->first][j].rankObj; - - if (periodicBoundary.boundary[rankIt->first][j].neighObj != recvedBound) { - int k = j + 1; - for (; k < static_cast(rankIt->second.size()); k++) - if (periodicBoundary.boundary[rankIt->first][k].neighObj == recvedBound) - break; - - // The element must always be found, because the list is just in - // another order. - TEST_EXIT_DBG(k < static_cast(rankIt->second.size())) - ("Should never happen!\n"); - - // Swap the current with the found element. - AtomicBoundary tmpBound = (rankIt->second)[k]; - (rankIt->second)[k] = (rankIt->second)[j]; - (rankIt->second)[j] = tmpBound; - } + for (unsigned int i = 0; i < it->second.size(); i++) { + AtomicBoundary& b = otherIntBoundary.getNewAtomic(it->first); + b.rankObj = it->second[i]; + b.neighObj = stdMpiObjs.getRecvData(it->first)[i]; } - } - } // periodicBoundary.boundary.size() > 0 - - // exit(0); + } + } } @@ -1494,13 +1583,15 @@ namespace AMDiS { dofIt != sendIt->second.end(); ++dofIt) sendDofs[sendIt->first].push_back(dofIt->first); - StdMpi > > stdMpi(mpiComm); + typedef std::vector > DofMapVec; + + StdMpi stdMpi(mpiComm); stdMpi.send(sendNewDofs); for (std::map::iterator recvIt = recvNewDofs.begin(); recvIt != recvNewDofs.end(); ++recvIt, i++) stdMpi.recv(recvIt->first, recvIt->second * 2); stdMpi.startCommunication(MPI_INT); - std::map > > &dofMap = stdMpi.getRecvData(); + std::map &dofMap = stdMpi.getRecvData(); // === Change dof indices at boundary from other ranks. === @@ -1726,7 +1817,7 @@ namespace AMDiS { // === Send new DOF indices. === #if 0 - StdMpi, std::vector > stdMpi(mpiComm); + StdMpi > stdMpi(mpiComm); for (RankToDofContainer::iterator sendIt = sendDofs.begin(); sendIt != sendDofs.end(); ++sendIt, i++) { stdMpi.getSendData(sendIt->first).resize(0); @@ -2435,7 +2526,7 @@ namespace AMDiS { // === the same on both corresponding ranks. === typedef std::vector > CoordsVec; - StdMpi stdMpi(mpiComm, true); + StdMpi stdMpi(mpiComm, true); stdMpi.send(sendCoords); stdMpi.recv(recvCoords); stdMpi.startCommunication(MPI_DOUBLE); diff --git a/AMDiS/src/ParallelDomainBase.h b/AMDiS/src/ParallelDomainBase.h index cdc9dceff388f30e072bb1edb74ea5801897dda5..d04c90a1bf78b1d65c68cf6add10246cbe8badc0 100644 --- a/AMDiS/src/ParallelDomainBase.h +++ b/AMDiS/src/ParallelDomainBase.h @@ -193,11 +193,27 @@ namespace AMDiS { protected: /** \brief - * Determine the interior boundaries, i.e. boundaries between ranks, and store + * Determines the interior boundaries, i.e. boundaries between ranks, and stores * all information about them in \ref interiorBoundary. */ void createInteriorBoundaryInfo(); + /** \brief + * Deterimes the interior boundaries between ranks, that are based on the + * neighbourhood information of the macro elements. That means that in 2d the + * function search for all edge based interior boundaries and in 3d for all face + * based interior boundaries. This function cannot find boundaries of substructure + * elements, i.e. vertex boundaries in 2d and vertex and edge boundaries in 3d. + */ + void createMacroElementInteriorBoundaryInfo(); + + /** \brief + * Determines all interior boundaries between rank that consist of element's + * substructures. In 2d these may be only vertices. In 3d there may be also + * interior boundaries consisting of either a whole edge or of a single vertex. + */ + void createSubstructureInteriorBoundaryInfo(); + /// Removes all macro elements from the mesh that are not part of ranks partition. void removeMacroElements(); @@ -342,16 +358,23 @@ namespace AMDiS { void writePartitioningMesh(std::string filename); /** \brief - * This function must be used if the values of a SystemVector must be synchronised + * This function must be used if the values of a DOFVector must be synchronised * over all ranks. That means, that each rank sends the values of the DOFs, which * are owned by the rank and lie on an interior bounday, to all other ranks also - * having this DOF. + * having these DOFs. * * This function must be used, for example, after the lineary system is solved, or - * after the SystemVector is set by some user defined functions, e.g., initial + * after the DOFVector is set by some user defined functions, e.g., initial * solution functions. + */ + void synchVector(DOFVector &vec); + + /** \brief + * Works in the same way as the function above defined for DOFVectors. Due to + * performance, this function does not call \ref synchVector for each DOFVector, + * but instead sends all values of all DOFVectors all at once. */ - void synchVectors(SystemVector &vec); + void synchVector(SystemVector &vec); /// Writes a vector of dof pointers to an output stream. void serialize(std::ostream &out, DofContainer &data); diff --git a/AMDiS/src/ParallelDomainVec.cc b/AMDiS/src/ParallelDomainVec.cc index 54fed89978d92bf171d87977a1ea2106de2ff792..286cbf9e80da94154af1af696c585eaa4ffaf7eb 100644 --- a/AMDiS/src/ParallelDomainVec.cc +++ b/AMDiS/src/ParallelDomainVec.cc @@ -103,7 +103,7 @@ namespace AMDiS { { ParallelDomainBase::solveInitialProblem(adaptInfo); - synchVectors(*(probVec->getSolution())); + synchVector(*(probVec->getSolution())); } void ParallelDomainVec::solve() diff --git a/AMDiS/src/ProblemVec.cc b/AMDiS/src/ProblemVec.cc index 8a65c066e81ae7324869ae5c1b40699db9b6dcfe..5b265a3123d05fb84fffea294b71b372a53dfe64 100644 --- a/AMDiS/src/ProblemVec.cc +++ b/AMDiS/src/ProblemVec.cc @@ -1,5 +1,5 @@ #include -#include "boost/lexical_cast.hpp" +#include #include "ProblemVec.h" #include "RecoveryEstimator.h" #include "Serializer.h" diff --git a/AMDiS/src/StdMpi.h b/AMDiS/src/StdMpi.h index dc29720abfeb497bc4f875c4ae17013d03799eda..41d72e13e22af127bf8be311eb1de31c5f206b2f 100644 --- a/AMDiS/src/StdMpi.h +++ b/AMDiS/src/StdMpi.h @@ -48,6 +48,16 @@ namespace AMDiS { return data.size() * Global::getGeo(WORLD); } + int intSizeOf(std::vector &data) + { + return data.size(); + } + + int intSizeOf(std::vector &data) + { + return data.size(); + } + int intSizeOf(std::vector > &data) { return data.size() * 2; @@ -146,6 +156,34 @@ namespace AMDiS { data[i][j] = buf[pos++]; } + void makeBuf(std::vector &data, int *buf) + { + for (unsigned int i = 0; i < data.size(); i++) + buf[i] = data[i]; + } + + void makeFromBuf(std::vector &data, int *buf, int bufSize) + { + data.resize(bufSize); + + for (int i = 0; i < bufSize; i++) + data[i] = buf[i]; + } + + void makeBuf(std::vector &data, double *buf) + { + for (unsigned int i = 0; i < data.size(); i++) + buf[i] = data[i]; + } + + void makeFromBuf(std::vector &data, double *buf, int bufSize) + { + data.resize(bufSize); + + for (int i = 0; i < bufSize; i++) + data[i] = buf[i]; + } + void makeBuf(std::vector > &data, int *buf) { for (unsigned int i = 0; i < data.size(); i++) { @@ -206,7 +244,7 @@ namespace AMDiS { } } - template + template class StdMpi { public: diff --git a/AMDiS/src/Tetrahedron.cc b/AMDiS/src/Tetrahedron.cc index 424516bd781f4c6f1f8a7d7362f02c62e88e94da..8a34ec25598b27e1184d3b64f6da7c8f60759934 100644 --- a/AMDiS/src/Tetrahedron.cc +++ b/AMDiS/src/Tetrahedron.cc @@ -193,6 +193,9 @@ namespace AMDiS { FUNCNAME("Tetrahedron::getVertexDofs()"); switch (bound.subObj) { + case VERTEX: + dofs.push_back(dof[bound.ithObj]); + break; case EDGE: getVertexDofsAtEdge(feSpace, bound, dofs, parentVertices); break; diff --git a/AMDiS/src/VtkWriter.cc b/AMDiS/src/VtkWriter.cc index b40b91070fd9dae7bcd1fa41ca3b471b781b033c..9bb7c0ac0968ecb8649b2a8b467ee2e1f6548960 100644 --- a/AMDiS/src/VtkWriter.cc +++ b/AMDiS/src/VtkWriter.cc @@ -5,6 +5,10 @@ #include #include #include +#include +#ifdef HAVE_PARALLEL_DOMAIN_AMDIS +#include +#endif #include "VtkWriter.h" #include "DataCollector.h" @@ -125,16 +129,27 @@ namespace AMDiS { void VtkWriter::writeFile(DOFVector *values, std::string filename) { - DataCollector *dc = new DataCollector(values->getFESpace(), values); + FUNCNAME("VtkWriter::writeFile()"); + DataCollector dc(values->getFESpace(), values); std::vector dcList(0); - dcList.push_back(dc); + dcList.push_back(&dc); + VtkWriter writer(&dcList); + +#ifdef HAVE_PARALLEL_DOMAIN_AMDIS + using boost::lexical_cast; + + int sPos = filename.find(".vtu"); + TEST_EXIT(sPos >= 0)("Failed to find file postfix!\n"); + std::string name = filename.substr(0, sPos); - VtkWriter *writer = new VtkWriter(&dcList); - writer->writeFile(filename); + if (MPI::COMM_WORLD.Get_rank() == 0) + writer.writeParallelFile(name + ".pvtu", MPI::COMM_WORLD.Get_size(), name, ".vtu"); + + filename = name + "-p" + lexical_cast(MPI::COMM_WORLD.Get_rank()) + "-.vtu"; +#endif - delete writer; - delete dc; + writer.writeFile(filename); } }