83 typedef CellIntegrator<dim, dim, Number> IntegratorCell;
84 typedef FaceIntegrator<dim, dim, Number> IntegratorFace;
86 typedef dealii::LinearAlgebra::distributed::Vector<Number> VectorType;
88 typedef dealii::VectorizedArray<Number> scalar;
89 typedef dealii::Tensor<1, dim, dealii::VectorizedArray<Number>> vector;
93 : matrix_free(
nullptr), dof_index(0), quad_index(0), array_penalty_parameter(0)
98 reinit(dealii::MatrixFree<dim, Number>
const & matrix_free,
99 unsigned int const dof_index,
100 unsigned int const quad_index,
103 this->matrix_free = &matrix_free;
105 this->dof_index = dof_index;
106 this->quad_index = quad_index;
110 unsigned int n_cells = matrix_free.n_cell_batches() + matrix_free.n_ghost_cell_batches();
111 array_penalty_parameter.resize(n_cells);
121 get_integrator_flags()
const
127 flags.face_evaluate = dealii::EvaluationFlags::values;
128 flags.face_integrate = dealii::EvaluationFlags::values;
140 flags.inner_faces = dealii::update_JxW_values | dealii::update_normal_vectors;
141 flags.boundary_faces =
142 dealii::update_JxW_values | dealii::update_normal_vectors | dealii::update_quadrature_points;
148 calculate_penalty_parameter(VectorType
const & velocity)
150 velocity.update_ghost_values();
152 IntegratorCell integrator(*matrix_free, dof_index, quad_index);
154 dealii::AlignedVector<scalar> JxW_values(integrator.n_q_points);
156 ElementType
const element_type =
157 get_element_type(matrix_free->get_dof_handler(dof_index).get_triangulation());
159 unsigned int n_cells = matrix_free->n_cell_batches() + matrix_free->n_ghost_cell_batches();
160 for(
unsigned int cell = 0; cell < n_cells; ++cell)
162 integrator.reinit(cell);
163 integrator.read_dof_values(velocity);
164 integrator.evaluate(dealii::EvaluationFlags::values);
165 scalar volume = dealii::make_vectorized_array<Number>(0.0);
166 scalar norm_U_mean = dealii::make_vectorized_array<Number>(0.0);
168 for(
unsigned int q = 0; q < integrator.n_q_points; ++q)
170 volume += integrator.JxW(q);
171 norm_U_mean += integrator.JxW(q) * integrator.get_value(q).norm();
174 norm_U_mean /= volume;
176 scalar tau_convective = norm_U_mean;
179 scalar tau_viscous = dealii::make_vectorized_array<Number>(data.viscosity) / h_eff;
181 if(data.type_penalty_parameter == TypePenaltyParameter::ConvectiveTerm)
183 array_penalty_parameter[cell] = data.penalty_factor * tau_convective;
185 else if(data.type_penalty_parameter == TypePenaltyParameter::ViscousTerm)
187 array_penalty_parameter[cell] = data.penalty_factor * tau_viscous;
189 else if(data.type_penalty_parameter == TypePenaltyParameter::ViscousAndConvectiveTerms)
191 array_penalty_parameter[cell] = data.penalty_factor * (tau_convective + tau_viscous);
197 reinit_face(IntegratorFace & integrator_m, IntegratorFace & integrator_p)
const
199 tau = 0.5 * (integrator_m.read_cell_data(array_penalty_parameter) +
200 integrator_p.read_cell_data(array_penalty_parameter));
204 reinit_boundary_face(IntegratorFace & integrator_m)
const
206 tau = integrator_m.read_cell_data(array_penalty_parameter);
210 reinit_face_cell_based(dealii::types::boundary_id
const boundary_id,
211 IntegratorFace & integrator_m,
212 IntegratorFace & integrator_p)
const
214 if(boundary_id == dealii::numbers::internal_face_boundary_id)
216 tau = 0.5 * (integrator_m.read_cell_data(array_penalty_parameter) +
217 integrator_p.read_cell_data(array_penalty_parameter));
221 tau = integrator_m.read_cell_data(array_penalty_parameter);
225 inline DEAL_II_ALWAYS_INLINE
227 calculate_flux(vector
const & u_m, vector
const & u_p, vector
const & normal_m)
const
229 vector jump_value = u_m - u_p;
233 if(data.which_components == ContinuityPenaltyComponents::All)
236 flux = tau * jump_value;
238 else if(data.which_components == ContinuityPenaltyComponents::Normal)
240 flux = tau * (jump_value * normal_m) * normal_m;
244 AssertThrow(
false, dealii::ExcMessage(
"not implemented."));
252 dealii::MatrixFree<dim, Number>
const * matrix_free;
254 unsigned int dof_index;
255 unsigned int quad_index;
259 dealii::AlignedVector<scalar> array_penalty_parameter;
288 typedef dealii::LinearAlgebra::distributed::Vector<Number> VectorType;
290 typedef dealii::Tensor<1, dim, dealii::VectorizedArray<Number>> vector;
292 typedef std::pair<unsigned int, unsigned int> Range;
294 typedef FaceIntegrator<dim, dim, Number> IntegratorFace;
302 initialize(dealii::MatrixFree<dim, Number>
const & matrix_free,
304 std::shared_ptr<Kernel>
const kernel);
307 update(VectorType
const & velocity);
311 apply(VectorType & dst, VectorType
const & src)
const;
314 apply_add(VectorType & dst, VectorType
const & src)
const;
318 rhs(VectorType & dst, Number
const evaluation_time)
const;
321 rhs_add(VectorType & dst, Number
const evaluation_time)
const;
325 evaluate(VectorType & dst, VectorType
const & src, Number
const evaluation_time)
const;
328 evaluate_add(VectorType & dst, VectorType
const & src, Number
const evaluation_time)
const;
332 cell_loop_empty(dealii::MatrixFree<dim, Number>
const & matrix_free,
334 VectorType
const & src,
335 Range
const & range)
const;
338 face_loop(dealii::MatrixFree<dim, Number>
const & matrix_free,
340 VectorType
const & src,
341 Range
const & face_range)
const;
344 face_loop_empty(dealii::MatrixFree<dim, Number>
const & matrix_free,
346 VectorType
const & src,
347 Range
const & face_range)
const;
350 boundary_face_loop_hom(dealii::MatrixFree<dim, Number>
const & matrix_free,
352 VectorType
const & src,
353 Range
const & face_range)
const;
356 boundary_face_loop_full(dealii::MatrixFree<dim, Number>
const & matrix_free,
358 VectorType
const & src,
359 Range
const & face_range)
const;
362 boundary_face_loop_inhom(dealii::MatrixFree<dim, Number>
const & matrix_free,
364 VectorType
const & src,
365 Range
const & face_range)
const;
368 do_face_integral(IntegratorFace & integrator_m, IntegratorFace & integrator_p)
const;
371 do_boundary_integral(IntegratorFace & integrator_m,
372 OperatorType
const & operator_type,
373 dealii::types::boundary_id
const & boundary_id)
const;
375 dealii::MatrixFree<dim, Number>
const * matrix_free;
381 std::shared_ptr<Operators::ContinuityPenaltyKernel<dim, Number>> kernel;