Newer
Older
#include <dune/common/test/testsuite.hh>
#include <dune/curvedgrid/curvedgrid.hh>
#include <dune/foamgrid/foamgrid.hh>
#include <dune/gmsh4/gmsh4reader.hh>
#include <dune/gmsh4/gridcreators/lagrangegridcreator.hh>
#include <dune/functions/functionspacebases/interpolate.hh>
#include <dune/functions/functionspacebases/lagrangebasis.hh>
#include <dune/functions/functionspacebases/powerbasis.hh>
#include <dune/functions/gridfunctions/analyticgridviewfunction.hh>
#include <dune/gfe/cosseratvtkwriter.hh>
#include <dune/gfe/assemblers/nonplanarcosseratshellenergy.hh>
#include <dune/gfe/spaces/productmanifold.hh>
#include <dune/gfe/spaces/realtuple.hh>
#include <dune/gfe/spaces/rotation.hh>
using namespace Dune;
template <typename FlatGridView, typename CurvedGridView, typename GridGeometry,
typename DeformationFunction, typename OrientationFunction>
double calculateEnergy(const FlatGridView& flatGridView,
const CurvedGridView& curvedGridView,
const GridGeometry curvedGridGeometry,
const DeformationFunction deformationFunction,
const OrientationFunction orientationFunction)
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";
using FlatFEBasis = Functions::LagrangeBasis<FlatGridView,2>;
FlatFEBasis flatFEBasis(flatGridView);
using namespace Dune::Functions::BasisFactory;
using namespace Dune::Indices;
TupleVector<std::vector<RealTuple<double,3> >, std::vector<Rotation<double,3> > > configuration;
configuration[_0].resize(flatFEBasis.size());
configuration[_1].resize(flatFEBasis.size());
/////////////////////////////////////////////////////////////////////////
// FE representation of the deformation field
/////////////////////////////////////////////////////////////////////////
auto curvedGridDeformationBasis = makeBasis(
curvedGridView,
power<3>(
lagrange<2>()
));
BlockVector<FieldVector<double,3> > deformationAsVector(flatFEBasis.size());
Functions::interpolate(curvedGridDeformationBasis, deformationAsVector, deformationFunction);
for (std::size_t i = 0; i < flatFEBasis.size(); i++)
configuration[_0][i].globalCoordinates() = deformationAsVector[i];
/////////////////////////////////////////////////////////////////////////
// FE representation of the microrotation field
/////////////////////////////////////////////////////////////////////////
auto curvedGridQuaternionBasis = makeBasis(
curvedGridView,
power<4>(
lagrange<2>()
));
// Turn matrix-valued function into quaternion-valued function
// TODO: There is no real need to turn orientationFunction into a
// GridViewFunction, but makeComposedGridFunction expects such
// a GridViewFunction.
auto orientationGridViewFunction = Functions::makeAnalyticGridViewFunction(orientationFunction, curvedGridView);
const auto orientationQuaternionFunction
= Functions::makeComposedGridFunction(Rotation<double,3>::matrixToQuaternion,
orientationGridViewFunction);
BlockVector<FieldVector<double,4> > orientationAsVector(flatFEBasis.size());
Functions::interpolate(curvedGridQuaternionBasis, orientationAsVector, orientationQuaternionFunction);
for (std::size_t i = 0; i < flatFEBasis.size(); i++)
configuration[_1][i] = orientationAsVector[i];
/////////////////////////////////////////////////////////////////////////
// Write the configuration to a file (just for debugging)
/////////////////////////////////////////////////////////////////////////
// The deformation and orientation functions need to become a GridViewFunctions,
// otherwise the current CosseratVTKWriter will not accept them.
auto deformationGridViewFunction = Functions::makeAnalyticGridViewFunction(deformationFunction, flatGridView);
// TODO: This is essentially the same as orientationQuaternionFunction
auto orientationQuaternionGridViewFunction = Functions::makeAnalyticGridViewFunction(orientationQuaternionFunction, flatGridView);
// TODO: Write the curved grid, not the flat one
// BUG: The second argument should be the displacement, not the deformation
CosseratVTKWriter<FlatGridView>::write(flatGridView,
deformationGridViewFunction,
orientationQuaternionGridViewFunction,
2, // VTK output element order
"nonplanarcosseratenergytest-result.vtu");
///////////////////////////////////////////////////
// Construct the energy functional
///////////////////////////////////////////////////
using GlobalCoordinate = typename FlatGridView::template Codim<0>::Entity::Geometry::GlobalCoordinate;
auto density = std::make_shared<GFE::CosseratShellDensity<GlobalCoordinate, double> >(materialParameters);
using ShellEnergy = NonplanarCosseratShellEnergy<FlatFEBasis,
3, // Dimension of the target space
double,
GridGeometry>;
ShellEnergy nonplanarCosseratShellEnergy(density, &curvedGridGeometry);
///////////////////////////////////////////////////
// Compute the energy
///////////////////////////////////////////////////
using TargetSpace = GFE::ProductManifold<RealTuple<double,3>,Rotation<double,3> >;
double energy = 0;
// A view on the FE basis on a single element
auto localView = flatFEBasis.localView();
// Loop over all elements
for (const auto& element : elements(flatGridView))
{
localView.bind(element);
// Number of degrees of freedom on this element
size_t nDofs = localView.tree().size();
std::vector<TargetSpace> localConfiguration(nDofs);
for (size_t i=0; i<nDofs; i++)
{
localConfiguration[i][_0] = configuration[_0][localView.index(i)[0]];
localConfiguration[i][_1] = configuration[_1][localView.index(i)[0]];
}
energy += nonplanarCosseratShellEnergy.energy(localView, localConfiguration);
}
return energy;
int main(int argc, char** argv)
{
MPIHelper::instance(argc, argv);
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
//////////////////////////////////
// Create the grid
//////////////////////////////////
using FlatGrid = FoamGrid<2,3>;
GridFactory<FlatGrid> factory;
Gmsh4::LagrangeGridCreator gridGeometry{factory};
Gmsh4Reader reader{gridGeometry};
reader.read(GRID_PATH "/sphere_order2.msh");
auto flatGrid = factory.createGrid();
// Wrap the flat grid to build a curved grid, with Lagrange elements of order 2
CurvedGrid curvedGrid{*flatGrid, gridGeometry, 2};
///////////////////////////////////////////////////////
// Create configurations and check their energies
///////////////////////////////////////////////////////
TestSuite test;
// The reference configuration
// ---------------------------------
auto deformationIdentity = [](FieldVector<double,3> x){
return x;
};
auto orientationIdentity = [](FieldVector<double,3> x) -> FieldMatrix<double,3,3>
{
return {{1,0,0}, {0,1,0}, {0,0,1}};
};
double energyIdentity = calculateEnergy(flatGrid->leafGridView(), curvedGrid.leafGridView(), gridGeometry,
deformationIdentity, orientationIdentity);
test.check(std::fabs(energyIdentity) < 1e-12, "reference configuration has zero energy");
// A translation
// ---------------------------------
auto deformationTranslation = [](FieldVector<double,3> x){
return x + FieldVector<double,3>{1.5,1.5,1.5};
};
double energyTranslation = calculateEnergy(flatGrid->leafGridView(), curvedGrid.leafGridView(), gridGeometry,
deformationTranslation, orientationIdentity);
test.check(std::fabs(energyTranslation) < 1e-12, "translated configuration has zero energy");
// A rotation -- for simplicity about an axis
// ---------------------------------
double angle = M_PI/4;
FieldMatrix<double,3,3> globalRotation = {{1,0,0},
{0,cos(angle),-sin(angle)},
{0,sin(angle),cos(angle)}};
// A translation
auto deformationRotation = [&globalRotation](FieldVector<double,3> x){
FieldVector<double,3> result;
globalRotation.mv(x,result);
return result;
};
auto orientationRotation = [&globalRotation](FieldVector<double,3> x) -> FieldMatrix<double,3,3>
{
return globalRotation;
double energyRotation = calculateEnergy(flatGrid->leafGridView(), curvedGrid.leafGridView(), gridGeometry,
deformationRotation, orientationRotation);
test.check(std::fabs(energyRotation) < 1e-12, "rotated configuration has zero energy");
// Stretching
// ---------------------------------
auto deformationStretchY = [](FieldVector<double,3> x) -> FieldVector<double,3>
{
return {x[0], 2*x[1], x[2]};
double energyStretchY = calculateEnergy(flatGrid->leafGridView(), curvedGrid.leafGridView(), gridGeometry,
deformationStretchY, orientationIdentity);
test.check(std::fabs(energyStretchY-357163.4280328181) < 1e-6,
"stretched configuration has energy 357163.4280328181");
// Something wild
// ----------------------------------------
auto deformationIrregular = [](FieldVector<double,3> x) -> FieldVector<double,3>
{
return {x[0] + sin(3*M_PI*(x[1]+x[2])),
x[1] + sin(5*M_PI*(x[0]+x[2])),
x[2] + sin(2*M_PI*(x[0]+x[1]))};
};
auto orientationIrregular = [](FieldVector<double,3> x) -> FieldMatrix<double,3>
{
double angle = sqrt(x[0]*x[0] + x[1]*x[1]);
return {{1,0,0},
{0,cos(angle),-sin(angle)},
{0,sin(angle),cos(angle)}};
};
double energyIrregular = calculateEnergy(flatGrid->leafGridView(), curvedGrid.leafGridView(), gridGeometry,
deformationIrregular, orientationIrregular);
test.check(std::fabs(energyIrregular-56331610.0823059) < 1e-6,
"irregular configuration has energy 56331610.0823059");