72 typedef CellIntegrator<dim, dim, Number> IntegratorCell;
74 typedef dealii::LinearAlgebra::distributed::Vector<Number> VectorType;
76 typedef dealii::VectorizedArray<Number> scalar;
80 : matrix_free(
nullptr), dof_index(0), quad_index(0), array_penalty_parameter(0)
85 reinit(dealii::MatrixFree<dim, Number>
const & matrix_free,
86 unsigned int const dof_index,
87 unsigned int const quad_index,
90 this->matrix_free = &matrix_free;
92 this->dof_index = dof_index;
93 this->quad_index = quad_index;
97 unsigned int n_cells = matrix_free.n_cell_batches() + matrix_free.n_ghost_cell_batches();
98 array_penalty_parameter.resize(n_cells);
108 get_integrator_flags()
const
112 flags.cell_evaluate = dealii::EvaluationFlags::gradients;
113 flags.cell_integrate = dealii::EvaluationFlags::gradients;
125 flags.cells = dealii::update_JxW_values | dealii::update_gradients;
133 calculate_penalty_parameter(VectorType
const & velocity)
135 velocity.update_ghost_values();
137 IntegratorCell integrator(*matrix_free, dof_index, quad_index);
139 dealii::AlignedVector<scalar> JxW_values(integrator.n_q_points);
141 ElementType
const element_type =
142 get_element_type(matrix_free->get_dof_handler(dof_index).get_triangulation());
144 unsigned int n_cells = matrix_free->n_cell_batches() + matrix_free->n_ghost_cell_batches();
145 for(
unsigned int cell = 0; cell < n_cells; ++cell)
147 scalar tau_convective = dealii::make_vectorized_array<Number>(0.0);
148 scalar tau_viscous = dealii::make_vectorized_array<Number>(data.viscosity);
150 if(data.type_penalty_parameter == TypePenaltyParameter::ConvectiveTerm or
151 data.type_penalty_parameter == TypePenaltyParameter::ViscousAndConvectiveTerms)
153 integrator.reinit(cell);
154 integrator.read_dof_values(velocity);
155 integrator.evaluate(dealii::EvaluationFlags::values);
157 scalar volume = dealii::make_vectorized_array<Number>(0.0);
158 scalar norm_U_mean = dealii::make_vectorized_array<Number>(0.0);
159 for(
unsigned int q = 0; q < integrator.n_q_points; ++q)
161 volume += integrator.JxW(q);
162 norm_U_mean += integrator.JxW(q) * integrator.get_value(q).norm();
164 norm_U_mean /= volume;
169 tau_convective = norm_U_mean * h_eff;
172 if(data.type_penalty_parameter == TypePenaltyParameter::ConvectiveTerm)
174 array_penalty_parameter[cell] = data.penalty_factor * tau_convective;
176 else if(data.type_penalty_parameter == TypePenaltyParameter::ViscousTerm)
178 array_penalty_parameter[cell] = data.penalty_factor * tau_viscous;
180 else if(data.type_penalty_parameter == TypePenaltyParameter::ViscousAndConvectiveTerms)
182 array_penalty_parameter[cell] = data.penalty_factor * (tau_convective + tau_viscous);
188 reinit_cell(IntegratorCell & integrator)
const
190 tau = integrator.read_cell_data(array_penalty_parameter);
196 inline DEAL_II_ALWAYS_INLINE
198 get_volume_flux(IntegratorCell
const & integrator,
unsigned int const q)
const
200 return tau * integrator.get_divergence(q);
204 dealii::MatrixFree<dim, Number>
const * matrix_free;
206 unsigned int dof_index;
207 unsigned int quad_index;
211 dealii::AlignedVector<scalar> array_penalty_parameter;
234 typedef dealii::LinearAlgebra::distributed::Vector<Number> VectorType;
236 typedef std::pair<unsigned int, unsigned int> Range;
238 typedef CellIntegrator<dim, dim, Number> IntegratorCell;
246 initialize(dealii::MatrixFree<dim, Number>
const & matrix_free,
248 std::shared_ptr<Kernel>
const kernel);
251 update(VectorType
const & velocity);
254 apply(VectorType & dst, VectorType
const & src)
const;
257 apply_add(VectorType & dst, VectorType
const & src)
const;
261 cell_loop(dealii::MatrixFree<dim, Number>
const & matrix_free,
263 VectorType
const & src,
264 Range
const & range)
const;
267 do_cell_integral(IntegratorCell & integrator)
const;
269 dealii::MatrixFree<dim, Number>
const * matrix_free;
273 std::shared_ptr<Operators::DivergencePenaltyKernel<dim, Number>> kernel;