16 #ifndef dealii_integrators_laplace_h
17 #define dealii_integrators_laplace_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
59 const double factor = 1.)
62 const unsigned int n_components = fe.
get_fe().n_components();
66 const double dx = fe.
JxW(k) * factor;
67 for (
unsigned int i = 0; i < n_dofs; ++i)
70 for (
unsigned int d = 0; d < n_components; ++d)
76 for (
unsigned int j = i + 1; j < n_dofs; ++j)
79 for (
unsigned int d = 0; d < n_components; ++d)
108 for (
unsigned int k = 0; k < nq; ++k)
110 const double dx = factor * fe.
JxW(k);
111 for (
unsigned int i = 0; i < n_dofs; ++i)
131 const unsigned int n_comp = fe.
get_fe().n_components();
137 for (
unsigned int k = 0; k < nq; ++k)
139 const double dx = factor * fe.
JxW(k);
140 for (
unsigned int i = 0; i < n_dofs; ++i)
141 for (
unsigned int d = 0; d < n_comp; ++d)
171 const unsigned int n_comp = fe.
get_fe().n_components();
178 const double dx = fe.
JxW(k) * factor;
180 for (
unsigned int i = 0; i < n_dofs; ++i)
181 for (
unsigned int j = 0; j < n_dofs; ++j)
182 for (
unsigned int d = 0; d < n_comp; ++d)
221 const double dx = fe.
JxW(k) * factor;
223 for (
unsigned int i = 0; i < n_dofs; ++i)
224 for (
unsigned int j = 0; j < n_dofs; ++j)
231 for (
unsigned int d = 0; d < dim; ++d)
239 for (
unsigned int d = 0; d < dim; ++d)
250 M(i, j) +=
dx * (2. * penalty * u_t * v_t - dnu_t * v_t -
277 const std::vector<double> & input,
279 const std::vector<double> & data,
290 const double dx = factor * fe.
JxW(k);
292 for (
unsigned int i = 0; i < n_dofs; ++i)
295 const double dnu = Dinput[k] * n;
297 const double u = input[k];
298 const double g = data[k];
301 dx * (2. * penalty * (u - g) * v - dnv * (u - g) - dnu * v);
327 const ArrayView<
const std::vector<double>> & input,
329 const ArrayView<
const std::vector<double>> & data,
334 const unsigned int n_comp = fe.
get_fe().n_components();
341 const double dx = factor * fe.
JxW(k);
343 for (
unsigned int i = 0; i < n_dofs; ++i)
344 for (
unsigned int d = 0; d < n_comp; ++d)
347 const double dnu = Dinput[d][k] * n;
349 const double u = input[d][k];
350 const double g = data[d][k];
353 dx * (2. * penalty * (u - g) * v - dnv * (u - g) - dnu * v);
387 double factor2 = -1.)
399 const double nui = factor1;
400 const double nue = (factor2 < 0) ? factor1 : factor2;
401 const double nu = .5 * (nui + nue);
405 const double dx = fe1.
JxW(k);
407 for (
unsigned int d = 0; d < fe1.
get_fe().n_components(); ++d)
409 for (
unsigned int i = 0; i < n_dofs; ++i)
411 for (
unsigned int j = 0; j < n_dofs; ++j)
422 dx * (-.5 * nui * dnvi * ui - .5 * nui * dnui * vi +
423 nu * penalty * ui * vi);
425 dx * (.5 * nui * dnvi * ue - .5 * nue * dnue * vi -
426 nu * penalty * vi * ue);
428 dx * (-.5 * nue * dnve * ui + .5 * nui * dnui * ve -
429 nu * penalty * ui * ve);
431 dx * (.5 * nue * dnve * ue + .5 * nue * dnue * ve +
432 nu * penalty * ue * ve);
463 double factor2 = -1.)
477 const double nui = factor1;
478 const double nue = (factor2 < 0) ? factor1 : factor2;
479 const double nu = .5 * (nui + nue);
483 const double dx = fe1.
JxW(k);
485 for (
unsigned int i = 0; i < n_dofs; ++i)
490 for (
unsigned int j = 0; j < n_dofs; ++j)
497 double ngradu1n = 0.;
498 double ngradv1n = 0.;
499 double ngradu2n = 0.;
500 double ngradv2n = 0.;
502 for (
unsigned int d = 0; d < dim; ++d)
518 for (
unsigned int d = 0; d < dim; ++d)
541 dx * (-.5 * nui * dnvi * ui - .5 * nui * dnui * vi +
542 nu * penalty * ui * vi);
544 dx * (.5 * nui * dnvi * ue - .5 * nue * dnue * vi -
545 nu * penalty * vi * ue);
547 dx * (-.5 * nue * dnve * ui + .5 * nui * dnui * ve -
548 nu * penalty * ui * ve);
550 dx * (.5 * nue * dnve * ue + .5 * nue * dnue * ve +
551 nu * penalty * ue * ve);
574 const std::vector<double> & input1,
576 const std::vector<double> & input2,
579 double int_factor = 1.,
580 double ext_factor = -1.)
587 const double nui = int_factor;
588 const double nue = (ext_factor < 0) ? int_factor : ext_factor;
589 const double penalty = .5 * pen * (nui + nue);
595 const double dx = fe1.
JxW(k);
598 for (
unsigned int i = 0; i < n_dofs; ++i)
602 const double dnvi = Dvi * n;
605 const double dnve = Dve * n;
607 const double ui = input1[k];
609 const double dnui = Dui * n;
610 const double ue = input2[k];
612 const double dnue = Due * n;
614 result1(i) +=
dx * (-.5 * nui * dnvi * ui - .5 * nui * dnui * vi +
616 result1(i) +=
dx * (.5 * nui * dnvi * ue - .5 * nue * dnue * vi -
618 result2(i) +=
dx * (-.5 * nue * dnve * ui + .5 * nui * dnui * ve -
620 result2(i) +=
dx * (.5 * nue * dnve * ue + .5 * nue * dnue * ve +
644 const ArrayView<
const std::vector<double>> & input1,
646 const ArrayView<
const std::vector<double>> & input2,
649 double int_factor = 1.,
650 double ext_factor = -1.)
652 const unsigned int n_comp = fe1.
get_fe().n_components();
660 const double nui = int_factor;
661 const double nue = (ext_factor < 0) ? int_factor : ext_factor;
662 const double penalty = .5 * pen * (nui + nue);
667 const double dx = fe1.
JxW(k);
670 for (
unsigned int i = 0; i < n1; ++i)
671 for (
unsigned int d = 0; d < n_comp; ++d)
675 const double dnvi = Dvi * n;
678 const double dnve = Dve * n;
680 const double ui = input1[d][k];
682 const double dnui = Dui * n;
683 const double ue = input2[d][k];
685 const double dnue = Due * n;
687 result1(i) +=
dx * (-.5 * nui * dnvi * ui -
688 .5 * nui * dnui * vi + penalty * ui * vi);
689 result1(i) +=
dx * (.5 * nui * dnvi * ue -
690 .5 * nue * dnue * vi - penalty * vi * ue);
691 result2(i) +=
dx * (-.5 * nue * dnve * ui +
692 .5 * nui * dnui * ve - penalty * ui * ve);
693 result2(i) +=
dx * (.5 * nue * dnve * ue +
694 .5 * nue * dnue * ve + penalty * ue * ve);
715 template <
int dim,
int spacedim,
typename number>
722 const unsigned int normal1 =
724 const unsigned int normal2 =
726 const unsigned int deg1sq = (deg1 == 0) ? 1 : deg1 * (deg1 + 1);
727 const unsigned int deg2sq = (deg2 == 0) ? 1 : deg2 * (deg2 + 1);
729 double penalty1 = deg1sq / dinfo1.
cell->extent_in_direction(normal1);
730 double penalty2 = deg2sq / dinfo2.
cell->extent_in_direction(normal2);
731 if (dinfo1.
cell->has_children() ^ dinfo2.
cell->has_children())
737 const double penalty = 0.5 * (penalty1 + penalty2);
744 DEAL_II_NAMESPACE_CLOSE