13#ifndef dealii_integrators_advection_h
14#define dealii_integrators_advection_h
75 const ArrayView<
const std::vector<double>> &velocity,
76 const double factor = 1.)
87 const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
99 const double dx = factor * fe.
JxW(k);
100 const unsigned int vindex = k * v_increment;
102 for (
unsigned j = 0; j < n_dofs; ++j)
103 for (
unsigned i = 0; i < t_dofs; ++i)
104 for (
unsigned int c = 0; c < n_components; ++c)
108 for (
unsigned int d = 1; d < dim; ++d)
131 const ArrayView<
const std::vector<double>> &velocity,
141 const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
142 if (v_increment == 1)
147 for (
unsigned k = 0; k < nq; ++k)
149 const double dx = factor * fe.
JxW(k);
150 for (
unsigned i = 0; i < n_dofs; ++i)
151 for (
unsigned int d = 0; d < dim; ++d)
152 result(i) += dx * input[k][d] * fe.
shape_value(i, k) *
153 velocity[d][k * v_increment];
174 const ArrayView<
const std::vector<double>> &velocity,
186 const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
187 if (v_increment == 1)
192 for (
unsigned k = 0; k < nq; ++k)
194 const double dx = factor * fe.
JxW(k);
195 for (
unsigned i = 0; i < n_dofs; ++i)
196 for (
unsigned int c = 0; c < n_comp; ++c)
197 for (
unsigned int d = 0; d < dim; ++d)
198 result(i) += dx * input[c][k][d] *
200 velocity[d][k * v_increment];
215 const std::vector<double> &input,
216 const ArrayView<
const std::vector<double>> &velocity,
226 const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
227 if (v_increment == 1)
232 for (
unsigned k = 0; k < nq; ++k)
234 const double dx = factor * fe.
JxW(k);
235 for (
unsigned i = 0; i < n_dofs; ++i)
236 for (
unsigned int d = 0; d < dim; ++d)
237 result(i) -= dx * input[k] * fe.
shape_grad(i, k)[d] *
238 velocity[d][k * v_increment];
255 const ArrayView<
const std::vector<double>> &input,
256 const ArrayView<
const std::vector<double>> &velocity,
268 const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
269 if (v_increment == 1)
274 for (
unsigned k = 0; k < nq; ++k)
276 const double dx = factor * fe.
JxW(k);
277 for (
unsigned i = 0; i < n_dofs; ++i)
278 for (
unsigned int c = 0; c < n_comp; ++c)
279 for (
unsigned int d = 0; d < dim; ++d)
280 result(i) -= dx * input[c][k] *
282 velocity[d][k * v_increment];
310 const ArrayView<
const std::vector<double>> &velocity,
320 const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
321 if (v_increment == 1)
328 const double dx = factor * fe.
JxW(k);
331 for (
unsigned int d = 0; d < dim; ++d)
336 for (
unsigned i = 0; i < t_dofs; ++i)
337 for (
unsigned j = 0; j < n_dofs; ++j)
343 for (
unsigned int c = 0; c < n_components; ++c)
382 const std::vector<double> &input,
383 const std::vector<double> &
data,
384 const ArrayView<
const std::vector<double>> &velocity,
393 const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
394 if (v_increment == 1)
402 const double dx = factor * fe.
JxW(k);
405 for (
unsigned int d = 0; d < dim; ++d)
409 const double val = (nv > 0.) ? input[k] : -
data[k];
411 for (
unsigned i = 0; i < n_dofs; ++i)
414 result(i) += dx * nv * val * v;
449 const ArrayView<
const std::vector<double>> &input,
451 const ArrayView<
const std::vector<double>> &velocity,
461 const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
462 if (v_increment == 1)
470 const double dx = factor * fe.
JxW(k);
473 for (
unsigned int d = 0; d < dim; ++d)
476 std::vector<double> val(n_comp);
478 for (
unsigned int d = 0; d < n_comp; ++d)
480 val[d] = (nv > 0.) ? input[d][k] : -
data[d][k];
481 for (
unsigned i = 0; i < n_dofs; ++i)
484 result(i) += dx * nv * val[d] * v;
522 const ArrayView<
const std::vector<double>> &velocity,
523 const double factor = 1.)
531 const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
532 if (v_increment == 1)
539 double nbeta = fe1.
normal_vector(k)[0] * velocity[0][k * v_increment];
540 for (
unsigned int d = 1; d < dim; ++d)
541 nbeta += fe1.
normal_vector(k)[d] * velocity[d][k * v_increment];
542 const double dx_nbeta = factor *
std::abs(nbeta) * fe1.
JxW(k);
548 for (
unsigned i = 0; i < n1; ++i)
549 for (
unsigned j = 0; j < n1; ++j)
560 for (
unsigned int d = 0; d < fe1.
get_fe().n_components();
563 M1(i, j) += dx_nbeta *
566 M2(i, j) -= dx_nbeta *
603 const std::vector<double> &input1,
604 const std::vector<double> &input2,
605 const ArrayView<
const std::vector<double>> &velocity,
606 const double factor = 1.)
619 const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
620 if (v_increment == 1)
627 double nbeta = fe1.
normal_vector(k)[0] * velocity[0][k * v_increment];
628 for (
unsigned int d = 1; d < dim; ++d)
629 nbeta += fe1.
normal_vector(k)[d] * velocity[d][k * v_increment];
630 const double dx_nbeta = factor * nbeta * fe1.
JxW(k);
632 for (
unsigned i = 0; i < n1; ++i)
636 const double u1 = input1[k];
637 const double u2 = input2[k];
640 result1(i) += dx_nbeta * u1 *
v1;
641 result2(i) -= dx_nbeta * u1 * v2;
645 result1(i) += dx_nbeta * u2 *
v1;
646 result2(i) -= dx_nbeta * u2 * v2;
680 const ArrayView<
const std::vector<double>> &input1,
681 const ArrayView<
const std::vector<double>> &input2,
682 const ArrayView<
const std::vector<double>> &velocity,
683 const double factor = 1.)
695 const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
696 if (v_increment == 1)
703 double nbeta = fe1.
normal_vector(k)[0] * velocity[0][k * v_increment];
704 for (
unsigned int d = 1; d < dim; ++d)
705 nbeta += fe1.
normal_vector(k)[d] * velocity[d][k * v_increment];
706 const double dx_nbeta = factor * nbeta * fe1.
JxW(k);
708 for (
unsigned i = 0; i < n1; ++i)
709 for (
unsigned int d = 0; d < n_comp; ++d)
713 const double u1 = input1[d][k];
714 const double u2 = input2[d][k];
717 result1(i) += dx_nbeta * u1 *
v1;
718 result2(i) -= dx_nbeta * u1 * v2;
722 result1(i) += dx_nbeta * u2 *
v1;
723 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 > &)