void SIPGLaplace<dim>::assemble_system()
{
const auto cell_worker =
- [&](const auto &cell, auto &scratch_data, auto ©_data) {
+ [&](const typename DoFHandler<dim>::active_cell_iterator &cell,
+ ScratchData &scratch_data,
+ CopyData ©_data) {
const FEValues<dim> &fe_v = scratch_data.reinit(cell);
const unsigned int dofs_per_cell = fe_v.dofs_per_cell;
copy_data.reinit(cell, dofs_per_cell);
- const auto &q_points = scratch_data.get_quadrature_points();
- const unsigned int n_q_points = q_points.size();
- const std::vector<double> &JxW = scratch_data.get_JxW_values();
+ const std::vector<Point<dim>> &q_points =
+ scratch_data.get_quadrature_points();
+ const unsigned int n_q_points = q_points.size();
+ const std::vector<double> &JxW = scratch_data.get_JxW_values();
std::vector<double> rhs(n_q_points);
rhs_function->value_list(q_points, rhs);
};
// Next, we need a function that assembles face integrals on the boundary:
- const auto boundary_worker = [&](const auto &cell,
- const unsigned int &face_no,
- auto &scratch_data,
- auto ©_data) {
- const FEFaceValuesBase<dim> &fe_fv = scratch_data.reinit(cell, face_no);
+ const auto boundary_worker =
+ [&](const typename DoFHandler<dim>::active_cell_iterator &cell,
+ const unsigned int &face_no,
+ ScratchData &scratch_data,
+ CopyData ©_data) {
+ const FEFaceValuesBase<dim> &fe_fv = scratch_data.reinit(cell, face_no);
- const auto &q_points = scratch_data.get_quadrature_points();
- const unsigned int n_q_points = q_points.size();
- const unsigned int dofs_per_cell = fe_fv.dofs_per_cell;
+ const std::vector<Point<dim>> &q_points =
+ scratch_data.get_quadrature_points();
+ const unsigned int n_q_points = q_points.size();
+ const unsigned int dofs_per_cell = fe_fv.dofs_per_cell;
- const std::vector<double> &JxW = scratch_data.get_JxW_values();
- const std::vector<Tensor<1, dim>> &normals =
- scratch_data.get_normal_vectors();
+ const std::vector<double> &JxW = scratch_data.get_JxW_values();
+ const std::vector<Tensor<1, dim>> &normals =
+ scratch_data.get_normal_vectors();
- std::vector<double> g(n_q_points);
- exact_solution->value_list(q_points, g);
+ std::vector<double> g(n_q_points);
+ exact_solution->value_list(q_points, g);
- const double extent1 = cell->measure() / cell->face(face_no)->measure();
- const double penalty = get_penalty_factor(degree, extent1, extent1);
+ const double extent1 = cell->measure() / cell->face(face_no)->measure();
+ const double penalty = get_penalty_factor(degree, extent1, extent1);
- for (unsigned int point = 0; point < n_q_points; ++point)
- {
- for (unsigned int i = 0; i < dofs_per_cell; ++i)
- for (unsigned int j = 0; j < dofs_per_cell; ++j)
- copy_data.cell_matrix(i, j) +=
- (-diffusion_coefficient * // - nu
- fe_fv.shape_value(i, point) * // v_h
- (fe_fv.shape_grad(j, point) * // (grad u_h .
- normals[point]) // n)
+ for (unsigned int point = 0; point < n_q_points; ++point)
+ {
+ for (unsigned int i = 0; i < dofs_per_cell; ++i)
+ for (unsigned int j = 0; j < dofs_per_cell; ++j)
+ copy_data.cell_matrix(i, j) +=
+ (-diffusion_coefficient * // - nu
+ fe_fv.shape_value(i, point) * // v_h
+ (fe_fv.shape_grad(j, point) * // (grad u_h .
+ normals[point]) // n)
- - diffusion_coefficient * // - nu
- (fe_fv.shape_grad(i, point) * // (grad v_h .
- normals[point]) * // n)
- fe_fv.shape_value(j, point) // u_h
+ - diffusion_coefficient * // - nu
+ (fe_fv.shape_grad(i, point) * // (grad v_h .
+ normals[point]) * // n)
+ fe_fv.shape_value(j, point) // u_h
- + diffusion_coefficient * penalty * // + nu sigma
- fe_fv.shape_value(i, point) * // v_h
- fe_fv.shape_value(j, point) // u_h
+ + diffusion_coefficient * penalty * // + nu sigma
+ fe_fv.shape_value(i, point) * // v_h
+ fe_fv.shape_value(j, point) // u_h
- ) *
- JxW[point]; // dx
+ ) *
+ JxW[point]; // dx
- for (unsigned int i = 0; i < dofs_per_cell; ++i)
- copy_data.cell_rhs(i) +=
- (-diffusion_coefficient * // - nu
- (fe_fv.shape_grad(i, point) * // (grad v_h .
- normals[point]) * // n)
- g[point] // g
+ for (unsigned int i = 0; i < dofs_per_cell; ++i)
+ copy_data.cell_rhs(i) +=
+ (-diffusion_coefficient * // - nu
+ (fe_fv.shape_grad(i, point) * // (grad v_h .
+ normals[point]) * // n)
+ g[point] // g
- + diffusion_coefficient * penalty * // + nu sigma
- fe_fv.shape_value(i, point) * g[point] // v_h g
+ + diffusion_coefficient * penalty * // + nu sigma
+ fe_fv.shape_value(i, point) * g[point] // v_h g
- ) *
- JxW[point]; // dx
- }
- };
+ ) *
+ JxW[point]; // dx
+ }
+ };
// Finally, a function that assembles face integrals on interior
// faces. To reinitialize FEInterfaceValues, we need to pass
// cells, face and subface indices (for adaptive refinement) to
// the reinit() function of FEInterfaceValues:
- const auto face_worker = [&](const auto &cell,
- const unsigned int &f,
- const unsigned int &sf,
- const auto &ncell,
- const unsigned int &nf,
- const unsigned int &nsf,
- auto &scratch_data,
- auto ©_data) {
- const FEInterfaceValues<dim> &fe_iv =
- scratch_data.reinit(cell, f, sf, ncell, nf, nsf);
-
- copy_data.face_data.emplace_back();
- CopyDataFace ©_data_face = copy_data.face_data.back();
- const unsigned int n_dofs_face = fe_iv.n_current_interface_dofs();
- copy_data_face.joint_dof_indices = fe_iv.get_interface_dof_indices();
- copy_data_face.cell_matrix.reinit(n_dofs_face, n_dofs_face);
-
- const std::vector<double> &JxW = fe_iv.get_JxW_values();
- const std::vector<Tensor<1, dim>> &normals = fe_iv.get_normal_vectors();
-
- const double extent1 = cell->measure() / cell->face(f)->measure();
- const double extent2 = ncell->measure() / ncell->face(nf)->measure();
- const double penalty = get_penalty_factor(degree, extent1, extent2);
-
- for (const unsigned int point : fe_iv.quadrature_point_indices())
- {
- for (const unsigned int i : fe_iv.dof_indices())
- for (const unsigned int j : fe_iv.dof_indices())
- copy_data_face.cell_matrix(i, j) +=
- (-diffusion_coefficient * // - nu
- fe_iv.jump_in_shape_values(i, point) * // [v_h]
- (fe_iv.average_of_shape_gradients(j, //
- point) * // ({grad u_h} .
- normals[point]) // n)
-
- -
- diffusion_coefficient * // - nu
- (fe_iv.average_of_shape_gradients(i, point) * // (grad v_h .
- normals[point]) * // n)
- fe_iv.jump_in_shape_values(j, point) // [u_h]
-
- + diffusion_coefficient * penalty * // + nu sigma
- fe_iv.jump_in_shape_values(i, point) * // [v_h]
- fe_iv.jump_in_shape_values(j, point) // [u_h]
-
- ) *
- JxW[point]; // dx
- }
- };
+ const auto face_worker =
+ [&](const typename DoFHandler<dim>::cell_iterator &cell,
+ const unsigned int &f,
+ const unsigned int &sf,
+ const typename DoFHandler<dim>::cell_iterator &ncell,
+ const unsigned int &nf,
+ const unsigned int &nsf,
+ ScratchData &scratch_data,
+ CopyData ©_data) {
+ const FEInterfaceValues<dim> &fe_iv =
+ scratch_data.reinit(cell, f, sf, ncell, nf, nsf);
+
+ copy_data.face_data.emplace_back();
+ CopyDataFace ©_data_face = copy_data.face_data.back();
+ const unsigned int n_dofs_face = fe_iv.n_current_interface_dofs();
+ copy_data_face.joint_dof_indices = fe_iv.get_interface_dof_indices();
+ copy_data_face.cell_matrix.reinit(n_dofs_face, n_dofs_face);
+
+ const std::vector<double> &JxW = fe_iv.get_JxW_values();
+ const std::vector<Tensor<1, dim>> &normals = fe_iv.get_normal_vectors();
+
+ const double extent1 = cell->measure() / cell->face(f)->measure();
+ const double extent2 = ncell->measure() / ncell->face(nf)->measure();
+ const double penalty = get_penalty_factor(degree, extent1, extent2);
+
+ for (const unsigned int point : fe_iv.quadrature_point_indices())
+ {
+ for (const unsigned int i : fe_iv.dof_indices())
+ for (const unsigned int j : fe_iv.dof_indices())
+ copy_data_face.cell_matrix(i, j) +=
+ (-diffusion_coefficient * // - nu
+ fe_iv.jump_in_shape_values(i, point) * // [v_h]
+ (fe_iv.average_of_shape_gradients(j, //
+ point) * // ({grad u_h} .
+ normals[point]) // n)
+
+ - diffusion_coefficient * // - nu
+ (fe_iv.average_of_shape_gradients(i,
+ point) * // (grad v_h .
+ normals[point]) * // n)
+ fe_iv.jump_in_shape_values(j, point) // [u_h]
+
+ + diffusion_coefficient * penalty * // + nu sigma
+ fe_iv.jump_in_shape_values(i, point) * // [v_h]
+ fe_iv.jump_in_shape_values(j, point) // [u_h]
+
+ ) *
+ JxW[point]; // dx
+ }
+ };
// The following lambda function will then copy data into the
// global matrix and right-hand side. Though there are no hanging
// AffineConstraints::distribute_local_to_global() functionality.
AffineConstraints<double> constraints;
constraints.close();
- const auto copier = [&](const auto &c) {
+ const auto copier = [&](const CopyData &c) {
constraints.distribute_local_to_global(c.cell_matrix,
c.cell_rhs,
c.local_dof_indices,
system_rhs);
// Copy data from interior face assembly to the global matrix.
- for (auto &cdf : c.face_data)
+ for (const CopyDataFace &cdf : c.face_data)
{
constraints.distribute_local_to_global(cdf.cell_matrix,
cdf.joint_dof_indices,
void SIPGLaplace<dim>::compute_error_estimate()
{
const auto cell_worker =
- [&](const auto &cell, auto &scratch_data, auto ©_data) {
+ [&](const typename DoFHandler<dim>::active_cell_iterator &cell,
+ ScratchData &scratch_data,
+ CopyData ©_data) {
const FEValues<dim> &fe_v = scratch_data.reinit(cell);
copy_data.cell_index = cell->active_cell_index();
- const auto &q_points = fe_v.get_quadrature_points();
- const unsigned int n_q_points = q_points.size();
- const std::vector<double> &JxW = fe_v.get_JxW_values();
+ const std::vector<Point<dim>> &q_points = fe_v.get_quadrature_points();
+ const unsigned int n_q_points = q_points.size();
+ const std::vector<double> &JxW = fe_v.get_JxW_values();
std::vector<Tensor<2, dim>> hessians(n_q_points);
fe_v.get_function_hessians(solution, hessians);
// Next compute boundary terms $\sum_{f\in \partial K \cap \partial \Omega}
// \sigma \left\| [ u_h-g_D ] \right\|_f^2 $:
- const auto boundary_worker = [&](const auto &cell,
- const unsigned int &face_no,
- auto &scratch_data,
- auto ©_data) {
- const FEFaceValuesBase<dim> &fe_fv = scratch_data.reinit(cell, face_no);
+ const auto boundary_worker =
+ [&](const typename DoFHandler<dim>::active_cell_iterator &cell,
+ const unsigned int &face_no,
+ ScratchData &scratch_data,
+ CopyData ©_data) {
+ const FEFaceValuesBase<dim> &fe_fv = scratch_data.reinit(cell, face_no);
- const auto &q_points = fe_fv.get_quadrature_points();
- const unsigned n_q_points = q_points.size();
+ const std::vector<Point<dim>> &q_points = fe_fv.get_quadrature_points();
+ const unsigned n_q_points = q_points.size();
- const std::vector<double> &JxW = fe_fv.get_JxW_values();
+ const std::vector<double> &JxW = fe_fv.get_JxW_values();
- std::vector<double> g(n_q_points);
- exact_solution->value_list(q_points, g);
+ std::vector<double> g(n_q_points);
+ exact_solution->value_list(q_points, g);
- std::vector<double> sol_u(n_q_points);
- fe_fv.get_function_values(solution, sol_u);
+ std::vector<double> sol_u(n_q_points);
+ fe_fv.get_function_values(solution, sol_u);
- const double extent1 = cell->measure() / cell->face(face_no)->measure();
- const double penalty = get_penalty_factor(degree, extent1, extent1);
+ const double extent1 = cell->measure() / cell->face(face_no)->measure();
+ const double penalty = get_penalty_factor(degree, extent1, extent1);
- double difference_norm_square = 0.;
- for (unsigned int point = 0; point < q_points.size(); ++point)
- {
- const double diff = (g[point] - sol_u[point]);
- difference_norm_square += diff * diff * JxW[point];
- }
- copy_data.value += penalty * difference_norm_square;
- };
+ double difference_norm_square = 0.;
+ for (unsigned int point = 0; point < q_points.size(); ++point)
+ {
+ const double diff = (g[point] - sol_u[point]);
+ difference_norm_square += diff * diff * JxW[point];
+ }
+ copy_data.value += penalty * difference_norm_square;
+ };
// And finally interior face terms $\sum_{f\in \partial K}\lbrace \sigma
// \left\| [u_h] \right\|_f^2 + h_f \left\| [\nu \nabla u_h \cdot
// \mathbf n ] \right\|_f^2 \rbrace$:
- const auto face_worker = [&](const auto &cell,
- const unsigned int &f,
- const unsigned int &sf,
- const auto &ncell,
- const unsigned int &nf,
- const unsigned int &nsf,
- auto &scratch_data,
- auto ©_data) {
- const FEInterfaceValues<dim> &fe_iv =
- scratch_data.reinit(cell, f, sf, ncell, nf, nsf);
+ const auto face_worker =
+ [&](const typename DoFHandler<dim>::cell_iterator &cell,
+ const unsigned int &f,
+ const unsigned int &sf,
+ const typename DoFHandler<dim>::cell_iterator &ncell,
+ const unsigned int &nf,
+ const unsigned int &nsf,
+ ScratchData &scratch_data,
+ CopyData ©_data) {
+ const FEInterfaceValues<dim> &fe_iv =
+ scratch_data.reinit(cell, f, sf, ncell, nf, nsf);
- copy_data.face_data.emplace_back();
- CopyDataFace ©_data_face = copy_data.face_data.back();
+ copy_data.face_data.emplace_back();
+ CopyDataFace ©_data_face = copy_data.face_data.back();
- copy_data_face.cell_indices[0] = cell->active_cell_index();
- copy_data_face.cell_indices[1] = ncell->active_cell_index();
+ copy_data_face.cell_indices[0] = cell->active_cell_index();
+ copy_data_face.cell_indices[1] = ncell->active_cell_index();
- const std::vector<double> &JxW = fe_iv.get_JxW_values();
- const std::vector<Tensor<1, dim>> &normals = fe_iv.get_normal_vectors();
+ const std::vector<double> &JxW = fe_iv.get_JxW_values();
+ const std::vector<Tensor<1, dim>> &normals = fe_iv.get_normal_vectors();
- const auto &q_points = fe_iv.get_quadrature_points();
- const unsigned int n_q_points = q_points.size();
+ const std::vector<Point<dim>> &q_points = fe_iv.get_quadrature_points();
+ const unsigned int n_q_points = q_points.size();
- std::vector<double> jump(n_q_points);
- fe_iv.get_jump_in_function_values(solution, jump);
+ std::vector<double> jump(n_q_points);
+ fe_iv.get_jump_in_function_values(solution, jump);
- std::vector<Tensor<1, dim>> grad_jump(n_q_points);
- fe_iv.get_jump_in_function_gradients(solution, grad_jump);
+ std::vector<Tensor<1, dim>> grad_jump(n_q_points);
+ fe_iv.get_jump_in_function_gradients(solution, grad_jump);
- const double h = cell->face(f)->diameter();
+ const double h = cell->face(f)->diameter();
- const double extent1 = cell->measure() / cell->face(f)->measure();
- const double extent2 = ncell->measure() / ncell->face(nf)->measure();
- const double penalty = get_penalty_factor(degree, extent1, extent2);
+ const double extent1 = cell->measure() / cell->face(f)->measure();
+ const double extent2 = ncell->measure() / ncell->face(nf)->measure();
+ const double penalty = get_penalty_factor(degree, extent1, extent2);
- double flux_jump_square = 0;
- double u_jump_square = 0;
- for (unsigned int point = 0; point < n_q_points; ++point)
- {
- u_jump_square += jump[point] * jump[point] * JxW[point];
- const double flux_jump = grad_jump[point] * normals[point];
- flux_jump_square +=
- diffusion_coefficient * flux_jump * flux_jump * JxW[point];
- }
- copy_data_face.values[0] =
- 0.5 * h * (flux_jump_square + penalty * u_jump_square);
- copy_data_face.values[1] = copy_data_face.values[0];
- };
+ double flux_jump_square = 0;
+ double u_jump_square = 0;
+ for (unsigned int point = 0; point < n_q_points; ++point)
+ {
+ u_jump_square += jump[point] * jump[point] * JxW[point];
+ const double flux_jump = grad_jump[point] * normals[point];
+ flux_jump_square +=
+ diffusion_coefficient * flux_jump * flux_jump * JxW[point];
+ }
+ copy_data_face.values[0] =
+ 0.5 * h * (flux_jump_square + penalty * u_jump_square);
+ copy_data_face.values[1] = copy_data_face.values[0];
+ };
// Having computed local contributions for each cell, we still
// need a way to copy these into the global vector that will hold
// the error estimators for all cells:
- const auto copier = [&](const auto ©_data) {
+ const auto copier = [&](const CopyData ©_data) {
if (copy_data.cell_index != numbers::invalid_unsigned_int)
estimated_error_square_per_cell[copy_data.cell_index] +=
copy_data.value;
- for (auto &cdf : copy_data.face_data)
+ for (const CopyDataFace &cdf : copy_data.face_data)
for (unsigned int j = 0; j < 2; ++j)
estimated_error_square_per_cell[cdf.cell_indices[j]] += cdf.values[j];
};
// Assemble $\sum_{K \in \Gamma_h} \nu\|\nabla (u_h - u) \|_K^2 $.
const auto cell_worker =
- [&](const auto &cell, auto &scratch_data, auto ©_data) {
+ [&](const typename DoFHandler<dim>::active_cell_iterator &cell,
+ ScratchData &scratch_data,
+ CopyData ©_data) {
const FEValues<dim> &fe_v = scratch_data.reinit(cell);
copy_data.cell_index = cell->active_cell_index();
- const auto &q_points = fe_v.get_quadrature_points();
- const unsigned int n_q_points = q_points.size();
- const std::vector<double> &JxW = fe_v.get_JxW_values();
+ const std::vector<Point<dim>> &q_points = fe_v.get_quadrature_points();
+ const unsigned int n_q_points = q_points.size();
+ const std::vector<double> &JxW = fe_v.get_JxW_values();
std::vector<Tensor<1, dim>> grad_u(n_q_points);
fe_v.get_function_gradients(solution, grad_u);
};
// Assemble $\sum_{f \in F_b}\sigma \|u_h-g_D\|_f^2$.
- const auto boundary_worker = [&](const auto &cell,
- const unsigned int &face_no,
- auto &scratch_data,
- auto ©_data) {
- const FEFaceValuesBase<dim> &fe_fv = scratch_data.reinit(cell, face_no);
+ const auto boundary_worker =
+ [&](const typename DoFHandler<dim>::active_cell_iterator &cell,
+ const unsigned int &face_no,
+ ScratchData &scratch_data,
+ CopyData ©_data) {
+ const FEFaceValuesBase<dim> &fe_fv = scratch_data.reinit(cell, face_no);
- const auto &q_points = fe_fv.get_quadrature_points();
- const unsigned n_q_points = q_points.size();
+ const std::vector<Point<dim>> &q_points = fe_fv.get_quadrature_points();
+ const unsigned n_q_points = q_points.size();
- const std::vector<double> &JxW = fe_fv.get_JxW_values();
+ const std::vector<double> &JxW = fe_fv.get_JxW_values();
- std::vector<double> g(n_q_points);
- exact_solution->value_list(q_points, g);
+ std::vector<double> g(n_q_points);
+ exact_solution->value_list(q_points, g);
- std::vector<double> sol_u(n_q_points);
- fe_fv.get_function_values(solution, sol_u);
+ std::vector<double> sol_u(n_q_points);
+ fe_fv.get_function_values(solution, sol_u);
- const double extent1 = cell->measure() / cell->face(face_no)->measure();
- const double penalty = get_penalty_factor(degree, extent1, extent1);
+ const double extent1 = cell->measure() / cell->face(face_no)->measure();
+ const double penalty = get_penalty_factor(degree, extent1, extent1);
- double difference_norm_square = 0.;
- for (unsigned int point = 0; point < q_points.size(); ++point)
- {
- const double diff = (g[point] - sol_u[point]);
- difference_norm_square += diff * diff * JxW[point];
- }
- copy_data.value += penalty * difference_norm_square;
- };
+ double difference_norm_square = 0.;
+ for (unsigned int point = 0; point < q_points.size(); ++point)
+ {
+ const double diff = (g[point] - sol_u[point]);
+ difference_norm_square += diff * diff * JxW[point];
+ }
+ copy_data.value += penalty * difference_norm_square;
+ };
// Assemble $\sum_{f \in F_i} \sigma \| [ u_h ] \|_f^2$.
- const auto face_worker = [&](const auto &cell,
- const unsigned int &f,
- const unsigned int &sf,
- const auto &ncell,
- const unsigned int &nf,
- const unsigned int &nsf,
- auto &scratch_data,
- auto ©_data) {
- const FEInterfaceValues<dim> &fe_iv =
- scratch_data.reinit(cell, f, sf, ncell, nf, nsf);
+ const auto face_worker =
+ [&](const typename DoFHandler<dim>::cell_iterator &cell,
+ const unsigned int &f,
+ const unsigned int &sf,
+ const typename DoFHandler<dim>::cell_iterator &ncell,
+ const unsigned int &nf,
+ const unsigned int &nsf,
+ ScratchData &scratch_data,
+ CopyData ©_data) {
+ const FEInterfaceValues<dim> &fe_iv =
+ scratch_data.reinit(cell, f, sf, ncell, nf, nsf);
- copy_data.face_data.emplace_back();
- CopyDataFace ©_data_face = copy_data.face_data.back();
+ copy_data.face_data.emplace_back();
+ CopyDataFace ©_data_face = copy_data.face_data.back();
- copy_data_face.cell_indices[0] = cell->active_cell_index();
- copy_data_face.cell_indices[1] = ncell->active_cell_index();
+ copy_data_face.cell_indices[0] = cell->active_cell_index();
+ copy_data_face.cell_indices[1] = ncell->active_cell_index();
- const std::vector<double> &JxW = fe_iv.get_JxW_values();
+ const std::vector<double> &JxW = fe_iv.get_JxW_values();
- const auto &q_points = fe_iv.get_quadrature_points();
- const unsigned int n_q_points = q_points.size();
+ const std::vector<Point<dim>> &q_points = fe_iv.get_quadrature_points();
+ const unsigned int n_q_points = q_points.size();
- std::vector<double> jump(n_q_points);
- fe_iv.get_jump_in_function_values(solution, jump);
+ std::vector<double> jump(n_q_points);
+ fe_iv.get_jump_in_function_values(solution, jump);
- const double extent1 = cell->measure() / cell->face(f)->measure();
- const double extent2 = ncell->measure() / ncell->face(nf)->measure();
- const double penalty = get_penalty_factor(degree, extent1, extent2);
+ const double extent1 = cell->measure() / cell->face(f)->measure();
+ const double extent2 = ncell->measure() / ncell->face(nf)->measure();
+ const double penalty = get_penalty_factor(degree, extent1, extent2);
- double u_jump_square = 0;
- for (unsigned int point = 0; point < n_q_points; ++point)
- {
- u_jump_square += jump[point] * jump[point] * JxW[point];
- }
- copy_data_face.values[0] = 0.5 * penalty * u_jump_square;
- copy_data_face.values[1] = copy_data_face.values[0];
- };
+ double u_jump_square = 0;
+ for (unsigned int point = 0; point < n_q_points; ++point)
+ {
+ u_jump_square += jump[point] * jump[point] * JxW[point];
+ }
+ copy_data_face.values[0] = 0.5 * penalty * u_jump_square;
+ copy_data_face.values[1] = copy_data_face.values[0];
+ };
- const auto copier = [&](const auto ©_data) {
+ const auto copier = [&](const CopyData ©_data) {
if (copy_data.cell_index != numbers::invalid_unsigned_int)
energy_norm_square_per_cell[copy_data.cell_index] += copy_data.value;
- for (auto &cdf : copy_data.face_data)
+ for (const CopyDataFace &cdf : copy_data.face_data)
for (unsigned int j = 0; j < 2; ++j)
energy_norm_square_per_cell[cdf.cell_indices[j]] += cdf.values[j];
};
projects / dealii.git / commitdiff