From ce6d8778d3af35c78ce3603f05efc1866b2697ad Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Mon, 20 Apr 2009 09:54:51 +0000
Subject: [PATCH] complete implementation, not tested beyond compilation yet,
though
[[Imported from SVN: r4059]]
---
src/harmonicenergystiffness.hh | 44 +++++++++++++++++++++-------------
1 file changed, 28 insertions(+), 16 deletions(-)
diff --git a/src/harmonicenergystiffness.hh b/src/harmonicenergystiffness.hh
index 4ae83e62..42ea57fd 100644
--- a/src/harmonicenergystiffness.hh
+++ b/src/harmonicenergystiffness.hh
@@ -1,18 +1,19 @@
#ifndef HARMONIC_ENERGY_LOCAL_STIFFNESS_HH
#define HARMONIC_ENERGY_LOCAL_STIFFNESS_HH
-#include <dune/istl/bcrsmatrix.hh>
#include <dune/common/fmatrix.hh>
-#include <dune/istl/matrixindexset.hh>
-#include <dune/istl/matrix.hh>
-#include <dune/disc/operators/localstiffness.hh>
+#include <dune/grid/common/quadraturerules.hh>
+
+#include "localgeodesicfestiffness.hh"
+#include "localgeodesicfefunction.hh"
template<class GridView, class TargetSpace>
class HarmonicEnergyLocalStiffness
- : public Dune::LocalGeodesicFEStiffness<GridView,TargetSpace>
+ : public LocalGeodesicFEStiffness<GridView,TargetSpace>
{
// grid types
typedef typename GridView::Grid::ctype DT;
+ typedef typename TargetSpace::ctype RT;
typedef typename GridView::template Codim<0>::Entity Entity;
// some other sizes
@@ -23,14 +24,18 @@ public:
//! Dimension of a tangent space
enum { blocksize = TargetSpace::TangentVector::size };
+#if 0
// types for matrics, vectors and boundary conditions
typedef Dune::FieldMatrix<RT,m,m> MBlockType; // one entry in the stiffness matrix
typedef Dune::FieldVector<RT,m> VBlockType; // one entry in the global vectors
typedef Dune::array<Dune::BoundaryConditions::Flags,m> BCBlockType; // componentwise boundary conditions
+#endif
+#if 0
//! Default Constructor
HarmonicEnergyLocalStiffness ()
{}
+#endif
/** \brief Assemble the energy for a single element */
RT energy (const Entity& e,
@@ -38,19 +43,19 @@ public:
};
-template <class GridType, class RT>
-RT RodLocalStiffness<GridType, RT>::
+template <class GridView, class TargetSpace>
+typename HarmonicEnergyLocalStiffness<GridView, TargetSpace>::RT HarmonicEnergyLocalStiffness<GridView, TargetSpace>::
energy(const Entity& element,
- const Dune::array<RigidBodyMotion<3>,2>& localSolution) const
+ const std::vector<TargetSpace>& localSolution) const
{
RT energy = 0;
- assert(element.type().isSimplex());
+ LocalGeodesicFEFunction<gridDim, double, TargetSpace> localGeodesicFEFunction(element.type(), localSolution);
int quadOrder = gridDim;
- const Dune::QuadratureRule<double, 1>& quad
- = Dune::QuadratureRules<double, 1>::rule(element.type(), quadOrder);
+ const Dune::QuadratureRule<double, gridDim>& quad
+ = Dune::QuadratureRules<double, gridDim>::rule(element.type(), quadOrder);
for (size_t pt=0; pt<quad.size(); pt++) {
@@ -58,17 +63,24 @@ energy(const Entity& element,
const Dune::FieldVector<double,gridDim>& quadPos = quad[pt].position();
const double integrationElement = element.geometry().integrationElement(quadPos);
+
+ const Dune::FieldMatrix<double,gridDim,gridDim>& jacobianInverseTransposed = element.geometry().jacobianInverseTransposed(quadPos);
double weight = quad[pt].weight() * integrationElement;
- for (int comp=0; comp<4; comp++) {
+ // The derivative of the local function defined on the reference element
+ Dune::FieldMatrix<double, TargetSpace::EmbeddedTangentVector::size, gridDim> referenceDerivative = localGeodesicFEFunction.evaluateDerivative(quadPos);
+
+ // The derivative of the function defined on the actual element
+ Dune::FieldMatrix<double, TargetSpace::EmbeddedTangentVector::size, gridDim> derivative(0);
- for (int dir=0; dir<gridDim; dir++) {
+ for (int comp=0; comp<4; comp++)
+ jacobianInverseTransposed.umv(referenceDerivative[comp], derivative[comp]);
- double derivative = localGeodesicFunction.evaluateDerivative(comp, dir);
- energy += weight * derivative * derivative;
+ for (int comp=0; comp<4; comp++) {
- }
+ for (int dir=0; dir<gridDim; dir++)
+ energy += weight * derivative[comp][dir] * derivative[comp][dir];
}
--
GitLab