This and that.

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

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]);

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);
 

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));
+  }
+
 }

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;
     }

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);
   }
 
 }

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;