From 699c6b19cda0b6da8d6e5fff242bc7ae2f39d2b5 Mon Sep 17 00:00:00 2001
From: Oliver Sander <oliver.sander@tu-dresden.de>
Date: Sun, 13 Oct 2024 10:39:09 +0200
Subject: [PATCH] rod3d: Get the Dirichlet boundary from the Python file
This is more flexible, and more consistent with what
cosserat-continuum.cc does.
---
problems/cosserat-rod-example.py | 19 +++++++++++++++----
src/rod3d.cc | 30 ++++++++++++++++++++----------
2 files changed, 35 insertions(+), 14 deletions(-)
diff --git a/problems/cosserat-rod-example.py b/problems/cosserat-rod-example.py
index d58d65fd7..477a9362b 100644
--- a/problems/cosserat-rod-example.py
+++ b/problems/cosserat-rod-example.py
@@ -5,6 +5,9 @@ class ParameterSet(dict):
parameterSet = ParameterSet()
+# Interpolation method
+parameterSet.interpolationMethod = "geodesic"
+
#############################################
# Grid parameters
#############################################
@@ -56,19 +59,27 @@ parameterSet.baseTolerance = 1e-8
parameterSet.instrumented = "no"
+
############################
-# Problem specifications
+# Material parameters
############################
-# Interpolation method
-parameterSet.interpolationMethod = "geodesic"
-
parameterSet.A = 1
parameterSet.J1 = 1
parameterSet.J2 = 1
parameterSet.E = 2.5e5
parameterSet.nu = 0.3
+#############################################
+# Boundary values
+#############################################
+
+### Python predicate specifying all Dirichlet grid vertices
+# x is the vertex coordinate
+parameterSet.dirichletVerticesPredicate = "[x[2] < 0.001 or x[2] > 0.999, x[2] < 0.001 or x[2] > 0.999, x[2] < 0.001 or x[2] > 0.999]"
+parameterSet.dirichletRotationVerticesPredicate = "x[2] < 0.001 or x[2] > 0.999"
+
+
parameterSet.dirichletValue = "1 0 0"
parameterSet.dirichletAxis = "1 0 0"
diff --git a/src/rod3d.cc b/src/rod3d.cc
index f28f30053..786c48138 100644
--- a/src/rod3d.cc
+++ b/src/rod3d.cc
@@ -27,8 +27,6 @@
#include <dune/vtk/vtkwriter.hh>
#include <dune/vtk/datacollectors/lagrangedatacollector.hh>
-#include <dune/fufem/boundarypatch.hh>
-#include <dune/fufem/functiontools/boundarydofs.hh>
#include <dune/fufem/dunepython.hh>
#include <dune/gfe/assemblers/cosseratrodenergy.hh>
@@ -154,15 +152,27 @@ int main (int argc, char *argv[]) try
blockedInterleaved()
));
- // Find all boundary dofs
- BoundaryPatch<GridView> dirichletBoundary(gridView,
- true); // true: The entire boundary is Dirichlet boundary
BitSetVector<TargetSpace::TangentVector::dimension> dirichletNodes(tangentBasis.size(), false);
-#if DUNE_VERSION_GTE(DUNE_FUFEM, 2, 10)
- Fufem::markBoundaryPatchDofs(dirichletBoundary,tangentBasis,dirichletNodes);
-#else
- constructBoundaryDofs(dirichletBoundary,tangentBasis,dirichletNodes);
-#endif
+
+ // Make Python function that computes which vertices are on the Dirichlet boundary, based on the vertex positions.
+ std::string lambda = std::string("lambda x: (") + parameterSet.get<std::string>("dirichletVerticesPredicate") + std::string(")");
+ auto pythonDirichletVertices = Python::make_function<FieldVector<bool,3> >(Python::evaluate(lambda));
+
+ lambda = std::string("lambda x: (") + parameterSet.get<std::string>("dirichletRotationVerticesPredicate") + std::string(")");
+ auto pythonOrientationDirichletVertices = Python::make_function<bool>(Python::evaluate(lambda));
+
+ for (size_t i=0; i<tangentBasis.size(); i++)
+ {
+ FieldVector<bool,3> isDirichlet = pythonDirichletVertices(referenceConfiguration[i][_0].globalCoordinates());
+ for (size_t j=0; j<3; j++)
+ dirichletNodes[i][j] = isDirichlet[j];
+
+ bool isDirichletOrientation = pythonOrientationDirichletVertices(referenceConfiguration[i][_0].globalCoordinates());
+ for (size_t j=0; j<3; j++)
+ dirichletNodes[i][j+3] = isDirichletOrientation;
+ }
+
+ std::cout << "Dirichlet boundary has " << dirichletNodes.count() << " degrees of freedom.\n";
// Find the dof on the right boundary
std::size_t rightBoundaryDof;
--
GitLab