Reference documentation for deal.II version GIT relicensing-439-g5fda5c893d 2024-04-20 06:50:02+00:00
\(\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// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2010 - 2023 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_integrators_elasticity_h
16#define dealii_integrators_elasticity_h
17
18
19#include <deal.II/base/config.h>
20
23
25#include <deal.II/fe/mapping.h>
26
28
30
32
33namespace LocalIntegrators
34{
40 namespace Elasticity
41 {
48 template <int dim>
49 inline void
51 const FEValuesBase<dim> &fe,
52 const double factor = 1.)
53 {
54 const unsigned int n_dofs = fe.dofs_per_cell;
55
57 AssertDimension(M.m(), n_dofs);
58 AssertDimension(M.n(), n_dofs);
59
60 for (unsigned int k = 0; k < fe.n_quadrature_points; ++k)
61 {
62 const double dx = factor * fe.JxW(k);
63 for (unsigned int i = 0; i < n_dofs; ++i)
64 for (unsigned int j = 0; j < n_dofs; ++j)
65 for (unsigned int d1 = 0; d1 < dim; ++d1)
66 for (unsigned int d2 = 0; d2 < dim; ++d2)
67 M(i, j) += dx * .25 *
68 (fe.shape_grad_component(j, k, d1)[d2] +
69 fe.shape_grad_component(j, k, d2)[d1]) *
70 (fe.shape_grad_component(i, k, d1)[d2] +
71 fe.shape_grad_component(i, k, d2)[d1]);
72 }
73 }
74
75
81 template <int dim, typename number>
82 inline void
84 const FEValuesBase<dim> &fe,
85 const ArrayView<const std::vector<Tensor<1, dim>>> &input,
86 double factor = 1.)
87 {
88 const unsigned int nq = fe.n_quadrature_points;
89 const unsigned int n_dofs = fe.dofs_per_cell;
91
93 Assert(result.size() == n_dofs,
94 ExcDimensionMismatch(result.size(), n_dofs));
95
96 for (unsigned int k = 0; k < nq; ++k)
97 {
98 const double dx = factor * fe.JxW(k);
99 for (unsigned int i = 0; i < n_dofs; ++i)
100 for (unsigned int d1 = 0; d1 < dim; ++d1)
101 for (unsigned int d2 = 0; d2 < dim; ++d2)
102 {
103 result(i) += dx * .25 *
104 (input[d1][k][d2] + input[d2][k][d1]) *
105 (fe.shape_grad_component(i, k, d1)[d2] +
106 fe.shape_grad_component(i, k, d2)[d1]);
107 }
108 }
109 }
110
111
120 template <int dim>
121 inline void
123 const FEValuesBase<dim> &fe,
124 double penalty,
125 double factor = 1.)
126 {
127 const unsigned int n_dofs = fe.dofs_per_cell;
128
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 Tensor<1, 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 *
148 fe.shape_grad_component(j, k, d1)[d2] * n[d2] *
149 v;
150 // v (nabla u)^T n
151 M(i, j) -= .5 * dx *
152 fe.shape_grad_component(j, k, d2)[d1] * n[d2] *
153 v;
154 // u nabla v n
155 M(i, j) -= .5 * dx *
156 fe.shape_grad_component(i, k, d1)[d2] * n[d2] *
157 u;
158 // u (nabla v)^T n
159 M(i, j) -= .5 * dx *
160 fe.shape_grad_component(i, k, d2)[d1] * n[d2] *
161 u;
162 }
163 }
164 }
165 }
166
175 template <int dim>
176 inline void
178 const FEValuesBase<dim> &fe,
179 double penalty,
180 double factor = 1.)
181 {
182 const unsigned int n_dofs = fe.dofs_per_cell;
183
185 AssertDimension(M.m(), n_dofs);
186 AssertDimension(M.n(), n_dofs);
187
188 for (unsigned int k = 0; k < fe.n_quadrature_points; ++k)
189 {
190 const double dx = factor * fe.JxW(k);
191 const Tensor<1, dim> n = fe.normal_vector(k);
192 for (unsigned int i = 0; i < n_dofs; ++i)
193 for (unsigned int j = 0; j < n_dofs; ++j)
194 {
195 double udotn = 0.;
196 double vdotn = 0.;
197 double ngradun = 0.;
198 double ngradvn = 0.;
199
200 for (unsigned int d = 0; d < dim; ++d)
201 {
202 udotn += n[d] * fe.shape_value_component(j, k, d);
203 vdotn += n[d] * fe.shape_value_component(i, k, d);
204 ngradun += n * fe.shape_grad_component(j, k, d) * n[d];
205 ngradvn += n * fe.shape_grad_component(i, k, d) * n[d];
206 }
207 for (unsigned int d1 = 0; d1 < dim; ++d1)
208 {
209 const double u =
210 fe.shape_value_component(j, k, d1) - udotn * n[d1];
211 const double v =
212 fe.shape_value_component(i, k, d1) - vdotn * n[d1];
213 M(i, j) += dx * 2. * penalty * u * v;
214 // Correct the gradients below and subtract normal component
215 M(i, j) += dx * (ngradun * v + ngradvn * u);
216 for (unsigned int d2 = 0; d2 < dim; ++d2)
217 {
218 // v . nabla u n
219 M(i, j) -= .5 * dx *
220 fe.shape_grad_component(j, k, d1)[d2] *
221 n[d2] * v;
222 // v (nabla u)^T n
223 M(i, j) -= .5 * dx *
224 fe.shape_grad_component(j, k, d2)[d1] *
225 n[d2] * v;
226 // u nabla v n
227 M(i, j) -= .5 * dx *
228 fe.shape_grad_component(i, k, d1)[d2] *
229 n[d2] * u;
230 // u (nabla v)^T n
231 M(i, j) -= .5 * dx *
232 fe.shape_grad_component(i, k, d2)[d1] *
233 n[d2] * u;
234 }
235 }
236 }
237 }
238 }
239
254 template <int dim, typename number>
255 void
257 const FEValuesBase<dim> &fe,
258 const ArrayView<const std::vector<double>> &input,
259 const ArrayView<const std::vector<Tensor<1, dim>>> &Dinput,
260 const ArrayView<const std::vector<double>> &data,
261 double penalty,
262 double factor = 1.)
263 {
264 const unsigned int n_dofs = fe.dofs_per_cell;
268
269 for (unsigned int k = 0; k < fe.n_quadrature_points; ++k)
270 {
271 const double dx = factor * fe.JxW(k);
272 const Tensor<1, dim> n = fe.normal_vector(k);
273 for (unsigned int i = 0; i < n_dofs; ++i)
274 for (unsigned int d1 = 0; d1 < dim; ++d1)
275 {
276 const double u = input[d1][k];
277 const double v = fe.shape_value_component(i, k, d1);
278 const double g = data[d1][k];
279 result(i) += dx * 2. * penalty * (u - g) * v;
280
281 for (unsigned int d2 = 0; d2 < dim; ++d2)
282 {
283 // v . nabla u n
284 result(i) -= .5 * dx * v * Dinput[d1][k][d2] * n[d2];
285 // v . (nabla u)^T n
286 result(i) -= .5 * dx * v * Dinput[d2][k][d1] * n[d2];
287 // u nabla v n
288 result(i) -= .5 * dx * (u - g) *
289 fe.shape_grad_component(i, k, d1)[d2] * n[d2];
290 // u (nabla v)^T n
291 result(i) -= .5 * dx * (u - g) *
292 fe.shape_grad_component(i, k, d2)[d1] * n[d2];
293 }
294 }
295 }
296 }
297
306 template <int dim, typename number>
307 inline void
309 Vector<number> &result,
310 const FEValuesBase<dim> &fe,
311 const ArrayView<const std::vector<double>> &input,
312 const ArrayView<const std::vector<Tensor<1, dim>>> &Dinput,
313 const ArrayView<const std::vector<double>> &data,
314 double penalty,
315 double factor = 1.)
316 {
317 const unsigned int n_dofs = fe.dofs_per_cell;
321
322 for (unsigned int k = 0; k < fe.n_quadrature_points; ++k)
323 {
324 const double dx = factor * fe.JxW(k);
325 const Tensor<1, dim> n = fe.normal_vector(k);
326 for (unsigned int i = 0; i < n_dofs; ++i)
327 {
328 double udotn = 0.;
329 double gdotn = 0.;
330 double vdotn = 0.;
331 double ngradun = 0.;
332 double ngradvn = 0.;
333
334 for (unsigned int d = 0; d < dim; ++d)
335 {
336 udotn += n[d] * input[d][k];
337 gdotn += n[d] * data[d][k];
338 vdotn += n[d] * fe.shape_value_component(i, k, d);
339 ngradun += n * Dinput[d][k] * n[d];
340 ngradvn += n * fe.shape_grad_component(i, k, d) * n[d];
341 }
342 for (unsigned int d1 = 0; d1 < dim; ++d1)
343 {
344 const double u = input[d1][k] - udotn * n[d1];
345 const double v =
346 fe.shape_value_component(i, k, d1) - vdotn * n[d1];
347 const double g = data[d1][k] - gdotn * n[d1];
348 result(i) += dx * 2. * penalty * (u - g) * v;
349 // Correct the gradients below and subtract normal component
350 result(i) += dx * (ngradun * v + ngradvn * (u - g));
351 for (unsigned int d2 = 0; d2 < dim; ++d2)
352 {
353 // v . nabla u n
354 result(i) -= .5 * dx * Dinput[d1][k][d2] * n[d2] * v;
355 // v (nabla u)^T n
356 result(i) -= .5 * dx * Dinput[d2][k][d1] * n[d2] * v;
357 // u nabla v n
358 result(i) -= .5 * dx * (u - g) *
359 fe.shape_grad_component(i, k, d1)[d2] *
360 n[d2];
361 // u (nabla v)^T n
362 result(i) -= .5 * dx * (u - g) *
363 fe.shape_grad_component(i, k, d2)[d1] *
364 n[d2];
365 }
366 }
367 }
368 }
369 }
370
384 template <int dim, typename number>
385 void
387 Vector<number> &result,
388 const FEValuesBase<dim> &fe,
389 const ArrayView<const std::vector<double>> &input,
390 const ArrayView<const std::vector<Tensor<1, dim>>> &Dinput,
391 double penalty,
392 double factor = 1.)
393 {
394 const unsigned int n_dofs = fe.dofs_per_cell;
397
398 for (unsigned int k = 0; k < fe.n_quadrature_points; ++k)
399 {
400 const double dx = factor * fe.JxW(k);
401 const Tensor<1, dim> n = fe.normal_vector(k);
402 for (unsigned int i = 0; i < n_dofs; ++i)
403 for (unsigned int d1 = 0; d1 < dim; ++d1)
404 {
405 const double u = input[d1][k];
406 const double v = fe.shape_value_component(i, k, d1);
407 result(i) += dx * 2. * penalty * u * v;
408
409 for (unsigned int d2 = 0; d2 < dim; ++d2)
410 {
411 // v . nabla u n
412 result(i) -= .5 * dx * v * Dinput[d1][k][d2] * n[d2];
413 // v . (nabla u)^T n
414 result(i) -= .5 * dx * v * Dinput[d2][k][d1] * n[d2];
415 // u nabla v n
416 result(i) -= .5 * dx * u *
417 fe.shape_grad_component(i, k, d1)[d2] * n[d2];
418 // u (nabla v)^T n
419 result(i) -= .5 * dx * u *
420 fe.shape_grad_component(i, k, d2)[d1] * n[d2];
421 }
422 }
423 }
424 }
425
429 template <int dim>
430 inline void
435 const FEValuesBase<dim> &fe1,
436 const FEValuesBase<dim> &fe2,
437 const double pen,
438 const double int_factor = 1.,
439 const double ext_factor = -1.)
440 {
441 const unsigned int n_dofs = fe1.dofs_per_cell;
442
443 AssertDimension(fe1.get_fe().n_components(), dim);
444 AssertDimension(fe2.get_fe().n_components(), dim);
445 AssertDimension(M11.m(), n_dofs);
446 AssertDimension(M11.n(), n_dofs);
447 AssertDimension(M12.m(), n_dofs);
448 AssertDimension(M12.n(), n_dofs);
449 AssertDimension(M21.m(), n_dofs);
450 AssertDimension(M21.n(), n_dofs);
451 AssertDimension(M22.m(), n_dofs);
452 AssertDimension(M22.n(), n_dofs);
453
454 const double nu1 = int_factor;
455 const double nu2 = (ext_factor < 0) ? int_factor : ext_factor;
456 const double penalty = .5 * pen * (nu1 + nu2);
457
458 for (unsigned int k = 0; k < fe1.n_quadrature_points; ++k)
459 {
460 const double dx = fe1.JxW(k);
461 const Tensor<1, dim> n = fe1.normal_vector(k);
462 for (unsigned int i = 0; i < n_dofs; ++i)
463 for (unsigned int j = 0; j < n_dofs; ++j)
464 for (unsigned int d1 = 0; d1 < dim; ++d1)
465 {
466 const double u1 = fe1.shape_value_component(j, k, d1);
467 const double u2 = fe2.shape_value_component(j, k, d1);
468 const double v1 = fe1.shape_value_component(i, k, d1);
469 const double v2 = fe2.shape_value_component(i, k, d1);
470
471 M11(i, j) += dx * penalty * u1 * v1;
472 M12(i, j) -= dx * penalty * u2 * v1;
473 M21(i, j) -= dx * penalty * u1 * v2;
474 M22(i, j) += dx * penalty * u2 * v2;
475
476 for (unsigned int d2 = 0; d2 < dim; ++d2)
477 {
478 // v . nabla u n
479 M11(i, j) -= .25 * dx * nu1 *
480 fe1.shape_grad_component(j, k, d1)[d2] *
481 n[d2] * v1;
482 M12(i, j) -= .25 * dx * nu2 *
483 fe2.shape_grad_component(j, k, d1)[d2] *
484 n[d2] * v1;
485 M21(i, j) += .25 * dx * nu1 *
486 fe1.shape_grad_component(j, k, d1)[d2] *
487 n[d2] * v2;
488 M22(i, j) += .25 * dx * nu2 *
489 fe2.shape_grad_component(j, k, d1)[d2] *
490 n[d2] * v2;
491 // v (nabla u)^T n
492 M11(i, j) -= .25 * dx * nu1 *
493 fe1.shape_grad_component(j, k, d2)[d1] *
494 n[d2] * v1;
495 M12(i, j) -= .25 * dx * nu2 *
496 fe2.shape_grad_component(j, k, d2)[d1] *
497 n[d2] * v1;
498 M21(i, j) += .25 * dx * nu1 *
499 fe1.shape_grad_component(j, k, d2)[d1] *
500 n[d2] * v2;
501 M22(i, j) += .25 * dx * nu2 *
502 fe2.shape_grad_component(j, k, d2)[d1] *
503 n[d2] * v2;
504 // u nabla v n
505 M11(i, j) -= .25 * dx * nu1 *
506 fe1.shape_grad_component(i, k, d1)[d2] *
507 n[d2] * u1;
508 M12(i, j) += .25 * dx * nu1 *
509 fe1.shape_grad_component(i, k, d1)[d2] *
510 n[d2] * u2;
511 M21(i, j) -= .25 * dx * nu2 *
512 fe2.shape_grad_component(i, k, d1)[d2] *
513 n[d2] * u1;
514 M22(i, j) += .25 * dx * nu2 *
515 fe2.shape_grad_component(i, k, d1)[d2] *
516 n[d2] * u2;
517 // u (nabla v)^T n
518 M11(i, j) -= .25 * dx * nu1 *
519 fe1.shape_grad_component(i, k, d2)[d1] *
520 n[d2] * u1;
521 M12(i, j) += .25 * dx * nu1 *
522 fe1.shape_grad_component(i, k, d2)[d1] *
523 n[d2] * u2;
524 M21(i, j) -= .25 * dx * nu2 *
525 fe2.shape_grad_component(i, k, d2)[d1] *
526 n[d2] * u1;
527 M22(i, j) += .25 * dx * nu2 *
528 fe2.shape_grad_component(i, k, d2)[d1] *
529 n[d2] * u2;
530 }
531 }
532 }
533 }
537 template <int dim, typename number>
538 void
540 Vector<number> &result2,
541 const FEValuesBase<dim> &fe1,
542 const FEValuesBase<dim> &fe2,
543 const ArrayView<const std::vector<double>> &input1,
544 const ArrayView<const std::vector<Tensor<1, dim>>> &Dinput1,
545 const ArrayView<const std::vector<double>> &input2,
546 const ArrayView<const std::vector<Tensor<1, dim>>> &Dinput2,
547 double pen,
548 double int_factor = 1.,
549 double ext_factor = -1.)
550 {
551 const unsigned int n1 = fe1.dofs_per_cell;
552
553 AssertDimension(fe1.get_fe().n_components(), dim);
554 AssertDimension(fe2.get_fe().n_components(), dim);
559
560 const double nu1 = int_factor;
561 const double nu2 = (ext_factor < 0) ? int_factor : ext_factor;
562 const double penalty = .5 * pen * (nu1 + nu2);
563
564
565 for (unsigned int k = 0; k < fe1.n_quadrature_points; ++k)
566 {
567 const double dx = fe1.JxW(k);
568 const Tensor<1, dim> n = fe1.normal_vector(k);
569
570 for (unsigned int i = 0; i < n1; ++i)
571 for (unsigned int d1 = 0; d1 < dim; ++d1)
572 {
573 const double v1 = fe1.shape_value_component(i, k, d1);
574 const double v2 = fe2.shape_value_component(i, k, d1);
575 const double u1 = input1[d1][k];
576 const double u2 = input2[d1][k];
577
578 result1(i) += dx * penalty * u1 * v1;
579 result1(i) -= dx * penalty * u2 * v1;
580 result2(i) -= dx * penalty * u1 * v2;
581 result2(i) += dx * penalty * u2 * v2;
582
583 for (unsigned int d2 = 0; d2 < dim; ++d2)
584 {
585 // v . nabla u n
586 result1(i) -=
587 .25 * dx *
588 (nu1 * Dinput1[d1][k][d2] + nu2 * Dinput2[d1][k][d2]) *
589 n[d2] * v1;
590 result2(i) +=
591 .25 * dx *
592 (nu1 * Dinput1[d1][k][d2] + nu2 * Dinput2[d1][k][d2]) *
593 n[d2] * v2;
594 // v . (nabla u)^T n
595 result1(i) -=
596 .25 * dx *
597 (nu1 * Dinput1[d2][k][d1] + nu2 * Dinput2[d2][k][d1]) *
598 n[d2] * v1;
599 result2(i) +=
600 .25 * dx *
601 (nu1 * Dinput1[d2][k][d1] + nu2 * Dinput2[d2][k][d1]) *
602 n[d2] * v2;
603 // u nabla v n
604 result1(i) -= .25 * dx * nu1 *
605 fe1.shape_grad_component(i, k, d1)[d2] *
606 n[d2] * (u1 - u2);
607 result2(i) -= .25 * dx * nu2 *
608 fe2.shape_grad_component(i, k, d1)[d2] *
609 n[d2] * (u1 - u2);
610 // u (nabla v)^T n
611 result1(i) -= .25 * dx * nu1 *
612 fe1.shape_grad_component(i, k, d2)[d1] *
613 n[d2] * (u1 - u2);
614 result2(i) -= .25 * dx * nu2 *
615 fe2.shape_grad_component(i, k, d2)[d1] *
616 n[d2] * (u1 - u2);
617 }
618 }
619 }
620 }
621 } // namespace Elasticity
622} // namespace LocalIntegrators
623
625
626#endif
const unsigned int dofs_per_cell
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
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
virtual size_type size() const override
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
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:177
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:431
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:256
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:386
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:539
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:308
void nitsche_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, double penalty, double factor=1.)
Definition elasticity.h:122
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:83
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const double factor=1.)
Definition elasticity.h:50
Library of integrals over cells and faces.
Definition advection.h:34