Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
elasticity.h
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
9 // Public License as published by the Free Software Foundation; either
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 
22 #include <deal.II/base/exceptions.h>
23 #include <deal.II/base/quadrature.h>
24 
25 #include <deal.II/fe/fe_values.h>
26 #include <deal.II/fe/mapping.h>
27 
28 #include <deal.II/lac/full_matrix.h>
29 
30 #include <deal.II/meshworker/dof_info.h>
31 
32 DEAL_II_NAMESPACE_OPEN
33 
34 namespace LocalIntegrators
35 {
43  namespace Elasticity
44  {
51  template <int dim>
52  inline void
54  const FEValuesBase<dim> &fe,
55  const double factor = 1.)
56  {
57  const unsigned int n_dofs = fe.dofs_per_cell;
58 
59  AssertDimension(fe.get_fe().n_components(), dim);
60  AssertDimension(M.m(), n_dofs);
61  AssertDimension(M.n(), n_dofs);
62 
63  for (unsigned int k = 0; k < fe.n_quadrature_points; ++k)
64  {
65  const double dx = factor * fe.JxW(k);
66  for (unsigned int i = 0; i < n_dofs; ++i)
67  for (unsigned int j = 0; j < n_dofs; ++j)
68  for (unsigned int d1 = 0; d1 < dim; ++d1)
69  for (unsigned int d2 = 0; d2 < dim; ++d2)
70  M(i, j) += dx * .25 *
71  (fe.shape_grad_component(j, k, d1)[d2] +
72  fe.shape_grad_component(j, k, d2)[d1]) *
73  (fe.shape_grad_component(i, k, d1)[d2] +
74  fe.shape_grad_component(i, k, d2)[d1]);
75  }
76  }
77 
78 
84  template <int dim, typename number>
85  inline void
87  const FEValuesBase<dim> & fe,
88  const ArrayView<const 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,
97  ExcDimensionMismatch(result.size(), n_dofs));
98 
99  for (unsigned int k = 0; k < nq; ++k)
100  {
101  const double dx = factor * fe.JxW(k);
102  for (unsigned int i = 0; i < n_dofs; ++i)
103  for (unsigned int d1 = 0; d1 < dim; ++d1)
104  for (unsigned int d2 = 0; d2 < dim; ++d2)
105  {
106  result(i) += dx * .25 *
107  (input[d1][k][d2] + input[d2][k][d1]) *
108  (fe.shape_grad_component(i, k, d1)[d2] +
109  fe.shape_grad_component(i, k, d2)[d1]);
110  }
111  }
112  }
113 
114 
123  template <int dim>
124  inline void
126  const FEValuesBase<dim> &fe,
127  double penalty,
128  double factor = 1.)
129  {
130  const unsigned int n_dofs = fe.dofs_per_cell;
131 
132  AssertDimension(fe.get_fe().n_components(), dim);
133  AssertDimension(M.m(), n_dofs);
134  AssertDimension(M.n(), n_dofs);
135 
136  for (unsigned int k = 0; k < fe.n_quadrature_points; ++k)
137  {
138  const double dx = factor * fe.JxW(k);
139  const Tensor<1, dim> n = fe.normal_vector(k);
140  for (unsigned int i = 0; i < n_dofs; ++i)
141  for (unsigned int j = 0; j < n_dofs; ++j)
142  for (unsigned int d1 = 0; d1 < dim; ++d1)
143  {
144  const double u = fe.shape_value_component(j, k, d1);
145  const double v = fe.shape_value_component(i, k, d1);
146  M(i, j) += dx * 2. * penalty * u * v;
147  for (unsigned int d2 = 0; d2 < dim; ++d2)
148  {
149  // v . nabla u n
150  M(i, j) -= .5 * dx *
151  fe.shape_grad_component(j, k, d1)[d2] * n[d2] *
152  v;
153  // v (nabla u)^T n
154  M(i, j) -= .5 * dx *
155  fe.shape_grad_component(j, k, d2)[d1] * n[d2] *
156  v;
157  // u nabla v n
158  M(i, j) -= .5 * dx *
159  fe.shape_grad_component(i, k, d1)[d2] * n[d2] *
160  u;
161  // u (nabla v)^T n
162  M(i, j) -= .5 * dx *
163  fe.shape_grad_component(i, k, d2)[d1] * n[d2] *
164  u;
165  }
166  }
167  }
168  }
169 
178  template <int dim>
179  inline void
181  const FEValuesBase<dim> &fe,
182  double penalty,
183  double factor = 1.)
184  {
185  const unsigned int n_dofs = fe.dofs_per_cell;
186 
187  AssertDimension(fe.get_fe().n_components(), dim);
188  AssertDimension(M.m(), n_dofs);
189  AssertDimension(M.n(), n_dofs);
190 
191  for (unsigned int k = 0; k < fe.n_quadrature_points; ++k)
192  {
193  const double dx = factor * fe.JxW(k);
194  const Tensor<1, dim> n = fe.normal_vector(k);
195  for (unsigned int i = 0; i < n_dofs; ++i)
196  for (unsigned int j = 0; j < n_dofs; ++j)
197  {
198  double udotn = 0.;
199  double vdotn = 0.;
200  double ngradun = 0.;
201  double ngradvn = 0.;
202 
203  for (unsigned int d = 0; d < dim; ++d)
204  {
205  udotn += n[d] * fe.shape_value_component(j, k, d);
206  vdotn += n[d] * fe.shape_value_component(i, k, d);
207  ngradun += n * fe.shape_grad_component(j, k, d) * n[d];
208  ngradvn += n * fe.shape_grad_component(i, k, d) * n[d];
209  }
210  for (unsigned int d1 = 0; d1 < dim; ++d1)
211  {
212  const double u =
213  fe.shape_value_component(j, k, d1) - udotn * n[d1];
214  const double v =
215  fe.shape_value_component(i, k, d1) - vdotn * n[d1];
216  M(i, j) += dx * 2. * penalty * u * v;
217  // Correct the gradients below and subtract normal component
218  M(i, j) += dx * (ngradun * v + ngradvn * u);
219  for (unsigned int d2 = 0; d2 < dim; ++d2)
220  {
221  // v . nabla u n
222  M(i, j) -= .5 * dx *
223  fe.shape_grad_component(j, k, d1)[d2] *
224  n[d2] * v;
225  // v (nabla u)^T n
226  M(i, j) -= .5 * dx *
227  fe.shape_grad_component(j, k, d2)[d1] *
228  n[d2] * v;
229  // u nabla v n
230  M(i, j) -= .5 * dx *
231  fe.shape_grad_component(i, k, d1)[d2] *
232  n[d2] * u;
233  // u (nabla v)^T n
234  M(i, j) -= .5 * dx *
235  fe.shape_grad_component(i, k, d2)[d1] *
236  n[d2] * u;
237  }
238  }
239  }
240  }
241  }
242 
260  template <int dim, typename number>
261  void
263  const FEValuesBase<dim> & fe,
264  const ArrayView<const std::vector<double>> & input,
265  const ArrayView<const std::vector<Tensor<1, dim>>> &Dinput,
266  const ArrayView<const std::vector<double>> & data,
267  double penalty,
268  double factor = 1.)
269  {
270  const unsigned int n_dofs = fe.dofs_per_cell;
274 
275  for (unsigned int k = 0; k < fe.n_quadrature_points; ++k)
276  {
277  const double dx = factor * fe.JxW(k);
278  const Tensor<1, dim> n = fe.normal_vector(k);
279  for (unsigned int i = 0; i < n_dofs; ++i)
280  for (unsigned int d1 = 0; d1 < dim; ++d1)
281  {
282  const double u = input[d1][k];
283  const double v = fe.shape_value_component(i, k, d1);
284  const double g = data[d1][k];
285  result(i) += dx * 2. * penalty * (u - g) * v;
286 
287  for (unsigned int d2 = 0; d2 < dim; ++d2)
288  {
289  // v . nabla u n
290  result(i) -= .5 * dx * v * Dinput[d1][k][d2] * n[d2];
291  // v . (nabla u)^T n
292  result(i) -= .5 * dx * v * Dinput[d2][k][d1] * n[d2];
293  // u nabla v n
294  result(i) -= .5 * dx * (u - g) *
295  fe.shape_grad_component(i, k, d1)[d2] * n[d2];
296  // u (nabla v)^T n
297  result(i) -= .5 * dx * (u - g) *
298  fe.shape_grad_component(i, k, d2)[d1] * n[d2];
299  }
300  }
301  }
302  }
303 
312  template <int dim, typename number>
313  inline void
315  Vector<number> & result,
316  const FEValuesBase<dim> & fe,
317  const ArrayView<const std::vector<double>> & input,
318  const ArrayView<const std::vector<Tensor<1, dim>>> &Dinput,
319  const ArrayView<const std::vector<double>> & data,
320  double penalty,
321  double factor = 1.)
322  {
323  const unsigned int n_dofs = fe.dofs_per_cell;
327 
328  for (unsigned int k = 0; k < fe.n_quadrature_points; ++k)
329  {
330  const double dx = factor * fe.JxW(k);
331  const Tensor<1, dim> n = fe.normal_vector(k);
332  for (unsigned int i = 0; i < n_dofs; ++i)
333  {
334  double udotn = 0.;
335  double gdotn = 0.;
336  double vdotn = 0.;
337  double ngradun = 0.;
338  double ngradvn = 0.;
339 
340  for (unsigned int d = 0; d < dim; ++d)
341  {
342  udotn += n[d] * input[d][k];
343  gdotn += n[d] * data[d][k];
344  vdotn += n[d] * fe.shape_value_component(i, k, d);
345  ngradun += n * Dinput[d][k] * n[d];
346  ngradvn += n * fe.shape_grad_component(i, k, d) * n[d];
347  }
348  for (unsigned int d1 = 0; d1 < dim; ++d1)
349  {
350  const double u = input[d1][k] - udotn * n[d1];
351  const double v =
352  fe.shape_value_component(i, k, d1) - vdotn * n[d1];
353  const double g = data[d1][k] - gdotn * n[d1];
354  result(i) += dx * 2. * penalty * (u - g) * v;
355  // Correct the gradients below and subtract normal component
356  result(i) += dx * (ngradun * v + ngradvn * (u - g));
357  for (unsigned int d2 = 0; d2 < dim; ++d2)
358  {
359  // v . nabla u n
360  result(i) -= .5 * dx * Dinput[d1][k][d2] * n[d2] * v;
361  // v (nabla u)^T n
362  result(i) -= .5 * dx * Dinput[d2][k][d1] * n[d2] * v;
363  // u nabla v n
364  result(i) -= .5 * dx * (u - g) *
365  fe.shape_grad_component(i, k, d1)[d2] *
366  n[d2];
367  // u (nabla v)^T n
368  result(i) -= .5 * dx * (u - g) *
369  fe.shape_grad_component(i, k, d2)[d1] *
370  n[d2];
371  }
372  }
373  }
374  }
375  }
376 
393  template <int dim, typename number>
394  void
396  Vector<number> & result,
397  const FEValuesBase<dim> & fe,
398  const ArrayView<const std::vector<double>> & input,
399  const ArrayView<const std::vector<Tensor<1, dim>>> &Dinput,
400  double penalty,
401  double factor = 1.)
402  {
403  const unsigned int n_dofs = fe.dofs_per_cell;
406 
407  for (unsigned int k = 0; k < fe.n_quadrature_points; ++k)
408  {
409  const double dx = factor * fe.JxW(k);
410  const Tensor<1, dim> n = fe.normal_vector(k);
411  for (unsigned int i = 0; i < n_dofs; ++i)
412  for (unsigned int d1 = 0; d1 < dim; ++d1)
413  {
414  const double u = input[d1][k];
415  const double v = fe.shape_value_component(i, k, d1);
416  result(i) += dx * 2. * penalty * u * v;
417 
418  for (unsigned int d2 = 0; d2 < dim; ++d2)
419  {
420  // v . nabla u n
421  result(i) -= .5 * dx * v * Dinput[d1][k][d2] * n[d2];
422  // v . (nabla u)^T n
423  result(i) -= .5 * dx * v * Dinput[d2][k][d1] * n[d2];
424  // u nabla v n
425  result(i) -= .5 * dx * u *
426  fe.shape_grad_component(i, k, d1)[d2] * n[d2];
427  // u (nabla v)^T n
428  result(i) -= .5 * dx * u *
429  fe.shape_grad_component(i, k, d2)[d1] * n[d2];
430  }
431  }
432  }
433  }
434 
438  template <int dim>
439  inline void
441  FullMatrix<double> & M12,
442  FullMatrix<double> & M21,
443  FullMatrix<double> & M22,
444  const FEValuesBase<dim> &fe1,
445  const FEValuesBase<dim> &fe2,
446  const double pen,
447  const double int_factor = 1.,
448  const double ext_factor = -1.)
449  {
450  const unsigned int n_dofs = fe1.dofs_per_cell;
451 
452  AssertDimension(fe1.get_fe().n_components(), dim);
453  AssertDimension(fe2.get_fe().n_components(), dim);
454  AssertDimension(M11.m(), n_dofs);
455  AssertDimension(M11.n(), n_dofs);
456  AssertDimension(M12.m(), n_dofs);
457  AssertDimension(M12.n(), n_dofs);
458  AssertDimension(M21.m(), n_dofs);
459  AssertDimension(M21.n(), n_dofs);
460  AssertDimension(M22.m(), n_dofs);
461  AssertDimension(M22.n(), n_dofs);
462 
463  const double nu1 = int_factor;
464  const double nu2 = (ext_factor < 0) ? int_factor : ext_factor;
465  const double penalty = .5 * pen * (nu1 + nu2);
466 
467  for (unsigned int k = 0; k < fe1.n_quadrature_points; ++k)
468  {
469  const double dx = fe1.JxW(k);
470  const Tensor<1, dim> n = fe1.normal_vector(k);
471  for (unsigned int i = 0; i < n_dofs; ++i)
472  for (unsigned int j = 0; j < n_dofs; ++j)
473  for (unsigned int d1 = 0; d1 < dim; ++d1)
474  {
475  const double u1 = fe1.shape_value_component(j, k, d1);
476  const double u2 = fe2.shape_value_component(j, k, d1);
477  const double v1 = fe1.shape_value_component(i, k, d1);
478  const double v2 = fe2.shape_value_component(i, k, d1);
479 
480  M11(i, j) += dx * penalty * u1 * v1;
481  M12(i, j) -= dx * penalty * u2 * v1;
482  M21(i, j) -= dx * penalty * u1 * v2;
483  M22(i, j) += dx * penalty * u2 * v2;
484 
485  for (unsigned int d2 = 0; d2 < dim; ++d2)
486  {
487  // v . nabla u n
488  M11(i, j) -= .25 * dx * nu1 *
489  fe1.shape_grad_component(j, k, d1)[d2] *
490  n[d2] * v1;
491  M12(i, j) -= .25 * dx * nu2 *
492  fe2.shape_grad_component(j, k, d1)[d2] *
493  n[d2] * v1;
494  M21(i, j) += .25 * dx * nu1 *
495  fe1.shape_grad_component(j, k, d1)[d2] *
496  n[d2] * v2;
497  M22(i, j) += .25 * dx * nu2 *
498  fe2.shape_grad_component(j, k, d1)[d2] *
499  n[d2] * v2;
500  // v (nabla u)^T n
501  M11(i, j) -= .25 * dx * nu1 *
502  fe1.shape_grad_component(j, k, d2)[d1] *
503  n[d2] * v1;
504  M12(i, j) -= .25 * dx * nu2 *
505  fe2.shape_grad_component(j, k, d2)[d1] *
506  n[d2] * v1;
507  M21(i, j) += .25 * dx * nu1 *
508  fe1.shape_grad_component(j, k, d2)[d1] *
509  n[d2] * v2;
510  M22(i, j) += .25 * dx * nu2 *
511  fe2.shape_grad_component(j, k, d2)[d1] *
512  n[d2] * v2;
513  // u nabla v n
514  M11(i, j) -= .25 * dx * nu1 *
515  fe1.shape_grad_component(i, k, d1)[d2] *
516  n[d2] * u1;
517  M12(i, j) += .25 * dx * nu1 *
518  fe1.shape_grad_component(i, k, d1)[d2] *
519  n[d2] * u2;
520  M21(i, j) -= .25 * dx * nu2 *
521  fe2.shape_grad_component(i, k, d1)[d2] *
522  n[d2] * u1;
523  M22(i, j) += .25 * dx * nu2 *
524  fe2.shape_grad_component(i, k, d1)[d2] *
525  n[d2] * u2;
526  // u (nabla v)^T n
527  M11(i, j) -= .25 * dx * nu1 *
528  fe1.shape_grad_component(i, k, d2)[d1] *
529  n[d2] * u1;
530  M12(i, j) += .25 * dx * nu1 *
531  fe1.shape_grad_component(i, k, d2)[d1] *
532  n[d2] * u2;
533  M21(i, j) -= .25 * dx * nu2 *
534  fe2.shape_grad_component(i, k, d2)[d1] *
535  n[d2] * u1;
536  M22(i, j) += .25 * dx * nu2 *
537  fe2.shape_grad_component(i, k, d2)[d1] *
538  n[d2] * u2;
539  }
540  }
541  }
542  }
549  template <int dim, typename number>
550  void
552  Vector<number> & result2,
553  const FEValuesBase<dim> & fe1,
554  const FEValuesBase<dim> & fe2,
555  const ArrayView<const std::vector<double>> & input1,
556  const ArrayView<const std::vector<Tensor<1, dim>>> &Dinput1,
557  const ArrayView<const std::vector<double>> & input2,
558  const ArrayView<const std::vector<Tensor<1, dim>>> &Dinput2,
559  double pen,
560  double int_factor = 1.,
561  double ext_factor = -1.)
562  {
563  const unsigned int n1 = fe1.dofs_per_cell;
564 
565  AssertDimension(fe1.get_fe().n_components(), dim);
566  AssertDimension(fe2.get_fe().n_components(), dim);
571 
572  const double nu1 = int_factor;
573  const double nu2 = (ext_factor < 0) ? int_factor : ext_factor;
574  const double penalty = .5 * pen * (nu1 + nu2);
575 
576 
577  for (unsigned int k = 0; k < fe1.n_quadrature_points; ++k)
578  {
579  const double dx = fe1.JxW(k);
580  const Tensor<1, dim> n = fe1.normal_vector(k);
581 
582  for (unsigned int i = 0; i < n1; ++i)
583  for (unsigned int d1 = 0; d1 < dim; ++d1)
584  {
585  const double v1 = fe1.shape_value_component(i, k, d1);
586  const double v2 = fe2.shape_value_component(i, k, d1);
587  const double u1 = input1[d1][k];
588  const double u2 = input2[d1][k];
589 
590  result1(i) += dx * penalty * u1 * v1;
591  result1(i) -= dx * penalty * u2 * v1;
592  result2(i) -= dx * penalty * u1 * v2;
593  result2(i) += dx * penalty * u2 * v2;
594 
595  for (unsigned int d2 = 0; d2 < dim; ++d2)
596  {
597  // v . nabla u n
598  result1(i) -=
599  .25 * dx *
600  (nu1 * Dinput1[d1][k][d2] + nu2 * Dinput2[d1][k][d2]) *
601  n[d2] * v1;
602  result2(i) +=
603  .25 * dx *
604  (nu1 * Dinput1[d1][k][d2] + nu2 * Dinput2[d1][k][d2]) *
605  n[d2] * v2;
606  // v . (nabla u)^T n
607  result1(i) -=
608  .25 * dx *
609  (nu1 * Dinput1[d2][k][d1] + nu2 * Dinput2[d2][k][d1]) *
610  n[d2] * v1;
611  result2(i) +=
612  .25 * dx *
613  (nu1 * Dinput1[d2][k][d1] + nu2 * Dinput2[d2][k][d1]) *
614  n[d2] * v2;
615  // u nabla v n
616  result1(i) -= .25 * dx * nu1 *
617  fe1.shape_grad_component(i, k, d1)[d2] *
618  n[d2] * (u1 - u2);
619  result2(i) -= .25 * dx * nu2 *
620  fe2.shape_grad_component(i, k, d1)[d2] *
621  n[d2] * (u1 - u2);
622  // u (nabla v)^T n
623  result1(i) -= .25 * dx * nu1 *
624  fe1.shape_grad_component(i, k, d2)[d1] *
625  n[d2] * (u1 - u2);
626  result2(i) -= .25 * dx * nu2 *
627  fe2.shape_grad_component(i, k, d2)[d1] *
628  n[d2] * (u1 - u2);
629  }
630  }
631  }
632  }
633  } // namespace Elasticity
634 } // namespace LocalIntegrators
635 
636 DEAL_II_NAMESPACE_CLOSE
637 
638 #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:86
size_type m() const
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1567
const unsigned int dofs_per_cell
Definition: fe_values.h:2106
const FiniteElement< dim, spacedim > & get_fe() const
#define AssertVectorVectorDimension(VEC, DIM1, DIM2)
Definition: exceptions.h:1595
void nitsche_tangential_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, double penalty, double factor=1.)
Definition: elasticity.h:180
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:314
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.
Definition: advection.h:34
size_type n() const
#define Assert(cond, exc)
Definition: exceptions.h:1407
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
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:440
const unsigned int n_quadrature_points
Definition: fe_values.h:2099
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:551
void nitsche_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, double penalty, double factor=1.)
Definition: elasticity.h:125
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
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:262
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:395
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:53