22#ifndef INCLUDE_EXADG_OPERATORS_INTERIOR_PENALTY_PARAMETER_H_
23#define INCLUDE_EXADG_OPERATORS_INTERIOR_PENALTY_PARAMETER_H_
26#include <deal.II/base/aligned_vector.h>
27#include <deal.II/base/vectorization.h>
28#include <deal.II/fe/mapping.h>
29#include <deal.II/matrix_free/fe_evaluation.h>
30#include <deal.II/matrix_free/matrix_free.h>
33#include <exadg/grid/grid_data.h>
43template<
int dim,
typename Number>
45calculate_penalty_parameter(
46 dealii::AlignedVector<dealii::VectorizedArray<Number>> & array_penalty_parameter,
47 dealii::MatrixFree<dim, Number>
const & matrix_free,
48 unsigned int const dof_index = 0)
50 unsigned int n_cells = matrix_free.n_cell_batches() + matrix_free.n_ghost_cell_batches();
51 array_penalty_parameter.resize(n_cells);
53 dealii::Mapping<dim>
const & mapping = *matrix_free.get_mapping_info().mapping;
54 dealii::FiniteElement<dim>
const & fe = matrix_free.get_dof_handler(dof_index).get_fe();
55 unsigned int const degree = fe.degree;
57 auto const reference_cells =
58 matrix_free.get_dof_handler(dof_index).get_triangulation().get_reference_cells();
59 AssertThrow(reference_cells.size() == 1, dealii::ExcMessage(
"No mixed meshes allowed."));
61 auto const quadrature = reference_cells[0].template get_gauss_type_quadrature<dim>(degree + 1);
62 dealii::FEValues<dim> fe_values(mapping, fe, quadrature, dealii::update_JxW_values);
64 auto const face_quadrature =
65 reference_cells[0].face_reference_cell(0).template get_gauss_type_quadrature<dim - 1>(degree +
67 dealii::FEFaceValues<dim> fe_face_values(mapping, fe, face_quadrature, dealii::update_JxW_values);
69 for(
unsigned int i = 0; i < n_cells; ++i)
71 for(
unsigned int v = 0; v < matrix_free.n_active_entries_per_cell_batch(i); ++v)
73 typename dealii::DoFHandler<dim>::cell_iterator cell =
74 matrix_free.get_cell_iterator(i, v, dof_index);
75 fe_values.reinit(cell);
79 for(
unsigned int q = 0; q < quadrature.size(); ++q)
81 volume += fe_values.JxW(q);
85 Number surface_area = 0;
86 for(
unsigned int const f : cell->face_indices())
88 fe_face_values.reinit(cell, f);
90 (cell->at_boundary(f) and not(cell->has_periodic_neighbor(f))) ? 1. : 0.5;
91 for(
unsigned int q = 0; q < face_quadrature.size(); ++q)
93 surface_area += fe_face_values.JxW(q) * factor;
97 array_penalty_parameter[i][v] = surface_area / volume;
109template<
int dim,
typename Number>
111get_penalty_factor(
unsigned int const degree, ElementType element_type, Number
const factor = 1.0)
113 if(element_type == ElementType::Simplex)
118 return factor * (degree + 1.0) * (degree + dim) / dim;
120 else if(element_type == ElementType::Hypercube)
125 return factor * (degree + 1.0) * (degree + 1.0);
128 AssertThrow(
false, dealii::ExcMessage(
"Only hypercube or simplex elements are supported."));