Skip to content
Snippets Groups Projects
Commit 8830565c authored by Oliver Sander's avatar Oliver Sander Committed by sander
Browse files

Major cleanup and decrufting 

[[Imported from SVN: r10114]]
parent 05a013e0
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
dune/gfe/henckyenergy.hh
+ 15
− 241
View file @ 8830565c
...@@ -2,6 +2,8 @@ ...@@ -2,6 +2,8 @@
#define DUNE_GFE_HENCKYENERGY_HH #define DUNE_GFE_HENCKYENERGY_HH
#include <dune/common/fmatrix.hh> #include <dune/common/fmatrix.hh>
#include <dune/common/fmatrixev.hh>
#include <dune/geometry/quadraturerules.hh> #include <dune/geometry/quadraturerules.hh>
#include <dune/fufem/functions/virtualgridfunction.hh> #include <dune/fufem/functions/virtualgridfunction.hh>
...@@ -11,12 +13,9 @@ ...@@ -11,12 +13,9 @@
#include <dune/gfe/localfestiffness.hh> #include <dune/gfe/localfestiffness.hh>
#include <dune/gfe/localgeodesicfefunction.hh> #include <dune/gfe/localgeodesicfefunction.hh>
#include <dune/gfe/realtuple.hh> #include <dune/gfe/realtuple.hh>
#include <dune/gfe/eigenvalues.hh>
namespace Dune { namespace Dune {
namespace Fufem {
template<class GridView, class LocalFiniteElement, class field_type=double> template<class GridView, class LocalFiniteElement, class field_type=double>
class HenckyEnergy class HenckyEnergy
: public LocalFEStiffness<GridView,LocalFiniteElement,std::vector<Dune::FieldVector<field_type,3> > > : public LocalFEStiffness<GridView,LocalFiniteElement,std::vector<Dune::FieldVector<field_type,3> > >
...@@ -48,8 +47,7 @@ public: // for testing ...@@ -48,8 +47,7 @@ public: // for testing
/** \brief Assemble the energy for a single element */ /** \brief Assemble the energy for a single element */
field_type energy (const Entity& e, field_type energy (const Entity& e,
const LocalFiniteElement& localFiniteElement, const LocalFiniteElement& localFiniteElement,
const std::vector<Dune::FieldVector<field_type,gridDim> >& localConfiguration, const std::vector<Dune::FieldVector<field_type,gridDim> >& localConfiguration) const;
const std::vector<Dune::FieldVector<double,gridDim> >& localPointLoads) const;
/** \brief Lame constants */ /** \brief Lame constants */
double mu_, lambda_; double mu_, lambda_;
...@@ -66,17 +64,13 @@ field_type ...@@ -66,17 +64,13 @@ field_type
HenckyEnergy<GridView,LocalFiniteElement,field_type>:: HenckyEnergy<GridView,LocalFiniteElement,field_type>::
energy(const Entity& element, energy(const Entity& element,
const LocalFiniteElement& localFiniteElement, const LocalFiniteElement& localFiniteElement,
const std::vector<Dune::FieldVector<field_type,gridDim> >& localConfiguration, const std::vector<Dune::FieldVector<field_type,gridDim> >& localConfiguration) const
const std::vector<Dune::FieldVector<double,gridDim> >& localPointLoads) const
{ {
assert(element.type() == localFiniteElement.type()); assert(element.type() == localFiniteElement.type());
typedef typename GridView::template Codim<0>::Entity::Geometry Geometry; typedef typename GridView::template Codim<0>::Entity::Geometry Geometry;
field_type energy = 0; field_type energy = 0;
// typedef LocalGeodesicFEFunction<gridDim, DT, LocalFiniteElement, TargetSpace> LocalGFEFunctionType;
// LocalGFEFunctionType localGeodesicFEFunction(localFiniteElement,localConfiguration);
// store gradients of shape functions and base functions // store gradients of shape functions and base functions
std::vector<Dune::FieldMatrix<DT,1,gridDim> > referenceGradients(localFiniteElement.localBasis().size()); std::vector<Dune::FieldMatrix<DT,1,gridDim> > referenceGradients(localFiniteElement.localBasis().size());
std::vector<Dune::FieldVector<DT,gridDim> > gradients(localFiniteElement.localBasis().size()); std::vector<Dune::FieldVector<DT,gridDim> > gradients(localFiniteElement.localBasis().size());
...@@ -98,28 +92,13 @@ energy(const Entity& element, ...@@ -98,28 +92,13 @@ energy(const Entity& element,
DT weight = quad[pt].weight() * integrationElement; DT weight = quad[pt].weight() * integrationElement;
#if 1
// get gradients of shape functions // get gradients of shape functions
localFiniteElement.localBasis().evaluateJacobian(quadPos, referenceGradients); localFiniteElement.localBasis().evaluateJacobian(quadPos, referenceGradients);
// compute gradients of base functions // compute gradients of base functions
for (size_t i=0; i<gradients.size(); ++i) { for (size_t i=0; i<gradients.size(); ++i)
// transform gradients
jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]); jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]);
}
#endif
#if 0
// The derivative of the local function defined on the reference element
typename LocalGFEFunctionType::DerivativeType referenceDerivative = localGeodesicFEFunction.evaluateDerivative(quadPos);
// The derivative of the function defined on the actual element
typename LocalGFEFunctionType::DerivativeType derivative(0);
for (size_t comp=0; comp<referenceDerivative.N(); comp++)
jacobianInverseTransposed.umv(referenceDerivative[comp], derivative[comp]);
#endif
Dune::FieldMatrix<field_type,gridDim,gridDim> derivative(0); Dune::FieldMatrix<field_type,gridDim,gridDim> derivative(0);
for (size_t i=0; i<gradients.size(); i++) for (size_t i=0; i<gradients.size(); i++)
for (int j=0; j<gridDim; j++) for (int j=0; j<gridDim; j++)
...@@ -138,13 +117,9 @@ energy(const Entity& element, ...@@ -138,13 +117,9 @@ energy(const Entity& element,
////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////
// Eigenvalues of FTF // Eigenvalues of FTF
////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////
#if 0 // hencky
std::array<field_type,dim> lambda = eigenValues(FTF);
////////////////////////////////////////////////////////// Dune::FieldVector<field_type,dim> lambda;
// Compute the derivative of the rotation FMatrixHelp::eigenValues(FTF, lambda);
// Note: we need it in matrix coordinates
//////////////////////////////////////////////////////////
// logarithms of the eigenvalues // logarithms of the eigenvalues
std::array<field_type,dim> ln; std::array<field_type,dim> ln;
...@@ -160,53 +135,35 @@ energy(const Entity& element, ...@@ -160,53 +135,35 @@ energy(const Entity& element,
trace += ln[i]; trace += ln[i];
energy += weight * 0.5 * lambda_ * trace * trace; energy += weight * 0.5 * lambda_ * trace * trace;
#else // St.Venant-Kirchhoff
Dune::FieldMatrix<field_type,dim,dim> E = FTF;
for (int i=0; i<dim; i++)
E[i][i] -= 1.0;
E *= 0.5;
field_type trE = E[0][0] + E[1][1] + E[2][2];
Dune::FieldMatrix<field_type,dim,dim> ESquared = E*E;
field_type trESquared = ESquared[0][0] + ESquared[1][1] + ESquared[2][2];
energy += weight * mu_ * trESquared + weight * 0.5 * lambda_ * trE * trE;
#endif
} }
////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////
// Assemble boundary contributions // Assemble boundary contributions
////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////
for (size_t i=0; i<localPointLoads.size(); i++)
for (size_t j=0; j<dim; j++)
energy -= localConfiguration[i][j] * localPointLoads[i][j];
if (not neumannFunction_) if (not neumannFunction_)
return energy; return energy;
for (typename Entity::LeafIntersectionIterator it = element.ileafbegin(); it != element.ileafend(); ++it) { for (auto&& it : intersections(neumannBoundary_->gridView(), element)) {
if (not neumannBoundary_ or not neumannBoundary_->contains(*it)) if (not neumannBoundary_ or not neumannBoundary_->contains(it))
continue; continue;
const Dune::QuadratureRule<DT, gridDim-1>& quad const Dune::QuadratureRule<DT, gridDim-1>& quad
= Dune::QuadratureRules<DT, gridDim-1>::rule(it->type(), quadOrder); = Dune::QuadratureRules<DT, gridDim-1>::rule(it.type(), quadOrder);
for (size_t pt=0; pt<quad.size(); pt++) { for (size_t pt=0; pt<quad.size(); pt++) {
// Local position of the quadrature point // Local position of the quadrature point
const Dune::FieldVector<DT,gridDim>& quadPos = it->geometryInInside().global(quad[pt].position()); const Dune::FieldVector<DT,gridDim>& quadPos = it.geometryInInside().global(quad[pt].position());
const DT integrationElement = it->geometry().integrationElement(quad[pt].position()); const DT integrationElement = it.geometry().integrationElement(quad[pt].position());
// The value of the local function // The value of the local function
//RealTuple<field_type,dim> value = localGeodesicFEFunction.evaluate(quadPos);
// get gradients of shape functions
std::vector<Dune::FieldVector<DT,1> > shapeFunctionValues; std::vector<Dune::FieldVector<DT,1> > shapeFunctionValues;
localFiniteElement.localBasis().evaluateFunction(quadPos, shapeFunctionValues); localFiniteElement.localBasis().evaluateFunction(quadPos, shapeFunctionValues);
Dune::FieldVector<field_type,dim> value(0); Dune::FieldVector<field_type,dim> value(0);
for (int i=0; i<localFiniteElement.size(); i++) for (int i=0; i<localFiniteElement.size(); i++)
for (int j=0; j<dim; j++) for (int j=0; j<dim; j++)
...@@ -218,7 +175,7 @@ energy(const Entity& element, ...@@ -218,7 +175,7 @@ energy(const Entity& element,
if (dynamic_cast<const VirtualGridViewFunction<GridView,Dune::FieldVector<double,3> >*>(neumannFunction_)) 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); dynamic_cast<const VirtualGridViewFunction<GridView,Dune::FieldVector<double,3> >*>(neumannFunction_)->evaluateLocal(element, quadPos, neumannValue);
else else
neumannFunction_->evaluate(it->geometry().global(quad[pt].position()), neumannValue); neumannFunction_->evaluate(it.geometry().global(quad[pt].position()), neumannValue);
// Only translational dofs are affected by the Neumann force // Only translational dofs are affected by the Neumann force
for (size_t i=0; i<neumannValue.size(); i++) for (size_t i=0; i<neumannValue.size(); i++)
...@@ -231,191 +188,8 @@ energy(const Entity& element, ...@@ -231,191 +188,8 @@ energy(const Entity& element,
return energy; return energy;
} }
} // namespace Fufem
} // namespace Dune } // namespace Dune
template<class GridView, class LocalFiniteElement, int dim, class field_type=double> #endif //#ifndef DUNE_GFE_HENCKYENERGY_HH
class HenckyEnergy
: public LocalGeodesicFEStiffness<GridView,LocalFiniteElement,RealTuple<field_type,dim> >
{
// grid types
typedef typename GridView::Grid::ctype DT;
typedef RealTuple<field_type,dim> TargetSpace;
typedef typename TargetSpace::ctype RT;
typedef typename GridView::template Codim<0>::Entity Entity;
// some other sizes
enum {gridDim=GridView::dimension};
public: // for testing
/** \brief Constructor with a set of material parameters
* \param parameters The material parameters
*/
HenckyEnergy(const Dune::ParameterTree& parameters,
const BoundaryPatch<GridView>* neumannBoundary,
const Dune::VirtualFunction<Dune::FieldVector<double,gridDim>, Dune::FieldVector<double,3> >* neumannFunction)
: neumannBoundary_(neumannBoundary),
neumannFunction_(neumannFunction)
{
// Lame constants
mu_ = parameters.template get<double>("mu");
lambda_ = parameters.template get<double>("lambda");
}
/** \brief Assemble the energy for a single element */
RT energy (const Entity& e,
const LocalFiniteElement& localFiniteElement,
const std::vector<TargetSpace>& localSolution) const;
/** \brief Lame constants */
double mu_, lambda_;
/** \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, class field_type>
typename HenckyEnergy<GridView,LocalFiniteElement,dim,field_type>::RT
HenckyEnergy<GridView,LocalFiniteElement,dim,field_type>::
energy(const Entity& element,
const LocalFiniteElement& localFiniteElement,
const std::vector<RealTuple<field_type,dim> >& localConfiguration) const
{
assert(element.type() == localFiniteElement.type());
typedef typename GridView::template Codim<0>::Entity::Geometry Geometry;
RT energy = 0;
typedef LocalGeodesicFEFunction<gridDim, DT, LocalFiniteElement, TargetSpace> LocalGFEFunctionType;
LocalGFEFunctionType localGeodesicFEFunction(localFiniteElement,localConfiguration);
int quadOrder = (element.type().isSimplex()) ? localFiniteElement.localBasis().order()
: localFiniteElement.localBasis().order() * gridDim;
const Dune::QuadratureRule<DT, gridDim>& quad
= Dune::QuadratureRules<DT, gridDim>::rule(element.type(), quadOrder);
for (size_t pt=0; pt<quad.size(); pt++) {
// Local position of the quadrature point
const Dune::FieldVector<DT,gridDim>& quadPos = quad[pt].position();
const DT integrationElement = element.geometry().integrationElement(quadPos);
const typename Geometry::JacobianInverseTransposed& jacobianInverseTransposed = element.geometry().jacobianInverseTransposed(quadPos);
DT weight = quad[pt].weight() * integrationElement;
// The derivative of the local function defined on the reference element
typename LocalGFEFunctionType::DerivativeType referenceDerivative = localGeodesicFEFunction.evaluateDerivative(quadPos);
// The derivative of the function defined on the actual element
typename LocalGFEFunctionType::DerivativeType derivative(0);
for (size_t comp=0; comp<referenceDerivative.N(); comp++)
jacobianInverseTransposed.umv(referenceDerivative[comp], derivative[comp]);
/////////////////////////////////////////////////////////
// compute F^T F
/////////////////////////////////////////////////////////
Dune::FieldMatrix<field_type,dim,dim> FTF(0);
for (int i=0; i<dim; i++)
for (int j=0; j<dim; j++)
for (int k=0; k<dim; k++)
FTF[i][j] += derivative[k][i] * derivative[k][j];
//////////////////////////////////////////////////////////
// Eigenvalues of FTF
//////////////////////////////////////////////////////////
#if 0 // hencky
std::array<field_type,dim> lambda = eigenValues(FTF);
//////////////////////////////////////////////////////////
// Compute the derivative of the rotation
// Note: we need it in matrix coordinates
//////////////////////////////////////////////////////////
// logarithms of the eigenvalues
std::array<field_type,dim> ln;
for (int i=0; i<dim; i++)
ln[i] = std::log(lambda[i]);
// Add the local energy density
for (int i=0; i<dim; i++)
energy += weight * mu_ * ln[i]*ln[i];
field_type trace = 0;
for (int i=0; i<dim; i++)
trace += ln[i];
energy += weight * 0.5 * lambda_ * trace * trace;
#else // St.Venant-Kirchhoff
Dune::FieldMatrix<field_type,dim,dim> E = FTF;
for (int i=0; i<dim; i++)
E[i][i] -= 1.0;
E *= 0.5;
field_type trE = E[0][0] + E[1][1] + E[2][2];
Dune::FieldMatrix<field_type,dim,dim> ESquared = E*E;
field_type trESquared = ESquared[0][0] + ESquared[1][1] + ESquared[2][2];
energy += weight * mu_ * trESquared + weight * 0.5 * lambda_ * trE * trE;
#endif
}
//////////////////////////////////////////////////////////////////////////////
// Assemble boundary contributions
//////////////////////////////////////////////////////////////////////////////
if (not neumannFunction_)
return energy;
for (typename Entity::LeafIntersectionIterator it = element.ileafbegin(); it != element.ileafend(); ++it) {
if (not neumannBoundary_ or not neumannBoundary_->contains(*it))
continue;
const Dune::QuadratureRule<DT, gridDim-1>& quad
= Dune::QuadratureRules<DT, gridDim-1>::rule(it->type(), quadOrder);
for (size_t pt=0; pt<quad.size(); pt++) {
// Local position of the quadrature point
const Dune::FieldVector<DT,gridDim>& quadPos = it->geometryInInside().global(quad[pt].position());
const DT integrationElement = it->geometry().integrationElement(quad[pt].position());
// The value of the local function
RealTuple<field_type,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);
// Only translational dofs are affected by the Neumann force
for (size_t i=0; i<neumannValue.size(); i++)
energy += (neumannValue[i] * value.globalCoordinates()[i]) * quad[pt].weight() * integrationElement;
}
}
return energy;
}
#endif //#ifndef COSSERAT_ENERGY_LOCAL_STIFFNESS_HH
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