From 73dfe50e45b7a05580b4cd0b5a7efd5ec02c559f Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@dune-project.org>
Date: Mon, 5 Sep 2005 08:24:57 +0000
Subject: [PATCH] added a preliminary method to insert linear boundary segments
 into the coarse grid

[[Imported from SVN: r2730]]
---
 grid/uggrid.hh               |  10 +++
 grid/uggrid/ugfunctions.hh   |  22 ++++++
 grid/uggrid/ugfunctions3d.hh |  21 ++++++
 grid/uggrid/uggrid.cc        | 141 +++++++++++++++++++++++++++++++++++
 4 files changed, 194 insertions(+)

diff --git a/grid/uggrid.hh b/grid/uggrid.hh
index fc40741b0..14d8d6236 100644
--- a/grid/uggrid.hh
+++ b/grid/uggrid.hh
@@ -364,6 +364,16 @@ namespace Dune {
     /** \brief End the coarse grid creation process */
     void createend();
 
+    /** \brief Preliminary method to insert a linear boundary segment into a UG coarse grid
+        \param index The index number of the segment
+        \param vertices The vertices of the segment
+        \param userData Pointer handed over to the segment implementation method.
+        To be removed quickly!
+     */
+    void insertLinearSegment(int index,
+                             const std::vector<int> vertices,
+                             void* userData);
+
     /** \brief Insert an element into the coarse grid
         \param type The GeometryType of the new element
         \param vertices The vertices of the new element, using the DUNE numbering
diff --git a/grid/uggrid/ugfunctions.hh b/grid/uggrid/ugfunctions.hh
index 2f901435e..2cccf4935 100644
--- a/grid/uggrid/ugfunctions.hh
+++ b/grid/uggrid/ugfunctions.hh
@@ -33,6 +33,8 @@ namespace Dune {
 
     typedef UG2d::UserProcPtr UserProcPtr;
 
+    typedef UG2d::BndSegFuncPtr BndSegFuncPtr;
+
     enum {GM_REFINE_NOT_CLOSED = UG2d::GM_REFINE_NOT_CLOSED};
 
     enum {GM_COPY_ALL = UG2d::GM_COPY_ALL};
@@ -399,6 +401,26 @@ namespace Dune {
       return UG2d::CreateFormatCmd(argc, argv);
     }
 
+    /** \todo Remove the const_casts */
+    static void* CreateBoundarySegment(const char *name, int left, int right,
+                                       int index, int res,
+                                       int *point,
+                                       const double *alpha, const double *beta,
+                                       UG2d::BndSegFuncPtr boundarySegmentFunction,
+                                       void *userData) {
+      return UG2d::CreateBoundarySegment(const_cast<char*>(name),            // internal name of the boundary segment
+                                         left,                      //  id of left subdomain
+                                         right,                      //  id of right subdomain
+                                         index,                  // Index of the segment
+                                         UG2d::NON_PERIODIC,     // I don't know what this means
+                                         res,                      // Resolution, only for the UG graphics
+                                         point,
+                                         const_cast<double*>(alpha),
+                                         const_cast<double*>(beta),
+                                         boundarySegmentFunction,
+                                         userData);
+    }
+
   };
   //! \todo Please doc me!
   template <int codim, int dimworld>
diff --git a/grid/uggrid/ugfunctions3d.hh b/grid/uggrid/ugfunctions3d.hh
index f6ab2bb49..30779d580 100644
--- a/grid/uggrid/ugfunctions3d.hh
+++ b/grid/uggrid/ugfunctions3d.hh
@@ -30,6 +30,8 @@ namespace Dune {
 
     typedef UG3d::UserProcPtr UserProcPtr;
 
+    typedef UG3d::BndSegFuncPtr BndSegFuncPtr;
+
     enum {GM_REFINE_NOT_CLOSED = UG3d::GM_REFINE_NOT_CLOSED};
 
     enum {GM_COPY_ALL = UG3d::GM_COPY_ALL};
@@ -401,6 +403,25 @@ namespace Dune {
       return UG3d::CreateFormatCmd(argc, argv);
     }
 
+    /** \todo Remove the const_casts */
+    static void* CreateBoundarySegment(const char *name, int left, int right,
+                                       int index, int res,
+                                       int *point,
+                                       const double *alpha, const double *beta,
+                                       UG3d::BndSegFuncPtr boundarySegmentFunction,
+                                       void *userData) {
+      return UG3d::CreateBoundarySegment(const_cast<char*>(name),            // internal name of the boundary segment
+                                         left,                      //  id of left subdomain
+                                         right,                      //  id of right subdomain
+                                         index,                  // Index of the segment
+                                         UG3d::NON_PERIODIC,     // I don't know what this means
+                                         res,                      // Resolution, only for the UG graphics
+                                         point,
+                                         const_cast<double*>(alpha),
+                                         const_cast<double*>(beta),
+                                         boundarySegmentFunction,
+                                         userData);
+    }
   };
 
   // Specializations for dimworld==3
diff --git a/grid/uggrid/uggrid.cc b/grid/uggrid/uggrid.cc
index 1eb87a8cd..4376aeb23 100644
--- a/grid/uggrid/uggrid.cc
+++ b/grid/uggrid/uggrid.cc
@@ -104,6 +104,79 @@ namespace Dune {
 
 using namespace Dune;
 
+/** This method implements a linear function in order to be able to
+ *  work with straight line boundaries.
+ *  We interpret data as a DOUBLE* to the world coordinates of the
+ *  two endpoints.
+ */
+static int linearSegmentDescription2d(void *data, double *param, double *result)
+{
+  double a[2], b[2];
+  a[0] = ((double*)data)[0];
+  a[1] = ((double*)data)[1];
+  b[0] = ((double*)data)[2];
+  b[1] = ((double*)data)[3];
+
+  // linear interpolation
+  result[0] = a[0] + (*param)*(b[0]-a[0]);
+  result[1] = a[1] + (*param)*(b[1]-a[1]);
+
+  return 0;
+}
+
+/** This method implements a linear function in order to be able to
+ *  work with straight line boundaries.
+ *  We interpret data as a DOUBLE* to the world coordinates of the
+ *  three endpoints.
+ *
+ * \todo This should actually be replaced by using LinearSegments
+ * instead of BoundarySegments.  But LinearSegments are buggy in UG.
+ */
+static int linearSegmentDescription3d(void *data, double *param, double *result)
+{
+  Dune::FieldVector<double, 3> a,b,c,d;
+  a[0] = ((double*)data)[0];
+  a[1] = ((double*)data)[1];
+  a[2] = ((double*)data)[2];
+  b[0] = ((double*)data)[3];
+  b[1] = ((double*)data)[4];
+  b[2] = ((double*)data)[5];
+  c[0] = ((double*)data)[6];
+  c[1] = ((double*)data)[7];
+  c[2] = ((double*)data)[8];
+  d[0] = ((double*)data)[9];
+  d[1] = ((double*)data)[10];
+  d[2] = ((double*)data)[11];
+
+  // linear interpolation
+  for (int i=0; i<3; i++)
+    result[i] = a[i] + param[0]*(b[i]-a[i]) + param[1]*(d[i]-a[i])
+                + param[0]*param[1]*(a[i]+c[i]-b[i]-d[i]);
+
+  return 0;
+}
+
+// The same thing for triangles
+static int linearSegmentDescription3d_tri(void *data, double *param, double *result)
+{
+  Dune::FieldVector<double, 3> a,b,c;
+  a[0] = ((double*)data)[0];
+  a[1] = ((double*)data)[1];
+  a[2] = ((double*)data)[2];
+  b[0] = ((double*)data)[3];
+  b[1] = ((double*)data)[4];
+  b[2] = ((double*)data)[5];
+  c[0] = ((double*)data)[6];
+  c[1] = ((double*)data)[7];
+  c[2] = ((double*)data)[8];
+
+  // linear interpolation
+  for (int i=0; i<3; i++)
+    result[i] = a[i] + param[0]*(b[i]-a[i]) + param[1]*(c[i]-b[i]);
+
+  return 0;
+}
+
 //***********************************************************************
 //
 // --UGGrid
@@ -858,6 +931,74 @@ void Dune::UGGrid < dim, dimworld >::createend()
 
 }
 
+template <int dim, int dimworld>
+void Dune::UGGrid<dim, dimworld>::insertLinearSegment(int index,
+                                                      const std::vector<int> vertices,
+                                                      void* userData)
+{
+  /** \todo Make sure this is the current multigrid, so that CreateBoundarySegment
+      really inserts boundary segments into this grid.*/
+
+  // Copy the vertices into a C-style array
+  /** \todo Due to some UG weirdness, in 3d, CreateBoundarySegment always expects
+      this array to have four entries, even if only a triangular segment is
+      inserted.  If not, undefined values are will be introduced. */
+  int vertices_c_style[4] = {-1, -1, -1, -1};
+  for (unsigned int i=0; i<vertices.size(); i++)
+    vertices_c_style[i] = vertices[i];
+
+  // Create dummy parameter ranges
+  const double alpha[2] = {0, 0};
+  const double beta[2]  = {1, 1};
+
+  // Create some boundary segment name
+  char segmentName[20];
+  if(sprintf(segmentName, "BS %d", index) < 0)
+    DUNE_THROW(GridError, "sprintf returned error code!");
+
+  // Choose the method which implements the shape of the boundary segment
+  typename UG_NS<dim>::BndSegFuncPtr boundarySegmentFunction;
+  if (dim==3) {
+    if (vertices.size()==3)
+      boundarySegmentFunction = &linearSegmentDescription3d_tri;
+    else
+      boundarySegmentFunction = &linearSegmentDescription3d;
+  } else
+    boundarySegmentFunction = &linearSegmentDescription2d;
+
+  // Actually create the segment
+  if (UG_NS<dim>::CreateBoundarySegment(segmentName,            // internal name of the boundary segment
+                                        1,                      //  id of left subdomain
+                                        2,                      //  id of right subdomain
+                                        index,                  // Index of the segment
+                                        1,                      // Resolution, only for the UG graphics
+                                        vertices_c_style,
+                                        alpha,
+                                        beta,
+                                        boundarySegmentFunction,
+                                        userData)==NULL) {
+    DUNE_THROW(GridError, "Calling UG" << dim << "d::CreateBoundarySegment failed!");
+  }
+
+#if 0
+  // It would be a lot smarter to use this way of describing
+  // boundary segments.  But as of yet, UG crashes when using
+  // linear segments.
+  double paramCoords[3][2] = {{0,0}, {1,0}, {0,1}};
+  if (UG3d::CreateLinearSegment(segmentName,
+                                left,               /*id of left subdomain */
+                                right,              /*id of right subdomain*/
+                                i,                  /*id of segment*/
+                                4,                  // Number of corners
+                                point,
+                                paramCoords
+                                )==NULL)
+    DUNE_THROW(IOError, "Error calling CreateLinearSegment");
+
+#endif
+
+}
+
 template <int dim, int dimworld>
 void Dune::UGGrid<dim, dimworld>::insertElement(GeometryType type,
                                                 const std::vector<unsigned int>& vertices)
-- 
GitLab