From c44c3229d41dac2c970e449e1729558659a50837 Mon Sep 17 00:00:00 2001
From: podlesny <podlesny@zedat.fu-berlin.de>
Date: Fri, 19 Mar 2021 15:33:20 +0100
Subject: [PATCH] new gridConstructor for cuboids, uniform refinement, towards
foam production
---
dune/tectonic/factories/twoblocksfactory.cc | 41 ++-
dune/tectonic/factories/twoblocksfactory.hh | 4 +-
.../tectonic/problem-data/grid/CMakeLists.txt | 2 +
.../grid/cuboidgridconstructor.hh | 234 ++++++++++++++++++
.../problem-data/grid/gridconstructor.hh | 9 +-
src/foam/CMakeLists.txt | 2 +-
src/foam/foam-2D.cfg | 8 +-
src/foam/foam.cc | 4 +-
8 files changed, 266 insertions(+), 38 deletions(-)
create mode 100644 dune/tectonic/problem-data/grid/cuboidgridconstructor.hh
diff --git a/dune/tectonic/factories/twoblocksfactory.cc b/dune/tectonic/factories/twoblocksfactory.cc
index a60ace65..9b7834c9 100644
--- a/dune/tectonic/factories/twoblocksfactory.cc
+++ b/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);
diff --git a/dune/tectonic/factories/twoblocksfactory.hh b/dune/tectonic/factories/twoblocksfactory.hh
index a02cd74a..39447a4a 100644
--- a/dune/tectonic/factories/twoblocksfactory.hh
+++ b/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_;
diff --git a/dune/tectonic/problem-data/grid/CMakeLists.txt b/dune/tectonic/problem-data/grid/CMakeLists.txt
index b7abfed5..6effa5a3 100644
--- a/dune/tectonic/problem-data/grid/CMakeLists.txt
+++ b/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
diff --git a/dune/tectonic/problem-data/grid/cuboidgridconstructor.hh b/dune/tectonic/problem-data/grid/cuboidgridconstructor.hh
new file mode 100644
index 00000000..ac2c09c2
--- /dev/null
+++ b/dune/tectonic/problem-data/grid/cuboidgridconstructor.hh
@@ -0,0 +1,234 @@
+#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
+
+
diff --git a/dune/tectonic/problem-data/grid/gridconstructor.hh b/dune/tectonic/problem-data/grid/gridconstructor.hh
index e93d3c01..85f97a21 100644
--- a/dune/tectonic/problem-data/grid/gridconstructor.hh
+++ b/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_;
};
diff --git a/src/foam/CMakeLists.txt b/src/foam/CMakeLists.txt
index b38552f3..c228fda9 100644
--- a/src/foam/CMakeLists.txt
+++ b/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
diff --git a/src/foam/foam-2D.cfg b/src/foam/foam-2D.cfg
index d5028e7f..c7e12533 100644
--- a/src/foam/foam-2D.cfg
+++ b/src/foam/foam-2D.cfg
@@ -1,9 +1,13 @@
# -*- 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
diff --git a/src/foam/foam.cc b/src/foam/foam.cc
index 21dc1cec..b606d65b 100644
--- a/src/foam/foam.cc
+++ b/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);
--
GitLab