diff --git a/AMDiS/src/Estimator.cc b/AMDiS/src/Estimator.cc
index 0354e8add638dd3213b48beebd760ea4a062316b..a6c2542b36a5c4f273a4353ae3baf99dd8df5582 100755
--- a/AMDiS/src/Estimator.cc
+++ b/AMDiS/src/Estimator.cc
@@ -100,6 +100,7 @@ namespace AMDiS {
TraverseStack stack;
ElInfo *elInfo = stack.traverseFirst(mesh, -1, traverseFlag);
+
while(elInfo) {
estimateElement(elInfo);
elInfo = stack.traverseNext(elInfo);
diff --git a/AMDiS/src/Operator.hh b/AMDiS/src/Operator.hh
new file mode 100644
index 0000000000000000000000000000000000000000..7e9c5b1143e1cda9f7b1baefcf87b86d0f0949b3
--- /dev/null
+++ b/AMDiS/src/Operator.hh
@@ -0,0 +1,49 @@
+namespace AMDiS {
+
+ template <typename T>
+ void Operator::addSecondOrderTerm(T *term)
+ {
+ secondOrder[0].push_back(term);
+ term->operat = this;
+
+ for (int i = 1; i < omp_get_max_threads(); i++) {
+ T *newTerm = NEW T(static_cast<const T>(*term));
+ secondOrder[i].push_back(newTerm);
+ }
+ }
+
+ template <typename T>
+ void Operator::addFirstOrderTerm(T *term,
+ FirstOrderType type)
+ {
+ if (type == GRD_PSI) {
+ firstOrderGrdPsi[0].push_back(term);
+ } else {
+ firstOrderGrdPhi[0].push_back(term);
+ }
+ term->operat = this;
+
+ for (int i = 1; i < omp_get_max_threads(); i++) {
+ T *newTerm = NEW T(static_cast<const T>(*term));
+ if (type == GRD_PSI) {
+ firstOrderGrdPsi[i].push_back(newTerm);
+ } else {
+ firstOrderGrdPhi[i].push_back(newTerm);
+ }
+ }
+ }
+
+ template <typename T>
+ void Operator::addZeroOrderTerm(T *term)
+ {
+ zeroOrder[0].push_back(term);
+ term->operat = this;
+
+ for (int i = 1; i < omp_get_max_threads(); i++) {
+ T *newTerm = NEW T(static_cast<const T>(*term));
+ zeroOrder[i].push_back(newTerm);
+ }
+ }
+
+
+}
diff --git a/AMDiS/src/PardisoSolver.cc b/AMDiS/src/PardisoSolver.cc
index f53ef71cb8e7c5eb34a326ca1fc64ec6f7756e51..907294ef3b2d72489fe97643c5081d5cce4b94f7 100644
--- a/AMDiS/src/PardisoSolver.cc
+++ b/AMDiS/src/PardisoSolver.cc
@@ -125,8 +125,8 @@ namespace AMDiS {
iparm[7] = 2; // Max numbers of iterative refinement steps
iparm[9] = 13; // Perturb the pivot elements with 1e-13
iparm[10] = 1; // Use nonsymmetric permutation and scaling MPS
- iparm[17] = -1; // Output: Number of nonzeros in the factor LU
- iparm[18] = -1; // Output: Mflops for LU factorization
+ iparm[17] = 0; // Output: Number of nonzeros in the factor LU
+ iparm[18] = 0; // Output: Mflops for LU factorization
// Maximum number of numerical factorizations
MKL_INT maxfct = 1;
@@ -135,7 +135,7 @@ namespace AMDiS {
MKL_INT mnum = 1;
// Print statistical information in file
- MKL_INT msglvl = 1;
+ MKL_INT msglvl = 0;
// Error flag
MKL_INT error = 0;
diff --git a/AMDiS/src/ResidualEstimator.cc b/AMDiS/src/ResidualEstimator.cc
new file mode 100644
index 0000000000000000000000000000000000000000..202762640d663a7e1cb92a976bdc32e080145f06
--- /dev/null
+++ b/AMDiS/src/ResidualEstimator.cc
@@ -0,0 +1,543 @@
+#include "ResidualEstimator.h"
+#include "Operator.h"
+#include "DOFMatrix.h"
+#include "DOFVector.h"
+#include "Assembler.h"
+#include "Traverse.h"
+#include "Parameters.h"
+
+namespace AMDiS {
+
+ ResidualEstimator::ResidualEstimator(::std::string name, int r)
+ : Estimator(name, r),
+ C0(1.0),
+ C1(1.0),
+ C2(1.0),
+ C3(1.0)
+ {
+ GET_PARAMETER(0, name + "->C0", "%f", &C0);
+ GET_PARAMETER(0, name + "->C1", "%f", &C1);
+ GET_PARAMETER(0, name + "->C2", "%f", &C2);
+ GET_PARAMETER(0, name + "->C3", "%f", &C3);
+
+ C0 = C0 > 1.e-25 ? sqr(C0) : 0.0;
+ C1 = C1 > 1.e-25 ? sqr(C1) : 0.0;
+ C2 = C2 > 1.e-25 ? sqr(C2) : 0.0;
+ C3 = C3 > 1.e-25 ? sqr(C3) : 0.0;
+ }
+
+ void ResidualEstimator::init(double ts)
+ {
+ FUNCNAME("ResidualEstimator::init()");
+
+ timestep = ts;
+
+ mesh = uh[row == -1 ? 0 : row]->getFESpace()->getMesh();
+
+ numSystems = static_cast<int>(uh.size());
+ TEST_EXIT_DBG(numSystems > 0)("no system set\n");
+
+ dim = mesh->getDim();
+ basFcts = GET_MEMORY(const BasisFunction*, numSystems);
+ quadFast = GET_MEMORY(FastQuadrature*, numSystems);
+
+ degree = 0;
+ for (int system = 0; system < numSystems; system++) {
+ basFcts[system] = uh[system]->getFESpace()->getBasisFcts();
+ degree = ::std::max(degree, basFcts[system]->getDegree());
+ }
+
+ degree *= 2;
+
+ quad = Quadrature::provideQuadrature(dim, degree);
+ numPoints = quad->getNumPoints();
+
+ Flag flag = INIT_PHI | INIT_GRD_PHI;
+ if (degree > 2) {
+ flag |= INIT_D2_PHI;
+ }
+
+ for (int system = 0; system < numSystems; system++) {
+ quadFast[system] = FastQuadrature::provideFastQuadrature(basFcts[system],
+ *quad,
+ flag);
+ }
+
+ uhEl = GET_MEMORY(double*, numSystems);
+ uhOldEl = timestep ? GET_MEMORY(double*, numSystems) : NULL;
+
+ for (int system = 0; system < numSystems; system++) {
+ uhEl[system] = GET_MEMORY(double, basFcts[system]->getNumber());
+ if (timestep)
+ uhOldEl[system] = GET_MEMORY(double, basFcts[system]->getNumber());
+ }
+
+ uhQP = timestep ? GET_MEMORY(double, numPoints) : NULL;
+ uhOldQP = timestep ? GET_MEMORY(double, numPoints) : NULL;
+
+ riq = GET_MEMORY(double, numPoints);
+
+ TraverseStack stack;
+ ElInfo *elInfo = NULL;
+
+ // clear error indicators and mark elements for jumpRes
+ elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
+ while (elInfo) {
+ elInfo->getElement()->setEstimation(0.0, row);
+ elInfo->getElement()->setMark(1);
+ elInfo = stack.traverseNext(elInfo);
+ }
+
+ est_sum = 0.0;
+ est_max = 0.0;
+ est_t_sum = 0.0;
+ est_t_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;
+ }
+
+ void ResidualEstimator::exit(bool output)
+ {
+ FUNCNAME("ResidualEstimator::exit()");
+
+ est_sum = sqrt(est_sum);
+ est_t_sum = sqrt(est_t_sum);
+
+ for (int system = 0; system < numSystems; system++) {
+ FREE_MEMORY(uhEl[system], double, basFcts[system]->getNumber());
+ if (timestep)
+ FREE_MEMORY(uhOldEl[system], double, basFcts[system]->getNumber());
+ }
+
+ FREE_MEMORY(uhEl, double*, numSystems);
+ if (timestep)
+ FREE_MEMORY(uhOldEl, double*, numSystems);
+
+ if (timestep) {
+ FREE_MEMORY(uhQP, double, numPoints);
+ FREE_MEMORY(uhOldQP, double, numPoints);
+ }
+
+ if (output) {
+ MSG("estimate = %.8e\n", est_sum);
+ if (C3) {
+ MSG("time estimate = %.8e\n", est_t_sum);
+ }
+ }
+
+ FREE_MEMORY(riq, double, numPoints);
+ FREE_MEMORY(basFcts, const BasisFunction*, numSystems);
+ FREE_MEMORY(quadFast, FastQuadrature*, numSystems);
+ }
+
+ void ResidualEstimator::estimateElement(ElInfo *elInfo)
+ {
+ FUNCNAME("ResidualEstimator::estimateElement()");
+
+ double est_el = 0.0;
+ double val = 0.0;
+ Element *el, *neigh;
+ int iq;
+
+ TEST_EXIT_DBG(numSystems > 0)("no system set\n");
+
+ ::std::vector<Operator*>::iterator it;
+
+ WorldVector<double> *grdUh_qp = NULL;
+ WorldMatrix<double> *D2uhqp = NULL;
+
+ el = elInfo->getElement();
+
+ double det = elInfo->getDet();
+ const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
+
+ est_el = el->getEstimation(row);
+
+ double h2 = h2_from_det(det, dim);
+
+ for (iq = 0; iq < numPoints; iq++) {
+ riq[iq] = 0.0;
+ }
+
+ for (int system = 0; system < numSystems; system++) {
+
+ if (matrix[system] == NULL)
+ continue;
+
+ // init assemblers
+ ::std::vector<Operator*>::iterator it;
+
+ for (it = const_cast<DOFMatrix*>(matrix[system])->getOperatorsBegin();
+ it != const_cast<DOFMatrix*>(matrix[system])->getOperatorsEnd();
+ ++it) {
+ (*it)->getAssembler(omp_get_thread_num())->initElement(elInfo, quad);
+ }
+
+ for (it = const_cast<DOFVector<double>*>(fh[system])->getOperatorsBegin();
+ it != const_cast<DOFVector<double>*>(fh[system])->getOperatorsEnd();
+ ++it) {
+ (*it)->getAssembler(omp_get_thread_num())->initElement(elInfo, quad);
+ }
+
+ if (timestep) {
+ TEST_EXIT_DBG(uhOld[system])("no uhOld\n");
+ uhOld[system]->getLocalVector(el, uhOldEl[system]);
+
+ // ===== time and element residuals
+ if (C0 || C3) {
+ uh[system]->getVecAtQPs(elInfo, NULL, quadFast[system], uhQP);
+ uhOld[system]->getVecAtQPs(elInfo, NULL, quadFast[system], uhOldQP);
+
+ if (C3 && uhOldQP && system == ::std::max(row, 0)) {
+ for (val = iq = 0; iq < numPoints; iq++) {
+ double tiq = (uhQP[iq] - uhOldQP[iq]);
+ val += quad->getWeight(iq) * tiq * tiq;
+ }
+ double v = C3 * det * val;
+ est_t_sum += v;
+ est_t_max = max(est_t_max, v);
+ }
+ }
+ }
+
+
+ if (C0) {
+ for (it = const_cast<DOFMatrix*>(matrix[system])->getOperatorsBegin();
+ it != const_cast<DOFMatrix*>(matrix[system])->getOperatorsEnd();
+ ++it) {
+ if ((*it)->zeroOrderTerms() && !uhQP) {
+ uhQP = GET_MEMORY(double, numPoints);
+ uh[system]->getVecAtQPs(elInfo, NULL, quadFast[system], uhQP);
+ }
+ if ((*it)->firstOrderTermsGrdPsi() || (*it)->firstOrderTermsGrdPhi()
+ && !grdUh_qp) {
+ grdUh_qp = new WorldVector<double>[numPoints];
+ uh[system]->getGrdAtQPs(elInfo, NULL, quadFast[system], grdUh_qp);
+ }
+ if ((*it)->secondOrderTerms() && !D2uhqp) {
+ if (degree > 2) {
+ D2uhqp = new WorldMatrix<double>[numPoints];
+ uh[system]->getD2AtQPs(elInfo, NULL, quadFast[system], D2uhqp);
+ }
+ }
+ }
+
+ r(elInfo,
+ numPoints,
+ uhQP,
+ grdUh_qp,
+ D2uhqp,
+ uhOldQP,
+ NULL, // grdUhOldQP
+ NULL, // D2UhOldQP
+ matrix[system],
+ fh[system],
+ quad,
+ riq);
+ }
+ }
+
+ // add integral over r square
+ for (val = iq = 0; iq < numPoints; iq++) {
+ val += quad->getWeight(iq) * riq[iq] * riq[iq];
+ }
+
+ if (timestep != 0.0 || norm == NO_NORM || norm == L2_NORM)
+ val = C0 * h2 * h2 * det * val;
+ else
+ val = C0 * h2 * det * val;
+
+ est_el += val;
+
+ // ===== jump residuals
+ if (C1 && (dim > 1)) {
+ int neighbours = Global::getGeo(NEIGH, dim);
+ int face;
+ Quadrature *surfaceQuad = Quadrature::provideQuadrature(dim - 1, degree);
+ int numPointsSurface = surfaceQuad->getNumPoints();
+ Vector<WorldVector<double> > jump(numPointsSurface);
+
+
+ for (face = 0; face < neighbours; face++) {
+ neigh = const_cast<Element*>(elInfo->getNeighbour(face));
+ if (neigh && neigh->getMark()) {
+ ElInfo *neighInfo = mesh->createNewElInfo();
+ WorldVector<int> faceIndEl, faceIndNeigh;
+ int oppV = elInfo->getOppVertex(face);
+ DimVec<WorldVector<double> > LambdaNeigh(dim, NO_INIT);
+ double detNeigh;
+ DimVec<double> lambda(dim, NO_INIT);
+
+ el->sortFaceIndices(face, &faceIndEl);
+ neigh->sortFaceIndices(oppV, &faceIndNeigh);
+
+ neighInfo->setElement(const_cast<Element*>(neigh));
+ neighInfo->setFillFlag(Mesh::FILL_COORDS);
+
+ int dow = Global::getGeo(WORLD);
+
+ int j, i1, i2;
+
+ for (int i = 0; i < dow; i++)
+ neighInfo->getCoord(oppV)[i] = elInfo->getOppCoord(face)[i];
+
+ // periodic leaf data ?
+ ElementData *ldp = el->getElementData()->getElementData(PERIODIC);
+
+ bool periodicCoords = false;
+
+ if (ldp) {
+ ::std::list<LeafDataPeriodic::PeriodicInfo>::iterator it;
+
+ for (it = dynamic_cast<LeafDataPeriodic*>(ldp)->getInfoList().begin();
+ it != dynamic_cast<LeafDataPeriodic*>(ldp)->getInfoList().end();
+ ++it) {
+ if (it->elementSide == face) {
+ for (int i = 0; i < dim; i++) {
+ i1 = faceIndEl[i];
+ i2 = faceIndNeigh[i];
+
+ for (j = 0; j < dim; j++) {
+ if (i1 == el->getVertexOfPosition(INDEX_OF_DIM(dim - 1,
+ dim),
+ face,
+ j)) {
+ break;
+ }
+ }
+
+ TEST_EXIT_DBG(j != dim)("vertex i1 not on face ???\n");
+
+ neighInfo->getCoord(i2) = (*(it->periodicCoords))[j];
+ }
+ periodicCoords = true;
+ break;
+ }
+ }
+ }
+
+ if (!periodicCoords) {
+ for (int i = 0; i < dim; i++) {
+ i1 = faceIndEl[i];
+ i2 = faceIndNeigh[i];
+ for (j = 0; j < dow; j++)
+ neighInfo->getCoord(i2)[j] = elInfo->getCoord(i1)[j];
+ }
+ }
+
+ Parametric *parametric = mesh->getParametric();
+ if (parametric) {
+ neighInfo = parametric->addParametricInfo(neighInfo);
+ }
+
+ detNeigh = abs(neighInfo->calcGrdLambda(LambdaNeigh));
+
+ Vector<WorldVector<double> > grdUhEl(numPointsSurface);
+ Vector<WorldVector<double> > grdUhNeigh(numPointsSurface);
+
+ for (iq = 0; iq < numPointsSurface; iq++) {
+ jump[iq].set(0.0);
+ }
+
+
+ for (int system = 0; system < numSystems; system++) {
+ if (matrix[system] == NULL)
+ continue;
+
+ uh[system]->getLocalVector(el, uhEl[system]);
+
+ double* uhNeigh;
+ uhNeigh = new double[basFcts[system]->getNumber()];
+ uh[system]->getLocalVector(neigh, uhNeigh);
+
+
+ for (iq = 0; iq < numPointsSurface; iq++) {
+ lambda[face] = 0.0;
+ for (int i = 0; i < dim; i++) {
+ lambda[faceIndEl[i]] = surfaceQuad->getLambda(iq, i);
+ }
+
+ basFcts[system]->evalGrdUh(lambda,
+ Lambda,
+ uhEl[system],
+ &grdUhEl[iq]);
+
+ lambda[oppV] = 0.0;
+ for (int i = 0; i < dim; i++) {
+ lambda[faceIndNeigh[i]] = surfaceQuad->getLambda(iq, i);
+ }
+
+ basFcts[system]->evalGrdUh(lambda,
+ LambdaNeigh,
+ uhNeigh,
+ &grdUhNeigh[iq]);
+
+ grdUhEl[iq] -= grdUhNeigh[iq];
+ }
+
+ delete [] uhNeigh;
+
+ ::std::vector<double*>::iterator fac;
+
+ for (it = const_cast<DOFMatrix*>(matrix[system])->getOperatorsBegin(),
+ fac = const_cast<DOFMatrix*>(matrix[system])->getOperatorEstFactorBegin();
+ it != const_cast<DOFMatrix*>(matrix[system])->getOperatorsEnd();
+ ++it, ++fac) {
+ Vector<WorldVector<double> > localJump(numPointsSurface);
+ for (iq = 0; iq < numPointsSurface; iq++) {
+ localJump[iq].set(0.0);
+ }
+
+ (*it)->weakEvalSecondOrder(numPointsSurface,
+ grdUhEl.getValArray(),
+ localJump.getValArray());
+ double factor = *fac ? **fac : 1.0;
+ if (factor != 1.0) {
+ for (int i = 0; i < numPointsSurface; i++) {
+ localJump[i] *= factor;
+ }
+ }
+
+ for (int i = 0; i < numPointsSurface; i++) {
+ jump[i] += localJump[i];
+ }
+ }
+ }
+
+
+ for (val = iq = 0; iq < numPointsSurface; 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;
+
+ if (parametric) {
+ neighInfo = parametric->removeParametricInfo(neighInfo);
+ }
+
+ DELETE neighInfo;
+
+ neigh->setEstimation(neigh->getEstimation(row) + val, row);
+
+ est_el += val;
+
+ }
+ }
+
+ val = fh[::std::max(row, 0)]->
+ getBoundaryManager()->
+ boundResidual(elInfo, matrix[::std::max(row, 0)], uh[::std::max(row, 0)]);
+ if (norm == NO_NORM || norm == L2_NORM)
+ val *= C1 * h2;
+ else
+ val *= C1;
+
+ est_el += val;
+ }
+
+
+ el->setEstimation(est_el, row);
+
+ est_sum += est_el;
+ est_max = max(est_max, est_el);
+
+ elInfo->getElement()->setMark(0);
+ }
+
+ void r(const ElInfo *elInfo,
+ int numPoints,
+ const double *uhIq,
+ const WorldVector<double> *grdUhIq,
+ const WorldMatrix<double> *D2UhIq,
+ const double *uhOldIq,
+ const WorldVector<double> *grdUhOldIq,
+ const WorldMatrix<double> *D2UhOldIq,
+ DOFMatrix *A,
+ DOFVector<double> *fh,
+ Quadrature *quad,
+ double *result)
+ {
+ ::std::vector<Operator*>::iterator it;
+ ::std::vector<double*>::iterator fac;
+ double factor;
+
+ // lhs
+ for (it = const_cast<DOFMatrix*>(A)->getOperatorsBegin(),
+ fac = const_cast<DOFMatrix*>(A)->getOperatorEstFactorBegin();
+ it != const_cast<DOFMatrix*>(A)->getOperatorsEnd();
+ ++it, ++fac) {
+
+ factor = *fac ? **fac : 1.0;
+ if (factor) {
+ if (D2UhIq) {
+ (*it)->evalSecondOrder(numPoints, uhIq, grdUhIq, D2UhIq, result, -factor);
+ }
+
+ if (grdUhIq) {
+ (*it)->evalFirstOrderGrdPsi(numPoints, uhIq, grdUhIq, D2UhIq, result, factor);
+ (*it)->evalFirstOrderGrdPhi(numPoints, uhIq, grdUhIq, D2UhIq, result, factor);
+ }
+
+ if (uhIq) {
+ (*it)->evalZeroOrder(numPoints, uhIq, grdUhIq, D2UhIq, result, factor);
+ }
+ }
+ }
+
+ // rhs
+ for (it = const_cast<DOFVector<double>*>(fh)->getOperatorsBegin(),
+ fac = const_cast<DOFVector<double>*>(fh)->getOperatorEstFactorBegin();
+ it != const_cast<DOFVector<double>*>(fh)->getOperatorsEnd();
+ ++it, ++fac) {
+
+ factor = *fac ? **fac : 1.0;
+ if (factor) {
+ if ((*it)->getUhOld()) {
+ if (D2UhOldIq) {
+ (*it)->evalSecondOrder(numPoints,
+ uhOldIq, grdUhOldIq, D2UhOldIq,
+ result, factor);
+ }
+ if (grdUhOldIq) {
+ (*it)->evalFirstOrderGrdPsi(numPoints,
+ uhOldIq, grdUhOldIq, D2UhOldIq,
+ result, -factor);
+ (*it)->evalFirstOrderGrdPhi(numPoints,
+ uhOldIq, grdUhOldIq, D2UhOldIq,
+ result, -factor);
+ }
+ if (uhOldIq) {
+ (*it)->evalZeroOrder(numPoints,
+ uhOldIq, grdUhOldIq, D2UhOldIq,
+ result, -factor);
+ }
+ } else {
+ double *fx = GET_MEMORY(double, numPoints);
+ for (int iq = 0; iq < numPoints; iq++) {
+ fx[iq] = 0.0;
+ }
+ (*it)->getC(elInfo, numPoints, fx);
+
+ for (int iq = 0; iq < numPoints; iq++) {
+ result[iq] -= factor * fx[iq];
+ }
+ FREE_MEMORY(fx, double, numPoints);
+ }
+ }
+ }
+ }
+
+
+}
diff --git a/AMDiS/src/ResidualEstimator.h b/AMDiS/src/ResidualEstimator.h
new file mode 100644
index 0000000000000000000000000000000000000000..d07db3880a593b96634719a514efa50fa75c121a
--- /dev/null
+++ b/AMDiS/src/ResidualEstimator.h
@@ -0,0 +1,154 @@
+// ============================================================================
+// == ==
+// == AMDiS - Adaptive multidimensional simulations ==
+// == ==
+// ============================================================================
+// == ==
+// == crystal growth group ==
+// == ==
+// == Stiftung caesar ==
+// == Ludwig-Erhard-Allee 2 ==
+// == 53175 Bonn ==
+// == germany ==
+// == ==
+// ============================================================================
+// == ==
+// == http://www.caesar.de/cg/AMDiS ==
+// == ==
+// ============================================================================
+
+/** \file ResidualEstimator.h */
+
+/** \defgroup Estimator Estimator module
+ * @{ <img src="estimator.png"> @}
+ */
+
+#ifndef AMDIS_RESIDUALESTIMATOR_H
+#define AMDIS_RESIDUALESTIMATOR_H
+
+#include "Estimator.h"
+#include "FixVec.h"
+
+namespace AMDiS {
+
+ /** \brief
+ * Returns residual square at quadrature point. Not Member of
+ * Estimator to avoid multiple instantiation.
+ */
+ void r(const ElInfo *elInfo,
+ int numPoints,
+ const double *uhIq,
+ const WorldVector<double> *grdUhIq,
+ const WorldMatrix<double> *D2UhIq,
+ const double *uhOldIq,
+ const WorldVector<double> *grdUhOldIq,
+ const WorldMatrix<double> *D2UhOldIq,
+ DOFMatrix *A,
+ DOFVector<double> *fh,
+ Quadrature *quad,
+ double *result);
+
+ /** \brief
+ * Returns pow(det,2.0/dim). Not Member of
+ * Estimator to avoid multiple instantiation.
+ */
+ inline double h2_from_det(double det, int dim) {
+ return pow(det, 2.0 / dim);
+ };
+
+ // ============================================================================
+ // ===== class ResidualEstimator ==============================================
+ // ============================================================================
+
+ /**
+ * \ingroup Estimator
+ *
+ * \brief
+ * Estimator for scalar problems.
+ */
+ class ResidualEstimator : public Estimator
+ {
+ public:
+ MEMORY_MANAGED(ResidualEstimator);
+
+ /** \brief
+ * Creator class used in the OEMSolverMap.
+ */
+ class Creator : public EstimatorCreator
+ {
+ public:
+ MEMORY_MANAGED(Creator);
+
+ Creator() : EstimatorCreator() {};
+
+ virtual ~Creator() {};
+
+ /** \brief
+ * Returns a new ODirSolver object.
+ */
+ Estimator* create() {
+ return NEW ResidualEstimator(name, row);
+ };
+ };
+
+
+ /** \brief
+ * Constructor.
+ */
+ ResidualEstimator(::std::string name, int r);
+
+ virtual void init(double timestep);
+
+ virtual void estimateElement(ElInfo *elInfo);
+
+ virtual void exit(bool output = true);
+
+ protected:
+ /** \brief
+ * Constant in front of element residual
+ */
+ double C0;
+
+ /** \brief
+ * Constant in front of edge/face residual
+ */
+ double C1;
+
+ /** \brief
+ * Constant in front of coarsening term
+ */
+ double C2;
+
+ /** \brief
+ * Constant in front of the time
+ */
+ double C3;
+
+
+ int numSystems;
+
+ int numPoints;
+
+ int dim;
+
+ int degree;
+
+ Quadrature *quad;
+
+ FastQuadrature **quadFast;
+
+ const BasisFunction **basFcts;
+
+ double **uhEl;
+
+ double **uhOldEl;
+
+ double *uhQP;
+
+ double *uhOldQP;
+
+ double *riq;
+ };
+}
+
+#endif // AMDIS_RESIDUALESTIMATOR_H
diff --git a/AMDiS/src/ResidualParallelEstimator.cc b/AMDiS/src/ResidualParallelEstimator.cc
new file mode 100644
index 0000000000000000000000000000000000000000..809608916a5680cb8690468575ef680b210f62f6
--- /dev/null
+++ b/AMDiS/src/ResidualParallelEstimator.cc
@@ -0,0 +1,108 @@
+#include "ResidualParallelEstimator.h"
+#include "ResidualEstimator.h"
+#include "TraverseParallel.h"
+#include "Parameters.h"
+#include "OpenMP.h"
+
+namespace AMDiS {
+ ResidualParallelEstimator::ResidualParallelEstimator(::std::string name, int r)
+ : Estimator(name, r),
+ C3(1.0)
+ {
+ GET_PARAMETER(0, name + "->C3", "%f", &C3);
+ C3 = C3 > 1.e-25 ? sqr(C3) : 0.0;
+
+ seqEstimators_.resize(omp_get_max_threads());
+
+ for (int i = 0; i < omp_get_max_threads(); i++) {
+ seqEstimators_[i] = NEW ResidualEstimator(name, r);
+ }
+ }
+
+ ResidualParallelEstimator::~ResidualParallelEstimator()
+ {
+ for (int i = 0; i < static_cast<int>(seqEstimators_.size()); i++) {
+ DELETE seqEstimators_[i];
+ }
+ }
+
+ void ResidualParallelEstimator::addSystem(DOFMatrix *matrix_,
+ DOFVector<double> *uh_,
+ DOFVector<double> *fh_,
+ DOFVector<double> *uhOld_)
+ {
+ Estimator::addSystem(matrix_, uh_, fh_, uhOld_);
+
+ for (int i = 0; i < static_cast<int>(seqEstimators_.size()); i++) {
+ seqEstimators_[i]->addSystem(matrix_, uh_, fh_, uhOld_);
+ }
+ }
+
+ void ResidualParallelEstimator::addUhOldToSystem(int system,
+ DOFVector<double> *uhOld_)
+ {
+ Estimator::addUhOldToSystem(system, uhOld_);
+
+ for (int i = 0; i < static_cast<int>(seqEstimators_.size()); i++) {
+ seqEstimators_[i]->addUhOldToSystem(system, uhOld_);
+ }
+ }
+
+ double ResidualParallelEstimator::estimate(double ts)
+ {
+ mesh = uh[row == -1 ? 0 : row]->getFESpace()->getMesh();
+
+ for (int i = 0; i < static_cast<int>(seqEstimators_.size()); i++) {
+ seqEstimators_[i]->init(ts);
+ }
+
+ TraverseParallelStack stack;
+
+#pragma omp parallel
+ {
+ ResidualEstimator *estimator = seqEstimators_[omp_get_thread_num()];
+
+ 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;
+
+ ElInfo *elInfo = stack.traverseFirst(mesh, -1, traverseFlag);
+
+ while (elInfo) {
+ estimator->estimateElement(elInfo);
+ elInfo = stack.traverseNext(elInfo);
+ }
+ }
+
+ est_sum = 0.0;
+ est_max = 0.0;
+ est_t_sum = 0.0;
+ est_t_max = 0.0;
+
+ for (int i = 0; i < static_cast<int>(seqEstimators_.size()); i++) {
+ seqEstimators_[i]->exit(false);
+ est_sum += pow(seqEstimators_[i]->getErrorSum(), 2);
+ est_max = max(est_max, seqEstimators_[i]->getErrorMax());
+ est_t_sum += pow(seqEstimators_[i]->getTimeEst(), 2);
+ est_t_max = max(est_t_max, seqEstimators_[i]->getTimeEstMax());
+ }
+
+ est_sum = sqrt(est_sum);
+ est_t_sum = sqrt(est_t_sum);
+
+ MSG("estimate = %.8e\n", est_sum);
+ if (C3) {
+ MSG("time estimate = %.8e\n", est_t_sum);
+ }
+
+ return est_sum;
+ }
+
+
+}
+
diff --git a/AMDiS/src/ResidualParallelEstimator.h b/AMDiS/src/ResidualParallelEstimator.h
new file mode 100644
index 0000000000000000000000000000000000000000..d68173ba8dbeb65195023b5c9fd45b22c2822012
--- /dev/null
+++ b/AMDiS/src/ResidualParallelEstimator.h
@@ -0,0 +1,109 @@
+// ============================================================================
+// == ==
+// == AMDiS - Adaptive multidimensional simulations ==
+// == ==
+// ============================================================================
+// == ==
+// == crystal growth group ==
+// == ==
+// == Stiftung caesar ==
+// == Ludwig-Erhard-Allee 2 ==
+// == 53175 Bonn ==
+// == germany ==
+// == ==
+// ============================================================================
+// == ==
+// == http://www.caesar.de/cg/AMDiS ==
+// == ==
+// ============================================================================
+
+/** \file ResidualParallelEstimator.h */
+
+/** \defgroup Estimator Estimator module
+ * @{ <img src="estimator.png"> @}
+ */
+
+#ifndef AMDIS_RESIDUALPARALLELESTIMATOR_H
+#define AMDIS_RESIDUALPARALLELESTIMATOR_H
+
+#include "ResidualEstimator.h"
+#include <vector>
+
+namespace AMDiS {
+
+ // ============================================================================
+ // ===== class ResidualParallelEstimator ======================================
+ // ============================================================================
+
+ /**
+ * \ingroup Estimator
+ *
+ * \brief
+ * Parallel Residual Estimator for scalar problems.
+ */
+ class ResidualParallelEstimator : public Estimator
+ {
+ public:
+ MEMORY_MANAGED(ResidualParallelEstimator);
+
+ /** \brief
+ * Creator class used in the OEMSolverMap.
+ */
+ class Creator : public EstimatorCreator
+ {
+ public:
+ MEMORY_MANAGED(Creator);
+
+ Creator() : EstimatorCreator() {};
+
+ virtual ~Creator() {};
+
+ /** \brief
+ * Returns a new ODirSolver object.
+ */
+ Estimator* create() {
+ return NEW ResidualParallelEstimator(name, row);
+ };
+ };
+
+ /** \brief
+ * Constructor.
+ */
+ ResidualParallelEstimator(::std::string name, int r);
+
+ /** \brief
+ * Destructor.
+ */
+ ~ResidualParallelEstimator();
+
+ /** \brief
+ *
+ */
+ virtual void addSystem(DOFMatrix *matrix_,
+ DOFVector<double> *uh_,
+ DOFVector<double> *fh_,
+ DOFVector<double> *uhOld_ = NULL);
+
+ /** \brief
+ *
+ */
+ virtual void addUhOldToSystem(int system, DOFVector<double> *uhOld_);
+
+ /** \brief
+ *
+ */
+ virtual double estimate(double timestep);
+
+ private:
+ ::std::vector<ResidualEstimator*> seqEstimators_;
+
+
+ /** \brief
+ * Constant in fromt of the time
+ */
+ double C3;
+ };
+
+}
+
+#endif // AMDIS_RESIDUALPARALLELESTIMATOR_H
diff --git a/AMDiS/src/Traverse.cc b/AMDiS/src/Traverse.cc
index 51786ff9e5fff95022d9ea2674069a1b11b87969..a5c7e3505caada3f7f00aa35dbcafe486875d1cc 100644
--- a/AMDiS/src/Traverse.cc
+++ b/AMDiS/src/Traverse.cc
@@ -340,7 +340,7 @@ namespace AMDiS {
el = elinfo_stack[stack_used]->getElement();
if ((el == NULL) || (el->getFirstChild() == NULL)) {
- return(elinfo_stack[stack_used]);
+ return (elinfo_stack[stack_used]);
}
} else {
el = elinfo_stack[stack_used]->getElement();
diff --git a/AMDiS/src/TraverseParallel.cc b/AMDiS/src/TraverseParallel.cc
new file mode 100644
index 0000000000000000000000000000000000000000..d6b99c7a4b85913d969efa1362950cda6903546d
--- /dev/null
+++ b/AMDiS/src/TraverseParallel.cc
@@ -0,0 +1,49 @@
+#include "TraverseParallel.h"
+#include "OpenMP.h"
+#include "Mesh.h"
+
+#ifdef _OPENMP
+
+namespace AMDiS {
+
+ TraverseParallelStack::TraverseParallelStack(int nThreads)
+ {
+ if (nThreads == 0) {
+ nThreads_ = omp_get_max_threads();
+ } else {
+ nThreads_ = nThreads;
+ }
+
+ stacks_.resize(nThreads_);
+ int i = 0;
+ for (::std::vector<TraverseStack*>::iterator it = stacks_.begin();
+ it != stacks_.end();
+ ++it) {
+
+ (*it) = NEW TraverseStack();
+
+ (*it)->setMyThreadId(i++);
+ (*it)->setMaxThreads(nThreads_);
+ }
+ }
+
+ TraverseParallelStack::~TraverseParallelStack()
+ {
+ for (::std::vector<TraverseStack*>::iterator it = stacks_.begin();
+ it != stacks_.end();
+ ++it) {
+
+ DELETE (*it);
+ }
+ }
+
+ ElInfo* TraverseParallelStack::traverseFirst(Mesh *mesh, int level, Flag fill_flag)
+ {
+ TEST_EXIT(fill_flag.isSet(Mesh::CALL_LEAF_EL))("not yet implemented");
+ TEST_EXIT(stacks_[omp_get_thread_num()])("something wrong with parallel stack traverse");
+
+ return stacks_[omp_get_thread_num()]->traverseFirst(mesh, level, fill_flag);
+ }
+}
+
+#endif // _OPENMP
diff --git a/AMDiS/src/TraverseParallel.h b/AMDiS/src/TraverseParallel.h
new file mode 100644
index 0000000000000000000000000000000000000000..83f5072b1945228375a4431978b7c7c8f98d4c3e
--- /dev/null
+++ b/AMDiS/src/TraverseParallel.h
@@ -0,0 +1,69 @@
+// ============================================================================
+// == ==
+// == AMDiS - Adaptive multidimensional simulations ==
+// == ==
+// ============================================================================
+// == ==
+// == crystal growth group ==
+// == ==
+// == Stiftung caesar ==
+// == Ludwig-Erhard-Allee 2 ==
+// == 53175 Bonn ==
+// == germany ==
+// == ==
+// ============================================================================
+// == ==
+// == http://www.caesar.de/cg/AMDiS ==
+// == ==
+// ============================================================================
+
+/** \file TraverseParallel.h */
+
+#ifndef AMDIS_TRAVERSEPARALLEL_H
+#define AMDIS_TRAVERSEPARALLEL_H
+
+#ifdef _OPENMP
+
+#include "OpenMP.h"
+#include "Traverse.h"
+#include <vector>
+
+namespace AMDiS {
+
+ // ===========================================================================
+ // ===== class TraverseParallelStack =========================================
+ // ===========================================================================
+
+ /** \ingroup Traverse
+ * \brief
+ */
+ class TraverseParallelStack
+ {
+ public:
+ MEMORY_MANAGED(TraverseParallelStack);
+
+ TraverseParallelStack(int nThreads = 0);
+
+ ~TraverseParallelStack();
+
+ ElInfo* traverseFirst(Mesh *mesh, int level, Flag fill_flag);
+
+ inline ElInfo* traverseNext(ElInfo* elinfo_old) {
+ return stacks_[omp_get_thread_num()]->traverseNext(elinfo_old);
+ }
+
+ private:
+ /** \brief
+ * Number of threads using the stack in parallel.
+ */
+ int nThreads_;
+
+ ::std::vector<TraverseStack*> stacks_;
+ };
+
+}
+
+#endif // _OPENMP
+
+#endif // AMDIS_TRAVERSEPARALLEL_H
+