From 0d66c91051011997db76b5398ecc77ebaddd9265 Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Thu, 9 Jan 2014 17:00:47 +0100
Subject: [PATCH] [Extend] Introduce geometry class
---
src/mygeometry.hh | 15 ++++++++++++
src/mygrid.hh | 34 +++++++++++++++++++++++++++
src/one-body-sample.cc | 53 +++++++++++-------------------------------
3 files changed, 63 insertions(+), 39 deletions(-)
create mode 100644 src/mygeometry.hh
create mode 100644 src/mygrid.hh
diff --git a/src/mygeometry.hh b/src/mygeometry.hh
new file mode 100644
index 00000000..9f73081e
--- /dev/null
+++ b/src/mygeometry.hh
@@ -0,0 +1,15 @@
+#ifndef MY_GEOMETRY_HH
+#define MY_GEOMETRY_HH
+
+struct MyGeometry {
+ MyGeometry() {}
+
+ Dune::FieldVector<double, 2> A = { 0, 0 };
+ Dune::FieldVector<double, 2> B = { 5, 0 };
+ Dune::FieldVector<double, 2> C = { 5, 1 };
+ Dune::FieldVector<double, 2> D = { 0, 1 };
+
+ Dune::FieldVector<double, 2> zenith = { 0, 1 };
+};
+
+#endif
diff --git a/src/mygrid.hh b/src/mygrid.hh
new file mode 100644
index 00000000..964c82dc
--- /dev/null
+++ b/src/mygrid.hh
@@ -0,0 +1,34 @@
+#include <dune/grid/common/gridfactory.hh>
+
+#include <dune/fufem/boundarypatch.hh>
+
+#include "mygeometry.hh"
+
+template <class Grid>
+std::shared_ptr<Grid> constructGrid(MyGeometry const &myGeometry) {
+ std::array<unsigned int, 2> elements = { { 5, 1 } };
+ return Dune::StructuredGridFactory<Grid>::createSimplexGrid(
+ myGeometry.A, myGeometry.C, elements);
+}
+
+template <class GridView, class MyGeometry> class MyFaces {
+private:
+ bool isClose(double a, double b) {
+ return std::abs(a - b) < 1e-14;
+ };
+
+public:
+ MyFaces(GridView const &gridView, MyGeometry const &myGeometry)
+ : lower(gridView), upper(gridView) {
+ lower.insertFacesByProperty([&](typename GridView::Intersection const &in) {
+ return isClose(myGeometry.A[1], in.geometry().center()[1]);
+ });
+
+ upper.insertFacesByProperty([&](typename GridView::Intersection const &in) {
+ return isClose(myGeometry.C[1], in.geometry().center()[1]);
+ });
+ }
+
+ BoundaryPatch<GridView> lower;
+ BoundaryPatch<GridView> upper;
+};
diff --git a/src/one-body-sample.cc b/src/one-body-sample.cc
index ef903594..2f21840c 100644
--- a/src/one-body-sample.cc
+++ b/src/one-body-sample.cc
@@ -72,6 +72,8 @@
#include "enum_verbosity.cc"
#include "enums.hh"
#include "friction_writer.hh"
+#include "mygeometry.hh"
+#include "mygrid.hh"
#include "solverfactory.hh"
#include "state.hh"
#include "timestepping.hh"
@@ -93,6 +95,8 @@ int main(int argc, char *argv[]) {
parset);
Dune::ParameterTreeParser::readOptions(argc, argv, parset);
+ MyGeometry const myGeometry;
+
auto const youngModulus = parset.get<double>("body.youngModulus");
auto const poissonRatio = parset.get<double>("body.poissonRatio");
auto const shearViscosity = parset.get<double>("body.shearViscosity");
@@ -102,23 +106,7 @@ int main(int argc, char *argv[]) {
// {{{ Set up grid
using Grid = Dune::ALUGrid<dims, dims, Dune::simplex, Dune::nonconforming>;
-
- Dune::FieldVector<size_t, dims> lowerLeft(0);
- Dune::FieldVector<size_t, dims> upperRight(1);
- upperRight[0] = 5;
- upperRight[1] = 1;
-
- std::shared_ptr<Grid> grid;
- {
- Dune::array<unsigned int, dims> elements;
- for (size_t i = 0; i < dims; ++i)
- elements[i] = upperRight[i]; // 1x1 elements
-
- grid = Dune::StructuredGridFactory<Grid>::createSimplexGrid(
- lowerLeft, upperRight, elements);
- }
- Dune::FieldVector<double, dims> zenith(0);
- zenith[1] = 1;
+ auto grid = constructGrid<Grid>(myGeometry); // FIXME
auto const refinements = parset.get<size_t>("grid.refinements");
grid->globalRefine(refinements);
@@ -128,28 +116,14 @@ int main(int argc, char *argv[]) {
GridView const leafView = grid->leafView();
// }}}
- // Set up faces
- BoundaryPatch<GridView> lowerFace(leafView);
- BoundaryPatch<GridView> upperFace(leafView);
- {
- auto const isClose = [](double a,
- double b) { return std::abs(a - b) < 1e-14; };
- lowerFace.insertFacesByProperty([&](
- typename GridView::Intersection const &in) {
- return isClose(lowerLeft[1], in.geometry().center()[1]);
- });
-
- upperFace.insertFacesByProperty([&](
- typename GridView::Intersection const &in) {
- return isClose(upperRight[1], in.geometry().center()[1]);
- });
- }
+ // Set up myFaces
+ MyFaces<GridView, MyGeometry> myFaces(leafView, myGeometry);
// Neumann boundary
BoundaryPatch<GridView> const neumannBoundary(leafView);
// Frictional Boundary
- BoundaryPatch<GridView> const &frictionalBoundary = lowerFace;
+ BoundaryPatch<GridView> const &frictionalBoundary = myFaces.lower;
Dune::BitSetVector<1> frictionalNodes(fineVertexCount);
frictionalBoundary.getVertices(frictionalNodes);
@@ -157,10 +131,10 @@ int main(int argc, char *argv[]) {
Dune::BitSetVector<dims> noNodes(fineVertexCount);
Dune::BitSetVector<dims> dirichletNodes(fineVertexCount);
for (size_t i = 0; i < fineVertexCount; ++i) {
- if (lowerFace.containsVertex(i))
+ if (myFaces.lower.containsVertex(i))
dirichletNodes[i][1] = true;
- if (upperFace.containsVertex(i))
+ if (myFaces.upper.containsVertex(i))
dirichletNodes[i] = true;
}
@@ -204,7 +178,8 @@ int main(int argc, char *argv[]) {
// Assemble forces
Vector gravityFunctional;
- myAssembler.assembleBodyForce(gravity, density, zenith, gravityFunctional);
+ myAssembler.assembleBodyForce(gravity, density, myGeometry.zenith,
+ gravityFunctional);
// Problem formulation: right-hand side
auto const computeExternalForces = [&](double _relativeTime, Vector &_ell) {
@@ -219,12 +194,12 @@ int main(int argc, char *argv[]) {
{
double volume = 1.0;
for (size_t i = 0; i < dims; ++i)
- volume *= (upperRight[i] - lowerLeft[i]);
+ volume *= (myGeometry.C[i] - myGeometry.A[i]);
double area = 1.0;
for (size_t i = 0; i < dims; ++i)
if (i != 1)
- area *= (upperRight[i] - lowerLeft[i]);
+ area *= (myGeometry.C[i] - myGeometry.A[i]);
// volume * gravity * density / area = normal stress
// V * g * rho / A = sigma_n
--
GitLab