Skip to content
Snippets Groups Projects
Commit 8633a47c authored by Sander, Oliver's avatar Sander, Oliver
Browse files

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.
parent d7bdb0f0
Branches
No related tags found
1 merge request!167Simplify interface of CosseratVTKWriter
Hide whitespace changes
Inline Side-by-side
dune/gfe/cosseratvtkwriter.hh
+ 23
− 24
View file @ 8633a47c
......@@ -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));
......
src/cosserat-continuum.cc
+ 11
− 11
View file @ 8633a47c
......@@ -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) {
......
src/cosserat-rod.cc
+ 7
− 9
View file @ 8633a47c
......@@ -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");
......
test/nonplanarcosseratenergytest.cc
+ 10
− 12
View file @ 8633a47c
......@@ -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");
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment