15#ifndef dealii_integrators_elasticity_h
16#define dealii_integrators_elasticity_h
52 const double factor = 1.)
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)
81 template <
int dim,
typename number>
96 for (
unsigned int k = 0; k < nq; ++k)
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)
103 result(i) += dx * .25 *
104 (input[d1][k][d2] + input[d2][k][d1]) *
135 const double dx = factor * fe.
JxW(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)
143 M(i, j) += dx * 2. * penalty * u * v;
144 for (
unsigned int d2 = 0; d2 < dim; ++d2)
190 const double dx = factor * fe.
JxW(k);
192 for (
unsigned int i = 0; i < n_dofs; ++i)
193 for (
unsigned int j = 0; j < n_dofs; ++j)
200 for (
unsigned int d = 0; d < dim; ++d)
207 for (
unsigned int d1 = 0; d1 < dim; ++d1)
213 M(i, j) += dx * 2. * penalty * u * v;
215 M(i, j) += dx * (ngradun * v + ngradvn * u);
216 for (
unsigned int d2 = 0; d2 < dim; ++d2)
254 template <
int dim,
typename number>
258 const ArrayView<
const std::vector<double>> &input,
271 const double dx = factor * fe.
JxW(k);
273 for (
unsigned int i = 0; i < n_dofs; ++i)
274 for (
unsigned int d1 = 0; d1 < dim; ++d1)
276 const double u = input[d1][k];
278 const double g =
data[d1][k];
279 result(i) += dx * 2. * penalty * (u - g) * v;
281 for (
unsigned int d2 = 0; d2 < dim; ++d2)
284 result(i) -= .5 * dx * v * Dinput[d1][k][d2] * n[d2];
286 result(i) -= .5 * dx * v * Dinput[d2][k][d1] * n[d2];
288 result(i) -= .5 * dx * (u - g) *
291 result(i) -= .5 * dx * (u - g) *
306 template <
int dim,
typename number>
311 const ArrayView<
const std::vector<double>> &input,
324 const double dx = factor * fe.
JxW(k);
326 for (
unsigned int i = 0; i < n_dofs; ++i)
334 for (
unsigned int d = 0; d < dim; ++d)
336 udotn += n[d] * input[d][k];
337 gdotn += n[d] *
data[d][k];
339 ngradun += n * Dinput[d][k] * n[d];
342 for (
unsigned int d1 = 0; d1 < dim; ++d1)
344 const double u = input[d1][k] - udotn * n[d1];
347 const double g =
data[d1][k] - gdotn * n[d1];
348 result(i) += dx * 2. * penalty * (u - g) * v;
350 result(i) += dx * (ngradun * v + ngradvn * (u - g));
351 for (
unsigned int d2 = 0; d2 < dim; ++d2)
354 result(i) -= .5 * dx * Dinput[d1][k][d2] * n[d2] * v;
356 result(i) -= .5 * dx * Dinput[d2][k][d1] * n[d2] * v;
358 result(i) -= .5 * dx * (u - g) *
362 result(i) -= .5 * dx * (u - g) *
384 template <
int dim,
typename number>
389 const ArrayView<
const std::vector<double>> &input,
400 const double dx = factor * fe.
JxW(k);
402 for (
unsigned int i = 0; i < n_dofs; ++i)
403 for (
unsigned int d1 = 0; d1 < dim; ++d1)
405 const double u = input[d1][k];
407 result(i) += dx * 2. * penalty * u * v;
409 for (
unsigned int d2 = 0; d2 < dim; ++d2)
412 result(i) -= .5 * dx * v * Dinput[d1][k][d2] * n[d2];
414 result(i) -= .5 * dx * v * Dinput[d2][k][d1] * n[d2];
416 result(i) -= .5 * dx * u *
419 result(i) -= .5 * dx * u *
438 const double int_factor = 1.,
439 const double ext_factor = -1.)
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);
460 const double dx = fe1.
JxW(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)
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;
476 for (
unsigned int d2 = 0; d2 < dim; ++d2)
479 M11(i, j) -= .25 * dx * nu1 *
482 M12(i, j) -= .25 * dx * nu2 *
485 M21(i, j) += .25 * dx * nu1 *
488 M22(i, j) += .25 * dx * nu2 *
492 M11(i, j) -= .25 * dx * nu1 *
495 M12(i, j) -= .25 * dx * nu2 *
498 M21(i, j) += .25 * dx * nu1 *
501 M22(i, j) += .25 * dx * nu2 *
505 M11(i, j) -= .25 * dx * nu1 *
508 M12(i, j) += .25 * dx * nu1 *
511 M21(i, j) -= .25 * dx * nu2 *
514 M22(i, j) += .25 * dx * nu2 *
518 M11(i, j) -= .25 * dx * nu1 *
521 M12(i, j) += .25 * dx * nu1 *
524 M21(i, j) -= .25 * dx * nu2 *
527 M22(i, j) += .25 * dx * nu2 *
537 template <
int dim,
typename number>
543 const ArrayView<
const std::vector<double>> &input1,
545 const ArrayView<
const std::vector<double>> &input2,
548 double int_factor = 1.,
549 double ext_factor = -1.)
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);
567 const double dx = fe1.
JxW(k);
570 for (
unsigned int i = 0; i < n1; ++i)
571 for (
unsigned int d1 = 0; d1 < dim; ++d1)
575 const double u1 = input1[d1][k];
576 const double u2 = input2[d1][k];
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;
583 for (
unsigned int d2 = 0; d2 < dim; ++d2)
588 (nu1 * Dinput1[d1][k][d2] + nu2 * Dinput2[d1][k][d2]) *
592 (nu1 * Dinput1[d1][k][d2] + nu2 * Dinput2[d1][k][d2]) *
597 (nu1 * Dinput1[d2][k][d1] + nu2 * Dinput2[d2][k][d1]) *
601 (nu1 * Dinput1[d2][k][d1] + nu2 * Dinput2[d2][k][d1]) *
604 result1(i) -= .25 * dx * nu1 *
607 result2(i) -= .25 * dx * nu2 *
611 result1(i) -= .25 * dx * nu1 *
614 result2(i) -= .25 * dx * nu2 *
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
virtual size_type size() const override
#define DEAL_II_DEPRECATED
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
#define Assert(cond, exc)
#define AssertVectorVectorDimension(VEC, DIM1, DIM2)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
std::vector< index_type > data
void nitsche_tangential_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, double penalty, double factor=1.)
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.)
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.)
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_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.)
void nitsche_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, double penalty, double factor=1.)
void cell_residual(Vector< number > &result, const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim > > > &input, double factor=1.)
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const double factor=1.)
Library of integrals over cells and faces.