Allow to create the tube segments in any Euclidean axis

dune/fufem/makering.hh

@@ -229,8 +229,7 @@ void makeRingSegment2D(const Dune::FieldVector<typename GridType::ctype,2>& cent
 
 
 /** \todo Make tube orientation arbitrary. */
-/** \todo Allow arbitrary fromAngle <= toAngle. */
-/** \brief Create a 3D tube segment prism grid along the x-axis.
+/** \brief Create a 3D tube segment prism grid along a prescribed Euclidean axis.
  *
  *  If the tube is very thin the boundary parametrization can require the improvement of the inner grid vertices.
  *  This is for example implemented in the dune-biomech/prismlayer.hh class.
@@ -246,6 +245,7 @@ void makeRingSegment2D(const Dune::FieldVector<typename GridType::ctype,2>& cent
  *  \param nElementRing - The number of initial elements in the ring direction.
  *  \param nElementLength - The number of initial elements in the length direction.
  *  \param closeTube - Close the tube, needs fromAngle == 0, toAngle == 2*PI
+ *  \param axis - The axis along the tube should be build 0 <= axis <= 2
  *  \param innerBoundary - If given, make a level 0 boundary patch containing the inner nodes
  *  \param outerBoundary - If given, make a level 0 boundary patch containing the outer nodes
  *
@@ -256,6 +256,7 @@ GridType* makeTubeSegment3D(const Dune::FieldVector<typename GridType::ctype,3>&
                        const typename GridType::ctype innerRadius,
                        const typename GridType::ctype fromAngle, const typename GridType::ctype toAngle,
                        unsigned int nElementRing, unsigned int nElementLength, bool closeTube,
+                       unsigned int axis = 0,
                        bool tetrahedra=true, bool parametrizedBoundary=true,
                        BoundaryPatch<typename GridType::LevelGridView>* innerBoundary = NULL,
                        BoundaryPatch<typename GridType::LevelGridView>* outerBoundary = NULL)
@@ -263,9 +264,8 @@ GridType* makeTubeSegment3D(const Dune::FieldVector<typename GridType::ctype,3>&
     static_assert(GridType::dimension==3, "makeTubeSegment3D can only be called for 3-dimensional grids");
     typedef typename GridType::ctype ctype;
 
-    ctype eps = 1e-10;
-    if ( (-eps >= fromAngle) || (fromAngle > toAngle) || (toAngle > 2*M_PI) )
-        DUNE_THROW(Dune::Exception, "Parameters do not fulfill 0 <= fromAngle <= toAngle <= 2*M_PI");
+    if (fromAngle >= toAngle)
+        DUNE_THROW(Dune::Exception, "Parameters do not fulfill fromAngle < toAngle");
 
     bool fullTube = closeTube && (fromAngle==0.0) && (toAngle == 2*M_PI);
 
@@ -281,8 +281,8 @@ GridType* makeTubeSegment3D(const Dune::FieldVector<typename GridType::ctype,3>&
 
     // project center on the yz-plane
     FieldVector<ctype,2> localCenter;
-    localCenter[0] = center[1];
-    localCenter[1] = center[2];
+    localCenter[0] = center[(axis+1)%3];
+    localCenter[1] = center[(axis+2)%3];
 
     ArcOfCircle innerArc(localCenter, innerRadius, fromAngle, toAngle);
 
@@ -295,12 +295,12 @@ GridType* makeTubeSegment3D(const Dune::FieldVector<typename GridType::ctype,3>&
         FieldVector<ctype,3> vertex(0);
 
         // x-coordinate is determined by the center x-coord
-        vertex[0] = center[0] + j*stepLength;
+        vertex[axis] = center[axis] + j*stepLength;
 
         for (size_t i=0; i<nRingNodes; i++) {
 
             FieldVector<ctype,2> local = innerArc(i*stepRing);
-            vertex[1] = local[0]; vertex[2] = local[1];
+            vertex[(axis+1)%3] = local[0]; vertex[(axis+2)%3] = local[1];
 
             factory.insertVertex(vertex);
         }
@@ -314,12 +314,12 @@ GridType* makeTubeSegment3D(const Dune::FieldVector<typename GridType::ctype,3>&
         FieldVector<ctype,3> vertex(0);
 
         // x-coordinate is determined by the center x-coord
-        vertex[0] = center[0] + j*stepLength;
+        vertex[axis] = center[axis] + j*stepLength;
 
         for (size_t i=0; i<nRingNodes; i++) {
 
             FieldVector<ctype,2> local = outerArc(i*stepRing);
-            vertex[1] = local[0]; vertex[2] = local[1];
+            vertex[(axis+1)%3] = local[0]; vertex[(axis+2)%3] = local[1];
 
             factory.insertVertex(vertex);
         }