16#ifndef dealii_integrators_advection_h
17#define dealii_integrators_advection_h
78 const ArrayView<
const std::vector<double>> &velocity,
79 const double factor = 1.)
90 const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
102 const double dx = factor * fe.
JxW(k);
103 const unsigned int vindex = k * v_increment;
105 for (
unsigned j = 0; j < n_dofs; ++j)
106 for (
unsigned i = 0; i < t_dofs; ++i)
107 for (
unsigned int c = 0; c < n_components; ++c)
111 for (
unsigned int d = 1;
d < dim; ++
d)
134 const ArrayView<
const std::vector<double>> &velocity,
144 const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
145 if (v_increment == 1)
150 for (
unsigned k = 0; k < nq; ++k)
152 const double dx = factor * fe.
JxW(k);
153 for (
unsigned i = 0; i < n_dofs; ++i)
154 for (
unsigned int d = 0;
d < dim; ++
d)
156 velocity[
d][k * v_increment];
177 const ArrayView<
const std::vector<double>> & velocity,
189 const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
190 if (v_increment == 1)
195 for (
unsigned k = 0; k < nq; ++k)
197 const double dx = factor * fe.
JxW(k);
198 for (
unsigned i = 0; i < n_dofs; ++i)
199 for (
unsigned int c = 0; c < n_comp; ++c)
200 for (
unsigned int d = 0;
d < dim; ++
d)
201 result(i) +=
dx * input[c][k][
d] *
203 velocity[
d][k * v_increment];
218 const std::vector<double> & input,
219 const ArrayView<
const std::vector<double>> &velocity,
229 const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
230 if (v_increment == 1)
235 for (
unsigned k = 0; k < nq; ++k)
237 const double dx = factor * fe.
JxW(k);
238 for (
unsigned i = 0; i < n_dofs; ++i)
239 for (
unsigned int d = 0;
d < dim; ++
d)
241 velocity[
d][k * v_increment];
258 const ArrayView<
const std::vector<double>> &input,
259 const ArrayView<
const std::vector<double>> &velocity,
271 const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
272 if (v_increment == 1)
277 for (
unsigned k = 0; k < nq; ++k)
279 const double dx = factor * fe.
JxW(k);
280 for (
unsigned i = 0; i < n_dofs; ++i)
281 for (
unsigned int c = 0; c < n_comp; ++c)
282 for (
unsigned int d = 0;
d < dim; ++
d)
283 result(i) -=
dx * input[c][k] *
285 velocity[
d][k * v_increment];
313 const ArrayView<
const std::vector<double>> &velocity,
323 const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
324 if (v_increment == 1)
331 const double dx = factor * fe.
JxW(k);
334 for (
unsigned int d = 0;
d < dim; ++
d)
339 for (
unsigned i = 0; i < t_dofs; ++i)
340 for (
unsigned j = 0; j < n_dofs; ++j)
346 for (
unsigned int c = 0; c < n_components; ++c)
385 const std::vector<double> & input,
386 const std::vector<double> & data,
387 const ArrayView<
const std::vector<double>> &velocity,
396 const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
397 if (v_increment == 1)
405 const double dx = factor * fe.
JxW(k);
408 for (
unsigned int d = 0;
d < dim; ++
d)
412 const double val = (nv > 0.) ? input[k] : -data[k];
414 for (
unsigned i = 0; i < n_dofs; ++i)
417 result(i) +=
dx * nv * val * v;
452 const ArrayView<
const std::vector<double>> &input,
453 const ArrayView<
const std::vector<double>> &data,
454 const ArrayView<
const std::vector<double>> &velocity,
464 const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
465 if (v_increment == 1)
473 const double dx = factor * fe.
JxW(k);
476 for (
unsigned int d = 0;
d < dim; ++
d)
479 std::vector<double> val(n_comp);
481 for (
unsigned int d = 0;
d < n_comp; ++
d)
483 val[
d] = (nv > 0.) ? input[
d][k] : -data[
d][k];
484 for (
unsigned i = 0; i < n_dofs; ++i)
487 result(i) +=
dx * nv * val[
d] * v;
525 const ArrayView<
const std::vector<double>> &velocity,
526 const double factor = 1.)
534 const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
535 if (v_increment == 1)
542 double nbeta = fe1.
normal_vector(k)[0] * velocity[0][k * v_increment];
543 for (
unsigned int d = 1;
d < dim; ++
d)
545 const double dx_nbeta = factor *
std::abs(nbeta) * fe1.
JxW(k);
551 for (
unsigned i = 0; i < n1; ++i)
552 for (
unsigned j = 0; j < n1; ++j)
563 for (
unsigned int d = 0;
d < fe1.
get_fe().n_components();
566 M1(i, j) += dx_nbeta *
569 M2(i, j) -= dx_nbeta *
606 const std::vector<double> & input1,
607 const std::vector<double> & input2,
608 const ArrayView<
const std::vector<double>> &velocity,
609 const double factor = 1.)
622 const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
623 if (v_increment == 1)
630 double nbeta = fe1.
normal_vector(k)[0] * velocity[0][k * v_increment];
631 for (
unsigned int d = 1;
d < dim; ++
d)
633 const double dx_nbeta = factor * nbeta * fe1.
JxW(k);
635 for (
unsigned i = 0; i < n1; ++i)
639 const double u1 = input1[k];
640 const double u2 = input2[k];
643 result1(i) += dx_nbeta * u1 *
v1;
644 result2(i) -= dx_nbeta * u1 * v2;
648 result1(i) += dx_nbeta * u2 *
v1;
649 result2(i) -= dx_nbeta * u2 * v2;
683 const ArrayView<
const std::vector<double>> &input1,
684 const ArrayView<
const std::vector<double>> &input2,
685 const ArrayView<
const std::vector<double>> &velocity,
686 const double factor = 1.)
698 const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
699 if (v_increment == 1)
706 double nbeta = fe1.
normal_vector(k)[0] * velocity[0][k * v_increment];
707 for (
unsigned int d = 1;
d < dim; ++
d)
709 const double dx_nbeta = factor * nbeta * fe1.
JxW(k);
711 for (
unsigned i = 0; i < n1; ++i)
712 for (
unsigned int d = 0;
d < n_comp; ++
d)
716 const double u1 = input1[
d][k];
717 const double u2 = input2[
d][k];
720 result1(i) += dx_nbeta * u1 *
v1;
721 result2(i) -= dx_nbeta * u1 * v2;
725 result1(i) += dx_nbeta * u2 *
v1;
726 result2(i) -= dx_nbeta * u2 * v2;
const unsigned int dofs_per_cell
const unsigned int n_quadrature_points
Tensor< 1, spacedim > shape_grad_component(const unsigned int function_no, const unsigned int point_no, const unsigned int component) const
double shape_value_component(const unsigned int function_no, const unsigned int point_no, const unsigned int component) const
const Tensor< 1, spacedim > & normal_vector(const unsigned int i) const
double JxW(const unsigned int quadrature_point) const
const FiniteElement< dim, spacedim > & get_fe() const
const Tensor< 1, spacedim > & shape_grad(const unsigned int function_no, const unsigned int quadrature_point) const
const double & shape_value(const unsigned int function_no, const unsigned int point_no) const
unsigned int n_components() const
bool is_primitive() const
#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)
void upwind_face_residual(Vector< double > &result1, Vector< double > &result2, const FEValuesBase< dim > &fe1, const FEValuesBase< dim > &fe2, const std::vector< double > &input1, const std::vector< double > &input2, const ArrayView< const std::vector< double > > &velocity, const double factor=1.)
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, const ArrayView< const std::vector< double > > &velocity, const double factor=1.)
void cell_residual(Vector< double > &result, const FEValuesBase< dim > &fe, const std::vector< Tensor< 1, dim > > &input, const ArrayView< const std::vector< double > > &velocity, double factor=1.)
void upwind_value_residual(Vector< double > &result, const FEValuesBase< dim > &fe, const std::vector< double > &input, const std::vector< double > &data, const ArrayView< const std::vector< double > > &velocity, double factor=1.)
void upwind_value_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, const ArrayView< const std::vector< double > > &velocity, double factor=1.)
Library of integrals over cells and faces.
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)