From 2eeed182c901809e802fe8c031861601b8ffe848 Mon Sep 17 00:00:00 2001 From: Lisa Julia Nebel <lisa_julia.nebel@tu-dresden.de> Date: Fri, 2 Jul 2021 08:09:42 +0200 Subject: [PATCH 1/2] Add nonplanarcosseratshellenergytest Add a simple test with a grid containing one element testing if the energy is invariant of the number of grid refinements. --- test/CMakeLists.txt | 1 + test/nonplanarcosseratenergytest.cc | 193 ++++++++++++++++++++++++++++ 2 files changed, 194 insertions(+) create mode 100644 test/nonplanarcosseratenergytest.cc diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 75384afe4..b62fb7a9d 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -8,6 +8,7 @@ set(TESTS localgeodesicfefunctiontest localgfetestfunctiontest localprojectedfefunctiontest + nonplanarcosseratenergytest orthogonalmatrixtest rigidbodymotiontest rotationtest diff --git a/test/nonplanarcosseratenergytest.cc b/test/nonplanarcosseratenergytest.cc new file mode 100644 index 000000000..68b055697 --- /dev/null +++ b/test/nonplanarcosseratenergytest.cc @@ -0,0 +1,193 @@ +#include "config.h" + +#include <dune/foamgrid/foamgrid.hh> + +#include <dune/geometry/type.hh> +#include <dune/geometry/quadraturerules.hh> + +#include <dune/functions/functionspacebases/interpolate.hh> +#include <dune/functions/functionspacebases/lagrangebasis.hh> +#include <dune/functions/functionspacebases/powerbasis.hh> + +#include <dune/gfe/cosseratvtkwriter.hh> +#include <dune/gfe/nonplanarcosseratshellenergy.hh> +#include <dune/gfe/rigidbodymotion.hh> + +#include "multiindex.hh" +#include "valuefactory.hh" + +using namespace Dune; + +static const int dim = 2; +static const int dimworld = 3; + +using GridType = FoamGrid<dim,dimworld>; +using TargetSpace = RigidBodyMotion<double,dimworld>; + +////////////////////////////////////////////////////////// +// Make a test grid consisting of a single triangle +////////////////////////////////////////////////////////// + +template <class GridType> +std::unique_ptr<GridType> makeSingleElementGrid() +{ + constexpr auto triangle = Dune::GeometryTypes::triangle; + GridFactory<GridType> factory; + + FieldVector<double,dimworld> vertex0{0,0,0}; + FieldVector<double,dimworld> vertex1{0,1,0}; + FieldVector<double,dimworld> vertex2{1,0,0}; + factory.insertVertex(vertex0); + factory.insertVertex(vertex1); + factory.insertVertex(vertex2); + + factory.insertElement(triangle, {0,1,2}); + + return std::unique_ptr<GridType>(factory.createGrid()); +} + + +////////////////////////////////////////////////////////////////////////////////////// +// Test energy computation for the same grid with different refinement levels +////////////////////////////////////////////////////////////////////////////////////// + +TargetSpace getConfiguration(const FieldVector<double, dimworld>& point) { + FieldVector<double, dimworld> displacementAtPoint(0); + FieldVector<double,4> rotationVectorAtPoint(0); + + if (point[0] == 0 and point[1] == 0 and point[2] == 0) { + //0 0 0 + displacementAtPoint = {0, 0, 0}; + rotationVectorAtPoint = {0, 0, 0, 1}; + } else if (point[0] == 1 and point[1] == 0 and point[2] == 0) { + //1 0 0 + displacementAtPoint = {0, 0, 1}; + rotationVectorAtPoint = {0, 0, 0, 1}; + } else if (point[0] == 0 and point[1] == 1 and point[2] == 0) { + //0 1 0 + displacementAtPoint = {0, 0, 1}; + rotationVectorAtPoint = {0, 0, 0, 1}; + } else if (point[0] == 0.5 and point[1] == 0 and point[2] == 0) { + //0.5 0 0 + displacementAtPoint = {0, 0, 0.5}; + rotationVectorAtPoint = {0, 0, 0, 1}; + } else if (point[0] == 0 and point[1] == 0.5 and point[2] == 0) { + //0 0.5 0 + displacementAtPoint = {0, 0, 0.5}; + rotationVectorAtPoint = {0, 0, 0, 1}; + } else if (point[0] == 0.5 and point[1] == 0.5 and point[2] == 0) { + //0.5 0.5 0 + displacementAtPoint = {0, 0, 1}; + rotationVectorAtPoint = {0, 0, 0, 1}; + } + + TargetSpace configuration; + for (int i = 0; i < dimworld; i++) + configuration.r[i] = point[i] + displacementAtPoint[i]; + + Rotation<double,dimworld> rotation(rotationVectorAtPoint); + FieldMatrix<double,dimworld,dimworld> rotationMatrix(0); + rotation.matrix(rotationMatrix); + configuration.q.set(rotationMatrix); + + return configuration; +} + +double calculateEnergy(const int numLevels) +{ + ParameterTree materialParameters; + materialParameters["thickness"] = "0.1"; + materialParameters["mu"] = "3.8462e+05"; + materialParameters["lambda"] = "2.7149e+05"; + materialParameters["mu_c"] = "0"; + materialParameters["L_c"] = "1e-3"; + materialParameters["q"] = "2.5"; + materialParameters["kappa"] = "0.1"; + materialParameters["b1"] = "1"; + materialParameters["b2"] = "1"; + materialParameters["b3"] = "1"; + + const std::unique_ptr<GridType> grid = makeSingleElementGrid<GridType>(); + grid->globalRefine(numLevels-1); + GridType::LeafGridView gridView = grid->leafGridView(); + + using FEBasis = Dune::Functions::LagrangeBasis<typename GridType::LeafGridView,1>; + FEBasis feBasis(gridView); + + using namespace Dune::Functions::BasisFactory; + using namespace Dune::TypeTree::Indices; + + auto deformationPowerBasis = makeBasis( + gridView, + power<dimworld>( + lagrange<1>() + )); + + BlockVector<FieldVector<double,3> > helperVector; + Dune::Functions::interpolate(deformationPowerBasis, helperVector, [](FieldVector<double,dimworld> x){ + auto out = x; + out[2] += x[0]; + return out; + }); + //Dune::Functions::interpolate(deformationPowerBasis, helperVector, [](FieldVector<double,dimworld> x){ return x; }); + + auto stressFreeConfiguration = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dimworld>>(deformationPowerBasis, helperVector); + + NonplanarCosseratShellEnergy<FEBasis, 3, double, decltype(stressFreeConfiguration)> localCosseratEnergyPlanar(materialParameters, + &stressFreeConfiguration, + nullptr, + nullptr, + nullptr); + + BlockVector<FieldVector<double,3> > id; + Dune::Functions::interpolate(deformationPowerBasis, id, [](FieldVector<double,dimworld> x){ return x; }); + BlockVector<TargetSpace> sol(feBasis.size()); + TupleVector<std::vector<RealTuple<double,3> >, + std::vector<Rotation<double,3> > > solTuple; + solTuple[_0].resize(feBasis.size()); + solTuple[_1].resize(feBasis.size()); + for (int i = 0; i < feBasis.size(); i++) { + sol[i] = getConfiguration(id[i]); + solTuple[_0][i] = sol[i].r; + solTuple[_1][i] = sol[i].q; + } + + CosseratVTKWriter<GridType>::write<FEBasis>(feBasis, solTuple, "configuration_l" + std::to_string(numLevels)); + + double energy = 0; + + // A view on the FE basis on a single element + auto localView = feBasis.localView(); + + // Loop over all elements + for (const auto& element : elements(feBasis.gridView(), Dune::Partitions::interior)) + { + localView.bind(element); + + // Number of degrees of freedom on this element + size_t nDofs = localView.tree().size(); + + std::vector<TargetSpace> localSolution(nDofs); + + for (size_t i=0; i<nDofs; i++) + localSolution[i] = sol[localView.index(i)[0]]; + + energy += localCosseratEnergyPlanar.energy(localView, localSolution); + + } + + return energy; +} + +int main(int argc, char** argv) +{ + MPIHelper::instance(argc, argv); + + double energyFine = calculateEnergy(2); + double energyCoarse = calculateEnergy(1); + std::cout << "energyFine: " << energyFine << std::endl; + std::cout << "energyCoarse: " << energyCoarse << std::endl; + + + assert(std::abs(energyFine - energyCoarse) < 1e-3); +} -- GitLab From f0f113e7f4b3b02725d3e8d0494892e17021cc9a Mon Sep 17 00:00:00 2001 From: Lisa Julia Nebel <lisa_julia.nebel@tu-dresden.de> Date: Fri, 2 Jul 2021 08:30:43 +0200 Subject: [PATCH 2/2] Fix error in nonplanarcosseratenergy The contravariant base vectors were not calculated correctly. The contravariant base vectors are the *columns* of the inverse of the covariant matrix, not the rows. To fix this, take the rows of the transpose of inverse of the covariant matrix. --- dune/gfe/nonplanarcosseratshellenergy.hh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/dune/gfe/nonplanarcosseratshellenergy.hh b/dune/gfe/nonplanarcosseratshellenergy.hh index 21a1f2a9e..0eed3750c 100644 --- a/dune/gfe/nonplanarcosseratshellenergy.hh +++ b/dune/gfe/nonplanarcosseratshellenergy.hh @@ -228,7 +228,7 @@ energy(const typename Basis::LocalView& localView, aCovariant[2] = Dune::MatrixVector::crossProduct(aCovariant[0], aCovariant[1]); aCovariant[2] /= aCovariant[2].two_norm(); - auto aContravariant = aCovariant; + auto aContravariant = Dune::GFE::transpose(aCovariant); aContravariant.invert(); Dune::FieldMatrix<double,3,3> a(0); -- GitLab