diff --git a/AMDiS/src/DOFVector.h b/AMDiS/src/DOFVector.h
index 7e91601236c775ce5f5796c1d0d9fc9a08575668..893ded7ebaddb3407577b047dbace66899417b74 100644
--- a/AMDiS/src/DOFVector.h
+++ b/AMDiS/src/DOFVector.h
@@ -471,7 +471,23 @@ namespace AMDiS {
}
/// Calculates Integral of this DOFVector
- double Int(Quadrature* q = NULL) const;
+ double Int(Quadrature* q = NULL) const
+ {
+ return Int(-1, q);
+ }
+
+ /** \brief
+ * Calculates Integral of this DOFVector.
+ *
+ * \param[in] meshlevel Then mesh level on which the integral should be
+ * calculated. Usually, only -1 for the leaf level
+ * makes sence, but other values can be used, e.g.,
+ * for debugging.
+ * \param[in] q Quadrature object. If not specified, the function
+ * creates a new quadrature object.
+ */
+ double Int(int meshLevel, Quadrature* q = NULL) const;
+
/** \brief
* Calculates Integral of this DOFVector over parts of the domain
diff --git a/AMDiS/src/DOFVector.hh b/AMDiS/src/DOFVector.hh
index e749f3075b86b2360639db8367a79eee92b245b7..30b165f72b9116d55bce55af9ab700d29b98512e 100644
--- a/AMDiS/src/DOFVector.hh
+++ b/AMDiS/src/DOFVector.hh
@@ -493,10 +493,10 @@ namespace AMDiS {
template<typename T>
- double DOFVector<T>::Int(Quadrature* q) const
+ double DOFVector<T>::Int(int meshLevel, Quadrature* q) const
{
FUNCNAME("DOFVector::Int()");
-
+
Mesh* mesh = this->feSpace->getMesh();
if (!q) {
@@ -507,15 +507,20 @@ namespace AMDiS {
FastQuadrature *quadFast =
FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), *q, INIT_PHI);
- T result; nullify(result);
+ Flag flag = Mesh::FILL_COORDS | Mesh::FILL_DET;
+ if (meshLevel == -1)
+ flag |= Mesh::CALL_LEAF_EL;
+ else
+ flag |= Mesh::CALL_EL_LEVEL;
+
+ double result = 0.0;
int nPoints = quadFast->getNumPoints();
mtl::dense_vector<T> uh_vec(nPoints);
TraverseStack stack;
- ElInfo *elInfo = stack.traverseFirst(mesh, -1,
- Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET);
+ ElInfo *elInfo = stack.traverseFirst(mesh, meshLevel, flag);
while (elInfo) {
double det = elInfo->getDet();
- T normT; nullify(normT);
+ double normT = 0.0;
this->getVecAtQPs(elInfo, NULL, quadFast, uh_vec);
for (int iq = 0; iq < nPoints; iq++)
normT += quadFast->getWeight(iq) * (uh_vec[iq]);
diff --git a/AMDiS/src/Error.h b/AMDiS/src/Error.h
index 61af95643e547b09fec66b9ab405005abea5384f..bb9e32ab592c4340ebc50ef9ec5e116d24d361f4 100644
--- a/AMDiS/src/Error.h
+++ b/AMDiS/src/Error.h
@@ -69,6 +69,14 @@ namespace AMDiS {
double* max,
bool writeLeafData = false,
int comp = 0);
+
+ static double L2Err_ElementWise(const AbstractFunction<T, WorldVector<double> >& u,
+ const DOFVector<T>& uh,
+ double* max,
+ bool writeLeafData,
+ int comp,
+ int level,
+ map<int, double>& estMap);
public:
static T errUFct(const DimVec<double>& lambda);
diff --git a/AMDiS/src/Error.hh b/AMDiS/src/Error.hh
index 860c99acf79a9a201cefad1fa99ca0bae033a7a6..74c5b2e54df43fb92ff045b1405773104c973c82 100644
--- a/AMDiS/src/Error.hh
+++ b/AMDiS/src/Error.hh
@@ -296,4 +296,92 @@ namespace AMDiS {
return (sqrt(l2Err2));
}
+ template<typename T>
+ double Error<T>::L2Err_ElementWise(const AbstractFunction<T, WorldVector<double> >& u,
+ const DOFVector<T>& uh,
+ double* max,
+ bool writeLeafData,
+ int comp,
+ int level,
+ map<int, double>& estMap)
+ {
+ FUNCNAME("Error<T>::L2Err_ElementWise()");
+
+ const FiniteElemSpace *fe_space;
+ Quadrature *q = NULL;
+ writeInLeafData = writeLeafData;
+ component = comp;
+
+ if (!(pU = &u)) {
+ ERROR("no function u specified; doing nothing\n");
+ return(0.0);
+ }
+
+ if (!(errUh = &uh) || !(fe_space = uh.getFeSpace())) {
+ ERROR("no discrete function or no fe_space for it; doing nothing\n");
+ return(0.0);
+ }
+
+ if (!(basFct = fe_space->getBasisFcts())) {
+ ERROR("no basis functions at discrete solution ; doing nothing\n");
+ return(0.0);
+ }
+
+ int dim = fe_space->getMesh()->getDim();
+ int deg = u.getDegree();
+ int degree = deg ? deg : 2 * fe_space->getBasisFcts()->getDegree() - 2;
+
+ q = Quadrature::provideQuadrature(dim, degree);
+ quadFast = FastQuadrature::provideFastQuadrature(basFct, *q, INIT_PHI);
+
+ double maxErr = 0.0, l2Err2 = 0.0, l2Norm2 = 0.0;
+ double *u_vec = new double[quadFast->getQuadrature()->getNumPoints()];
+ int nPoints = quadFast->getNumPoints();
+ mtl::dense_vector<double> uh_vec(nPoints);
+
+ Flag flag = Mesh::FILL_COORDS | Mesh::FILL_DET | Mesh::FILL_GRD_LAMBDA;
+ if (level == -1)
+ flag |= Mesh::CALL_LEAF_EL;
+ else
+ flag |= Mesh::CALL_EL_LEVEL;
+
+ TraverseStack stack;
+ ElInfo *elInfo = stack.traverseFirst(fe_space->getMesh(), level, flag);
+
+ while (elInfo) {
+ elinfo = elInfo;
+
+ quadFast->getQuadrature()->fAtQp(errU, u_vec);
+ errUh->getVecAtQPs(elInfo, NULL, quadFast, uh_vec);
+ double det = elInfo->getDet();
+ double l2_err_el = 0.0, err = 0.0;
+
+ double int_u_vec = 0.0;
+ double int_uh_vec = 0.0;
+
+ for (int i = 0; i < nPoints; i++) {
+ int_u_vec += quadFast->getWeight(i) * u_vec[i];
+ int_uh_vec += quadFast->getWeight(i) * uh_vec[i];
+ }
+
+ l2_err_el = pow(fabs(det * int_u_vec - det * int_uh_vec), 2.0);
+ l2Err2 += l2_err_el;
+ maxErr = std::max(maxErr, l2_err_el);
+
+// if (writeInLeafData)
+// elInfo->getElement()->setEstimation(l2_err_el, component);
+
+ estMap[elInfo->getElement()->getIndex()] = l2_err_el;
+
+ elInfo = stack.traverseNext(elInfo);
+ }
+
+ delete [] u_vec;
+
+ if (max)
+ *max = maxErr;
+
+ return (sqrt(l2Err2));
+ }
+
}
diff --git a/AMDiS/src/nonlin/NonLinSolver.h b/AMDiS/src/nonlin/NonLinSolver.h
index bb1d8b306d020d6aded233f37e22dcd357a12603..0857d870f112194af68a3470e4a994d41740a910 100644
--- a/AMDiS/src/nonlin/NonLinSolver.h
+++ b/AMDiS/src/nonlin/NonLinSolver.h
@@ -86,7 +86,6 @@ namespace AMDiS {
init();
solution.set(0.0);
int result = nlsolve(mat, solution, rhs, adaptInfo, prob);
- exit();
return result;
}
diff --git a/AMDiS/src/nonlin/ProblemNonLin.cc b/AMDiS/src/nonlin/ProblemNonLin.cc
index 655e76a048fbb07549be817ca24c725c0497e446..e5c5f3125cf48763d31920c80eabc7b0f48b3505 100644
--- a/AMDiS/src/nonlin/ProblemNonLin.cc
+++ b/AMDiS/src/nonlin/ProblemNonLin.cc
@@ -60,49 +60,13 @@ namespace AMDiS {
}
- void ProblemNonLin::solve(AdaptInfo *adaptInfo, bool)
+ void ProblemNonLin::solve(AdaptInfo *adaptInfo, bool b0, bool b1)
{
TEST_EXIT(nonLinSolver)("no non-linear solver!\n");
- nonLinSolver->solve(solverMatrix, *solution, *rhs, adaptInfo, this);
- }
+ MSG("HERE A\n");
-
- void ProblemNonLin::buildAfterCoarsen(AdaptInfo *adaptInfo, Flag f, bool b0, bool b1)
- {
- FUNCNAME("ProblemNonLin::buildAfterCoarsen()");
-
- ProblemStat::buildAfterCoarsen(adaptInfo, f, b0, b1);
-
-#if 0
- int nMeshes = static_cast<int>(meshes.size());
-
- for (int i = 0; i < nMeshes; i++)
- meshes[i]->dofCompress();
-
- for (int i = 0; i < nComponents; i++)
- MSG("%d DOFs for %s\n",
- componentSpaces[i]->getAdmin()->getUsedSize(),
- componentSpaces[i]->getName().c_str());
-
- // for all elements ...
- for (int i = 0; i < nComponents; i++) {
- TraverseStack stack;
- ElInfo *elInfo = stack.traverseFirst(componentMeshes[i], -1,
- Mesh::CALL_LEAF_EL |
- Mesh::FILL_BOUND |
- Mesh::FILL_COORDS |
- Mesh::FILL_DET |
- Mesh::FILL_GRD_LAMBDA);
- while (elInfo) {
- if (solution->getDOFVector(i)->getBoundaryManager())
- solution->getDOFVector(i)->
- getBoundaryManager()->fillBoundaryConditions(elInfo, solution->getDOFVector(i));
-
- elInfo = stack.traverseNext(elInfo);
- }
- }
-#endif
+ nonLinSolver->solve(solverMatrix, *solution, *rhs, adaptInfo, this);
}
}
diff --git a/AMDiS/src/nonlin/ProblemNonLin.h b/AMDiS/src/nonlin/ProblemNonLin.h
index 0158d4cb6c8970561f8b8dfebf32ab96651a5ee8..f370b41d0ff98d9833de55a80b41e6d5a9ae9f27 100644
--- a/AMDiS/src/nonlin/ProblemNonLin.h
+++ b/AMDiS/src/nonlin/ProblemNonLin.h
@@ -27,7 +27,7 @@
#include "NonLinSolver.h"
#include "DOFVector.h"
#include "SystemVector.h"
-#include "MatrixVector.h"
+#include <vector>
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
#ifdef HAVE_PARALLEL_MTL4
@@ -37,6 +37,8 @@
#endif
#endif
+using namespace std;
+
namespace AMDiS {
/**
@@ -54,7 +56,8 @@ namespace AMDiS {
nonLinSolver(NULL)
{
u0.resize(nComponents);
- u0.set(NULL);
+ for (int i = 0; i < nComponents; i++)
+ u0[i] = NULL;
}
/// Sets \ref u0 and interpolates it to \ref solution.
@@ -84,12 +87,9 @@ namespace AMDiS {
* Overrides ProblemStat::solve(). Uses the non linear solver
* \ref nonLinSolver.
*/
- virtual void solve(AdaptInfo *adaptInfo, bool fixedMatrix = false);
-
- /// Overrides ProblemStat::buildAfterCoarsen().
- virtual void buildAfterCoarsen(AdaptInfo *adaptInfo, Flag,
- bool assembleMatrix = true,
- bool assembleVector = true);
+ void solve(AdaptInfo *adaptInfo,
+ bool createMatrixData = true,
+ bool storeMatrixData = false);
/// Returns \ref nonLinSolver.
inline NonLinSolver *getNonLinSolver()
@@ -105,7 +105,7 @@ namespace AMDiS {
protected:
/// Initial guess for the solution.
- Vector<AbstractFunction<double, WorldVector<double> >*> u0;
+ vector<AbstractFunction<double, WorldVector<double> >*> u0;
/// Non linear solver object. Used in \ref solve().
NonLinSolver *nonLinSolver;