Reference documentation for deal.II version 9.5.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
elasticity.h
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2010 - 2020 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
24
26#include <deal.II/fe/mapping.h>
27
29
31
33
34namespace 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
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 *
69 (fe.shape_grad_component(j, k, d1)[d2] +
70 fe.shape_grad_component(j, k, d2)[d1]) *
71 (fe.shape_grad_component(i, k, d1)[d2] +
72 fe.shape_grad_component(i, k, d2)[d1]);
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;
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]) *
106 (fe.shape_grad_component(i, k, d1)[d2] +
107 fe.shape_grad_component(i, k, d2)[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
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
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.;
198 double ngradun = 0.;
199 double ngradvn = 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);
205 ngradun += n * fe.shape_grad_component(j, k, d) * n[d];
206 ngradvn += n * fe.shape_grad_component(i, k, d) * n[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 *
221 fe.shape_grad_component(j, k, d1)[d2] *
222 n[d2] * v;
223 // v (nabla u)^T n
224 M(i, j) -= .5 * dx *
225 fe.shape_grad_component(j, k, d2)[d1] *
226 n[d2] * v;
227 // u nabla v n
228 M(i, j) -= .5 * dx *
229 fe.shape_grad_component(i, k, d1)[d2] *
230 n[d2] * u;
231 // u (nabla v)^T n
232 M(i, j) -= .5 * dx *
233 fe.shape_grad_component(i, k, d2)[d1] *
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.;
332 double ngradun = 0.;
333 double ngradvn = 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];
341 ngradvn += n * fe.shape_grad_component(i, k, d) * 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) *
360 fe.shape_grad_component(i, k, d1)[d2] *
361 n[d2];
362 // u (nabla v)^T n
363 result(i) -= .5 * dx * (u - g) *
364 fe.shape_grad_component(i, k, d2)[d1] *
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 *
481 fe1.shape_grad_component(j, k, d1)[d2] *
482 n[d2] * v1;
483 M12(i, j) -= .25 * dx * nu2 *
484 fe2.shape_grad_component(j, k, d1)[d2] *
485 n[d2] * v1;
486 M21(i, j) += .25 * dx * nu1 *
487 fe1.shape_grad_component(j, k, d1)[d2] *
488 n[d2] * v2;
489 M22(i, j) += .25 * dx * nu2 *
490 fe2.shape_grad_component(j, k, d1)[d2] *
491 n[d2] * v2;
492 // v (nabla u)^T n
493 M11(i, j) -= .25 * dx * nu1 *
494 fe1.shape_grad_component(j, k, d2)[d1] *
495 n[d2] * v1;
496 M12(i, j) -= .25 * dx * nu2 *
497 fe2.shape_grad_component(j, k, d2)[d1] *
498 n[d2] * v1;
499 M21(i, j) += .25 * dx * nu1 *
500 fe1.shape_grad_component(j, k, d2)[d1] *
501 n[d2] * v2;
502 M22(i, j) += .25 * dx * nu2 *
503 fe2.shape_grad_component(j, k, d2)[d1] *
504 n[d2] * v2;
505 // u nabla v n
506 M11(i, j) -= .25 * dx * nu1 *
507 fe1.shape_grad_component(i, k, d1)[d2] *
508 n[d2] * u1;
509 M12(i, j) += .25 * dx * nu1 *
510 fe1.shape_grad_component(i, k, d1)[d2] *
511 n[d2] * u2;
512 M21(i, j) -= .25 * dx * nu2 *
513 fe2.shape_grad_component(i, k, d1)[d2] *
514 n[d2] * u1;
515 M22(i, j) += .25 * dx * nu2 *
516 fe2.shape_grad_component(i, k, d1)[d2] *
517 n[d2] * u2;
518 // u (nabla v)^T n
519 M11(i, j) -= .25 * dx * nu1 *
520 fe1.shape_grad_component(i, k, d2)[d1] *
521 n[d2] * u1;
522 M12(i, j) += .25 * dx * nu1 *
523 fe1.shape_grad_component(i, k, d2)[d1] *
524 n[d2] * u2;
525 M21(i, j) -= .25 * dx * nu2 *
526 fe2.shape_grad_component(i, k, d2)[d1] *
527 n[d2] * u1;
528 M22(i, j) += .25 * dx * nu2 *
529 fe2.shape_grad_component(i, k, d2)[d1] *
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 *
606 fe1.shape_grad_component(i, k, d1)[d2] *
607 n[d2] * (u1 - u2);
608 result2(i) -= .25 * dx * nu2 *
609 fe2.shape_grad_component(i, k, d1)[d2] *
610 n[d2] * (u1 - u2);
611 // u (nabla v)^T n
612 result1(i) -= .25 * dx * nu1 *
613 fe1.shape_grad_component(i, k, d2)[d1] *
614 n[d2] * (u1 - u2);
615 result2(i) -= .25 * dx * nu2 *
616 fe2.shape_grad_component(i, k, d2)[d1] *
617 n[d2] * (u1 - u2);
618 }
619 }
620 }
621 }
622 } // namespace Elasticity
623} // namespace LocalIntegrators
624
626
627#endif
const unsigned int dofs_per_cell
Definition fe_values.h:2451
const Tensor< 1, spacedim > & normal_vector(const unsigned int q_point) const
double shape_value_component(const unsigned int i, const unsigned int q_point, const unsigned int component) const
const unsigned int n_quadrature_points
Definition fe_values.h:2433
Tensor< 1, spacedim > shape_grad_component(const unsigned int i, const unsigned int q_point, const unsigned int component) const
const FiniteElement< dim, spacedim > & get_fe() const
double JxW(const unsigned int q_point) const
unsigned int n_components() const
size_type n() const
size_type m() const
size_type size() const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
const unsigned int v1
#define Assert(cond, exc)
#define AssertVectorVectorDimension(VEC, DIM1, DIM2)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
void nitsche_tangential_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, double penalty, double factor=1.)
Definition elasticity.h:178
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
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
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_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
void nitsche_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, double penalty, double factor=1.)
Definition elasticity.h:123
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
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const double factor=1.)
Definition elasticity.h:51
Library of integrals over cells and faces.
Definition advection.h:35