Integration over boundary parts

AMDiS/src/DOFVector.cc

@@ -160,6 +160,146 @@ namespace AMDiS {
   };
 
 
+  template<>
+  double DOFVector<double>::IntOnBoundary(BoundaryType boundaryType, 
+                                          Quadrature* q) const
+  {
+    FUNCNAME("DOFVector<T>::IntOnBoundary()");
+  
+    Mesh* mesh = this->feSpace->getMesh();
+    int dim = mesh->getDim();
+  
+    if (!q) {
+      int deg = 2 * this->feSpace->getBasisFcts()->getDegree();
+      q = Quadrature::provideQuadrature(dim - 1, deg);
+    } else {
+      TEST_EXIT(q->getDim() == dim-1)
+        ("Provided quadrature has wrong dimension for surfaceQuadrature!\n");
+    }
+  
+    // create barycentric coords for each vertex of each side
+    const Element *refElement = Global::getReferenceElement(dim);
+    VectorOfFixVecs<DimVec<double> >**coords;
+    coords = new VectorOfFixVecs<DimVec<double> >*[dim + 1];
+  
+    // for all element sides
+    for (int i = 0; i < dim + 1; i++) {
+      coords[i] = new VectorOfFixVecs<DimVec<double> >(dim, dim, DEFAULT_VALUE,
+        DimVec<double>(dim, DEFAULT_VALUE, 0.0));
+      // for each vertex of the side
+      for (int k = 0; k < dim; k++) {
+        int index = refElement->getVertexOfPosition(INDEX_OF_DIM(dim - 1, dim), 
+                                                    i, k);
+        (*coords[i])[k][index] = 1.0;
+      }
+    }
+  
+    std::vector<SurfaceQuadrature*> quadSurfaces(dim + 1);
+    for (int i = 0; i < dim + 1; i++)
+      quadSurfaces[i] = new SurfaceQuadrature(q, *coords[i]);
+  
+    double result = 0.0;
+    TraverseStack stack;
+    ElInfo *elInfo = stack.traverseFirst(mesh, -1, 
+      Mesh::CALL_LEAF_EL | Mesh::FILL_BOUND | Mesh::FILL_COORDS);
+    while (elInfo) {
+      for (int face = 0; face < dim + 1; face++) {
+        if (elInfo->getBoundary(face) == boundaryType) {
+          int nPoints = quadSurfaces[face]->getNumPoints();
+          mtl::dense_vector<double> uh_vec(nPoints);
+          double det = elInfo->calcSurfaceDet(*coords[face]);
+          double normT = 0.0;
+          this->getVecAtQPs(elInfo, quadSurfaces[face], NULL, uh_vec);
+          for (int iq = 0; iq < nPoints; iq++)
+            normT += quadSurfaces[face]->getWeight(iq) * (uh_vec[iq]);
+          result += det * normT;
+        }
+      }
+  
+      elInfo = stack.traverseNext(elInfo);
+    }
+  
+    #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
+    double localResult = result;
+    MPI::COMM_WORLD.Allreduce(&localResult, &result, 1, MPI_DOUBLE, MPI_SUM);
+    #endif
+    
+    return result;
+  }
+
+
+  template<>
+  double DOFVector<WorldVector<double> >::IntOnBoundaryNormal(
+    BoundaryType boundaryType, Quadrature* q) const
+  {
+    FUNCNAME("DOFVector<T>::IntOnBoundary()");
+ 
+    Mesh* mesh = this->feSpace->getMesh();
+    int dim = mesh->getDim();
+  
+    if (!q) {
+      int deg = 2 * this->feSpace->getBasisFcts()->getDegree();
+      q = Quadrature::provideQuadrature(dim - 1, deg);
+    } else {
+      TEST_EXIT(q->getDim() == dim-1)
+        ("Provided quadrature has wrong dimension for surfaceQuadrature!\n");
+    }
+  
+    // create barycentric coords for each vertex of each side
+    const Element *refElement = Global::getReferenceElement(dim);
+    VectorOfFixVecs<DimVec<double> >**coords;
+    coords = new VectorOfFixVecs<DimVec<double> >*[dim + 1];
+  
+    // for all element sides
+    for (int i = 0; i < dim + 1; i++) {
+      coords[i] = new VectorOfFixVecs<DimVec<double> >(dim, dim, DEFAULT_VALUE,
+        DimVec<double>(dim, DEFAULT_VALUE, 0.0));
+      // for each vertex of the side
+      for (int k = 0; k < dim; k++) {
+        int index = refElement->getVertexOfPosition(INDEX_OF_DIM(dim - 1, dim), 
+                                                    i, k);
+        (*coords[i])[k][index] = 1.0;
+      }
+    }
+  
+    std::vector<SurfaceQuadrature*> quadSurfaces(dim + 1);
+    for (int i = 0; i < dim + 1; i++)
+      quadSurfaces[i] = new SurfaceQuadrature(q, *coords[i]);
+  
+    double result = 0.0;
+    WorldVector<double> normal;
+    TraverseStack stack;
+    ElInfo *elInfo = stack.traverseFirst(mesh, -1, 
+      Mesh::CALL_LEAF_EL | Mesh::FILL_BOUND | Mesh::FILL_COORDS);
+    while (elInfo) {
+      for (int face = 0; face < dim + 1; face++) {
+        if (elInfo->getBoundary(face) == boundaryType) {
+          elInfo->getNormal(face, normal);
+          double det = elInfo->calcSurfaceDet(*coords[face]);
+  
+          int nPoints = quadSurfaces[face]->getNumPoints();
+          mtl::dense_vector<WorldVector<double> > uh_vec(nPoints);
+          WorldVector<double> normT; normT.set(0.0);
+          this->getVecAtQPs(elInfo, quadSurfaces[face], NULL, uh_vec);
+          for (int iq = 0; iq < nPoints; iq++)
+            normT += quadSurfaces[face]->getWeight(iq) * (uh_vec[iq]);
+          // scalar product between vector-valued solution and normal vector
+          result += det * (normT * normal); 
+        }
+      }
+  
+      elInfo = stack.traverseNext(elInfo);
+    }
+  
+    #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
+    double localResult = result;
+    MPI::COMM_WORLD.Allreduce(&localResult, &result, 1, MPI_DOUBLE, MPI_SUM);
+    #endif
+    
+    return result;
+  }
+
+
   template<>
   DOFVector<WorldVector<double> >*
   DOFVector<double>::getGradient(DOFVector<WorldVector<double> > *grad) const 

AMDiS/src/DOFVector.h

@@ -41,6 +41,7 @@
 #include "DOFMatrix.h" 
 #include "BasisFunction.h"
 #include "FiniteElemSpace.h"
+#include "SurfaceQuadrature.h"
 
 namespace AMDiS {
 
@@ -446,6 +447,29 @@ namespace AMDiS {
     /// Calculates Integral of this DOFVector
     double Int(Quadrature* q = NULL) const;
 
+    /** \brief
+     * Calculates Integral of this DOFVector over parts of the domain
+     * boundary, indicated by boundaryType. Implemented for DOFVector<double>
+     **/
+    T IntOnBoundary(BoundaryType boundary, Quadrature* q = NULL) const
+    {
+      FUNCNAME("DOFVector::IntOnBoundary())");
+      TEST_EXIT(false)("Please implement your integration\n");
+      return 0.0;
+    }
+
+    /** \brief
+     * Calculates Integral of this DOFVector times normal vector over parts 
+     * of the domain boundary, indicated by boundaryType. Implemented for 
+     * DOFVector<WorldVector<double> >
+     **/
+    double IntOnBoundaryNormal(BoundaryType boundary, Quadrature* q = NULL) const
+    {
+      FUNCNAME("DOFVector::IntOnBoundaryNormal())");
+      TEST_EXIT(false)("Please implement your integration\n");
+      return 0.0;
+    }
+
     /// Calculates L1 norm of this DOFVector
     double L1Norm(Quadrature* q = NULL) const;
  
@@ -544,7 +568,8 @@ namespace AMDiS {
      *  starts from oldElInfo. 
      *  implemented for: double, WorldVector< double >
      */
-    inline const T& evalAtPoint(WorldVector<double> &p, ElInfo *oldElInfo = NULL, T* value = NULL) const 
+    inline const T& evalAtPoint(
+      WorldVector<double> &p, ElInfo *oldElInfo = NULL, T* value = NULL) const 
     {
       FUNCNAME("DOFVector::evalAtPoint())");
       TEST_EXIT(false)("Please implement your evaluation\n");
@@ -575,6 +600,7 @@ namespace AMDiS {
     DOFVector<WorldVector<T> >* getRecoveryGradient(DOFVector<WorldVector<T> >*) const;
 
   protected: 
+
     /// Data container
     std::vector<T> vec; 
  
@@ -590,10 +616,20 @@ namespace AMDiS {
 
 
   template<>
-  const double& DOFVector<double>::evalAtPoint(WorldVector<double> &p, ElInfo *oldElInfo, double* value) const;
+  double DOFVector<double>::IntOnBoundary(
+    BoundaryType boundaryType, Quadrature* q) const;
+
+  template<>
+  double DOFVector<WorldVector<double> >::IntOnBoundaryNormal(
+    BoundaryType boundaryType, Quadrature* q) const;
+
+  template<>
+  const double& DOFVector<double>::evalAtPoint(
+    WorldVector<double> &p, ElInfo *oldElInfo, double* value) const;
 
   template<>
-  const WorldVector<double>& DOFVector<WorldVector<double> >::evalAtPoint(WorldVector<double> &p, ElInfo *oldElInfo, WorldVector<double>* value) const;
+  const WorldVector<double>& DOFVector<WorldVector<double> >::evalAtPoint(
+    WorldVector<double> &p, ElInfo *oldElInfo, WorldVector<double>* value) const;
 
   template<>
   void DOFVector<double>::refineInterpol(RCNeighbourList&, int);

AMDiS/src/DOFVector.hh

@@ -520,9 +520,8 @@ namespace AMDiS {
     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, -1, 
+      Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET);
     while (elInfo) {
       double det = elInfo->getDet();
       double normT = 0.0;

AMDiS/src/ElInfo.cc

@@ -156,6 +156,24 @@ namespace AMDiS {
   }
 
 
+  double ElInfo::calcSurfaceDet(VectorOfFixVecs<DimVec<double> > &surfVert) const
+  {
+    double surfDet;
+    int dim = surfVert[0].getSize() - 1;
+    FixVec<WorldVector<double>, VERTEX> worldCoords(dim - 1, NO_INIT);
+  
+    // transform barycentric coords to world coords
+    for (int i = 0; i < dim; i++) {
+      coordToWorld(surfVert[i], worldCoords[i]);
+    }
+  
+    // calculate determinant for surface
+    surfDet = calcDet(worldCoords);
+
+    return surfDet;
+  }
+
+
   void ElInfo::testFlag(const Flag& flags) const
   {    
     TEST_EXIT_DBG(fillFlag.isSet(flags))("flag not set\n");

AMDiS/src/ElInfo.h

@@ -368,6 +368,9 @@ namespace AMDiS {
      */
     double calcDet(const FixVec<WorldVector<double>, VERTEX> &coords) const;
 
+    /// from CFE_Integration
+    double calcSurfaceDet(VectorOfFixVecs<DimVec<double> > &surfVert) const;
+
     /// Checks whether flag is set in ElInfo's \ref fillFlag. If not, the program exits.
     void testFlag(const Flag& flag) const;
 

AMDiS/src/Initfile.cc

@@ -68,7 +68,7 @@ namespace AMDiS {
   /// Fill an initfile from an input stream
   void Initfile::read(istream& in) 
   {
-    unsigned line_length = 256;
+    unsigned line_length = 512;
     char swap[line_length];
     in.getline(swap, line_length);
     while (in.good()) {
@@ -76,11 +76,14 @@ namespace AMDiS {
       std::string sw(swap);
       size_t pos0 = sw.find_first_not_of(whitespaces);
 
-      if (pos0 != std::string::npos && sw[pos0] != '%' 
-          && sw[pos0] != '#' && sw[pos0] != 0) {
+      if (pos0 != std::string::npos
+          && sw[pos0] != '%' 
+          && sw[pos0] != '#'
+          && sw[pos0] != 0) {
         // parse line and extract map: tag->value
         Parser parser(sw);
         
+        // add parameter to map after variable replacement
         operator[](variableReplacement(InitfileInternal::trim(parser.name))) = variableReplacement(InitfileInternal::trim(parser.value));
       } else if (sw[pos0] == '#' 
           && static_cast<size_t>(sw.find("#include")) == pos0) {
@@ -130,7 +133,7 @@ namespace AMDiS {
       std::string varName = inputSwap.substr(posVarBegin + 1 , posVarEnd - posVarBegin - 1);
       std::string varParam = checkedGet(varName);
       // if varname is found in parameter list then replace variable by value
-      // otherwise throw tarNotFound exception
+      // otherwise throw tagNotFound exception
       std::string replaceName = inputSwap.substr(posVar , posVarEnd - posVar + (posVarBegin - posVar));
       inputSwap.replace(inputSwap.find(replaceName), replaceName.length(), varParam);
 
@@ -181,7 +184,7 @@ namespace AMDiS {
   {
     initIntern();
     for (Initfile::iterator it = singlett->begin(); it != singlett->end(); it++)
-    std::cout << (*it).first << " => " << (*it).second << std::endl;
+      std::cout << (*it).first << " => " << (*it).second << std::endl;
   };
 
 
@@ -189,7 +192,7 @@ namespace AMDiS {
   void Initfile::write(ostream& out)
   {
     for (Initfile::iterator it = begin() ; it!=end(); it++)
-    out << (*it).first << ": " << (*it).second << std::endl;	
+      out << (*it).first << ": " << (*it).second << std::endl;	
   }
 
 

AMDiS/src/ProblemTimeInterface.h

@@ -45,6 +45,9 @@ namespace AMDiS {
   public:
     virtual ~ProblemTimeInterface() {};
 
+    /// Called at the beginning of the adaption loop before any other call
+    virtual void initTimeInterface() {};
+
     /// Executes all needed operations when the simulation time changes.
     virtual void setTime(AdaptInfo *adaptInfo) = 0;