diff --git a/AMDiS/other/include/Makefile_AMDiS.mk b/AMDiS/other/include/Makefile_AMDiS.mk index 322b8fcd04264d262e3223a05ba4612ff3526112..091c1fc70d0020f81d88f57ad4b1726711de67f3 100644 --- a/AMDiS/other/include/Makefile_AMDiS.mk +++ b/AMDiS/other/include/Makefile_AMDiS.mk @@ -20,6 +20,7 @@ AMDIS_LIB = -L$(AMDIS_DIR)/lib -lamdis -lcompositeFEM PNG_LIB = -lpng UMFPACK_LIB = -L$(AMDIS_DIR)/lib/AMD/Lib -L$(AMDIS_DIR)/lib/UMFPACK/Lib +#UMFPACK_LIB = -L/u/witkowski/test/AMD/Lib -L/u/witkowski/test/UMFPACK/Lib ifeq ($(strip $(USE_SERVER)), mars) UMFPACK_LIB += -lmkl -lumfpack -lamd else diff --git a/AMDiS/src/Assembler.cc b/AMDiS/src/Assembler.cc index 3e24e16673168f66d2126658f9052b524de2bca8..0216b10e78a4a3deb68f0662d6c797d0e7460da1 100644 --- a/AMDiS/src/Assembler.cc +++ b/AMDiS/src/Assembler.cc @@ -83,7 +83,7 @@ namespace AMDiS { { FUNCNAME("Assembler::calculateElementMatrix()"); - if (remember && ((factor != 1.0) || (operat->uhOld))) + if (remember && (factor != 1.0 || operat->uhOld)) rememberElMat = true; Element *el = smallElInfo->getElement(); @@ -103,10 +103,11 @@ namespace AMDiS { return; } } - + ElementMatrix& mat = rememberElMat ? elementMatrix : userMat; if (secondOrderAssembler) { + ERROR_EXIT("Da muss i noch ma testen!\n"); secondOrderAssembler->calculateElementMatrix(smallElInfo, mat); ElementMatrix m(nRow, nRow); @@ -134,17 +135,17 @@ namespace AMDiS { smallElInfo->getSubElemCoordsMat(m, rowFESpace->getBasisFcts()->getDegree()); ElementMatrix tmpMat(nRow, nRow); - - if (smallElInfo == colElInfo) - tmpMat = m * mat; - else + + if (smallElInfo == colElInfo) + tmpMat = m * mat; + else tmpMat = mat * trans(m); - mat = tmpMat; + mat = tmpMat; } - if (rememberElMat) - userMat += factor * elementMatrix; + if (rememberElMat && &userMat != &elementMatrix) + userMat += factor * elementMatrix; } @@ -178,7 +179,7 @@ namespace AMDiS { if (operat->uhOld && remember) { matVecAssemble(elInfo, vec); if (rememberElVec) - userVec += factor * elementVector; + userVec += factor * elementVector; return; } @@ -189,9 +190,10 @@ namespace AMDiS { zeroOrderAssembler->calculateElementVector(elInfo, vec); if (rememberElVec) - userVec += factor * elementVector; + userVec += factor * elementVector; } + void Assembler::calculateElementVector(const ElInfo *mainElInfo, const ElInfo *auxElInfo, const ElInfo *smallElInfo, @@ -207,8 +209,8 @@ namespace AMDiS { Element *el = mainElInfo->getElement(); if ((el != lastMatEl && el != lastVecEl) || !operat->isOptimized()) - initElement(auxElInfo); - + initElement(smallElInfo, largeElInfo); + if (el != lastVecEl || !operat->isOptimized()) { if (rememberElVec) set_to_zero(elementVector); @@ -223,32 +225,30 @@ namespace AMDiS { ElementVector& vec = rememberElVec ? elementVector : userVec; if (operat->uhOld && remember) { - - if (smallElInfo->getLevel() == largeElInfo->getLevel()) { + if (smallElInfo->getLevel() == largeElInfo->getLevel()) matVecAssemble(auxElInfo, vec); - } else { - matVecAssemble(mainElInfo, auxElInfo, smallElInfo, largeElInfo, vec); - } + else + matVecAssemble(mainElInfo, auxElInfo, smallElInfo, largeElInfo, vec); if (rememberElVec) - userVec += factor * elementVector; + userVec += factor * elementVector; return; } ERROR_EXIT("Not yet implemented!\n"); -// if (firstOrderAssemblerGrdPsi) { -// firstOrderAssemblerGrdPsi->calculateElementVector(elInfo, vec); -// } -// if (zeroOrderAssembler) { -// zeroOrderAssembler->calculateElementVector(elInfo, vec); -// } -// if (rememberElVec) { -// axpy(factor, *elementVector, *userVec); -// } +#if 0 + if (firstOrderAssemblerGrdPsi) + firstOrderAssemblerGrdPsi->calculateElementVector(elInfo, vec); + if (zeroOrderAssembler) + zeroOrderAssembler->calculateElementVector(elInfo, vec); + if (rememberElVec) + axpy(factor, *elementVector, *userVec); +#endif } + void Assembler::matVecAssemble(const ElInfo *elInfo, ElementVector& vec) { FUNCNAME("Assembler::matVecAssemble()"); @@ -258,20 +258,23 @@ namespace AMDiS { operat->uhOld->getLocalVector(el, uhOldLoc); - if (el != lastMatEl) + if (el != lastMatEl) { + set_to_zero(elementMatrix); calculateElementMatrix(elInfo, elementMatrix); - + } + for (int i = 0; i < nRow; i++) { double val = 0.0; - for (int j = 0; j < nRow; j++) { + for (int j = 0; j < nRow; j++) val += elementMatrix[i][j] * uhOldLoc[j]; - } + vec[i] += val; } delete [] uhOldLoc; } + void Assembler::matVecAssemble(const ElInfo *mainElInfo, const ElInfo *auxElInfo, const ElInfo *smallElInfo, const ElInfo *largeElInfo, ElementVector& vec) @@ -289,35 +292,22 @@ namespace AMDiS { const BasisFunction *basFcts = rowFESpace->getBasisFcts(); int nBasFcts = basFcts->getNumber(); - double *uhOldLoc = new double[nBasFcts]; - double *uhOldLoc2 = new double[nBasFcts]; - - operat->uhOld->getLocalVector(auxEl, uhOldLoc); + std::vector<double> uhOldLoc(nBasFcts); - mtl::dense2D<double> m(nRow, nRow); - smallElInfo->getSubElemCoordsMat(m, rowFESpace->getBasisFcts()->getDegree()); - - for (int i = 0; i < nBasFcts; i++) { - uhOldLoc2[i] = 0.0; - for (int j = 0; j < nBasFcts; j++) - uhOldLoc2[i] += m[j][i] * uhOldLoc[i]; - } + operat->uhOld->getLocalVector(auxEl, &(uhOldLoc[0])); - if (mainEl != lastMatEl) { + set_to_zero(elementMatrix); calculateElementMatrix(mainElInfo, auxElInfo, smallElInfo, largeElInfo, - elementMatrix); + elementMatrix); } - + for (int i = 0; i < nBasFcts; i++) { double val = 0.0; for (int j = 0; j < nBasFcts; j++) - val += elementMatrix[i][j] * uhOldLoc2[j]; + val += elementMatrix[i][j] * uhOldLoc[j]; vec[i] += val; - } - - delete [] uhOldLoc; - delete [] uhOldLoc2; + } } diff --git a/AMDiS/src/ComponentTraverseInfo.cc b/AMDiS/src/ComponentTraverseInfo.cc index 772b39fc2fd7f81f2d0274061d74d111e34d5198..3ee7cb140f843548f834aaf2c5073a45a85de764 100644 --- a/AMDiS/src/ComponentTraverseInfo.cc +++ b/AMDiS/src/ComponentTraverseInfo.cc @@ -22,8 +22,10 @@ namespace AMDiS { status = SingleComponentInfo::EQ_SPACES_NO_AUX; } else { status = SingleComponentInfo::EQ_SPACES_WITH_AUX; - for (int i = 0; i < static_cast<int>(auxFeSpaces.size()); i++) { - if (auxFeSpaces[i]->getMesh() != rowFeSpace->getMesh()) { + + for (std::set<const FiniteElemSpace*>::iterator it = auxFeSpaces.begin(); + it != auxFeSpaces.end(); ++it) { + if ((*it)->getMesh() != rowFeSpace->getMesh()) { status = SingleComponentInfo::EQ_SPACES_WITH_DIF_AUX; break; } @@ -34,9 +36,10 @@ namespace AMDiS { status = SingleComponentInfo::DIF_SPACES_NO_AUX; } else { status = SingleComponentInfo::DIF_SPACES_WITH_AUX; - for (int i = 0; i < static_cast<int>(auxFeSpaces.size()); i++) { - if (auxFeSpaces[i]->getMesh() != rowFeSpace->getMesh() && - auxFeSpaces[i]->getMesh() != colFeSpace->getMesh()) { + for (std::set<const FiniteElemSpace*>::iterator it = auxFeSpaces.begin(); + it != auxFeSpaces.end(); ++it) { + if ((*it)->getMesh() != rowFeSpace->getMesh() && + (*it)->getMesh() != colFeSpace->getMesh()) { status = SingleComponentInfo::DIF_SPACES_WITH_DIF_AUX; break; } @@ -72,12 +75,12 @@ namespace AMDiS { for (int i = 0; i < nComponents; i++) { if (matrixComponents[row][i].getColFeSpace() != rowFeSpace) return matrixComponents[row][i].getColFeSpace(); - if (matrixComponents[row][i].getAuxFeSpace(0) != rowFeSpace) - return matrixComponents[row][i].getAuxFeSpace(0); + if (matrixComponents[row][i].getAuxFeSpace() != rowFeSpace) + return matrixComponents[row][i].getAuxFeSpace(); } - if (vectorComponents[row].getAuxFeSpace(0) != rowFeSpace) - return vectorComponents[row].getAuxFeSpace(0); + if (vectorComponents[row].getAuxFeSpace() != rowFeSpace) + return vectorComponents[row].getAuxFeSpace(); return NULL; } diff --git a/AMDiS/src/ComponentTraverseInfo.h b/AMDiS/src/ComponentTraverseInfo.h index 0a7d8a27bf7c086319c9ae9a0c33a1d93db3e9b7..d029e3587f874f79451881892a6731b3a265c710 100644 --- a/AMDiS/src/ComponentTraverseInfo.h +++ b/AMDiS/src/ComponentTraverseInfo.h @@ -22,7 +22,7 @@ #ifndef AMDIS_COMPONENTTRAVERSEINFO_H #define AMDIS_COMPONENTTRAVERSEINFO_H -#include <vector> +#include <set> #include "Global.h" #include "FiniteElemSpace.h" @@ -34,7 +34,6 @@ namespace AMDiS { SingleComponentInfo() : rowFeSpace(NULL), colFeSpace(NULL), - auxFeSpaces(0), status(0) {} @@ -44,14 +43,14 @@ namespace AMDiS { colFeSpace = col; } - void setAuxFeSpaces(std::vector<const FiniteElemSpace*> feSpaces) + void setAuxFeSpaces(std::set<const FiniteElemSpace*> feSpaces) { auxFeSpaces = feSpaces; } void addAuxFeSpace(const FiniteElemSpace *fe) { - auxFeSpaces.push_back(fe); + auxFeSpaces.insert(fe); } bool hasFeSpace() @@ -76,14 +75,16 @@ namespace AMDiS { return colFeSpace; } - const FiniteElemSpace *getAuxFeSpace(int i) + const FiniteElemSpace *getAuxFeSpace() { - return ((i < static_cast<int>(auxFeSpaces.size())) ? auxFeSpaces[i] : NULL); - } + FUNCNAME("SingleComponentInfo::getAuxFeSpace()"); - void setAuxFeSpace(const FiniteElemSpace* fe, int pos) - { - auxFeSpaces[pos] = fe; + TEST_EXIT_DBG(auxFeSpaces.size() <= 1)("More than one aux FE space!\n"); + + if (auxFeSpaces.size() == 1) + return (*(auxFeSpaces.begin())); + + return NULL; } int getStatus() @@ -96,7 +97,7 @@ namespace AMDiS { FiniteElemSpace *colFeSpace; - std::vector<const FiniteElemSpace*> auxFeSpaces; + std::set<const FiniteElemSpace*> auxFeSpaces; /// Status of the component, see the possible status flags below. int status; @@ -177,12 +178,12 @@ namespace AMDiS { const FiniteElemSpace* getAuxFeSpace(int row, int col) { - return matrixComponents[row][col].getAuxFeSpace(0); + return matrixComponents[row][col].getAuxFeSpace(); } const FiniteElemSpace* getAuxFeSpace(int row) { - return vectorComponents[row].getAuxFeSpace(0); + return vectorComponents[row].getAuxFeSpace(); } /** \brief diff --git a/AMDiS/src/DOFVector.cc b/AMDiS/src/DOFVector.cc index 2b41bfe3a7f5a361db804cdf2a3b0537976ba746..b90a5d152a7ad47a41199569a307772ab100e9b9 100644 --- a/AMDiS/src/DOFVector.cc +++ b/AMDiS/src/DOFVector.cc @@ -808,9 +808,11 @@ namespace AMDiS { return; set_to_zero(this->elementVector); + bool addVector = false; if (op) { op->getElementVector(elInfo, this->elementVector); + addVector = true; } else { std::vector<Operator*>::iterator it; std::vector<double*>::iterator factorIt; @@ -818,12 +820,15 @@ namespace AMDiS { for (it = this->operators.begin(), factorIt = this->operatorFactor.begin(); it != this->operators.end(); ++it, ++factorIt) - if ((*it)->getNeedDualTraverse() == false) + if ((*it)->getNeedDualTraverse() == false) { (*it)->getElementVector(elInfo, this->elementVector, *factorIt ? **factorIt : 1.0); + addVector = true; + } } - addElementVector(factor, this->elementVector, bound, elInfo); + if (addVector) + addElementVector(factor, this->elementVector, bound, elInfo); } template<> @@ -839,27 +844,25 @@ namespace AMDiS { return; set_to_zero(this->elementVector); - - if (op) { - ERROR_EXIT("TODO"); - // op->getElementVector(mainElInfo, this->elementVector); - } else { - std::vector<Operator*>::iterator it; - std::vector<double*>::iterator factorIt; - for (it = this->operators.begin(), factorIt = this->operatorFactor.begin(); - it != this->operators.end(); - ++it, ++factorIt) { - if ((*it)->getNeedDualTraverse()) { - (*it)->getElementVector(mainElInfo, auxElInfo, - smallElInfo, largeElInfo, - this->elementVector, - *factorIt ? **factorIt : 1.0); - } + bool addVector = false; + + TEST_EXIT(!op)("Not yet implemented!\n"); + + std::vector<Operator*>::iterator it; + std::vector<double*>::iterator factorIt; + for (it = this->operators.begin(), factorIt = this->operatorFactor.begin(); + it != this->operators.end(); + ++it, ++factorIt) + if ((*it)->getNeedDualTraverse()) { + (*it)->getElementVector(mainElInfo, auxElInfo, + smallElInfo, largeElInfo, + this->elementVector, + *factorIt ? **factorIt : 1.0); + addVector = true; } - } - - addElementVector(factor, this->elementVector, bound, mainElInfo); + if (addVector) + addElementVector(factor, this->elementVector, bound, mainElInfo); } diff --git a/AMDiS/src/ElInfo2d.cc b/AMDiS/src/ElInfo2d.cc index ae87da7a8ae1b82d9f5dc7a5af90f737c7900a7b..44ebec462cb9a194abd04976af1f6240556469b6 100644 --- a/AMDiS/src/ElInfo2d.cc +++ b/AMDiS/src/ElInfo2d.cc @@ -19,28 +19,16 @@ namespace AMDiS { {0.0, 0.0, 1.0}}; mtl::dense2D<double> ElInfo2d::mat_d1(mat_d1_val); - /* + double ElInfo2d::mat_d1_left_val[3][3] = {{0.0, 1.0, 0.5}, {0.0, 0.0, 0.5}, - {1.0, 0.0, 0.0}}; - */ - - double ElInfo2d::mat_d1_left_val[3][3] = {{0.0, 0.0, 1.0}, - {1.0, 0.0, 0.0}, - {0.5, 0.5, 0.0}}; - + {1.0, 0.0, 0.0}}; mtl::dense2D<double> ElInfo2d::mat_d1_left(mat_d1_left_val); - /* + double ElInfo2d::mat_d1_right_val[3][3] = {{0.0, 0.0, 0.5}, {1.0, 0.0, 0.5}, - {0.0, 1.0, 0.0}}; - */ - - double ElInfo2d::mat_d1_right_val[3][3] = {{0.0, 1.0, 0.0}, - {0.0, 0.0, 1.0}, - {0.5, 0.5, 0.0}}; - + {0.0, 1.0, 0.0}}; mtl::dense2D<double> ElInfo2d::mat_d1_right(mat_d1_right_val); diff --git a/AMDiS/src/Estimator.cc b/AMDiS/src/Estimator.cc index c534c2227697107a1835d77e85b81a5e07b8bce9..ba4ca4bca896ee53816c8a632c6d2453376687ed 100644 --- a/AMDiS/src/Estimator.cc +++ b/AMDiS/src/Estimator.cc @@ -24,6 +24,8 @@ namespace AMDiS { FUNCNAME("Estimator::estimate()"); bool dualTraverse = false; + + /* for (unsigned int i = 0; i < matrix.size(); i++) { TEST_EXIT(traverseInfo.getStatus(row, i) != SingleComponentInfo::DIF_SPACES_WITH_DIF_AUX) ("Not yet implemented!\n"); @@ -51,7 +53,11 @@ namespace AMDiS { auxFeSpace->getBasisFcts()->getDegree()) ("Mh, do you really want to do this? Think about it ...\n"); } - + */ + + mesh = uh[row == -1 ? 0 : row]->getFESpace()->getMesh(); + auxMesh = NULL; + init(ts); if (!dualTraverse) diff --git a/AMDiS/src/Operator.cc b/AMDiS/src/Operator.cc index 9d3bd3d74291062b7e60cbddf13e939ba8434729..07ee6458db018f26a2e6da9bd358ba03caabe68c 100644 --- a/AMDiS/src/Operator.cc +++ b/AMDiS/src/Operator.cc @@ -248,9 +248,11 @@ namespace AMDiS { nCol = colFESpace->getBasisFcts()->getNumber(); } - void Operator::setUhOld(const DOFVectorBase<double> *uhOld_) + + void Operator::setUhOld(const DOFVectorBase<double> *vec) { - uhOld = uhOld_; + uhOld = vec; + auxFeSpaces.insert(vec->getFESpace()); } @@ -344,7 +346,7 @@ namespace AMDiS { { TEST_EXIT(dv)("No vector!\n"); - auxFESpaces.push_back(dv->getFESpace()); + auxFeSpaces.insert(dv->getFESpace()); } MultVecAtQP_ZOT::MultVecAtQP_ZOT(DOFVectorBase<double> *dv1, @@ -357,8 +359,8 @@ namespace AMDiS { TEST_EXIT(dv1)("No first vector!\n"); TEST_EXIT(dv2)("No second vector!\n"); - auxFESpaces.push_back(dv1->getFESpace()); - auxFESpaces.push_back(dv2->getFESpace()); + auxFeSpaces.insert(dv1->getFESpace()); + auxFeSpaces.insert(dv2->getFESpace()); } Vec2AtQP_ZOT::Vec2AtQP_ZOT(DOFVectorBase<double> *dv1, @@ -369,8 +371,8 @@ namespace AMDiS { TEST_EXIT(dv1)("No first vector!\n"); TEST_EXIT(dv2)("No second vector!\n"); - auxFESpaces.push_back(dv1->getFESpace()); - auxFESpaces.push_back(dv2->getFESpace()); + auxFeSpaces.insert(dv1->getFESpace()); + auxFeSpaces.insert(dv2->getFESpace()); } Vec3AtQP_ZOT::Vec3AtQP_ZOT(DOFVectorBase<double> *dv1, @@ -383,9 +385,9 @@ namespace AMDiS { TEST_EXIT(dv2)("No second vector!\n"); TEST_EXIT(dv3)("No thierd vector!\n"); - auxFESpaces.push_back(dv1->getFESpace()); - auxFESpaces.push_back(dv2->getFESpace()); - auxFESpaces.push_back(dv3->getFESpace()); + auxFeSpaces.insert(dv1->getFESpace()); + auxFeSpaces.insert(dv2->getFESpace()); + auxFeSpaces.insert(dv3->getFESpace()); } FctGradientCoords_ZOT::FctGradientCoords_ZOT(DOFVectorBase<double> *dv, @@ -394,7 +396,7 @@ namespace AMDiS { { TEST_EXIT(dv)("No vector!\n"); - auxFESpaces.push_back(dv->getFESpace()); + auxFeSpaces.insert(dv->getFESpace()); } VecGradCoordsAtQP_ZOT::VecGradCoordsAtQP_ZOT(DOFVectorBase<double> *dv, @@ -403,7 +405,7 @@ namespace AMDiS { { TEST_EXIT(dv)("No vector!\n"); - auxFESpaces.push_back(dv->getFESpace()); + auxFeSpaces.insert(dv->getFESpace()); } VecAndCoordsAtQP_ZOT::VecAndCoordsAtQP_ZOT(DOFVectorBase<double> *dv, @@ -412,7 +414,7 @@ namespace AMDiS { { TEST_EXIT(dv)("No vector!\n"); - auxFESpaces.push_back(dv->getFESpace()); + auxFeSpaces.insert(dv->getFESpace()); } Vec2AndGradAtQP_ZOT::Vec2AndGradAtQP_ZOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2, @@ -422,8 +424,8 @@ namespace AMDiS { TEST_EXIT(dv1)("No first vector!\n"); TEST_EXIT(dv2)("No second vector!\n"); - auxFESpaces.push_back(dv1->getFESpace()); - auxFESpaces.push_back(dv2->getFESpace()); + auxFeSpaces.insert(dv1->getFESpace()); + auxFeSpaces.insert(dv2->getFESpace()); } FctGradient_ZOT::FctGradient_ZOT(DOFVectorBase<double> *dv, @@ -432,7 +434,7 @@ namespace AMDiS { { TEST_EXIT(dv)("No vector!\n"); - auxFESpaces.push_back(dv->getFESpace()); + auxFeSpaces.insert(dv->getFESpace()); } VecAndGradAtQP_ZOT::VecAndGradAtQP_ZOT(DOFVectorBase<double> *dv, @@ -441,7 +443,7 @@ namespace AMDiS { { TEST_EXIT(dv)("No vector!\n"); - auxFESpaces.push_back(dv->getFESpace()); + auxFeSpaces.insert(dv->getFESpace()); } VecAndGradVecAtQP_ZOT::VecAndGradVecAtQP_ZOT(DOFVectorBase<double> *dv, @@ -452,8 +454,8 @@ namespace AMDiS { TEST_EXIT(dv)("No vector!\n"); TEST_EXIT(dGrd)("No gradient vector!\n"); - auxFESpaces.push_back(dv->getFESpace()); - auxFESpaces.push_back(dGrd->getFESpace()); + auxFeSpaces.insert(dv->getFESpace()); + auxFeSpaces.insert(dGrd->getFESpace()); } VecAndGradVec2AtQP_ZOT::VecAndGradVec2AtQP_ZOT(DOFVectorBase<double> *dv, @@ -466,9 +468,9 @@ namespace AMDiS { TEST_EXIT(dGrd1)("No first gradient vector!\n"); TEST_EXIT(dGrd2)("No second gradient vector!\n"); - auxFESpaces.push_back(dv->getFESpace()); - auxFESpaces.push_back(dGrd1->getFESpace()); - auxFESpaces.push_back(dGrd2->getFESpace()); + auxFeSpaces.insert(dv->getFESpace()); + auxFeSpaces.insert(dGrd1->getFESpace()); + auxFeSpaces.insert(dGrd2->getFESpace()); } VecOfDOFVecsAtQP_ZOT::VecOfDOFVecsAtQP_ZOT(const std::vector<DOFVectorBase<double>*>& dv, @@ -480,7 +482,7 @@ namespace AMDiS { for (int i = 0; i < static_cast<int>(dv.size()); i++) { TEST_EXIT(dv[i])("One vector is NULL!\n"); - auxFESpaces.push_back(dv[i]->getFESpace()); + auxFeSpaces.insert(dv[i]->getFESpace()); } } @@ -493,7 +495,7 @@ namespace AMDiS { for (int i = 0; i < static_cast<int>(dv.size()); i++) { TEST_EXIT(dv[i])("One vector is NULL!\n"); - auxFESpaces.push_back(dv[i]->getFESpace()); + auxFeSpaces.insert(dv[i]->getFESpace()); } } @@ -509,11 +511,11 @@ namespace AMDiS { vecs[1] = vec1; vecs[2] = vec2; - auxFESpaces.push_back(vec0->getFESpace()); + auxFeSpaces.insert(vec0->getFESpace()); if (vec1) - auxFESpaces.push_back(vec1->getFESpace()); + auxFeSpaces.insert(vec1->getFESpace()); if (vec2) - auxFESpaces.push_back(vec2->getFESpace()); + auxFeSpaces.insert(vec2->getFESpace()); } @@ -526,11 +528,11 @@ namespace AMDiS { TEST_EXIT(v)("No vector!\n"); - auxFESpaces.push_back(v->getFESpace()); + auxFeSpaces.insert(v->getFESpace()); for (int i = 0; i < static_cast<int>(dv.size()); i++) { TEST_EXIT(dv[i])("One gradient vector is NULL!\n"); - auxFESpaces.push_back(dv[i]->getFESpace()); + auxFeSpaces.insert(dv[i]->getFESpace()); } } @@ -543,8 +545,8 @@ namespace AMDiS { TEST_EXIT(dv1)("No first vector!\n"); TEST_EXIT(dv2)("No second vector!\n"); - auxFESpaces.push_back(dv1->getFESpace()); - auxFESpaces.push_back(dv2->getFESpace()); + auxFeSpaces.insert(dv1->getFESpace()); + auxFeSpaces.insert(dv2->getFESpace()); } void Vec2AndGrad2AtQP_ZOT::initElement(const ElInfo* elInfo, @@ -589,9 +591,9 @@ namespace AMDiS { TEST_EXIT(dv2)("No vector!\n"); TEST_EXIT(dGrd)("No gradient vector!\n"); - auxFESpaces.push_back(dv1->getFESpace()); - auxFESpaces.push_back(dv2->getFESpace()); - auxFESpaces.push_back(dGrd->getFESpace()); + auxFeSpaces.insert(dv1->getFESpace()); + auxFeSpaces.insert(dv2->getFESpace()); + auxFeSpaces.insert(dGrd->getFESpace()); } void Vec2AndGradVecAtQP_ZOT::initElement(const ElInfo* elInfo, @@ -633,13 +635,13 @@ namespace AMDiS { for (int i = 0; i < static_cast<int>(vecs.size()); i++) { TEST_EXIT(vecs[i])("One vector is NULL!\n"); - auxFESpaces.push_back(vecs[i]->getFESpace()); + auxFeSpaces.insert(vecs[i]->getFESpace()); } for (int i = 0; i < static_cast<int>(grads.size()); i++) { TEST_EXIT(grads[i])("One gradient vector is NULL!\n"); - auxFESpaces.push_back(grads[i]->getFESpace()); + auxFeSpaces.insert(grads[i]->getFESpace()); } } @@ -658,13 +660,13 @@ namespace AMDiS { for (int i = 0; i < static_cast<int>(vecs.size()); i++) { TEST_EXIT(vecs[i])("One vector is NULL!\n"); - auxFESpaces.push_back(vecs[i]->getFESpace()); + auxFeSpaces.insert(vecs[i]->getFESpace()); } for (int i = 0; i < static_cast<int>(grads.size()); i++) { TEST_EXIT(grads[i])("One gradient vector is NULL!\n"); - auxFESpaces.push_back(grads[i]->getFESpace()); + auxFeSpaces.insert(grads[i]->getFESpace()); } } @@ -675,7 +677,7 @@ namespace AMDiS { { TEST_EXIT(dv)("No vector!\n"); - auxFESpaces.push_back(dv->getFESpace()); + auxFeSpaces.insert(dv->getFESpace()); } VecAtQP_FOT::VecAtQP_FOT(DOFVectorBase<double> *dv, @@ -686,7 +688,7 @@ namespace AMDiS { TEST_EXIT(dv)("No vector!\n"); bOne = bIdx; - auxFESpaces.push_back(dv->getFESpace()); + auxFeSpaces.insert(dv->getFESpace()); } VectorGradient_FOT::VectorGradient_FOT(DOFVectorBase<double> *dv, @@ -695,7 +697,7 @@ namespace AMDiS { { TEST_EXIT(dv)("No vector!\n"); - auxFESpaces.push_back(dv->getFESpace()); + auxFeSpaces.insert(dv->getFESpace()); } VectorFct_FOT::VectorFct_FOT(DOFVectorBase<double> *dv, @@ -704,7 +706,7 @@ namespace AMDiS { { TEST_EXIT(dv)("No vector!\n"); - auxFESpaces.push_back(dv->getFESpace()); + auxFeSpaces.insert(dv->getFESpace()); } VecGrad_FOT::VecGrad_FOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2, @@ -714,8 +716,8 @@ namespace AMDiS { TEST_EXIT(dv1)("No first vector!\n"); TEST_EXIT(dv2)("No second vector!\n"); - auxFESpaces.push_back(dv1->getFESpace()); - auxFESpaces.push_back(dv2->getFESpace()); + auxFeSpaces.insert(dv1->getFESpace()); + auxFeSpaces.insert(dv2->getFESpace()); } Vec2AndGradAtQP_FOT::Vec2AndGradAtQP_FOT(DOFVectorBase<double> *dv1, @@ -727,8 +729,8 @@ namespace AMDiS { f(f_), b(b_) { - auxFESpaces.push_back(dv1->getFESpace()); - auxFESpaces.push_back(dv2->getFESpace()); + auxFeSpaces.insert(dv1->getFESpace()); + auxFeSpaces.insert(dv2->getFESpace()); } void Vec2AndGradAtQP_FOT::initElement(const ElInfo* elInfo, @@ -770,7 +772,7 @@ namespace AMDiS { WorldVector<double> *b_) : FirstOrderTerm(f_->getDegree()), vec(dv), f(f_), b(b_) { - auxFESpaces.push_back(dv->getFESpace()); + auxFeSpaces.insert(dv->getFESpace()); } @@ -812,8 +814,8 @@ namespace AMDiS { WorldVector<double> *b_) : FirstOrderTerm(0), vec1(dv1), vec2(dv2), b(b_) { - auxFESpaces.push_back(dv1->getFESpace()); - auxFESpaces.push_back(dv2->getFESpace()); + auxFeSpaces.insert(dv1->getFESpace()); + auxFeSpaces.insert(dv2->getFESpace()); } Vec2AtQP_FOT::Vec2AtQP_FOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2, @@ -821,8 +823,8 @@ namespace AMDiS { : FirstOrderTerm(0), vec1(dv1), vec2(dv2) { bOne = bIdx; - auxFESpaces.push_back(dv1->getFESpace()); - auxFESpaces.push_back(dv2->getFESpace()); + auxFeSpaces.insert(dv1->getFESpace()); + auxFeSpaces.insert(dv2->getFESpace()); } void Vec2AtQP_FOT::initElement(const ElInfo* elInfo, @@ -872,9 +874,9 @@ namespace AMDiS { f(f_), b(bvec) { - auxFESpaces.push_back(dv1->getFESpace()); - auxFESpaces.push_back(dv2->getFESpace()); - auxFESpaces.push_back(dv3->getFESpace()); + auxFeSpaces.insert(dv1->getFESpace()); + auxFeSpaces.insert(dv2->getFESpace()); + auxFeSpaces.insert(dv3->getFESpace()); } Vec3FctAtQP_FOT::Vec3FctAtQP_FOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2,DOFVectorBase<double> *dv3, @@ -887,9 +889,9 @@ namespace AMDiS { f(f_) { bOne = b; - auxFESpaces.push_back(dv1->getFESpace()); - auxFESpaces.push_back(dv2->getFESpace()); - auxFESpaces.push_back(dv3->getFESpace()); + auxFeSpaces.insert(dv1->getFESpace()); + auxFeSpaces.insert(dv2->getFESpace()); + auxFeSpaces.insert(dv3->getFESpace()); } void Vec3FctAtQP_FOT::initElement(const ElInfo* elInfo, @@ -948,13 +950,13 @@ namespace AMDiS { for (int i = 0; i < static_cast<int>(vecs.size()); i++) { TEST_EXIT(vecs[i])("One vector is NULL!\n"); - auxFESpaces.push_back(vecs[i]->getFESpace()); + auxFeSpaces.insert(vecs[i]->getFESpace()); } for (int i = 0; i < static_cast<int>(grads.size()); i++) { TEST_EXIT(grads[i])("One gradient vector is NULL!\n"); - auxFESpaces.push_back(grads[i]->getFESpace()); + auxFeSpaces.insert(grads[i]->getFESpace()); } } @@ -973,13 +975,13 @@ namespace AMDiS { for (int i = 0; i < static_cast<int>(vecs.size()); i++) { TEST_EXIT(vecs[i])("One vector is NULL!\n"); - auxFESpaces.push_back(vecs[i]->getFESpace()); + auxFeSpaces.insert(vecs[i]->getFESpace()); } for (int i = 0; i < static_cast<int>(grads.size()); i++) { TEST_EXIT(grads[i])("One gradient vector is NULL!\n"); - auxFESpaces.push_back(grads[i]->getFESpace()); + auxFeSpaces.insert(grads[i]->getFESpace()); } } @@ -997,7 +999,7 @@ namespace AMDiS { TEST_EXIT(dv)("No vector!\n"); - auxFESpaces.push_back(dv->getFESpace()); + auxFeSpaces.insert(dv->getFESpace()); } VecAtQP_SOT::VecAtQP_SOT(DOFVectorBase<double> *dv, @@ -1008,7 +1010,7 @@ namespace AMDiS { TEST_EXIT(dv)("No vector!\n"); - auxFESpaces.push_back(dv->getFESpace()); + auxFeSpaces.insert(dv->getFESpace()); } Vec2AtQP_SOT::Vec2AtQP_SOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2, @@ -1020,8 +1022,8 @@ namespace AMDiS { TEST_EXIT(dv1)("No first vector!\n"); TEST_EXIT(dv2)("No second vector!\n"); - auxFESpaces.push_back(dv1->getFESpace()); - auxFESpaces.push_back(dv2->getFESpace()); + auxFeSpaces.insert(dv1->getFESpace()); + auxFeSpaces.insert(dv2->getFESpace()); } MatrixGradient_SOT::MatrixGradient_SOT(DOFVectorBase<double> *dv, @@ -1039,7 +1041,7 @@ namespace AMDiS { TEST_EXIT(dv)("No vector!\n"); - auxFESpaces.push_back(dv->getFESpace()); + auxFeSpaces.insert(dv->getFESpace()); } FctGradient_SOT::FctGradient_SOT(DOFVectorBase<double> *dv, @@ -1048,7 +1050,7 @@ namespace AMDiS { { TEST_EXIT(dv)("No vector!\n"); - auxFESpaces.push_back(dv->getFESpace()); + auxFeSpaces.insert(dv->getFESpace()); } VecAndGradientAtQP_SOT::VecAndGradientAtQP_SOT(DOFVectorBase<double> *dv, @@ -1057,7 +1059,7 @@ namespace AMDiS { { TEST_EXIT(dv)("No vector!\n"); - auxFESpaces.push_back(dv->getFESpace()); + auxFeSpaces.insert(dv->getFESpace()); } VecMatrixGradientAtQP_SOT::VecMatrixGradientAtQP_SOT(DOFVectorBase<double> *dv, @@ -1071,7 +1073,7 @@ namespace AMDiS { TEST_EXIT(dv)("No vector!\n"); - auxFESpaces.push_back(dv->getFESpace()); + auxFeSpaces.insert(dv->getFESpace()); } VecGradCoordsAtQP_SOT::VecGradCoordsAtQP_SOT(DOFVectorBase<double> *dv, @@ -1081,7 +1083,7 @@ namespace AMDiS { { TEST_EXIT(dv)("No vector!\n"); - auxFESpaces.push_back(dv->getFESpace()); + auxFeSpaces.insert(dv->getFESpace()); } VecAndCoordsAtQP_SOT::VecAndCoordsAtQP_SOT(DOFVectorBase<double> *dv, @@ -1092,7 +1094,7 @@ namespace AMDiS { TEST_EXIT(dv)("No vector!\n"); - auxFESpaces.push_back(dv->getFESpace()); + auxFeSpaces.insert(dv->getFESpace()); } MatrixGradientAndCoords_SOT::MatrixGradientAndCoords_SOT(DOFVectorBase<double> *dv, @@ -1107,7 +1109,7 @@ namespace AMDiS { TEST_EXIT(dv)("No vector!\n"); - auxFESpaces.push_back(dv->getFESpace()); + auxFeSpaces.insert(dv->getFESpace()); } @@ -1129,8 +1131,8 @@ namespace AMDiS { TEST_EXIT_DBG(dv1)("No vector!\n"); TEST_EXIT_DBG(dv2)("No vector!\n"); - auxFESpaces.push_back(dv1->getFESpace()); - auxFESpaces.push_back(dv2->getFESpace()); + auxFeSpaces.insert(dv1->getFESpace()); + auxFeSpaces.insert(dv2->getFESpace()); } void MatrixVec2_SOT::initElement(const ElInfo* elInfo, @@ -1200,13 +1202,13 @@ namespace AMDiS { for (int i = 0; i < static_cast<int>(vecs.size()); i++) { TEST_EXIT(vecs[i])("One vector is NULL!\n"); - auxFESpaces.push_back(vecs[i]->getFESpace()); + auxFeSpaces.insert(vecs[i]->getFESpace()); } for (int i = 0; i < static_cast<int>(grads.size()); i++) { TEST_EXIT(grads[i])("One gradient vector is NULL!\n"); - auxFESpaces.push_back(grads[i]->getFESpace()); + auxFeSpaces.insert(grads[i]->getFESpace()); } } @@ -1233,13 +1235,13 @@ namespace AMDiS { for (int i = 0; i < static_cast<int>(vecs.size()); i++) { TEST_EXIT(vecs[i])("One vector is NULL!\n"); - auxFESpaces.push_back(vecs[i]->getFESpace()); + auxFeSpaces.insert(vecs[i]->getFESpace()); } for (int i = 0; i < static_cast<int>(grads.size()); i++) { TEST_EXIT(grads[i])("One gradient vector is NULL!\n"); - auxFESpaces.push_back(grads[i]->getFESpace()); + auxFeSpaces.insert(grads[i]->getFESpace()); } } @@ -1249,7 +1251,7 @@ namespace AMDiS { { TEST_EXIT(dv)("No vector!\n"); - auxFESpaces.push_back(dv->getFESpace()); + auxFeSpaces.insert(dv->getFESpace()); } VecAtQP_IJ_SOT::VecAtQP_IJ_SOT(DOFVectorBase<double> *dv, @@ -1261,7 +1263,7 @@ namespace AMDiS { TEST_EXIT(dv)("No vector!\n"); - auxFESpaces.push_back(dv->getFESpace()); + auxFeSpaces.insert(dv->getFESpace()); } void MatrixFct_SOT::initElement(const ElInfo* elInfo, diff --git a/AMDiS/src/Operator.h b/AMDiS/src/Operator.h index 3ed2aa4cbcbbeeb46316657ed3423286ff4b8ef7..93482fbfbf82966ab1286bd2b215a71127293a4e 100644 --- a/AMDiS/src/Operator.h +++ b/AMDiS/src/Operator.h @@ -55,7 +55,6 @@ namespace AMDiS { : properties(0), degree(deg), dimOfWorld(Global::getGeo(WORLD)), - auxFESpaces(0), bOne(-1) {} @@ -68,22 +67,17 @@ namespace AMDiS { * and coordinates at quadrature points can be calculated here. */ virtual void initElement(const ElInfo*, SubAssembler*, Quadrature *quad = NULL) - { - FUNCNAME("OperatorTerm::initElement()"); - ERROR_EXIT("This function has not been implemented in the operator!\n"); - } + {} virtual void initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo, SubAssembler*, Quadrature *quad = NULL) - { - FUNCNAME("OperatorTerm::initElement()"); - ERROR_EXIT("The function initElement has not been implemented for dual mesh traverse!\n"); - } + {} + - /// Returs \auxFESpaces, the list of all aux fe spaces the operator makes use off. - std::vector<const FiniteElemSpace*>& getAuxFeSpaces() + /// Returs \auxFeSpaces, the list of all aux fe spaces the operator makes use off. + inline std::set<const FiniteElemSpace*>& getAuxFeSpaces() { - return auxFESpaces; + return auxFeSpaces; } /// Specifies whether the matrix of the term is symmetric @@ -254,7 +248,7 @@ namespace AMDiS { int dimOfWorld; /// List off all fe spaces, the operator term makes use off. - std::vector<const FiniteElemSpace*> auxFESpaces; + std::set<const FiniteElemSpace*> auxFeSpaces; /// Pointer to the Operator this OperatorTerm belongs to. Operator* operat; @@ -3224,10 +3218,10 @@ namespace AMDiS { return colFESpace; } - /// Returns \ref auxFESpaces. - inline std::vector<const FiniteElemSpace*> getAuxFeSpaces() + /// Returns \ref auxFeSpaces. + inline std::set<const FiniteElemSpace*>& getAuxFeSpaces() { - return auxFESpaces; + return auxFeSpaces; } /// Sets \ref uhOld. @@ -3463,7 +3457,7 @@ namespace AMDiS { const FiniteElemSpace *colFESpace; /// List of aux fe spaces, e.g., if a term is multiplied with a DOF vector - std::vector<const FiniteElemSpace*> auxFESpaces; + std::set<const FiniteElemSpace*> auxFeSpaces; /// Number of rows in the element matrix int nRow; diff --git a/AMDiS/src/Operator.hh b/AMDiS/src/Operator.hh index c1b36002be752e22acc5bd80f1a557037bfa63f9..96dacf076857d74e8255ccc69c9a8c039935430f 100644 --- a/AMDiS/src/Operator.hh +++ b/AMDiS/src/Operator.hh @@ -5,9 +5,9 @@ namespace AMDiS { { secondOrder[0].push_back(term); term->operat = this; - auxFESpaces.insert(auxFESpaces.end(), - term->getAuxFeSpaces().begin(), - term->getAuxFeSpaces().end()); + for (std::set<const FiniteElemSpace*>::iterator it = term->getAuxFeSpaces().begin(); + it != term->getAuxFeSpaces().end(); ++it) + auxFeSpaces.insert(*it); for (int i = 1; i < omp_get_overall_max_threads(); i++) { T *newTerm = new T(static_cast<const T>(*term)); @@ -25,17 +25,16 @@ namespace AMDiS { firstOrderGrdPhi[0].push_back(term); } term->operat = this; - auxFESpaces.insert(auxFESpaces.end(), - term->getAuxFeSpaces().begin(), - term->getAuxFeSpaces().end()); + for (std::set<const FiniteElemSpace*>::iterator it = term->getAuxFeSpaces().begin(); + it != term->getAuxFeSpaces().end(); ++it) + auxFeSpaces.insert(*it); for (int i = 1; i < omp_get_overall_max_threads(); i++) { T *newTerm = new T(static_cast<const T>(*term)); - if (type == GRD_PSI) { + if (type == GRD_PSI) firstOrderGrdPsi[i].push_back(newTerm); - } else { - firstOrderGrdPhi[i].push_back(newTerm); - } + else + firstOrderGrdPhi[i].push_back(newTerm); } } @@ -44,9 +43,9 @@ namespace AMDiS { { zeroOrder[0].push_back(term); term->operat = this; - auxFESpaces.insert(auxFESpaces.end(), - term->getAuxFeSpaces().begin(), - term->getAuxFeSpaces().end()); + for (std::set<const FiniteElemSpace*>::iterator it = term->getAuxFeSpaces().begin(); + it != term->getAuxFeSpaces().end(); ++it) + auxFeSpaces.insert(*it); for (int i = 1; i < omp_get_overall_max_threads(); i++) { T *newTerm = new T(static_cast<const T>(*term)); diff --git a/AMDiS/src/ProblemVec.cc b/AMDiS/src/ProblemVec.cc index 5668230f6145569e497299e1198a48f0325ab41b..9816d17f4b37da90967d95b76516417c5646e6a7 100644 --- a/AMDiS/src/ProblemVec.cc +++ b/AMDiS/src/ProblemVec.cc @@ -216,7 +216,8 @@ namespace AMDiS { } // === do global refinements === - for (int i = 0; i < static_cast<int>(meshes.size()); i++) + + for (unsigned int i = 0; i < meshes.size(); i++) if (initFlag.isSet(INIT_MESH) && meshes[i]) refinementManager->globalRefine(meshes[i], globalRefinements); } @@ -662,7 +663,8 @@ namespace AMDiS { for (int j = 0; j < nComponents; j++) { - // std::cout << "-------" << i << " " << j << "----------------" << std::endl; + if (writeAsmInfo) + std::cout << "-------" << i << " " << j << "----------------" << std::endl; // Only if this variable is true, the current matrix will be assembled. bool assembleMatrix = true; @@ -710,6 +712,9 @@ namespace AMDiS { if (traverseInfo.getStatus(i) == SingleComponentInfo::EQ_SPACES_NO_AUX || traverseInfo.getStatus(i) == SingleComponentInfo::EQ_SPACES_WITH_AUX) { + if (writeAsmInfo) + std::cout << "ASM on one mesh!" << std::endl; + // The simplest case: either the right hand side has no operaters, no aux // fe spaces, or all aux fe spaces are equal to the row and col fe space. assembleOnOneMesh(componentSpaces[i], @@ -721,27 +726,41 @@ namespace AMDiS { // Row fe space and col fe space are both equal, but right hand side has at // least one another aux fe space. - if (assembleMatrix) - assembleOnOneMesh(componentSpaces[i], - assembleFlag, - assembleMatrix ? matrix : NULL, - ((i == j) && asmVector) ? rhs->getDOFVector(i) : NULL); + + if (writeAsmInfo) + std::cout << "ASM on one mesh!" << std::endl; + + assembleOnOneMesh(componentSpaces[i], + assembleFlag, + assembleMatrix ? matrix : NULL, + ((i == j) && asmVector) ? rhs->getDOFVector(i) : NULL); - if (i == j && asmVector) + if (i == j && asmVector) { + if (writeAsmInfo) + std::cout << "ASM vector on two meshes!" << std::endl; + assembleOnDifMeshes2(componentSpaces[i], traverseInfo.getAuxFeSpace(i), assembleFlag, NULL, rhs->getDOFVector(i)); + } } else { ERROR_EXIT("Possible? If yes, not yet implemented!\n"); } } else if (traverseInfo.getStatus(i, j) == SingleComponentInfo::EQ_SPACES_WITH_DIF_AUX) { + if (writeAsmInfo) + std::cout << "ASM on one mesh!" << std::endl; + assembleOnOneMesh(componentSpaces[i], assembleFlag, assembleMatrix ? matrix : NULL, ((i == j) && asmVector) ? rhs->getDOFVector(i) : NULL); + + if (writeAsmInfo) + std::cout << "ASM on two meshes!" << std::endl; + assembleOnDifMeshes2(componentSpaces[i], traverseInfo.getAuxFeSpace(i, j), assembleFlag, @@ -750,6 +769,10 @@ namespace AMDiS { } else if (traverseInfo.getStatus(i, j) == SingleComponentInfo::DIF_SPACES_NO_AUX || traverseInfo.getStatus(i, j) == SingleComponentInfo::DIF_SPACES_WITH_AUX) { + + if (writeAsmInfo) + std::cout << "ASM on two meshes!" << std::endl; + assembleOnDifMeshes(componentSpaces[i], componentSpaces[j], assembleFlag, assembleMatrix ? matrix : NULL, @@ -878,10 +901,11 @@ namespace AMDiS { (*systemMatrix)[i][j]->addOperator(op, factor, estFactor); traverseInfo.getMatrix(i, j).setAuxFeSpaces(op->getAuxFeSpaces()); - - for (int k = 0; k < static_cast<int>(op->getAuxFeSpaces().size()); k++) { - if ((op->getAuxFeSpaces())[k]->getMesh() != componentSpaces[i]->getMesh() || - (op->getAuxFeSpaces())[k]->getMesh() != componentSpaces[j]->getMesh()) { + + for (std::set<const FiniteElemSpace*>::iterator it = op->getAuxFeSpaces().begin(); + it != op->getAuxFeSpaces().end(); ++it) { + if ((*it)->getMesh() != componentSpaces[i]->getMesh() || + (*it)->getMesh() != componentSpaces[j]->getMesh()) { op->setNeedDualTraverse(true); break; } @@ -911,9 +935,10 @@ namespace AMDiS { rhs->getDOFVector(i)->addOperator(op, factor, estFactor); traverseInfo.getVector(i).setAuxFeSpaces(op->getAuxFeSpaces()); - - for (int j = 0; j < static_cast<int>(op->getAuxFeSpaces().size()); j++) { - if ((op->getAuxFeSpaces())[j]->getMesh() != componentSpaces[i]->getMesh()) { + + for (std::set<const FiniteElemSpace*>::iterator it = op->getAuxFeSpaces().begin(); + it != op->getAuxFeSpaces().end(); ++it) { + if ((*it)->getMesh() != componentSpaces[i]->getMesh()) { op->setNeedDualTraverse(true); break; } @@ -937,6 +962,9 @@ namespace AMDiS { { FUNCNAME("ProblemVec::addDirichletBC()"); + TEST_EXIT(row >= 0 && row < nComponents)("Wrong row number: %d\n", row); + TEST_EXIT(col >= 0 && col < nComponents)("Wrong col number: %d\n", col); + boundaryConditionSet = true; DirichletBC *dirichletApply = diff --git a/AMDiS/src/ProblemVec.h b/AMDiS/src/ProblemVec.h index 6199f9ace6ce3b1e189eceac0a58951e5788fc7a..3017c89d47e71cb42f3b63bfebfcb70db60201d1 100644 --- a/AMDiS/src/ProblemVec.h +++ b/AMDiS/src/ProblemVec.h @@ -71,7 +71,8 @@ namespace AMDiS { info(10), allowFirstRef(false), computeExactError(false), - boundaryConditionSet(false) + boundaryConditionSet(false), + writeAsmInfo(false) { GET_PARAMETER(0, name + "->components", "%d", &nComponents); TEST_EXIT(nComponents > 0)("components not set!\n"); @@ -470,6 +471,12 @@ namespace AMDiS { computeExactError = v; } + /// Sets \ref writeAsmInfo; + void setWriteAsmInfo(bool b) + { + writeAsmInfo = b; + } + /** \} */ /** \brief @@ -619,6 +626,9 @@ namespace AMDiS { std::map<Operator*, std::vector<OperatorPos> > operators; std::map<Operator*, Flag> opFlags; + + /// If true, AMDiS prints information about the assembling procedure to the screen. + bool writeAsmInfo; }; } diff --git a/AMDiS/src/RefinementManager2d.cc b/AMDiS/src/RefinementManager2d.cc index 4db3aaa9d3fa8045af6ad7ef62d151bb0ec230de..55b86eb51a26bf7951ff3592ec8b4448e9c78566 100644 --- a/AMDiS/src/RefinementManager2d.cc +++ b/AMDiS/src/RefinementManager2d.cc @@ -235,7 +235,7 @@ namespace AMDiS { { FUNCNAME("RefinementManager2d::bisectTriangle()"); TEST_EXIT_DBG(mesh)("no mesh!\n"); - + Triangle *child[2]; child[0] = dynamic_cast<Triangle*>(mesh->createNewElement(el)); diff --git a/AMDiS/src/ResidualEstimator.cc b/AMDiS/src/ResidualEstimator.cc index e81f334d7ad5b6ec5f1b8298b230ac66f44e9480..271a6ae980475956ab84be329035aa4fea36b73a 100644 --- a/AMDiS/src/ResidualEstimator.cc +++ b/AMDiS/src/ResidualEstimator.cc @@ -17,7 +17,8 @@ namespace AMDiS { C0(1.0), C1(1.0), C2(0.0), - C3(1.0) + C3(1.0), + jumpResidualOnly(false) { FUNCNAME("ResidualEstimator::ResidualEstimator()"); @@ -31,6 +32,10 @@ namespace AMDiS { C2 = C2 > 1.e-25 ? sqr(C2) : 0.0; C3 = C3 > 1.e-25 ? sqr(C3) : 0.0; + if (C1 != 0.0 && C0 == 0.0 && C3 == 0.0) + jumpResidualOnly = true; + + TEST_EXIT(C2 == 0.0)("C2 is not used! Please remove it or set it to 0.0!\n"); } @@ -124,6 +129,9 @@ namespace AMDiS { for (int system = 0; system < nSystems; system++) { secondOrderTerms[system] = false; + if (matrix[system] == NULL) + continue; + for (std::vector<Operator*>::iterator it = matrix[system]->getOperators().begin(); it != matrix[system]->getOperators().end(); ++it) secondOrderTerms[system] = secondOrderTerms[system] || (*it)->secondOrderTerms(); @@ -215,23 +223,28 @@ namespace AMDiS { for (it = dofMat->getOperatorsBegin(), itfac = dofMat->getOperatorEstFactorBegin(); it != dofMat->getOperatorsEnd(); ++it, ++itfac) - if (*itfac == NULL || **itfac != 0.0) { - if (dualElInfo) { + if (*itfac == NULL || **itfac != 0.0) { + // If the estimator must only compute the jump residual but there are no + // second order terms in the operator, it can be skipped. + if (jumpResidualOnly && (*it)->secondOrderTerms() == false) + continue; + + if (dualElInfo) (*it)->getAssembler(omp_get_thread_num())->initElement(dualElInfo->smallElInfo, dualElInfo->largeElInfo, quad); - } else - (*it)->getAssembler(omp_get_thread_num())->initElement(elInfo, NULL, quad); + else + (*it)->getAssembler(omp_get_thread_num())->initElement(elInfo, NULL, quad); } if (C0 > 0.0) for (it = dofVec->getOperatorsBegin(); it != dofVec->getOperatorsEnd(); ++it) { - if (dualElInfo) { + if (dualElInfo) (*it)->getAssembler(omp_get_thread_num())->initElement(dualElInfo->smallElInfo, dualElInfo->largeElInfo, quad); - } else - (*it)->getAssembler(omp_get_thread_num())->initElement(elInfo, NULL, quad); + else + (*it)->getAssembler(omp_get_thread_num())->initElement(elInfo, NULL, quad); } } @@ -257,7 +270,7 @@ namespace AMDiS { { FUNCNAME("ResidualEstimator::computeElementResidual()"); - TEST_EXIT(dualElInfo)("Not yet implemented!\n"); + TEST_EXIT(!dualElInfo)("Not yet implemented!\n"); std::vector<Operator*>::iterator it; std::vector<double*>::iterator itfac; diff --git a/AMDiS/src/ResidualEstimator.h b/AMDiS/src/ResidualEstimator.h index ba737d6db20010e60734e0d9c0972ef630ff36e7..79c3112633e9e8474a68f20823cd563ee89d73d9 100644 --- a/AMDiS/src/ResidualEstimator.h +++ b/AMDiS/src/ResidualEstimator.h @@ -111,6 +111,13 @@ namespace AMDiS { /// Constant in front of the time double C3; + /** \brief + * Is true, if C1 != 0, and C0 and C3 = 0, hence only the jump residual must be + * calculated. In this case, some optimizations can be done, because the jump + * residual is calculated only on second order terms. + */ + bool jumpResidualOnly; + /// Number of systems, e.g., number of variables in the equation. int nSystems; diff --git a/AMDiS/src/UmfPackSolver.h b/AMDiS/src/UmfPackSolver.h index d41d0bbd687b9d97fcc3d729926be7d8e0844fc1..c1cf67042f526c783f8a7dd4ab1bcc3079a2ce6b 100644 --- a/AMDiS/src/UmfPackSolver.h +++ b/AMDiS/src/UmfPackSolver.h @@ -50,9 +50,7 @@ namespace AMDiS { public: virtual ~Creator() {} - /** \brief - * Returns a new UmfPackSolver object. - */ + /// Returns a new UmfPackSolver object. OEMSolver* create() { return new UmfPackSolver(this->name); @@ -63,10 +61,13 @@ namespace AMDiS { UmfPackSolver(::std::string name) : OEMSolver(name), solver(0), - store_symbolic(0) + store_symbolic(0), + symmetric_strategy(0) { - GET_PARAMETER(0, name + "->store symbolic", "%d", &store_symbolic); - // GET_PARAMETER(0, name + "->multiple rhs", "%d", &multipleRhs); + GET_PARAMETER(0, name + "->store symbolic", "%d", &store_symbolic); + GET_PARAMETER(0, name + "->symmetric strategy", "%d", &symmetric_strategy); + + // GET_PARAMETER(0, name + "->multiple rhs", "%d", &multipleRhs); } /// Destructor @@ -81,24 +82,31 @@ namespace AMDiS { mtl::dense_vector<value_type>& x, const mtl::dense_vector<value_type>& b) { - if (!solver) - solver= new mtl::matrix::umfpack::solver<matrix_type>(A); - else if (store_symbolic) - solver->update_numeric(); + if (!solver) { + if (!symmetric_strategy) + solver = new mtl::matrix::umfpack::solver<matrix_type>(A); + else + solver = new mtl::matrix::umfpack::solver<matrix_type>(A, UMFPACK_STRATEGY_SYMMETRIC); + } else { + if (store_symbolic) + solver->update_numeric(); else - solver->update(); - - int code= (*solver)(x, b); - mtl::dense_vector<value_type> r(b); r-= A * x; residual= two_norm(r); - std::cout << "UmfPackSolver: ||b-Ax|| = " << residual << "\n"; - return code; + solver->update(); + } + + int code= (*solver)(x, b); + mtl::dense_vector<value_type> r(b); r-= A * x; residual= two_norm(r); + std::cout << "UmfPackSolver: ||b-Ax|| = " << residual << "\n"; + return code; } private: mtl::matrix::umfpack::solver<matrix_type> *solver; - int store_symbolic; + int store_symbolic; + + int symmetric_strategy; }; }