diff --git a/dune/gfe/cosseratvtkwriter.hh b/dune/gfe/cosseratvtkwriter.hh
index 4fe9fffe1b0d05f7b1cf2b5f2e0290b86f747a0a..69eb2dfbac78c7a5f3d76d0be9fbcb85de243b0c 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 856287f544a5c79dc634d4493b0345b870eb99c3..f656ae22e382ada3b830110d66e850be3571be60 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 6dcd7930c8c9fd86238e1e16b133cd20cde15fc3..c7f937da5bd087ab9c0aa3363a8f94298fb60778 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 350a544e0d59661814c5043a4b4435e7fd44a0e5..e44b1c24c8880c535065f4c1658b24b42c2b70d0 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");