From 543146ef8130ed4cf4707e59d9e9cfef5dfeaedb Mon Sep 17 00:00:00 2001
From: Oliver Sander <oliver.sander@tu-dresden.de>
Date: Tue, 26 Jan 2016 06:44:10 +0100
Subject: [PATCH] Implement general surface normal fields for shells with
nonplanar reference geometry
Previously, the NonplanarCosseratShellEnergy had the normal field of the unit sphere
hard-wired in the code. With this patch, using general triangulated surface as
reference geometries is possible.
---
dune/gfe/nonplanarcosseratshellenergy.hh | 15 ++++---
dune/gfe/vertexnormals.hh | 53 ++++++++++++++++++++++++
src/cosserat-continuum.cc | 3 ++
3 files changed, 65 insertions(+), 6 deletions(-)
create mode 100644 dune/gfe/vertexnormals.hh
diff --git a/dune/gfe/nonplanarcosseratshellenergy.hh b/dune/gfe/nonplanarcosseratshellenergy.hh
index a0b80b63..b995ccaa 100644
--- a/dune/gfe/nonplanarcosseratshellenergy.hh
+++ b/dune/gfe/nonplanarcosseratshellenergy.hh
@@ -181,10 +181,12 @@ public:
* \param parameters The material parameters
*/
NonplanarCosseratShellEnergy(const Dune::ParameterTree& parameters,
+ const std::vector<UnitVector<double,3> >& vertexNormals,
const BoundaryPatch<GridView>* neumannBoundary,
const std::shared_ptr<Dune::VirtualFunction<Dune::FieldVector<double,dimworld>, Dune::FieldVector<double,3> > > neumannFunction,
const std::shared_ptr<Dune::VirtualFunction<Dune::FieldVector<double,dimworld>, Dune::FieldVector<double,3> > > volumeLoad)
- : neumannBoundary_(neumannBoundary),
+ : vertexNormals_(vertexNormals),
+ neumannBoundary_(neumannBoundary),
neumannFunction_(neumannFunction),
volumeLoad_(volumeLoad)
{
@@ -248,6 +250,9 @@ public:
/** \brief Curvature parameters */
double b1_, b2_, b3_;
+ /** \brief The normal vectors at the grid vertices. This are used to compute the reference surface curvature. */
+ std::vector<UnitVector<double,3> > vertexNormals_;
+
/** \brief The Neumann boundary */
const BoundaryPatch<GridView>* neumannBoundary_;
@@ -274,14 +279,12 @@ energy(const typename Basis::LocalView& localView,
////////////////////////////////////////////////////////////////////////////////////
// Construct a linear (i.e., non-constant!) normal field on the surface
////////////////////////////////////////////////////////////////////////////////////
+ auto gridView = localView.globalBasis().gridView();
+ 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++)
- {
- // TODO: WARNING WARNING WARNING: Shape is hard-coded here!
- // The sphere:
- cornerNormals[i] = geometry.corner(i) / geometry.corner(i).two_norm();
- }
+ cornerNormals[i] = vertexNormals_[gridView.indexSet().subIndex(element,i,2)];
typedef typename Dune::PQkLocalFiniteElementCache<DT, double, gridDim, 1> P1FiniteElementCache;
typedef typename P1FiniteElementCache::FiniteElementType P1LocalFiniteElement;
diff --git a/dune/gfe/vertexnormals.hh b/dune/gfe/vertexnormals.hh
new file mode 100644
index 00000000..fb899e0e
--- /dev/null
+++ b/dune/gfe/vertexnormals.hh
@@ -0,0 +1,53 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_GFE_VERTEXNORMALS_HH
+#define DUNE_GFE_VERTEXNORMALS_HH
+
+#include <vector>
+
+#include <dune/common/fvector.hh>
+
+#include <dune/geometry/referenceelements.hh>
+
+#include <dune/gfe/unitvector.hh>
+
+/** \brief Compute averaged vertex normals for a 2d-in-3d grid
+ */
+template <class GridView>
+std::vector<UnitVector<typename GridView::ctype,3> > computeVertexNormals(const GridView& gridView)
+{
+ const auto& indexSet = gridView.indexSet();
+
+ assert(GridView::dimension == 2);
+
+ std::vector<Dune::FieldVector<typename GridView::ctype,3> > unscaledNormals(indexSet.size(2));
+ std::fill(unscaledNormals.begin(), unscaledNormals.end(), Dune::FieldVector<typename GridView::ctype,3>(0.0));
+
+ for (const auto& element : elements(gridView))
+ {
+ for (int i=0; i<element.subEntities(2); i++)
+ {
+ auto cornerPos = Dune::ReferenceElements<double,2>::general(element.type()).position(i,2);
+ auto tangent = element.geometry().jacobianTransposed(cornerPos);
+ auto cornerNormal = Arithmetic::crossProduct(tangent[0], tangent[1]);
+ cornerNormal /= cornerNormal.two_norm();
+
+ unscaledNormals[indexSet.subIndex(element,i,2)] += cornerNormal;
+ }
+ }
+
+ // Normalize
+ for (auto& normal : unscaledNormals)
+ normal /= normal.two_norm();
+
+ // Convert array of vectors to array of UnitVectors
+ std::vector<UnitVector<typename GridView::ctype,3> > result(indexSet.size(2));
+
+ for (std::size_t i=0; i<unscaledNormals.size(); i++)
+ result[i] = unscaledNormals[i];
+
+ return result;
+}
+
+
+#endif
\ No newline at end of file
diff --git a/src/cosserat-continuum.cc b/src/cosserat-continuum.cc
index 58aafbd8..ec7e405d 100644
--- a/src/cosserat-continuum.cc
+++ b/src/cosserat-continuum.cc
@@ -43,6 +43,7 @@
#include <dune/gfe/vtkreader.hh>
#include <dune/gfe/geodesicfeassembler.hh>
#include <dune/gfe/riemanniantrsolver.hh>
+#include <dune/gfe/vertexnormals.hh>
// grid dimension
const int dim = 2;
@@ -311,7 +312,9 @@ int main (int argc, char *argv[]) try
}
else
{
+ std::vector<UnitVector<double,3> > vertexNormals = computeVertexNormals(gridView);
cosseratEnergyADOLCLocalStiffness = std::make_shared<NonplanarCosseratShellEnergy<FEBasis,3,adouble> >(materialParameters,
+ std::move(vertexNormals),
&neumannBoundary,
neumannFunction,
volumeLoad);
--
GitLab