neohookeoperatorassembler.hh 8.36 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
#ifndef NEO_HOOKE_OPERATOR_ASSEMBLER_HH
#define NEO_HOOKE_OPERATOR_ASSEMBLER_HH


#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>

#include <dune/istl/matrix.hh>

#include <dune/fufem/quadraturerules/quadraturerulecache.hh>
#include <dune/fufem/assemblers/localoperatorassembler.hh>
#include <dune/fufem/staticmatrixtools.hh>
#include <dune/fufem/symmetrictensor.hh>
#include <dune/fufem/functions/virtualgridfunction.hh>

#include <dune/elasticity/common/elasticityhelpers.hh>

/** \brief Assemble the second derivative of a neo-hookean material.
 *
 *  Laursen:
 *    \f[
 *      W(u)= \frac{\lambda}4(J^2-1)-(\frac{\lambda}2+\mu)log(J)+\mu tr(E(u))
 *    \f]
 *
 *  Fischer/Wriggers:
 *    \f[
 *      W(u)= \frac{\lambda}2(J-1)^2-2\mu)log(J)+\mu tr(E(u))
 *    \f]
 *
 *
 *      - \f$D\varphi\f$: the deformation gradient
 *      - \f$E\f$: the nonlinear strain tensor
 *      - \f$E'_u(v)\f$: the linearized strain tensor at the point u in direction v
 *      - \f$\sigma(u) = C:E(u)\f$: the second Piola-Kirchhoff stress tensor
 *      - \f$C\f$: the Hooke tensor
 *      - \f$Dv\f$: the gradient of the test function
 *
 */
template <class GridType, class TrialLocalFE, class AnsatzLocalFE>
class NeoHookeOperatorAssembler 
    : public LocalOperatorAssembler < GridType, TrialLocalFE, AnsatzLocalFE, Dune::FieldMatrix<typename GridType::ctype ,GridType::dimension,GridType::dimension> >
{


    static const int dim = GridType::dimension;
    typedef typename GridType::ctype ctype;
    
    typedef VirtualGridFunction<GridType, Dune::FieldVector<ctype,GridType::dimension> > GridFunction;

public:
    typedef typename Dune::FieldMatrix<ctype,GridType::dimension,GridType::dimension> T;

    typedef typename LocalOperatorAssembler < GridType, TrialLocalFE, AnsatzLocalFE, T >::Element Element;
    typedef typename LocalOperatorAssembler < GridType, TrialLocalFE, AnsatzLocalFE,T >::BoolMatrix BoolMatrix;
    typedef typename LocalOperatorAssembler < GridType, TrialLocalFE, AnsatzLocalFE,T >::LocalMatrix LocalMatrix;


    NeoHookeOperatorAssembler(ctype lambda, ctype mu, const Dune::shared_ptr<GridFunction> displacement):
        lambda_(lambda), mu_(mu), displacement_(displacement)
    {}

    NeoHookeOperatorAssembler(ctype lambda, ctype mu):
        lambda_(lambda), mu_(mu)
    {}

    void indices(const Element& element, BoolMatrix& isNonZero, const TrialLocalFE& tFE, const AnsatzLocalFE& aFE) const
    {
        isNonZero = true;
    }

    void assemble(const Element& element, LocalMatrix& localMatrix, const TrialLocalFE& tFE, const AnsatzLocalFE& aFE) const
    {
        typedef Dune::FieldVector<ctype,dim> FVdim;
        typedef Dune::FieldMatrix<ctype,dim,dim> FMdimdim;
        typedef typename TrialLocalFE::Traits::LocalBasisType::Traits::JacobianType JacobianType;

        int rows = tFE.localBasis().size();
        int cols = aFE.localBasis().size();

        localMatrix.setSize(rows,cols);
        localMatrix = 0.0;

        // get quadrature rule
        const int order = (element.type().isSimplex()) 
            ? 4*(tFE.localBasis().order())
            : 4*(tFE.localBasis().order()*dim);
        const Dune::template QuadratureRule<ctype, dim>& quad = QuadratureRuleCache<ctype, dim>::rule(element.type(), order, IsRefinedLocalFiniteElement<TrialLocalFE>::value(tFE) );

        // store gradients of shape functions and base functions
        std::vector<JacobianType> referenceGradients(tFE.localBasis().size());
        std::vector<FVdim> gradients(tFE.localBasis().size());

        // store entity geometry mapping
        const typename Element::Geometry geometry = element.geometry();

        // loop over quadrature points
        for (size_t pt=0; pt < quad.size(); ++pt) {

            // get quadrature point
            const FVdim& quadPos = quad[pt].position();

            // get transposed inverse of Jacobian of transformation
            const FMdimdim& invJacobian = geometry.jacobianInverseTransposed(quadPos);

            // get integration factor
            const ctype integrationElement = geometry.integrationElement(quadPos);

            // get gradients of shape functions
            tFE.localBasis().evaluateJacobian(quadPos, referenceGradients);

            // compute gradients of base functions
            for (size_t i=0; i<gradients.size(); ++i) {

                // transform gradients
                gradients[i] = 0.0;
                invJacobian.umv(referenceGradients[i][0], gradients[i]);

            }

            // evaluate displacement gradients the configuration at the quadrature point
            typename GridFunction::DerivativeType localConfGrad;
            if (displacement_->isDefinedOn(element))
                displacement_->evaluateDerivativeLocal(element, quadPos, localConfGrad);
            else
                displacement_->evaluateDerivative(geometry.global(quadPos),localConfGrad);

            // compute linearization of the determinante of the deformation gradient
            FMdimdim linDefDet;
            Dune::Elasticity::linearisedDefDet(localConfGrad,linDefDet);

            // compute deformation gradient of the configuration
            for (int i=0;i<dim;i++)
                localConfGrad[i][i] += 1;                

            // evaluate the derminante of the deformation gradient
            const ctype J = localConfGrad.determinant();

            // /////////////////////////////////////////////////
            //   Assemble matrix
            // /////////////////////////////////////////////////
            ctype z = quad[pt].weight() * integrationElement;

            ctype coeff,coeff2;

#ifdef LAURSEN
            coeff= (0.5*lambda_*J-(lambda_*0.5+mu_)/J);
            coeff2 = (0.5*lambda_ + (lambda_*0.5+mu_)/(J*J));
#else
            coeff = lambda_*(J-1) -2*mu_/J;
            coeff2 = lambda_ + 2*mu_/(J*J);
#endif

            for (int row=0; row<rows; row++)
                    for (int col=0; col<cols; col++) {
                        
                        // second derivative of the determinat term 
                        localMatrix[row][col].axpy(z*coeff,hessDefDet(localConfGrad,gradients[row],gradients[col]));
                      
                        FVdim rowTerm(0),colTerm(0);
                        linDefDet.mv(gradients[row],rowTerm); 
                        linDefDet.mv(gradients[col],colTerm); 

                        // derivative of the coefficient term
                        for (int rcomp=0; rcomp<dim; rcomp++)
                            localMatrix[row][col][rcomp].axpy(coeff2*z*rowTerm[rcomp],colTerm);
 
                        // second derivative of the trace term
                        StaticMatrix::addToDiagonal(localMatrix[row][col],z*mu_*((gradients[row]*gradients[col])));

            }
        }

    }

    /** \brief Set new configuration to be assembled at. Needed e.g. during Newton iteration.*/
    void setConfiguration(Dune::shared_ptr<GridFunction> newConfig) {
        displacement_ = newConfig;
    }

private:
    /** \brief Lame constant */
    ctype lambda_;

    /** \brief Lame constant */
    ctype mu_;

    /** \brief The configuration at which the linearized operator is evaluated.*/
    Dune::shared_ptr<GridFunction> displacement_;


    /** \brief Compute the matrix entry that arises from the second derivative of the determinant. 
     *  Note that it is assumed the defGrad is the deformation gradient and not the displacement gradient   
     */
    Dune::FieldMatrix<ctype,dim,dim> hessDefDet(const Dune::FieldMatrix<ctype,dim,dim>& defGrad,
akbib's avatar
akbib committed
195
                            const Dune::FieldVector<ctype,dim>& testGrad, const Dune::FieldVector<ctype,dim>& ansatzGrad) const
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
    {
        Dune::FieldMatrix<ctype, dim, dim> res(0);
        
        if (dim==2) {
            res[0][1] = testGrad[0]*ansatzGrad[1]-ansatzGrad[0]*testGrad[1];
            res[1][0] = -res[0][1];
            return res;
        }

        //compute the cross product
        Dune::FieldVector<ctype,dim> cross(0);    
        for (int k=0; k<dim; k++)
            cross[k]=testGrad[(k+1)%dim]*ansatzGrad[(k+2)%dim] -testGrad[(k+2)%dim]*ansatzGrad[(k+1)%dim];

        // the resulting matrix is skew symmetric with entries <cross,degGrad[i]>
        for (int i=0; i<dim; i++)
            for (int j=i+1; j<dim; j++) {
akbib's avatar
akbib committed
213
                int k= (dim-(i+j))%dim;
214
215
216
217
218
219
220
221
222
223
224
            res[i][j] = (cross*defGrad[k])*std::pow(-1,k);
            res[j][i] = -res[i][j];
        }
        
        return res;
    }
};


#endif