16 #ifndef dealii_integrators_grad_div_h 17 #define dealii_integrators_grad_div_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
69 const double dx = factor * fe.
JxW(k);
70 for (
unsigned int i = 0; i < n_dofs; ++i)
71 for (
unsigned int j = 0; j < n_dofs; ++j)
78 M(i, j) += dx * divu * divv;
92 template <
int dim,
typename number>
98 const double factor = 1.)
107 const double dx = factor * fetest.
JxW(k);
108 for (
unsigned int i = 0; i < n_dofs; ++i)
113 for (
unsigned int d = 0; d < dim; ++d)
114 du += input[d][k][d];
116 result(i) += dx * du * divv;
144 const double dx = factor * fe.
JxW(k);
146 for (
unsigned int i = 0; i < n_dofs; ++i)
147 for (
unsigned int j = 0; j < n_dofs; ++j)
153 double un = 0., vn = 0.;
154 for (
unsigned int d = 0; d < dim; ++d)
160 M(i, j) += dx * 2. * penalty * un * vn;
161 M(i, j) -= dx * (divu * vn + divv * un);
189 const VectorSlice<
const std::vector<std::vector<double>>> & input,
191 const VectorSlice<
const std::vector<std::vector<double>>> & data,
203 const double dx = factor * fe.
JxW(k);
208 for (
unsigned int d = 0; d < dim; ++d)
210 umgn += (input[d][k] - data[d][k]) * n[d];
211 divu += Dinput[d][k][d];
214 for (
unsigned int i = 0; i < n_dofs; ++i)
219 for (
unsigned int d = 0; d < dim; ++d)
223 dx * (2. * penalty * umgn * vn - divv * umgn - divu * vn);
246 double factor2 = -1.)
258 const double fi = factor1;
259 const double fe = (factor2 < 0) ? factor1 : factor2;
260 const double f = .5 * (fi + fe);
264 const double dx = fe1.
JxW(k);
266 for (
unsigned int i = 0; i < n_dofs; ++i)
267 for (
unsigned int j = 0; j < n_dofs; ++j)
282 for (
unsigned int d = 0; d < dim; ++d)
290 dx * (-.5 * fi * divvi * uni - .5 * fi * divui * vni +
291 f * penalty * uni * vni);
293 dx * (.5 * fi * divvi * une - .5 * fe * divue * vni -
294 f * penalty * vni * une);
296 dx * (-.5 * fe * divve * uni + .5 * fi * divui * vne -
297 f * penalty * uni * vne);
299 dx * (.5 * fe * divve * une + .5 * fe * divue * vne +
300 f * penalty * une * vne);
326 const VectorSlice<
const std::vector<std::vector<double>>> &input1,
329 const VectorSlice<
const std::vector<std::vector<double>>> &input2,
333 double int_factor = 1.,
334 double ext_factor = -1.)
344 const double fi = int_factor;
345 const double fe = (ext_factor < 0) ? int_factor : ext_factor;
346 const double penalty = .5 * pen * (fi + fe);
351 const double dx = fe1.
JxW(k);
357 for (
unsigned int d = 0; d < dim; ++d)
359 uni += input1[d][k] * n[d];
360 une += input2[d][k] * n[d];
361 divui += Dinput1[d][k][d];
362 divue += Dinput2[d][k][d];
365 for (
unsigned int i = 0; i < n1; ++i)
373 for (
unsigned int d = 0; d < dim; ++d)
379 result1(i) += dx * (-.5 * fi * divvi * uni -
380 .5 * fi * divui * vni + penalty * uni * vni);
381 result1(i) += dx * (.5 * fi * divvi * une -
382 .5 * fe * divue * vni - penalty * vni * une);
383 result2(i) += dx * (-.5 * fe * divve * uni +
384 .5 * fi * divui * vne - penalty * uni * vne);
385 result2(i) += dx * (.5 * fe * divve * une +
386 .5 * fe * divue * vne + penalty * une * vne);
393 DEAL_II_NAMESPACE_CLOSE
double shape_value_component(const unsigned int function_no, const unsigned int point_no, const unsigned int component) const
void nitsche_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, double penalty, double factor=1.)
#define AssertDimension(dim1, dim2)
const unsigned int dofs_per_cell
const Tensor< 1, spacedim > & normal_vector(const unsigned int i) const
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, 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.)
#define AssertVectorVectorDimension(vec, dim1, dim2)
Library of integrals over cells and faces.
void cell_residual(Vector< number > &result, const FEValuesBase< dim > &fetest, const VectorSlice< const std::vector< std::vector< Tensor< 1, dim >>>> &input, const double factor=1.)
void nitsche_residual(Vector< double > &result, const FEValuesBase< dim > &fe, const VectorSlice< const std::vector< std::vector< double >>> &input, const VectorSlice< const std::vector< std::vector< Tensor< 1, dim >>>> &Dinput, const VectorSlice< const std::vector< std::vector< double >>> &data, double penalty, double factor=1.)
const FiniteElement< dim, spacedim > & get_fe() const
const unsigned int n_quadrature_points
void ip_residual(Vector< double > &result1, Vector< double > &result2, const FEValuesBase< dim > &fe1, const FEValuesBase< dim > &fe2, const VectorSlice< const std::vector< std::vector< double >>> &input1, const VectorSlice< const std::vector< std::vector< Tensor< 1, dim >>>> &Dinput1, const VectorSlice< const std::vector< std::vector< double >>> &input2, const VectorSlice< const std::vector< std::vector< Tensor< 1, dim >>>> &Dinput2, double pen, double int_factor=1., double ext_factor=-1.)
double JxW(const unsigned int quadrature_point) const