16 #ifndef dealii_integrators_advection_h
17 #define dealii_integrators_advection_h
83 const ArrayView<
const std::vector<double>> &velocity,
84 const double factor = 1.)
95 const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
107 const double dx = factor * fe.
JxW(k);
108 const unsigned int vindex = k * v_increment;
110 for (
unsigned j = 0; j < n_dofs; ++j)
111 for (
unsigned i = 0; i < t_dofs; ++i)
116 for (
unsigned int d = 1;
d < dim; ++
d)
139 const ArrayView<
const std::vector<double>> &velocity,
149 const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
150 if (v_increment == 1)
155 for (
unsigned k = 0; k < nq; ++k)
157 const double dx = factor * fe.
JxW(k);
158 for (
unsigned i = 0; i < n_dofs; ++i)
159 for (
unsigned int d = 0;
d < dim; ++
d)
161 velocity[
d][k * v_increment];
182 const ArrayView<
const std::vector<double>> & velocity,
187 const unsigned int n_comp = fe.
get_fe().n_components();
194 const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
195 if (v_increment == 1)
200 for (
unsigned k = 0; k < nq; ++k)
202 const double dx = factor * fe.
JxW(k);
203 for (
unsigned i = 0; i < n_dofs; ++i)
204 for (
unsigned int c = 0; c < n_comp; ++c)
205 for (
unsigned int d = 0;
d < dim; ++
d)
206 result(i) +=
dx * input[c][k][
d] *
208 velocity[
d][k * v_increment];
223 const std::vector<double> & input,
224 const ArrayView<
const std::vector<double>> &velocity,
234 const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
235 if (v_increment == 1)
240 for (
unsigned k = 0; k < nq; ++k)
242 const double dx = factor * fe.
JxW(k);
243 for (
unsigned i = 0; i < n_dofs; ++i)
244 for (
unsigned int d = 0;
d < dim; ++
d)
246 velocity[
d][k * v_increment];
263 const ArrayView<
const std::vector<double>> &input,
264 const ArrayView<
const std::vector<double>> &velocity,
269 const unsigned int n_comp = fe.
get_fe().n_components();
276 const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
277 if (v_increment == 1)
282 for (
unsigned k = 0; k < nq; ++k)
284 const double dx = factor * fe.
JxW(k);
285 for (
unsigned i = 0; i < n_dofs; ++i)
286 for (
unsigned int c = 0; c < n_comp; ++c)
287 for (
unsigned int d = 0;
d < dim; ++
d)
288 result(i) -=
dx * input[c][k] *
290 velocity[
d][k * v_increment];
318 const ArrayView<
const std::vector<double>> &velocity,
328 const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
329 if (v_increment == 1)
336 const double dx = factor * fe.
JxW(k);
339 for (
unsigned int d = 0;
d < dim; ++
d)
344 for (
unsigned i = 0; i < t_dofs; ++i)
345 for (
unsigned j = 0; j < n_dofs; ++j)
347 if (fe.
get_fe().is_primitive())
390 const std::vector<double> & input,
391 const std::vector<double> & data,
392 const ArrayView<
const std::vector<double>> &velocity,
401 const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
402 if (v_increment == 1)
410 const double dx = factor * fe.
JxW(k);
413 for (
unsigned int d = 0;
d < dim; ++
d)
417 const double val = (nv > 0.) ? input[k] : -data[k];
419 for (
unsigned i = 0; i < n_dofs; ++i)
422 result(i) +=
dx * nv * val * v;
457 const ArrayView<
const std::vector<double>> &input,
458 const ArrayView<
const std::vector<double>> &data,
459 const ArrayView<
const std::vector<double>> &velocity,
463 const unsigned int n_comp = fe.
get_fe().n_components();
469 const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
470 if (v_increment == 1)
478 const double dx = factor * fe.
JxW(k);
481 for (
unsigned int d = 0;
d < dim; ++
d)
484 std::vector<double> val(n_comp);
486 for (
unsigned int d = 0;
d < n_comp; ++
d)
488 val[
d] = (nv > 0.) ? input[
d][k] : -data[
d][k];
489 for (
unsigned i = 0; i < n_dofs; ++i)
492 result(i) +=
dx * nv * val[
d] * v;
530 const ArrayView<
const std::vector<double>> &velocity,
531 const double factor = 1.)
539 const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
540 if (v_increment == 1)
547 double nbeta = fe1.
normal_vector(k)[0] * velocity[0][k * v_increment];
548 for (
unsigned int d = 1;
d < dim; ++
d)
550 const double dx_nbeta = factor * std::abs(nbeta) * fe1.
JxW(k);
556 for (
unsigned i = 0; i < n1; ++i)
557 for (
unsigned j = 0; j < n1; ++j)
559 if (fe1.
get_fe().is_primitive())
571 M1(i, j) += dx_nbeta *
574 M2(i, j) -= dx_nbeta *
611 const std::vector<double> & input1,
612 const std::vector<double> & input2,
613 const ArrayView<
const std::vector<double>> &velocity,
614 const double factor = 1.)
627 const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
628 if (v_increment == 1)
635 double nbeta = fe1.
normal_vector(k)[0] * velocity[0][k * v_increment];
636 for (
unsigned int d = 1;
d < dim; ++
d)
638 const double dx_nbeta = factor * nbeta * fe1.
JxW(k);
640 for (
unsigned i = 0; i < n1; ++i)
644 const double u1 = input1[k];
645 const double u2 = input2[k];
648 result1(i) += dx_nbeta * u1 * v1;
649 result2(i) -= dx_nbeta * u1 * v2;
653 result1(i) += dx_nbeta * u2 * v1;
654 result2(i) -= dx_nbeta * u2 * v2;
688 const ArrayView<
const std::vector<double>> &input1,
689 const ArrayView<
const std::vector<double>> &input2,
690 const ArrayView<
const std::vector<double>> &velocity,
691 const double factor = 1.)
693 const unsigned int n_comp = fe1.
get_fe().n_components();
703 const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
704 if (v_increment == 1)
711 double nbeta = fe1.
normal_vector(k)[0] * velocity[0][k * v_increment];
712 for (
unsigned int d = 1;
d < dim; ++
d)
714 const double dx_nbeta = factor * nbeta * fe1.
JxW(k);
716 for (
unsigned i = 0; i < n1; ++i)
717 for (
unsigned int d = 0;
d < n_comp; ++
d)
721 const double u1 = input1[
d][k];
722 const double u2 = input2[
d][k];
725 result1(i) += dx_nbeta * u1 * v1;
726 result2(i) -= dx_nbeta * u1 * v2;
730 result1(i) += dx_nbeta * u2 * v1;
731 result2(i) -= dx_nbeta * u2 * v2;