elasticityhelpers.hh 9.71 KB
Newer Older
1
2
3
#ifndef DUNE_ELASTICITY_HELPERS_HH
#define DUNE_ELASTICITY_HELPERS_HH

Oliver Sander's avatar
Oliver Sander committed
4
5
6
7
#include <cmath>

#include <dune/common/fmatrix.hh>

8
9
#include <dune/fufem/symmetrictensor.hh>

10
11
12
13
namespace Dune {

    namespace Elasticity {

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
        /** \brief Compute nonlinear strain tensor from a displacement gradient.
         *
         *  \param grad The gradient of the direction in which the linearisation is computed.
         *  \param strain The tensor to store the strain in. 
         */
        template <int dim>
        void computeNonlinearStrain(const Dune::FieldMatrix<double, dim, dim>& grad, SymmetricTensor<dim>& strain) {
            strain = 0;
            for (int i=0; i<dim ; ++i) {
                strain(i,i) +=grad[i][i];

                for (int k=0;k<dim;++k)
                    strain(i,i) += 0.5*grad[k][i]*grad[k][i];

                for (int j=i+1; j<dim; ++j) {
                    strain(i,j) += 0.5*(grad[i][j] + grad[j][i]);
                    for (int k=0;k<dim;++k)
                        strain(i,j) += 0.5*grad[k][i]*grad[k][j];
                }
            }
        }

        /** \brief Compute linearised strain at the identity from a given displacement gradient. 
         *
         *  \param grad The gradient of the direction in which the linearisation is computed.
         *  \param strain The tensor to store the strain in. 
         */
        template <int dim>
        void computeLinearisedStrain(const Dune::FieldMatrix<double, dim, dim>& grad, SymmetricTensor<dim>& strain) {
            for (int i=0; i<dim ; ++i)
            {
                strain(i,i) = grad[i][i];
                for (int j=i+1; j<dim; ++j)
                    strain(i,j) = 0.5*(grad[i][j] + grad[j][i]);
            }
        }


        /** \brief Compute linearised strain at a prescribed configuration in the direction of a given displacement gradient. 
         *
         *  \param grad The gradient of the direction in which the linearisation is computed.
         *  \param conf The deformation gradient(!) of the configuration at which the linearisation is evaluated.
         *  \param strain The tensor to store the strain in. 
         */
        template <int dim>
        void computeLinearisedStrain(const Dune::FieldMatrix<double, dim, dim>& grad, const Dune::FieldMatrix<double, dim, dim>& conf,
                                     SymmetricTensor<dim>& strain) {
            strain = 0;
            for (int i=0;i<dim; i++)
                for (int j=i;j<dim;j++)
                    for (int k=0;k<dim;k++)
                        strain(i,j)+= grad[k][i]*conf[k][j]+conf[k][i]*grad[k][j];
            strain *= 0.5;
        }

69
70
        /** \brief The penalty term \f$ \Gamma \f$ */
        double Gamma(double x) {
Oliver Sander's avatar
Oliver Sander committed
71
            return x*x - std::log(x);
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
        }
        
        /** \brief First derivative of the penalty term \f$ \Gamma \f$ */
        double Gamma_x(double x) {
            return 2*x - 1/x;
        }
        
        /** \brief Second derivative of the penalty term \f$ \Gamma \f$ */
        double Gamma_xx(double x) {
            return 2 + 1/(x*x);
        }
        
        double delta(int i, int j) { return i==j? 1.0: 0.0; }
        
        template <int dim>
        double E_du(const Dune::FieldMatrix<double,dim,dim>& u, int i, int j, int p, int q) 
        {
            double result = delta(i,p)*delta(j,q) + delta(j,p)*delta(i,q);
            
            result += ( delta(i,q)*u[p][j] + u[p][i]*delta(j,q) );
            return 0.5*result;
        }
        
        // template <int i, int j, int p, int q, int r, int s>
        // static inline double e_dudu() 
        // {
        //     return 0.5 * Delta<p,r>::delta
        //         * ( Delta<i,q>::delta*Delta<j,s>::delta + Delta<i,s>::delta*Delta<j,q>::delta );
        // }
        
        static inline double e_dudu(int i, int j, int p, int q, int r, int s) 
        {
            return 0.5 * delta(p,r)* ( delta(i,q)*delta(j,s) + delta(i,s)*delta(j,q) );
        }
        
        // The Kronecker delta function
        template <int i, int j>
        struct Delta {
            static int const delta = 0;
        };
        
        template <int i>
        struct Delta<i,i> {
            static int const delta = 1;
        };
        
        // \parder {det I + u}{ u_pq }
        template <int p, int qq>
        double det_du_tmpl(const Dune::FieldMatrix<double,3,3>& u) {
            int const q = qq+3;
            
            return (Delta<(p+1)%3,(q+1)%3>::delta + u[(p+1)%3][(q+1)%3])
                * (Delta<(p+2)%3,(q+2)%3>::delta+u[(p+2)%3][(q+2)%3]) 
                - (Delta<(p+1)%3,(q-1)%3>::delta+u[(p+1)%3][(q-1)%3])
                * (Delta<(p+2)%3,(q-2)%3>::delta+u[(p+2)%3][(q-2)%3]);
        }
128
         
akbib's avatar
akbib committed
129
        /** \brief Compute linearization of the determinant of a deformation gradient. Multiplication with a test function(scalar) gradient, gives the linearization in the direction of the test function.
130
131
132
         *
         *  \param u The displacement gradient(!) at which the determinant is evaluated
         */
akbib's avatar
akbib committed
133
        void linearisedDefDet(const Dune::FieldMatrix<double,3,3>& u, Dune::FieldMatrix<double,3,3>& linDet) {
134
135
136
137

            linDet = 0;
            for (int i=0; i<2; i++)                       
                for (int j=(i+1)%3; j<3; j++) {
138
                    int k=(-(i+j))%3;
139
140
141
142
143
144
145
146
                    linDet[i][j] = u[j][k]*u[k][i] - u[j][i]*(1+u[k][k]);
                    linDet[j][i] = u[k][j]*u[i][k] - u[i][j]*(1+u[k][k]);
            }    

            // the diagonal parts
            for (int i=0; i<3; i++)
                linDet[i][i] = (u[(i+1)%3][(i+1)%3]+1)*(u[(i+2)%3][(i+2)%3]+1)-u[(i+1)%3][(i+2)%3]*u[(i+2)%3][(i+1)%3]
        }
147
        
akbib's avatar
akbib committed
148
        /** \brief Compute linearization of the determinant of a deformation gradient. Multiplication with a test function(scalar) gradient, gives the linearization in the direction of the test function.
149
150
151
         *
         *  \param u The displacement gradient(!) at which the determinant is evaluated
         */
akbib's avatar
akbib committed
152
        void linearisedDefDet(const Dune::FieldMatrix<double,2,2>& u, Dune::FieldMatrix<double,2,2>& linDet) {
153
154
155
156
157
158
159
160

            linDet = 0;
            for (int i=0; i<2; i++) {
                linDet[i][i] = 1+u[(i+1)%2][(i+1)%2];
                linDet[i][(i+1)%2] = -u[(i+1)%2][i];
            }
        }

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
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
        // \parder {det I + u}{ u_pq }
        template <int p, int q>
        double det_du_tmpl(const Dune::FieldMatrix<double,2,2>& u) {
            return Delta<p,q>::delta * (1+u[(p+1)%2][(q+1)%2]) 
                - Delta<p,(q+1)%2>::delta * u[(p+1)%2][(q+1)%2];
            
            // p = 0, q = 0;  (1+u[1][1]) 
            // p = 0, q = 1; -u[1][0] 
            // p = 1, q = 0; -u[0][1]
            // p = 1, q = 1;  (1+u[0][0])
            //
        }
        
        // template <int p, int qq, int r, int s>
        // double det_dudu_tmpl(const Dune::FieldMatrix<double,3,3>& u) {
        //     int const q = qq+3;
        
        //     return
        //         Delta<(p+1)%3,r>::delta*Delta<(q+1)%3,s>::delta * (Delta<(p+2)%3,(q+2)%3>::delta+u[(p+2)%3][(q+2)%3])
        //         + Delta<(p+2)%3,r>::delta*Delta<(q+2)%3,s>::delta * (Delta<(p+1)%3,(q+1)%3>::delta+u[(p+1)%3][(q+1)%3])
        //         - Delta<(p+1)%3,r>::delta*Delta<(q-1)%3,s>::delta * (Delta<(p+2)%3,(q-2)%3>::delta+u[(p+2)%3][(q-2)%3])
        //         - Delta<(p+2)%3,r>::delta*Delta<(q-2)%3,s>::delta * (Delta<(p+1)%3,(q-1)%3>::delta+u[(p+1)%3][(q-1)%3]);
        // }
        
        double det_dudu_tmpl(const Dune::FieldMatrix<double,3,3>& u, int p, int qq, int r, int s) {
            int const q = qq+3;
            
            return
                delta((p+1)%3,r) * delta((q+1)%3,s) * (delta((p+2)%3,(q+2)%3)+u[(p+2)%3][(q+2)%3])
                + delta((p+2)%3,r)*delta((q+2)%3,s) * (delta((p+1)%3,(q+1)%3)+u[(p+1)%3][(q+1)%3])
                - delta((p+1)%3,r)*delta((q-1)%3,s) * (delta((p+2)%3,(q-2)%3)+u[(p+2)%3][(q-2)%3])
                - delta((p+2)%3,r)*delta((q-2)%3,s) * (delta((p+1)%3,(q-1)%3)+u[(p+1)%3][(q-1)%3]);
        }
        
        // parder det (I + u) partial u_pq \partial u_rs
        //template <int p, int q, int r, int s>
        double det_dudu_tmpl(const Dune::FieldMatrix<double,2,2>& u, int p, int q, int r, int s) {
            return delta(p,q)*delta(r,(p+1)%2)*delta(s,(q+1)%2)
                - delta(p,(q+1)%2)*delta(r,(p+1)%2)*delta(s,(q+1)%2);
            
            // 0 0 0 0   -->  0  - 0
            // 0 0 0 1   -->  0  - 0
            // 0 0 1 0   -->  0  - 0
            // 0 0 1 1   -->  1  - 0  = 1
            // 0 1 0 0   -->  0  - 0
            // 0 1 0 1   -->  0  - 0
            // 0 1 1 0   -->  0  - 1  = -1
            // 0 1 1 1   -->  0  - 0
            // 1 0 0 0   -->  0  - 0
            // 1 0 0 1   -->  0  - 1  = -1
            // 1 0 1 0   -->  0  - 0
            // 1 0 1 1   -->  0  - 0
            // 1 1 0 0   -->  1  - 0  = 1
            // 1 1 0 1   -->  0  - 0
            // 1 1 1 0   -->  0  - 0
            // 1 1 1 1   -->  0  - 0
        }
218

219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
        double det_val(const Dune::FieldMatrix<double,3,3>& u) {
            return (1+u[0][0])*(1+u[1][1])*(1+u[2][2])
                + u[0][1]*u[1][2]*u[2][0]
                + u[0][2]*u[1][0]*u[2][1]
                - u[2][0]*(1+u[1][1])*u[0][2]
                - u[2][1]*u[1][2]*(1+u[0][0])
                - (1+u[2][2])*u[1][0]*u[0][1];
        }
        
        double det_val(const Dune::FieldMatrix<double,2,2>& u) {
            return (1+u[0][0])*(1+u[1][1]) - u[1][0]*u[0][1];
        }
        
        template <int dim>
        double E_val(const Dune::FieldMatrix<double,dim,dim>& u, int i, int j) 
        {
            double result = u[i][j] + u[j][i];
            
            for (int k=0; k<dim; ++k)
                result += u[k][i]*u[k][j];  
            return 0.5*result;
        }
        
        

    }

}

#endif