From 1409a3f07d92f022ea68f3ecbe44df4dd58ccc27 Mon Sep 17 00:00:00 2001
From: Oliver Sander <oliver.sander@tu-dresden.de>
Date: Wed, 21 Aug 2024 10:01:36 +0200
Subject: [PATCH] Use the DirichletValues Python class also for the initial
iterate
That way, we can also set the initial microrotation -- previously
it was hard-wired to the identity matrix.
Being able to explicitly choose the initial microrotation allows
to start, in particular, from a twisted strip.
---
problems/cosserat-continuum-cantilever.py | 4 +--
problems/cosserat-continuum-twisted-strip.py | 4 +--
...rat-continuum-wong-pellegrino-wrinkling.py | 8 +++--
.../cosserat-continuum-wriggers-l-shape.py | 2 ++
src/cosserat-continuum.cc | 30 ++++++++++++++-----
5 files changed, 32 insertions(+), 16 deletions(-)
diff --git a/problems/cosserat-continuum-cantilever.py b/problems/cosserat-continuum-cantilever.py
index f99abbfe..0dfec2fe 100644
--- a/problems/cosserat-continuum-cantilever.py
+++ b/problems/cosserat-continuum-cantilever.py
@@ -112,6 +112,7 @@ parameterSet.dirichletVerticesPredicate = "[x[0] < 0.01, x[0] < 0.01, x[0] < 0.0
parameterSet.dirichletRotationVerticesPredicate = "x[0] < 0.01"
### The actual Dirichlet values
+# With 'homotopyParameter==0', this class gives the initial iterate
class DirichletValues:
def __init__(self, homotopyParameter):
self.homotopyParameter = homotopyParameter
@@ -134,8 +135,5 @@ parameterSet.neumannVerticesPredicate = "x[0] > 99.99"
### Neumann values
parameterSet.neumannValues = "0 0 3"
-# Initial deformation
-parameterSet.initialDeformation = "[x[0], x[1], 0]"
-
#parameterSet.startFromFile = yes
#parameterSet.initialIterateFilename = initial_iterate.vtu
diff --git a/problems/cosserat-continuum-twisted-strip.py b/problems/cosserat-continuum-twisted-strip.py
index d32cf33e..8439f111 100644
--- a/problems/cosserat-continuum-twisted-strip.py
+++ b/problems/cosserat-continuum-twisted-strip.py
@@ -109,6 +109,7 @@ parameterSet.dirichletVerticesPredicate = "[x[0] < 0.001 or x[0] > 0.0999, x[0]
parameterSet.dirichletRotationVerticesPredicate = "x[0] < 0.001 or x[0] > 0.0999"
### The actual Dirichlet values
+# With 'homotopyParameter==0', this class gives the initial iterate
class DirichletValues:
def __init__(self, homotopyParameter):
self.homotopyParameter = homotopyParameter
@@ -135,7 +136,4 @@ class DirichletValues:
rotation = [[1,0,0], [0, math.cos(angle), -math.sin(angle)], [0, math.sin(angle), math.cos(angle)]]
return rotation
-# Initial deformation
-parameterSet.initialDeformation = "[x[0], x[1], 0]"
-
diff --git a/problems/cosserat-continuum-wong-pellegrino-wrinkling.py b/problems/cosserat-continuum-wong-pellegrino-wrinkling.py
index 7b805648..8bce680b 100644
--- a/problems/cosserat-continuum-wong-pellegrino-wrinkling.py
+++ b/problems/cosserat-continuum-wong-pellegrino-wrinkling.py
@@ -108,12 +108,17 @@ parameterSet.dirichletVerticesPredicate = "[x[1] < 0.0001 or x[1] > 0.128 - 0.00
parameterSet.dirichletRotationVerticesPredicate = "x[1] < 0.0001 or x[1] > 0.128 - 0.0001"
### The actual Dirichlet values
+# With 'homotopyParameter==0', this class gives the initial iterate
class DirichletValues:
def __init__(self, homotopyParameter):
self.homotopyParameter = homotopyParameter
def deformation(self, x):
- out = [x[0], x[1], 0]
+ if self.homotopyParameter<0.01:
+ out = [x[0] + 0.003*x[1] / 0.128, x[1], 0.002*math.cos(1e4*x[0])]
+ else:
+ out = [x[0], x[1], 0]
+
if x[1] > 0.128-1e-4 :
out[0] += 0.003 * self.homotopyParameter
out[1] += 0.0005
@@ -126,4 +131,3 @@ class DirichletValues:
# Initial deformation
#parameterSet.startFromFile = True
parameterSet.initialIterateFilename = "cosserat_iterate_2.vtu"
-parameterSet.initialDeformation = "[x[0] + 0.003*x[1] / 0.128, x[1], 0.002*math.cos(1e4*x[0])]"
diff --git a/problems/cosserat-continuum-wriggers-l-shape.py b/problems/cosserat-continuum-wriggers-l-shape.py
index cecbc1ae..edd8a706 100644
--- a/problems/cosserat-continuum-wriggers-l-shape.py
+++ b/problems/cosserat-continuum-wriggers-l-shape.py
@@ -104,6 +104,7 @@ parameterSet.dirichletVerticesPredicate = "[x[0] < 1, x[0] < 1, x[0] < 1]"
parameterSet.dirichletRotationVerticesPredicate = "x[0] < 1"
### The actual Dirichlet values
+# With 'homotopyParameter==0', this class gives the initial iterate
class DirichletValues:
def __init__(self, homotopyParameter):
self.homotopyParameter = homotopyParameter
@@ -134,6 +135,7 @@ parameterSet.neumannValues = "0.09 0 0"
#parameterSet.initialDeformation = "[x[0], x[1], 0 if (x[0] < 225 or x[1] < -15) else 0.001*(x[0]-225)*(x[1]+15)]"
##parameterSet.initialDeformation = "[x[0], x[1], 0]"
+# Overrides the configuration given by the 'DirichletValues' class
parameterSet.startFromFile = True
parameterSet.initialIterateGridFilename = "wriggers-L-shape_99_mm.msh"
parameterSet.initialIterateFilename = "initial-wriggers-l-shape.vtu"
diff --git a/src/cosserat-continuum.cc b/src/cosserat-continuum.cc
index a94d0a66..cdfe2623 100644
--- a/src/cosserat-continuum.cc
+++ b/src/cosserat-continuum.cc
@@ -356,19 +356,33 @@ int main (int argc, char *argv[]) try
// //////////////////////////
SolutionType x;
-
x[_0].resize(compositeBasis.size({0}));
+ x[_1].resize(compositeBasis.size({1}));
- lambda = std::string("lambda x: (") + parameterSet.get<std::string>("initialDeformation") + std::string(")");
- auto pythonInitialDeformation = Python::make_function<FieldVector<double,3> >(Python::evaluate(lambda));
+ // Load the initial configuration from the Python options file
+ // Yes, the Python class really is called 'DirichletValues'.
+ // We use that class both for the initial iterate and the Dirichlet values
+ Python::Callable initialConfigurationPythonClass = pyModule.get("DirichletValues");
- std::vector<FieldVector<double,3> > v;
- Functions::interpolate(deformationPowerBasis, v, pythonInitialDeformation);
+ // Construct an object
+ Python::Reference initialConfigurationPythonObject = initialConfigurationPythonClass(0.0 /* homotopy parameter*/);
- for (size_t i=0; i<x[_0].size(); i++)
- x[_0][i] = v[i];
+ // Extract object member functions as Dune functions
+ auto initialDeformationFunction = Python::make_function<FieldVector<double,3> > (initialConfigurationPythonObject.get("deformation"));
+ auto initialOrientationFunction = Python::make_function<FieldMatrix<double,3,3> > (initialConfigurationPythonObject.get("orientation"));
+
+ BlockVector<FieldVector<double,3> > ddV;
+ Functions::interpolate(deformationPowerBasis, ddV, initialDeformationFunction);
+
+ BlockVector<FieldMatrix<double,3,3> > dOV;
+ Functions::interpolate(orientationMatrixBasis, dOV, initialOrientationFunction);
+
+ for (std::size_t i = 0; i < compositeBasis.size({0}); i++)
+ x[_0][i] = ddV[i];
+
+ for (std::size_t i = 0; i < compositeBasis.size({1}); i++)
+ x[_1][i].set(dOV[i]);
- x[_1].resize(compositeBasis.size({1}));
#if !MIXED_SPACE
if (parameterSet.hasKey("startFromFile"))
{
--
GitLab