From 51efe788e48a26f116d6c9f82e577aea3b1b4d3f Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Fri, 12 Jul 2013 21:36:46 +0000
Subject: [PATCH] remove trailing whitespace
[[Imported from SVN: r9308]]
---
dune/gfe/cosseratenergystiffness.hh | 306 ++++++++++++++--------------
1 file changed, 153 insertions(+), 153 deletions(-)
diff --git a/dune/gfe/cosseratenergystiffness.hh b/dune/gfe/cosseratenergystiffness.hh
index 4a216dce..3157e8fb 100644
--- a/dune/gfe/cosseratenergystiffness.hh
+++ b/dune/gfe/cosseratenergystiffness.hh
@@ -20,7 +20,7 @@
template<class GridView, class LocalFiniteElement, int dim>
-class CosseratEnergyLocalStiffness
+class CosseratEnergyLocalStiffness
: public LocalGeodesicFEStiffness<GridView,LocalFiniteElement,RigidBodyMotion<double,dim> >
{
// grid types
@@ -28,11 +28,11 @@ class CosseratEnergyLocalStiffness
typedef RigidBodyMotion<double,dim> TargetSpace;
typedef typename TargetSpace::ctype RT;
typedef typename GridView::template Codim<0>::Entity Entity;
-
+
// some other sizes
enum {gridDim=GridView::dimension};
-
-
+
+
/** \brief Compute the symmetric part of a matrix A, i.e. \f$ \frac 12 (A + A^T) \f$ */
static Dune::FieldMatrix<double,dim,dim> sym(const Dune::FieldMatrix<double,dim,dim>& A)
{
@@ -52,7 +52,7 @@ class CosseratEnergyLocalStiffness
result[i][j] = 0.5 * (A[i][j] - A[j][i]);
return result;
}
-
+
/** \brief Return the square of the trace of a matrix */
template <int N>
static double traceSquared(const Dune::FieldMatrix<double,N,N>& A)
@@ -63,30 +63,30 @@ class CosseratEnergyLocalStiffness
return trace*trace;
}
- /** \brief Compute the (row-wise) curl of a matrix R \f$
+ /** \brief Compute the (row-wise) curl of a matrix R \f$
\param DR The partial derivatives of the matrix R
*/
static Dune::FieldMatrix<double,dim,dim> curl(const Tensor3<double,dim,dim,dim>& DR)
{
Dune::FieldMatrix<double,dim,dim> result;
-
+
for (int i=0; i<dim; i++) {
result[i][0] = DR[i][2][1] - DR[i][1][2];
result[i][1] = DR[i][0][2] - DR[i][2][0];
result[i][2] = DR[i][1][0] - DR[i][0][1];
}
-
+
return result;
}
/** \brief Return the adjugate of the strain matrix U
- *
+ *
* Optimized for U, as we know a bit about its structure
*/
static Dune::FieldMatrix<double,dim,dim> adjugateU(const Dune::FieldMatrix<double,dim,dim>& U)
{
Dune::FieldMatrix<double,dim,dim> result;
-
+
result[0][0] = U[1][1];
result[0][1] = -U[0][1];
result[0][2] = 0;
@@ -101,15 +101,15 @@ class CosseratEnergyLocalStiffness
return result;
}
-
-
+
+
public: // for testing
/** \brief Compute the derivative of the rotation (with respect to x), but wrt matrix coordinates
\param value Value of the gfe function at a certain point
\param derivative First derivative of the gfe function wrt x at that point, in quaternion coordinates
\param DR First derivative of the gfe function wrt x at that point, in matrix coordinates
*/
- static void computeDR(const RigidBodyMotion<double,3>& value,
+ static void computeDR(const RigidBodyMotion<double,3>& value,
const Dune::FieldMatrix<double,7,gridDim>& derivative,
Tensor3<double,3,3,3>& DR)
{
@@ -124,7 +124,7 @@ public: // for testing
// So, if I am not mistaken, DR[i][j][k] contains \partial R_ij / \partial k
Tensor3<double,3 , 3, 4> dd_dq;
value.q.getFirstDerivativesOfDirectors(dd_dq);
-
+
DR = 0;
for (int i=0; i<3; i++)
for (int j=0; j<3; j++)
@@ -138,7 +138,7 @@ public: // for testing
\param value Value of the gfe function at a certain point
\param derivative First derivative of the gfe function wrt the coefficients at that point, in quaternion coordinates
*/
- static void computeDR_dv(const RigidBodyMotion<double,3>& value,
+ static void computeDR_dv(const RigidBodyMotion<double,3>& value,
const Dune::FieldMatrix<double,7,7>& derivative,
Tensor3<double,3,3,4>& dR_dv)
{
@@ -153,7 +153,7 @@ public: // for testing
// So, if I am not mistaken, DR[i][j][k] contains \partial R_ij / \partial k
Tensor3<double,3 , 3, 4> dd_dq;
value.q.getFirstDerivativesOfDirectors(dd_dq);
-
+
dR_dv = 0;
for (int i=0; i<3; i++)
for (int j=0; j<3; j++)
@@ -165,11 +165,11 @@ public: // for testing
/** \brief Compute the derivative of the gradient of the rotation with respect to the corner coefficients,
* but in matrix coordinates
- *
+ *
\param value Value of the gfe function at a certain point
\param derivative First derivative of the gfe function wrt the coefficients at that point, in quaternion coordinates
*/
- static void compute_dDR_dv(const RigidBodyMotion<double,3>& value,
+ static void compute_dDR_dv(const RigidBodyMotion<double,3>& value,
const Dune::FieldMatrix<double,7,gridDim>& derOfValueWRTx,
const Dune::FieldMatrix<double,7,7>& derOfValueWRTCoefficient,
const Tensor3<double,7,7,gridDim>& derOfGradientWRTCoefficient,
@@ -187,14 +187,14 @@ public: // for testing
// So, if I am not mistaken, DR[i][j][k] contains \partial R_ij / \partial k
Tensor3<double,3 , 3, 4> dd_dq;
value.q.getFirstDerivativesOfDirectors(dd_dq);
-
+
/** \todo This is a constant -- don't recompute it every time */
Dune::array<Tensor3<double,3,4,4>, 3> dd_dq_dq;
Rotation<double,3>::getSecondDerivativesOfDirectors(dd_dq_dq);
-
+
for (size_t i=0; i<dDR_dv.size(); i++)
dDR_dv[i] = 0;
-
+
for (int i=0; i<3; i++)
for (int j=0; j<3; j++)
for (int k=0; k<gridDim; k++)
@@ -214,10 +214,10 @@ public: // for testing
// Compute covariant derivative for the third director
// by projecting onto the tangent space at S^2. Yes, that is something different!
///////////////////////////////////////////////////////////////////////////////////////////
-
+
Dune::FieldMatrix<double,3,3> Mtmp;
value.q.matrix(Mtmp);
- OrthogonalMatrix<double,3> M(Mtmp);
+ OrthogonalMatrix<double,3> M(Mtmp);
Dune::FieldVector<double,3> tmpDirector;
for (int i=0; i<3; i++)
@@ -226,13 +226,13 @@ public: // for testing
for (int k=0; k<gridDim; k++)
for (int v_i=0; v_i<4; v_i++) {
-
+
Dune::FieldVector<double,3> unprojected;
for (int i=0; i<3; i++)
unprojected[i] = dDR_dv[i][2][k][v_i];
-
+
Dune::FieldVector<double,3> projected = director.projectOntoTangentSpace(unprojected);
-
+
for (int i=0; i<3; i++)
dDR3_dv[i][k][v_i] = projected[i];
}
@@ -244,23 +244,23 @@ public: // for testing
for (int k=0; k<gridDim; k++)
for (int v_i=0; v_i<4; v_i++) {
-
+
Dune::FieldMatrix<double,3,3> unprojected;
for (int i=0; i<3; i++)
for (int j=0; j<3; j++)
unprojected[i][j] = dDR_dv[i][j][k][v_i];
-
+
Dune::FieldMatrix<double,3,3> projected = M.projectOntoTangentSpace(unprojected);
-
+
for (int i=0; i<3; i++)
for (int j=0; j<3; j++)
dDR_dv[i][j][k][v_i] = projected[i][j];
}
-
+
}
public:
-
+
/** \brief Constructor with a set of material parameters
* \param parameters The material parameters
*/
@@ -272,20 +272,20 @@ public:
{
// The shell thickness
thickness_ = parameters.template get<double>("thickness");
-
+
// Lame constants
mu_ = parameters.template get<double>("mu");
lambda_ = parameters.template get<double>("lambda");
// Cosserat couple modulus, preferably 0
mu_c_ = parameters.template get<double>("mu_c");
-
+
// Length scale parameter
L_c_ = parameters.template get<double>("L_c");
-
+
// Curvature exponent
q_ = parameters.template get<double>("q");
-
+
// Shear correction factor
kappa_ = parameters.template get<double>("kappa");
}
@@ -309,7 +309,7 @@ public:
Dune::FieldMatrix<double,3,3> UMinus1 = U;
for (int i=0; i<dim; i++)
UMinus1[i][i] -= 1;
-
+
return mu_ * sym(UMinus1).frobenius_norm2()
+ mu_c_ * skew(UMinus1).frobenius_norm2()
+ (mu_*lambda_)/(2*mu_ + lambda_) * traceSquared(sym(UMinus1));
@@ -321,7 +321,7 @@ public:
RT longQuadraticMembraneEnergy(const Dune::FieldMatrix<double,3,3>& U) const
{
RT result = 0;
-
+
// shear-stretch energy
Dune::FieldMatrix<double,dim-1,dim-1> sym2x2;
for (int i=0; i<dim-1; i++)
@@ -329,7 +329,7 @@ public:
sym2x2[i][j] = 0.5 * (U[i][j] + U[j][i]) - (i==j);
result += mu_ * sym2x2.frobenius_norm2();
-
+
// first order drill energy
Dune::FieldMatrix<double,dim-1,dim-1> skew2x2;
for (int i=0; i<dim-1; i++)
@@ -337,17 +337,17 @@ public:
skew2x2[i][j] = 0.5 * (U[i][j] - U[j][i]);
result += mu_c_ * skew2x2.frobenius_norm2();
-
-
+
+
// classical transverse shear energy
result += kappa_ * (mu_ + mu_c_)/2 * (U[2][0]*U[2][0] + U[2][1]*U[2][1]);
-
+
// elongational stretch energy
result += mu_*lambda_ / (2*mu_ + lambda_) * traceSquared(sym2x2);
-
+
return result;
}
-
+
/** \brief Energy for large-deformation problems (private communication by Patrizio Neff)
*/
RT nonquadraticMembraneEnergy(const Dune::FieldMatrix<double,3,3>& U) const
@@ -355,9 +355,9 @@ public:
Dune::FieldMatrix<double,3,3> UMinus1 = U;
for (int i=0; i<dim; i++)
UMinus1[i][i] -= 1;
-
+
RT detU = U.determinant();
-
+
return mu_ * sym(UMinus1).frobenius_norm2()
+ (mu_*lambda_)/(2*mu_ + lambda_) * 0.5 * ((detU-1)*(detU-1) + (1.0/detU -1)*(1.0/detU -1));
}
@@ -381,9 +381,9 @@ public:
for (int k=0; k<3; k++)
RT_DR3[i][j] += R[k][i] * DR[k][2][j];
}
-
-
-
+
+
+
return mu_ * sym(RT_DR3).frobenius_norm2()
+ mu_c_ * skew(RT_DR3).frobenius_norm2()
+ mu_*lambda_/(2*mu_+lambda_) * traceSquared(RT_DR3);
@@ -410,35 +410,35 @@ public:
const Dune::FieldMatrix<double,3,3>& R,
const Tensor3<double,3,3,3>& DR,
const Dune::array<Tensor3<double,3,3,4>, 3>& dDR_dv) const;
-
+
void bendingEnergyGradient(typename TargetSpace::EmbeddedTangentVector& embeddedLocalGradient,
const Dune::FieldMatrix<double,3,3>& R,
const Tensor3<double,3,3,4>& dR_dv,
const Tensor3<double,3,3,3>& DR,
const Tensor3<double,3,gridDim,4>& dDR3_dv) const;
-
+
/** \brief The shell thickness */
double thickness_;
-
+
/** \brief Lame constants */
double mu_, lambda_;
/** \brief Cosserat couple modulus, preferably 0 */
double mu_c_;
-
+
/** \brief Length scale parameter */
double L_c_;
-
+
/** \brief Curvature exponent */
double q_;
/** \brief Shear correction factor */
double kappa_;
-
+
/** \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_;
};
@@ -458,20 +458,20 @@ energy(const Entity& element,
int quadOrder = (element.type().isSimplex()) ? localFiniteElement.localBasis().order()
: localFiniteElement.localBasis().order() * gridDim;
- const Dune::QuadratureRule<double, gridDim>& quad
+ const Dune::QuadratureRule<double, gridDim>& quad
= Dune::QuadratureRules<double, gridDim>::rule(element.type(), quadOrder);
-
+
for (size_t pt=0; pt<quad.size(); pt++) {
-
+
// Local position of the quadrature point
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;
-
+
// The value of the local function
RigidBodyMotion<double,dim> value = localGeodesicFEFunction.evaluate(quadPos);
@@ -483,16 +483,16 @@ energy(const Entity& element,
for (size_t comp=0; comp<referenceDerivative.N(); comp++)
jacobianInverseTransposed.umv(referenceDerivative[comp], derivative[comp]);
-
+
/////////////////////////////////////////////////////////
// compute U, the Cosserat strain
/////////////////////////////////////////////////////////
dune_static_assert(dim>=gridDim, "Codim of the grid must be nonnegative");
-
+
//
Dune::FieldMatrix<double,dim,dim> R;
value.q.matrix(R);
-
+
/* Compute F, the deformation gradient.
In the continuum case this is
\f$ \hat{F} = \nabla m \f$
@@ -503,11 +503,11 @@ energy(const Entity& element,
for (int i=0; i<dim; i++)
for (int j=0; j<gridDim; j++)
F[i][j] = derivative[i][j];
-
+
for (int i=0; i<dim; i++)
for (int j=gridDim; j<dim; j++)
F[i][j] = R[i][j];
-
+
// U = R^T F
Dune::FieldMatrix<double,dim,dim> U;
for (int i=0; i<dim; i++)
@@ -516,15 +516,15 @@ energy(const Entity& element,
for (int k=0; k<dim; k++)
U[i][j] += R[k][i] * F[k][j];
}
-
+
//////////////////////////////////////////////////////////
// Compute the derivative of the rotation
// Note: we need it in matrix coordinates
//////////////////////////////////////////////////////////
-
+
Tensor3<double,3,3,3> DR;
computeDR(value, derivative, DR);
-
+
// Add the local energy density
if (gridDim==2) {
#ifdef QUADRATIC_MEMBRANE_ENERGY
@@ -542,27 +542,27 @@ 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
+
+ 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
@@ -570,7 +570,7 @@ energy(const Entity& element,
// 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
@@ -578,11 +578,11 @@ energy(const Entity& element,
// Only translational dofs are affected by the Neumann force
energy += thickness_ * (neumannValue * value.r) * quad[pt].weight() * integrationElement;
-
+
}
-
+
}
-
+
return energy;
}
@@ -598,36 +598,36 @@ nonquadraticMembraneEnergyGradient(typename TargetSpace::EmbeddedTangentVector&
{
// Compute partial/partial v ((R_1|R_2)^T\nablam)
Tensor3<double,3,3,7> dU_dv(0);
-
+
for (size_t i=0; i<3; i++) {
for (size_t j=0; j<2; j++) {
-
+
for (size_t k=0; k<3; k++) {
-
+
// Translational dofs of the coefficient
for (size_t v_i=0; v_i<3; v_i++)
dU_dv[i][j][v_i] += R[k][i] * derOfGradientWRTCoefficient[k][v_i][j];
-
+
// Rotational dofs of the coefficient
for (size_t v_i=0; v_i<4; v_i++)
dU_dv[i][j][v_i+3] += dR_dv[k][i][v_i] * derivative[k][j];
-
+
}
-
+
}
-
+
dU_dv[i][2] = 0;
-
+
}
-
-
+
+
////////////////////////////////////////////////////////////////////////////////////
// Quadratic part
////////////////////////////////////////////////////////////////////////////////////
for (size_t v_i=0; v_i<7; v_i++) {
-
+
for (size_t i=0; i<3; i++)
for (size_t j=0; j<3; j++) {
double symUMinusI = 0.5*U[i][j] + 0.5*U[j][i] - (i==j);
@@ -642,12 +642,12 @@ nonquadraticMembraneEnergyGradient(typename TargetSpace::EmbeddedTangentVector&
/////////////////////////////////////////////////////////////////////////////////////
RT detU = U.determinant(); // todo: use structure of U
Dune::FieldMatrix<RT,3,3> adjU = adjugateU(U); // the transpose of the derivative of detU
-
+
for (size_t v_i=0; v_i<7; v_i++)
for (size_t i=0; i<3; i++)
for (size_t j=0; j<3; j++)
embeddedLocalGradient[v_i] += mu_*lambda_/(2*mu_+lambda_) * (detU - 1 - (1.0/detU -1)/(detU*detU)) * adjU[j][i] * dU_dv[i][j][v_i];
-
+
}
@@ -663,34 +663,34 @@ longQuadraticMembraneEnergyGradient(typename TargetSpace::EmbeddedTangentVector&
{
// Compute partial/partial v ((R_1|R_2)^T\nablam)
Tensor3<double,2,2,7> dR1R2Tnablam_dv(0);
-
+
for (size_t i=0; i<2; i++) {
for (size_t j=0; j<2; j++) {
-
+
for (size_t k=0; k<3; k++) {
-
+
// Translational dofs of the coefficient
for (size_t v_i=0; v_i<3; v_i++)
dR1R2Tnablam_dv[i][j][v_i] += R[k][i] * derOfGradientWRTCoefficient[k][v_i][j];
-
+
// Rotational dofs of the coefficient
for (size_t v_i=0; v_i<4; v_i++)
dR1R2Tnablam_dv[i][j][v_i+3] += dR_dv[k][i][v_i] * derivative[k][j];
-
+
}
-
+
}
-
+
}
-
-
+
+
////////////////////////////////////////////////////////////////////////////////////
// Shear--Stretch Energy
////////////////////////////////////////////////////////////////////////////////////
for (size_t v_i=0; v_i<7; v_i++) {
-
+
for (size_t i=0; i<2; i++)
for (size_t j=0; j<2; j++) {
double symUMinusI = 0.5*U[i][j] + 0.5*U[j][i] - (i==j);
@@ -705,33 +705,33 @@ longQuadraticMembraneEnergyGradient(typename TargetSpace::EmbeddedTangentVector&
////////////////////////////////////////////////////////////////////////////////////
for (size_t v_i=0; v_i<7; v_i++) {
-
+
for (size_t i=0; i<2; i++)
for (size_t j=0; j<2; j++)
embeddedLocalGradient[v_i] += mu_c_ * 0.5* (U[i][j]-U[j][i]) * (dR1R2Tnablam_dv[i][j][v_i] - dR1R2Tnablam_dv[j][i][v_i]);
-
+
}
////////////////////////////////////////////////////////////////////////////////////
// Classical Transverse Shear Energy
////////////////////////////////////////////////////////////////////////////////////
-
+
for (size_t v_i=0; v_i<7; v_i++) {
-
+
for (size_t j=0; j<2; j++) {
double sp = 0;
for (size_t k=0; k<3; k++)
sp += R[k][2]*derivative[k][j];
double spDer = 0;
-
+
if (v_i < 3)
for (size_t k=0; k<3; k++)
spDer += R[k][2] * derOfGradientWRTCoefficient[k][v_i][j];
else
for (size_t k=0; k<3; k++)
spDer += dR_dv[k][2][v_i-3] * derivative[k][j];
-
+
embeddedLocalGradient[v_i] += 0.5 * kappa_ * (mu_ + mu_c_) * 2 * sp * spDer;
}
@@ -741,9 +741,9 @@ longQuadraticMembraneEnergyGradient(typename TargetSpace::EmbeddedTangentVector&
////////////////////////////////////////////////////////////////////////////////////
// Elongational Stretch Energy
////////////////////////////////////////////////////////////////////////////////////
-
+
double trace = U[0][0] + U[1][1] - 2;
-
+
for (size_t v_i=0; v_i<7; v_i++)
for (size_t i=0; i<2; i++)
embeddedLocalGradient[v_i] += 2 * mu_*lambda_ / (2*mu_+lambda_)* trace * dR1R2Tnablam_dv[i][i][v_i];
@@ -762,14 +762,14 @@ curvatureEnergyGradient(typename TargetSpace::EmbeddedTangentVector& embeddedLoc
#error curvatureEnergyGradient not implemented for the curl curvature energy
#endif
embeddedLocalGradient = 0;
-
+
for (size_t v_i=0; v_i<4; v_i++) {
-
+
for (size_t i=0; i<3; i++)
for (size_t j=0; j<3; j++)
for (size_t k=0; k<3; k++)
embeddedLocalGradient[v_i+3] += DR[i][j][k] * dDR_dv[i][j][k][v_i];
-
+
}
// multiply with constant pre-factor
@@ -786,7 +786,7 @@ bendingEnergyGradient(typename TargetSpace::EmbeddedTangentVector& embeddedLocal
const Tensor3<double,3,gridDim,4>& dDR3_dv) const
{
embeddedLocalGradient = 0;
-
+
// left-multiply the derivative of the third director (in DR[][2][]) with R^T
Dune::FieldMatrix<double,3,3> RT_DR3(0);
for (int i=0; i<3; i++)
@@ -796,9 +796,9 @@ bendingEnergyGradient(typename TargetSpace::EmbeddedTangentVector& embeddedLocal
for (int i=0; i<gridDim; i++)
assert(std::fabs(RT_DR3[2][i]) < 1e-7);
-
+
// -----------------------------------------------------------------------------
-
+
Tensor3<double,3,3,4> d_RT_DR3(0);
for (size_t v_i=0; v_i<4; v_i++)
for (int i=0; i<3; i++)
@@ -812,7 +812,7 @@ bendingEnergyGradient(typename TargetSpace::EmbeddedTangentVector& embeddedLocal
for (int i=0; i<gridDim; i++)
for (size_t v_i=0; v_i<4; v_i++)
d_RT_DR3[2][i][v_i] = 0;
-
+
////////////////////////////////////////////////////////////////////////////////
// "Shear--Stretch Energy"
////////////////////////////////////////////////////////////////////////////////
@@ -839,7 +839,7 @@ bendingEnergyGradient(typename TargetSpace::EmbeddedTangentVector& embeddedLocal
for (size_t v_i=0; v_i<4; v_i++)
for (int i=0; i<3; i++)
embeddedLocalGradient[v_i+3] += factor * d_RT_DR3[i][i][v_i];
-
+
}
template <class GridView, class LocalFiniteElement, int dim>
@@ -852,7 +852,7 @@ assembleGradient(const Entity& element,
std::vector<typename TargetSpace::EmbeddedTangentVector> embeddedLocalGradient(localSolution.size());
std::fill(embeddedLocalGradient.begin(), embeddedLocalGradient.end(), 0);
-
+
LocalGeodesicFEFunction<gridDim, double, LocalFiniteElement, TargetSpace> localGeodesicFEFunction(localFiniteElement,
localSolution);
@@ -860,20 +860,20 @@ assembleGradient(const Entity& element,
int quadOrder = (element.type().isSimplex()) ? localFiniteElement.localBasis().order()
: localFiniteElement.localBasis().order() * gridDim;
- const Dune::QuadratureRule<double, gridDim>& quad
+ const Dune::QuadratureRule<double, gridDim>& quad
= Dune::QuadratureRules<double, gridDim>::rule(element.type(), quadOrder);
-
+
for (size_t pt=0; pt<quad.size(); pt++) {
-
+
// Local position of the quadrature point
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;
-
+
// The value of the local function
RigidBodyMotion<double,dim> value = localGeodesicFEFunction.evaluate(quadPos);
@@ -885,16 +885,16 @@ assembleGradient(const Entity& element,
for (size_t comp=0; comp<referenceDerivative.N(); comp++)
jacobianInverseTransposed.umv(referenceDerivative[comp], derivative[comp]);
-
+
/////////////////////////////////////////////////////////
// compute U, the Cosserat strain
/////////////////////////////////////////////////////////
dune_static_assert(dim>=gridDim, "Codim of the grid must be nonnegative");
-
+
//
Dune::FieldMatrix<double,dim,dim> R;
value.q.matrix(R);
-
+
/* Compute F, the deformation gradient.
In the continuum case this is
\f$ \hat{F} = \nabla m \f$
@@ -905,11 +905,11 @@ assembleGradient(const Entity& element,
for (int i=0; i<dim; i++)
for (int j=0; j<gridDim; j++)
F[i][j] = derivative[i][j];
-
+
for (int i=0; i<dim; i++)
for (int j=gridDim; j<dim; j++)
F[i][j] = R[i][j];
-
+
// U = R^T F
Dune::FieldMatrix<double,dim,dim> U;
for (int i=0; i<dim; i++)
@@ -918,18 +918,18 @@ assembleGradient(const Entity& element,
for (int k=0; k<dim; k++)
U[i][j] += R[k][i] * F[k][j];
}
-
+
//////////////////////////////////////////////////////////
// Compute the derivative of the rotation
// Note: we need it in matrix coordinates
//////////////////////////////////////////////////////////
-
+
Tensor3<double,3,3,3> DR;
computeDR(value, derivative, DR);
// loop over all the element's degrees of freedom and compute the gradient wrt it
for (size_t i=0; i<localSolution.size(); i++) {
-
+
// --------------------------------------------------------------------
Dune::FieldMatrix<double,7,7> derOfValueWRTCoefficient;
localGeodesicFEFunction.evaluateDerivativeOfValueWRTCoefficient(quadPos, i, derOfValueWRTCoefficient);
@@ -940,7 +940,7 @@ assembleGradient(const Entity& element,
// --------------------------------------------------------------------
Tensor3<double, TargetSpace::EmbeddedTangentVector::dimension,TargetSpace::EmbeddedTangentVector::dimension,gridDim> referenceDerivativeDerivative;
localGeodesicFEFunction.evaluateDerivativeOfGradientWRTCoefficient(quadPos, i, referenceDerivativeDerivative);
-
+
// multiply the transformation from the reference element to the actual element
Tensor3<double, TargetSpace::EmbeddedTangentVector::dimension,TargetSpace::EmbeddedTangentVector::dimension,gridDim> derOfGradientWRTCoefficient;
for (int ii=0; ii<TargetSpace::EmbeddedTangentVector::dimension; ii++)
@@ -955,7 +955,7 @@ assembleGradient(const Entity& element,
Dune::array<Tensor3<double,3,3,4>, 3> dDR_dv;
Tensor3<double,3,gridDim,4> dDR3_dv;
compute_dDR_dv(value, derivative, derOfValueWRTCoefficient, derOfGradientWRTCoefficient, dDR_dv, dDR3_dv);
-
+
// Add the local energy density
if (gridDim==2) {
typename TargetSpace::EmbeddedTangentVector tmp(0);
@@ -965,11 +965,11 @@ assembleGradient(const Entity& element,
nonquadraticMembraneEnergyGradient(tmp,R,dR_dv,derivative,derOfGradientWRTCoefficient,U);
#endif
embeddedLocalGradient[i].axpy(weight * thickness_, tmp);
-
+
tmp = 0;
curvatureEnergyGradient(tmp,R,DR,dDR_dv);
embeddedLocalGradient[i].axpy(weight * thickness_, tmp);
-
+
tmp = 0;
bendingEnergyGradient(tmp,R,dR_dv,DR,dDR3_dv);
embeddedLocalGradient[i].axpy(weight * std::pow(thickness_,3) / 12.0, tmp);
@@ -981,33 +981,33 @@ assembleGradient(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
+
+ 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());
// 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
@@ -1015,7 +1015,7 @@ assembleGradient(const Entity& element,
// loop over all the element's degrees of freedom and compute the gradient wrt it
for (size_t i=0; i<localSolution.size(); i++) {
-
+
// --------------------------------------------------------------------
Dune::FieldMatrix<double,7,7> derOfValueWRTCoefficient;
localGeodesicFEFunction.evaluateDerivativeOfValueWRTCoefficient(quadPos, i, derOfValueWRTCoefficient);
@@ -1025,10 +1025,10 @@ assembleGradient(const Entity& element,
for (size_t j=0; j<3; j++)
#warning Try whether the arguments of derOfValueWRTCoefficient are swapped
embeddedLocalGradient[i][v_i] += thickness_ * (neumannValue[j] * derOfValueWRTCoefficient[j][v_i]) * quad[pt].weight() * integrationElement;
-
+
}
}
-
+
}
--
GitLab