From 8633a47c85b3810db1824fa07f2379f44d86907f Mon Sep 17 00:00:00 2001
From: Oliver Sander <oliver.sander@tu-dresden.de>
Date: Mon, 11 Nov 2024 13:49:41 +0100
Subject: [PATCH] Simplify interface of CosseratVTKWriter
Previously, the writer requested an FE basis to represent the director
vectors. With the help of ComposedGridViewFunction this patch allows
to avoid this extra basis, which makes for a much nicer interface
for CosseratVTKWriter.
---
dune/gfe/cosseratvtkwriter.hh | 47 ++++++++++++++---------------
src/cosserat-continuum.cc | 22 +++++++-------
src/cosserat-rod.cc | 16 +++++-----
test/nonplanarcosseratenergytest.cc | 22 ++++++--------
4 files changed, 51 insertions(+), 56 deletions(-)
diff --git a/dune/gfe/cosseratvtkwriter.hh b/dune/gfe/cosseratvtkwriter.hh
index 4fe9fffe..69eb2dfb 100644
--- a/dune/gfe/cosseratvtkwriter.hh
+++ b/dune/gfe/cosseratvtkwriter.hh
@@ -11,6 +11,7 @@
#include <dune/functions/functionspacebases/lagrangebasis.hh>
#include <dune/functions/functionspacebases/interpolate.hh>
#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
+#include <dune/functions/gridfunctions/composedgridfunction.hh>
#include <dune/vtk/vtkwriter.hh>
#include <dune/vtk/datacollectors/lagrangedatacollector.hh>
@@ -395,11 +396,10 @@ public:
* \param directorBasis The basis that will be used to represent director fields
* \param order The polynomial order that the data will have in the VTK file
*/
- template <typename DisplacementFunction, typename DirectorBasis>
+ template <typename DisplacementFunction, typename OrientationFunction>
static void write(const GridView& gridView,
const DisplacementFunction& displacement,
- const DirectorBasis& directorBasis,
- const std::vector<Rotation<double,3> >& orientationConfiguration,
+ const OrientationFunction& orientation,
int order,
const std::string& filename)
{
@@ -417,27 +417,26 @@ public:
vtkWriter.addPointData(displacement, Dune::VTK::FieldInfo("displacement", Dune::VTK::FieldInfo::Type::vector, 3));
// Attach the director fields
- BlockVector<FieldVector<double,3> > director0(orientationConfiguration.size());
- BlockVector<FieldVector<double,3> > director1(orientationConfiguration.size());
- BlockVector<FieldVector<double,3> > director2(orientationConfiguration.size());
-
- for (size_t i=0; i<orientationConfiguration.size(); i++)
- {
- FieldMatrix<double,3,3> rotationMatrix;
- orientationConfiguration[i].matrix(rotationMatrix);
-
- // Extract the directors, i.e., the columns of the matrix
- for (size_t j=0; j<3; j++)
- {
- director0[i][j] = rotationMatrix[j][0];
- director1[i][j] = rotationMatrix[j][1];
- director2[i][j] = rotationMatrix[j][2];
- }
- }
-
- auto director0Function = Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,3> >(directorBasis, director0);
- auto director1Function = Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,3> >(directorBasis, director1);
- auto director2Function = Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,3> >(directorBasis, director2);
+ // This lambda takes a unit quaternion and extracts one column
+ // of the corresponding rotation matrix
+ auto directorExtractor = [](FieldVector<double,4> q,int columnNumber) -> FieldVector<double,3>
+ {
+ FieldMatrix<double,3,3> matrix;
+ Rotation<double,3>(q).matrix(matrix);
+ FieldVector<double,3> column;
+ for (size_t i=0; i<3; ++i)
+ column[i] = matrix[i][columnNumber];
+ return column;
+ };
+
+ auto director0Function = Functions::makeComposedGridFunction(std::bind(directorExtractor,std::placeholders::_1,0),
+ orientation);
+
+ auto director1Function = Functions::makeComposedGridFunction(std::bind(directorExtractor,std::placeholders::_1,1),
+ orientation);
+
+ auto director2Function = Functions::makeComposedGridFunction(std::bind(directorExtractor,std::placeholders::_1,2),
+ orientation);
vtkWriter.addPointData(director0Function, Dune::VTK::FieldInfo("director0", Dune::VTK::FieldInfo::Type::vector, 3));
vtkWriter.addPointData(director1Function, Dune::VTK::FieldInfo("director1", Dune::VTK::FieldInfo::Type::vector, 3));
diff --git a/src/cosserat-continuum.cc b/src/cosserat-continuum.cc
index 856287f5..f656ae22 100644
--- a/src/cosserat-continuum.cc
+++ b/src/cosserat-continuum.cc
@@ -287,13 +287,6 @@ int main (int argc, char *argv[]) try
lagrange<displacementOrder>()
));
- // Vector-valued basis to represent directors
- auto directorBasis = makeBasis(
- gridView,
- power<3>(
- lagrange<rotationOrder>()
- ));
-
// Matrix-valued basis for treating the microrotation as a matrix field
auto orientationMatrixBasis = makeBasis(
gridView,
@@ -483,11 +476,19 @@ int main (int argc, char *argv[]) try
auto displacementFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,3> >(deformationPowerBasis, displacement);
+#ifdef PROJECTED_INTERPOLATION
+ using RotationInterpolationRule = GFE::LocalProjectedFEFunction<dim, double, OrientationFEBasis::LocalView::Tree::FiniteElement, Rotation<double,3> >;
+#else
+ using RotationInterpolationRule = LocalGeodesicFEFunction<dim, double, OrientationFEBasis::LocalView::Tree::FiniteElement, Rotation<double,3> >;
+#endif
+
+ GFE::EmbeddedGlobalGFEFunction<OrientationFEBasis, RotationInterpolationRule,Rotation<double,3> > orientationFunction(orientationFEBasis,
+ x[_1]);
+
if (dim == dimworld) {
CosseratVTKWriter<GridView>::write(gridView,
displacementFunction,
- directorBasis,
- x[_1],
+ orientationFunction,
std::max(LFE_ORDER, GFE_ORDER),
resultPath + "cosserat_homotopy_0_l" + std::to_string(numLevels));
} else if (dim == 2 && dimworld == 3) {
@@ -753,8 +754,7 @@ int main (int argc, char *argv[]) try
if (dim == dimworld) {
CosseratVTKWriter<GridView>::write(gridView,
displacementFunction,
- directorBasis,
- x[_1],
+ orientationFunction,
std::max(LFE_ORDER, GFE_ORDER),
resultPath + "cosserat_homotopy_" + std::to_string(i+1) + "_l" + std::to_string(numLevels));
} else if (dim == 2 && dimworld == 3) {
diff --git a/src/cosserat-rod.cc b/src/cosserat-rod.cc
index 6dcd7930..c7f937da 100644
--- a/src/cosserat-rod.cc
+++ b/src/cosserat-rod.cc
@@ -31,6 +31,7 @@
#include <dune/gfe/assemblers/geodesicfeassembler.hh>
#include <dune/gfe/assemblers/localgeodesicfeadolcstiffness.hh>
#include <dune/gfe/cosseratvtkwriter.hh>
+#include <dune/gfe/embeddedglobalgfefunction.hh>
#include <dune/gfe/localgeodesicfefunction.hh>
#include <dune/gfe/localprojectedfefunction.hh>
#include <dune/gfe/riemanniantrsolver.hh>
@@ -288,13 +289,6 @@ int main (int argc, char *argv[]) try
power<3>(lagrange<order>())
);
- // Vector-valued basis to represent directors
- auto directorBasis = makeBasis(
- gridView,
- power<3>(
- lagrange<order>()
- ));
-
// Compute the displacement from the deformation, because that's more easily visualized
// in ParaView
BlockVector<FieldVector<double,3> > displacement(worldBasis.size());
@@ -319,10 +313,14 @@ int main (int argc, char *argv[]) try
for (size_t i=0; i<x.size(); ++i)
orientationConfiguration[i] = x[i][_1];
+ using RotationInterpolationRule = LocalGeodesicFEFunction<1, double, ScalarBasis::LocalView::Tree::FiniteElement, Rotation<double,3> >;
+
+ GFE::EmbeddedGlobalGFEFunction<ScalarBasis, RotationInterpolationRule,Rotation<double,3> > orientationFunction(scalarBasis,
+ orientationConfiguration);
+
CosseratVTKWriter<GridView>::write(gridView,
displacementFunction,
- directorBasis,
- orientationConfiguration,
+ orientationFunction,
order,
resultPath + "cosserat-rod-result");
diff --git a/test/nonplanarcosseratenergytest.cc b/test/nonplanarcosseratenergytest.cc
index 350a544e..e44b1c24 100644
--- a/test/nonplanarcosseratenergytest.cc
+++ b/test/nonplanarcosseratenergytest.cc
@@ -81,11 +81,11 @@ double calculateEnergy(const FlatGridView& flatGridView,
// Turn matrix-valued function into quaternion-valued function
auto orientationQuaternionFunction
= [&orientationFunction](FieldVector<double,3> x) -> FieldVector<double,4>
- {
- Rotation<double,3> rotation;
- rotation.set(orientationFunction(x));
- return rotation;
- };
+ {
+ Rotation<double,3> rotation;
+ rotation.set(orientationFunction(x));
+ return rotation;
+ };
BlockVector<FieldVector<double,4> > orientationAsVector(flatFEBasis.size());
Functions::interpolate(curvedGridQuaternionBasis, orientationAsVector, orientationQuaternionFunction);
@@ -96,17 +96,15 @@ double calculateEnergy(const FlatGridView& flatGridView,
// Write the configuration to a file (just for debugging)
/////////////////////////////////////////////////////////////////////////
- auto directorBasis = makeBasis(
- flatGridView,
- power<3>(
- lagrange<2>()
- ));
+ // The orientation function needs to become a GridViewFunction,
+ // otherwise the current CosseratVTKWriter will not accept it.
+ auto orientationQuaternionGridViewFunction = Functions::makeAnalyticGridViewFunction(orientationQuaternionFunction, flatGridView);
+
// TODO: Write the curved grid, not the flat one
CosseratVTKWriter<FlatGridView>::write(flatGridView,
curvedGridGeometry,
- directorBasis,
- configuration[_1],
+ orientationQuaternionGridViewFunction,
2, // VTK output element order
"nonplanarcosseratenergytest-result.vtu");
--
GitLab