51 typedef dealii::VectorizedArray<Number> scalar;
52 typedef dealii::Tensor<1, dim, dealii::VectorizedArray<Number>> vector;
53 typedef dealii::Tensor<2, dim, dealii::VectorizedArray<Number>> tensor;
55 typedef CellIntegrator<dim, dim, Number> IntegratorCell;
56 typedef FaceIntegrator<dim, dim, Number> IntegratorFace;
59 ViscousKernel() : quad_index(0), degree(1), tau(dealii::make_vectorized_array<Number>(0.0))
64 reinit(dealii::MatrixFree<dim, Number>
const & matrix_free,
66 unsigned int const dof_index,
67 unsigned int const quad_index)
71 this->quad_index = quad_index;
73 dealii::FiniteElement<dim>
const & fe = matrix_free.get_dof_handler(dof_index).get_fe();
74 this->degree = fe.degree;
76 calculate_penalty_parameter(matrix_free, dof_index);
78 AssertThrow(data.viscosity >= 0.0, dealii::ExcMessage(
"Viscosity is not set!"));
80 if(data.viscosity_is_variable)
83 viscosity_coefficients.
initialize(matrix_free, quad_index,
true,
false);
89 calculate_penalty_parameter(dealii::MatrixFree<dim, Number>
const & matrix_free,
90 unsigned int const dof_index)
92 IP::calculate_penalty_parameter<dim, Number>(array_penalty_parameter, matrix_free, dof_index);
102 get_quad_index()
const
104 return this->quad_index;
114 set_constant_coefficient(Number
const & constant_coefficient)
120 set_coefficient_cell(
unsigned int const cell,
unsigned int const q, scalar
const & value)
126 get_coefficient_face(
unsigned int const face,
unsigned int const q)
132 set_coefficient_face(
unsigned int const face,
unsigned int const q, scalar
const & value)
138 get_coefficient_face_neighbor(
unsigned int const face,
unsigned int const q)
144 set_coefficient_face_neighbor(
unsigned int const face,
unsigned int const q, scalar
const & value)
150 get_integrator_flags()
const
154 flags.cell_evaluate = dealii::EvaluationFlags::gradients;
155 flags.cell_integrate = dealii::EvaluationFlags::gradients;
157 flags.face_evaluate = dealii::EvaluationFlags::values | dealii::EvaluationFlags::gradients;
158 flags.face_integrate = dealii::EvaluationFlags::values | dealii::EvaluationFlags::gradients;
164 get_mapping_flags(
bool const compute_interior_face_integrals,
165 bool const compute_boundary_face_integrals)
169 flags.cells = dealii::update_JxW_values | dealii::update_gradients;
170 if(compute_interior_face_integrals)
172 dealii::update_JxW_values | dealii::update_gradients | dealii::update_normal_vectors;
173 if(compute_boundary_face_integrals)
174 flags.boundary_faces = dealii::update_JxW_values | dealii::update_gradients |
175 dealii::update_normal_vectors | dealii::update_quadrature_points;
181 reinit_face(IntegratorFace & integrator_m,
182 IntegratorFace & integrator_p,
183 unsigned int const dof_index)
const
185 tau = std::max(integrator_m.read_cell_data(array_penalty_parameter),
186 integrator_p.read_cell_data(array_penalty_parameter)) *
187 IP::get_penalty_factor<dim, Number>(
190 integrator_m.get_matrix_free().get_dof_handler(dof_index).get_triangulation()),
195 reinit_boundary_face(IntegratorFace & integrator_m,
unsigned int const dof_index)
const
197 tau = integrator_m.read_cell_data(array_penalty_parameter) *
198 IP::get_penalty_factor<dim, Number>(
201 integrator_m.get_matrix_free().get_dof_handler(dof_index).get_triangulation()),
206 reinit_face_cell_based(dealii::types::boundary_id
const boundary_id,
207 IntegratorFace & integrator_m,
208 IntegratorFace & integrator_p,
209 unsigned int const dof_index)
const
211 if(boundary_id == dealii::numbers::internal_face_boundary_id)
213 tau = std::max(integrator_m.read_cell_data(array_penalty_parameter),
214 integrator_p.read_cell_data(array_penalty_parameter)) *
215 IP::get_penalty_factor<dim, Number>(
218 integrator_m.get_matrix_free().get_dof_handler(dof_index).get_triangulation()),
223 tau = integrator_m.read_cell_data(array_penalty_parameter) *
224 IP::get_penalty_factor<dim, Number>(
227 integrator_m.get_matrix_free().get_dof_handler(dof_index).get_triangulation()),
235 inline DEAL_II_ALWAYS_INLINE
237 get_viscosity_cell(
unsigned int const cell,
unsigned int const q)
const
239 scalar viscosity = dealii::make_vectorized_array<Number>(data.viscosity);
241 if(data.viscosity_is_variable)
252 inline DEAL_II_ALWAYS_INLINE
254 calculate_average_viscosity(
unsigned int const face,
unsigned int const q)
const
256 scalar average_viscosity = dealii::make_vectorized_array<Number>(0.0);
259 scalar coefficient_face_neighbor =
263 average_viscosity = 2.0 * coefficient_face * coefficient_face_neighbor /
264 (coefficient_face + coefficient_face_neighbor);
272 return average_viscosity;
278 inline DEAL_II_ALWAYS_INLINE
280 get_viscosity_interior_face(
unsigned int const face,
unsigned int const q)
const
282 scalar viscosity = dealii::make_vectorized_array<Number>(data.viscosity);
284 if(data.viscosity_is_variable)
286 viscosity = calculate_average_viscosity(face, q);
295 inline DEAL_II_ALWAYS_INLINE
297 get_viscosity_boundary_face(
unsigned int const face,
unsigned int const q)
const
299 scalar viscosity = dealii::make_vectorized_array<Number>(data.viscosity);
301 if(data.viscosity_is_variable)
312 inline DEAL_II_ALWAYS_INLINE
314 get_volume_flux(tensor
const & gradient, scalar
const & viscosity)
const
316 if(data.formulation_viscous_term == FormulationViscousTerm::DivergenceFormulation)
318 return viscosity * (gradient + transpose(gradient));
320 else if(data.formulation_viscous_term == FormulationViscousTerm::LaplaceFormulation)
322 return viscosity * gradient;
326 AssertThrow(data.formulation_viscous_term == FormulationViscousTerm::DivergenceFormulation or
327 data.formulation_viscous_term == FormulationViscousTerm::LaplaceFormulation,
328 dealii::ExcMessage(
"Specified formulation of viscous term is not implemented."));
338 inline DEAL_II_ALWAYS_INLINE
341 vector
const & value_p,
342 vector
const & normal,
343 scalar
const & viscosity)
const
347 vector jump_value = value_m - value_p;
348 tensor jump_tensor = outer_product(jump_value, normal);
350 if(data.formulation_viscous_term == FormulationViscousTerm::LaplaceFormulation)
352 if(data.IP_formulation == InteriorPenaltyFormulation::NIPG)
354 flux = 0.5 * viscosity * jump_tensor;
356 else if(data.IP_formulation == InteriorPenaltyFormulation::SIPG)
358 flux = -0.5 * viscosity * jump_tensor;
363 false, dealii::ExcMessage(
"Specified interior penalty formulation is not implemented."));
366 else if(data.formulation_viscous_term == FormulationViscousTerm::DivergenceFormulation)
368 if(data.IP_formulation == InteriorPenaltyFormulation::NIPG)
370 flux = 0.5 * viscosity * (jump_tensor + transpose(jump_tensor));
372 else if(data.IP_formulation == InteriorPenaltyFormulation::SIPG)
374 flux = -0.5 * viscosity * (jump_tensor + transpose(jump_tensor));
379 false, dealii::ExcMessage(
"Specified interior penalty formulation is not implemented."));
385 dealii::ExcMessage(
"Specified formulation of viscous term is not implemented."));
395 template<
typename Integrator>
396 inline DEAL_II_ALWAYS_INLINE
398 calculate_normal_gradient(
unsigned int const q, Integrator & integrator)
const
402 if(data.formulation_viscous_term == FormulationViscousTerm::DivergenceFormulation)
408 gradient = dealii::make_vectorized_array<Number>(2.0) * integrator.get_symmetric_gradient(q);
410 else if(data.formulation_viscous_term == FormulationViscousTerm::LaplaceFormulation)
415 gradient = integrator.get_gradient(q);
419 AssertThrow(data.formulation_viscous_term == FormulationViscousTerm::DivergenceFormulation or
420 data.formulation_viscous_term == FormulationViscousTerm::LaplaceFormulation,
421 dealii::ExcMessage(
"Specified formulation of viscous term is not implemented."));
424 vector normal_gradient = gradient * integrator.get_normal_vector(q);
426 return normal_gradient;
436 inline DEAL_II_ALWAYS_INLINE
439 vector
const & normal_gradient_p,
440 vector
const & value_m,
441 vector
const & value_p,
442 vector
const & normal,
443 scalar
const & viscosity)
const
447 vector jump_value = value_m - value_p;
448 vector average_normal_gradient = 0.5 * (normal_gradient_m + normal_gradient_p);
450 if(data.formulation_viscous_term == FormulationViscousTerm::DivergenceFormulation)
452 if(data.penalty_term_div_formulation == PenaltyTermDivergenceFormulation::Symmetrized)
454 flux = viscosity * average_normal_gradient -
455 viscosity * tau * (jump_value + (jump_value * normal) * normal);
457 else if(data.penalty_term_div_formulation == PenaltyTermDivergenceFormulation::NotSymmetrized)
459 flux = viscosity * average_normal_gradient - viscosity * tau * jump_value;
464 data.penalty_term_div_formulation == PenaltyTermDivergenceFormulation::Symmetrized or
465 data.penalty_term_div_formulation == PenaltyTermDivergenceFormulation::NotSymmetrized,
466 dealii::ExcMessage(
"Specified formulation of viscous term is not implemented."));
469 else if(data.formulation_viscous_term == FormulationViscousTerm::LaplaceFormulation)
471 flux = viscosity * average_normal_gradient - viscosity * tau * jump_value;
475 AssertThrow(data.formulation_viscous_term == FormulationViscousTerm::DivergenceFormulation or
476 data.formulation_viscous_term == FormulationViscousTerm::LaplaceFormulation,
477 dealii::ExcMessage(
"Specified formulation of viscous term is not implemented."));
495 template<
typename Integrator>
496 inline DEAL_II_ALWAYS_INLINE
498 calculate_interior_normal_gradient(
unsigned int const q,
499 Integrator
const & integrator,
500 OperatorType
const & operator_type)
const
502 vector normal_gradient_m;
504 if(operator_type == OperatorType::full or operator_type == OperatorType::homogeneous)
506 normal_gradient_m = calculate_normal_gradient(q, integrator);
508 else if(operator_type == OperatorType::inhomogeneous)
514 AssertThrow(
false, dealii::ExcMessage(
"Specified OperatorType is not implemented!"));
517 return normal_gradient_m;
521 ViscousKernelData data;
523 unsigned int quad_index;
527 dealii::AlignedVector<scalar> array_penalty_parameter;
531 VariableCoefficients<dealii::VectorizedArray<Number>> viscosity_coefficients;
552 typedef dealii::VectorizedArray<Number> scalar;
553 typedef dealii::Tensor<1, dim, dealii::VectorizedArray<Number>> vector;
554 typedef dealii::Tensor<2, dim, dealii::VectorizedArray<Number>> tensor;
558 typedef typename Base::VectorType VectorType;
559 typedef typename Base::Range Range;
560 typedef typename Base::IntegratorCell IntegratorCell;
561 typedef typename Base::IntegratorFace IntegratorFace;
564 initialize(dealii::MatrixFree<dim, Number>
const & matrix_free,
565 dealii::AffineConstraints<Number>
const & affine_constraints,
574 reinit_face_derived(IntegratorFace & integrator_m,
575 IntegratorFace & integrator_p,
576 unsigned int const face)
const final;
579 reinit_boundary_face_derived(IntegratorFace & integrator_m,
unsigned int const face)
const final;
582 reinit_face_cell_based_derived(IntegratorFace & integrator_m,
583 IntegratorFace & integrator_p,
584 unsigned int const cell,
585 unsigned int const face,
586 dealii::types::boundary_id
const boundary_id)
const final;
589 do_cell_integral(IntegratorCell & integrator)
const final;
592 do_face_integral(IntegratorFace & integrator_m, IntegratorFace & integrator_p)
const final;
595 do_face_int_integral(IntegratorFace & integrator_m, IntegratorFace & integrator_p)
const final;
598 do_face_ext_integral(IntegratorFace & integrator_m, IntegratorFace & integrator_p)
const final;
601 do_boundary_integral(IntegratorFace & integrator,
602 OperatorType
const & operator_type,
603 dealii::types::boundary_id
const & boundary_id)
const final;
607 std::shared_ptr<Operators::ViscousKernel<dim, Number>> kernel;