diff --git a/AMDiS/src/Assembler.cc b/AMDiS/src/Assembler.cc
index 0216b10e78a4a3deb68f0662d6c797d0e7460da1..a9a1078e99c6db24f7a03aad9b443d5c31c77c13 100644
--- a/AMDiS/src/Assembler.cc
+++ b/AMDiS/src/Assembler.cc
@@ -23,15 +23,17 @@ namespace AMDiS {
rememberElVec(false),
elementMatrix(nRow, nCol),
elementVector(nRow),
+ tmpMat(nRow, nCol),
lastMatEl(NULL),
lastVecEl(NULL),
lastTraverseId(-1)
-
{}
+
Assembler::~Assembler()
{}
+
void Assembler::calculateElementMatrix(const ElInfo *elInfo,
ElementMatrix& userMat,
double factor)
@@ -110,10 +112,9 @@ namespace AMDiS {
ERROR_EXIT("Da muss i noch ma testen!\n");
secondOrderAssembler->calculateElementMatrix(smallElInfo, mat);
- ElementMatrix m(nRow, nRow);
- smallElInfo->getSubElemCoordsMat_so(m, rowFESpace->getBasisFcts()->getDegree());
+ ElementMatrix &m =
+ smallElInfo->getSubElemCoordsMat_so(rowFESpace->getBasisFcts()->getDegree());
- ElementMatrix tmpMat(nRow, nRow);
tmpMat = m * mat;
mat = tmpMat;
}
@@ -131,10 +132,8 @@ namespace AMDiS {
if (zeroOrderAssembler) {
zeroOrderAssembler->calculateElementMatrix(smallElInfo, mat);
- ElementMatrix m(nRow, nRow);
- smallElInfo->getSubElemCoordsMat(m, rowFESpace->getBasisFcts()->getDegree());
-
- ElementMatrix tmpMat(nRow, nRow);
+ ElementMatrix &m =
+ smallElInfo->getSubElemCoordsMat(rowFESpace->getBasisFcts()->getDegree());
if (smallElInfo == colElInfo)
tmpMat = m * mat;
@@ -366,22 +365,24 @@ namespace AMDiS {
}
}
+
void Assembler::finishAssembling()
{
lastVecEl = NULL;
lastMatEl = NULL;
}
+
OptimizedAssembler::OptimizedAssembler(Operator *op,
Quadrature *quad2,
Quadrature *quad1GrdPsi,
Quadrature *quad1GrdPhi,
Quadrature *quad0,
- const FiniteElemSpace *row,
- const FiniteElemSpace *col)
- : Assembler(op, row, col)
+ const FiniteElemSpace *rowFeSpace,
+ const FiniteElemSpace *colFeSpace)
+ : Assembler(op, rowFeSpace, colFeSpace)
{
- bool opt = (row == col);
+ bool opt = (rowFeSpace->getBasisFcts() == colFeSpace->getBasisFcts());
// create sub assemblers
secondOrderAssembler =
@@ -396,14 +397,15 @@ namespace AMDiS {
checkQuadratures();
}
+
StandardAssembler::StandardAssembler(Operator *op,
Quadrature *quad2,
Quadrature *quad1GrdPsi,
Quadrature *quad1GrdPhi,
Quadrature *quad0,
- const FiniteElemSpace *row,
- const FiniteElemSpace *col)
- : Assembler(op, row, col)
+ const FiniteElemSpace *rowFeSpace,
+ const FiniteElemSpace *colFeSpace)
+ : Assembler(op, rowFeSpace, colFeSpace)
{
remember = false;
diff --git a/AMDiS/src/Assembler.h b/AMDiS/src/Assembler.h
index 37a84e27e67f91a10edb6d3c07e9bdae3e5fe12c..1328cc219b45aad292e2278acd6e98050ef806da 100644
--- a/AMDiS/src/Assembler.h
+++ b/AMDiS/src/Assembler.h
@@ -247,6 +247,9 @@ namespace AMDiS {
/// Locally stored element vector
ElementVector elementVector;
+
+ ///
+ ElementMatrix tmpMat;
/** \brief
* Used to check whether \ref initElement() must be called, because
diff --git a/AMDiS/src/ComponentTraverseInfo.h b/AMDiS/src/ComponentTraverseInfo.h
index d029e3587f874f79451881892a6731b3a265c710..b7123c935900a7e6ee8c89fa7ac4b6b49c3dbe98 100644
--- a/AMDiS/src/ComponentTraverseInfo.h
+++ b/AMDiS/src/ComponentTraverseInfo.h
@@ -87,7 +87,7 @@ namespace AMDiS {
return NULL;
}
- int getStatus()
+ inline int getStatus() const
{
return status;
}
@@ -166,16 +166,30 @@ namespace AMDiS {
return vectorComponents[row];
}
- int getStatus(int row, int col)
+ inline int getStatus(int row, int col) const
{
return matrixComponents[row][col].getStatus();
}
- int getStatus(int row)
+ inline int getStatus(int row) const
{
return vectorComponents[row].getStatus();
}
+ /** \brief
+ * Returns true, if for the given matrix component the row and the col FE spaces
+ * are equal. Note that this is also the case, if there is another aux FE space,
+ * that may be different from both, the row and the col FE spaces.
+ */
+ inline bool eqSpaces(int row, int col) const
+ {
+ int status = matrixComponents[row][col].getStatus();
+
+ return (status == SingleComponentInfo::EQ_SPACES_NO_AUX ||
+ status == SingleComponentInfo::EQ_SPACES_WITH_AUX ||
+ status == SingleComponentInfo::EQ_SPACES_WITH_DIF_AUX);
+ }
+
const FiniteElemSpace* getAuxFeSpace(int row, int col)
{
return matrixComponents[row][col].getAuxFeSpace();
@@ -185,6 +199,24 @@ namespace AMDiS {
{
return vectorComponents[row].getAuxFeSpace();
}
+
+ /// Returns true if there is an aux FE that is different from row and col FE spaces.
+ inline bool difAuxSpace(int row, int col) const
+ {
+ int status = matrixComponents[row][col].getStatus();
+
+ return (status == SingleComponentInfo::EQ_SPACES_WITH_DIF_AUX ||
+ status == SingleComponentInfo::DIF_SPACES_WITH_DIF_AUX);
+ }
+
+ /// Returns true if there is an aux FE that is different from row and col FE spaces.
+ inline bool difAuxSpace(int row) const
+ {
+ int status = vectorComponents[row].getStatus();
+
+ return (status == SingleComponentInfo::EQ_SPACES_WITH_DIF_AUX ||
+ status == SingleComponentInfo::DIF_SPACES_WITH_DIF_AUX);
+ }
/** \brief
* Returns the row FE space for a given row number, i.e., the FE space
diff --git a/AMDiS/src/DOFMatrix.cc b/AMDiS/src/DOFMatrix.cc
index 0ee1b5a45a0be2e6d81c0f8e41b82ed221a56f87..a8d6589e64f49fce08de7d50c282847c9af2c2a2 100644
--- a/AMDiS/src/DOFMatrix.cc
+++ b/AMDiS/src/DOFMatrix.cc
@@ -221,6 +221,7 @@ namespace AMDiS {
DegreeOfFreedom col = colIndices[j];
double entry = elMat[i][j];
+ // std::cout << "ADD at " << row << " " << col << std::endl;
ins[row][col] += entry;
}
}
diff --git a/AMDiS/src/DOFVector.cc b/AMDiS/src/DOFVector.cc
index b90a5d152a7ad47a41199569a307772ab100e9b9..04ae730c246561c245317174612301ffa6b6a81f 100644
--- a/AMDiS/src/DOFVector.cc
+++ b/AMDiS/src/DOFVector.cc
@@ -337,8 +337,7 @@ namespace AMDiS {
getLocalVector(largeElInfo->getElement(), localVec);
const BasisFunction *basFcts = feSpace->getBasisFcts();
- mtl::dense2D<double> m(nBasFcts, nBasFcts);
- smallElInfo->getSubElemCoordsMat(m, basFcts->getDegree());
+ mtl::dense2D<double> &m = smallElInfo->getSubElemCoordsMat(basFcts->getDegree());
DimVec<double> &grd1 = *grdTmp[myRank];
int parts = Global::getGeo(PARTS, dim);
@@ -779,8 +778,7 @@ namespace AMDiS {
double *localVec = localVectors[omp_get_thread_num()];
getLocalVector(largeElInfo->getElement(), localVec);
- mtl::dense2D<double> m(nBasFcts, nBasFcts);
- smallElInfo->getSubElemCoordsMat(m, basFcts->getDegree());
+ mtl::dense2D<double> &m = smallElInfo->getSubElemCoordsMat(basFcts->getDegree());
for (int i = 0; i < nPoints; i++) {
result[i] = 0.0;
diff --git a/AMDiS/src/ElInfo.cc b/AMDiS/src/ElInfo.cc
index fe12e9df36a0633aa5cb0b2d9a326cad663b41b7..6ed47f8feba5295f5c95a7036070d8355152f096 100644
--- a/AMDiS/src/ElInfo.cc
+++ b/AMDiS/src/ElInfo.cc
@@ -14,7 +14,7 @@
namespace AMDiS {
- std::map<unsigned long, mtl::dense2D<double> > ElInfo::subElemMatrices;
+ std::vector<std::map<unsigned long, mtl::dense2D<double> > > ElInfo::subElemMatrices(4);
ElInfo::ElInfo(Mesh *aMesh)
: mesh(aMesh),
diff --git a/AMDiS/src/ElInfo.h b/AMDiS/src/ElInfo.h
index a60d70b194657b6f101ffe3e6edf4fecc1baecb1..bd0f395530d0fc945357db966860a02568d4afc3 100644
--- a/AMDiS/src/ElInfo.h
+++ b/AMDiS/src/ElInfo.h
@@ -239,12 +239,15 @@ namespace AMDiS {
return refinementPathLength;
}
+ virtual mtl::dense2D<double>& getSubElemCoordsMat(int degree) const
+ {
+ return subElemMatrices[degree][refinementPath];
+ }
- virtual void getSubElemCoordsMat(mtl::dense2D<double>& mat,
- int basisFctsDegree) const {}
-
- virtual void getSubElemCoordsMat_so(mtl::dense2D<double>& mat,
- int basisFctsDegree) const {}
+ virtual mtl::dense2D<double>& getSubElemCoordsMat_so(int degree) const
+ {
+ return subElemMatrices[degree][refinementPath];
+ }
/** \} */
@@ -544,7 +547,7 @@ namespace AMDiS {
mtl::dense2D<double> subElemCoordsMat_so;
public:
- static std::map<unsigned long, mtl::dense2D<double> > subElemMatrices;
+ static std::vector<std::map<unsigned long, mtl::dense2D<double> > > subElemMatrices;
/** \brief
* child_vertex[el_type][child][i] = father's local vertex index of new
diff --git a/AMDiS/src/ElInfo1d.cc b/AMDiS/src/ElInfo1d.cc
index 4038ff857016df73c087325b39d5c649b299d4b0..0bf6f774d40057ee3b33a5f99c3d272d493334ad 100644
--- a/AMDiS/src/ElInfo1d.cc
+++ b/AMDiS/src/ElInfo1d.cc
@@ -295,45 +295,69 @@ namespace AMDiS {
}
- void ElInfo1d::getSubElemCoordsMat(mtl::dense2D<double>& mat, int degree) const
+ mtl::dense2D<double>& ElInfo1d::getSubElemCoordsMat(int degree) const
{
FUNCNAME("ElInfo1d::getSubElemCoordsMat()");
- switch (degree) {
- case 1:
- mat = mat_d1;
+ using namespace mtl;
+
+ if (subElemMatrices[degree].count(refinementPath) == 0) {
+ switch (degree) {
+ case 1:
+ {
+ dense2D<double> mat(mat_d1);
+ dense2D<double> tmpMat(num_rows(mat), num_rows(mat));
+
+ for (int i = 0; i < refinementPathLength; i++) {
+ if (refinementPath & (1 << i)) {
+ tmpMat = mat_d1_right * mat;
+ mat = tmpMat;
+ } else {
+ tmpMat = mat_d1_left * mat;
+ mat = tmpMat;
+ }
+ }
- for (int i = 0; i < refinementPathLength; i++) {
- if (refinementPath & (1 << i))
- mat *= mat_d1_left;
- else
- mat *= mat_d1_right;
- }
+ subElemMatrices[1][refinementPath] = mat;
+ }
+
+ break;
+ case 2:
+ {
+ dense2D<double> mat(mat_d2);
+ dense2D<double> tmpMat(num_rows(mat), num_rows(mat));
+
+ for (int i = 0; i < refinementPathLength; i++) {
+ if (refinementPath & (1 << i)) {
+ tmpMat = mat_d2_right * mat;
+ mat = tmpMat;
+ } else {
+ tmpMat = mat_d2_left * mat;
+ mat = tmpMat;
+ }
+ }
- break;
-
- case 2:
- mat = mat_d2;
-
- for (int i = 0; i < refinementPathLength; i++) {
- if (refinementPath & (1 << i))
- mat *= mat_d2_left;
- else
- mat *= mat_d2_right;
+ subElemMatrices[2][refinementPath] = mat;
+ }
+ break;
+ default:
+ ERROR_EXIT("Not supported for basis function degree: %d\n", degree);
}
-
- break;
-
- default:
- ERROR_EXIT("Not supported for basis function degree: %d\n", degree);
}
+
+ return subElemMatrices[degree][refinementPath];
}
- void ElInfo1d::getSubElemCoordsMat_so(mtl::dense2D<double>& mat, int degree) const
+ mtl::dense2D<double>& ElInfo1d::getSubElemCoordsMat_so(int degree) const
{
FUNCNAME("ElInfo1d::getSubElemCoordsMat_so()");
+ ERROR_EXIT("Reimplement!\n");
+
+ return subElemMatrices[degree][refinementPath];
+
+#if 0
switch (degree) {
case 1:
mat = test2_val;
@@ -345,6 +369,7 @@ namespace AMDiS {
default:
ERROR_EXIT("Not supported for basis function degree: %d\n", degree);
}
+#endif
}
diff --git a/AMDiS/src/ElInfo1d.h b/AMDiS/src/ElInfo1d.h
index 35c9ff28420d466549e54a289a7392a2a262a712..e2636f12ed99fbc68f49251005b4e491d7ad0934 100644
--- a/AMDiS/src/ElInfo1d.h
+++ b/AMDiS/src/ElInfo1d.h
@@ -62,9 +62,9 @@ namespace AMDiS {
return (i + 1) % 2;
}
- void getSubElemCoordsMat(mtl::dense2D<double>& mat, int degree) const;
+ mtl::dense2D<double>& getSubElemCoordsMat(int degree) const;
- void getSubElemCoordsMat_so(mtl::dense2D<double>& mat, int degree) const;
+ mtl::dense2D<double>& getSubElemCoordsMat_so(int degree) const;
protected:
static double mat_d1_val[2][2];
diff --git a/AMDiS/src/ElInfo2d.cc b/AMDiS/src/ElInfo2d.cc
index 9fc6000d7bee21654f0c276673846f9012dcf4b7..2b0ffbee2209872e7c6aac8a3f8a4e36fbb9b561 100644
--- a/AMDiS/src/ElInfo2d.cc
+++ b/AMDiS/src/ElInfo2d.cc
@@ -688,46 +688,47 @@ namespace AMDiS {
}
- void ElInfo2d::getSubElemCoordsMat(mtl::dense2D<double>& mat, int degree) const
+ mtl::dense2D<double>& ElInfo2d::getSubElemCoordsMat(int degree) const
{
FUNCNAME("ElInfo2d::getSubElemCoordsMat()");
using namespace mtl;
- switch (degree) {
- case 1:
- {
- if (subElemMatrices.count(refinementPath) > 0) {
- mat = subElemMatrices[refinementPath];
- return;
- }
-
- mat = mat_d1;
- dense2D<double> tmpMat(num_rows(mat), num_rows(mat));
-
- for (int i = 0; i < refinementPathLength; i++) {
- if (refinementPath & (1 << i)) {
- tmpMat = mat_d1_right * mat;
- mat = tmpMat;
- } else {
- tmpMat = mat_d1_left * mat;
- mat = tmpMat;
+ if (subElemMatrices[degree].count(refinementPath) == 0) {
+ switch (degree) {
+ case 1:
+ {
+ dense2D<double> mat(mat_d1);
+ dense2D<double> tmpMat(num_rows(mat), num_rows(mat));
+
+ for (int i = 0; i < refinementPathLength; i++) {
+ if (refinementPath & (1 << i)) {
+ tmpMat = mat_d1_right * mat;
+ mat = tmpMat;
+ } else {
+ tmpMat = mat_d1_left * mat;
+ mat = tmpMat;
+ }
}
- }
- subElemMatrices[refinementPath] = mat;
+ subElemMatrices[1][refinementPath] = mat;
+ }
+ break;
+ default:
+ ERROR_EXIT("Not supported for basis function degree: %d\n", degree);
}
- break;
- default:
- ERROR_EXIT("Not supported for basis function degree: %d\n", degree);
}
+
+ return subElemMatrices[degree][refinementPath];
}
- void ElInfo2d::getSubElemCoordsMat_so(mtl::dense2D<double>& mat, int degree) const
+ mtl::dense2D<double>& ElInfo2d::getSubElemCoordsMat_so(int degree) const
{
FUNCNAME("ElInfo2d::getSubElemCoordsMat_so()");
ERROR_EXIT("Not yet implemented!\n");
+
+ return subElemMatrices[degree][refinementPath];
}
}
diff --git a/AMDiS/src/ElInfo2d.h b/AMDiS/src/ElInfo2d.h
index 4d530ab4b5e5e25617424000624212a439bdf955..155221a3479a85631c24c941cc791bfa8540b347 100644
--- a/AMDiS/src/ElInfo2d.h
+++ b/AMDiS/src/ElInfo2d.h
@@ -58,9 +58,9 @@ namespace AMDiS {
/// 2-dimensional realisation of ElInfo's getElementNormal method.
double getElementNormal(WorldVector<double> &normal) const;
- void getSubElemCoordsMat(mtl::dense2D<double>& mat, int degree) const;
+ mtl::dense2D<double>& getSubElemCoordsMat(int degree) const;
- void getSubElemCoordsMat_so(mtl::dense2D<double>& mat, int degree) const;
+ mtl::dense2D<double>& getSubElemCoordsMat_so(int degree) const;
protected:
/// Temp vectors for function \ref calcGrdLambda.
diff --git a/AMDiS/src/ElInfo3d.cc b/AMDiS/src/ElInfo3d.cc
index 9d72a1056a2c7177a317d21452809e38b7a8d3f4..3e8ace91d9330bdd2c5587121d0acc63631384d4 100644
--- a/AMDiS/src/ElInfo3d.cc
+++ b/AMDiS/src/ElInfo3d.cc
@@ -591,19 +591,23 @@ namespace AMDiS {
}
- void ElInfo3d::getSubElemCoordsMat(mtl::dense2D<double>& mat, int degree) const
+ mtl::dense2D<double>& ElInfo3d::getSubElemCoordsMat(int degree) const
{
FUNCNAME("ElInfo3d::getSubElemCoordsMat()");
ERROR_EXIT("Not yet implemented!\n");
+
+ return subElemMatrices[degree][refinementPath];
}
- void ElInfo3d::getSubElemCoordsMat_so(mtl::dense2D<double>& mat, int degree) const
+ mtl::dense2D<double>& ElInfo3d::getSubElemCoordsMat_so(int degree) const
{
FUNCNAME("ElInfo3d::getSubElemCoordsMat_so()");
ERROR_EXIT("Not yet implemented!\n");
+
+ return subElemMatrices[degree][refinementPath];
}
}
diff --git a/AMDiS/src/ElInfo3d.h b/AMDiS/src/ElInfo3d.h
index f7364b39c27983f679d7f6e1a0e4a6896a578fa5..d08866f37b46368deeee4f7ba4c53f11adcf03ce 100644
--- a/AMDiS/src/ElInfo3d.h
+++ b/AMDiS/src/ElInfo3d.h
@@ -80,9 +80,9 @@ namespace AMDiS {
orientation = o;
}
- void getSubElemCoordsMat(mtl::dense2D<double>& mat, int degree) const;
+ mtl::dense2D<double>& getSubElemCoordsMat(int degree) const;
- void getSubElemCoordsMat_so(mtl::dense2D<double>& mat, int degree) const;
+ mtl::dense2D<double>& getSubElemCoordsMat_so(int degree) const;
protected:
/** \brief
diff --git a/AMDiS/src/FirstOrderAssembler.cc b/AMDiS/src/FirstOrderAssembler.cc
index ce304f0af370e0782c328922741b23ce3b9ff1c2..f07620c2b986b560d2b70081f7dea2fbd096520d 100644
--- a/AMDiS/src/FirstOrderAssembler.cc
+++ b/AMDiS/src/FirstOrderAssembler.cc
@@ -30,12 +30,12 @@ namespace AMDiS {
tmpLb.resize(omp_get_overall_max_threads(), Lb);
}
- FirstOrderAssembler*
- FirstOrderAssembler::getSubAssembler(Operator* op,
- Assembler *assembler,
- Quadrature *quad,
- FirstOrderType type,
- bool optimized)
+
+ FirstOrderAssembler* FirstOrderAssembler::getSubAssembler(Operator* op,
+ Assembler *assembler,
+ Quadrature *quad,
+ FirstOrderType type,
+ bool optimized)
{
std::vector<SubAssembler*> *subAssemblers =
optimized ?
@@ -101,13 +101,17 @@ namespace AMDiS {
return newAssembler;
}
+
Stand10::Stand10(Operator *op, Assembler *assembler, Quadrature *quad)
: FirstOrderAssembler(op, assembler, quad, false, GRD_PSI)
{
+ name = "standard first order assembler";
+
psi = owner->getRowFESpace()->getBasisFcts();
phi = owner->getColFESpace()->getBasisFcts();
}
+
void Stand10::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
{
DimVec<double> grdPsi(dim, DEFAULT_VALUE, 0.0);
@@ -137,6 +141,7 @@ namespace AMDiS {
}
}
+
void Stand10::calculateElementVector(const ElInfo *elInfo, ElementVector& vec)
{
DimVec<double> grdPsi(dim, DEFAULT_VALUE, 0.0);
@@ -163,7 +168,9 @@ namespace AMDiS {
Quad10::Quad10(Operator *op, Assembler *assembler, Quadrature *quad)
: FirstOrderAssembler(op, assembler, quad, true, GRD_PSI)
- {}
+ {
+ name = "fast quadrature first order assembler";
+ }
void Quad10::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
@@ -208,6 +215,7 @@ namespace AMDiS {
}
}
+
void Quad10::calculateElementVector(const ElInfo *elInfo, ElementVector& vec)
{
VectorOfFixVecs<DimVec<double> > *grdPsi;
@@ -244,9 +252,11 @@ namespace AMDiS {
}
}
+
Pre10::Pre10(Operator *op, Assembler *assembler, Quadrature *quad)
: FirstOrderAssembler(op, assembler, quad, true, GRD_PSI)
{
+ name = "precalculated first order assembler";
}
@@ -331,9 +341,13 @@ namespace AMDiS {
}
}
+
Quad01::Quad01(Operator *op, Assembler *assembler, Quadrature *quad)
: FirstOrderAssembler(op, assembler, quad, true, GRD_PHI)
- {}
+ {
+ name = "fast quadrature first order assembler";
+ }
+
void Quad01::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
{
@@ -374,9 +388,13 @@ namespace AMDiS {
}
}
+
Pre01::Pre01(Operator *op, Assembler *assembler, Quadrature *quad)
: FirstOrderAssembler(op, assembler, quad, true, GRD_PHI)
- {}
+ {
+ name = "precalculated first order assembler";
+ }
+
void Pre01::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
{
@@ -422,6 +440,7 @@ namespace AMDiS {
}
}
+
void Pre10::calculateElementVector(const ElInfo *elInfo, ElementVector& vec)
{
const int *k;
diff --git a/AMDiS/src/Operator.cc b/AMDiS/src/Operator.cc
index ea95dfbd80e7b7b938a58ccb469c836055cf30a5..354b97ff7155c098906d8391deb29f9095e493be 100644
--- a/AMDiS/src/Operator.cc
+++ b/AMDiS/src/Operator.cc
@@ -1173,7 +1173,7 @@ namespace AMDiS {
}
void MatrixVec2_SOT::weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
- std::vector<WorldVector<double> > &result) const
+ std::vector<WorldVector<double> > &result)
{
int nPoints = grdUhAtQP.size();
for (int iq = 0; iq < nPoints; iq++)
@@ -1953,7 +1953,7 @@ namespace AMDiS {
}
void MatrixFct_SOT::weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
- std::vector<WorldVector<double> > &result) const
+ std::vector<WorldVector<double> > &result)
{
int nPoints = grdUhAtQP.size();
WorldMatrix<double> A;
@@ -1986,7 +1986,7 @@ namespace AMDiS {
}
void Matrix_SOT::weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
- std::vector<WorldVector<double> > &result) const
+ std::vector<WorldVector<double> > &result)
{
int nPoints = grdUhAtQP.size();
for (int iq = 0; iq < nPoints; iq++)
@@ -2021,7 +2021,7 @@ namespace AMDiS {
void MatrixGradient_SOT::weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
- std::vector<WorldVector<double> > &result) const
+ std::vector<WorldVector<double> > &result)
{
FUNCNAME("MatrixGradient_SOT::weakEval()");
@@ -2029,10 +2029,9 @@ namespace AMDiS {
TEST_EXIT_DBG(gradAtQPs)("Operator was not initialized correctly!\n");
int nPoints = grdUhAtQP.size();
- WorldMatrix<double> A;
for (int iq = 0; iq < nPoints; iq++) {
- A = (*f)(gradAtQPs[iq]);
- result[iq] += A * grdUhAtQP[iq];
+ tmpMat = (*f)(gradAtQPs[iq]);
+ result[iq] += tmpMat * grdUhAtQP[iq];
}
}
@@ -2058,7 +2057,7 @@ namespace AMDiS {
}
void VecAtQP_SOT::weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
- std::vector<WorldVector<double> > &result) const
+ std::vector<WorldVector<double> > &result)
{
int nPoints = grdUhAtQP.size();
for (int iq = 0; iq < nPoints; iq++) {
@@ -2088,7 +2087,7 @@ namespace AMDiS {
}
void Vec2AtQP_SOT::weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
- std::vector<WorldVector<double> > &result) const
+ std::vector<WorldVector<double> > &result)
{
int nPoints = grdUhAtQP.size();
for (int iq = 0; iq < nPoints; iq++) {
@@ -2118,7 +2117,7 @@ namespace AMDiS {
}
void CoordsAtQP_SOT::weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
- std::vector<WorldVector<double> > &result) const
+ std::vector<WorldVector<double> > &result)
{
int nPoints = grdUhAtQP.size();
for (int iq = 0; iq < nPoints; iq++) {
@@ -2148,7 +2147,7 @@ namespace AMDiS {
}
void FctGradient_SOT::weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
- std::vector<WorldVector<double> > &result) const
+ std::vector<WorldVector<double> > &result)
{
int nPoints = grdUhAtQP.size();
for (int iq = 0; iq < nPoints; iq++) {
@@ -2178,7 +2177,7 @@ namespace AMDiS {
}
void VecAndGradientAtQP_SOT::weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
- std::vector<WorldVector<double> > &result) const
+ std::vector<WorldVector<double> > &result)
{
int nPoints = grdUhAtQP.size();
for (int iq = 0; iq < nPoints; iq++) {
@@ -2214,7 +2213,7 @@ namespace AMDiS {
}
void VecMatrixGradientAtQP_SOT::weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
- std::vector<WorldVector<double> > &result) const
+ std::vector<WorldVector<double> > &result)
{
int nPoints = grdUhAtQP.size();
WorldMatrix<double> A;
@@ -2245,7 +2244,7 @@ namespace AMDiS {
}
void VecAndCoordsAtQP_SOT::weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
- std::vector<WorldVector<double> > &result) const
+ std::vector<WorldVector<double> > &result)
{
int nPoints = grdUhAtQP.size();
for (int iq = 0; iq < nPoints; iq++) {
@@ -2281,7 +2280,7 @@ namespace AMDiS {
}
void MatrixGradientAndCoords_SOT::weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
- std::vector<WorldVector<double> > &result) const
+ std::vector<WorldVector<double> > &result)
{
int nPoints = grdUhAtQP.size();
WorldMatrix<double> A;
@@ -2312,7 +2311,7 @@ namespace AMDiS {
}
void VecGradCoordsAtQP_SOT::weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
- std::vector<WorldVector<double> > &result) const
+ std::vector<WorldVector<double> > &result)
{
int nPoints = grdUhAtQP.size();
for (int iq = 0; iq < nPoints; iq++) {
@@ -2482,7 +2481,7 @@ namespace AMDiS {
}
void CoordsAtQP_IJ_SOT::weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
- std::vector<WorldVector<double> > &result) const
+ std::vector<WorldVector<double> > &result)
{
int nPoints = grdUhAtQP.size();
for (int iq = 0; iq < nPoints; iq++) {
@@ -2524,7 +2523,7 @@ namespace AMDiS {
}
void VecAtQP_IJ_SOT::weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
- std::vector<WorldVector<double> > &result) const
+ std::vector<WorldVector<double> > &result)
{
int nPoints = grdUhAtQP.size();
for (int iq = 0; iq < nPoints; iq++) {
@@ -2674,7 +2673,7 @@ namespace AMDiS {
}
void VecAndGradAtQP_SOT::weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
- std::vector<WorldVector<double> > &result) const
+ std::vector<WorldVector<double> > &result)
{
int nPoints = grdUhAtQP.size();
for (int iq = 0; iq < nPoints; iq++) {
@@ -2775,7 +2774,7 @@ namespace AMDiS {
}
void General_SOT::weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
- std::vector<WorldVector<double> > &result) const
+ std::vector<WorldVector<double> > &result)
{
int nPoints = grdUhAtQP.size();
int nVecs = static_cast<int>(vecs_.size());
@@ -2992,7 +2991,7 @@ namespace AMDiS {
}
void GeneralParametric_SOT::weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
- std::vector<WorldVector<double> > &result) const
+ std::vector<WorldVector<double> > &result)
{
int nPoints = grdUhAtQP.size();
int nVecs = static_cast<int>(vecs_.size());
@@ -3242,7 +3241,7 @@ namespace AMDiS {
}
void VecGrad_SOT::weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
- std::vector<WorldVector<double> > &result) const
+ std::vector<WorldVector<double> > &result)
{
int nPoints = grdUhAtQP.size();
for (int iq = 0; iq < nPoints; iq++) {
diff --git a/AMDiS/src/Operator.h b/AMDiS/src/Operator.h
index 2989692821b1a7c97128dbd5b16f0b2e22abceee..54fb1667896338da9d1c333be4c24831018b17aa 100644
--- a/AMDiS/src/Operator.h
+++ b/AMDiS/src/Operator.h
@@ -300,7 +300,7 @@ namespace AMDiS {
/// Evaluation of \f$ A \nabla u(\vec{x}) \f$ at all quadrature points.
virtual void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
- std::vector<WorldVector<double> > &result) const = 0;
+ std::vector<WorldVector<double> > &result) = 0;
};
@@ -351,7 +351,7 @@ namespace AMDiS {
/// Implenetation of SecondOrderTerm::weakEval().
void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
- std::vector<WorldVector<double> > &result) const
+ std::vector<WorldVector<double> > &result)
{
int nPoints = grdUhAtQP.size();
for (int iq = 0; iq < nPoints; iq++)
@@ -418,7 +418,7 @@ namespace AMDiS {
/// Implenetation of SecondOrderTerm::weakEval().
void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
- std::vector<WorldVector<double> > &result) const
+ std::vector<WorldVector<double> > &result)
{
int nPoints = grdUhAtQP.size();
for (int iq = 0; iq < nPoints; iq++)
@@ -464,7 +464,7 @@ namespace AMDiS {
/// Implenetation of SecondOrderTerm::weakEval().
void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
- std::vector<WorldVector<double> > &result) const;
+ std::vector<WorldVector<double> > &result);
protected:
/// DOFVector to be evaluated at quadrature points.
@@ -518,7 +518,7 @@ namespace AMDiS {
/// Implenetation of SecondOrderTerm::weakEval().
void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
- std::vector<WorldVector<double> > &result) const;
+ std::vector<WorldVector<double> > &result);
protected:
/// Matrix stroring A.
@@ -579,7 +579,7 @@ namespace AMDiS {
/// Implenetation of SecondOrderTerm::weakEval().
void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
- std::vector<WorldVector<double> > &result) const
+ std::vector<WorldVector<double> > &result)
{
int nPoints = grdUhAtQP.size();
for (int iq = 0; iq < nPoints; iq++)
@@ -630,7 +630,7 @@ namespace AMDiS {
/// Implenetation of SecondOrderTerm::weakEval().
void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
- std::vector<WorldVector<double> > &result) const;
+ std::vector<WorldVector<double> > &result);
protected:
/// DOFVector to be evaluated at quadrature points.
@@ -674,7 +674,7 @@ namespace AMDiS {
/// Implenetation of SecondOrderTerm::weakEval().
void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
- std::vector<WorldVector<double> > &result) const;
+ std::vector<WorldVector<double> > &result);
protected:
/// DOFVector to be evaluated at quadrature points.
@@ -722,7 +722,7 @@ namespace AMDiS {
/// Implenetation of SecondOrderTerm::weakEval().
void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
- std::vector<WorldVector<double> > &result) const;
+ std::vector<WorldVector<double> > &result);
protected:
@@ -773,7 +773,7 @@ namespace AMDiS {
/// Implenetation of SecondOrderTerm::weakEval().
void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
- std::vector<WorldVector<double> > &result) const;
+ std::vector<WorldVector<double> > &result);
protected:
DOFVectorBase<double>* vec;
@@ -791,6 +791,10 @@ namespace AMDiS {
/// True, if \ref f provides always a symmetric WorldMatrix<double>.
bool symmetric;
+
+ private:
+ /// Temporary matrix used in \ref MatrixGradient_SOT::weakEval.
+ WorldMatrix<double> tmpMat;
};
/**
@@ -825,7 +829,7 @@ namespace AMDiS {
/// Implenetation of SecondOrderTerm::weakEval().
void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
- std::vector<WorldVector<double> > &result) const;
+ std::vector<WorldVector<double> > &result);
protected:
DOFVectorBase<double>* vec;
@@ -873,7 +877,7 @@ namespace AMDiS {
/// Implenetation of SecondOrderTerm::weakEval().
void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
- std::vector<WorldVector<double> > &result) const;
+ std::vector<WorldVector<double> > &result);
protected:
DOFVectorBase<double>* vec;
@@ -923,7 +927,7 @@ namespace AMDiS {
/// Implenetation of SecondOrderTerm::weakEval().
void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
- std::vector<WorldVector<double> > &result) const;
+ std::vector<WorldVector<double> > &result);
protected:
DOFVectorBase<double>* vec;
@@ -977,7 +981,7 @@ namespace AMDiS {
/// Implenetation of SecondOrderTerm::weakEval().
void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
- std::vector<WorldVector<double> > &result) const;
+ std::vector<WorldVector<double> > &result);
protected:
DOFVectorBase<double>* vec;
@@ -1026,7 +1030,7 @@ namespace AMDiS {
double factor);
void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
- std::vector<WorldVector<double> > &result) const;
+ std::vector<WorldVector<double> > &result);
protected:
/// DOFVector to be evaluated at quadrature points.
@@ -1078,7 +1082,7 @@ namespace AMDiS {
/// Implenetation of SecondOrderTerm::weakEval().
void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
- std::vector<WorldVector<double> > &result) const;
+ std::vector<WorldVector<double> > &result);
protected:
DOFVectorBase<double>* vec;
@@ -1127,7 +1131,7 @@ namespace AMDiS {
/// Implenetation of SecondOrderTerm::weakEval().
void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
- std::vector<WorldVector<double> > &result) const;
+ std::vector<WorldVector<double> > &result);
protected:
/// DOFVector to be evaluated at quadrature points.
@@ -1179,7 +1183,7 @@ namespace AMDiS {
/// Implenetation of SecondOrderTerm::weakEval().
void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
- std::vector<WorldVector<double> > &result) const;
+ std::vector<WorldVector<double> > &result);
protected:
std::vector<DOFVectorBase<double>*> vecs_;
@@ -1236,7 +1240,7 @@ namespace AMDiS {
/// Implenetation of SecondOrderTerm::weakEval().
void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
- std::vector<WorldVector<double> > &result) const;
+ std::vector<WorldVector<double> > &result);
protected:
std::vector<DOFVectorBase<double>*> vecs_;
@@ -2529,7 +2533,7 @@ namespace AMDiS {
/// Implements SecondOrderTerm::weakEval().
void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
- std::vector<WorldVector<double> > &result) const;
+ std::vector<WorldVector<double> > &result);
protected:
/// DOFVector to be evaluated at quadrature points.
@@ -2579,7 +2583,7 @@ namespace AMDiS {
/// Implements SecondOrderTerm::weakEval().
void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
- std::vector<WorldVector<double> > &result) const;
+ std::vector<WorldVector<double> > &result);
protected:
/// DOFVectors to be evaluated at quadrature points.
@@ -2671,7 +2675,7 @@ namespace AMDiS {
/// Implements SecondOrderTerm::weakEval().
void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
- std::vector<WorldVector<double> > &result) const;
+ std::vector<WorldVector<double> > &result);
private:
/// Stores coordinates at quadrature points. Set in \ref initElement().
@@ -2717,7 +2721,7 @@ namespace AMDiS {
/// Implements SecondOrderTerm::weakEval().
void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
- std::vector<WorldVector<double> > &result) const;
+ std::vector<WorldVector<double> > &result);
protected:
/// DOFVector to be evaluated at quadrature points.
@@ -3240,7 +3244,7 @@ namespace AMDiS {
}
/// Returns \ref assembler
- Assembler *getAssembler(int rank);
+ Assembler *getAssembler(int rank = 0);
/// Sets \ref assembler
void setAssembler(int rank, Assembler *ass);
diff --git a/AMDiS/src/ProblemVec.cc b/AMDiS/src/ProblemVec.cc
index db21cd06e7254450d2b70ee39eb9cdf8fc1dab3f..822cd79bad16c18a0f08387b77541b2a6ea2be01 100644
--- a/AMDiS/src/ProblemVec.cc
+++ b/AMDiS/src/ProblemVec.cc
@@ -625,8 +625,12 @@ namespace AMDiS {
{
FUNCNAME("ProblemVec::buildAfterCoarsen()");
- // printOpenmpTraverseInfo(this, true);
+ if (dualMeshTraverseRequired()) {
+ dualAssemble(adaptInfo, flag);
+ return;
+ }
+ // printOpenmpTraverseInfo(this, true);
// std::cout << "ElInfo = " << ElInfo::subElemMatrices.size() << std::endl;
clock_t first = clock();
@@ -674,6 +678,21 @@ namespace AMDiS {
// will be set to false).
DOFMatrix *matrix = (*systemMatrix)[i][j];
+ if (writeAsmInfo && matrix) {
+ for (std::vector<Operator*>::iterator it = matrix->getOperatorsBegin();
+ it != matrix->getOperatorsEnd(); ++it) {
+ Assembler *assembler = (*it)->getAssembler();
+ if (assembler->getZeroOrderAssembler())
+ std::cout << "ZOA: " << assembler->getZeroOrderAssembler()->getName() << std::endl;
+ if (assembler->getFirstOrderAssembler(GRD_PSI))
+ std::cout << "FOA GRD_PSI: " << assembler->getFirstOrderAssembler(GRD_PSI)->getName() << std::endl;
+ if (assembler->getFirstOrderAssembler(GRD_PHI))
+ std::cout << "FOA GRD_PHI: " << assembler->getFirstOrderAssembler(GRD_PHI)->getName() << std::endl;
+ if (assembler->getSecondOrderAssembler())
+ std::cout << "SOA: " << assembler->getSecondOrderAssembler()->getName() << std::endl;
+ }
+ }
+
if (matrix)
matrix->calculateNnz();
@@ -706,93 +725,290 @@ namespace AMDiS {
if (assembleMatrix && matrix->getBoundaryManager())
matrix->getBoundaryManager()->initMatrix(matrix);
- if (traverseInfo.getStatus(i, j) == SingleComponentInfo::EQ_SPACES_NO_AUX ||
- traverseInfo.getStatus(i, j) == SingleComponentInfo::EQ_SPACES_WITH_AUX) {
+ // 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],
+ assembleFlag,
+ assembleMatrix ? matrix : NULL,
+ ((i == j) && asmVector) ? rhs->getDOFVector(i) : NULL);
+
+ assembledMatrix[i][j] = true;
- // Row fe space and col fe space are both equal
+ if (assembleMatrix)
+ matrix->finishInsertion();
- if (traverseInfo.getStatus(i) == SingleComponentInfo::EQ_SPACES_NO_AUX ||
- traverseInfo.getStatus(i) == SingleComponentInfo::EQ_SPACES_WITH_AUX) {
+ if (assembleMatrix && matrix->getBoundaryManager())
+ matrix->getBoundaryManager()->exitMatrix(matrix);
+
+ if (matrix)
+ nnz += matrix->getBaseMatrix().nnz();
+ }
- if (writeAsmInfo)
- std::cout << "ASM on one mesh!" << std::endl;
+ // And now assemble boundary conditions on the vectors
+ assembleBoundaryConditions(rhs->getDOFVector(i),
+ solution->getDOFVector(i),
+ componentMeshes[i],
+ assembleFlag);
+ }
- // 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],
- assembleFlag,
- assembleMatrix ? matrix : NULL,
- ((i == j) && asmVector) ? rhs->getDOFVector(i) : NULL);
+ solverMatrix.setMatrix(*systemMatrix);
- } else if (traverseInfo.getStatus(i) == SingleComponentInfo::EQ_SPACES_WITH_DIF_AUX) {
+ createPrecon();
- // Row fe space and col fe space are both equal, but right hand side has at
- // least one another aux fe space.
+ INFO(info, 8)("fillin of assembled matrix: %d\n", nnz);
- if (writeAsmInfo)
- std::cout << "ASM on one mesh!" << std::endl;
+#ifdef _OPENMP
+ INFO(info, 8)("buildAfterCoarsen needed %.5f seconds system time / %.5f seconds wallclock time\n",
+ TIME_USED(first, clock()), omp_get_wtime() - wtime);
+#else
+ INFO(info, 8)("buildAfterCoarsen needed %.5f seconds\n",
+ TIME_USED(first, clock()));
+#endif
+ }
- assembleOnOneMesh(componentSpaces[i],
- assembleFlag,
- assembleMatrix ? matrix : NULL,
- ((i == j) && asmVector) ? rhs->getDOFVector(i) : NULL);
-
- 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");
+
+ bool ProblemVec::dualMeshTraverseRequired()
+ {
+ FUNCNAME("ProblemVec::dualMeshTraverseRequired()");
+
+ TEST_EXIT(meshes.size() <= 2)("More than two mneshes are not supported yet!\n");
+
+ return (meshes.size() == 2);
+ }
+
+
+ void ProblemVec::dualAssemble(AdaptInfo *adaptInfo, Flag flag)
+ {
+ FUNCNAME("ProblemVec::dualAssemble()");
+
+ clock_t first = clock();
+#ifdef _OPENMP
+ double wtime = omp_get_wtime();
+#endif
+
+ for (int i = 0; i < static_cast<int>(meshes.size()); i++)
+ meshes[i]->dofCompress();
+
+ Flag assembleFlag =
+ flag |
+ (*systemMatrix)[0][0]->getAssembleFlag() |
+ rhs->getDOFVector(0)->getAssembleFlag() |
+ Mesh::CALL_LEAF_EL |
+ Mesh::FILL_COORDS |
+ Mesh::FILL_DET |
+ Mesh::FILL_GRD_LAMBDA |
+ Mesh::FILL_NEIGH;
+
+ if (useGetBound)
+ assembleFlag |= Mesh::FILL_BOUND;
+
+ traverseInfo.updateStatus();
+
+ // Used to calculate the overall number of non zero entries.
+ int nnz = 0;
+
+ for (int i = 0; i < nComponents; i++) {
+
+ MSG("%d DOFs for %s\n",
+ componentSpaces[i]->getAdmin()->getUsedSize(),
+ componentSpaces[i]->getName().c_str());
+
+ rhs->getDOFVector(i)->set(0.0);
+
+ for (int j = 0; j < nComponents; j++) {
+
+ if (writeAsmInfo)
+ std::cout << "-------" << i << " " << j << "----------------" << std::endl;
+
+ // Only if this variable is true, the current matrix will be assembled.
+ bool assembleMatrix = true;
+ // The DOFMatrix which should be assembled (or not, if assembleMatrix
+ // will be set to false).
+ DOFMatrix *matrix = (*systemMatrix)[i][j];
+
+ if (writeAsmInfo && matrix) {
+ for (std::vector<Operator*>::iterator it = matrix->getOperatorsBegin();
+ it != matrix->getOperatorsEnd(); ++it) {
+ Assembler *assembler = (*it)->getAssembler();
+ if (assembler->getZeroOrderAssembler())
+ std::cout << "ZOA: " << assembler->getZeroOrderAssembler()->getName() << std::endl;
+ if (assembler->getFirstOrderAssembler(GRD_PSI))
+ std::cout << "FOA GRD_PSI: " << assembler->getFirstOrderAssembler(GRD_PSI)->getName() << std::endl;
+ if (assembler->getFirstOrderAssembler(GRD_PHI))
+ std::cout << "FOA GRD_PHI: " << assembler->getFirstOrderAssembler(GRD_PHI)->getName() << std::endl;
+ if (assembler->getSecondOrderAssembler())
+ std::cout << "SOA: " << assembler->getSecondOrderAssembler()->getName() << std::endl;
}
+ }
- } else if (traverseInfo.getStatus(i, j) == SingleComponentInfo::EQ_SPACES_WITH_DIF_AUX) {
- if (writeAsmInfo)
- std::cout << "ASM on one mesh!" << std::endl;
+ if (matrix)
+ matrix->calculateNnz();
+
+ // If the matrix was assembled before and it is marked to be assembled
+ // only once, it will not be assembled.
+ if (assembleMatrixOnlyOnce[i][j] && assembledMatrix[i][j]) {
+ assembleMatrix = false;
+ } else if (matrix) {
+ matrix->getBaseMatrix().
+ change_dim(componentSpaces[i]->getAdmin()->getUsedSize(),
+ componentSpaces[j]->getAdmin()->getUsedSize());
- assembleOnOneMesh(componentSpaces[i],
- assembleFlag,
- assembleMatrix ? matrix : NULL,
- ((i == j) && asmVector) ? rhs->getDOFVector(i) : NULL);
+ set_to_zero(matrix->getBaseMatrix());
+ matrix->startInsertion(matrix->getNnz());
- if (writeAsmInfo)
- std::cout << "ASM on two meshes!" << std::endl;
+ if (matrix->getBoundaryManager())
+ matrix->getBoundaryManager()->initMatrix(matrix);
+ }
- assembleOnDifMeshes2(componentSpaces[i],
- traverseInfo.getAuxFeSpace(i, j),
- assembleFlag,
- assembleMatrix ? matrix : NULL,
- ((i == j) && asmVector) ? rhs->getDOFVector(i) : NULL);
+ // If there is no DOFMatrix, e.g., if it is completly 0, do not assemble.
+ if (!matrix || !assembleMatrix)
+ assembleMatrix = false;
- } else if (traverseInfo.getStatus(i, j) == SingleComponentInfo::DIF_SPACES_NO_AUX ||
- traverseInfo.getStatus(i, j) == SingleComponentInfo::DIF_SPACES_WITH_AUX) {
+ // If the matrix should not be assembled, the rhs vector has to be considered.
+ // This will be only done, if i == j. So, if both is not true, we can jump
+ // to the next matrix.
+ if (!assembleMatrix && i != j)
+ if (matrix)
+ nnz += matrix->getBaseMatrix().nnz();
+ }
+ }
- if (writeAsmInfo)
- std::cout << "ASM on two meshes!" << std::endl;
+ TEST_EXIT(meshes.size() == 2)("There is something wrong!\n");
- assembleOnDifMeshes(componentSpaces[i], componentSpaces[j],
- assembleFlag,
- assembleMatrix ? matrix : NULL,
- ((i == j) && asmVector) ? rhs->getDOFVector(i) : NULL);
- } else {
- ERROR_EXIT("Not yet implemented!\n");
+ const BasisFunction *basisFcts = componentSpaces[0]->getBasisFcts();
+ BoundaryType *bound =
+ useGetBound ? new BoundaryType[basisFcts->getNumber()] : NULL;
+
+ DualTraverse dualTraverse;
+ DualElInfo dualElInfo;
+ int oldElIndex0 = -1;
+ int oldElIndex1 = -1;
+ dualTraverse.setFillSubElemMat(true, basisFcts);
+ bool cont = dualTraverse.traverseFirst(meshes[0], meshes[1], -1, -1,
+ assembleFlag, assembleFlag, dualElInfo);
+
+ while (cont) {
+ bool newEl0 = (dualElInfo.rowElInfo->getElement()->getIndex() != oldElIndex0);
+ bool newEl1 = (dualElInfo.colElInfo->getElement()->getIndex() != oldElIndex1);
+ oldElIndex0 = dualElInfo.rowElInfo->getElement()->getIndex();
+ oldElIndex1 = dualElInfo.colElInfo->getElement()->getIndex();
+
+ for (int i = 0; i < nComponents; i++) {
+ for (int j = 0; j < nComponents; j++) {
+ DOFMatrix *matrix = (*systemMatrix)[i][j];
+
+ if (!matrix)
+ continue;
+
+ if (traverseInfo.eqSpaces(i, j)) {
+ ElInfo *elInfo = NULL;
+ if (componentMeshes[i] == meshes[0] && newEl0)
+ elInfo = dualElInfo.rowElInfo;
+ if (componentMeshes[i] == meshes[1] && newEl1)
+ elInfo = dualElInfo.colElInfo;
+
+ if (elInfo != NULL) {
+ if (useGetBound)
+ basisFcts->getBound(elInfo, bound);
+
+ if (matrix) {
+ matrix->assemble(1.0, elInfo, bound);
+ // Take the matrix boundary manager from the public matrix,
+ // but assemble the boundary conditions on the thread private matrix.
+ if (matrix->getBoundaryManager())
+ matrix->getBoundaryManager()->fillBoundaryConditions(elInfo, matrix);
+ }
+
+ if (i == j)
+ rhs->getDOFVector(i)->assemble(1.0, elInfo, bound);
+ }
+
+ if (traverseInfo.difAuxSpace(i) && i == j) {
+ ElInfo *mainElInfo, *auxElInfo;
+ if (traverseInfo.getRowFeSpace(i)->getMesh() == meshes[0]) {
+ mainElInfo = dualElInfo.rowElInfo;
+ auxElInfo = dualElInfo.colElInfo;
+ } else {
+ mainElInfo = dualElInfo.colElInfo;
+ auxElInfo = dualElInfo.rowElInfo;
+ }
+
+ if (useGetBound)
+ basisFcts->getBound(mainElInfo, bound);
+
+ rhs->getDOFVector(i)->assemble2(1.0, mainElInfo, auxElInfo,
+ dualElInfo.smallElInfo, dualElInfo.largeElInfo, bound);
+ }
+
+ if (traverseInfo.difAuxSpace(i, j) && matrix) {
+ ElInfo *mainElInfo, *auxElInfo;
+ if (traverseInfo.getRowFeSpace(i)->getMesh() == meshes[0]) {
+ mainElInfo = dualElInfo.rowElInfo;
+ auxElInfo = dualElInfo.colElInfo;
+ } else {
+ mainElInfo = dualElInfo.colElInfo;
+ auxElInfo = dualElInfo.rowElInfo;
+ }
+
+ if (useGetBound)
+ basisFcts->getBound(mainElInfo, bound);
+
+ matrix->assemble2(1.0, mainElInfo, auxElInfo,
+ dualElInfo.smallElInfo, dualElInfo.largeElInfo, bound);
+
+ if (matrix->getBoundaryManager())
+ matrix->getBoundaryManager()->fillBoundaryConditions(mainElInfo, matrix);
+ }
+ } else {
+ TEST_EXIT_DBG(traverseInfo.getStatus(i, j) != SingleComponentInfo::DIF_SPACES_WITH_DIF_AUX)
+ ("Not yet supported!\n");
+
+ ElInfo *rowElInfo, *colElInfo;
+ if (componentMeshes[i] == meshes[0]) {
+ rowElInfo = dualElInfo.rowElInfo;
+ colElInfo = dualElInfo.colElInfo;
+ } else {
+ rowElInfo = dualElInfo.colElInfo;
+ colElInfo = dualElInfo.rowElInfo;
+ }
+
+ if (useGetBound)
+ basisFcts->getBound(rowElInfo, bound);
+
+ if (matrix) {
+ matrix->assemble(1.0, rowElInfo, colElInfo,
+ dualElInfo.smallElInfo, dualElInfo.largeElInfo, bound);
+
+ if (matrix->getBoundaryManager())
+ matrix->getBoundaryManager()->fillBoundaryConditions(rowElInfo, matrix);
+ }
+
+ if (i == j)
+ rhs->getDOFVector(i)->assemble(1.0, rowElInfo, bound);
+
+ }
}
+ }
- assembledMatrix[i][j] = true;
+ cont = dualTraverse.traverseNext(dualElInfo);
+ }
- if (assembleMatrix)
+
+ for (int i = 0; i < nComponents; i++) {
+ for (int j = 0; j < nComponents; j++) {
+ DOFMatrix *matrix = (*systemMatrix)[i][j];
+
+ if (matrix)
+ matrix->removeRowsWithDBC(matrix->getApplyDBCs());
+
+ if (matrix)
matrix->finishInsertion();
- if (assembleMatrix && matrix->getBoundaryManager())
+ if (matrix && matrix->getBoundaryManager())
matrix->getBoundaryManager()->exitMatrix(matrix);
-
+
if (matrix)
- nnz += matrix->getBaseMatrix().nnz();
+ nnz += matrix->getBaseMatrix().nnz();
}
// And now assemble boundary conditions on the vectors
@@ -1191,112 +1407,6 @@ namespace AMDiS {
}
- void ProblemVec::assembleOnDifMeshes(FiniteElemSpace *rowFeSpace,
- FiniteElemSpace *colFeSpace,
- Flag assembleFlag,
- DOFMatrix *matrix,
- DOFVector<double> *vector)
- {
- FUNCNAME("ProblemVec::assembleOnDifMeshes()");
-
- const BasisFunction *basisFcts = rowFeSpace->getBasisFcts();
- BoundaryType *bound = NULL;
- if (useGetBound)
- bound = new BoundaryType[basisFcts->getNumber()];
-
- if (matrix)
- matrix->startInsertion();
-
- DualTraverse dualTraverse;
- ElInfo *rowElInfo, *colElInfo;
- ElInfo *largeElInfo, *smallElInfo;
-
- dualTraverse.setFillSubElemMat(true, basisFcts);
- bool cont = dualTraverse.traverseFirst(rowFeSpace->getMesh(),
- colFeSpace->getMesh(),
- -1, -1,
- assembleFlag, assembleFlag,
- &rowElInfo, &colElInfo,
- &smallElInfo, &largeElInfo);
- while (cont) {
- if (useGetBound)
- basisFcts->getBound(rowElInfo, bound);
-
- if (matrix) {
- matrix->assemble(1.0, rowElInfo, colElInfo, smallElInfo, largeElInfo, bound);
-
- if (matrix->getBoundaryManager())
- matrix->getBoundaryManager()->fillBoundaryConditions(rowElInfo, matrix);
- }
-
- if (vector)
- vector->assemble(1.0, rowElInfo, bound);
-
- cont = dualTraverse.traverseNext(&rowElInfo, &colElInfo,
- &smallElInfo, &largeElInfo);
- }
-
- if (matrix)
- matrix->removeRowsWithDBC(matrix->getApplyDBCs());
-
- if (useGetBound)
- delete [] bound;
- }
-
-
- void ProblemVec::assembleOnDifMeshes2(const FiniteElemSpace *mainFeSpace,
- const FiniteElemSpace *auxFeSpace,
- Flag assembleFlag,
- DOFMatrix *matrix,
- DOFVector<double> *vector)
- {
- FUNCNAME("ProblemVec::assembleOnDifMeshes2()");
-
- TEST_EXIT(mainFeSpace)("No main FE space!\n");
- TEST_EXIT(auxFeSpace)("No aux FE space!\n");
-
- Mesh *mainMesh = mainFeSpace->getMesh();
- Mesh *auxMesh = auxFeSpace->getMesh();
-
- const BasisFunction *basisFcts = mainFeSpace->getBasisFcts();
- BoundaryType *bound = NULL;
- if (useGetBound)
- bound = new BoundaryType[basisFcts->getNumber()];
-
- if (matrix)
- matrix->startInsertion();
-
- DualTraverse dualTraverse;
- ElInfo *mainElInfo, *auxElInfo;
- ElInfo *largeElInfo, *smallElInfo;
-
- dualTraverse.setFillSubElemMat(true, basisFcts);
- bool cont = dualTraverse.traverseFirst(mainMesh, auxMesh, -1, -1,
- assembleFlag, assembleFlag,
- &mainElInfo, &auxElInfo,
- &smallElInfo, &largeElInfo);
- while (cont) {
- if (useGetBound)
- basisFcts->getBound(mainElInfo, bound);
-
- if (matrix)
- matrix->assemble2(1.0, mainElInfo, auxElInfo, smallElInfo, largeElInfo, bound);
-
- if (vector)
- vector->assemble2(1.0, mainElInfo, auxElInfo, smallElInfo, largeElInfo, bound);
-
- cont = dualTraverse.traverseNext(&mainElInfo, &auxElInfo,
- &smallElInfo, &largeElInfo);
- }
-
- if (matrix)
- matrix->removeRowsWithDBC(matrix->getApplyDBCs());
-
- if (useGetBound)
- delete [] bound;
- }
-
-
void ProblemVec::assembleBoundaryConditions(DOFVector<double> *rhs,
DOFVector<double> *solution,
Mesh *mesh,
@@ -1364,7 +1474,7 @@ namespace AMDiS {
ElInfo *elInfo = stack.traverseFirst(mesh, -1, assembleFlag);
while (elInfo) {
if (rhs->getBoundaryManager())
- rhs->getBoundaryManager()-> fillBoundaryConditions(elInfo, rhs);
+ rhs->getBoundaryManager()->fillBoundaryConditions(elInfo, rhs);
if (solution->getBoundaryManager())
solution->getBoundaryManager()->fillBoundaryConditions(elInfo, solution);
diff --git a/AMDiS/src/ProblemVec.h b/AMDiS/src/ProblemVec.h
index 3017c89d47e71cb42f3b63bfebfcb70db60201d1..9aac3317410c068c354b57915198dde8c1255de8 100644
--- a/AMDiS/src/ProblemVec.h
+++ b/AMDiS/src/ProblemVec.h
@@ -180,6 +180,10 @@ namespace AMDiS {
bool assembleMatrix = true,
bool assembleVector = true);
+ bool dualMeshTraverseRequired();
+
+ void dualAssemble(AdaptInfo *adaptInfo, Flag flag);
+
void createPrecon();
@@ -266,19 +270,6 @@ namespace AMDiS {
Flag assembleFlag,
DOFMatrix *matrix, DOFVector<double> *vector);
- ///
- void assembleOnDifMeshes(FiniteElemSpace *rowFeSpace,
- FiniteElemSpace *colFeSpace,
- Flag assembleFlag,
- DOFMatrix *matrix,
- DOFVector<double> *vector);
-
- ///
- void assembleOnDifMeshes2(const FiniteElemSpace *mainFeSpace,
- const FiniteElemSpace *auxFeSpace,
- Flag assembleFlag,
- DOFMatrix *matrix,
- DOFVector<double> *vector);
///
void assembleBoundaryConditions(DOFVector<double> *rhs,
diff --git a/AMDiS/src/SecondOrderAssembler.cc b/AMDiS/src/SecondOrderAssembler.cc
index 0d6fd931da0e0732abb9e54722cbd039cb6b0718..2bd7858f14170760ff7b843c59984f43f33d2096 100644
--- a/AMDiS/src/SecondOrderAssembler.cc
+++ b/AMDiS/src/SecondOrderAssembler.cc
@@ -21,11 +21,11 @@ namespace AMDiS {
: SubAssembler(op, assembler, quad, 2, optimized)
{}
- SecondOrderAssembler*
- SecondOrderAssembler::getSubAssembler(Operator* op,
- Assembler *assembler,
- Quadrature *quad,
- bool optimized)
+
+ SecondOrderAssembler* SecondOrderAssembler::getSubAssembler(Operator* op,
+ Assembler *assembler,
+ Quadrature *quad,
+ bool optimized)
{
int myRank = omp_get_thread_num();
@@ -78,9 +78,12 @@ namespace AMDiS {
return newAssembler;
}
+
Pre2::Pre2(Operator *op, Assembler *assembler, Quadrature *quad)
: SecondOrderAssembler(op, assembler, quad, true)
{
+ name = "precalculated second order assembler";
+
q11 = Q11PsiPhi::provideQ11PsiPhi(owner->getRowFESpace()->getBasisFcts(),
owner->getColFESpace()->getBasisFcts(),
quadrature);
@@ -91,6 +94,7 @@ namespace AMDiS {
}
}
+
Pre2::~Pre2()
{
for (int i = 0; i < static_cast<int>(tmpLALt.size()); i++) {
@@ -99,6 +103,7 @@ namespace AMDiS {
}
}
+
void Pre2::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
{
const int **nEntries;
@@ -158,13 +163,17 @@ namespace AMDiS {
}
}
+
Quad2::Quad2(Operator *op, Assembler *assembler, Quadrature *quad)
: SecondOrderAssembler(op, assembler, quad, true)
{
+ name = "fast quadrature second order assembler";
+
tmpVec.resize(omp_get_overall_max_threads());
tmpLALt.resize(omp_get_overall_max_threads());
}
+
Quad2::~Quad2()
{
if (!firstCall) {
@@ -178,6 +187,7 @@ namespace AMDiS {
}
}
+
void Quad2::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
{
int nPoints = quadrature->getNumPoints();
@@ -247,9 +257,13 @@ namespace AMDiS {
}
}
+
Stand2::Stand2(Operator *op, Assembler *assembler, Quadrature *quad)
: SecondOrderAssembler(op, assembler, quad, false)
- {}
+ {
+ name = "standard second order assembler";
+ }
+
void Stand2::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
{
diff --git a/AMDiS/src/SubAssembler.cc b/AMDiS/src/SubAssembler.cc
index e780d6e2a6f404f943a09c72b47443f9b2b33341..2e41ca904073685a397688e2e27561497854f83f 100644
--- a/AMDiS/src/SubAssembler.cc
+++ b/AMDiS/src/SubAssembler.cc
@@ -26,7 +26,8 @@ namespace AMDiS {
owner(assembler),
symmetric(true),
opt(optimized),
- firstCall(true)
+ firstCall(true),
+ name("")
{
const BasisFunction *psi = assembler->rowFESpace->getBasisFcts();
const BasisFunction *phi = assembler->colFESpace->getBasisFcts();
diff --git a/AMDiS/src/SubAssembler.h b/AMDiS/src/SubAssembler.h
index 4f5ae9bedf832350b821858d2e403d4a876f7f63..8e66d1ab1b433f8f84304c58ddf5e53ff919a8af 100644
--- a/AMDiS/src/SubAssembler.h
+++ b/AMDiS/src/SubAssembler.h
@@ -131,6 +131,12 @@ namespace AMDiS {
return phiFast;
}
+ /// Returns \ref name.
+ std::string getName() const
+ {
+ return name;
+ }
+
protected:
/// Updates \ref psiFast and \ref phiFast.
FastQuadrature *updateFastQuadrature(FastQuadrature *quadFast,
@@ -209,6 +215,9 @@ namespace AMDiS {
///
bool firstCall;
+ /// Name of the assembler. Is used to print information about used assembler.
+ std::string name;
+
friend class Assembler;
};
diff --git a/AMDiS/src/ZeroOrderAssembler.cc b/AMDiS/src/ZeroOrderAssembler.cc
index 736b07b80d151de2b5220fd90c3bec440f702e53..270b15a4fe91a6fab2d0425474eb99bea0ed493b 100644
--- a/AMDiS/src/ZeroOrderAssembler.cc
+++ b/AMDiS/src/ZeroOrderAssembler.cc
@@ -21,6 +21,7 @@ namespace AMDiS {
: SubAssembler(op, assembler, quad, 0, optimized)
{}
+
ZeroOrderAssembler* ZeroOrderAssembler::getSubAssembler(Operator* op,
Assembler *assembler,
Quadrature *quad,
@@ -28,10 +29,9 @@ namespace AMDiS {
{
int myRank = omp_get_thread_num();
- // check if a assembler is needed at all
- if (op->zeroOrder[myRank].size() == 0) {
+ // check if an assembler is needed at all
+ if (op->zeroOrder[myRank].size() == 0)
return NULL;
- }
ZeroOrderAssembler *newAssembler;
@@ -77,9 +77,13 @@ namespace AMDiS {
return newAssembler;
}
+
StandardZOA::StandardZOA(Operator *op, Assembler *assembler, Quadrature *quad)
- : ZeroOrderAssembler(op, assembler, quad, false)
- {}
+ : ZeroOrderAssembler(op, assembler, quad, false)
+ {
+ name = "standard zero order assembler";
+ }
+
void StandardZOA::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
{
@@ -132,6 +136,7 @@ namespace AMDiS {
}
}
+
void StandardZOA::calculateElementVector(const ElInfo *elInfo, ElementVector& vec)
{
int nPoints = quadrature->getNumPoints();
@@ -153,12 +158,15 @@ namespace AMDiS {
}
}
+
FastQuadZOA::FastQuadZOA(Operator *op, Assembler *assembler, Quadrature *quad)
: ZeroOrderAssembler(op, assembler, quad, true)
{
+ name = "fast quadrature zero order assembler";
tmpC.resize(omp_get_overall_max_threads());
}
+
void FastQuadZOA::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
{
int nPoints = quadrature->getNumPoints();
@@ -217,6 +225,7 @@ namespace AMDiS {
}
}
+
void FastQuadZOA::calculateElementVector(const ElInfo *elInfo, ElementVector& vec)
{
int myRank = omp_get_thread_num();
@@ -251,11 +260,14 @@ namespace AMDiS {
}
}
+
PrecalcZOA::PrecalcZOA(Operator *op, Assembler *assembler, Quadrature *quad)
: ZeroOrderAssembler(op, assembler, quad, true)
{
+ name = "precalculated zero order assembler";
}
+
void PrecalcZOA::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
{
if (firstCall) {
@@ -300,6 +312,7 @@ namespace AMDiS {
}
}
+
void PrecalcZOA::calculateElementVector(const ElInfo *elInfo, ElementVector& vec)
{
if (firstCall) {