15#ifndef dealii_integrators_laplace_h
16#define dealii_integrators_laplace_h
53 const double factor = 1.)
60 const double dx = fe.
JxW(k) * factor;
61 for (
unsigned int i = 0; i < n_dofs; ++i)
64 for (
unsigned int d = 0; d < n_components; ++d)
70 for (
unsigned int j = i + 1; j < n_dofs; ++j)
73 for (
unsigned int d = 0; d < n_components; ++d)
102 for (
unsigned int k = 0; k < nq; ++k)
104 const double dx = factor * fe.
JxW(k);
105 for (
unsigned int i = 0; i < n_dofs; ++i)
106 result(i) += dx * (input[k] * fe.
shape_grad(i, k));
131 for (
unsigned int k = 0; k < nq; ++k)
133 const double dx = factor * fe.
JxW(k);
134 for (
unsigned int i = 0; i < n_dofs; ++i)
135 for (
unsigned int d = 0; d < n_comp; ++d)
169 const double dx = fe.
JxW(k) * factor;
171 for (
unsigned int i = 0; i < n_dofs; ++i)
172 for (
unsigned int j = 0; j < n_dofs; ++j)
173 for (
unsigned int d = 0; d < n_comp; ++d)
209 const double dx = fe.
JxW(k) * factor;
211 for (
unsigned int i = 0; i < n_dofs; ++i)
212 for (
unsigned int j = 0; j < n_dofs; ++j)
219 for (
unsigned int d = 0; d < dim; ++d)
227 for (
unsigned int d = 0; d < dim; ++d)
238 M(i, j) += dx * (2. * penalty * u_t * v_t - dnu_t * v_t -
262 const std::vector<double> &input,
264 const std::vector<double> &
data,
275 const double dx = factor * fe.
JxW(k);
277 for (
unsigned int i = 0; i < n_dofs; ++i)
280 const double dnu = Dinput[k] * n;
282 const double u = input[k];
283 const double g =
data[k];
286 dx * (2. * penalty * (u - g) * v - dnv * (u - g) - dnu * v);
309 const ArrayView<
const std::vector<double>> &input,
323 const double dx = factor * fe.
JxW(k);
325 for (
unsigned int i = 0; i < n_dofs; ++i)
326 for (
unsigned int d = 0; d < n_comp; ++d)
329 const double dnu = Dinput[d][k] * n;
331 const double u = input[d][k];
332 const double g =
data[d][k];
335 dx * (2. * penalty * (u - g) * v - dnv * (u - g) - dnu * v);
366 double factor2 = -1.)
378 const double nui = factor1;
379 const double nue = (factor2 < 0) ? factor1 : factor2;
380 const double nu = .5 * (nui + nue);
384 const double dx = fe1.
JxW(k);
388 for (
unsigned int i = 0; i < n_dofs; ++i)
390 for (
unsigned int j = 0; j < n_dofs; ++j)
401 dx * (-.5 * nui * dnvi * ui - .5 * nui * dnui * vi +
402 nu * penalty * ui * vi);
404 dx * (.5 * nui * dnvi * ue - .5 * nue * dnue * vi -
405 nu * penalty * vi * ue);
407 dx * (-.5 * nue * dnve * ui + .5 * nui * dnui * ve -
408 nu * penalty * ui * ve);
410 dx * (.5 * nue * dnve * ue + .5 * nue * dnue * ve +
411 nu * penalty * ue * ve);
439 double factor2 = -1.)
453 const double nui = factor1;
454 const double nue = (factor2 < 0) ? factor1 : factor2;
455 const double nu = .5 * (nui + nue);
459 const double dx = fe1.
JxW(k);
461 for (
unsigned int i = 0; i < n_dofs; ++i)
466 for (
unsigned int j = 0; j < n_dofs; ++j)
473 double ngradu1n = 0.;
474 double ngradv1n = 0.;
475 double ngradu2n = 0.;
476 double ngradv2n = 0.;
478 for (
unsigned int d = 0; d < dim; ++d)
494 for (
unsigned int d = 0; d < dim; ++d)
517 dx * (-.5 * nui * dnvi * ui - .5 * nui * dnui * vi +
518 nu * penalty * ui * vi);
520 dx * (.5 * nui * dnvi * ue - .5 * nue * dnue * vi -
521 nu * penalty * vi * ue);
523 dx * (-.5 * nue * dnve * ui + .5 * nui * dnui * ve -
524 nu * penalty * ui * ve);
526 dx * (.5 * nue * dnve * ue + .5 * nue * dnue * ve +
527 nu * penalty * ue * ve);
547 const std::vector<double> &input1,
549 const std::vector<double> &input2,
552 double int_factor = 1.,
553 double ext_factor = -1.)
560 const double nui = int_factor;
561 const double nue = (ext_factor < 0) ? int_factor : ext_factor;
562 const double penalty = .5 * pen * (nui + nue);
568 const double dx = fe1.
JxW(k);
571 for (
unsigned int i = 0; i < n_dofs; ++i)
575 const double dnvi = Dvi * n;
578 const double dnve = Dve * n;
580 const double ui = input1[k];
582 const double dnui = Dui * n;
583 const double ue = input2[k];
585 const double dnue = Due * n;
587 result1(i) += dx * (-.5 * nui * dnvi * ui - .5 * nui * dnui * vi +
589 result1(i) += dx * (.5 * nui * dnvi * ue - .5 * nue * dnue * vi -
591 result2(i) += dx * (-.5 * nue * dnve * ui + .5 * nui * dnui * ve -
593 result2(i) += dx * (.5 * nue * dnve * ue + .5 * nue * dnue * ve +
614 const ArrayView<
const std::vector<double>> &input1,
616 const ArrayView<
const std::vector<double>> &input2,
619 double int_factor = 1.,
620 double ext_factor = -1.)
630 const double nui = int_factor;
631 const double nue = (ext_factor < 0) ? int_factor : ext_factor;
632 const double penalty = .5 * pen * (nui + nue);
637 const double dx = fe1.
JxW(k);
640 for (
unsigned int i = 0; i < n1; ++i)
641 for (
unsigned int d = 0; d < n_comp; ++d)
645 const double dnvi = Dvi * n;
648 const double dnve = Dve * n;
650 const double ui = input1[d][k];
652 const double dnui = Dui * n;
653 const double ue = input2[d][k];
655 const double dnue = Due * n;
657 result1(i) += dx * (-.5 * nui * dnvi * ui -
658 .5 * nui * dnui * vi + penalty * ui * vi);
659 result1(i) += dx * (.5 * nui * dnvi * ue -
660 .5 * nue * dnue * vi - penalty * vi * ue);
661 result2(i) += dx * (-.5 * nue * dnve * ui +
662 .5 * nui * dnui * ve - penalty * ui * ve);
663 result2(i) += dx * (.5 * nue * dnve * ue +
664 .5 * nue * dnue * ve + penalty * ue * ve);
682 template <
int dim,
int spacedim,
typename number>
689 const unsigned int normal1 =
691 const unsigned int normal2 =
693 const unsigned int deg1sq = (deg1 == 0) ? 1 : deg1 * (deg1 + 1);
694 const unsigned int deg2sq = (deg2 == 0) ? 1 : deg2 * (deg2 + 1);
696 double penalty1 = deg1sq / dinfo1.
cell->extent_in_direction(normal1);
697 double penalty2 = deg2sq / dinfo2.
cell->extent_in_direction(normal2);
698 if (dinfo1.
cell->has_children() ^ dinfo2.
cell->has_children())
704 const double penalty = 0.5 * (penalty1 + penalty2);
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 Tensor< 1, spacedim > & shape_grad(const unsigned int i, const unsigned int q_point) const
const FiniteElement< dim, spacedim > & get_fe() const
double JxW(const unsigned int q_point) const
const double & shape_value(const unsigned int i, const unsigned int q_point) const
unsigned int n_components() const
Triangulation< dim, spacedim >::face_iterator face
The current face.
Triangulation< dim, spacedim >::cell_iterator cell
The current cell.
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 & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
std::vector< index_type > data
void ip_tangential_matrix(FullMatrix< double > &M11, FullMatrix< double > &M12, FullMatrix< double > &M21, FullMatrix< double > &M22, const FEValuesBase< dim > &fe1, const FEValuesBase< dim > &fe2, double penalty, double factor1=1., double factor2=-1.)
void cell_residual(Vector< double > &result, const FEValuesBase< dim > &fe, const std::vector< Tensor< 1, dim > > &input, double factor=1.)
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const double factor=1.)
double compute_penalty(const MeshWorker::DoFInfo< dim, spacedim, number > &dinfo1, const MeshWorker::DoFInfo< dim, spacedim, number > &dinfo2, unsigned int deg1, unsigned int deg2)
void nitsche_residual(Vector< double > &result, const FEValuesBase< dim > &fe, const std::vector< double > &input, const std::vector< Tensor< 1, dim > > &Dinput, const std::vector< double > &data, 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, double penalty, double factor1=1., double factor2=-1.)
void nitsche_tangential_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, double penalty, double factor=1.)
void ip_residual(Vector< double > &result1, Vector< double > &result2, const FEValuesBase< dim > &fe1, const FEValuesBase< dim > &fe2, const std::vector< double > &input1, const std::vector< Tensor< 1, dim > > &Dinput1, const std::vector< double > &input2, 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.)
Library of integrals over cells and faces.