Skip to content
Snippets Groups Projects
Commit 68f8cf14 authored by Lisa Julia Nebel's avatar Lisa Julia Nebel
Browse files

Use the calculation of the vertexnomals from dune-curvedgeometry in film-on-substrate

parent 08f8f940
No related branches found
No related tags found
1 merge request!63Feature/remove vertex normals
Hide whitespace changes
Inline Side-by-side
dune/gfe/surfacecosseratenergy.hh
+ 39
− 38
View file @ 68f8cf14
......@@ -15,11 +15,18 @@
#include <dune/gfe/orthogonalmatrix.hh>
#include <dune/gfe/rigidbodymotion.hh>
#include <dune/gfe/tensor3.hh>
#include <dune/gfe/vertexnormals.hh>
namespace Dune::GFE {
#include <dune/curvedgeometry/curvedgeometry.hh>
#include <dune/localfunctions/lagrange/lfecache.hh>
template<class Basis, class... TargetSpaces>
namespace Dune::GFE {
/** \brief Assembles the cosserat energy on the given boundary for a single element.
*
* \tparam CurvedGeometryGridFunction Type of the grid function that gives the geometry of the deformed surface
* \tparam Basis Type of the Basis used for assembling
* \tparam TargetSpaces The List of TargetSpaces - SurfaceCosseratEnergy needs exactly two TargetSpaces!
*/
template<class CurvedGeometryGridFunction, class Basis, class... TargetSpaces>
class SurfaceCosseratEnergy
: public Dune::GFE::LocalEnergy<Basis, TargetSpaces...>
{
......@@ -70,16 +77,19 @@ public:
/** \brief Constructor with a set of material parameters
* \param parameters The material parameters
* \param shellBoundary The shellBoundary contains the faces where the cosserat energy is assembled
* \param curvedGeometryGridFunction The curvedGeometryGridFunction gives the geometry of the shell in stress-free state.
When assembling, we deform the intersections using the curvedGeometryGridFunction and then use the deformed geometries.
* \param thicknessF The shell thickness parameter, given as a function and evaluated at each quadrature point
* \param lameF The Lame parameters, given as a function and evaluated at each quadrature point
*/
SurfaceCosseratEnergy(const Dune::ParameterTree& parameters,
const std::vector<UnitVector<double,dimWorld> >& vertexNormals,
const BoundaryPatch<GridView>* shellBoundary,
const std::unordered_map<typename GridView::Grid::GlobalIdSet::IdType,Dune::MultiLinearGeometry<double, dimWorld-1, dimWorld>>& geometriesOnShellBoundary,
const CurvedGeometryGridFunction& curvedGeometryGridFunction,
const std::function<double(Dune::FieldVector<double,dimWorld>)> thicknessF,
const std::function<Dune::FieldVector<double,2>(Dune::FieldVector<double,dimWorld>)> lameF)
: shellBoundary_(shellBoundary),
vertexNormals_(vertexNormals),
geometriesOnShellBoundary_(geometriesOnShellBoundary),
curvedGeometryGridFunction_(curvedGeometryGridFunction),
thicknessF_(thicknessF),
lameF_(lameF)
{
......@@ -98,7 +108,7 @@ public:
RT energy(const typename Basis::LocalView& localView,
const std::vector<TargetSpaces>&... localSolutions) const
{
static_assert(sizeof...(TargetSpaces) == 2, "SurfaceCosseratEnergy needs exactly two TargetSpace!");
static_assert(sizeof...(TargetSpaces) == 2, "SurfaceCosseratEnergy needs exactly two TargetSpaces!");
using namespace Dune::Indices;
using TargetSpace0 = typename std::tuple_element<0, std::tuple<TargetSpaces...> >::type;
......@@ -110,21 +120,6 @@ RT energy(const typename Basis::LocalView& localView,
auto element = localView.element();
auto gridView = localView.globalBasis().gridView();
////////////////////////////////////////////////////////////////////////////////////
// Construct a linear (i.e., non-constant!) normal field on each element
////////////////////////////////////////////////////////////////////////////////////
typedef typename Dune::PQkLocalFiniteElementCache<DT, double, gridDim, 1> P1FiniteElementCache;
typedef typename P1FiniteElementCache::FiniteElementType P1LocalFiniteElement;
//it.type()?
P1FiniteElementCache p1FiniteElementCache;
const auto& p1LocalFiniteElement = p1FiniteElementCache.get(element.type());
assert(vertexNormals_.size() == gridView.indexSet().size(gridDim));
std::vector<UnitVector<double,3> > cornerNormals(element.subEntities(gridDim));
for (size_t i=0; i<cornerNormals.size(); i++)
cornerNormals[i] = vertexNormals_[gridView.indexSet().subIndex(element,i,gridDim)];
Dune::GFE::LocalProjectedFEFunction<gridDim, DT, P1LocalFiniteElement, UnitVector<double,3> > unitNormals(p1LocalFiniteElement, cornerNormals);
////////////////////////////////////////////////////////////////////////////////////
// Set up the local nonlinear finite element function
////////////////////////////////////////////////////////////////////////////////////
......@@ -159,14 +154,25 @@ RT energy(const typename Basis::LocalView& localView,
RT energy = 0;
auto& idSet = gridView.grid().globalIdSet();
for (auto&& it : intersections(shellBoundary_->gridView(), element)) {
if (not shellBoundary_->contains(it))
continue;
auto id = idSet.subId(it.inside(), it.indexInInside(), 1);
auto boundaryGeometry = geometriesOnShellBoundary_.at(id);
auto localGridFunction = localFunction(curvedGeometryGridFunction_);
auto curvedGeometryGridFunctionOrder = deformationLocalFiniteElement.localBasis().order();//curvedGeometryGridFunction_.basis().localView().tree().child(0).finiteElement().localBasis().order();
localGridFunction.bind(element);
auto referenceElement = Dune::referenceElement<DT,boundaryDim>(it.type());
// Construct the geometry on the boundary using the map lGF(localGeometry.global(local)):
// The variable local holds the local coordinates in the 2D reference element, localGeometry.global maps them to the 3D reference element.
// The function lGF is the gridfunction bound to the current element, so lGF(localGeometry.global(local)) is the value of curvedGeometryGridFunction_ at
// the point on the intersection face.
using BoundaryGeometry = Dune::CurvedGeometry<DT, boundaryDim, dimWorld, Dune::CurvedGeometryTraits<DT, Dune::LagrangeLFECache<DT,DT,boundaryDim>>>;
BoundaryGeometry boundaryGeometry(referenceElement,
[localGridFunction, localGeometry=it.geometryInInside()](const auto& local) {
return localGridFunction(localGeometry.global(local));
}, curvedGeometryGridFunctionOrder);
auto quadOrder = (it.type().isSimplex()) ? deformationLocalFiniteElement.localBasis().order()
: deformationLocalFiniteElement.localBasis().order() * boundaryDim;
......@@ -243,7 +249,6 @@ RT energy(const typename Basis::LocalView& localView,
// If dimWorld==3, then the first two lines of aCovariant are simply the jacobianTransposed
// of the element. If dimWorld<3 (i.e., ==2), we have to explicitly enters 0.0 in the last column.
const auto jacobianTransposed = boundaryGeometry.jacobianTransposed(quad[pt].position());
// auto jacobianTransposed = geometry.jacobianTransposed(quadPos);
for (int i=0; i<2; i++)
{
......@@ -253,7 +258,7 @@ RT energy(const typename Basis::LocalView& localView,
aCovariant[i][j] = 0.0;
}
aCovariant[2] = Dune::MatrixVector::crossProduct(aCovariant[0], aCovariant[1]);
aCovariant[2] = Dune::FMatrixHelp::Impl::crossProduct(aCovariant[0], aCovariant[1]);
aCovariant[2] /= aCovariant[2].two_norm();
auto aContravariant = aCovariant;
......@@ -277,9 +282,8 @@ RT energy(const typename Basis::LocalView& localView,
for (int beta=0; beta<2; beta++)
c += aScalar * eps[alpha][beta] * Dune::GFE::dyadicProduct(aContravariant[alpha], aContravariant[beta]);
// Second fundamental form
// The derivative of the normal field
auto normalDerivative = unitNormals.evaluateDerivative(quadPos);
// Second fundamental form: The derivative of the normal field
auto normalDerivative = boundaryGeometry.normalGradient(quad[pt].position());
Dune::FieldMatrix<double,3,3> b(0);
for (int alpha=0; alpha<boundaryDim; alpha++)
......@@ -375,11 +379,8 @@ private:
/** \brief The shell boundary */
const BoundaryPatch<GridView>* shellBoundary_;
/** \brief Stress-free geometries of the shell elements*/
const std::unordered_map<typename GridView::Grid::GlobalIdSet::IdType, Dune::MultiLinearGeometry<double, dimWorld-1, dimWorld>> geometriesOnShellBoundary_;
/** \brief The normal vectors at the grid vertices. They are used to compute the reference surface curvature. */
std::vector<UnitVector<double,3> > vertexNormals_;
/** \brief The function used to create the Geometries used for assembling */
const CurvedGeometryGridFunction curvedGeometryGridFunction_;
/** \brief The shell thickness as a function*/
std::function<double(Dune::FieldVector<double,dimWorld>)> thicknessF_;
......
problems/film-on-substrate.parset
+ 1
− 0
View file @ 68f8cf14
......@@ -18,6 +18,7 @@ numLevels = 1
startFromFile = false
pathToGridDeformationFile = ./
gridDeformationFile = deformation
writeOutStressFreeData = true
# When not starting from a file, deformation of the surface shell part can be given here using the gridDeformation function
gridDeformation="[1.3*x[0], x[1], x[2]]"
......
src/film-on-substrate.cc
+ 40
− 65
View file @ 68f8cf14
......@@ -47,7 +47,6 @@
#include <dune/gfe/neumannenergy.hh>
#include <dune/gfe/surfacecosseratenergy.hh>
#include <dune/gfe/sumenergy.hh>
#include <dune/gfe/vertexnormals.hh>
#if MIXED_SPACE
#include <dune/gfe/mixedriemanniantrsolver.hh>
......@@ -63,6 +62,7 @@
#include <dune/solvers/solvers/iterativesolver.hh>
#include <dune/solvers/norms/energynorm.hh>
// grid dimension
#ifndef WORLD_DIM
# define WORLD_DIM 3
......@@ -74,6 +74,8 @@ const int targetDim = WORLD_DIM;
const int displacementOrder = 2;
const int rotationOrder = 2;
const int stressFreeDataOrder = 2;
#if !MIXED_SPACE
static_assert(displacementOrder==rotationOrder, "displacement and rotation order do not match!");
#endif
......@@ -343,32 +345,28 @@ int main (int argc, char *argv[]) try
vtkWriter.write(resultPath + "finite-strain_homotopy_" + parameterSet.get<std::string>("energy") + "_0");
/////////////////////////////////////////////////////////////
// INITIAL SURFACE SHELL DATA
// STRESS-FREE SURFACE SHELL DATA
/////////////////////////////////////////////////////////////
auto stressFreeFEBasis = makeBasis(
gridView,
power<dim>(
lagrange<stressFreeDataOrder>(),
blockedInterleaved()
));
typedef MultiLinearGeometry<double, dim-1, dim> ML;
std::unordered_map<GridType::GlobalIdSet::IdType, ML> geometriesOnShellBoundary;
auto& idSet = grid->globalIdSet();
GlobalIndexSet<GridView> globalVertexIndexSet(gridView,dim);
BlockVector<FieldVector<double,dim> > stressFreeShellVector(stressFreeFEBasis.size());
// Read in the grid deformation
if (startFromFile) {
// Create a basis of order 1 in order to deform the geometries on the surface shell boundary
const std::string pathToGridDeformationFile = parameterSet.get("pathToGridDeformationFile", "");
// for this, we create a basis of order 1 in order to deform the geometries on the surface shell boundary
auto feBasisOrder1 = makeBasis(
gridView,
power<dim>(
lagrange<1>(),
blockedInterleaved()
));
GlobalIndexSet<GridView> globalVertexIndexSet(gridView,dim);
std::unordered_map<std::string, FieldVector<double,3>> deformationMap;
std::string line, displacement, entry;
if (mpiHelper.rank() == 0)
std::cout << "Reading in deformation file: " << pathToGridDeformationFile + parameterSet.get<std::string>("gridDeformationFile") << std::endl;
std::cout << "Reading in deformation file (" << "order is " << stressFreeDataOrder << "): " << pathToGridDeformationFile + parameterSet.get<std::string>("gridDeformationFile") << std::endl;
// Read grid deformation information from the file specified in the parameter set via gridDeformationFile
std::ifstream file(pathToGridDeformationFile + parameterSet.get<std::string>("gridDeformationFile"), std::ios::in);
if (file.is_open()) {
while (std::getline(file, line)) {
......@@ -384,62 +382,40 @@ int main (int argc, char *argv[]) try
}
if (mpiHelper.rank() == 0)
std::cout << "... done: The grid has " << globalVertexIndexSet.size(dim) << " vertices and the defomation file has " << deformationMap.size() << " entries." << std::endl;
if (deformationMap.size() != globalVertexIndexSet.size(dim))
if (stressFreeDataOrder == 1 && deformationMap.size() != globalVertexIndexSet.size(dim))
DUNE_THROW(Exception, "Error: Grid and deformation vector do not match!");
file.close();
} else {
DUNE_THROW(Exception, "Error: Could not open the file containing the deformation vector!");
}
BlockVector<FieldVector<double,dim>> gridDeformationFromFile;
Dune::Functions::interpolate(feBasisOrder1, gridDeformationFromFile, [](FieldVector<double,dim> x){ return x; });
Dune::Functions::interpolate(stressFreeFEBasis, stressFreeShellVector, [](FieldVector<double,dim> x){ return x; });
for (auto& entry : gridDeformationFromFile) {
for (auto& entry : stressFreeShellVector) {
std::stringstream stream;
stream << entry;
entry = deformationMap.at(stream.str()); //Look up the deformation for this vertex in the deformationMap
}
auto gridDeformationFromFileFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim>>(feBasisOrder1, gridDeformationFromFile);
auto localGridDeformationFromFileFunction = localFunction(gridDeformationFromFileFunction);
//Write out the stress-free geometries that were read in
SubsamplingVTKWriter<GridView> vtkWriter(gridView, Dune::refinementLevels(0));
vtkWriter.addVertexData(localGridDeformationFromFileFunction, VTK::FieldInfo("displacement", VTK::FieldInfo::Type::scalar, dim));
vtkWriter.write("stress-free-geometries");
//Iterate over boundary, each facet on the boundary has an element (boundaryElement.inside()) with a unique global id (idSet.subId);
//we store the new geometry in the map with this id as reference
for (auto boundaryElement : surfaceShellBoundary) {
localGridDeformationFromFileFunction.bind(boundaryElement.inside());
std::vector<Dune::FieldVector<double,dim>> corners;
for (int i = 0; i < boundaryElement.geometry().corners(); i++) {
auto corner = boundaryElement.geometry().corner(i);
corner += localGridDeformationFromFileFunction(boundaryElement.inside().geometry().local(boundaryElement.geometry().corner(i)));
corners.push_back(corner);
}
localGridDeformationFromFileFunction.unbind();
GridType::GlobalIdSet::IdType id = idSet.subId(boundaryElement.inside(), boundaryElement.indexInInside(), 1);
ML boundaryGeometry(boundaryElement.geometry().type(), corners);
geometriesOnShellBoundary.insert({id, boundaryGeometry});
entry += deformationMap.at(stream.str()); //Look up the displacement for this vertex in the deformationMap
}
} else {
// Read grid deformation from deformation function
auto gridDeformationLambda = std::string("lambda x: (") + parameterSet.get<std::string>("gridDeformation") + std::string(")");
auto gridDeformation = Python::make_function<FieldVector<double,dim> >(Python::evaluate(gridDeformationLambda));
Dune::Functions::interpolate(stressFreeFEBasis, stressFreeShellVector, gridDeformation);
}
//Iterate over boundary, each facet on the boundary has an element (boundaryElement.inside()) with a unique global id (idSet.subId);
//we store the new geometry in the map with this id as reference
for (auto boundaryElement : surfaceShellBoundary) {
std::vector<Dune::FieldVector<double,dim>> corners;
for (int i = 0; i < boundaryElement.geometry().corners(); i++) {
auto corner = gridDeformation(boundaryElement.geometry().corner(i));
corners.push_back(corner);
}
GridType::GlobalIdSet::IdType id = idSet.subId(boundaryElement.inside(), boundaryElement.indexInInside(), 1);
ML boundaryGeometry(boundaryElement.geometry().type(), corners);
geometriesOnShellBoundary.insert({id, boundaryGeometry});
auto stressFreeShellFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim>>(stressFreeFEBasis, stressFreeShellVector);
if (parameterSet.hasKey("writeOutStressFreeData") && parameterSet.get<bool>("writeOutStressFreeData")) {
BlockVector<FieldVector<double,dim> > stressFreeDisplacement(stressFreeFEBasis.size());
Dune::Functions::interpolate(stressFreeFEBasis, stressFreeDisplacement, [](FieldVector<double,dim> x){ return (-1.0)*x; });
for (int i = 0; i < stressFreeFEBasis.size(); i++) {
stressFreeDisplacement[i] += stressFreeShellVector[i];
}
auto stressFreeDisplacementFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim>>(stressFreeFEBasis, stressFreeDisplacement);
//Write out the stress-free shell function that was read in
SubsamplingVTKWriter<GridView> vtkWriterStressFree(gridView, Dune::refinementLevels(1));
vtkWriterStressFree.addVertexData(stressFreeDisplacementFunction, VTK::FieldInfo("displacement", VTK::FieldInfo::Type::scalar, dim));
vtkWriterStressFree.write("stress-free-shell-function");
}
/////////////////////////////////////////////////////////////
......@@ -450,13 +426,6 @@ int main (int argc, char *argv[]) try
neumannValues = parameterSet.get<FieldVector<double,dim> >("neumannValues");
std::cout << "Neumann values: " << neumannValues << std::endl;
// Vertex Normals for the 3D-Part
std::vector<UnitVector<double,dim> > vertexNormals(gridView.size(dim));
Dune::FieldVector<double,dim> vertexNormalRaw = {0,0,1};
for (int i = 0; i< vertexNormals.size(); i++) {
UnitVector vertexNormal(vertexNormalRaw);
vertexNormals[i] = vertexNormal;
}
//Function for the Cosserat material parameters
const ParameterTree& materialParameters = parameterSet.sub("materialParameters");
Python::Reference surfaceShellClass = Python::import(materialParameters.get<std::string>("surfaceShellParameters"));
......@@ -512,7 +481,13 @@ int main (int argc, char *argv[]) try
auto elasticEnergy = std::make_shared<GFE::LocalIntegralEnergy<CompositeBasis, RealTuple<ValueType,targetDim>, Rotation<ValueType,dim>>>(elasticDensity);
auto neumannEnergy = std::make_shared<GFE::NeumannEnergy<CompositeBasis, RealTuple<ValueType,targetDim>, Rotation<ValueType,dim>>>(neumannBoundary,*neumannFunctionPtr);
auto surfaceCosseratEnergy = std::make_shared<GFE::SurfaceCosseratEnergy<CompositeBasis, RealTuple<ValueType,dim>, Rotation<ValueType,dim> >>(materialParameters, std::move(vertexNormals), &surfaceShellBoundary, std::move(geometriesOnShellBoundary), fThickness, fLame);
auto surfaceCosseratEnergy = std::make_shared<GFE::SurfaceCosseratEnergy<
decltype(stressFreeShellFunction), CompositeBasis, RealTuple<ValueType,dim>, Rotation<ValueType,dim> >>(
materialParameters,
&surfaceShellBoundary,
stressFreeShellFunction,
fThickness,
fLame);
GFE::SumEnergy<CompositeBasis, RealTuple<ValueType,targetDim>, Rotation<ValueType,targetDim>> sumEnergy;
sumEnergy.addLocalEnergy(neumannEnergy);
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment