Reference documentation for deal.II version 8.1.0
elasticity.h
1 // ---------------------------------------------------------------------
2 // @f$Id: elasticity.h 30963 2013-09-26 18:42:28Z kanschat @f$
3 //
4 // Copyright (C) 2010 - 2013 by the deal.II authors
5 //
6 // This file is part of the deal.II library.
7 //
8 // The deal.II library is free software; you can use it, redistribute
9 // it, and/or modify it under the terms of the GNU Lesser General
10 // Public License as published by the Free Software Foundation; either
11 // version 2.1 of the License, or (at your option) any later version.
12 // The full text of the license can be found in the file LICENSE at
13 // the top level of the deal.II distribution.
14 //
15 // ---------------------------------------------------------------------
16 
17 #ifndef __deal2__integrators_elasticity_h
18 #define __deal2__integrators_elasticity_h
19 
20 
21 #include <deal.II/base/config.h>
23 #include <deal.II/base/quadrature.h>
24 #include <deal.II/lac/full_matrix.h>
25 #include <deal.II/fe/mapping.h>
26 #include <deal.II/fe/fe_values.h>
27 #include <deal.II/meshworker/dof_info.h>
28 
30 
31 namespace LocalIntegrators
32 {
40  namespace Elasticity
41  {
50  template <int dim>
51  inline void cell_matrix (
53  const FEValuesBase<dim> &fe,
54  const double factor = 1.)
55  {
56  const unsigned int n_dofs = fe.dofs_per_cell;
57 
58  AssertDimension(fe.get_fe().n_components(), dim);
59  AssertDimension(M.m(), n_dofs);
60  AssertDimension(M.n(), n_dofs);
61 
62  for (unsigned int k=0; k<fe.n_quadrature_points; ++k)
63  {
64  const double dx = factor * fe.JxW(k);
65  for (unsigned int i=0; i<n_dofs; ++i)
66  for (unsigned int j=0; j<n_dofs; ++j)
67  for (unsigned int d1=0; d1<dim; ++d1)
68  for (unsigned int d2=0; d2<dim; ++d2)
69  M(i,j) += dx * .25 *
70  (fe.shape_grad_component(j,k,d1)[d2] + fe.shape_grad_component(j,k,d2)[d1]) *
71  (fe.shape_grad_component(i,k,d1)[d2] + fe.shape_grad_component(i,k,d2)[d1]);
72  }
73  }
74 
75 
83  template <int dim, typename number>
84  inline void
86  Vector<number> &result,
87  const FEValuesBase<dim> &fe,
88  const VectorSlice<const std::vector<std::vector<Tensor<1,dim> > > > &input,
89  double factor = 1.)
90  {
91  const unsigned int nq = fe.n_quadrature_points;
92  const unsigned int n_dofs = fe.dofs_per_cell;
93  AssertDimension(fe.get_fe().n_components(), dim);
94 
96  Assert(result.size() == n_dofs, ExcDimensionMismatch(result.size(), n_dofs));
97 
98  for (unsigned int k=0; k<nq; ++k)
99  {
100  const double dx = factor * fe.JxW(k);
101  for (unsigned int i=0; i<n_dofs; ++i)
102  for (unsigned int d1=0; d1<dim; ++d1)
103  for (unsigned int d2=0; d2<dim; ++d2)
104  {
105  result(i) += dx * .25 *
106  (input[d1][k][d2] + input[d2][k][d1]) *
107  (fe.shape_grad_component(i,k,d1)[d2] + fe.shape_grad_component(i,k,d2)[d1]);
108  }
109  }
110  }
111 
112 
120  template <int dim>
121  inline void nitsche_matrix (
123  const FEValuesBase<dim> &fe,
124  double penalty,
125  double factor = 1.)
126  {
127  const unsigned int n_dofs = fe.dofs_per_cell;
128 
129  AssertDimension(fe.get_fe().n_components(), dim);
130  AssertDimension(M.m(), n_dofs);
131  AssertDimension(M.n(), n_dofs);
132 
133  for (unsigned int k=0; k<fe.n_quadrature_points; ++k)
134  {
135  const double dx = factor * fe.JxW(k);
136  const Point<dim> &n = fe.normal_vector(k);
137  for (unsigned int i=0; i<n_dofs; ++i)
138  for (unsigned int j=0; j<n_dofs; ++j)
139  for (unsigned int d1=0; d1<dim; ++d1)
140  {
141  const double u = fe.shape_value_component(j,k,d1);
142  const double v = fe.shape_value_component(i,k,d1);
143  M(i,j) += dx * 2. * penalty * u * v;
144  for (unsigned int d2=0; d2<dim; ++d2)
145  {
146  // v . nabla u n
147  M(i,j) -= .5*dx* fe.shape_grad_component(j,k,d1)[d2] *n(d2)* v;
148  // v (nabla u)^T n
149  M(i,j) -= .5*dx* fe.shape_grad_component(j,k,d2)[d1] *n(d2)* v;
150  // u nabla v n
151  M(i,j) -= .5*dx* fe.shape_grad_component(i,k,d1)[d2] *n(d2)* u;
152  // u (nabla v)^T n
153  M(i,j) -= .5*dx* fe.shape_grad_component(i,k,d2)[d1] *n(d2)* u;
154  }
155  }
156  }
157  }
158 
175  template <int dim, typename number>
177  Vector<number> &result,
178  const FEValuesBase<dim> &fe,
179  const VectorSlice<const std::vector<std::vector<double> > > &input,
180  const VectorSlice<const std::vector<std::vector<Tensor<1,dim> > > > &Dinput,
181  const VectorSlice<const std::vector<std::vector<double> > > &data,
182  double penalty,
183  double factor = 1.)
184  {
185  const unsigned int n_dofs = fe.dofs_per_cell;
189 
190  for (unsigned int k=0; k<fe.n_quadrature_points; ++k)
191  {
192  const double dx = factor * fe.JxW(k);
193  const Point<dim> &n = fe.normal_vector(k);
194  for (unsigned int i=0; i<n_dofs; ++i)
195  for (unsigned int d1=0; d1<dim; ++d1)
196  {
197  const double u= input[d1][k];
198  const double v= fe.shape_value_component(i,k,d1);
199  const double g= data[d1][k];
200  result(i) += dx * 2.*penalty * (u-g) * v;
201 
202  for (unsigned int d2=0; d2<dim; ++d2)
203  {
204  // v . nabla u n
205  result(i) -= .5*dx* v * Dinput[d1][k][d2] * n(d2);
206  // v . (nabla u)^T n
207  result(i) -= .5*dx* v * Dinput[d2][k][d1] * n(d2);
208  // u nabla v n
209  result(i) -= .5*dx * (u-g) * fe.shape_grad_component(i,k,d1)[d2] * n(d2);
210  // u (nabla v)^T n
211  result(i) -= .5*dx * (u-g) * fe.shape_grad_component(i,k,d2)[d1] * n(d2);
212  }
213  }
214  }
215  }
216 
232  template <int dim, typename number>
234  Vector<number> &result,
235  const FEValuesBase<dim> &fe,
236  const VectorSlice<const std::vector<std::vector<double> > > &input,
237  const VectorSlice<const std::vector<std::vector<Tensor<1,dim> > > > &Dinput,
238  double penalty,
239  double factor = 1.)
240  {
241  const unsigned int n_dofs = fe.dofs_per_cell;
244 
245  for (unsigned int k=0; k<fe.n_quadrature_points; ++k)
246  {
247  const double dx = factor * fe.JxW(k);
248  const Point<dim> &n = fe.normal_vector(k);
249  for (unsigned int i=0; i<n_dofs; ++i)
250  for (unsigned int d1=0; d1<dim; ++d1)
251  {
252  const double u= input[d1][k];
253  const double v= fe.shape_value_component(i,k,d1);
254  result(i) += dx * 2.*penalty * u * v;
255 
256  for (unsigned int d2=0; d2<dim; ++d2)
257  {
258  // v . nabla u n
259  result(i) -= .5*dx* v * Dinput[d1][k][d2] * n(d2);
260  // v . (nabla u)^T n
261  result(i) -= .5*dx* v * Dinput[d2][k][d1] * n(d2);
262  // u nabla v n
263  result(i) -= .5*dx * u * fe.shape_grad_component(i,k,d1)[d2] * n(d2);
264  // u (nabla v)^T n
265  result(i) -= .5*dx * u * fe.shape_grad_component(i,k,d2)[d1] * n(d2);
266  }
267  }
268  }
269  }
270 
275  template <int dim>
276  inline void ip_matrix (
277  FullMatrix<double> &M11,
278  FullMatrix<double> &M12,
279  FullMatrix<double> &M21,
280  FullMatrix<double> &M22,
281  const FEValuesBase<dim> &fe1,
282  const FEValuesBase<dim> &fe2,
283  const double pen,
284  const double int_factor = 1.,
285  const double ext_factor = -1.)
286  {
287  const unsigned int n_dofs = fe1.dofs_per_cell;
288 
289  AssertDimension(fe1.get_fe().n_components(), dim);
290  AssertDimension(fe2.get_fe().n_components(), dim);
291  AssertDimension(M11.m(), n_dofs);
292  AssertDimension(M11.n(), n_dofs);
293  AssertDimension(M12.m(), n_dofs);
294  AssertDimension(M12.n(), n_dofs);
295  AssertDimension(M21.m(), n_dofs);
296  AssertDimension(M21.n(), n_dofs);
297  AssertDimension(M22.m(), n_dofs);
298  AssertDimension(M22.n(), n_dofs);
299 
300  const double nu1 = int_factor;
301  const double nu2 = (ext_factor < 0) ? int_factor : ext_factor;
302  const double penalty = .5 * pen * (nu1 + nu2);
303 
304  for (unsigned int k=0; k<fe1.n_quadrature_points; ++k)
305  {
306  const double dx = fe1.JxW(k);
307  const Point<dim> &n = fe1.normal_vector(k);
308  for (unsigned int i=0; i<n_dofs; ++i)
309  for (unsigned int j=0; j<n_dofs; ++j)
310  for (unsigned int d1=0; d1<dim; ++d1)
311  {
312  const double u1 = fe1.shape_value_component(j,k,d1);
313  const double u2 = fe2.shape_value_component(j,k,d1);
314  const double v1 = fe1.shape_value_component(i,k,d1);
315  const double v2 = fe2.shape_value_component(i,k,d1);
316 
317  M11(i,j) += dx * penalty * u1*v1;
318  M12(i,j) -= dx * penalty * u2*v1;
319  M21(i,j) -= dx * penalty * u1*v2;
320  M22(i,j) += dx * penalty * u2*v2;
321 
322  for (unsigned int d2=0; d2<dim; ++d2)
323  {
324  // v . nabla u n
325  M11(i,j) -= .25 * dx * nu1 * fe1.shape_grad_component(j,k,d1)[d2] * n(d2) * v1;
326  M12(i,j) -= .25 * dx * nu2 * fe2.shape_grad_component(j,k,d1)[d2] * n(d2) * v1;
327  M21(i,j) += .25 * dx * nu1 * fe1.shape_grad_component(j,k,d1)[d2] * n(d2) * v2;
328  M22(i,j) += .25 * dx * nu2 * fe2.shape_grad_component(j,k,d1)[d2] * n(d2) * v2;
329  // v (nabla u)^T n
330  M11(i,j) -= .25 * dx * nu1 * fe1.shape_grad_component(j,k,d2)[d1] * n(d2) * v1;
331  M12(i,j) -= .25 * dx * nu2 * fe2.shape_grad_component(j,k,d2)[d1] * n(d2) * v1;
332  M21(i,j) += .25 * dx * nu1 * fe1.shape_grad_component(j,k,d2)[d1] * n(d2) * v2;
333  M22(i,j) += .25 * dx * nu2 * fe2.shape_grad_component(j,k,d2)[d1] * n(d2) * v2;
334  // u nabla v n
335  M11(i,j) -= .25 * dx * nu1 * fe1.shape_grad_component(i,k,d1)[d2] * n(d2) * u1;
336  M12(i,j) += .25 * dx * nu1 * fe1.shape_grad_component(i,k,d1)[d2] * n(d2) * u2;
337  M21(i,j) -= .25 * dx * nu2 * fe2.shape_grad_component(i,k,d1)[d2] * n(d2) * u1;
338  M22(i,j) += .25 * dx * nu2 * fe2.shape_grad_component(i,k,d1)[d2] * n(d2) * u2;
339  // u (nabla v)^T n
340  M11(i,j) -= .25 * dx * nu1 * fe1.shape_grad_component(i,k,d2)[d1] * n(d2) * u1;
341  M12(i,j) += .25 * dx * nu1 * fe1.shape_grad_component(i,k,d2)[d1] * n(d2) * u2;
342  M21(i,j) -= .25 * dx * nu2 * fe2.shape_grad_component(i,k,d2)[d1] * n(d2) * u1;
343  M22(i,j) += .25 * dx * nu2 * fe2.shape_grad_component(i,k,d2)[d1] * n(d2) * u2;
344  }
345  }
346  }
347  }
354  template<int dim, typename number>
355  void
357  Vector<number> &result1,
358  Vector<number> &result2,
359  const FEValuesBase<dim> &fe1,
360  const FEValuesBase<dim> &fe2,
361  const VectorSlice<const std::vector<std::vector<double> > > &input1,
362  const VectorSlice<const std::vector<std::vector<Tensor<1,dim> > > > &Dinput1,
363  const VectorSlice<const std::vector<std::vector<double> > > &input2,
364  const VectorSlice<const std::vector<std::vector<Tensor<1,dim> > > > &Dinput2,
365  double pen,
366  double int_factor = 1.,
367  double ext_factor = -1.)
368  {
369  const unsigned int n1 = fe1.dofs_per_cell;
370 
371  AssertDimension(fe1.get_fe().n_components(), dim);
372  AssertDimension(fe2.get_fe().n_components(), dim);
377 
378  const double nu1 = int_factor;
379  const double nu2 = (ext_factor < 0) ? int_factor : ext_factor;
380  const double penalty = .5 * pen * (nu1 + nu2);
381 
382 
383  for (unsigned int k=0; k<fe1.n_quadrature_points; ++k)
384  {
385  const double dx = fe1.JxW(k);
386  const Point<dim> &n = fe1.normal_vector(k);
387 
388  for (unsigned int i=0; i<n1; ++i)
389  for (unsigned int d1=0; d1<dim; ++d1)
390  {
391  const double v1 = fe1.shape_value_component(i,k,d1);
392  const double v2 = fe2.shape_value_component(i,k,d1);
393  const double u1 = input1[d1][k];
394  const double u2 = input2[d1][k];
395 
396  result1(i) += dx * penalty * u1*v1;
397  result1(i) -= dx * penalty * u2*v1;
398  result2(i) -= dx * penalty * u1*v2;
399  result2(i) += dx * penalty * u2*v2;
400 
401  for (unsigned int d2=0; d2<dim; ++d2)
402  {
403  // v . nabla u n
404  result1(i) -= .25*dx* (nu1*Dinput1[d1][k][d2]+nu2*Dinput2[d1][k][d2]) * n(d2) * v1;
405  result2(i) += .25*dx* (nu1*Dinput1[d1][k][d2]+nu2*Dinput2[d1][k][d2]) * n(d2) * v2;
406  // v . (nabla u)^T n
407  result1(i) -= .25*dx* (nu1*Dinput1[d2][k][d1]+nu2*Dinput2[d2][k][d1]) * n(d2) * v1;
408  result2(i) += .25*dx* (nu1*Dinput1[d2][k][d1]+nu2*Dinput2[d2][k][d1]) * n(d2) * v2;
409  // u nabla v n
410  result1(i) -= .25*dx* nu1*fe1.shape_grad_component(i,k,d1)[d2] * n(d2) * (u1-u2);
411  result2(i) -= .25*dx* nu2*fe2.shape_grad_component(i,k,d1)[d2] * n(d2) * (u1-u2);
412  // u (nabla v)^T n
413  result1(i) -= .25*dx* nu1*fe1.shape_grad_component(i,k,d2)[d1] * n(d2) * (u1-u2);
414  result2(i) -= .25*dx* nu2*fe2.shape_grad_component(i,k,d2)[d1] * n(d2) * (u1-u2);
415  }
416  }
417  }
418  }
419  }
420 }
421 
422 DEAL_II_NAMESPACE_CLOSE
423 
424 #endif
double shape_value_component(const unsigned int function_no, const unsigned int point_no, const unsigned int component) const
void nitsche_residual_homogeneous(Vector< number > &result, const FEValuesBase< dim > &fe, const VectorSlice< const std::vector< std::vector< double > > > &input, const VectorSlice< const std::vector< std::vector< Tensor< 1, dim > > > > &Dinput, double penalty, double factor=1.)
Definition: elasticity.h:233
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:858
Tensor< 1, spacedim > shape_grad_component(const unsigned int function_no, const unsigned int point_no, const unsigned int component) const
const unsigned int dofs_per_cell
Definition: fe_values.h:1418
const Point< spacedim > & normal_vector(const unsigned int i) const
#define AssertVectorVectorDimension(vec, dim1, dim2)
Definition: exceptions.h:870
Library of integrals over cells and faces.
Definition: maxwell.h:31
void cell_residual(Vector< number > &result, const FEValuesBase< dim > &fe, const VectorSlice< const std::vector< std::vector< Tensor< 1, dim > > > > &input, double factor=1.)
Definition: elasticity.h:85
size_type n() const
#define Assert(cond, exc)
Definition: exceptions.h:299
std::size_t size() const
void ip_matrix(FullMatrix< double > &M11, FullMatrix< double > &M12, FullMatrix< double > &M21, FullMatrix< double > &M22, const FEValuesBase< dim > &fe1, const FEValuesBase< dim > &fe2, const double pen, const double int_factor=1., const double ext_factor=-1.)
Definition: elasticity.h:276
size_type m() const
const unsigned int n_quadrature_points
Definition: fe_values.h:1411
void nitsche_residual(Vector< number > &result, const FEValuesBase< dim > &fe, const VectorSlice< const std::vector< std::vector< double > > > &input, const VectorSlice< const std::vector< std::vector< Tensor< 1, dim > > > > &Dinput, const VectorSlice< const std::vector< std::vector< double > > > &data, double penalty, double factor=1.)
Definition: elasticity.h:176
void nitsche_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, double penalty, double factor=1.)
Definition: elasticity.h:121
const FiniteElement< dim, spacedim > & get_fe() const
::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
double JxW(const unsigned int quadrature_point) const
void ip_residual(Vector< number > &result1, Vector< number > &result2, const FEValuesBase< dim > &fe1, const FEValuesBase< dim > &fe2, const VectorSlice< const std::vector< std::vector< double > > > &input1, const VectorSlice< const std::vector< std::vector< Tensor< 1, dim > > > > &Dinput1, const VectorSlice< const std::vector< std::vector< double > > > &input2, const VectorSlice< const std::vector< std::vector< Tensor< 1, dim > > > > &Dinput2, double pen, double int_factor=1., double ext_factor=-1.)
Definition: elasticity.h:356
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const double factor=1.)
Definition: elasticity.h:51