From 477f9510d75ae3d9f61c678d4c69da3a13f237c4 Mon Sep 17 00:00:00 2001
From: Oliver Sander <oliver.sander@tu-dresden.de>
Date: Fri, 22 Jan 2016 06:01:53 +0100
Subject: [PATCH] Implement support for volume loads
---
dune/gfe/cosseratenergystiffness.hh | 24 +++++++++++++++++++++++-
src/cosserat-continuum.cc | 26 +++++++++++++++++++++++++-
2 files changed, 48 insertions(+), 2 deletions(-)
diff --git a/dune/gfe/cosseratenergystiffness.hh b/dune/gfe/cosseratenergystiffness.hh
index 0f33ad4b..989e5e50 100644
--- a/dune/gfe/cosseratenergystiffness.hh
+++ b/dune/gfe/cosseratenergystiffness.hh
@@ -189,7 +189,8 @@ public:
*/
CosseratEnergyLocalStiffness(const Dune::ParameterTree& parameters,
const BoundaryPatch<GridView>* neumannBoundary,
- const Dune::VirtualFunction<Dune::FieldVector<double,dimworld>, Dune::FieldVector<double,3> >* neumannFunction)
+ const 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),
neumannFunction_(neumannFunction)
{
@@ -363,6 +364,9 @@ public:
/** \brief The function implementing the Neumann data */
const Dune::VirtualFunction<Dune::FieldVector<double,dimworld>, Dune::FieldVector<double,3> >* neumannFunction_;
+
+ /** \brief The function implementing a volume load */
+ const std::shared_ptr<Dune::VirtualFunction<Dune::FieldVector<double,dimworld>, Dune::FieldVector<double,3> > > volumeLoad_;
};
template <class Basis, int dim, class field_type>
@@ -448,6 +452,24 @@ energy(const typename Basis::LocalView& localView,
} else
DUNE_THROW(Dune::NotImplemented, "CosseratEnergyStiffness for 1d grids");
+ ///////////////////////////////////////////////////////////
+ // Volume load contribution
+ ///////////////////////////////////////////////////////////
+
+ if (not volumeLoad_)
+ continue;
+
+ // Value of the volume load density at the current position
+ Dune::FieldVector<double,3> volumeLoadDensity;
+
+ if (std::dynamic_pointer_cast<const VirtualGridViewFunction<GridView,Dune::FieldVector<double,3> > >(volumeLoad_))
+ std::dynamic_pointer_cast<const VirtualGridViewFunction<GridView,Dune::FieldVector<double,3> > >(volumeLoad_)->evaluateLocal(element, quadPos, volumeLoadDensity);
+ else
+ volumeLoad_->evaluate(element.geometry().global(quad[pt].position()), volumeLoadDensity);
+
+ // Only translational dofs are affected by the volume load
+ for (size_t i=0; i<volumeLoadDensity.size(); i++)
+ energy -= thickness_ * (volumeLoadDensity[i] * value.r[i]) * quad[pt].weight() * integrationElement;
}
diff --git a/src/cosserat-continuum.cc b/src/cosserat-continuum.cc
index c2126be1..75eea015 100644
--- a/src/cosserat-continuum.cc
+++ b/src/cosserat-continuum.cc
@@ -76,6 +76,25 @@ struct NeumannFunction
double homotopyParameter_;
};
+/** \brief A constant vector-valued function, for simple volume loads */
+struct VolumeLoad
+ : public Dune::VirtualFunction<FieldVector<double,dimworld>, FieldVector<double,3> >
+{
+ VolumeLoad(const FieldVector<double,3> values,
+ double homotopyParameter)
+ : values_(values),
+ homotopyParameter_(homotopyParameter)
+ {}
+
+ void evaluate(const FieldVector<double, dimworld>& x, FieldVector<double,3>& out) const {
+ out = 0;
+ out.axpy(homotopyParameter_, values_);
+ }
+
+ FieldVector<double,3> values_;
+ double homotopyParameter_;
+};
+
int main (int argc, char *argv[]) try
{
@@ -260,6 +279,10 @@ int main (int argc, char *argv[]) try
neumannFunction = make_shared<NeumannFunction>(parameterSet.get<FieldVector<double,3> >("neumannValues"),
homotopyParameter);
+ shared_ptr<VolumeLoad> volumeLoad;
+ if (parameterSet.hasKey("volumeLoad"))
+ volumeLoad = make_shared<VolumeLoad>(parameterSet.get<FieldVector<double,3> >("volumeLoad"),
+ homotopyParameter);
if (mpiHelper.rank() == 0) {
std::cout << "Material parameters:" << std::endl;
@@ -270,7 +293,8 @@ int main (int argc, char *argv[]) try
CosseratEnergyLocalStiffness<FEBasis,
3,adouble> cosseratEnergyADOLCLocalStiffness(materialParameters,
&neumannBoundary,
- neumannFunction.get());
+ neumannFunction.get(),
+ volumeLoad);
LocalGeodesicFEADOLCStiffness<FEBasis,
TargetSpace> localGFEADOLCStiffness(&cosseratEnergyADOLCLocalStiffness);
--
GitLab