From 5da22230c29d3576494abf96c7c00dd800f4fb19 Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Wed, 28 May 2014 09:39:39 +0000
Subject: [PATCH] Start treating the twisted-strip example again
Make its DirichletValues function a proper Dune VirtualFunction.
Prescribe Dirichlet values for the transverse director, as Neff
does it in his papers.
[[Imported from SVN: r9773]]
---
cosserat-continuum.cc | 90 ++++++++++++++++++++++++++++++++-----------
1 file changed, 67 insertions(+), 23 deletions(-)
diff --git a/cosserat-continuum.cc b/cosserat-continuum.cc
index 44752746..018c4437 100644
--- a/cosserat-continuum.cc
+++ b/cosserat-continuum.cc
@@ -65,15 +65,15 @@ public:
};
-#if 1
// Dirichlet boundary data for the shear/wrinkling example
-class DirichletValues
+class WrinklingDirichletValues
: public Dune::VirtualFunction<FieldVector<double,dim>, FieldVector<double,3> >
{
+ FieldVector<double,2> upper_;
double homotopy_;
public:
- DirichletValues(double homotopy)
- : homotopy_(homotopy)
+ WrinklingDirichletValues(FieldVector<double,2> upper, double homotopy)
+ : upper_(upper), homotopy_(homotopy)
{}
void evaluate(const FieldVector<double,dim>& in, FieldVector<double,3>& out) const
@@ -83,22 +83,28 @@ public:
out[i] = in[i];
// if (out[1] > 1-1e-3)
- if (out[1] > 0.128-1e-4)
+ if (out[1] > upper_[1]-1e-4)
out[0] += 0.003*homotopy_;
}
};
-#endif
-#if 0
+
+
// Dirichlet boundary data for the 'twisted-strip' example
-void dirichletValues(const FieldVector<double,dim>& in, FieldVector<double,3>& out,
- double homotopy
-)
+class TwistedStripDeformationDirichletValues
+: public Dune::VirtualFunction<FieldVector<double,dim>, FieldVector<double,3> >
{
- if (not vIt->geometry().corner(0)[0] > upper[0]-1e-3)
- return;
+ FieldVector<double,2> upper_;
+ double homotopy_;
+public:
+
+ TwistedStripDeformationDirichletValues(FieldVector<double,2> upper, double homotopy)
+ : upper_(upper), homotopy_(homotopy)
+ {}
- double angle = 8*M_PI;
- angle *= homotopy;
+ void evaluate(const FieldVector<double,dim>& in, FieldVector<double,3>& out) const
+ {
+ double angle = 0.5*M_PI * in[0]/upper_[0];
+ angle *= homotopy_;
// center of rotation
FieldVector<double,3> center(0);
@@ -119,8 +125,38 @@ void dirichletValues(const FieldVector<double,dim>& in, FieldVector<double,3>& o
rotation.mv(inEmbedded, out);
out += center;
-}
-#endif
+ }
+};
+
+// Dirichlet boundary data for the 'twisted-strip' example
+class TwistedStripOrientationDirichletValues
+: public Dune::VirtualFunction<FieldVector<double,dim>, FieldMatrix<double,3,3> >
+{
+ FieldVector<double,2> upper_;
+ double homotopy_;
+public:
+
+ TwistedStripOrientationDirichletValues(FieldVector<double,2> upper, double homotopy)
+ : upper_(upper), homotopy_(homotopy)
+ {}
+
+ void evaluate(const FieldVector<double,dim>& in, FieldMatrix<double,3,3>& out) const
+ {
+ double angle = 0.5*M_PI * in[0]/upper_[0];
+ angle *= homotopy_;
+
+ // center of rotation
+ FieldMatrix<double,3,3> rotation(0);
+ rotation[0][0] = 1;
+ rotation[1][1] = std::cos(angle);
+ rotation[1][2] = -std::sin(angle);
+ rotation[2][1] = std::sin(angle);
+ rotation[2][2] = std::cos(angle);
+
+ out = rotation;
+ }
+};
+
/** \brief A constant vector-valued function, for simple Neumann boundary values */
struct NeumannFunction
@@ -228,11 +264,11 @@ int main (int argc, char *argv[]) try
if (vIt->geometry().corner(0)[0] > upper[0]-1e-3 )
neumannNodes[indexSet.index(*vIt)] = true;
#endif
-#if 1 // Boundary conditions for the shearing/wrinkling example
+#if 0 // Boundary conditions for the shearing/wrinkling example
if (vIt->geometry().corner(0)[1] < 1e-4 or vIt->geometry().corner(0)[1] > upper[1]-1e-4 )
dirichletVertices[indexSet.index(*vIt)] = true;
#endif
-#if 0 // Boundary conditions for the twisted-strip example
+#if 1 // Boundary conditions for the twisted-strip example
if (vIt->geometry().corner(0)[0] < lower[0]+1e-3 or vIt->geometry().corner(0)[0] > upper[0]-1e-3 )
dirichletVertices[indexSet.index(*vIt)] = true;
#endif
@@ -333,13 +369,21 @@ int main (int argc, char *argv[]) try
// Set Dirichlet values
////////////////////////////////////////////////////////
- DirichletValues dirichletValues(homotopyParameter);
- std::vector<FieldVector<double,3> > dV;
- Functions::interpolate(feBasis, dV, dirichletValues, dirichletDofs);
+ TwistedStripDeformationDirichletValues deformationDirichletValues(upper,homotopyParameter);
+ std::vector<FieldVector<double,3> > ddV;
+ Functions::interpolate(feBasis, ddV, deformationDirichletValues, dirichletDofs);
+
+
+ TwistedStripOrientationDirichletValues orientationDirichletValues(upper,homotopyParameter);
+ std::vector<FieldMatrix<double,3,3> > dOV;
+ Functions::interpolate(feBasis, dOV, orientationDirichletValues, dirichletDofs);
for (size_t j=0; j<x.size(); j++)
- if (dirichletDofs[j][0])
- x[j].r = dV[j];
+ if (dirichletNodes[j][0])
+ {
+ x[j].r = ddV[j];
+ x[j].q.set(dOV[j]);
+ }
// /////////////////////////////////////////////////////
// Solve!
--
GitLab