diff --git a/AMDiS/CMakeLists.txt b/AMDiS/CMakeLists.txt
index 356e28fa81bc68b0465dcdd97a72325f0913ee5a..f72db5fdb2dabc9c68fd0c3382c2947422b1bbff 100644
--- a/AMDiS/CMakeLists.txt
+++ b/AMDiS/CMakeLists.txt
@@ -92,7 +92,6 @@ SET(AMDIS_SRC ${SOURCE_DIR}/AdaptBase.cc
${SOURCE_DIR}/Element.cc
${SOURCE_DIR}/ElementData.cc
${SOURCE_DIR}/ElementDofIterator.cc
- ${SOURCE_DIR}/Estimator.cc
${SOURCE_DIR}/FiniteElemSpace.cc
${SOURCE_DIR}/FirstOrderAssembler.cc
${SOURCE_DIR}/FirstOrderTerm.cc
@@ -123,12 +122,10 @@ SET(AMDIS_SRC ${SOURCE_DIR}/AdaptBase.cc
${SOURCE_DIR}/Quadrature.cc
${SOURCE_DIR}/RCNeighbourList.cc
${SOURCE_DIR}/Recovery.cc
- ${SOURCE_DIR}/RecoveryEstimator.cc
${SOURCE_DIR}/RefinementManager.cc
${SOURCE_DIR}/RefinementManager1d.cc
${SOURCE_DIR}/RefinementManager2d.cc
${SOURCE_DIR}/RefinementManager3d.cc
- ${SOURCE_DIR}/ResidualEstimator.cc
${SOURCE_DIR}/RobinBC.cc
${SOURCE_DIR}/ScalableQuadrature.cc
${SOURCE_DIR}/SecondOrderAssembler.cc
@@ -147,6 +144,10 @@ SET(AMDIS_SRC ${SOURCE_DIR}/AdaptBase.cc
${SOURCE_DIR}/VertexVector.cc
${SOURCE_DIR}/ZeroOrderAssembler.cc
${SOURCE_DIR}/ZeroOrderTerm.cc
+ ${SOURCE_DIR}/est/Estimator.cc
+ ${SOURCE_DIR}/est/RecoveryEstimator.cc
+ ${SOURCE_DIR}/est/ResidualEstimator.cc
+ ${SOURCE_DIR}/est/SimpleResidualEstimator.cc
${SOURCE_DIR}/io/ArhReader.cc
${SOURCE_DIR}/io/ArhWriter.cc
${SOURCE_DIR}/io/DataCollector.cc
@@ -345,6 +346,11 @@ INSTALL(FILES ${HEADERS}
DESTINATION include/amdis/nonlin/)
list(APPEND deb_add_dirs "include/amdis/nonlin")
+FILE(GLOB HEADERS "${SOURCE_DIR}/est/*.h")
+INSTALL(FILES ${HEADERS}
+ DESTINATION include/amdis/est/)
+list(APPEND deb_add_dirs "include/amdis/est")
+
FILE(GLOB HEADERS "${SOURCE_DIR}/time/*.h")
INSTALL(FILES ${HEADERS}
DESTINATION include/amdis/time/)
diff --git a/AMDiS/src/AMDiS.h b/AMDiS/src/AMDiS.h
index 16224a2c36d3b1d7ab6b4408d968873629a12956..964027f31bad430f5f2331bc53762428ef72cf55 100644
--- a/AMDiS/src/AMDiS.h
+++ b/AMDiS/src/AMDiS.h
@@ -59,7 +59,6 @@
#include "Element.h"
#include "ElementDofIterator.h"
#include "Error.h"
-#include "Estimator.h"
#include "FiniteElemSpace.h"
#include "FirstOrderTerm.h"
#include "FixVec.h"
@@ -112,6 +111,8 @@
#include "VertexVector.h"
#include "ZeroOrderTerm.h"
+#include "est/Estimator.h"
+
#include "io/ArhReader.h"
#include "io/ArhWriter.h"
#include "io/FileWriter.h"
diff --git a/AMDiS/src/AdaptInstationary.cc b/AMDiS/src/AdaptInstationary.cc
index c402888a53b95f732781f327c1de15973ba85ca6..15839cd14b5184cfac55d779830b25546667d66d 100644
--- a/AMDiS/src/AdaptInstationary.cc
+++ b/AMDiS/src/AdaptInstationary.cc
@@ -12,7 +12,7 @@
#include "AdaptInstationary.h"
#include "Initfile.h"
-#include "Estimator.h"
+#include "est/Estimator.h"
#include "ProblemIterationInterface.h"
#include "ProblemTimeInterface.h"
#include "Serializer.h"
diff --git a/AMDiS/src/AdaptStationary.cc b/AMDiS/src/AdaptStationary.cc
index d6a817c20d3d4715fdb7dc30bcd0dde2c5ffd3cc..378ce85c69d5d852f5cecf5d7e3257ee041538bc 100644
--- a/AMDiS/src/AdaptStationary.cc
+++ b/AMDiS/src/AdaptStationary.cc
@@ -12,7 +12,7 @@
#include "AdaptStationary.h"
#include "Initfile.h"
-#include "Estimator.h"
+#include "est/Estimator.h"
#include "ProblemIterationInterface.h"
#include <math.h>
diff --git a/AMDiS/src/Assembler.cc b/AMDiS/src/Assembler.cc
index 2d5d44499ff7bb76e39976ca5fb2339f7a144a0a..6f2143ad4a0a8d0c8226463c62f42830c001720a 100644
--- a/AMDiS/src/Assembler.cc
+++ b/AMDiS/src/Assembler.cc
@@ -69,7 +69,7 @@ namespace AMDiS {
return;
}
}
-
+
ElementMatrix& mat = rememberElMat ? elementMatrix : userMat;
if (secondOrderAssembler)
@@ -315,8 +315,8 @@ namespace AMDiS {
FUNCNAME("Assembler::matVecAssemble()");
Element *el = elInfo->getElement();
- ElementVector uhOldLoc(operat->uhOld->getFeSpace() == rowFeSpace ? nRow : nCol);
-
+ ElementVector uhOldLoc(operat->uhOld->getFeSpace() == rowFeSpace ?
+ nRow : nCol);
operat->uhOld->getLocalVector(el, uhOldLoc);
if (el != lastMatEl) {
diff --git a/AMDiS/src/CreatorMap.cc b/AMDiS/src/CreatorMap.cc
index 830a1321f343921074fa975e50ee6c530659e206..1a0b10e503ae883f5a81439dac2651570d952840 100644
--- a/AMDiS/src/CreatorMap.cc
+++ b/AMDiS/src/CreatorMap.cc
@@ -17,14 +17,15 @@
#include "ITL_Preconditioner.h"
#include "MatrixVector.h"
#include "SystemVector.h"
-#include "Estimator.h"
-#include "RecoveryEstimator.h"
-#include "ResidualEstimator.h"
+#include "est/Estimator.h"
#include "LeafData.h"
#include "SurfaceRegion_ED.h"
#include "ElementRegion_ED.h"
#include "DOFMatrix.h"
#include "UmfPackSolver.h"
+#include "est/RecoveryEstimator.h"
+#include "est/ResidualEstimator.h"
+#include "est/SimpleResidualEstimator.h"
#include "time/RosenbrockMethod.h"
#include "nonlin/NonLinSolver.h"
@@ -102,6 +103,9 @@ namespace AMDiS {
creator = new ResidualEstimator::Creator;
addCreator("residual", creator);
+ creator = new SimpleResidualEstimator::Creator;
+ addCreator("simple-residual", creator);
+
creator = new RecoveryEstimator::Creator;
addCreator("recovery", creator);
}
diff --git a/AMDiS/src/ElInfo.cc b/AMDiS/src/ElInfo.cc
index 942862dfd9dcb07668dc6175d6327b220f71205b..3b6f9af4cc48d88c85b292a19a3f06d0e2049541 100644
--- a/AMDiS/src/ElInfo.cc
+++ b/AMDiS/src/ElInfo.cc
@@ -56,10 +56,12 @@ namespace AMDiS {
dimOfWorld = Global::getGeo(WORLD);
}
+
ElInfo::~ElInfo()
{}
+
void ElInfo::coordToWorld(const DimVec<double>& l,
WorldVector<double>& w) const
@@ -80,12 +82,14 @@ namespace AMDiS {
}
}
+
double ElInfo::calcDet() const
{
testFlag(Mesh::FILL_COORDS);
return calcDet(coord);
}
+
double ElInfo::calcDet(const FixVec<WorldVector<double>, VERTEX> &coords) const
{
FUNCNAME("ElInfo::calcDet()");
@@ -114,10 +118,8 @@ namespace AMDiS {
det = norm(&e1);
} else {
-
det = (coords[1][0] - coords[0][0]) * (coords[2][1] - coords[0][1]) -
(coords[1][1] - coords[0][1]) * (coords[2][0] - coords[0][0]);
-
}
break;
case 3:
diff --git a/AMDiS/src/ElInfo2d.cc b/AMDiS/src/ElInfo2d.cc
index 56dcc741a3e49f1babeba66fe2938c0cef1cb11d..95985f1dcb80d1d56c47c1ad9f1ceb38e9c4d583 100644
--- a/AMDiS/src/ElInfo2d.cc
+++ b/AMDiS/src/ElInfo2d.cc
@@ -735,13 +735,13 @@ namespace AMDiS {
vectorProduct(elementNormal, e0, normal);
}
- double det = norm(&normal);
+ double detn = norm(&normal);
- TEST_EXIT_DBG(det > 1.e-30)("det = 0 on face %d\n", side);
+ TEST_EXIT_DBG(detn > 1.e-30)("det = 0 on face %d\n", side);
- normal *= 1.0 / det;
+ normal *= 1.0 / detn;
- return det;
+ return detn;
}
@@ -762,13 +762,13 @@ namespace AMDiS {
vectorProduct(e0, e1, elementNormal);
- double det = norm(&elementNormal);
+ double detn = norm(&elementNormal);
- TEST_EXIT_DBG(det > 1.e-30)("det = 0");
+ TEST_EXIT_DBG(detn > 1.e-30)("det = 0");
- elementNormal *= 1.0 / det;
+ elementNormal *= 1.0 / detn;
- return det;
+ return detn;
}
diff --git a/AMDiS/src/ProblemInstat.cc b/AMDiS/src/ProblemInstat.cc
index af79b2f4c474873e0b648a0cadef016ac4e3b0f9..e4f07c84b5c13d836f5b309e39cc11ecf9590c29 100644
--- a/AMDiS/src/ProblemInstat.cc
+++ b/AMDiS/src/ProblemInstat.cc
@@ -16,8 +16,8 @@
#include "ProblemInstat.h"
#include "AdaptStationary.h"
#include "AdaptInstationary.h"
-#include "Estimator.h"
#include "StandardProblemIteration.h"
+#include "est/Estimator.h"
#include "io/FileWriter.h"
namespace AMDiS {
diff --git a/AMDiS/src/ProblemStat.cc b/AMDiS/src/ProblemStat.cc
index de7f0c07029c5d8ff62222910bfdf72d3cda64c3..ce774e3328c179ab143f3c6b77868d4b7c048c0a 100644
--- a/AMDiS/src/ProblemStat.cc
+++ b/AMDiS/src/ProblemStat.cc
@@ -14,14 +14,12 @@
#include <boost/lexical_cast.hpp>
#include "ProblemStat.h"
-#include "RecoveryEstimator.h"
#include "Serializer.h"
#include "AbstractFunction.h"
#include "Operator.h"
#include "SystemVector.h"
#include "DOFMatrix.h"
#include "FiniteElemSpace.h"
-#include "Estimator.h"
#include "Marker.h"
#include "AdaptInfo.h"
#include "io/FileWriter.h"
@@ -35,6 +33,7 @@
#include "PeriodicBC.h"
#include "Lagrange.h"
#include "Flag.h"
+#include "est/Estimator.h"
#include "io/VtkWriter.h"
#include "io/ValueReader.h"
#include "ProblemStatDbg.h"
@@ -790,6 +789,8 @@ namespace AMDiS {
if (matrix)
nnz += matrix->getBaseMatrix().nnz();
}
+
+
// And now assemble boundary conditions on the vectors
assembleBoundaryConditions(rhs->getDOFVector(i),
diff --git a/AMDiS/src/RobinBC.cc b/AMDiS/src/RobinBC.cc
index a5a736c7b3efcdd8f80e93c657e9e8a5c51770d9..4cb69498d6fe0ec9297c4fadb2d16b548a34e469 100644
--- a/AMDiS/src/RobinBC.cc
+++ b/AMDiS/src/RobinBC.cc
@@ -11,11 +11,12 @@
#include "RobinBC.h"
-#include "Estimator.h"
#include "Assembler.h"
#include "DOFVector.h"
#include "DOFMatrix.h"
#include "SurfaceOperator.h"
+#include "est/Estimator.h"
+
#include <math.h>
namespace AMDiS {
diff --git a/AMDiS/src/RobinBC.hh b/AMDiS/src/RobinBC.hh
index f60d96d9d8cb980bed7f8fb72e591ab376f5ebcb..fb4fa85fd54b88c874926677ee532e503e2fce08 100644
--- a/AMDiS/src/RobinBC.hh
+++ b/AMDiS/src/RobinBC.hh
@@ -10,8 +10,8 @@
// See also license.opensource.txt in the distribution.
-#include "Estimator.h"
#include "Assembler.h"
+#include "est/Estimator.h"
template<typename T>
void RobinBC<T>::fillLocalBC(DOFVector<T>* vector,
diff --git a/AMDiS/src/SecondOrderAssembler.cc b/AMDiS/src/SecondOrderAssembler.cc
index 618d5e3f9a3cde3703e48303e631ac0723c5b198..9dd029e05bcd89121c7049f90a58113b6884e271 100644
--- a/AMDiS/src/SecondOrderAssembler.cc
+++ b/AMDiS/src/SecondOrderAssembler.cc
@@ -77,9 +77,9 @@ namespace AMDiS {
newAssembler = new Stand2(op, assembler, quad);
} else {
if (pwConst) {
- newAssembler = new Pre2(op, assembler, quad);
+ newAssembler = new Pre2(op, assembler, quad);
} else {
- newAssembler = new Quad2(op, assembler, quad);
+ newAssembler = new Quad2(op, assembler, quad);
}
}
@@ -330,7 +330,7 @@ namespace AMDiS {
for (int i = 0; i < nRow; i++) {
(*(psi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPsi);
for (int j = 0; j < nCol; j++) {
- tmpVec = (LALt[iq] * grdPhi[j]);
+ tmpVec = (LALt[iq] * grdPhi[j]);
mat[i][j] += quadrature->getWeight(iq) * dot(grdPsi, tmpVec);
}
}
diff --git a/AMDiS/src/SecondOrderTerm.h b/AMDiS/src/SecondOrderTerm.h
index c35455728248b2c80b478a2b43efd7c687dee7dd..202dc781dbeaff921cd819b3348e807d28d075f3 100644
--- a/AMDiS/src/SecondOrderTerm.h
+++ b/AMDiS/src/SecondOrderTerm.h
@@ -132,7 +132,8 @@ namespace AMDiS {
}
/// Implements SecondOrderTerm::getLALt().
- inline void getLALt(const ElInfo *elInfo, vector<mtl::dense2D<double> > &LALt) const
+ inline void getLALt(const ElInfo *elInfo,
+ vector<mtl::dense2D<double> > &LALt) const
{
const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
const int nPoints = static_cast<int>(LALt.size());
diff --git a/AMDiS/src/Triangle.cc b/AMDiS/src/Triangle.cc
index fec3579d7ba9b12dc5b9311301e82fc0fbfbe06e..e7ee0bc22bcb1ad33d79947e74b24c9c627eed0a 100644
--- a/AMDiS/src/Triangle.cc
+++ b/AMDiS/src/Triangle.cc
@@ -52,16 +52,16 @@ namespace AMDiS {
}
const FixVec<int, WORLD>& Triangle::sortFaceIndices(int face,
- FixVec<int,WORLD> *vec) const
+ FixVec<int, WORLD> *vec) const
{
- static MatrixOfFixVecs<FixVec<int,WORLD> > *sorted_2d = NULL;
+ static MatrixOfFixVecs<FixVec<int, WORLD> > *sorted_2d = NULL;
int no = 0;
- FixVec<int,WORLD> *val = NULL;
+ FixVec<int, WORLD> *val = NULL;
const int *vof = vertexOfEdge[face];
if (NULL == sorted_2d) {
- sorted_2d = new MatrixOfFixVecs<FixVec<int,WORLD> >(2, 3, 2, NO_INIT);
+ sorted_2d = new MatrixOfFixVecs<FixVec<int, WORLD> >(2, 3, 2, NO_INIT);
(*sorted_2d)[1][0][1] = (*sorted_2d)[1][1][0] =
(*sorted_2d)[2][0][0] = (*sorted_2d)[2][1][1] = 0;
@@ -83,7 +83,7 @@ namespace AMDiS {
val = &((*sorted_2d)[face][no]);
}
- return *(const_cast<const FixVec<int,WORLD>* >(val));
+ return *(const_cast<const FixVec<int, WORLD>* >(val));
}
diff --git a/AMDiS/src/Triangle.h b/AMDiS/src/Triangle.h
index 4af74f260a22e46b550c76bbffc5bd786410f36e..86e8c3b694e6b72db55da74d3d795c65f66a92eb 100644
--- a/AMDiS/src/Triangle.h
+++ b/AMDiS/src/Triangle.h
@@ -95,7 +95,7 @@ namespace AMDiS {
bool hasSide(Element* sideElem) const;
/// implements Element::sortFaceIndices
- const FixVec<int,WORLD>& sortFaceIndices(int face, FixVec<int,WORLD> *vec) const;
+ const FixVec<int, WORLD>& sortFaceIndices(int face, FixVec<int, WORLD> *vec) const;
/// implements Element::isLine. Returns false because this element is a Triangle
inline bool isLine() const
diff --git a/AMDiS/src/Estimator.cc b/AMDiS/src/est/Estimator.cc
similarity index 100%
rename from AMDiS/src/Estimator.cc
rename to AMDiS/src/est/Estimator.cc
diff --git a/AMDiS/src/Estimator.h b/AMDiS/src/est/Estimator.h
similarity index 100%
rename from AMDiS/src/Estimator.h
rename to AMDiS/src/est/Estimator.h
diff --git a/AMDiS/src/RecoveryEstimator.cc b/AMDiS/src/est/RecoveryEstimator.cc
similarity index 100%
rename from AMDiS/src/RecoveryEstimator.cc
rename to AMDiS/src/est/RecoveryEstimator.cc
diff --git a/AMDiS/src/RecoveryEstimator.h b/AMDiS/src/est/RecoveryEstimator.h
similarity index 100%
rename from AMDiS/src/RecoveryEstimator.h
rename to AMDiS/src/est/RecoveryEstimator.h
diff --git a/AMDiS/src/ResidualEstimator.cc b/AMDiS/src/est/ResidualEstimator.cc
similarity index 98%
rename from AMDiS/src/ResidualEstimator.cc
rename to AMDiS/src/est/ResidualEstimator.cc
index 703b22aef4b9603e8d2a62ad84435903aea22112..c54bae2aee02e4e4dc827c8cd081c60cc0f8d0f3 100644
--- a/AMDiS/src/ResidualEstimator.cc
+++ b/AMDiS/src/est/ResidualEstimator.cc
@@ -322,7 +322,7 @@ namespace AMDiS {
}
if (D2uhqp == NULL && degree > 2 && (*it)->secondOrderTerms()) {
D2uhqp = new WorldMatrix<double>[nPoints];
- uh[system]->getD2AtQPs(elInfo, NULL, quadFast[system], D2uhqp);
+ uh[system]->getD2AtQPs(elInfo, NULL, quadFast[system], D2uhqp);
}
}
}
@@ -533,11 +533,13 @@ namespace AMDiS {
if (factor) {
if (D2UhIq)
- (*it)->evalSecondOrder(nPoints, uhIq, grdUhIq, D2UhIq, result, -factor);
+ (*it)->evalSecondOrder(nPoints, uhIq, grdUhIq, D2UhIq, result, -factor);
if (grdUhIq) {
- (*it)->evalFirstOrderGrdPsi(nPoints, uhIq, grdUhIq, D2UhIq, result, factor);
- (*it)->evalFirstOrderGrdPhi(nPoints, uhIq, grdUhIq, D2UhIq, result, factor);
+ (*it)->evalFirstOrderGrdPsi(nPoints, uhIq, grdUhIq, D2UhIq,
+ result, factor);
+ (*it)->evalFirstOrderGrdPhi(nPoints, uhIq, grdUhIq, D2UhIq,
+ result, factor);
}
if (num_rows(uhIq) > 0)
diff --git a/AMDiS/src/ResidualEstimator.h b/AMDiS/src/est/ResidualEstimator.h
similarity index 97%
rename from AMDiS/src/ResidualEstimator.h
rename to AMDiS/src/est/ResidualEstimator.h
index e0e977c88f45e2148b05919ea3f53b18a47ee1d1..0dd39e1941881a2a82a800b9c99f1c3320bd378f 100644
--- a/AMDiS/src/ResidualEstimator.h
+++ b/AMDiS/src/est/ResidualEstimator.h
@@ -59,7 +59,7 @@ namespace AMDiS {
* \ingroup Estimator
*
* \brief
- * Estimator for scalar problems.
+ * Residual estimator.
*/
class ResidualEstimator : public Estimator
{
@@ -88,8 +88,8 @@ namespace AMDiS {
virtual void init(double timestep);
/** \brief
- * Estimates the error on an element. For more information about the parameter,
- * see the description \ref Estimator::estimateElement.
+ * Estimates the error on an element. For more information about the
+ * parameter, see the description \ref Estimator::estimateElement.
*/
virtual void estimateElement(ElInfo *elInfo, DualElInfo *dualElInfo = NULL);
diff --git a/AMDiS/src/est/SimpleResidualEstimator.cc b/AMDiS/src/est/SimpleResidualEstimator.cc
new file mode 100644
index 0000000000000000000000000000000000000000..588d1aae2b50501f0947d57b1c67ba779a303b7e
--- /dev/null
+++ b/AMDiS/src/est/SimpleResidualEstimator.cc
@@ -0,0 +1,396 @@
+//
+// Software License for AMDiS
+//
+// Copyright (c) 2010 Dresden University of Technology
+// All rights reserved.
+// Authors: Simon Vey, Thomas Witkowski et al.
+//
+// This file is part of AMDiS
+//
+// See also license.opensource.txt in the distribution.
+
+
+#include "est/SimpleResidualEstimator.h"
+#include "Operator.h"
+#include "DOFMatrix.h"
+#include "DOFVector.h"
+#include "Assembler.h"
+#include "Traverse.h"
+#include "Initfile.h"
+
+namespace AMDiS {
+
+ SimpleResidualEstimator::SimpleResidualEstimator(std::string name)
+ : Estimator(name, 0),
+ C0(0.0),
+ C1(0.0)
+ {
+ FUNCNAME("SimpleResidualEstimator::SimpleResidualEstimator()");
+
+ // === Read parameters C0 and C1 from init file. ===
+
+ Parameters::get(name + "->C0", C0);
+ Parameters::get(name + "->C1", C1);
+
+ C0 = C0 > 1.e-25 ? sqr(C0) : 0.0;
+ C1 = C1 > 1.e-25 ? sqr(C1) : 0.0;
+ }
+
+
+ void SimpleResidualEstimator::init(double)
+ {
+ FUNCNAME("SimpleResidualEstimator::init()");
+
+ // === Create data structures. ===
+
+ basFcts = uh[0]->getFeSpace()->getBasisFcts();
+ dim = mesh->getDim();
+ degree = basFcts->getDegree() * 2;
+ quad = Quadrature::provideQuadrature(dim, degree);
+ nPoints = quad->getNumPoints();
+
+ Flag flag = INIT_PHI;
+ if (degree > 2)
+ flag |= INIT_D2_PHI;
+ quadFast = FastQuadrature::provideFastQuadrature(basFcts, *quad, flag);
+
+ uhEl.change_dim(basFcts->getNumber());
+ uhNeigh.change_dim(basFcts->getNumber());
+ riq.change_dim(nPoints);
+ D2uhqp = NULL;
+
+ // === Clear error indicators and mark elements for jumpRes. ===
+
+ TraverseStack stack;
+ ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
+ while (elInfo) {
+ elInfo->getElement()->setEstimation(0.0, 0);
+ elInfo->getElement()->setMark(1);
+ elInfo = stack.traverseNext(elInfo);
+ }
+
+ est_sum = 0.0;
+ est_max = 0.0;
+
+ traverseFlag =
+ Mesh::FILL_NEIGH |
+ Mesh::FILL_COORDS |
+ Mesh::FILL_OPP_COORDS |
+ Mesh::FILL_BOUND |
+ Mesh::FILL_GRD_LAMBDA |
+ Mesh::FILL_DET |
+ Mesh::CALL_LEAF_EL;
+ neighInfo = mesh->createNewElInfo();
+
+
+ // === Prepare date for computing jump residual. ===
+
+ if (C1 > 0.0 && dim > 1) {
+ surfaceQuad = Quadrature::provideQuadrature(dim - 1, degree);
+ nPointsSurface = surfaceQuad->getNumPoints();
+ grdUhEl.resize(nPointsSurface);
+ grdUhNeigh.resize(nPointsSurface);
+ jump.resize(nPointsSurface);
+ localJump.resize(nPointsSurface);
+ nNeighbours = Global::getGeo(NEIGH, dim);
+ lambdaNeigh = new DimVec<WorldVector<double> >(dim, NO_INIT);
+ lambda = new DimVec<double>(dim, NO_INIT);
+ }
+ }
+
+
+ void SimpleResidualEstimator::exit(bool output)
+ {
+ FUNCNAME("SimpleResidualEstimator::exit()");
+
+ // === Calculate the root of the estimations and make output. ===
+
+ est_sum = sqrt(est_sum);
+
+ if (output)
+ MSG("estimate = %.8e\n", est_sum);
+
+ /// === Delete data structures. ===
+
+ if (D2uhqp != NULL)
+ delete [] D2uhqp;
+
+ if (C1 && (dim > 1)) {
+ delete lambdaNeigh;
+ delete lambda;
+ }
+
+ delete neighInfo;
+ }
+
+
+ void SimpleResidualEstimator::estimateElement(ElInfo *elInfo, DualElInfo *)
+ {
+ FUNCNAME("SimpleResidualEstimator::estimateElement()");
+
+ TEST_EXIT(matrix[0])("Should not happen!\n");
+
+ // Get pointer to element object.
+ Element *el = elInfo->getElement();
+ // Get error estimation of this element.
+ double est_el = el->getEstimation(0);
+ // Shortcut for the system matrix.
+ DOFMatrix *dofMat = const_cast<DOFMatrix*>(matrix[0]);
+ // Shortcut for the right hand side DOF vector
+ DOFVector<double> *dofVec = const_cast<DOFVector<double>*>(fh[0]);
+
+ // === Init assembler ===
+
+ std::vector<Operator*>::iterator it;
+ std::vector<double*>::iterator itfac;
+ // Matrix assembler are only initialized with the corresponding multiply
+ // factors are != 0
+ for (it = dofMat->getOperatorsBegin(), itfac = dofMat->getOperatorEstFactorBegin();
+ it != dofMat->getOperatorsEnd(); ++it, ++itfac)
+ if (*itfac == NULL || **itfac != 0.0)
+ (*it)->getAssembler()->initElement(elInfo, NULL, quad);
+
+
+ // Vector assembler are only initialized if C0 is set. Note that the jump
+ // residual (thus C1) does not contain the right hand side.
+ if (C0 > 0.0)
+ for (it = dofVec->getOperatorsBegin(); it != dofVec->getOperatorsEnd(); ++it)
+ (*it)->getAssembler()->initElement(elInfo, NULL, quad);
+
+
+ // === Compute element residuals and time error estimation. ===
+ if (C0 > 0.0)
+ est_el += computeElementResidual(elInfo);
+
+ // === Compute jump residuals. ===
+ if (C1 > 0.0 && dim > 1)
+ est_el += computeJumpResidual(elInfo);
+
+ // === Update global residual variables. ===
+ el->setEstimation(est_el, 0);
+ el->setMark(0);
+ est_sum += est_el;
+ est_max = std::max(est_max, est_el);
+ }
+
+
+ double SimpleResidualEstimator::computeElementResidual(ElInfo *elInfo)
+ {
+ FUNCNAME("SimpleResidualEstimator::computeElementResidual()");
+
+ double det = elInfo->getDet();
+ double h2 = h2_from_det(det, dim);
+ riq = 0.0;
+ DOFMatrix *dofMat = const_cast<DOFMatrix*>(matrix[0]);
+ DOFVector<double> *dofVec = const_cast<DOFVector<double>*>(fh[0]);
+
+ // === If there is a valid left hand side operator get the solution ===
+ // === vector or its derivations on the quadrature points of the ===
+ // === element. ===
+ std::vector<Operator*>::iterator it;
+ std::vector<double*>::iterator itfac;
+ for (it = dofMat->getOperatorsBegin(), itfac = dofMat->getOperatorEstFactorBegin();
+ it != dofMat->getOperatorsEnd(); ++it, ++itfac) {
+ if (*itfac == NULL || **itfac != 0.0) {
+ if (num_rows(uhQP) == 0 && (*it)->zeroOrderTerms()) {
+ uhQP.change_dim(nPoints);
+ uh[0]->getVecAtQPs(elInfo, NULL, quadFast, uhQP);
+ }
+
+ if (D2uhqp == NULL && degree > 2 && (*it)->secondOrderTerms()) {
+ D2uhqp = new WorldMatrix<double>[nPoints];
+ uh[0]->getD2AtQPs(elInfo, NULL, quadFast, D2uhqp);
+ }
+ }
+ }
+
+ // === Compute the element residual and store it in irq. ===
+ r(elInfo, nPoints, uhQP, D2uhqp, dofMat, dofVec, quad, riq);
+
+
+ // === Add integral over r square. ===
+ double result = 0.0;
+ for (int iq = 0; iq < nPoints; iq++)
+ result += quad->getWeight(iq) * riq[iq] * riq[iq];
+
+ if (norm == NO_NORM || norm == L2_NORM)
+ result = C0 * h2 * h2 * det * result;
+ else
+ result = C0 * h2 * det * result;
+
+ return result;
+ }
+
+
+ double SimpleResidualEstimator::computeJumpResidual(ElInfo *elInfo)
+ {
+ FUNCNAME("SimpleResidualEstimator::computeJumpResidual()");
+
+ // === Init temporary variables. ===
+ double result = 0.0;
+ int dow = Global::getGeo(WORLD);
+ Element *el = elInfo->getElement();
+ const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
+ double det = elInfo->getDet();
+ double h2 = h2_from_det(det, dim);
+
+ /// === Compute jump on all faces of the current element. ===
+ for (int face = 0; face < nNeighbours; face++) {
+ // Pointer to neighbouring element.
+ Element *neigh = const_cast<Element*>(elInfo->getNeighbour(face));
+
+ // If there is no neighbour, or we have already visited this pair,
+ // then continue with next one.
+ if (!(neigh && neigh->getMark()))
+ continue;
+
+ // === The next lines are only to set all information about the ===
+ // === neighbouring element. ===
+ int oppV = elInfo->getOppVertex(face);
+ el->sortFaceIndices(face, &faceIndEl);
+ neigh->sortFaceIndices(oppV, &faceIndNeigh);
+ neighInfo->setElement(const_cast<Element*>(neigh));
+ neighInfo->setFillFlag(Mesh::FILL_COORDS);
+
+ for (int i = 0; i < dow; i++)
+ neighInfo->getCoord(oppV)[i] = elInfo->getOppCoord(face)[i];
+
+ for (int i = 0; i < dim; i++) {
+ int i1 = faceIndEl[i];
+ int i2 = faceIndNeigh[i];
+ for (int j = 0; j < dow; j++)
+ neighInfo->getCoord(i2)[j] = elInfo->getCoord(i1)[j];
+ }
+
+ // Compute determinant of the neighbouring element.
+ double detNeigh = abs(neighInfo->calcGrdLambda(*lambdaNeigh));
+
+ // Reset data.
+ for (int iq = 0; iq < nPointsSurface; iq++)
+ jump[iq].set(0.0);
+
+ // Get solution vector on the nodes of the current element and
+ // its neighbour.
+ uh[0]->getLocalVector(el, uhEl);
+ uh[0]->getLocalVector(neigh, uhNeigh);
+
+ // Compute the jump of the gradients of the solution over the face.
+ for (int iq = 0; iq < nPointsSurface; iq++) {
+ (*lambda)[face] = 0.0;
+ for (int i = 0; i < dim; i++)
+ (*lambda)[faceIndEl[i]] = surfaceQuad->getLambda(iq, i);
+
+ basFcts->evalGrdUh(*lambda, grdLambda, uhEl, &grdUhEl[iq]);
+
+ (*lambda)[oppV] = 0.0;
+ for (int i = 0; i < dim; i++)
+ (*lambda)[faceIndNeigh[i]] = surfaceQuad->getLambda(iq, i);
+
+ basFcts->evalGrdUh(*lambda, *lambdaNeigh, uhNeigh, &grdUhNeigh[iq]);
+
+ grdUhEl[iq] -= grdUhNeigh[iq];
+ } // for iq
+
+
+ // === Compute the action of the second order operators on the jump. ===
+ std::vector<double*>::iterator fac;
+ std::vector<Operator*>::iterator it;
+ DOFMatrix *mat = const_cast<DOFMatrix*>(matrix[0]);
+ for (it = mat->getOperatorsBegin(), fac = mat->getOperatorEstFactorBegin();
+ it != mat->getOperatorsEnd(); ++it, ++fac) {
+
+ if (*fac == NULL || **fac != 0.0) {
+ for (int iq = 0; iq < nPointsSurface; iq++)
+ localJump[iq].set(0.0);
+
+ (*it)->weakEvalSecondOrder(grdUhEl, localJump);
+
+ double factor = *fac ? **fac : 1.0;
+ if (factor != 1.0)
+ for (int i = 0; i < nPointsSurface; i++)
+ localJump[i] *= factor;
+
+ for (int i = 0; i < nPointsSurface; i++)
+ jump[i] += localJump[i];
+ }
+ } // for (it = ...
+
+
+ // === Compute the squared integral of the jump and multiply with ===
+ // === the appropriate constants. ===
+
+ double val = 0.0;
+ for (int iq = 0; iq < nPointsSurface; iq++)
+ val += surfaceQuad->getWeight(iq) * (jump[iq] * jump[iq]);
+
+ double d = 0.5 * (det + detNeigh);
+
+ if (norm == NO_NORM || norm == L2_NORM)
+ val *= C1 * h2_from_det(d, dim) * d;
+ else
+ val *= C1 * d;
+
+ neigh->setEstimation(neigh->getEstimation(0) + val, 0);
+ result += val;
+ } // for face
+
+
+ double val = fh[0]->getBoundaryManager()->boundResidual(elInfo, matrix[0], uh[0]);
+ if (norm == NO_NORM || norm == L2_NORM)
+ val *= C1 * h2;
+ else
+ val *= C1;
+ result += val;
+
+ return result;
+ }
+
+
+ void SimpleResidualEstimator::r(const ElInfo *elInfo,
+ int nPoints,
+ const ElementVector& uhIq,
+ const WorldMatrix<double> *D2UhIq,
+ DOFMatrix *A,
+ DOFVector<double> *fh,
+ Quadrature *quad,
+ ElementVector& result)
+ {
+ FUNCNAME("SimpleResidualEstimator::r()");
+
+ std::vector<Operator*>::iterator it;
+ std::vector<double*>::iterator fac;
+
+ // === Get left hand side on quadrature points. ===
+ for (it = A->getOperatorsBegin(), fac = A->getOperatorEstFactorBegin();
+ it != A->getOperatorsEnd(); ++it, ++fac) {
+
+ double factor = *fac ? **fac : 1.0;
+
+ if (factor) {
+ if (D2UhIq)
+ (*it)->evalSecondOrder(nPoints, uhIq, NULL, D2UhIq, result, -factor);
+
+ if (num_rows(uhIq) > 0)
+ (*it)->evalZeroOrder(nPoints, uhIq, NULL, D2UhIq, result, factor);
+ }
+ }
+
+
+ // === Get right hand side on quadrature points. ===
+ for (it = fh->getOperatorsBegin(), fac = fh->getOperatorEstFactorBegin();
+ it != fh->getOperatorsEnd(); ++it, ++fac) {
+
+ double factor = *fac ? **fac : 1.0;
+
+ if (factor) {
+ ElementVector fx(nPoints, 0.0);
+ (*it)->getC(elInfo, nPoints, fx);
+
+ for (int iq = 0; iq < nPoints; iq++)
+ result[iq] -= factor * fx[iq];
+ }
+ }
+ }
+
+
+}
diff --git a/AMDiS/src/est/SimpleResidualEstimator.h b/AMDiS/src/est/SimpleResidualEstimator.h
new file mode 100644
index 0000000000000000000000000000000000000000..cb0695672946d2079682d678887ff2aefa62bd14
--- /dev/null
+++ b/AMDiS/src/est/SimpleResidualEstimator.h
@@ -0,0 +1,197 @@
+// ============================================================================
+// == ==
+// == AMDiS - Adaptive multidimensional simulations ==
+// == ==
+// == http://www.amdis-fem.org ==
+// == ==
+// ============================================================================
+//
+// Software License for AMDiS
+//
+// Copyright (c) 2010 Dresden University of Technology
+// All rights reserved.
+// Authors: Simon Vey, Thomas Witkowski et al.
+//
+// This file is part of AMDiS
+//
+// See also license.opensource.txt in the distribution.
+
+
+
+/** \file SimpleResidualEstimator.h */
+
+/** \defgroup Estimator Estimator module
+ * @{ <img src="estimator.png"> @}
+ */
+
+#ifndef AMDIS_SIMPLE_RESIDUAL_ESTIMATOR_H
+#define AMDIS_SIMPLE_RESIDUAL_ESTIMATOR_H
+
+#include "Estimator.h"
+#include "FixVec.h"
+
+namespace AMDiS {
+
+ /**
+ * \ingroup Estimator
+ *
+ * \brief
+ * Simple Residual estimator.
+ *
+ * The restrictions for using the simple residual estimator are the following:
+ * - no system of PDEs, thus a PDE with only one unknown
+ * - no first order terms
+ * - no time dependent terms
+ */
+ class SimpleResidualEstimator : public Estimator
+ {
+ public:
+ /// Creator class used in the OEMSolverMap.
+ class Creator : public EstimatorCreator
+ {
+ public:
+ Creator()
+ : EstimatorCreator()
+ {}
+
+ virtual ~Creator()
+ {}
+
+ /// Returns a new ODirSolver object.
+ Estimator* create()
+ {
+ return new SimpleResidualEstimator(name);
+ }
+ };
+
+ /// Constructor.
+ SimpleResidualEstimator(std::string name);
+
+ /// Init the error estimator.
+ void init(double);
+
+ /// Run estimator on one element.
+ /// \param[in] elInfo Info object for the element to be estimated.
+ /// \param[in] dualElInfo Not used here. In general, this may be used for
+ /// estimating with the multi-mesh technique.
+ void estimateElement(ElInfo *elInfo, DualElInfo *dualElInfo = NULL);
+
+ /// Finalize the error estimator, i.e., delete all temporary data structures.
+ void exit(bool output = true);
+
+ /// Returns pow(det,2.0/dim).
+ inline double h2_from_det(double det, int dim)
+ {
+ return pow(det, 2.0 / dim);
+ }
+
+ protected:
+ /// Computes the element residual for a given element.
+ double computeElementResidual(ElInfo *elInfo);
+
+ /// Computes the jump residual for a given element.
+ double computeJumpResidual(ElInfo *elInfo);
+
+ /** \brief
+ * Returns residual square at quadrature point. Not Member of
+ * Estimator to avoid multiple instantiation.
+ *
+ * \param[in] elInfo Current element information object.
+ * \param[in] nPoints Number of quadrature points on this element.
+ * \param[in] uhIq Solution vector on element quadrature points.
+ * \param[in] D2UhIq Second derivations of the solution vector on
+ * element quadrature points.
+ * \param[in] A Matrix that contains the left hand side operators.
+ * \param[in] fh Vector that contains the right hand side operators.
+ * \param[in] quad Object for numerical quadrature.
+ * \param[out] result Vector containing the residual on the quadrature
+ * points.
+ */
+ void r(const ElInfo *elInfo,
+ int nPoints,
+ const ElementVector& uhIq,
+ const WorldMatrix<double> *D2UhIq,
+ DOFMatrix *A,
+ DOFVector<double> *fh,
+ Quadrature *quad,
+ ElementVector& result);
+
+ protected:
+ /// Constant in front of element residual.
+ double C0;
+
+ /// Constant in front of edge/face residual.
+ double C1;
+
+ /// Number of quadrature points.
+ int nPoints;
+
+ /// Dimension of the mesh.
+ int dim;
+
+ /// Polynomial degree.
+ int degree;
+
+ /// Object for numerical quadrature
+ Quadrature *quad;
+
+ /// Object for fast numerical quadrature
+ FastQuadrature *quadFast;
+
+ /// Pointer to the basis functions of the FE space.
+ const BasisFunction *basFcts;
+
+ /// Vector that stores all global DOFs of one element.
+ ElementVector uhEl;
+
+ /// Vector that stores all global DOFs of one element.
+ ElementVector uhNeigh;
+
+ /// Vector that stores values on all quadrature points (QP) on one element.
+ ElementVector uhQP;
+
+ /// Matrix that stores the second derivations on all quadrature points
+ /// on one element.
+ WorldMatrix<double> *D2uhqp;
+
+ /// Stores the element residual computed at the quadrature points of
+ /// the element.
+ ElementVector riq;
+
+ /// Pointer to the information object of some neighbouring element.
+ ElInfo *neighInfo;
+
+ /// Surface quadrature object, used for computing the jump residual.
+ Quadrature *surfaceQuad;
+
+ /// Number of quadrature points of the surface quadrature object.
+ int nPointsSurface;
+
+ /// Stores on all surface quadrature points the gradient of a function.
+ std::vector<WorldVector<double> > grdUhEl;
+
+ /// Stores on all surface quadrature points the gradient of a function.
+ std::vector<WorldVector<double> > grdUhNeigh;
+
+ /// Stores on all surface quadrature points the jump of a function.
+ std::vector<WorldVector<double> > jump;
+
+ /// Stores on all surface quadrature points the jump of a function.
+ std::vector<WorldVector<double> > localJump;
+
+ /// Vector to store, for a given element, all its neighbouring element indices.
+ WorldVector<int> faceIndEl;
+
+ /// Vector to store, for a given element, all its neighbouring element indices.
+ WorldVector<int> faceIndNeigh;
+
+ DimVec<WorldVector<double> > *lambdaNeigh;
+
+ DimVec<double> *lambda;
+
+ /// Maximal number of neighbours an element may have in the used dimension.
+ int nNeighbours;
+ };
+}
+
+#endif // AMDIS_SIMPLE_RESIDUAL_ESTIMATOR_H
diff --git a/AMDiS/src/io/MacroReader.cc b/AMDiS/src/io/MacroReader.cc
index 80343b35c1cb711ba8ed770b1ebf3c3edd924f9d..6846853ac310d1b0b10471f872836e6a75bc6bb4 100644
--- a/AMDiS/src/io/MacroReader.cc
+++ b/AMDiS/src/io/MacroReader.cc
@@ -1270,6 +1270,7 @@ namespace AMDiS {
}
}
+
int MacroReader::basicCheckFct(ElInfo* elinfo, Mesh *mesh)
{
FUNCNAME("MacroReader::basicCheckFct()");