diff --git a/AMDiS/src/Assembler.cc b/AMDiS/src/Assembler.cc
index adc1b6423b249b472d5e7d76074d39dd2d7fd41b..3e24e16673168f66d2126658f9052b524de2bca8 100644
--- a/AMDiS/src/Assembler.cc
+++ b/AMDiS/src/Assembler.cc
@@ -299,9 +299,8 @@ namespace AMDiS {
for (int i = 0; i < nBasFcts; i++) {
uhOldLoc2[i] = 0.0;
- for (int j = 0; j < nBasFcts; j++) {
+ for (int j = 0; j < nBasFcts; j++)
uhOldLoc2[i] += m[j][i] * uhOldLoc[i];
- }
}
@@ -312,9 +311,8 @@ namespace AMDiS {
for (int i = 0; i < nBasFcts; i++) {
double val = 0.0;
- for (int j = 0; j < nBasFcts; j++) {
+ for (int j = 0; j < nBasFcts; j++)
val += elementMatrix[i][j] * uhOldLoc2[j];
- }
vec[i] += val;
}
@@ -322,6 +320,7 @@ namespace AMDiS {
delete [] uhOldLoc2;
}
+
void Assembler::initElement(const ElInfo *smallElInfo,
const ElInfo *largeElInfo,
Quadrature *quad)
@@ -336,6 +335,7 @@ namespace AMDiS {
zeroOrderAssembler->initElement(smallElInfo, largeElInfo, quad);
}
+
void Assembler::checkQuadratures()
{
if (secondOrderAssembler) {
diff --git a/AMDiS/src/DualTraverse.cc b/AMDiS/src/DualTraverse.cc
index de4ccd81d86b2e012f4adc93ebe04e385d39d069..709f34c1c945e2ed7801b688bff46342ae7ba61d 100644
--- a/AMDiS/src/DualTraverse.cc
+++ b/AMDiS/src/DualTraverse.cc
@@ -165,6 +165,7 @@ namespace AMDiS {
}
}
+
void DualTraverse::fillSubElInfo(ElInfo *elInfo1,
ElInfo *elInfo2,
ElInfo *elInfoSmall,
@@ -173,19 +174,36 @@ namespace AMDiS {
if (!fillSubElemMat)
return;
- // mtl::dense2D<double>& subCoordsMat =
- // const_cast<mtl::dense2D<double>&>(elInfoSmall->getSubElemCoordsMat());
- // mtl::dense2D<double>& subCoordsMat_so =
- // const_cast<mtl::dense2D<double>&>(elInfoSmall->getSubElemCoordsMat_so());
-
- if (elInfo1 == elInfoSmall) {
- // stack1.getCoordsInElem(elInfo2, basisFcts, subCoordsMat);
- // stack1.getCoordsInElem_so(elInfo2, basisFcts, subCoordsMat_so);
+ if (elInfo1 == elInfoSmall)
stack1.fillRefinementPath(*elInfoSmall, *elInfo2);
- } else {
- // stack2.getCoordsInElem(elInfo1, basisFcts, subCoordsMat);
- // stack2.getCoordsInElem_so(elInfo1, basisFcts, subCoordsMat_so);
- stack2.fillRefinementPath(*elInfoSmall, *elInfo1);
- }
+ else
+ stack2.fillRefinementPath(*elInfoSmall, *elInfo1);
+ }
+
+
+ int DualTraverse::getFace(DualElInfo *dualElInfo, int largeFace)
+ {
+ FUNCNAME("DualTraverse::getFace()");
+
+ TEST_EXIT_DBG(dualElInfo)("No dual element info object!\n");
+
+ ElInfo *largeElInfo = dualElInfo->largeElInfo;
+ ElInfo *smallElInfo = dualElInfo->smallElInfo;
+
+ TEST_EXIT_DBG(largeElInfo->getLevel() <= smallElInfo->getLevel())
+ ("Should not happen!\n");
+
+ if (largeElInfo->getLevel() == smallElInfo->getLevel())
+ return largeFace;
+
+ int refinementPathLength = smallElInfo->getRefinementPathLength();
+ unsigned long refinementPath = smallElInfo->getRefinementPath();
+
+ TEST_EXIT_DBG(refinementPathLength != 0)
+ ("Refinement path is empty. This should not happen!\n");
+
+ std::cout << "REFINEMENTPATHLENGTH = " << refinementPathLength << std::endl;
+
+ return -1;
}
}
diff --git a/AMDiS/src/DualTraverse.h b/AMDiS/src/DualTraverse.h
index 0c4059b2039bf2fbe7873325a4b58e9d65fb4c31..8226414014cc2828246ae0517b84eea2e475b3ca 100644
--- a/AMDiS/src/DualTraverse.h
+++ b/AMDiS/src/DualTraverse.h
@@ -115,6 +115,16 @@ namespace AMDiS {
fillSubElemMat = b;
basisFcts = fcts;
}
+
+ /** \brief
+ * Checks if the small element has an edge/face which is part of a given edge/face
+ * of the large element. If this is the case, it returns the local number of the
+ * small edge/face, and -1 otherwise.
+ *
+ * \param[in] dualElInfo Dual element info with large and small element infos.
+ * \param[in] largeFace A local edge/face number on the large element.
+ */
+ static int getFace(DualElInfo *dualElInfo, int largeFace);
protected:
/** \brief
diff --git a/AMDiS/src/ElInfo.h b/AMDiS/src/ElInfo.h
index b61dc357befd0d580a996c230169090cef728cc4..2b75dcfd335911ec5c42af923d04c941af293519 100644
--- a/AMDiS/src/ElInfo.h
+++ b/AMDiS/src/ElInfo.h
@@ -227,6 +227,19 @@ namespace AMDiS {
return parametric;
}
+ /// Returns \ref refinementPath
+ inline unsigned long getRefinementPath() const
+ {
+ return refinementPath;
+ }
+
+ /// Get \ref refinementPathLength
+ inline int getRefinementPathLength() const
+ {
+ return refinementPathLength;
+ }
+
+
virtual void getSubElemCoordsMat(mtl::dense2D<double>& mat,
int basisFctsDegree) const {}
diff --git a/AMDiS/src/ElInfo2d.cc b/AMDiS/src/ElInfo2d.cc
index 9a837edc5a99bfc25b706fa09e30f4915fbf1878..ae87da7a8ae1b82d9f5dc7a5af90f737c7900a7b 100644
--- a/AMDiS/src/ElInfo2d.cc
+++ b/AMDiS/src/ElInfo2d.cc
@@ -716,12 +716,9 @@ namespace AMDiS {
if (refinementPath & (1 << i)) {
tmpMat = mat_d1_right * mat;
mat = tmpMat;
- // mat *= mat_d1_right;
} else {
tmpMat = mat_d1_left * mat;
mat = tmpMat;
-
- // mat *= mat_d1_left;
}
}
}
diff --git a/AMDiS/src/Operator.cc b/AMDiS/src/Operator.cc
index 7d8136d79c1135522b13c07b5b6cc89df006ccf8..9d3bd3d74291062b7e60cbddf13e939ba8434729 100644
--- a/AMDiS/src/Operator.cc
+++ b/AMDiS/src/Operator.cc
@@ -132,6 +132,8 @@ namespace AMDiS {
{
FUNCNAME("OperatorTerm::getGradientsAtQPs()");
+ ERROR_EXIT("Not yet tested!\n");
+
TEST_EXIT(smallElInfo->getMesh() == vec->getFESpace()->getMesh() ||
largeElInfo->getMesh() == vec->getFESpace()->getMesh())
("There is something wrong!\n");
@@ -1026,7 +1028,12 @@ namespace AMDiS {
AbstractFunction<WorldMatrix<double>, WorldVector<double> > *af,
AbstractFunction<WorldVector<double>, WorldMatrix<double> > *divAf,
bool symm)
- : SecondOrderTerm(af->getDegree()), vec(dv), f(af), divFct(divAf), symmetric(symm)
+ : SecondOrderTerm(af->getDegree()),
+ vec(dv),
+ f(af),
+ divFct(divAf),
+ gradAtQPs(NULL),
+ symmetric(symm)
{
setSymmetric(symmetric);
@@ -1163,13 +1170,12 @@ namespace AMDiS {
}
}
- void MatrixVec2_SOT::weakEval(int nPoints,
- const WorldVector<double> *grdUhAtQP,
- WorldVector<double> *result) const
+ void MatrixVec2_SOT::weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
+ std::vector<WorldVector<double> > &result) const
{
- if (grdUhAtQP)
- for (int iq = 0; iq < nPoints; iq++)
- result[iq] += A * grdUhAtQP[iq] * (*fct)(vec1AtQPs[iq], vec2AtQPs[iq]);
+ int nPoints = grdUhAtQP.size();
+ for (int iq = 0; iq < nPoints; iq++)
+ result[iq] += A * grdUhAtQP[iq] * (*fct)(vec1AtQPs[iq], vec2AtQPs[iq]);
}
General_SOT::General_SOT(std::vector<DOFVectorBase<double>*> vecs,
@@ -1302,6 +1308,14 @@ namespace AMDiS {
gradAtQPs = getGradientsAtQPs(vec, elInfo, subAssembler, quad);
}
+ void MatrixGradient_SOT::initElement(const ElInfo* smallElInfo,
+ const ElInfo* largeElInfo,
+ SubAssembler* subAssembler,
+ Quadrature *quad)
+ {
+ gradAtQPs = getGradientsAtQPs(vec, smallElInfo, largeElInfo, subAssembler, quad);
+ }
+
void FctGradient_SOT::initElement(const ElInfo* elInfo,
SubAssembler* subAssembler,
Quadrature *quad)
@@ -1944,17 +1958,15 @@ namespace AMDiS {
}
}
- void MatrixFct_SOT::weakEval(int nPoints,
- const WorldVector<double> *grdUhAtQP,
- WorldVector<double> *result) const
+ void MatrixFct_SOT::weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
+ std::vector<WorldVector<double> > &result) const
{
- if (grdUhAtQP) {
- WorldMatrix<double> A;
- for (int iq = 0; iq < nPoints; iq++) {
- A = (*matrixFct)(vecAtQPs[iq]);
- result[iq] += A * grdUhAtQP[iq];
- }
- }
+ int nPoints = grdUhAtQP.size();
+ WorldMatrix<double> A;
+ for (int iq = 0; iq < nPoints; iq++) {
+ A = (*matrixFct)(vecAtQPs[iq]);
+ result[iq] += A * grdUhAtQP[iq];
+ }
}
@@ -1979,13 +1991,12 @@ namespace AMDiS {
}
}
- void Matrix_SOT::weakEval(int nPoints,
- const WorldVector<double> *grdUhAtQP,
- WorldVector<double> *result) const
+ void Matrix_SOT::weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
+ std::vector<WorldVector<double> > &result) const
{
- if (grdUhAtQP)
- for (int iq = 0; iq < nPoints; iq++)
- result[iq] += matrix * grdUhAtQP[iq];
+ int nPoints = grdUhAtQP.size();
+ for (int iq = 0; iq < nPoints; iq++)
+ result[iq] += matrix * grdUhAtQP[iq];
}
void MatrixGradient_SOT::eval(int nPoints,
@@ -2014,19 +2025,24 @@ namespace AMDiS {
}
}
- void MatrixGradient_SOT::weakEval(int nPoints,
- const WorldVector<double> *grdUhAtQP,
- WorldVector<double> *result) const
+
+ void MatrixGradient_SOT::weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
+ std::vector<WorldVector<double> > &result) const
{
- if (grdUhAtQP) {
- WorldMatrix<double> A;
- for (int iq = 0; iq < nPoints; iq++) {
- A = (*f)(gradAtQPs[iq]);
- result[iq] += A * grdUhAtQP[iq];
- }
+ FUNCNAME("MatrixGradient_SOT::weakEval()");
+
+ TEST_EXIT_DBG(f)("No function f!\n");
+ 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];
}
}
+
void VecAtQP_SOT::eval(int nPoints,
const double *uhAtQP,
const WorldVector<double> *grdUhAtQP,
@@ -2047,15 +2063,13 @@ namespace AMDiS {
}
}
- void VecAtQP_SOT::weakEval(int nPoints,
- const WorldVector<double> *grdUhAtQP,
- WorldVector<double> *result) const
+ void VecAtQP_SOT::weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
+ std::vector<WorldVector<double> > &result) const
{
- if (grdUhAtQP) {
- for (int iq = 0; iq < nPoints; iq++) {
- double factor = (*f)(vecAtQPs[iq]);
- axpy(factor, grdUhAtQP[iq], result[iq]);
- }
+ int nPoints = grdUhAtQP.size();
+ for (int iq = 0; iq < nPoints; iq++) {
+ double factor = (*f)(vecAtQPs[iq]);
+ axpy(factor, grdUhAtQP[iq], result[iq]);
}
}
@@ -2079,15 +2093,13 @@ namespace AMDiS {
}
}
- void Vec2AtQP_SOT::weakEval(int nPoints,
- const WorldVector<double> *grdUhAtQP,
- WorldVector<double> *result) const
+ void Vec2AtQP_SOT::weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
+ std::vector<WorldVector<double> > &result) const
{
- if (grdUhAtQP) {
- for (int iq = 0; iq < nPoints; iq++) {
- double factor = (*f)(vecAtQPs1[iq], vecAtQPs2[iq]);
- axpy(factor, grdUhAtQP[iq], result[iq]);
- }
+ int nPoints = grdUhAtQP.size();
+ for (int iq = 0; iq < nPoints; iq++) {
+ double factor = (*f)(vecAtQPs1[iq], vecAtQPs2[iq]);
+ axpy(factor, grdUhAtQP[iq], result[iq]);
}
}
@@ -2111,15 +2123,13 @@ namespace AMDiS {
}
}
- void CoordsAtQP_SOT::weakEval(int nPoints,
- const WorldVector<double> *grdUhAtQP,
- WorldVector<double> *result) const
+ void CoordsAtQP_SOT::weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
+ std::vector<WorldVector<double> > &result) const
{
- if (grdUhAtQP) {
- for (int iq = 0; iq < nPoints; iq++) {
- double factor = (*g)(coordsAtQPs[iq]);
- axpy(factor, grdUhAtQP[iq], result[iq]);
- }
+ int nPoints = grdUhAtQP.size();
+ for (int iq = 0; iq < nPoints; iq++) {
+ double factor = (*g)(coordsAtQPs[iq]);
+ axpy(factor, grdUhAtQP[iq], result[iq]);
}
}
@@ -2143,15 +2153,13 @@ namespace AMDiS {
}
}
- void FctGradient_SOT::weakEval(int nPoints,
- const WorldVector<double> *grdUhAtQP,
- WorldVector<double> *result) const
+ void FctGradient_SOT::weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
+ std::vector<WorldVector<double> > &result) const
{
- if (grdUhAtQP) {
- for (int iq = 0; iq < nPoints; iq++) {
- double factor = (*f)(gradAtQPs[iq]);
- axpy(factor, grdUhAtQP[iq], result[iq]);
- }
+ int nPoints = grdUhAtQP.size();
+ for (int iq = 0; iq < nPoints; iq++) {
+ double factor = (*f)(gradAtQPs[iq]);
+ axpy(factor, grdUhAtQP[iq], result[iq]);
}
}
@@ -2175,15 +2183,13 @@ namespace AMDiS {
}
}
- void VecAndGradientAtQP_SOT::weakEval(int nPoints,
- const WorldVector<double> *grdUhAtQP,
- WorldVector<double> *result) const
+ void VecAndGradientAtQP_SOT::weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
+ std::vector<WorldVector<double> > &result) const
{
- if (grdUhAtQP) {
- for (int iq = 0; iq < nPoints; iq++) {
- double factor = (*f)(vecAtQPs[iq], gradAtQPs[iq]);
- axpy(factor, grdUhAtQP[iq], result[iq]);
- }
+ int nPoints = grdUhAtQP.size();
+ for (int iq = 0; iq < nPoints; iq++) {
+ double factor = (*f)(vecAtQPs[iq], gradAtQPs[iq]);
+ axpy(factor, grdUhAtQP[iq], result[iq]);
}
}
@@ -2213,16 +2219,14 @@ namespace AMDiS {
}
}
- void VecMatrixGradientAtQP_SOT::weakEval(int nPoints,
- const WorldVector<double> *grdUhAtQP,
- WorldVector<double> *result) const
+ void VecMatrixGradientAtQP_SOT::weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
+ std::vector<WorldVector<double> > &result) const
{
- if (grdUhAtQP) {
- WorldMatrix<double> A;
- for (int iq = 0; iq < nPoints; iq++) {
- A = (*f)(vecAtQPs[iq], gradAtQPs[iq]);
- result[iq] += A * grdUhAtQP[iq];
- }
+ int nPoints = grdUhAtQP.size();
+ WorldMatrix<double> A;
+ for (int iq = 0; iq < nPoints; iq++) {
+ A = (*f)(vecAtQPs[iq], gradAtQPs[iq]);
+ result[iq] += A * grdUhAtQP[iq];
}
}
@@ -2246,15 +2250,13 @@ namespace AMDiS {
}
}
- void VecAndCoordsAtQP_SOT::weakEval(int nPoints,
- const WorldVector<double> *grdUhAtQP,
- WorldVector<double> *result) const
+ void VecAndCoordsAtQP_SOT::weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
+ std::vector<WorldVector<double> > &result) const
{
- if (grdUhAtQP) {
- for (int iq = 0; iq < nPoints; iq++) {
- double factor = (*f)(vecAtQPs[iq], coordsAtQPs[iq]);
- axpy(factor, grdUhAtQP[iq], result[iq]);
- }
+ int nPoints = grdUhAtQP.size();
+ for (int iq = 0; iq < nPoints; iq++) {
+ double factor = (*f)(vecAtQPs[iq], coordsAtQPs[iq]);
+ axpy(factor, grdUhAtQP[iq], result[iq]);
}
}
@@ -2284,16 +2286,14 @@ namespace AMDiS {
}
}
- void MatrixGradientAndCoords_SOT::weakEval(int nPoints,
- const WorldVector<double> *grdUhAtQP,
- WorldVector<double> *result) const
+ void MatrixGradientAndCoords_SOT::weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
+ std::vector<WorldVector<double> > &result) const
{
- if (grdUhAtQP) {
- WorldMatrix<double> A;
- for (int iq = 0; iq < nPoints; iq++) {
- A = (*f)(gradAtQPs[iq], coordsAtQPs[iq]);
- result[iq] += A * grdUhAtQP[iq];
- }
+ int nPoints = grdUhAtQP.size();
+ WorldMatrix<double> A;
+ for (int iq = 0; iq < nPoints; iq++) {
+ A = (*f)(gradAtQPs[iq], coordsAtQPs[iq]);
+ result[iq] += A * grdUhAtQP[iq];
}
}
@@ -2317,15 +2317,13 @@ namespace AMDiS {
}
}
- void VecGradCoordsAtQP_SOT::weakEval(int nPoints,
- const WorldVector<double> *grdUhAtQP,
- WorldVector<double> *result) const
+ void VecGradCoordsAtQP_SOT::weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
+ std::vector<WorldVector<double> > &result) const
{
- if (grdUhAtQP) {
- for (int iq = 0; iq < nPoints; iq++) {
- double factor = (*f)(vecAtQPs[iq], gradAtQPs[iq], coordsAtQPs[iq]);
- axpy(factor, grdUhAtQP[iq], result[iq]);
- }
+ int nPoints = grdUhAtQP.size();
+ for (int iq = 0; iq < nPoints; iq++) {
+ double factor = (*f)(vecAtQPs[iq], gradAtQPs[iq], coordsAtQPs[iq]);
+ axpy(factor, grdUhAtQP[iq], result[iq]);
}
}
@@ -2489,15 +2487,13 @@ namespace AMDiS {
}
}
- void CoordsAtQP_IJ_SOT::weakEval(int nPoints,
- const WorldVector<double> *grdUhAtQP,
- WorldVector<double> *result) const
+ void CoordsAtQP_IJ_SOT::weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
+ std::vector<WorldVector<double> > &result) const
{
- if (grdUhAtQP) {
- for (int iq = 0; iq < nPoints; iq++) {
- double factor = (*g)(coordsAtQPs[iq]);
- result[iq][xi] += grdUhAtQP[iq][xj] * factor;
- }
+ int nPoints = grdUhAtQP.size();
+ for (int iq = 0; iq < nPoints; iq++) {
+ double factor = (*g)(coordsAtQPs[iq]);
+ result[iq][xi] += grdUhAtQP[iq][xj] * factor;
}
}
@@ -2533,15 +2529,13 @@ namespace AMDiS {
}
}
- void VecAtQP_IJ_SOT::weakEval(int nPoints,
- const WorldVector<double> *grdUhAtQP,
- WorldVector<double> *result) const
+ void VecAtQP_IJ_SOT::weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
+ std::vector<WorldVector<double> > &result) const
{
- if (grdUhAtQP) {
- for (int iq = 0; iq < nPoints; iq++) {
- double factor = (*f)(vecAtQPs[iq]);
- result[iq][xi] += grdUhAtQP[iq][xj] * factor;
- }
+ int nPoints = grdUhAtQP.size();
+ for (int iq = 0; iq < nPoints; iq++) {
+ double factor = (*f)(vecAtQPs[iq]);
+ result[iq][xi] += grdUhAtQP[iq][xj] * factor;
}
}
@@ -2688,15 +2682,13 @@ namespace AMDiS {
}
}
- void VecAndGradAtQP_SOT::weakEval(int nPoints,
- const WorldVector<double> *grdUhAtQP,
- WorldVector<double> *result) const
+ void VecAndGradAtQP_SOT::weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
+ std::vector<WorldVector<double> > &result) const
{
- if (grdUhAtQP) {
- for (int iq = 0; iq < nPoints; iq++) {
- double factor = (*f)(vecAtQPs[iq], gradAtQPs[iq]);
- axpy(factor, grdUhAtQP[iq], result[iq]);
- }
+ int nPoints = grdUhAtQP.size();
+ for (int iq = 0; iq < nPoints; iq++) {
+ double factor = (*f)(vecAtQPs[iq], gradAtQPs[iq]);
+ axpy(factor, grdUhAtQP[iq], result[iq]);
}
}
@@ -2791,27 +2783,25 @@ namespace AMDiS {
}
}
- void General_SOT::weakEval(int nPoints,
- const WorldVector<double> *grdUhAtQP,
- WorldVector<double> *result) const
+ void General_SOT::weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
+ std::vector<WorldVector<double> > &result) const
{
+ int nPoints = grdUhAtQP.size();
int nVecs = static_cast<int>(vecs_.size());
int nGrads = static_cast<int>(grads_.size());
std::vector<double> vecsArg(nVecs);
std::vector<WorldVector<double> > gradsArg(nGrads);
- if (grdUhAtQP) {
- WorldMatrix<double> A;
- for (int iq = 0; iq < nPoints; iq++) {
- for (int i = 0; i < nVecs; i++)
- vecsArg[i] = vecsAtQPs_[i][iq];
- for (int i = 0; i < nGrads; i++)
- gradsArg[i] = gradsAtQPs_[i][iq];
-
- A = (*f_)(coordsAtQPs_[iq], vecsArg, gradsArg);
- result[iq] += A * grdUhAtQP[iq];
- }
+ WorldMatrix<double> A;
+ for (int iq = 0; iq < nPoints; iq++) {
+ for (int i = 0; i < nVecs; i++)
+ vecsArg[i] = vecsAtQPs_[i][iq];
+ for (int i = 0; i < nGrads; i++)
+ gradsArg[i] = gradsAtQPs_[i][iq];
+
+ A = (*f_)(coordsAtQPs_[iq], vecsArg, gradsArg);
+ result[iq] += A * grdUhAtQP[iq];
}
}
@@ -3015,27 +3005,25 @@ namespace AMDiS {
}
}
- void GeneralParametric_SOT::weakEval(int nPoints,
- const WorldVector<double> *grdUhAtQP,
- WorldVector<double> *result) const
+ void GeneralParametric_SOT::weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
+ std::vector<WorldVector<double> > &result) const
{
+ int nPoints = grdUhAtQP.size();
int nVecs = static_cast<int>(vecs_.size());
int nGrads = static_cast<int>(grads_.size());
std::vector<double> vecsArg(nVecs);
std::vector<WorldVector<double> > gradsArg(nGrads);
- if (grdUhAtQP) {
- WorldMatrix<double> A;
- for (int iq = 0; iq < nPoints; iq++) {
- for (int i = 0; i < nVecs; i++)
- vecsArg[i] = vecsAtQPs_[i][iq];
- for (int i = 0; i < nGrads; i++)
- gradsArg[i] = gradsAtQPs_[i][iq];
-
- A = (*f_)(coordsAtQPs_[iq], elementNormal_, vecsArg, gradsArg);
- result[iq] += A * grdUhAtQP[iq];
- }
+ WorldMatrix<double> A;
+ for (int iq = 0; iq < nPoints; iq++) {
+ for (int i = 0; i < nVecs; i++)
+ vecsArg[i] = vecsAtQPs_[i][iq];
+ for (int i = 0; i < nGrads; i++)
+ gradsArg[i] = gradsAtQPs_[i][iq];
+
+ A = (*f_)(coordsAtQPs_[iq], elementNormal_, vecsArg, gradsArg);
+ result[iq] += A * grdUhAtQP[iq];
}
}
@@ -3269,15 +3257,13 @@ namespace AMDiS {
}
}
- void VecGrad_SOT::weakEval(int nPoints,
- const WorldVector<double> *grdUhAtQP,
- WorldVector<double> *result) const
+ void VecGrad_SOT::weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
+ std::vector<WorldVector<double> > &result) const
{
- if (grdUhAtQP) {
- for (int iq = 0; iq < nPoints; iq++) {
- double factor = (*f)(vecAtQPs[iq], gradAtQPs[iq]);
- axpy(factor, grdUhAtQP[iq], result[iq]);
- }
+ int nPoints = grdUhAtQP.size();
+ for (int iq = 0; iq < nPoints; iq++) {
+ double factor = (*f)(vecAtQPs[iq], gradAtQPs[iq]);
+ axpy(factor, grdUhAtQP[iq], result[iq]);
}
}
diff --git a/AMDiS/src/Operator.h b/AMDiS/src/Operator.h
index 0c3f4cfcac368d402b1017715a8c2c7f6ab3b258..3ed2aa4cbcbbeeb46316657ed3423286ff4b8ef7 100644
--- a/AMDiS/src/Operator.h
+++ b/AMDiS/src/Operator.h
@@ -68,12 +68,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)
- {}
+ 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()
@@ -300,9 +305,8 @@ namespace AMDiS {
DimMat<double> **result) const = 0;
/// Evaluation of \f$ A \nabla u(\vec{x}) \f$ at all quadrature points.
- virtual void weakEval(int nPoints,
- const WorldVector<double> *grdUhAtQP,
- WorldVector<double> *result) const = 0;
+ virtual void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
+ std::vector<WorldVector<double> > &result) const = 0;
};
@@ -352,13 +356,12 @@ namespace AMDiS {
}
/// Implenetation of SecondOrderTerm::weakEval().
- void weakEval(int nPoints,
- const WorldVector<double> *grdUhAtQP,
- WorldVector<double> *result) const
+ void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
+ std::vector<WorldVector<double> > &result) const
{
- if (grdUhAtQP)
- for (int iq = 0; iq < nPoints; iq++)
- result[iq] += grdUhAtQP[iq];
+ int nPoints = grdUhAtQP.size();
+ for (int iq = 0; iq < nPoints; iq++)
+ result[iq] += grdUhAtQP[iq];
}
};
@@ -420,13 +423,12 @@ namespace AMDiS {
}
/// Implenetation of SecondOrderTerm::weakEval().
- void weakEval(int nPoints,
- const WorldVector<double> *grdUhAtQP,
- WorldVector<double> *result) const
+ void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
+ std::vector<WorldVector<double> > &result) const
{
- if (grdUhAtQP)
- for (int iq = 0; iq < nPoints; iq++)
- axpy(*factor, grdUhAtQP[iq], result[iq]);
+ int nPoints = grdUhAtQP.size();
+ for (int iq = 0; iq < nPoints; iq++)
+ axpy(*factor, grdUhAtQP[iq], result[iq]);
}
private:
@@ -467,9 +469,8 @@ namespace AMDiS {
double factor) const;
/// Implenetation of SecondOrderTerm::weakEval().
- void weakEval(int nPoints,
- const WorldVector<double> *grdUhAtQP,
- WorldVector<double> *result) const;
+ void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
+ std::vector<WorldVector<double> > &result) const;
protected:
/// DOFVector to be evaluated at quadrature points.
@@ -522,9 +523,8 @@ namespace AMDiS {
double factor) const;
/// Implenetation of SecondOrderTerm::weakEval().
- void weakEval(int nPoints,
- const WorldVector<double> *grdUhAtQP,
- WorldVector<double> *result) const;
+ void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
+ std::vector<WorldVector<double> > &result) const;
protected:
/// Matrix stroring A.
@@ -584,13 +584,12 @@ namespace AMDiS {
}
/// Implenetation of SecondOrderTerm::weakEval().
- void weakEval(int nPoints,
- const WorldVector<double> *grdUhAtQP,
- WorldVector<double> *result) const
+ void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
+ std::vector<WorldVector<double> > &result) const
{
- if (grdUhAtQP)
- for (int iq = 0; iq < nPoints; iq++)
- result[iq][xi] += (*factor) * grdUhAtQP[iq][xj];
+ int nPoints = grdUhAtQP.size();
+ for (int iq = 0; iq < nPoints; iq++)
+ result[iq][xi] += (*factor) * grdUhAtQP[iq][xj];
}
private:
@@ -636,9 +635,8 @@ namespace AMDiS {
double factor) const;
/// Implenetation of SecondOrderTerm::weakEval().
- void weakEval(int nPoints,
- const WorldVector<double> *grdUhAtQP,
- WorldVector<double> *result) const;
+ void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
+ std::vector<WorldVector<double> > &result) const;
protected:
/// DOFVector to be evaluated at quadrature points.
@@ -681,9 +679,8 @@ namespace AMDiS {
double factor) const;
/// Implenetation of SecondOrderTerm::weakEval().
- void weakEval(int nPoints,
- const WorldVector<double> *grdUhAtQP,
- WorldVector<double> *result) const;
+ void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
+ std::vector<WorldVector<double> > &result) const;
protected:
/// DOFVector to be evaluated at quadrature points.
@@ -730,9 +727,8 @@ namespace AMDiS {
double factor) const;
/// Implenetation of SecondOrderTerm::weakEval().
- void weakEval(int nPoints,
- const WorldVector<double> *grdUhAtQP,
- WorldVector<double> *result) const;
+ void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
+ std::vector<WorldVector<double> > &result) const;
protected:
@@ -764,6 +760,12 @@ namespace AMDiS {
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
Quadrature *quad = NULL);
+ /// Implementation of \ref OperatorTerm::initElement() for multilpe meshes.
+ void initElement(const ElInfo* smallElInfo,
+ const ElInfo* largeElInfo,
+ SubAssembler* subAssembler,
+ Quadrature *quad = NULL);
+
/// Implements SecondOrderTerm::getLALt().
void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;
@@ -776,9 +778,8 @@ namespace AMDiS {
double factor) const;
/// Implenetation of SecondOrderTerm::weakEval().
- void weakEval(int nPoints,
- const WorldVector<double> *grdUhAtQP,
- WorldVector<double> *result) const;
+ void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
+ std::vector<WorldVector<double> > &result) const;
protected:
DOFVectorBase<double>* vec;
@@ -829,9 +830,8 @@ namespace AMDiS {
double factor) const;
/// Implenetation of SecondOrderTerm::weakEval().
- void weakEval(int nPoints,
- const WorldVector<double> *grdUhAtQP,
- WorldVector<double> *result) const;
+ void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
+ std::vector<WorldVector<double> > &result) const;
protected:
DOFVectorBase<double>* vec;
@@ -878,9 +878,8 @@ namespace AMDiS {
double factor) const;
/// Implenetation of SecondOrderTerm::weakEval().
- void weakEval(int nPoints,
- const WorldVector<double> *grdUhAtQP,
- WorldVector<double> *result) const;
+ void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
+ std::vector<WorldVector<double> > &result) const;
protected:
DOFVectorBase<double>* vec;
@@ -929,9 +928,8 @@ namespace AMDiS {
double factor) const;
/// Implenetation of SecondOrderTerm::weakEval().
- void weakEval(int nPoints,
- const WorldVector<double> *grdUhAtQP,
- WorldVector<double> *result) const;
+ void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
+ std::vector<WorldVector<double> > &result) const;
protected:
DOFVectorBase<double>* vec;
@@ -984,9 +982,8 @@ namespace AMDiS {
double factor) const;
/// Implenetation of SecondOrderTerm::weakEval().
- void weakEval(int nPoints,
- const WorldVector<double> *grdUhAtQP,
- WorldVector<double> *result) const;
+ void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
+ std::vector<WorldVector<double> > &result) const;
protected:
DOFVectorBase<double>* vec;
@@ -1034,9 +1031,8 @@ namespace AMDiS {
double *result,
double factor) const;
- void weakEval(int nPoints,
- const WorldVector<double> *grdUhAtQP,
- WorldVector<double> *result) const;
+ void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
+ std::vector<WorldVector<double> > &result) const;
protected:
/// DOFVector to be evaluated at quadrature points.
@@ -1087,9 +1083,8 @@ namespace AMDiS {
double factor) const;
/// Implenetation of SecondOrderTerm::weakEval().
- void weakEval(int nPoints,
- const WorldVector<double> *grdUhAtQP,
- WorldVector<double> *result) const;
+ void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
+ std::vector<WorldVector<double> > &result) const;
protected:
DOFVectorBase<double>* vec;
@@ -1137,9 +1132,8 @@ namespace AMDiS {
double factor) const;
/// Implenetation of SecondOrderTerm::weakEval().
- void weakEval(int nPoints,
- const WorldVector<double> *grdUhAtQP,
- WorldVector<double> *result) const;
+ void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
+ std::vector<WorldVector<double> > &result) const;
protected:
/// DOFVector to be evaluated at quadrature points.
@@ -1190,9 +1184,8 @@ namespace AMDiS {
double factor) const;
/// Implenetation of SecondOrderTerm::weakEval().
- void weakEval(int nPoints,
- const WorldVector<double> *grdUhAtQP,
- WorldVector<double> *result) const;
+ void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
+ std::vector<WorldVector<double> > &result) const;
protected:
std::vector<DOFVectorBase<double>*> vecs_;
@@ -1248,9 +1241,8 @@ namespace AMDiS {
double factor) const;
/// Implenetation of SecondOrderTerm::weakEval().
- void weakEval(int nPoints,
- const WorldVector<double> *grdUhAtQP,
- WorldVector<double> *result) const;
+ void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
+ std::vector<WorldVector<double> > &result) const;
protected:
std::vector<DOFVectorBase<double>*> vecs_;
@@ -2540,9 +2532,8 @@ namespace AMDiS {
double factor) const;
/// Implements SecondOrderTerm::weakEval().
- void weakEval(int nPoints,
- const WorldVector<double> *grdUhAtQP,
- WorldVector<double> *result) const;
+ void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
+ std::vector<WorldVector<double> > &result) const;
protected:
/// DOFVector to be evaluated at quadrature points.
@@ -2591,9 +2582,8 @@ namespace AMDiS {
double factor) const;
/// Implements SecondOrderTerm::weakEval().
- void weakEval(int nPoints,
- const WorldVector<double> *grdUhAtQP,
- WorldVector<double> *result) const;
+ void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
+ std::vector<WorldVector<double> > &result) const;
protected:
/// DOFVectors to be evaluated at quadrature points.
@@ -2684,9 +2674,8 @@ namespace AMDiS {
double factor) const;
/// Implements SecondOrderTerm::weakEval().
- void weakEval(int nPoints,
- const WorldVector<double> *grdUhAtQP,
- WorldVector<double> *result) const;
+ void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
+ std::vector<WorldVector<double> > &result) const;
private:
/// Stores coordinates at quadrature points. Set in \ref initElement().
@@ -2731,9 +2720,8 @@ namespace AMDiS {
double factor) const;
/// Implements SecondOrderTerm::weakEval().
- void weakEval(int nPoints,
- const WorldVector<double> *grdUhAtQP,
- WorldVector<double> *result) const;
+ void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
+ std::vector<WorldVector<double> > &result) const;
protected:
/// DOFVector to be evaluated at quadrature points.
@@ -3373,17 +3361,15 @@ namespace AMDiS {
}
/// Weak evaluation of all terms in \ref secondOrder.
- void weakEvalSecondOrder(int nPoints,
- const WorldVector<double> *grdUhAtQP,
- WorldVector<double> *result) const
+ void weakEvalSecondOrder(const std::vector<WorldVector<double> > &grdUhAtQP,
+ std::vector<WorldVector<double> > &result) const
{
int myRank = omp_get_thread_num();
- std::vector<OperatorTerm*>::const_iterator termIt;
- for (termIt = secondOrder[myRank].begin();
+ for (std::vector<OperatorTerm*>::const_iterator termIt = secondOrder[myRank].begin();
termIt != secondOrder[myRank].end();
++termIt)
- static_cast<SecondOrderTerm*>(*termIt)->weakEval(nPoints, grdUhAtQP, result);
+ static_cast<SecondOrderTerm*>(*termIt)->weakEval(grdUhAtQP, result);
}
/// Calls getLALt() for each term in \ref secondOrder and adds the results to LALt.
diff --git a/AMDiS/src/ResidualEstimator.cc b/AMDiS/src/ResidualEstimator.cc
index 13ae3f63f22a80eb1097282b5ab88cc90c0f25cc..e81f334d7ad5b6ec5f1b8298b230ac66f44e9480 100644
--- a/AMDiS/src/ResidualEstimator.cc
+++ b/AMDiS/src/ResidualEstimator.cc
@@ -215,12 +215,24 @@ namespace AMDiS {
for (it = dofMat->getOperatorsBegin(), itfac = dofMat->getOperatorEstFactorBegin();
it != dofMat->getOperatorsEnd(); ++it, ++itfac)
- if (*itfac == NULL || **itfac != 0.0)
- (*it)->getAssembler(omp_get_thread_num())->initElement(elInfo, NULL, quad);
+ if (*itfac == NULL || **itfac != 0.0) {
+ 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);
+ }
if (C0 > 0.0)
- for (it = dofVec->getOperatorsBegin(); it != dofVec->getOperatorsEnd(); ++it)
- (*it)->getAssembler(omp_get_thread_num())->initElement(elInfo, NULL, quad);
+ for (it = dofVec->getOperatorsBegin(); it != dofVec->getOperatorsEnd(); ++it) {
+ 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);
+ }
}
@@ -352,6 +364,10 @@ namespace AMDiS {
for (int face = 0; face < nNeighbours; face++) {
Element *neigh = const_cast<Element*>(elInfo->getNeighbour(face));
+ if (dualElInfo) {
+ int smallFace = DualTraverse::getFace(dualElInfo, face);
+ }
+
if (!(neigh && neigh->getMark()))
continue;
@@ -445,7 +461,7 @@ namespace AMDiS {
for (int iq = 0; iq < nPointsSurface; iq++)
localJump[iq].set(0.0);
- (*it)->weakEvalSecondOrder(nPointsSurface, &(grdUhEl[0]), &(localJump[0]));
+ (*it)->weakEvalSecondOrder(grdUhEl, localJump);
double factor = *fac ? **fac : 1.0;
if (factor != 1.0)
diff --git a/AMDiS/src/RobinBC.cc b/AMDiS/src/RobinBC.cc
index def19d38f562cc26d2255f903ed6aedf142024df..9df83ac1838735235480daf5b74b8970979b2f93 100644
--- a/AMDiS/src/RobinBC.cc
+++ b/AMDiS/src/RobinBC.cc
@@ -223,25 +223,19 @@ namespace AMDiS {
getZeroOrderAssembler()->getQuadrature()->getNumPoints() >
(*robinOperators)[0]->getAssembler(omp_get_thread_num())->
getZeroOrderAssembler()->getQuadrature()->getNumPoints())
- {
- neumannQuad = true;
- }
+ neumannQuad = true;
}
}
- Quadrature *quadrature = NULL;
- int face;
-
std::vector<Operator*>::iterator op;
- for (op = matrix->getOperatorsBegin(); op != matrix->getOperatorsEnd(); ++op) {
- (*op)->getAssembler(omp_get_thread_num())->initElement(elInfo);
- }
+ for (op = matrix->getOperatorsBegin(); op != matrix->getOperatorsEnd(); ++op)
+ (*op)->getAssembler(omp_get_thread_num())->initElement(elInfo);
- for (face = 0; face < dim + 1; face++) {
+ for (int face = 0; face < dim + 1; face++) {
elInfo->getNormal(face, normal);
- quadrature =
+ Quadrature *quadrature =
neumannQuad ?
(*neumannOperators)[face]->getAssembler(omp_get_thread_num())->
getZeroOrderAssembler()->getQuadrature() :
@@ -269,8 +263,8 @@ namespace AMDiS {
&f[0],
1.0);
- WorldVector<double> *grdUh = new WorldVector<double>[numPoints];
- WorldVector<double> *A_grdUh = new WorldVector<double>[numPoints];
+ std::vector<WorldVector<double> > grdUh(numPoints);
+ std::vector<WorldVector<double> > A_grdUh(numPoints);
for (iq = 0; iq < numPoints; iq++) {
A_grdUh[iq].set(0.0);
@@ -278,12 +272,8 @@ namespace AMDiS {
basFcts->evalGrdUh(lambda, Lambda, uhEl, &grdUh[iq]);
}
- for (op = matrix->getOperatorsBegin(); op != matrix->getOperatorsEnd(); ++op) {
- (*op)->weakEvalSecondOrder(numPoints,
- grdUh,
- A_grdUh);
- }
-
+ for (op = matrix->getOperatorsBegin(); op != matrix->getOperatorsEnd(); ++op)
+ (*op)->weakEvalSecondOrder(grdUh, A_grdUh);
if (neumannOperators)
(*neumannOperators)[face]->getC(elInfo, numPoints, f);
@@ -292,9 +282,6 @@ namespace AMDiS {
n_A_grdUh = (normal*A_grdUh[iq]) - f[iq];
val += quadrature->getWeight(iq)*sqr(n_A_grdUh);
}
-
- delete [] grdUh;
- delete [] A_grdUh;
}
}