New function for traversing sub boundaries.

AMDiS/src/ElInfo.h

@@ -424,9 +424,10 @@ namespace AMDiS {
     void fillElInfo(const MacroElement *mel, int refinementPathLength, unsigned long refinementPath);
 
     /** \brief
-     * calculates the Jacobian of the barycentric coordinates on \element and stores
-     * the matrix in grd_lam. The return value of the function is the absolute
-     * value of the determinant of the affine linear paraetrization's Jacobian.
+     * Calculates the Jacobian of the barycentric coordinates on \element and
+     * stores the matrix in grd_lam. The return value of the function is the
+     * absolute value of the determinant of the affine linear paraetrization's 
+     * Jacobian.
      * pure virtual => must be overriden in sub-class.
      */
     virtual double calcGrdLambda(DimVec<WorldVector<double> >& grd_lam) = 0;

AMDiS/src/Element.h

@@ -340,10 +340,9 @@ namespace AMDiS {
      * \{ 
      */
 
-    /// orient the vertices of edges/faces.
-    /// Used by Estimator for the jumps => same quadrature nodes from both sides!
-    virtual const FixVec<int,WORLD>& 
-      sortFaceIndices(int face, FixVec<int, WORLD> *vec) const = 0;
+    /// Orient the vertices of edges/faces. Used
+    /// by Estimator for the jumps => same quadrature nodes from both sides!
+    virtual void sortFaceIndices(int face, FixVec<int, WORLD> &vec) const = 0;
 
     /// Returns a copy of itself. Needed by Mesh to create Elements by 
     /// a prototype. 
@@ -432,11 +431,15 @@ namespace AMDiS {
      *                         identifing the DOF to be a vertex, edge, face or
      *                         center DOF.
      */
-    virtual void getHigherOrderDofs(const FiniteElemSpace* feSpace, 
-				    BoundaryObject bound,
-				    DofContainer& dofs,
-				    bool baseDofPtr = false,
-				    vector<GeoIndex>* dofGeoIndex = NULL) const = 0;
+    virtual void 
+    getHigherOrderDofs(const FiniteElemSpace* feSpace, 
+		       BoundaryObject bound,
+		       DofContainer& dofs,
+		       bool baseDofPtr = false,
+		       vector<GeoIndex>* dofGeoIndex = NULL) const = 0;
+
+    virtual void getSubBoundary(BoundaryObject bound, 
+				vector<BoundaryObject> &subBound) const = 0;
 
     /// Combines \ref getNodeDofs and \ref getHigherOrderDofs to one function. 
     /// See parameter description there.

AMDiS/src/Line.cc

@@ -45,13 +45,9 @@ namespace AMDiS {
     }  
   }
 
-  const FixVec<int, WORLD>& Line::sortFaceIndices(int face, 
-						  FixVec<int,WORLD> *vec) const
+  void Line::sortFaceIndices(int face, FixVec<int,WORLD> &vec) const
   {
-    static FixVec<int, WORLD> result(1, NO_INIT);
-    FixVec<int, WORLD> *v = vec ? vec : &result;
-    (*v)[0] = face;
-    return *v;
+    vec[0] = face;
   }
 
 }

AMDiS/src/Line.h

@@ -102,8 +102,7 @@ namespace AMDiS {
     }
 
     /// implements Element::sortFaceIndices
-    const FixVec<int,WORLD>& sortFaceIndices(int face, 
-					     FixVec<int,WORLD> *vec) const;
+    void sortFaceIndices(int face, FixVec<int,WORLD> &vec) const;
   
 
     /// implements Element::clone
@@ -187,6 +186,14 @@ namespace AMDiS {
       ERROR_EXIT("Not yet implemented!\n");
     }
 
+
+    void getSubBoundary(BoundaryObject bound, 
+			vector<BoundaryObject> &subBound) const
+    {
+      FUNCNAME("Line::getSubBoundary()");
+      ERROR_EXIT("Not yet implemented!\n");
+    }
+
   protected:
     /// vertexOfEdge[i][j] is the local number of the j-th vertex of the i-th 
     /// edge of this element.

AMDiS/src/Tetrahedron.cc

@@ -120,14 +120,11 @@ namespace AMDiS {
   }
 
 
-  const FixVec<int, WORLD>& Tetrahedron::sortFaceIndices(int face, 
-							 FixVec<int,WORLD> *vec) const
+  void Tetrahedron::sortFaceIndices(int face, FixVec<int,WORLD> &vec) const
   {
     static MatrixOfFixVecs<FixVec<int,WORLD> > *sorted_3d = NULL;
-    FixVec<int,WORLD> *val = NULL;
-    const int *vof = vertexOfFace[face];
 
-    if (NULL == sorted_3d) {
+    if (sorted_3d == NULL) {
       sorted_3d = new MatrixOfFixVecs<FixVec<int,WORLD> >(3, 4, 7, NO_INIT);
  
       (*sorted_3d)[0][0][0] = (*sorted_3d)[0][0][1] =
@@ -177,6 +174,7 @@ namespace AMDiS {
 	(*sorted_3d)[2][5][0] = (*sorted_3d)[2][6][1] = 3;
     }
 
+    const int *vof = vertexOfFace[face];
     int no = 0;
     if (dof[vof[0]][0] < dof[vof[1]][0])
       no++;
@@ -185,18 +183,10 @@ namespace AMDiS {
     if (dof[vof[2]][0] < dof[vof[0]][0])
       no += 4;
 
-    if (!(no >= 1  &&  no <= 6))
-      ERROR_EXIT("Cannot sort face indices of element %d at face %d\n", 
-		 getIndex(), face);
-
-    if (vec) {
-      val = vec;
-      *val = (*sorted_3d)[face][no];
-    } else {
-      val = &((*sorted_3d)[face][no]);
-    }
+    TEST_EXIT_DBG(no >= 1 && no <= 6)
+      ("Cannot sort face indices of element %d at face %d\n", getIndex(), face);
 
-    return(*(const_cast<const FixVec<int,WORLD>* >(val)));
+    vec = (*sorted_3d)[face][no];
   }
 
 
@@ -347,32 +337,6 @@ namespace AMDiS {
   }
 
 
-  void Tetrahedron::prepareNextBound(BoundaryObject &bound, int ithChild) const
-  {
-    FUNCNAME("Tetrahedron::prepareNextBound()");
-
-    switch (bound.subObj) {
-    case FACE:      
-      for (ExcludeList::iterator it = bound.excludedSubstructures.begin(); 
-   	   it != bound.excludedSubstructures.end(); ++it) {
-	if (it->first == EDGE && it->second != -1)
-	  it->second = edgeOfChild[bound.elType][ithChild][it->second];	
-      }
-      
-      bound.ithObj = sideOfChild[bound.elType][ithChild][bound.ithObj];
-      bound.elType = (bound.elType + 1) % 3;
-
-      break;
-    case EDGE:
-      bound.ithObj = edgeOfChild[bound.elType][ithChild][bound.ithObj];
-      bound.elType = (bound.elType + 1) % 3;
-      break;
-    default:
-      ERROR_EXIT("Should not happen!\n");
-    }
-  }
-
-
   void Tetrahedron::getHigherOrderDofs(const FiniteElemSpace* feSpace,
 				       BoundaryObject bound,
 				       DofContainer& dofs,
@@ -491,4 +455,36 @@ namespace AMDiS {
     }
   }
 
+
+  void Tetrahedron::getSubBoundary(BoundaryObject bound, 
+				   vector<BoundaryObject> &subBound) const
+  {
+  }
+
+
+  void Tetrahedron::prepareNextBound(BoundaryObject &bound, int ithChild) const
+  {
+    FUNCNAME("Tetrahedron::prepareNextBound()");
+
+    switch (bound.subObj) {
+    case FACE:      
+      for (ExcludeList::iterator it = bound.excludedSubstructures.begin(); 
+   	   it != bound.excludedSubstructures.end(); ++it) {
+	if (it->first == EDGE && it->second != -1)
+	  it->second = edgeOfChild[bound.elType][ithChild][it->second];	
+      }
+      
+      bound.ithObj = sideOfChild[bound.elType][ithChild][bound.ithObj];
+      bound.elType = (bound.elType + 1) % 3;
+
+      break;
+    case EDGE:
+      bound.ithObj = edgeOfChild[bound.elType][ithChild][bound.ithObj];
+      bound.elType = (bound.elType + 1) % 3;
+      break;
+    default:
+      ERROR_EXIT("Should not happen!\n");
+    }
+  }
+
 }

AMDiS/src/Tetrahedron.h

@@ -104,8 +104,7 @@ namespace AMDiS {
     bool hasSide(Element* sideElem) const;
 
     /// implements Element::sortFaceIndices
-    const FixVec<int,WORLD>& sortFaceIndices(int face, 
-					     FixVec<int,WORLD> *vec) const;
+    void sortFaceIndices(int face, FixVec<int,WORLD> &vec) const;
 
     /// Returns false because this element is a Tetrahedron.
     inline bool isLine() const 
@@ -157,6 +156,9 @@ namespace AMDiS {
 			    bool baseDofPtr = false,
 			    vector<GeoIndex>* dofGeoIndex = NULL) const;
 
+    void getSubBoundary(BoundaryObject bound, 
+			vector<BoundaryObject> &subBound) const;
+
     void prepareNextBound(BoundaryObject &bound, int ithChild) const;
 
   public:

AMDiS/src/Triangle.cc

@@ -51,16 +51,11 @@ namespace AMDiS {
     }
   }
 
-  const FixVec<int, WORLD>& Triangle::sortFaceIndices(int face, 
-						      FixVec<int, WORLD> *vec) const
+  void Triangle::sortFaceIndices(int face, FixVec<int, WORLD> &vec) const
   {
     static MatrixOfFixVecs<FixVec<int, WORLD> > *sorted_2d = NULL;
 
-    int no = 0;
-    FixVec<int, WORLD> *val = NULL;
-    const int *vof = vertexOfEdge[face];
-
-    if (NULL == sorted_2d) {
+    if (sorted_2d == NULL) {
       sorted_2d = new MatrixOfFixVecs<FixVec<int, WORLD> >(2, 3, 2, NO_INIT);
 
       (*sorted_2d)[1][0][1] = (*sorted_2d)[1][1][0] =
@@ -71,19 +66,9 @@ namespace AMDiS {
 	(*sorted_2d)[1][0][0] = (*sorted_2d)[1][1][1] = 2;
     }
 
-    if (dof[vof[0]][0] < dof[vof[1]][0])
-      no = 0;
-    else
-      no = 1;
-
-    if (vec) {
-      val = vec;
-      *val = (*sorted_2d)[face][no];
-    } else {
-      val = &((*sorted_2d)[face][no]);
-    }
-
-    return *(const_cast<const FixVec<int, WORLD>* >(val));
+    const int *vof = vertexOfEdge[face];
+    int no = ((dof[vof[0]][0] < dof[vof[1]][0]) ? 0 : 1);
+    vec = (*sorted_2d)[face][no];
   }
 
 
@@ -282,4 +267,49 @@ namespace AMDiS {
   }
 
 
+  void Triangle::getSubBoundary(BoundaryObject bound, 
+				vector<BoundaryObject> &subBound) const
+  {
+    TEST_EXIT_DBG(bound.subObj == EDGE)("This should not happen!\n");
+
+    if (!child[0]) {
+      subBound.push_back(bound);
+      return;
+    }
+
+    BoundaryObject nextBound0 = bound;
+    BoundaryObject nextBound1 = bound;
+    prepareNextBound(nextBound0, 0);
+    prepareNextBound(nextBound1, 1);
+
+    if (nextBound0.ithObj >= 0 && nextBound1.ithObj >= 0) {
+      if (bound.reverseMode) {
+	child[1]->getSubBoundary(nextBound1, subBound);
+	child[0]->getSubBoundary(nextBound0, subBound);
+      } else {
+	child[0]->getSubBoundary(nextBound0, subBound);
+	child[1]->getSubBoundary(nextBound1, subBound);
+      }
+    } else {
+      if (nextBound0.ithObj >= 0)
+	child[0]->getSubBoundary(nextBound0, subBound);
+
+      if (nextBound1.ithObj >= 0)
+	child[1]->getSubBoundary(nextBound1, subBound);
+    }
+  }
+
+  
+  void Triangle::prepareNextBound(BoundaryObject &bound, int ithChild) const
+  {
+    FUNCNAME("Triangle::prepareNextBound()");
+
+    TEST_EXIT_DBG(bound.el == this)("Wrong element!\n");
+    TEST_EXIT_DBG(child[0])("Has no child!\n");
+
+    bound.ithObj = sideOfChild[ithChild][bound.ithObj];
+    bound.el = child[ithChild];
+    bound.elIndex =  child[ithChild]->getIndex();
+  }
+
 }

AMDiS/src/Triangle.h

@@ -95,7 +95,7 @@ namespace AMDiS {
     bool hasSide(Element* sideElem) const;
 
     /// implements Element::sortFaceIndices
-    const FixVec<int, WORLD>& sortFaceIndices(int face, FixVec<int, WORLD> *vec) const;
+    void sortFaceIndices(int face, FixVec<int, WORLD> &vec) const;
 
     /// implements Element::isLine. Returns false because this element is a Triangle
     inline bool isLine() const 
@@ -205,6 +205,11 @@ namespace AMDiS {
 			    bool baseDofPtr = false,
 			    vector<GeoIndex>* dofGeoIndex = NULL) const;
 
+    void prepareNextBound(BoundaryObject &bound, int ithChild) const;
+
+    void getSubBoundary(BoundaryObject bound, 
+			vector<BoundaryObject> &subBound) const;
+
   protected:
     /// vertexOfEdge[i][j] is the local number of the j-th vertex of the i-th 
     /// edge of this element.

AMDiS/src/est/ResidualEstimator.cc

@@ -370,8 +370,8 @@ namespace AMDiS {
 
       int oppV = elInfo->getOppVertex(face);
 	
-      el->sortFaceIndices(face, &faceIndEl);
-      neigh->sortFaceIndices(oppV, &faceIndNeigh);	
+      el->sortFaceIndices(face, faceIndEl);
+      neigh->sortFaceIndices(oppV, faceIndNeigh);	
       neighInfo->setElement(const_cast<Element*>(neigh));
       neighInfo->setFillFlag(Mesh::FILL_COORDS);
       

AMDiS/src/est/ResidualEstimator.h

@@ -87,10 +87,8 @@ namespace AMDiS {
 
     virtual void init(double timestep);
 
-    /** \brief
-     * Estimates the error on an element. For more information about the
-     * parameter, see the description \ref Estimator::estimateElement.
-     */
+    /// Estimates the error on an element. For more information about the
+    /// parameter, see the description \ref Estimator::estimateElement.
     virtual void estimateElement(ElInfo *elInfo, DualElInfo *dualElInfo = NULL);
 
     virtual void exit(bool output = true);
@@ -114,12 +112,10 @@ namespace AMDiS {
 
     /// Constant in front of the time
     double C3;
-
-    /** \brief
-     * Is true, if C1 != 0, and C0 and C3 = 0, hence only the jump residual must be
-     * calculated. In this case, some optimizations can be done, because the jump
-     * residual is calculated only on second order terms.
-     */
+    
+    /// Is true, if C1 != 0, and C0 and C3 = 0, hence only the jump residual must be
+    /// calculated. In this case, some optimizations can be done, because the jump
+    /// residual is calculated only on second order terms.
     bool jumpResidualOnly;
 
     /// Number of systems, e.g., number of variables in the equation.
@@ -184,11 +180,9 @@ namespace AMDiS {
     /// Maximal number of neighbours an element may have in the used dimension.
     int nNeighbours;
 
-    /** \brief 
-     * Defines for every system if there are second order terms. These values
-     * are used to ommit computations of the jump residual that is defined
-     * only on second order terms.
-     */    
+    /// Defines for every system if there are second order terms. These values
+    /// are used to ommit computations of the jump residual that is defined
+    /// only on second order terms.
     std::vector<bool> secondOrderTerms;
   };
 }

AMDiS/src/est/SimpleResidualEstimator.cc

@@ -248,8 +248,8 @@ namespace AMDiS {
       // === The next lines are only to set all information about the ===
       // === neighbouring element.                                    ===
       int oppV = elInfo->getOppVertex(face);	
-      el->sortFaceIndices(face, &faceIndEl);
-      neigh->sortFaceIndices(oppV, &faceIndNeigh);	
+      el->sortFaceIndices(face, faceIndEl);
+      neigh->sortFaceIndices(oppV, faceIndNeigh);	
       neighInfo->setElement(const_cast<Element*>(neigh));
       neighInfo->setFillFlag(Mesh::FILL_COORDS);
       

AMDiS/src/parallel/InteriorBoundary.h

@@ -60,6 +60,21 @@ namespace AMDiS {
       return periodic;
     }
 
+    vector<AtomicBoundary>& getOwn(int rank)
+    {
+      return own[rank];
+    }
+ 
+    vector<AtomicBoundary>& getOther(int rank)
+    {
+      return other[rank];
+    }
+
+    vector<AtomicBoundary>& getPeriodic(int rank)
+    {
+      return periodic[rank];
+    }
+
     bool hasPeriodic()
     {
       return static_cast<bool>(periodic.size());

test/mpi/data/data0007/bound_rank_0.dat

+1885 2 1
+1884 2 0
+1883 2 1
+1882 2 0
+1877 2 1
+1876 2 0
+1875 2 1
+1874 2 0
+1841 2 1
+1840 2 0
+1839 2 1
+1838 2 0
+1833 2 1
+1832 2 0
+1831 2 1
+1830 2 0
+1673 2 1
+1672 2 0
+1671 2 1
+1670 2 0
+1665 2 1
+1664 2 0
+1663 2 1
+1662 2 0
+1629 2 1
+1628 2 0
+1627 2 1
+1626 2 0
+1621 2 1
+1620 2 0
+1619 2 1
+1618 2 0

test/mpi/data/data0007/bound_rank_1.dat

+642 2 0
+643 2 1
+644 2 0
+645 2 1
+650 2 0
+651 2 1
+652 2 0
+653 2 1
+686 2 0
+687 2 1
+688 2 0
+689 2 1
+694 2 0
+695 2 1
+696 2 0
+697 2 1
+888 2 0
+889 2 1
+890 2 0
+891 2 1
+897 2 1
+898 2 0
+899 2 1
+942 2 0
+943 2 1
+944 2 0
+945 2 1
+952 2 0
+953 2 1
+954 2 0
+955 2 1

test/mpi/data/data0007/bound_rank_2.dat

+1885 2 1
+1884 2 0
+1883 2 1
+1882 2 0
+1877 2 1
+1876 2 0
+1875 2 1
+1874 2 0
+1841 2 1
+1840 2 0
+1839 2 1
+1838 2 0
+1833 2 1
+1832 2 0
+1831 2 1
+1830 2 0
+1673 2 1
+1672 2 0
+1671 2 1
+1670 2 0
+1665 2 1
+1664 2 0
+1663 2 1
+1662 2 0
+1629 2 1
+1628 2 0
+1627 2 1
+1626 2 0
+1621 2 1
+1620 2 0
+1619 2 1
+1618 2 0

test/mpi/data/data0007/bound_rank_3.dat

+642 2 0
+643 2 1
+644 2 0
+645 2 1
+650 2 0
+651 2 1
+652 2 0
+653 2 1
+686 2 0
+687 2 1
+688 2 0
+689 2 1
+694 2 0
+695 2 1
+696 2 0
+697 2 1
+888 2 0
+889 2 1
+890 2 0
+891 2 1
+896 2 0
+897 2 1
+898 2 0
+899 2 1
+942 2 0
+943 2 1
+944 2 0
+945 2 1
+952 2 0
+953 2 1
+954 2 0
+955 2 1

test/mpi/init/test0007.dat.2d

+dimension of world:             2
+
+elliptMesh->macro file name:    ./macro/macro.stand.2d
+elliptMesh->global refinements: 10
+
+ellipt->mesh:                   elliptMesh
+ellipt->dim:                    2
+ellipt->components:             1
+ellipt->polynomial degree[0]:   1
+ 
+ellipt->solver:                 cg
+ellipt->solver->max iteration:  10
+ellipt->solver->tolerance:      1.e-8
+ellipt->solver->info:           10
+ellipt->solver->left precon:    diag
+ellipt->solver->right precon:   no
+
+ellipt->estimator[0]:              0
+ellipt->estimator[0]->error norm:  1   % 1: H1_NORM, 2: L2_NORM
+ellipt->estimator[0]->C0:          0.1 % constant of element residual
+ellipt->estimator[0]->C1:          0.1 % constant of jump residual
+
+ellipt->marker[0]->strategy:       0   % 0: no adaption 1: GR 2: MS 3: ES 4:GERS
+
+ellipt->output->filename:       output/ellipt.2d
+ellipt->output->ParaView format: 1
+
+parallel->log main rank:     0
+parallel->pre refine:        0
+parallel->partitioner:       simple

test/mpi/src/test0007.cc

+#define BOOST_TEST_DYN_LINK
+#define BOOST_TEST_MODULE 0007
+#define BOOST_TEST_NO_MAIN
+
+#include <boost/test/unit_test.hpp>
+#include <AMDiS.h>
+
+using namespace AMDiS;
+using namespace std;
+
+/*
+  - Test must run with 4 ranks.
+
+  - Test for function Triangle::getSubBoundary()
+  
+  - Test on edge boundary between rank 0/1 and 2/3
+*/
+
+BOOST_AUTO_TEST_CASE(amdis_mpi_triangle_getsubboundary)
+{
+  BOOST_REQUIRE(MPI::COMM_WORLD.Get_size() == 4);
+