15#ifndef dealii_integrators_advection_h
16#define dealii_integrators_advection_h
77 const ArrayView<
const std::vector<double>> &velocity,
78 const double factor = 1.)
89 const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
101 const double dx = factor * fe.
JxW(k);
102 const unsigned int vindex = k * v_increment;
104 for (
unsigned j = 0; j < n_dofs; ++j)
105 for (
unsigned i = 0; i < t_dofs; ++i)
106 for (
unsigned int c = 0; c < n_components; ++c)
110 for (
unsigned int d = 1; d < dim; ++d)
133 const ArrayView<
const std::vector<double>> &velocity,
143 const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
144 if (v_increment == 1)
149 for (
unsigned k = 0; k < nq; ++k)
151 const double dx = factor * fe.
JxW(k);
152 for (
unsigned i = 0; i < n_dofs; ++i)
153 for (
unsigned int d = 0; d < dim; ++d)
154 result(i) += dx * input[k][d] * fe.
shape_value(i, k) *
155 velocity[d][k * v_increment];
176 const ArrayView<
const std::vector<double>> &velocity,
188 const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
189 if (v_increment == 1)
194 for (
unsigned k = 0; k < nq; ++k)
196 const double dx = factor * fe.
JxW(k);
197 for (
unsigned i = 0; i < n_dofs; ++i)
198 for (
unsigned int c = 0; c < n_comp; ++c)
199 for (
unsigned int d = 0; d < dim; ++d)
200 result(i) += dx * input[c][k][d] *
202 velocity[d][k * v_increment];
217 const std::vector<double> &input,
218 const ArrayView<
const std::vector<double>> &velocity,
228 const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
229 if (v_increment == 1)
234 for (
unsigned k = 0; k < nq; ++k)
236 const double dx = factor * fe.
JxW(k);
237 for (
unsigned i = 0; i < n_dofs; ++i)
238 for (
unsigned int d = 0; d < dim; ++d)
239 result(i) -= dx * input[k] * fe.
shape_grad(i, k)[d] *
240 velocity[d][k * v_increment];
257 const ArrayView<
const std::vector<double>> &input,
258 const ArrayView<
const std::vector<double>> &velocity,
270 const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
271 if (v_increment == 1)
276 for (
unsigned k = 0; k < nq; ++k)
278 const double dx = factor * fe.
JxW(k);
279 for (
unsigned i = 0; i < n_dofs; ++i)
280 for (
unsigned int c = 0; c < n_comp; ++c)
281 for (
unsigned int d = 0; d < dim; ++d)
282 result(i) -= dx * input[c][k] *
284 velocity[d][k * v_increment];
312 const ArrayView<
const std::vector<double>> &velocity,
322 const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
323 if (v_increment == 1)
330 const double dx = factor * fe.
JxW(k);
333 for (
unsigned int d = 0; d < dim; ++d)
338 for (
unsigned i = 0; i < t_dofs; ++i)
339 for (
unsigned j = 0; j < n_dofs; ++j)
345 for (
unsigned int c = 0; c < n_components; ++c)
384 const std::vector<double> &input,
385 const std::vector<double> &
data,
386 const ArrayView<
const std::vector<double>> &velocity,
395 const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
396 if (v_increment == 1)
404 const double dx = factor * fe.
JxW(k);
407 for (
unsigned int d = 0; d < dim; ++d)
411 const double val = (nv > 0.) ? input[k] : -
data[k];
413 for (
unsigned i = 0; i < n_dofs; ++i)
416 result(i) += dx * nv * val * v;
451 const ArrayView<
const std::vector<double>> &input,
453 const ArrayView<
const std::vector<double>> &velocity,
463 const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
464 if (v_increment == 1)
472 const double dx = factor * fe.
JxW(k);
475 for (
unsigned int d = 0; d < dim; ++d)
478 std::vector<double> val(n_comp);
480 for (
unsigned int d = 0; d < n_comp; ++d)
482 val[d] = (nv > 0.) ? input[d][k] : -
data[d][k];
483 for (
unsigned i = 0; i < n_dofs; ++i)
486 result(i) += dx * nv * val[d] * v;
524 const ArrayView<
const std::vector<double>> &velocity,
525 const double factor = 1.)
533 const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
534 if (v_increment == 1)
541 double nbeta = fe1.
normal_vector(k)[0] * velocity[0][k * v_increment];
542 for (
unsigned int d = 1; d < dim; ++d)
543 nbeta += fe1.
normal_vector(k)[d] * velocity[d][k * v_increment];
544 const double dx_nbeta = factor *
std::abs(nbeta) * fe1.
JxW(k);
550 for (
unsigned i = 0; i < n1; ++i)
551 for (
unsigned j = 0; j < n1; ++j)
562 for (
unsigned int d = 0; d < fe1.
get_fe().n_components();
565 M1(i, j) += dx_nbeta *
568 M2(i, j) -= dx_nbeta *
605 const std::vector<double> &input1,
606 const std::vector<double> &input2,
607 const ArrayView<
const std::vector<double>> &velocity,
608 const double factor = 1.)
621 const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
622 if (v_increment == 1)
629 double nbeta = fe1.
normal_vector(k)[0] * velocity[0][k * v_increment];
630 for (
unsigned int d = 1; d < dim; ++d)
631 nbeta += fe1.
normal_vector(k)[d] * velocity[d][k * v_increment];
632 const double dx_nbeta = factor * nbeta * fe1.
JxW(k);
634 for (
unsigned i = 0; i < n1; ++i)
638 const double u1 = input1[k];
639 const double u2 = input2[k];
642 result1(i) += dx_nbeta * u1 *
v1;
643 result2(i) -= dx_nbeta * u1 * v2;
647 result1(i) += dx_nbeta * u2 *
v1;
648 result2(i) -= dx_nbeta * u2 * v2;
682 const ArrayView<
const std::vector<double>> &input1,
683 const ArrayView<
const std::vector<double>> &input2,
684 const ArrayView<
const std::vector<double>> &velocity,
685 const double factor = 1.)
697 const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
698 if (v_increment == 1)
705 double nbeta = fe1.
normal_vector(k)[0] * velocity[0][k * v_increment];
706 for (
unsigned int d = 1; d < dim; ++d)
707 nbeta += fe1.
normal_vector(k)[d] * velocity[d][k * v_increment];
708 const double dx_nbeta = factor * nbeta * fe1.
JxW(k);
710 for (
unsigned i = 0; i < n1; ++i)
711 for (
unsigned int d = 0; d < n_comp; ++d)
715 const double u1 = input1[d][k];
716 const double u2 = input2[d][k];
719 result1(i) += dx_nbeta * u1 *
v1;
720 result2(i) -= dx_nbeta * u1 * v2;
724 result1(i) += dx_nbeta * u2 *
v1;
725 result2(i) -= dx_nbeta * u2 * v2;
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
bool is_primitive() const
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 & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
std::vector< index_type > data
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.
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)