new gridConstructor for cuboids, uniform refinement, towards foam production

dune/tectonic/factories/twoblocksfactory.cc

@@ -2,6 +2,8 @@
 #include "config.h"
 #endif
 
+#include <string>
+
 #include <dune/fufem/geometry/convexpolyhedron.hh>
 
 #include <dune/contact/projections/normalprojection.hh>
@@ -48,39 +50,23 @@ void TwoBlocksFactory<HostGridType, VectorTEMPLATE>::setBodies() {
 #error CuboidGeometry only supports 2D and 3D!"
 #endif
 
-    // set up reference grids
-    gridConstructor_ = std::make_unique<GridsConstructor<HostGridType>>(cuboidGeometries_);
-    auto& grids = gridConstructor_->getGrids();
-
-    std::array<double, 2> smallestDiameter = {this->parset_.template get<double>("body0.smallestDiameter"), this->parset_.template get<double>("body1.smallestDiameter")};
+    const auto gravity = this->parset_.template get<double>("general.gravity");
 
     for (size_t i=0; i<this->bodyCount_; i++) {
+        std::string body_str = "body" + std::to_string(i);
         const auto& cuboidGeometry = *cuboidGeometries_[i];
 
-        // define weak patch and refine grid
-        const auto& weakeningRegions = cuboidGeometry.weakeningRegions();
-        for (size_t j=0; j<weakeningRegions.size(); j++) {
-            refine(*grids[i], weakeningRegions[j], smallestDiameter[i], CuboidGeometry::lengthScale());
-        }
+        std::array<unsigned int, 2> initialElements = {
+            { this->parset_.template get<size_t>("body" + std::to_string(i) + ".initialXElements"),
+              this->parset_.template get<size_t>("body" + std::to_string(i) + ".initialYElements") } };
 
-        // determine minDiameter and maxDiameter
-        double minDiameter = std::numeric_limits<double>::infinity();
-        double maxDiameter = 0.0;
-        for (auto &&e : elements(grids[i]->leafGridView())) {
-          auto const geometry = e.geometry();
-          auto const diam = diameter(geometry);
-          minDiameter = std::min(minDiameter, diam);
-          maxDiameter = std::max(maxDiameter, diam);
-        }
-        std::cout << "Grid" << i << " min diameter: " << minDiameter << std::endl;
-        std::cout << "Grid" << i << " max diameter: " << maxDiameter << std::endl;
-    }
-
-    bodyData_[0] = std::make_shared<MyBodyData<dim>>(this->parset_.sub("body0"), this->parset_.template get<double>("general.gravity"), zenith_());
-    this->bodies_[0] = std::make_shared<typename Base::LeafBody>(bodyData_[0], grids[0]);
+        // set up reference grid
+        DefaultCuboidGridConstructor<HostGridType> gridConstructor(cuboidGeometry, initialElements);
+        gridConstructor.refine(this->parset_.template get<size_t>(body_str + ".refinements"));
 
-    bodyData_[1] = std::make_shared<MyBodyData<dim>>(this->parset_.sub("body1"), this->parset_.template get<double>("general.gravity"), zenith_());
-    this->bodies_[1] = std::make_shared<typename Base::LeafBody>(bodyData_[1], grids[1]);
+        bodyData_[i] = std::make_shared<MyBodyData<dim>>(this->parset_.sub(body_str), gravity, zenith_());
+        this->bodies_[i] = std::make_shared<typename Base::LeafBody>(bodyData_[i], gridConstructor.grid());
+    }
 }
 
 template <class HostGridType, class VectorTEMPLATE>
@@ -121,6 +107,7 @@ void TwoBlocksFactory<HostGridType, VectorTEMPLATE>::setBoundaryConditions() {
 
     using Function = Dune::VirtualFunction<double, double>;
     std::shared_ptr<Function> neumannFunction = std::make_shared<NeumannCondition>();
+    //std::shared_ptr<Function> velocityDirichletFunction = std::make_shared<VelocityDirichletConditionLinearLoading>(this->parset_.template get<double>("boundary.dirichlet.finalVelocity"), 0.25);
     std::shared_ptr<Function> velocityDirichletFunction = std::make_shared<VelocityDirichletCondition>(this->parset_.template get<double>("boundary.dirichlet.finalVelocity"), 0.25);
     //std::shared_ptr<Function> velocityDirichletFunction = std::make_shared<VelocityStepDirichletCondition>(this->parset_.template get<double>("boundary.dirichlet.finalVelocity"), 10*this->parset_.template get<double>("boundary.dirichlet.finalVelocity"), 0.25, 0.5);
 

dune/tectonic/factories/twoblocksfactory.hh

@@ -10,8 +10,8 @@
 #include "contactnetworkfactory.hh"
 
 #include "../problem-data/mybody.hh"
-#include "../problem-data/grid/mygrids.hh"
 #include "../problem-data/grid/cuboidgeometry.hh"
+#include "../problem-data/grid/cuboidgridconstructor.hh"
 
 template <class HostGridType, class VectorType> class TwoBlocksFactory : public ContactNetworkFactory<HostGridType, VectorType>{
 private:
@@ -39,8 +39,6 @@ template <class HostGridType, class VectorType> class TwoBlocksFactory : public
 
     std::vector<std::shared_ptr<MyBodyData<dim>>> bodyData_;     // material properties of bodies
 
-    std::unique_ptr<GridsConstructor<HostGridType>> gridConstructor_;
-
     std::vector<std::shared_ptr<CuboidGeometry>> cuboidGeometries_;
 
     std::vector<std::shared_ptr<LeafFaces>> leafFaces_;

dune/tectonic/problem-data/grid/CMakeLists.txt

@@ -6,6 +6,7 @@ add_custom_target(tectonic_dune_problem-data_grid SOURCES
   cubegridconstructor.hh
   cuboidgeometry.hh
   cuboidgeometry.cc
+  cuboidgridconstructor.hh
   gridconstructor.hh
   mygrids.hh
   mygrids.cc
@@ -23,6 +24,7 @@ install(FILES
   cubefaces.hh
   cubegridconstructor.hh
   cuboidgeometry.hh
+  cuboidgridconstructor.hh
   gridconstructor.hh
   mygrids.hh
   simplexmanager.hh

dune/tectonic/problem-data/grid/cuboidgridconstructor.hh

+#ifndef SRC_CUBOIDGRIDCONSTRUCTOR_HH
+#define SRC_CUBOIDGRIDCONSTRUCTOR_HH
+
+#include <dune/common/fmatrix.hh>
+#include <dune/grid/common/gridfactory.hh>
+#include <dune/grid/utility/structuredgridfactory.hh>
+
+#include <dune/fufem/boundarypatch.hh>
+#include <dune/fufem/geometry/convexpolyhedron.hh>
+
+#include "cuboidgeometry.hh"
+#include "gridconstructor.hh"
+#include "simplexmanager.hh"
+
+// -------------------------------------------------
+// MyFaces
+// -------------------------------------------------
+
+template <class GridView> struct MyFaces {
+  BoundaryPatch<GridView> lower;
+  BoundaryPatch<GridView> right;
+  BoundaryPatch<GridView> upper;
+  BoundaryPatch<GridView> left;
+
+#if MY_DIM == 3
+  BoundaryPatch<GridView> front;
+  BoundaryPatch<GridView> back;
+#endif
+
+  MyFaces(GridView const &gridView, const CuboidGeometry<typename GridView::ctype>& cuboidGeometry_) :
+    #if MY_DIM == 3
+          lower(gridView),
+          right(gridView),
+          upper(gridView),
+          left(gridView),
+          front(gridView),
+          back(gridView),
+    #else
+          lower(gridView),
+          right(gridView),
+          upper(gridView),
+          left(gridView),
+    #endif
+          cuboidGeometry(cuboidGeometry_)
+    {
+      lower.insertFacesByProperty([&](typename GridView::Intersection const &in) {
+        return xyBetween(cuboidGeometry.lowerLeft(), cuboidGeometry.lowerRight(), in.geometry().center());
+      });
+
+      right.insertFacesByProperty([&](typename GridView::Intersection const &in) {
+        return xyBetween(cuboidGeometry.lowerRight(), cuboidGeometry.upperRight(), in.geometry().center());
+      });
+
+      upper.insertFacesByProperty([&](typename GridView::Intersection const &in) {
+        return xyBetween(cuboidGeometry.upperLeft(), cuboidGeometry.upperRight(), in.geometry().center());
+      });
+
+      left.insertFacesByProperty([&](typename GridView::Intersection const &in) {
+        return xyBetween(cuboidGeometry.lowerLeft(), cuboidGeometry.upperLeft(), in.geometry().center());
+      });
+
+    #if MY_DIM == 3
+      front.insertFacesByProperty([&](typename GridView::Intersection const &in) {
+        return isClose(cuboidGeometry.depth() / 2.0, in.geometry().center()[2]);
+      });
+
+      back.insertFacesByProperty([&](typename GridView::Intersection const &in) {
+        return isClose(-cuboidGeometry.depth() / 2.0, in.geometry().center()[2]);
+      });
+    #endif
+  }
+
+private:
+  const CuboidGeometry<typename GridView::ctype>& cuboidGeometry;
+
+  bool isClose(double a, double b) {
+    return std::abs(a - b) < 1e-14 * cuboidGeometry.lengthScale();
+  }
+
+  bool isClose2(double a, double b) {
+    return std::abs(a - b) <
+           1e-14 * cuboidGeometry.lengthScale() * cuboidGeometry.lengthScale();
+  }
+
+  template <class Vector>
+  bool xyBoxed(Vector const &v1, Vector const &v2, Vector const &x) {
+      auto const minmax0 = std::minmax(v1[0], v2[0]);
+      auto const minmax1 = std::minmax(v1[1], v2[1]);
+
+      if (minmax0.first - 1e-14 * cuboidGeometry.lengthScale() > x[0] or
+          x[0] > minmax0.second + 1e-14 * cuboidGeometry.lengthScale())
+        return false;
+      if (minmax1.first - 1e-14 * cuboidGeometry.lengthScale() > x[1] or
+          x[1] > minmax1.second + 1e-14 * cuboidGeometry.lengthScale())
+        return false;
+
+      return true;
+  }
+
+  template <class Vector>
+  bool xyCollinear(Vector const &a, Vector const &b, Vector const &c) {
+      return isClose2((b[0] - a[0]) * (c[1] - a[1]), (b[1] - a[1]) * (c[0] - a[0]));
+  }
+
+  template <class Vector>
+  bool xyBetween(Vector const &v1, Vector const &v2, Vector const &x) {
+      return xyCollinear(v1, v2, x) && xyBoxed(v1, v2, x);
+  }
+};
+
+
+template <class GridView>
+MyFaces<GridView> constructFaces(const GridView& gridView, const CuboidGeometry<typename GridView::ctype>& cuboidGeometry) {
+  return MyFaces<GridView>(gridView, cuboidGeometry);
+}
+
+// -------------------------------------------------
+// CuboidGridConstructor
+// -------------------------------------------------
+
+template <class Grid>
+class CuboidGridConstructor : public GridConstructor<Grid> {
+public:
+  CuboidGridConstructor(const CuboidGeometry<typename Grid::ctype>& cuboidGeometry) :
+    cuboidGeometry_(cuboidGeometry) {
+
+    createGrid();
+  }
+
+protected:
+  void createGrid() {
+    Dune::GridFactory<Grid> gridFactory;
+
+    const auto& A = cuboidGeometry_.lowerLeft();
+    const auto& B = cuboidGeometry_.lowerRight();
+            const auto& C = cuboidGeometry_.upperRight();
+            const auto& D = cuboidGeometry_.upperLeft();
+
+            unsigned int const vc = 4;
+
+            #if MY_DIM == 3
+                Dune::FieldMatrix<double, 2 * vc, MY_DIM> vertices;
+            #else
+                Dune::FieldMatrix<double, vc, MY_DIM> vertices;
+            #endif
+                for (size_t i = 0; i < 2; ++i) {
+            #if MY_DIM == 3
+                  size_t numXYplanes = 2;
+            #else
+                  size_t numXYplanes = 1;
+            #endif
+                  size_t k = 0;
+                  for (size_t j = 1; j <= numXYplanes; ++j) {
+                    vertices[k++][i] = A[i];
+                    vertices[k++][i] = B[i];
+                    vertices[k++][i] = C[i];
+                    vertices[k++][i] = D[i];
+                    assert(k == j * vc);
+                  }
+                }
+
+            #if MY_DIM == 3
+                for (size_t k = 0; k < vc; ++k) {
+                  vertices[k][2] = -cuboidGeometry_.depth() / 2.0;
+                  vertices[k + vc][2] = cuboidGeometry_.depth() / 2.0;
+                }
+            #endif
+
+                for (size_t i = 0; i < vertices.N(); ++i)
+                  gridFactory.insertVertex(vertices[i]);
+
+                Dune::GeometryType cell;
+            #if MY_DIM == 3
+                cell = Dune::GeometryTypes::tetrahedron;
+            #else
+                cell = Dune::GeometryTypes::triangle;
+            #endif
+
+            #if MY_DIM == 3
+                SimplexManager sm(vc);
+            #else
+                SimplexManager sm;
+            #endif
+                sm.addFromVerticesFFB(1, 2, 0);
+                sm.addFromVerticesFFB(2, 3, 0);
+                auto const &simplices = sm.getSimplices();
+
+                // sanity-check choices of simplices
+                for (size_t i = 0; i < simplices.size(); ++i) {
+                  Dune::FieldMatrix<double, MY_DIM, MY_DIM> check;
+                  for (size_t j = 0; j < MY_DIM; ++j)
+                    check[j] = vertices[simplices[i][j + 1]] - vertices[simplices[i][j]];
+                  assert(check.determinant() > 0);
+                  gridFactory.insertElement(cell, simplices[i]);
+                }
+
+                this->grid_ = std::shared_ptr<Grid>(gridFactory.createGrid());
+                this->grid_->setRefinementType(Grid::RefinementType::COPY);
+              }
+
+    protected:
+        const CuboidGeometry<typename Grid::ctype>& cuboidGeometry_;
+};
+
+
+// -------------------------------------------------
+// DefaultCuboidGridConstructor
+// -------------------------------------------------
+template <class Grid>
+class DefaultCuboidGridConstructor : public GridConstructor<Grid> {
+  static const int dim = Grid::dimension;
+
+public:
+  DefaultCuboidGridConstructor(const CuboidGeometry<typename Grid::ctype>& cuboidGeometry, std::array<unsigned int, dim> elements) :
+    cuboidGeometry_(cuboidGeometry),
+    elements_(elements) {
+
+    createGrid();
+  }
+
+protected:
+  void createGrid() {
+       this->grid_ = Dune::StructuredGridFactory<Grid>::createSimplexGrid(
+           cuboidGeometry_.lowerLeft(), cuboidGeometry_.upperRight(), elements_);
+    }
+
+    protected:
+        const CuboidGeometry<typename Grid::ctype>& cuboidGeometry_;
+        std::array<unsigned int, dim> elements_;
+};
+
+#endif
+
+

dune/tectonic/problem-data/grid/gridconstructor.hh

@@ -8,6 +8,7 @@
 
 #include "../../utils/diameter.hh"
 
+
 template <class GridType>
 class GridConstructor {
     private:
@@ -50,10 +51,14 @@ class GridConstructor {
             }
         }
 
+        void refine(size_t globalRefinements) {
+            grid_->globalRefine(globalRefinements);
+        }
+
+    protected:
         virtual void createGrid() = 0;
 
-    private:
-        Dune::GridFactory<GridType> gridFactory_;
+    protected:
         std::shared_ptr<GridType> grid_;
 };
 

src/foam/CMakeLists.txt

@@ -18,7 +18,7 @@ set(FOAM_SOURCE_FILES
   #../../dune/tectonic/io/hdf5/surface-writer.cc
   #../../dune/tectonic/io/hdf5/time-writer.cc
   ../../dune/tectonic/problem-data/grid/cuboidgeometry.cc
-  ../../dune/tectonic/problem-data/grid/mygrids.cc
+  #../../dune/tectonic/problem-data/grid/mygrids.cc
   ../../dune/tectonic/problem-data/grid/simplexmanager.cc
   ../../dune/tectonic/spatial-solving/solverfactory.cc
   ../../dune/tectonic/spatial-solving/fixedpointiterator.cc

src/foam/foam-2D.cfg

 # -*- mode:conf -*-
 [body0]
-smallestDiameter = 6.25e-2  # 2e-3 [m]
+refinements = 4  # 2e-3 [m]
+initialXElements = 5
+initialYElements = 1
 
 [body1]
-smallestDiameter = 6.25e-2 # 2e-3 [m]
+refinements = 4 # 2e-3 [m]
+initialXElements = 5
+initialYElements = 1
 
 [timeSteps]
 refinementTolerance = 1e-5 # 1e-5

src/foam/foam.cc

@@ -66,7 +66,7 @@
 
 #include <dune/tectonic/problem-data/bc.hh>
 #include <dune/tectonic/problem-data/mybody.hh>
-#include <dune/tectonic/problem-data/grid/mygrids.hh>
+//#include <dune/tectonic/problem-data/grid/mygrids.hh>
 
 #include <dune/tectonic/spatial-solving/tnnmg/functional.hh>
 //#include <dune/tectonic/spatial-solving/preconditioners/multilevelpatchpreconditioner.hh>
@@ -473,8 +473,6 @@ int main(int argc, char *argv[]) {
     typename ContactNetwork::ExternalForces externalForces;
     contactNetwork.externalForces(externalForces);
 
-    auto&& noFriction = ZeroNonlinearity();
-
     StepBase<NonlinearFactory, std::decay_t<decltype(contactNetwork)>, Updaters, std::decay_t<decltype(stateEnergyNorms)>>
         stepBase(parset, contactNetwork, totalDirichletNodes, globalFriction, frictionNodes,
                  externalForces, stateEnergyNorms);