* Exact error computation

AMDiS/src/ProblemNonLin.cc

@@ -188,7 +188,7 @@ namespace AMDiS {
     
     // for all elements ...
     for (int i = 0; i < nComponents; i++) {
-      elInfo = stack.traverseFirst(componentMeshes_[i], -1, 
+      elInfo = stack.traverseFirst(componentMeshes[i], -1, 
 				   Mesh::CALL_LEAF_EL | 
 				   Mesh::FILL_BOUND |
 				   Mesh::FILL_COORDS |

AMDiS/src/ProblemVec.cc

@@ -20,6 +20,7 @@
 #include "RobinBC.h"
 #include "PeriodicBC.h"
 #include "Lagrange.h"
+#include "Flag.h"
 
 namespace AMDiS {
 
@@ -45,7 +46,7 @@ namespace AMDiS {
 	   adoptFlag.isSet(INIT_SYSTEM) ||
 	   adoptFlag.isSet(INIT_FE_SPACE))) {
 	meshes_ = adoptProblem->getMeshes();
-	componentMeshes_ = adoptProblem->componentMeshes_;
+	componentMeshes = adoptProblem->componentMeshes;
 	refinementManager_ = adoptProblem->refinementManager_;
 	coarseningManager_ = adoptProblem->coarseningManager_;
       }
@@ -180,7 +181,7 @@ namespace AMDiS {
   {
     FUNCNAME("ProblemVec::createMesh()");
 
-    componentMeshes_.resize(nComponents);
+    componentMeshes.resize(nComponents);
     std::map<int, Mesh*> meshForRefinementSet;
     char number[3];
 
@@ -203,7 +204,7 @@ namespace AMDiS {
 	meshForRefinementSet[refSet] = newMesh;
 	meshes_.push_back(newMesh);
       }
-      componentMeshes_[i] = meshForRefinementSet[refSet];
+      componentMeshes[i] = meshForRefinementSet[refSet];
     }
     switch(dim) {
     case 1:
@@ -243,17 +244,17 @@ namespace AMDiS {
 
       TEST_EXIT(componentSpaces[i] == NULL)("feSpace already created\n");
 
-      if (feSpaceMap[std::pair<Mesh*, int>(componentMeshes_[i], degree)] == NULL) {
+      if (feSpaceMap[std::pair<Mesh*, int>(componentMeshes[i], degree)] == NULL) {
 	FiniteElemSpace *newFESpace = 
 	  FiniteElemSpace::provideFESpace(NULL,
 					  Lagrange::getLagrange(dim, degree),
-					  componentMeshes_[i],
+					  componentMeshes[i],
 					  name_ + "->feSpace");
-	feSpaceMap[std::pair<Mesh*, int>(componentMeshes_[i], degree)] = newFESpace;
+	feSpaceMap[std::pair<Mesh*, int>(componentMeshes[i], degree)] = newFESpace;
 	feSpaces.push_back(newFESpace);
       }
       componentSpaces[i] = 
-	feSpaceMap[std::pair<Mesh*, int>(componentMeshes_[i], degree)];
+	feSpaceMap[std::pair<Mesh*, int>(componentMeshes[i], degree)];
     }
 
     // create dof admin for vertex dofs if neccessary
@@ -459,14 +460,14 @@ namespace AMDiS {
       std::vector< DOFVector<double>* > solutionList(nComponents);
 
       for (int i = 0; i < nComponents; i++) {
-	TEST_EXIT(componentMeshes_[0] == componentMeshes_[i])
+	TEST_EXIT(componentMeshes[0] == componentMeshes[i])
 	  ("All Meshes have to be equal to write a vector file.\n");
 
 	solutionList[i] = solution_->getDOFVector(i);
       }
 
       fileWriters_.push_back(NEW FileWriter(numberedName,
-					    componentMeshes_[0],
+					    componentMeshes[0],
 					    solutionList));
     }
 
@@ -481,7 +482,7 @@ namespace AMDiS {
 
       if (filename != "") {
 	fileWriters_.push_back(NEW FileWriter(numberedName, 
-					      componentMeshes_[i], 
+					      componentMeshes[i], 
 					      solution_->getDOFVector(i)));
       }
     }
@@ -549,18 +550,7 @@ namespace AMDiS {
 #endif
 
     if (computeExactError) {
-      for (int i = 0; i < nComponents; i++) {	
-	DOFVector<double> *tmp = NEW DOFVector<double>(componentSpaces[i], "tmp");
-	tmp->interpol(exactSolutionFcts[i]);
-	//	double t = max(abs(tmp->max()), abs(tmp->min()));
-	double t = tmp->absMax();
-	*tmp -= *(solution_->getDOFVector(i));
-	double l2Error = tmp->L2Norm();
-	double maxError = tmp->absMax() / t;
-	MSG("L2    error = %.8e\n", l2Error);
-	MSG("L-inf error = %.8e\n", maxError);
-	DELETE tmp;	
-      }						       
+      computeError(adaptInfo);
     } else {
       for (int i = 0; i < nComponents; i++) {
 	Estimator *scalEstimator = estimator_[i];
@@ -600,7 +590,7 @@ namespace AMDiS {
     Flag markFlag = 0;
     for (int i = 0; i < nComponents; i++) {
       if (marker[i]) {
-	markFlag |= marker[i]->markMesh(adaptInfo, componentMeshes_[i]);
+	markFlag |= marker[i]->markMesh(adaptInfo, componentMeshes[i]);
       } else {
 	WARNING("no marker for component %d\n", i);
       }
@@ -745,7 +735,7 @@ namespace AMDiS {
 	}
 	
 	TraverseStack stack;
-	ElInfo *elInfo = stack.traverseFirst(componentMeshes_[i], -1, assembleFlag);
+	ElInfo *elInfo = stack.traverseFirst(componentMeshes[i], -1, assembleFlag);
 	
 	while (elInfo) {
 	  if (useGetBound_) {
@@ -787,7 +777,7 @@ namespace AMDiS {
       	solution_->getDOFVector(i)->getBoundaryManager()->initVector(solution_->getDOFVector(i));
 
       TraverseStack stack;
-      ElInfo *elInfo = stack.traverseFirst(componentMeshes_[i], -1, assembleFlag);
+      ElInfo *elInfo = stack.traverseFirst(componentMeshes[i], -1, assembleFlag);
       while (elInfo) {
 	if (rhs_->getDOFVector(i)->getBoundaryManager())
 	  rhs_->getDOFVector(i)->getBoundaryManager()->
@@ -1018,5 +1008,61 @@ namespace AMDiS {
 
     solution_->deserialize(in);
   }
+
+  void ProblemVec::computeError(AdaptInfo *adaptInfo) 
+  {
+    FUNCNAME("ProblemVec::computeError()");
+
+    for (int i = 0; i < nComponents; i++) {		
+      TEST_EXIT(exactSolutionFcts[i])("No solution function given!\n");
+
+      // Compute the difference between exact and computed solution
+      DOFVector<double> *tmp = NEW DOFVector<double>(componentSpaces[i], "tmp");
+      tmp->interpol(exactSolutionFcts[i]);
+      double solMax = tmp->absMax();
+      *tmp -= *(solution_->getDOFVector(i));
+      
+      MSG("L2    error = %.8e\n", tmp->L2Norm());
+      MSG("L-inf error = %.8e\n", tmp->absMax() / solMax);
+      
+      adaptInfo->setEstSum(tmp->absMax() / solMax, i);
+      adaptInfo->setEstMax(tmp->absMax() / solMax, i);
+      
+      // To set element estimates, compute a vector with the difference
+      // between exact and computed solution for each DOF.
+      DOFVector<double> *sol = NEW DOFVector<double>(componentSpaces[i], "tmp");
+      sol->interpol(exactSolutionFcts[i]);
+      DOFVector<double>::Iterator it1(sol, USED_DOFS);
+      DOFVector<double>::Iterator it2(tmp, USED_DOFS);
+      for (it1.reset(), it2.reset(); !it1.end(); ++it1, ++it2) {
+	if ((abs(*it1) <= DBL_TOL) || (abs(*it2) <= DBL_TOL)) {
+	  *it2 = 0.0;
+	} else {
+	  *it2 = abs(*it2 / *it1);
+	}
+      }
+
+      // Compute estimate for every mesh element
+      Vector<DegreeOfFreedom> locInd(componentSpaces[i]->getBasisFcts()->getNumber());
+      TraverseStack stack;
+      ElInfo *elInfo = stack.traverseFirst(componentMeshes[i], -1, Mesh::CALL_LEAF_EL);
+      while (elInfo) {
+	componentSpaces[i]->getBasisFcts()->getLocalIndicesVec(elInfo->getElement(),
+							       componentSpaces[i]->getAdmin(),
+							       &locInd);
+	double estimate = 0.0;
+	for (int j = 0; j < componentSpaces[i]->getBasisFcts()->getNumber(); j++) {
+	  estimate += (*tmp)[locInd[j]];
+	}
+	elInfo->getElement()->setEstimation(estimate, i);
+	elInfo->getElement()->setMark(0);
+								
+	elInfo = stack.traverseNext(elInfo);
+      }  
+      
+      DELETE tmp;	
+      DELETE sol;
+    }						           
+  }
 }
  

AMDiS/src/ProblemVec.h

@@ -322,9 +322,9 @@ namespace AMDiS {
      */
     inline Mesh* getMesh(int comp) {
       FUNCNAME("ProblemVec::getMesh()");
-      TEST_EXIT(comp < static_cast<int>(componentMeshes_.size()) && comp >= 0)
+      TEST_EXIT(comp < static_cast<int>(componentMeshes.size()) && comp >= 0)
 	("invalid component number\n");
-      return componentMeshes_[comp]; 
+      return componentMeshes[comp]; 
     };
 
     /** \brief
@@ -523,6 +523,13 @@ namespace AMDiS {
       return fileWriters_;
     };
 
+  protected:
+    /** \brief
+     * If the exact solution is known, the problem can compute the exact
+     * error instead of the error estimation. This is done in this function.
+     */
+    void computeError(AdaptInfo *adaptInfo);
+
   protected:
     /** \brief
      * Name of this problem.
@@ -552,7 +559,7 @@ namespace AMDiS {
     /** \brief
      * Pointer to the meshes for the different problem components
      */
-    std::vector<Mesh*> componentMeshes_;
+    std::vector<Mesh*> componentMeshes;
 
     /** \brief
      * Responsible for element marking.