17 #ifndef __deal2__integrators_laplace_h
18 #define __deal2__integrators_laplace_h
21 #include <deal.II/base/config.h>
23 #include <deal.II/base/quadrature.h>
24 #include <deal.II/lac/full_matrix.h>
25 #include <deal.II/fe/mapping.h>
26 #include <deal.II/fe/fe_values.h>
27 #include <deal.II/meshworker/dof_info.h>
29 DEAL_II_NAMESPACE_OPEN
31 namespace LocalIntegrators
59 const double factor = 1.)
66 const double dx = fe.
JxW(k) * factor;
67 for (
unsigned int i=0; i<n_dofs; ++i)
69 for (
unsigned int j=0; j<n_dofs; ++j)
70 for (
unsigned int d=0; d<n_components; ++d)
97 for (
unsigned int k=0; k<nq; ++k)
99 const double dx = factor * fe.
JxW(k);
100 for (
unsigned int i=0; i<n_dofs; ++i)
101 result(i) += dx * (input[k] * fe.
shape_grad(i,k));
128 for (
unsigned int k=0; k<nq; ++k)
130 const double dx = factor * fe.
JxW(k);
131 for (
unsigned int i=0; i<n_dofs; ++i)
132 for (
unsigned int d=0; d<n_comp; ++d)
168 const double dx = fe.
JxW(k) * factor;
170 for (
unsigned int i=0; i<n_dofs; ++i)
171 for (
unsigned int j=0; j<n_dofs; ++j)
172 for (
unsigned int d=0; d<n_comp; ++d)
200 const std::vector<double> &input,
202 const std::vector<double> &data,
213 const double dx = factor * fe.
JxW(k);
215 for (
unsigned int i=0; i<n_dofs; ++i)
218 const double dnu = Dinput[k] * n;
220 const double u= input[k];
221 const double g= data[k];
223 result(i) += dx*(2.*penalty*(u-g)*v - dnv*(u-g) - dnu*v);
251 const VectorSlice<
const std::vector<std::vector<double> > > &input,
253 const VectorSlice<
const std::vector<std::vector<double> > > &data,
265 const double dx = factor * fe.
JxW(k);
267 for (
unsigned int i=0; i<n_dofs; ++i)
268 for (
unsigned int d=0; d<n_comp; ++d)
271 const double dnu = Dinput[d][k] * n;
273 const double u= input[d][k];
274 const double g= data[d][k];
276 result(i) += dx*(2.*penalty*(u-g)*v - dnv*(u-g) - dnu*v);
310 double factor2 = -1.)
322 const double nui = factor1;
323 const double nue = (factor2 < 0) ? factor1 : factor2;
324 const double nu = .5*(nui+nue);
328 const double dx = fe1.
JxW(k);
332 for (
unsigned int i=0; i<n_dofs; ++i)
334 for (
unsigned int j=0; j<n_dofs; ++j)
344 M11(i,j) += dx*(-.5*nui*dnvi*ui-.5*nui*dnui*vi+nu*penalty*ui*vi);
345 M12(i,j) += dx*( .5*nui*dnvi*ue-.5*nue*dnue*vi-nu*penalty*vi*ue);
346 M21(i,j) += dx*(-.5*nue*dnve*ui+.5*nui*dnui*ve-nu*penalty*ui*ve);
347 M22(i,j) += dx*( .5*nue*dnve*ue+.5*nue*dnue*ve+nu*penalty*ue*ve);
378 double factor2 = -1.)
392 const double nui = factor1;
393 const double nue = (factor2 < 0) ? factor1 : factor2;
394 const double nu = .5*(nui+nue);
398 const double dx = fe1.
JxW(k);
400 for (
unsigned int i=0; i<n_dofs; ++i)
402 for (
unsigned int j=0; j<n_dofs; ++j)
409 double ngradu1n = 0.;
410 double ngradv1n = 0.;
411 double ngradu2n = 0.;
412 double ngradv2n = 0.;
414 for (
unsigned int d=0; d<dim; ++d)
441 M11(i,j) += dx*(-.5*nui*dnvi*ui-.5*nui*dnui*vi+nu*penalty*ui*vi);
442 M12(i,j) += dx*( .5*nui*dnvi*ue-.5*nue*dnue*vi-nu*penalty*vi*ue);
443 M21(i,j) += dx*(-.5*nue*dnve*ui+.5*nui*dnui*ve-nu*penalty*ui*ve);
444 M22(i,j) += dx*( .5*nue*dnve*ue+.5*nue*dnue*ve+nu*penalty*ue*ve);
468 const std::vector<double> &input1,
470 const std::vector<double> &input2,
473 double int_factor = 1.,
474 double ext_factor = -1.)
481 const double nui = int_factor;
482 const double nue = (ext_factor < 0) ? int_factor : ext_factor;
483 const double penalty = .5 * pen * (nui + nue);
489 const double dx = fe1.
JxW(k);
492 for (
unsigned int i=0; i<n_dofs; ++i)
496 const double dnvi = Dvi * n;
499 const double dnve = Dve * n;
501 const double ui = input1[k];
503 const double dnui = Dui * n;
504 const double ue = input2[k];
506 const double dnue = Due * n;
508 result1(i) += dx*(-.5*nui*dnvi*ui-.5*nui*dnui*vi+penalty*ui*vi);
509 result1(i) += dx*( .5*nui*dnvi*ue-.5*nue*dnue*vi-penalty*vi*ue);
510 result2(i) += dx*(-.5*nue*dnve*ui+.5*nui*dnui*ve-penalty*ui*ve);
511 result2(i) += dx*( .5*nue*dnve*ue+.5*nue*dnue*ve+penalty*ue*ve);
535 const VectorSlice<
const std::vector<std::vector<double> > > &input1,
537 const VectorSlice<
const std::vector<std::vector<double> > > &input2,
540 double int_factor = 1.,
541 double ext_factor = -1.)
551 const double nui = int_factor;
552 const double nue = (ext_factor < 0) ? int_factor : ext_factor;
553 const double penalty = .5 * pen * (nui + nue);
558 const double dx = fe1.
JxW(k);
561 for (
unsigned int i=0; i<n1; ++i)
562 for (
unsigned int d=0; d<n_comp; ++d)
566 const double dnvi = Dvi * n;
569 const double dnve = Dve * n;
571 const double ui = input1[d][k];
573 const double dnui = Dui * n;
574 const double ue = input2[d][k];
576 const double dnue = Due * n;
578 result1(i) += dx*(-.5*nui*dnvi*ui-.5*nui*dnui*vi+penalty*ui*vi);
579 result1(i) += dx*( .5*nui*dnvi*ue-.5*nue*dnue*vi-penalty*vi*ue);
580 result2(i) += dx*(-.5*nue*dnve*ui+.5*nui*dnui*ve-penalty*ui*ve);
581 result2(i) += dx*( .5*nue*dnve*ue+.5*nue*dnue*ve+penalty*ue*ve);
602 template <
int dim,
int spacedim,
typename number>
611 const unsigned int deg1sq = (deg1 == 0) ? 1 : deg1 * (deg1+1);
612 const unsigned int deg2sq = (deg2 == 0) ? 1 : deg2 * (deg2+1);
614 double penalty1 = deg1sq / dinfo1.
cell->extent_in_direction(normal1);
615 double penalty2 = deg2sq / dinfo2.
cell->extent_in_direction(normal2);
616 if (dinfo1.
cell->has_children() ^ dinfo2.
cell->has_children())
622 const double penalty = 0.5*(penalty1 + penalty2);
629 DEAL_II_NAMESPACE_CLOSE
Triangulation< dim, spacedim >::cell_iterator cell
The current cell.
#define AssertDimension(dim1, dim2)
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.)
const unsigned int n_quadrature_points
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const double factor=1.)
const double & shape_value(const unsigned int function_no, const unsigned int point_no) const
#define AssertVectorVectorDimension(vec, dim1, dim2)
const unsigned int dofs_per_cell
void cell_residual(Vector< double > &result, const FEValuesBase< dim > &fe, const std::vector< Tensor< 1, dim > > &input, double factor=1.)
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_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.)
const Tensor< 1, spacedim > & shape_grad(const unsigned int function_no, const unsigned int quadrature_point) const
const FiniteElement< dim, spacedim > & get_fe() const
#define Assert(cond, exc)
double shape_value_component(const unsigned int function_no, const unsigned int point_no, const unsigned int component) const
const Point< spacedim > & normal_vector(const unsigned int i) const
unsigned int n_components() const
void nitsche_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, double penalty, double factor=1.)
double JxW(const unsigned int quadrature_point) const
double compute_penalty(const MeshWorker::DoFInfo< dim, spacedim, number > &dinfo1, const MeshWorker::DoFInfo< dim, spacedim, number > &dinfo2, unsigned int deg1, unsigned int deg2)
::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
Tensor< 1, spacedim > shape_grad_component(const unsigned int function_no, const unsigned int point_no, const unsigned int component) const
::ExceptionBase & ExcInternalError()
Triangulation< dim, spacedim >::face_iterator face
The current face.
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.)