From 71c325be56ff37db2dab80169b04f38b4d83d1ae Mon Sep 17 00:00:00 2001 From: Oliver Sander <sander@igpm.rwth-aachen.de> Date: Wed, 12 Oct 2011 16:39:21 +0000 Subject: [PATCH] Support for translation Neumann values [[Imported from SVN: r7888]] --- dune/gfe/cosseratenergystiffness.hh | 52 ++++++++++++++++++++++++++++- 1 file changed, 51 insertions(+), 1 deletion(-) diff --git a/dune/gfe/cosseratenergystiffness.hh b/dune/gfe/cosseratenergystiffness.hh index 1ab4a2dc..067120cb 100644 --- a/dune/gfe/cosseratenergystiffness.hh +++ b/dune/gfe/cosseratenergystiffness.hh @@ -5,6 +5,8 @@ #include <dune/common/parametertree.hh> #include <dune/grid/common/quadraturerules.hh> +#include <dune/fufem/functions/virtualgridfunction.hh> + #include "localgeodesicfestiffness.hh" #include "localgeodesicfefunction.hh" #include <dune/gfe/rigidbodymotion.hh> @@ -101,7 +103,11 @@ public: /** \brief Constructor with a set of material parameters * \param parameters The material parameters */ - CosseratEnergyLocalStiffness(const Dune::ParameterTree& parameters) + CosseratEnergyLocalStiffness(const Dune::ParameterTree& parameters, + const BoundaryPatch<GridView>* neumannBoundary, + const Dune::VirtualFunction<Dune::FieldVector<double,gridDim>, Dune::FieldVector<double,3> >* neumannFunction) + : neumannBoundary_(neumannBoundary), + neumannFunction_(neumannFunction) { // The shell thickness thickness_ = parameters.template get<double>("thickness"); @@ -182,7 +188,11 @@ public: /** \brief Curvature exponent */ double q_; + /** \brief The Neumann boundary */ + const BoundaryPatch<GridView>* neumannBoundary_; + /** \brief The function implementing the Neumann data */ + const Dune::VirtualFunction<Dune::FieldVector<double,gridDim>, Dune::FieldVector<double,3> >* neumannFunction_; }; template <class GridView, class LocalFiniteElement, int dim> @@ -278,7 +288,47 @@ energy(const Entity& element, DUNE_THROW(Dune::NotImplemented, "CosseratEnergyStiffness for 1d grids"); } + + + ////////////////////////////////////////////////////////////////////////////// + // Assemble boundary contributions + ////////////////////////////////////////////////////////////////////////////// + + assert(neumannFunction_); + + for (typename Entity::LeafIntersectionIterator it = element.ileafbegin(); it != element.ileafend(); ++it) { + + if (not neumannBoundary_ or not neumannBoundary_->contains(*it)) + continue; + + const Dune::QuadratureRule<double, gridDim-1>& quad + = Dune::QuadratureRules<double, gridDim-1>::rule(it->type(), quadOrder); + + for (size_t pt=0; pt<quad.size(); pt++) { + + // Local position of the quadrature point + const Dune::FieldVector<double,gridDim>& quadPos = it->geometryInInside().global(quad[pt].position()); + + const double integrationElement = it->geometry().integrationElement(quad[pt].position()); + + // The value of the local function + RigidBodyMotion<dim> value = localGeodesicFEFunction.evaluate(quadPos); + // Value of the Neumann data at the current position + Dune::FieldVector<double,3> neumannValue; + + if (dynamic_cast<const VirtualGridViewFunction<GridView,Dune::FieldVector<double,3> >*>(neumannFunction_)) + dynamic_cast<const VirtualGridViewFunction<GridView,Dune::FieldVector<double,3> >*>(neumannFunction_)->evaluateLocal(element, quadPos, neumannValue); + else + neumannFunction_->evaluate(it->geometry().global(quad[pt].position()), neumannValue); + + // Constant Neumann force in z-direction + energy += thickness_ * (neumannValue * value.r) * quad[pt].weight() * integrationElement; + + } + + } + return energy; } -- GitLab