Reference documentation for deal.II version Git 6d63218887 2020-10-30 17:17:53 -0400
elasticity.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2010 - 2019 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15
16 #ifndef dealii_integrators_elasticity_h
17 #define dealii_integrators_elasticity_h
18
19
20 #include <deal.II/base/config.h>
21
24
25 #include <deal.II/fe/fe_values.h>
26 #include <deal.II/fe/mapping.h>
27
29
31
33
34 namespace LocalIntegrators
35 {
41  namespace Elasticity
42  {
49  template <int dim>
50  inline void
52  const FEValuesBase<dim> &fe,
53  const double factor = 1.)
54  {
55  const unsigned int n_dofs = fe.dofs_per_cell;
56
57  AssertDimension(fe.get_fe().n_components(), dim);
58  AssertDimension(M.m(), n_dofs);
59  AssertDimension(M.n(), n_dofs);
60
61  for (unsigned int k = 0; k < fe.n_quadrature_points; ++k)
62  {
63  const double dx = factor * fe.JxW(k);
64  for (unsigned int i = 0; i < n_dofs; ++i)
65  for (unsigned int j = 0; j < n_dofs; ++j)
66  for (unsigned int d1 = 0; d1 < dim; ++d1)
67  for (unsigned int d2 = 0; d2 < dim; ++d2)
68  M(i, j) += dx * .25 *
73  }
74  }
75
76
82  template <int dim, typename number>
83  inline void
85  const FEValuesBase<dim> & fe,
86  const ArrayView<const std::vector<Tensor<1, dim>>> &input,
87  double factor = 1.)
88  {
89  const unsigned int nq = fe.n_quadrature_points;
90  const unsigned int n_dofs = fe.dofs_per_cell;
91  AssertDimension(fe.get_fe().n_components(), dim);
92
94  Assert(result.size() == n_dofs,
95  ExcDimensionMismatch(result.size(), n_dofs));
96
97  for (unsigned int k = 0; k < nq; ++k)
98  {
99  const double dx = factor * fe.JxW(k);
100  for (unsigned int i = 0; i < n_dofs; ++i)
101  for (unsigned int d1 = 0; d1 < dim; ++d1)
102  for (unsigned int d2 = 0; d2 < dim; ++d2)
103  {
104  result(i) += dx * .25 *
105  (input[d1][k][d2] + input[d2][k][d1]) *
108  }
109  }
110  }
111
112
121  template <int dim>
122  inline void
124  const FEValuesBase<dim> &fe,
125  double penalty,
126  double factor = 1.)
127  {
128  const unsigned int n_dofs = fe.dofs_per_cell;
129
130  AssertDimension(fe.get_fe().n_components(), dim);
131  AssertDimension(M.m(), n_dofs);
132  AssertDimension(M.n(), n_dofs);
133
134  for (unsigned int k = 0; k < fe.n_quadrature_points; ++k)
135  {
136  const double dx = factor * fe.JxW(k);
137  const Tensor<1, dim> n = fe.normal_vector(k);
138  for (unsigned int i = 0; i < n_dofs; ++i)
139  for (unsigned int j = 0; j < n_dofs; ++j)
140  for (unsigned int d1 = 0; d1 < dim; ++d1)
141  {
142  const double u = fe.shape_value_component(j, k, d1);
143  const double v = fe.shape_value_component(i, k, d1);
144  M(i, j) += dx * 2. * penalty * u * v;
145  for (unsigned int d2 = 0; d2 < dim; ++d2)
146  {
147  // v . nabla u n
148  M(i, j) -= .5 * dx *
149  fe.shape_grad_component(j, k, d1)[d2] * n[d2] *
150  v;
151  // v (nabla u)^T n
152  M(i, j) -= .5 * dx *
153  fe.shape_grad_component(j, k, d2)[d1] * n[d2] *
154  v;
155  // u nabla v n
156  M(i, j) -= .5 * dx *
157  fe.shape_grad_component(i, k, d1)[d2] * n[d2] *
158  u;
159  // u (nabla v)^T n
160  M(i, j) -= .5 * dx *
161  fe.shape_grad_component(i, k, d2)[d1] * n[d2] *
162  u;
163  }
164  }
165  }
166  }
167
176  template <int dim>
177  inline void
179  const FEValuesBase<dim> &fe,
180  double penalty,
181  double factor = 1.)
182  {
183  const unsigned int n_dofs = fe.dofs_per_cell;
184
185  AssertDimension(fe.get_fe().n_components(), dim);
186  AssertDimension(M.m(), n_dofs);
187  AssertDimension(M.n(), n_dofs);
188
189  for (unsigned int k = 0; k < fe.n_quadrature_points; ++k)
190  {
191  const double dx = factor * fe.JxW(k);
192  const Tensor<1, dim> n = fe.normal_vector(k);
193  for (unsigned int i = 0; i < n_dofs; ++i)
194  for (unsigned int j = 0; j < n_dofs; ++j)
195  {
196  double udotn = 0.;
197  double vdotn = 0.;
200
201  for (unsigned int d = 0; d < dim; ++d)
202  {
203  udotn += n[d] * fe.shape_value_component(j, k, d);
204  vdotn += n[d] * fe.shape_value_component(i, k, d);
207  }
208  for (unsigned int d1 = 0; d1 < dim; ++d1)
209  {
210  const double u =
211  fe.shape_value_component(j, k, d1) - udotn * n[d1];
212  const double v =
213  fe.shape_value_component(i, k, d1) - vdotn * n[d1];
214  M(i, j) += dx * 2. * penalty * u * v;
215  // Correct the gradients below and subtract normal component
216  M(i, j) += dx * (ngradun * v + ngradvn * u);
217  for (unsigned int d2 = 0; d2 < dim; ++d2)
218  {
219  // v . nabla u n
220  M(i, j) -= .5 * dx *
222  n[d2] * v;
223  // v (nabla u)^T n
224  M(i, j) -= .5 * dx *
226  n[d2] * v;
227  // u nabla v n
228  M(i, j) -= .5 * dx *
230  n[d2] * u;
231  // u (nabla v)^T n
232  M(i, j) -= .5 * dx *
234  n[d2] * u;
235  }
236  }
237  }
238  }
239  }
240
255  template <int dim, typename number>
256  void
258  const FEValuesBase<dim> & fe,
259  const ArrayView<const std::vector<double>> & input,
260  const ArrayView<const std::vector<Tensor<1, dim>>> &Dinput,
261  const ArrayView<const std::vector<double>> & data,
262  double penalty,
263  double factor = 1.)
264  {
265  const unsigned int n_dofs = fe.dofs_per_cell;
269
270  for (unsigned int k = 0; k < fe.n_quadrature_points; ++k)
271  {
272  const double dx = factor * fe.JxW(k);
273  const Tensor<1, dim> n = fe.normal_vector(k);
274  for (unsigned int i = 0; i < n_dofs; ++i)
275  for (unsigned int d1 = 0; d1 < dim; ++d1)
276  {
277  const double u = input[d1][k];
278  const double v = fe.shape_value_component(i, k, d1);
279  const double g = data[d1][k];
280  result(i) += dx * 2. * penalty * (u - g) * v;
281
282  for (unsigned int d2 = 0; d2 < dim; ++d2)
283  {
284  // v . nabla u n
285  result(i) -= .5 * dx * v * Dinput[d1][k][d2] * n[d2];
286  // v . (nabla u)^T n
287  result(i) -= .5 * dx * v * Dinput[d2][k][d1] * n[d2];
288  // u nabla v n
289  result(i) -= .5 * dx * (u - g) *
290  fe.shape_grad_component(i, k, d1)[d2] * n[d2];
291  // u (nabla v)^T n
292  result(i) -= .5 * dx * (u - g) *
293  fe.shape_grad_component(i, k, d2)[d1] * n[d2];
294  }
295  }
296  }
297  }
298
307  template <int dim, typename number>
308  inline void
310  Vector<number> & result,
311  const FEValuesBase<dim> & fe,
312  const ArrayView<const std::vector<double>> & input,
313  const ArrayView<const std::vector<Tensor<1, dim>>> &Dinput,
314  const ArrayView<const std::vector<double>> & data,
315  double penalty,
316  double factor = 1.)
317  {
318  const unsigned int n_dofs = fe.dofs_per_cell;
322
323  for (unsigned int k = 0; k < fe.n_quadrature_points; ++k)
324  {
325  const double dx = factor * fe.JxW(k);
326  const Tensor<1, dim> n = fe.normal_vector(k);
327  for (unsigned int i = 0; i < n_dofs; ++i)
328  {
329  double udotn = 0.;
330  double gdotn = 0.;
331  double vdotn = 0.;
334
335  for (unsigned int d = 0; d < dim; ++d)
336  {
337  udotn += n[d] * input[d][k];
338  gdotn += n[d] * data[d][k];
339  vdotn += n[d] * fe.shape_value_component(i, k, d);
340  ngradun += n * Dinput[d][k] * n[d];
342  }
343  for (unsigned int d1 = 0; d1 < dim; ++d1)
344  {
345  const double u = input[d1][k] - udotn * n[d1];
346  const double v =
347  fe.shape_value_component(i, k, d1) - vdotn * n[d1];
348  const double g = data[d1][k] - gdotn * n[d1];
349  result(i) += dx * 2. * penalty * (u - g) * v;
350  // Correct the gradients below and subtract normal component
351  result(i) += dx * (ngradun * v + ngradvn * (u - g));
352  for (unsigned int d2 = 0; d2 < dim; ++d2)
353  {
354  // v . nabla u n
355  result(i) -= .5 * dx * Dinput[d1][k][d2] * n[d2] * v;
356  // v (nabla u)^T n
357  result(i) -= .5 * dx * Dinput[d2][k][d1] * n[d2] * v;
358  // u nabla v n
359  result(i) -= .5 * dx * (u - g) *
361  n[d2];
362  // u (nabla v)^T n
363  result(i) -= .5 * dx * (u - g) *
365  n[d2];
366  }
367  }
368  }
369  }
370  }
371
385  template <int dim, typename number>
386  void
388  Vector<number> & result,
389  const FEValuesBase<dim> & fe,
390  const ArrayView<const std::vector<double>> & input,
391  const ArrayView<const std::vector<Tensor<1, dim>>> &Dinput,
392  double penalty,
393  double factor = 1.)
394  {
395  const unsigned int n_dofs = fe.dofs_per_cell;
398
399  for (unsigned int k = 0; k < fe.n_quadrature_points; ++k)
400  {
401  const double dx = factor * fe.JxW(k);
402  const Tensor<1, dim> n = fe.normal_vector(k);
403  for (unsigned int i = 0; i < n_dofs; ++i)
404  for (unsigned int d1 = 0; d1 < dim; ++d1)
405  {
406  const double u = input[d1][k];
407  const double v = fe.shape_value_component(i, k, d1);
408  result(i) += dx * 2. * penalty * u * v;
409
410  for (unsigned int d2 = 0; d2 < dim; ++d2)
411  {
412  // v . nabla u n
413  result(i) -= .5 * dx * v * Dinput[d1][k][d2] * n[d2];
414  // v . (nabla u)^T n
415  result(i) -= .5 * dx * v * Dinput[d2][k][d1] * n[d2];
416  // u nabla v n
417  result(i) -= .5 * dx * u *
418  fe.shape_grad_component(i, k, d1)[d2] * n[d2];
419  // u (nabla v)^T n
420  result(i) -= .5 * dx * u *
421  fe.shape_grad_component(i, k, d2)[d1] * n[d2];
422  }
423  }
424  }
425  }
426
430  template <int dim>
431  inline void
433  FullMatrix<double> & M12,
434  FullMatrix<double> & M21,
435  FullMatrix<double> & M22,
436  const FEValuesBase<dim> &fe1,
437  const FEValuesBase<dim> &fe2,
438  const double pen,
439  const double int_factor = 1.,
440  const double ext_factor = -1.)
441  {
442  const unsigned int n_dofs = fe1.dofs_per_cell;
443
444  AssertDimension(fe1.get_fe().n_components(), dim);
445  AssertDimension(fe2.get_fe().n_components(), dim);
446  AssertDimension(M11.m(), n_dofs);
447  AssertDimension(M11.n(), n_dofs);
448  AssertDimension(M12.m(), n_dofs);
449  AssertDimension(M12.n(), n_dofs);
450  AssertDimension(M21.m(), n_dofs);
451  AssertDimension(M21.n(), n_dofs);
452  AssertDimension(M22.m(), n_dofs);
453  AssertDimension(M22.n(), n_dofs);
454
455  const double nu1 = int_factor;
456  const double nu2 = (ext_factor < 0) ? int_factor : ext_factor;
457  const double penalty = .5 * pen * (nu1 + nu2);
458
459  for (unsigned int k = 0; k < fe1.n_quadrature_points; ++k)
460  {
461  const double dx = fe1.JxW(k);
462  const Tensor<1, dim> n = fe1.normal_vector(k);
463  for (unsigned int i = 0; i < n_dofs; ++i)
464  for (unsigned int j = 0; j < n_dofs; ++j)
465  for (unsigned int d1 = 0; d1 < dim; ++d1)
466  {
467  const double u1 = fe1.shape_value_component(j, k, d1);
468  const double u2 = fe2.shape_value_component(j, k, d1);
469  const double v1 = fe1.shape_value_component(i, k, d1);
470  const double v2 = fe2.shape_value_component(i, k, d1);
471
472  M11(i, j) += dx * penalty * u1 * v1;
473  M12(i, j) -= dx * penalty * u2 * v1;
474  M21(i, j) -= dx * penalty * u1 * v2;
475  M22(i, j) += dx * penalty * u2 * v2;
476
477  for (unsigned int d2 = 0; d2 < dim; ++d2)
478  {
479  // v . nabla u n
480  M11(i, j) -= .25 * dx * nu1 *
482  n[d2] * v1;
483  M12(i, j) -= .25 * dx * nu2 *
485  n[d2] * v1;
486  M21(i, j) += .25 * dx * nu1 *
488  n[d2] * v2;
489  M22(i, j) += .25 * dx * nu2 *
491  n[d2] * v2;
492  // v (nabla u)^T n
493  M11(i, j) -= .25 * dx * nu1 *
495  n[d2] * v1;
496  M12(i, j) -= .25 * dx * nu2 *
498  n[d2] * v1;
499  M21(i, j) += .25 * dx * nu1 *
501  n[d2] * v2;
502  M22(i, j) += .25 * dx * nu2 *
504  n[d2] * v2;
505  // u nabla v n
506  M11(i, j) -= .25 * dx * nu1 *
508  n[d2] * u1;
509  M12(i, j) += .25 * dx * nu1 *
511  n[d2] * u2;
512  M21(i, j) -= .25 * dx * nu2 *
514  n[d2] * u1;
515  M22(i, j) += .25 * dx * nu2 *
517  n[d2] * u2;
518  // u (nabla v)^T n
519  M11(i, j) -= .25 * dx * nu1 *
521  n[d2] * u1;
522  M12(i, j) += .25 * dx * nu1 *
524  n[d2] * u2;
525  M21(i, j) -= .25 * dx * nu2 *
527  n[d2] * u1;
528  M22(i, j) += .25 * dx * nu2 *
530  n[d2] * u2;
531  }
532  }
533  }
534  }
538  template <int dim, typename number>
539  void
541  Vector<number> & result2,
542  const FEValuesBase<dim> & fe1,
543  const FEValuesBase<dim> & fe2,
544  const ArrayView<const std::vector<double>> & input1,
545  const ArrayView<const std::vector<Tensor<1, dim>>> &Dinput1,
546  const ArrayView<const std::vector<double>> & input2,
547  const ArrayView<const std::vector<Tensor<1, dim>>> &Dinput2,
548  double pen,
549  double int_factor = 1.,
550  double ext_factor = -1.)
551  {
552  const unsigned int n1 = fe1.dofs_per_cell;
553
554  AssertDimension(fe1.get_fe().n_components(), dim);
555  AssertDimension(fe2.get_fe().n_components(), dim);
560
561  const double nu1 = int_factor;
562  const double nu2 = (ext_factor < 0) ? int_factor : ext_factor;
563  const double penalty = .5 * pen * (nu1 + nu2);
564
565
566  for (unsigned int k = 0; k < fe1.n_quadrature_points; ++k)
567  {
568  const double dx = fe1.JxW(k);
569  const Tensor<1, dim> n = fe1.normal_vector(k);
570
571  for (unsigned int i = 0; i < n1; ++i)
572  for (unsigned int d1 = 0; d1 < dim; ++d1)
573  {
574  const double v1 = fe1.shape_value_component(i, k, d1);
575  const double v2 = fe2.shape_value_component(i, k, d1);
576  const double u1 = input1[d1][k];
577  const double u2 = input2[d1][k];
578
579  result1(i) += dx * penalty * u1 * v1;
580  result1(i) -= dx * penalty * u2 * v1;
581  result2(i) -= dx * penalty * u1 * v2;
582  result2(i) += dx * penalty * u2 * v2;
583
584  for (unsigned int d2 = 0; d2 < dim; ++d2)
585  {
586  // v . nabla u n
587  result1(i) -=
588  .25 * dx *
589  (nu1 * Dinput1[d1][k][d2] + nu2 * Dinput2[d1][k][d2]) *
590  n[d2] * v1;
591  result2(i) +=
592  .25 * dx *
593  (nu1 * Dinput1[d1][k][d2] + nu2 * Dinput2[d1][k][d2]) *
594  n[d2] * v2;
595  // v . (nabla u)^T n
596  result1(i) -=
597  .25 * dx *
598  (nu1 * Dinput1[d2][k][d1] + nu2 * Dinput2[d2][k][d1]) *
599  n[d2] * v1;
600  result2(i) +=
601  .25 * dx *
602  (nu1 * Dinput1[d2][k][d1] + nu2 * Dinput2[d2][k][d1]) *
603  n[d2] * v2;
604  // u nabla v n
605  result1(i) -= .25 * dx * nu1 *
607  n[d2] * (u1 - u2);
608  result2(i) -= .25 * dx * nu2 *
610  n[d2] * (u1 - u2);
611  // u (nabla v)^T n
612  result1(i) -= .25 * dx * nu1 *
614  n[d2] * (u1 - u2);
615  result2(i) -= .25 * dx * nu2 *
617  n[d2] * (u1 - u2);
618  }
619  }
620  }
621  }
622  } // namespace Elasticity
623 } // namespace LocalIntegrators
624
626
627 #endif
void cell_residual(Vector< number > &result, const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim >>> &input, double factor=1.)
Definition: elasticity.h:84
size_type m() const
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1580
const unsigned int dofs_per_cell
Definition: fe_values.h:2097
const FiniteElement< dim, spacedim > & get_fe() const
#define AssertVectorVectorDimension(VEC, DIM1, DIM2)
Definition: exceptions.h:1607
void nitsche_tangential_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, double penalty, double factor=1.)
Definition: elasticity.h:178
void nitsche_tangential_residual(Vector< number > &result, const FEValuesBase< dim > &fe, const ArrayView< const std::vector< double >> &input, const ArrayView< const std::vector< Tensor< 1, dim >>> &Dinput, const ArrayView< const std::vector< double >> &data, double penalty, double factor=1.)
Definition: elasticity.h:309
double shape_value_component(const unsigned int function_no, const unsigned int point_no, const unsigned int component) const
Library of integrals over cells and faces.
size_type n() const
#define Assert(cond, exc)
Definition: exceptions.h:1423
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:369
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:432
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Definition: fe_values.h:2090
void ip_residual(Vector< number > &result1, Vector< number > &result2, const FEValuesBase< dim > &fe1, const FEValuesBase< dim > &fe2, const ArrayView< const std::vector< double >> &input1, const ArrayView< const std::vector< Tensor< 1, dim >>> &Dinput1, const ArrayView< const std::vector< double >> &input2, const ArrayView< const std::vector< Tensor< 1, dim >>> &Dinput2, double pen, double int_factor=1., double ext_factor=-1.)
Definition: elasticity.h:540
void nitsche_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, double penalty, double factor=1.)
Definition: elasticity.h:123
Tensor< 1, spacedim > shape_grad_component(const unsigned int function_no, const unsigned int point_no, const unsigned int component) const
double JxW(const unsigned int quadrature_point) const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:368
size_type size() const
void nitsche_residual(Vector< number > &result, const FEValuesBase< dim > &fe, const ArrayView< const std::vector< double >> &input, const ArrayView< const std::vector< Tensor< 1, dim >>> &Dinput, const ArrayView< const std::vector< double >> &data, double penalty, double factor=1.)
Definition: elasticity.h:257
void nitsche_residual_homogeneous(Vector< number > &result, const FEValuesBase< dim > &fe, const ArrayView< const std::vector< double >> &input, const ArrayView< const std::vector< Tensor< 1, dim >>> &Dinput, double penalty, double factor=1.)
Definition: elasticity.h:387
const Tensor< 1, spacedim > & normal_vector(const unsigned int i) const
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const double factor=1.)
Definition: elasticity.h:51