16 #ifndef dealii_integrators_elasticity_h 17 #define dealii_integrators_elasticity_h 20 #include <deal.II/base/config.h> 22 #include <deal.II/base/exceptions.h> 23 #include <deal.II/base/quadrature.h> 25 #include <deal.II/fe/fe_values.h> 26 #include <deal.II/fe/mapping.h> 28 #include <deal.II/lac/full_matrix.h> 30 #include <deal.II/meshworker/dof_info.h> 32 DEAL_II_NAMESPACE_OPEN
55 const double factor = 1.)
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)
84 template <
int dim,
typename number>
99 for (
unsigned int k = 0; k < nq; ++k)
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)
106 result(i) += dx * .25 *
107 (input[d1][k][d2] + input[d2][k][d1]) *
138 const double dx = factor * fe.
JxW(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)
146 M(i, j) += dx * 2. * penalty * u * v;
147 for (
unsigned int d2 = 0; d2 < dim; ++d2)
193 const double dx = factor * fe.
JxW(k);
195 for (
unsigned int i = 0; i < n_dofs; ++i)
196 for (
unsigned int j = 0; j < n_dofs; ++j)
203 for (
unsigned int d = 0; d < dim; ++d)
210 for (
unsigned int d1 = 0; d1 < dim; ++d1)
216 M(i, j) += dx * 2. * penalty * u * v;
218 M(i, j) += dx * (ngradun * v + ngradvn * u);
219 for (
unsigned int d2 = 0; d2 < dim; ++d2)
260 template <
int dim,
typename number>
264 const ArrayView<
const std::vector<double>> & input,
266 const ArrayView<
const std::vector<double>> & data,
277 const double dx = factor * fe.
JxW(k);
279 for (
unsigned int i = 0; i < n_dofs; ++i)
280 for (
unsigned int d1 = 0; d1 < dim; ++d1)
282 const double u = input[d1][k];
284 const double g = data[d1][k];
285 result(i) += dx * 2. * penalty * (u - g) * v;
287 for (
unsigned int d2 = 0; d2 < dim; ++d2)
290 result(i) -= .5 * dx * v * Dinput[d1][k][d2] * n[d2];
292 result(i) -= .5 * dx * v * Dinput[d2][k][d1] * n[d2];
294 result(i) -= .5 * dx * (u - g) *
297 result(i) -= .5 * dx * (u - g) *
312 template <
int dim,
typename number>
317 const ArrayView<
const std::vector<double>> & input,
319 const ArrayView<
const std::vector<double>> & data,
330 const double dx = factor * fe.
JxW(k);
332 for (
unsigned int i = 0; i < n_dofs; ++i)
340 for (
unsigned int d = 0; d < dim; ++d)
342 udotn += n[d] * input[d][k];
343 gdotn += n[d] * data[d][k];
345 ngradun += n * Dinput[d][k] * n[d];
348 for (
unsigned int d1 = 0; d1 < dim; ++d1)
350 const double u = input[d1][k] - udotn * n[d1];
353 const double g = data[d1][k] - gdotn * n[d1];
354 result(i) += dx * 2. * penalty * (u - g) * v;
356 result(i) += dx * (ngradun * v + ngradvn * (u - g));
357 for (
unsigned int d2 = 0; d2 < dim; ++d2)
360 result(i) -= .5 * dx * Dinput[d1][k][d2] * n[d2] * v;
362 result(i) -= .5 * dx * Dinput[d2][k][d1] * n[d2] * v;
364 result(i) -= .5 * dx * (u - g) *
368 result(i) -= .5 * dx * (u - g) *
393 template <
int dim,
typename number>
398 const ArrayView<
const std::vector<double>> & input,
409 const double dx = factor * fe.
JxW(k);
411 for (
unsigned int i = 0; i < n_dofs; ++i)
412 for (
unsigned int d1 = 0; d1 < dim; ++d1)
414 const double u = input[d1][k];
416 result(i) += dx * 2. * penalty * u * v;
418 for (
unsigned int d2 = 0; d2 < dim; ++d2)
421 result(i) -= .5 * dx * v * Dinput[d1][k][d2] * n[d2];
423 result(i) -= .5 * dx * v * Dinput[d2][k][d1] * n[d2];
425 result(i) -= .5 * dx * u *
428 result(i) -= .5 * dx * u *
447 const double int_factor = 1.,
448 const double ext_factor = -1.)
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);
469 const double dx = fe1.
JxW(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)
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;
485 for (
unsigned int d2 = 0; d2 < dim; ++d2)
488 M11(i, j) -= .25 * dx * nu1 *
491 M12(i, j) -= .25 * dx * nu2 *
494 M21(i, j) += .25 * dx * nu1 *
497 M22(i, j) += .25 * dx * nu2 *
501 M11(i, j) -= .25 * dx * nu1 *
504 M12(i, j) -= .25 * dx * nu2 *
507 M21(i, j) += .25 * dx * nu1 *
510 M22(i, j) += .25 * dx * nu2 *
514 M11(i, j) -= .25 * dx * nu1 *
517 M12(i, j) += .25 * dx * nu1 *
520 M21(i, j) -= .25 * dx * nu2 *
523 M22(i, j) += .25 * dx * nu2 *
527 M11(i, j) -= .25 * dx * nu1 *
530 M12(i, j) += .25 * dx * nu1 *
533 M21(i, j) -= .25 * dx * nu2 *
536 M22(i, j) += .25 * dx * nu2 *
549 template <
int dim,
typename number>
555 const ArrayView<
const std::vector<double>> & input1,
557 const ArrayView<
const std::vector<double>> & input2,
560 double int_factor = 1.,
561 double ext_factor = -1.)
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);
579 const double dx = fe1.
JxW(k);
582 for (
unsigned int i = 0; i < n1; ++i)
583 for (
unsigned int d1 = 0; d1 < dim; ++d1)
587 const double u1 = input1[d1][k];
588 const double u2 = input2[d1][k];
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;
595 for (
unsigned int d2 = 0; d2 < dim; ++d2)
600 (nu1 * Dinput1[d1][k][d2] + nu2 * Dinput2[d1][k][d2]) *
604 (nu1 * Dinput1[d1][k][d2] + nu2 * Dinput2[d1][k][d2]) *
609 (nu1 * Dinput1[d2][k][d1] + nu2 * Dinput2[d2][k][d1]) *
613 (nu1 * Dinput1[d2][k][d1] + nu2 * Dinput2[d2][k][d1]) *
616 result1(i) -= .25 * dx * nu1 *
619 result2(i) -= .25 * dx * nu2 *
623 result1(i) -= .25 * dx * nu1 *
626 result2(i) -= .25 * dx * nu2 *
636 DEAL_II_NAMESPACE_CLOSE
void cell_residual(Vector< number > &result, const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim >>> &input, double factor=1.)
#define AssertDimension(dim1, dim2)
const unsigned int dofs_per_cell
const FiniteElement< dim, spacedim > & get_fe() const
#define AssertVectorVectorDimension(VEC, DIM1, DIM2)
void nitsche_tangential_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, double penalty, double factor=1.)
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.)
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.
#define Assert(cond, exc)
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.)
const unsigned int n_quadrature_points
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.)
void nitsche_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, double penalty, double factor=1.)
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
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.)
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.)
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.)