diff --git a/doc/news/changes/incompatibilities/20210604DanielArndt b/doc/news/changes/incompatibilities/20210604DanielArndt
new file mode 100644
index 000000000000..9e5d6d39e9b8
--- /dev/null
+++ b/doc/news/changes/incompatibilities/20210604DanielArndt
@@ -0,0 +1,4 @@
+Removed: The deprecated specializations of VectorTools::integrate_difference()
+and VectorTools::project_boundary_values_curl_conforming() have been removed.
+
+(Daniel Arndt, 2021/06/04)
diff --git a/include/deal.II/numerics/vector_tools_boundary.h b/include/deal.II/numerics/vector_tools_boundary.h
index 273a1d0a2b5d..056d53c9158b 100644
--- a/include/deal.II/numerics/vector_tools_boundary.h
+++ b/include/deal.II/numerics/vector_tools_boundary.h
@@ -385,8 +385,8 @@ namespace VectorTools
* where $\vec{n}$ is an outward normal vector.
*
* This function throws an exception if used with $H_\text{curl}$ conforming
- * elements, so the project_boundary_values_curl_conforming() should be used
- * instead.
+ * elements, so the project_boundary_values_curl_conforming_l2() should be
+ * used instead.
*
* @param[in] mapping The mapping that will be used in the transformations
* necessary to integrate along the boundary.
@@ -548,96 +548,6 @@ namespace VectorTools
AffineConstraints &constraints,
std::vector component_mapping = {});
-
- /**
- * Compute constraints that correspond to boundary conditions of the form
- * $\vec{n}\times\vec{u}=\vec{n}\times\vec{f}$, i.e. the tangential
- * components of $u$ and $f$ shall coincide.
- *
- * If the AffineConstraints @p constraints contained values or other
- * constraints before, the new ones are added or the old ones overwritten,
- * if a node of the boundary part to be used was already in the list of
- * constraints. This is handled by using inhomogeneous constraints. Please
- * note that when combining adaptive meshes and this kind of constraints,
- * the Dirichlet conditions should be set first, and then completed by
- * hanging node constraints, in order to make sure that the discretization
- * remains consistent. See the discussion on conflicting constraints in the
- * module on
- * @ref constraints.
- *
- * This function is explicitly written to use with the FE_Nedelec elements.
- * Thus it throws an exception, if it is called with other finite elements.
- *
- * The second argument of this function denotes the first vector component
- * in the finite element that corresponds to the vector function that you
- * want to constrain. For example, if we want to solve Maxwell's equations
- * in 3d and the finite element has components $(E_x,E_y,E_z,B_x,B_y,B_z)$
- * and we want the boundary conditions
- * $\vec{n}\times\vec{B}=\vec{n}\times\vec{f}$, then @p
- * first_vector_component would be 3. Vectors are implicitly assumed to have
- * exactly dim
components that are ordered in the same way as
- * we usually order the coordinate directions, i.e. $x$-, $y$-, and finally
- * $z$-component.
- *
- * The parameter @p boundary_component corresponds to the number @p
- * boundary_id of the face. numbers::internal_face_boundary_id is an illegal
- * value, since it is reserved for interior faces.
- *
- * The last argument is denoted to compute the normal vector $\vec{n}$ at
- * the boundary points.
- *
- * Computing constraints
- *
- * To compute the constraints we use projection-based interpolation as
- * proposed in Solin, Segeth and Dolezel (Higher order finite elements,
- * Chapman&Hall, 2004) on every face located at the boundary.
- *
- * First one projects $\vec{f}$ on the lowest-order edge shape functions.
- * Then the remaining part $(I-P_0)\vec{f}$ of the function is projected on
- * the remaining higher-order edge shape functions. In the last step we
- * project $(I-P_0-P_e)\vec{f}$ on the bubble shape functions defined on the
- * face.
- *
- * @deprecated Use the project_boundary_values_curl_conforming_l2() function
- * instead of this one.
- *
- * @ingroup constraints
- *
- * @see
- * @ref GlossBoundaryIndicator "Glossary entry on boundary indicators"
- */
- template
- DEAL_II_DEPRECATED void
- project_boundary_values_curl_conforming(
- const DoFHandler & dof_handler,
- const unsigned int first_vector_component,
- const Function &boundary_function,
- const types::boundary_id boundary_component,
- AffineConstraints & constraints,
- const Mapping & mapping);
-
- /**
- * Same as above for the hp-namespace.
- *
- * @deprecated Use the project_boundary_values_curl_conforming_l2() function
- * instead of this one.
- *
- * @ingroup constraints
- *
- * @see
- * @ref GlossBoundaryIndicator "Glossary entry on boundary indicators"
- */
- template
- DEAL_II_DEPRECATED void
- project_boundary_values_curl_conforming(
- const DoFHandler & dof_handler,
- const unsigned int first_vector_component,
- const Function & boundary_function,
- const types::boundary_id boundary_component,
- AffineConstraints & constraints,
- const hp::MappingCollection &mapping_collection =
- hp::StaticMappingQ1::mapping_collection);
-
/**
* This function is an updated version of the
* project_boundary_values_curl_conforming function. The intention is to fix
diff --git a/include/deal.II/numerics/vector_tools_boundary.templates.h b/include/deal.II/numerics/vector_tools_boundary.templates.h
index 6fb5ec1c587a..e7e456aedd55 100644
--- a/include/deal.II/numerics/vector_tools_boundary.templates.h
+++ b/include/deal.II/numerics/vector_tools_boundary.templates.h
@@ -986,1117 +986,6 @@ namespace VectorTools
}
- namespace internals
- {
- // This function computes the
- // projection of the boundary
- // function on edges for 3D.
- template
- void
- compute_edge_projection(const cell_iterator &cell,
- const unsigned int face,
- const unsigned int line,
- hp::FEValues<3> & hp_fe_values,
- const Function<3> & boundary_function,
- const unsigned int first_vector_component,
- std::vector &dof_values,
- std::vector & dofs_processed)
- {
- const double tol =
- 0.5 * cell->face(face)->line(line)->diameter() / cell->get_fe().degree;
- const unsigned int dim = 3;
- const unsigned int spacedim = 3;
-
- hp_fe_values.reinit(cell,
- (cell->active_fe_index() * cell->n_faces() + face) *
- GeometryInfo::lines_per_face +
- line);
-
- // Initialize the required
- // objects.
- const FEValues &fe_values = hp_fe_values.get_present_fe_values();
- const FiniteElement & fe = cell->get_fe();
- const std::vector> &jacobians =
- fe_values.get_jacobians();
- const std::vector> &quadrature_points =
- fe_values.get_quadrature_points();
-
- std::vector> tangentials(fe_values.n_quadrature_points);
- std::vector> values(fe_values.n_quadrature_points,
- Vector(fe.n_components()));
-
- // Get boundary function values
- // at quadrature points.
- AssertDimension(boundary_function.n_components, fe.n_components());
- boundary_function.vector_value_list(quadrature_points, values);
-
- const std::vector> &reference_quadrature_points =
- fe_values.get_quadrature().get_points();
- std::pair base_indices(0, 0);
-
- if (dynamic_cast *>(&cell->get_fe()) != nullptr)
- {
- unsigned int fe_index = 0;
- unsigned int fe_index_old = 0;
- unsigned int i = 0;
-
- for (; i < fe.n_base_elements(); ++i)
- {
- fe_index_old = fe_index;
- fe_index +=
- fe.element_multiplicity(i) * fe.base_element(i).n_components();
-
- if (fe_index > first_vector_component)
- break;
- }
-
- base_indices.first = i;
- base_indices.second = (first_vector_component - fe_index_old) /
- fe.base_element(i).n_components();
- }
-
- // coordinate directions of
- // the edges of the face.
- const unsigned int
- edge_coordinate_direction[GeometryInfo::faces_per_cell]
- [GeometryInfo::lines_per_face] = {
- {2, 2, 1, 1},
- {2, 2, 1, 1},
- {0, 0, 2, 2},
- {0, 0, 2, 2},
- {1, 1, 0, 0},
- {1, 1, 0, 0}};
- const FEValuesExtractors::Vector vec(first_vector_component);
-
- // The interpolation for the
- // lowest order edge shape
- // functions is just the mean
- // value of the tangential
- // components of the boundary
- // function on the edge.
- for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points;
- ++q_point)
- {
- // Therefore compute the
- // tangential of the edge at
- // the quadrature point.
- Point shifted_reference_point_1 =
- reference_quadrature_points[q_point];
- Point shifted_reference_point_2 =
- reference_quadrature_points[q_point];
-
- shifted_reference_point_1(edge_coordinate_direction[face][line]) +=
- tol;
- shifted_reference_point_2(edge_coordinate_direction[face][line]) -=
- tol;
- tangentials[q_point] =
- (0.5 *
- (fe_values.get_mapping().transform_unit_to_real_cell(
- cell, shifted_reference_point_1) -
- fe_values.get_mapping().transform_unit_to_real_cell(
- cell, shifted_reference_point_2)) /
- tol);
- tangentials[q_point] /= tangentials[q_point].norm();
-
- // Compute the degrees of
- // freedom.
- for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i)
- if (((dynamic_cast *>(&fe) != nullptr) &&
- (fe.system_to_base_index(fe.face_to_cell_index(i, face))
- .first == base_indices) &&
- (fe.base_element(base_indices.first)
- .face_to_cell_index(line * fe.degree, face) <=
- fe.system_to_base_index(fe.face_to_cell_index(i, face))
- .second) &&
- (fe.system_to_base_index(fe.face_to_cell_index(i, face))
- .second <=
- fe.base_element(base_indices.first)
- .face_to_cell_index((line + 1) * fe.degree - 1, face))) ||
- ((dynamic_cast *>(&fe) != nullptr) &&
- (line * fe.degree <= i) && (i < (line + 1) * fe.degree)))
- {
- const double tangential_solution_component =
- (values[q_point](first_vector_component) *
- tangentials[q_point][0] +
- values[q_point](first_vector_component + 1) *
- tangentials[q_point][1] +
- values[q_point](first_vector_component + 2) *
- tangentials[q_point][2]);
- dof_values[i] +=
- (fe_values.JxW(q_point) * tangential_solution_component *
- (fe_values[vec].value(fe.face_to_cell_index(i, face),
- q_point) *
- tangentials[q_point]) /
- std::sqrt(
- jacobians[q_point][0]
- [edge_coordinate_direction[face][line]] *
- jacobians[q_point][0]
- [edge_coordinate_direction[face][line]] +
- jacobians[q_point][1]
- [edge_coordinate_direction[face][line]] *
- jacobians[q_point][1]
- [edge_coordinate_direction[face][line]] +
- jacobians[q_point][2]
- [edge_coordinate_direction[face][line]] *
- jacobians[q_point][2]
- [edge_coordinate_direction[face][line]]));
-
- if (q_point == 0)
- dofs_processed[i] = true;
- }
- }
- }
-
- // dummy implementation of above
- // function for all other
- // dimensions
- template
- void
- compute_edge_projection(const cell_iterator &,
- const unsigned int,
- const unsigned int,
- hp::FEValues &,
- const Function &,
- const unsigned int,
- std::vector &,
- std::vector &)
- {
- Assert(false, ExcInternalError());
- }
-
- // This function computes the
- // projection of the boundary
- // function on the interior of
- // faces.
- template
- void
- compute_face_projection_curl_conforming(
- const cell_iterator & cell,
- const unsigned int face,
- hp::FEValues & hp_fe_values,
- const Function &boundary_function,
- const unsigned int first_vector_component,
- std::vector & dof_values,
- std::vector & dofs_processed)
- {
- const unsigned int spacedim = dim;
- hp_fe_values.reinit(cell,
- cell->active_fe_index() * cell->n_faces() + face);
- // Initialize the required
- // objects.
- const FEValues &fe_values = hp_fe_values.get_present_fe_values();
- const FiniteElement & fe = cell->get_fe();
- const std::vector> &jacobians =
- fe_values.get_jacobians();
- const std::vector> &quadrature_points =
- fe_values.get_quadrature_points();
- const unsigned int degree = fe.degree - 1;
- std::pair base_indices(0, 0);
-
- if (dynamic_cast *>(&cell->get_fe()) != nullptr)
- {
- unsigned int fe_index = 0;
- unsigned int fe_index_old = 0;
- unsigned int i = 0;
-
- for (; i < fe.n_base_elements(); ++i)
- {
- fe_index_old = fe_index;
- fe_index +=
- fe.element_multiplicity(i) * fe.base_element(i).n_components();
-
- if (fe_index > first_vector_component)
- break;
- }
-
- base_indices.first = i;
- base_indices.second = (first_vector_component - fe_index_old) /
- fe.base_element(i).n_components();
- }
-
- std::vector> values(fe_values.n_quadrature_points,
- Vector(fe.n_components()));
-
- // Get boundary function
- // values at quadrature
- // points.
- AssertDimension(boundary_function.n_components, fe.n_components());
- boundary_function.vector_value_list(quadrature_points, values);
-
- switch (dim)
- {
- case 2:
- {
- const double tol =
- 0.5 * cell->face(face)->diameter() / cell->get_fe().degree;
- std::vector> tangentials(
- fe_values.n_quadrature_points);
-
- const std::vector> &reference_quadrature_points =
- fe_values.get_quadrature().get_points();
-
- // coordinate directions
- // of the face.
- const unsigned int
- face_coordinate_direction[GeometryInfo::faces_per_cell] = {
- 1, 1, 0, 0};
- const FEValuesExtractors::Vector vec(first_vector_component);
-
- // The interpolation for
- // the lowest order face
- // shape functions is just
- // the mean value of the
- // tangential components
- // of the boundary function
- // on the edge.
- for (unsigned int q_point = 0;
- q_point < fe_values.n_quadrature_points;
- ++q_point)
- {
- // Therefore compute the
- // tangential of the
- // face at the quadrature
- // point.
- Point shifted_reference_point_1 =
- reference_quadrature_points[q_point];
- Point shifted_reference_point_2 =
- reference_quadrature_points[q_point];
-
- shifted_reference_point_1(face_coordinate_direction[face]) +=
- tol;
- shifted_reference_point_2(face_coordinate_direction[face]) -=
- tol;
- tangentials[q_point] =
- (fe_values.get_mapping().transform_unit_to_real_cell(
- cell, shifted_reference_point_1) -
- fe_values.get_mapping().transform_unit_to_real_cell(
- cell, shifted_reference_point_2)) /
- tol;
- tangentials[q_point] /= tangentials[q_point].norm();
-
- // Compute the degrees
- // of freedom.
- for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i)
- if (((dynamic_cast *>(&fe) !=
- nullptr) &&
- (fe.system_to_base_index(
- fe.face_to_cell_index(i, face))
- .first == base_indices)) ||
- (dynamic_cast *>(&fe) != nullptr))
- {
- dof_values[i] +=
- fe_values.JxW(q_point) *
- (values[q_point](first_vector_component) *
- tangentials[q_point][0] +
- values[q_point](first_vector_component + 1) *
- tangentials[q_point][1]) *
- (fe_values[vec].value(fe.face_to_cell_index(i, face),
- q_point) *
- tangentials[q_point]);
-
- if (q_point == 0)
- dofs_processed[i] = true;
- }
- }
-
- break;
- }
-
- case 3:
- {
- const FEValuesExtractors::Vector vec(first_vector_component);
- FullMatrix assembling_matrix(
- degree * fe.degree, dim * fe_values.n_quadrature_points);
- Vector assembling_vector(assembling_matrix.n());
- Vector cell_rhs(assembling_matrix.m());
- FullMatrix cell_matrix(assembling_matrix.m(),
- assembling_matrix.m());
- FullMatrix cell_matrix_inv(assembling_matrix.m(),
- assembling_matrix.m());
- Vector solution(cell_matrix.m());
-
- // Get coordinate directions
- // of the face.
- const unsigned int global_face_coordinate_directions
- [GeometryInfo<3>::faces_per_cell][2] = {
- {1, 2}, {1, 2}, {2, 0}, {2, 0}, {0, 1}, {0, 1}};
-
- // The projection is divided into two steps. In the first step we
- // project the boundary function on the horizontal shape
- // functions. Then the boundary function is projected on the
- // vertical shape functions. We begin with the horizontal shape
- // functions and set up a linear system of equations to get the
- // values for degrees of freedom associated with the interior of
- // the face.
- for (unsigned int q_point = 0;
- q_point < fe_values.n_quadrature_points;
- ++q_point)
- {
- // The right hand
- // side of the
- // corresponding problem
- // is the residual
- // of the boundary
- // function and
- // the already
- // interpolated part
- // on the edges.
- Tensor<1, dim> tmp;
-
- for (unsigned int d = 0; d < dim; ++d)
- tmp[d] = values[q_point](first_vector_component + d);
-
- for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i)
- if (((dynamic_cast *>(&fe) !=
- nullptr) &&
- (fe.system_to_base_index(
- fe.face_to_cell_index(i, face))
- .first == base_indices) &&
- (fe.base_element(base_indices.first)
- .face_to_cell_index(2 * fe.degree, face) <=
- fe.system_to_base_index(
- fe.face_to_cell_index(i, face))
- .second) &&
- (fe.system_to_base_index(
- fe.face_to_cell_index(i, face))
- .second <=
- fe.base_element(base_indices.first)
- .face_to_cell_index(4 * fe.degree - 1, face))) ||
- ((dynamic_cast *>(&fe) !=
- nullptr) &&
- (2 * fe.degree <= i) && (i < 4 * fe.degree)))
- tmp -=
- dof_values[i] *
- fe_values[vec].value(fe.face_to_cell_index(i, face),
- q_point);
-
- const double JxW = std::sqrt(
- fe_values.JxW(q_point) /
- ((jacobians[q_point][0]
- [global_face_coordinate_directions[face][0]] *
- jacobians[q_point][0]
- [global_face_coordinate_directions[face][0]] +
- jacobians[q_point][1]
- [global_face_coordinate_directions[face][0]] *
- jacobians[q_point][1]
- [global_face_coordinate_directions[face][0]] +
- jacobians[q_point][2]
- [global_face_coordinate_directions[face][0]] *
- jacobians[q_point][2]
- [global_face_coordinate_directions[face][0]]) *
- (jacobians[q_point][0]
- [global_face_coordinate_directions[face][1]] *
- jacobians[q_point][0]
- [global_face_coordinate_directions[face][1]] +
- jacobians[q_point][1]
- [global_face_coordinate_directions[face][1]] *
- jacobians[q_point][1]
- [global_face_coordinate_directions[face][1]] +
- jacobians[q_point][2]
- [global_face_coordinate_directions[face][1]] *
- jacobians[q_point][2]
- [global_face_coordinate_directions[face]
- [1]])));
-
- // In the weak form
- // the right hand
- // side function
- // is multiplicated
- // by the horizontal
- // shape functions
- // defined in the
- // interior of
- // the face.
- for (unsigned int d = 0; d < dim; ++d)
- assembling_vector(dim * q_point + d) = JxW * tmp[d];
-
- unsigned int index = 0;
-
- for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i)
- if (((dynamic_cast *>(&fe) !=
- nullptr) &&
- (fe.system_to_base_index(
- fe.face_to_cell_index(i, face))
- .first == base_indices) &&
- (fe.base_element(base_indices.first)
- .face_to_cell_index(
- GeometryInfo::lines_per_face * fe.degree,
- face) <=
- fe.system_to_base_index(
- fe.face_to_cell_index(i, face))
- .second) &&
- (fe.system_to_base_index(
- fe.face_to_cell_index(i, face))
- .second <
- fe.base_element(base_indices.first)
- .face_to_cell_index(
- (degree + GeometryInfo::lines_per_face) *
- fe.degree,
- face))) ||
- ((dynamic_cast *>(&fe) !=
- nullptr) &&
- (GeometryInfo::lines_per_face * fe.degree <= i) &&
- (i < (degree + GeometryInfo::lines_per_face) *
- fe.degree)))
- {
- const Tensor<1, dim> shape_value =
- (JxW *
- fe_values[vec].value(fe.face_to_cell_index(i, face),
- q_point));
-
- for (unsigned int d = 0; d < dim; ++d)
- assembling_matrix(index, dim * q_point + d) =
- shape_value[d];
-
- ++index;
- }
- }
-
- // Create the system matrix by multiplying the assembling matrix
- // with its transposed and the right hand side vector by
- // multiplying the assembling matrix with the assembling vector.
- // Invert the system matrix.
- assembling_matrix.mTmult(cell_matrix, assembling_matrix);
- cell_matrix_inv.invert(cell_matrix);
- assembling_matrix.vmult(cell_rhs, assembling_vector);
- cell_matrix_inv.vmult(solution, cell_rhs);
-
- // Store the computed
- // values.
- {
- unsigned int index = 0;
-
- for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i)
- if (((dynamic_cast *>(&fe) != nullptr) &&
- (fe.system_to_base_index(fe.face_to_cell_index(i, face))
- .first == base_indices) &&
- (fe.base_element(base_indices.first)
- .face_to_cell_index(
- GeometryInfo::lines_per_face * fe.degree,
- face) <=
- fe.system_to_base_index(fe.face_to_cell_index(i, face))
- .second) &&
- (fe.system_to_base_index(fe.face_to_cell_index(i, face))
- .second <
- fe.base_element(base_indices.first)
- .face_to_cell_index(
- (degree + GeometryInfo::lines_per_face) *
- fe.degree,
- face))) ||
- ((dynamic_cast *>(&fe) !=
- nullptr) &&
- (GeometryInfo::lines_per_face * fe.degree <= i) &&
- (i < (degree + GeometryInfo::lines_per_face) *
- fe.degree)))
- {
- dof_values[i] = solution(index);
- dofs_processed[i] = true;
- ++index;
- }
- }
-
- // Now we do the same as above with the vertical shape functions
- // instead of the horizontal ones.
- for (unsigned int q_point = 0;
- q_point < fe_values.n_quadrature_points;
- ++q_point)
- {
- Tensor<1, dim> tmp;
-
- for (unsigned int d = 0; d < dim; ++d)
- tmp[d] = values[q_point](first_vector_component + d);
-
- for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i)
- if (((dynamic_cast *>(&fe) !=
- nullptr) &&
- (fe.system_to_base_index(
- fe.face_to_cell_index(i, face))
- .first == base_indices) &&
- (fe.system_to_base_index(
- fe.face_to_cell_index(i, face))
- .second <=
- fe.base_element(base_indices.first)
- .face_to_cell_index(2 * fe.degree - 1, face)) &&
- (fe.system_to_base_index(
- fe.face_to_cell_index(i, face))
- .second >= fe.base_element(base_indices.first)
- .face_to_cell_index(0, face))) ||
- ((dynamic_cast *>(&fe) !=
- nullptr) &&
- (i < 2 * fe.degree)))
- tmp -=
- dof_values[i] *
- fe_values[vec].value(fe.face_to_cell_index(i, face),
- q_point);
-
- const double JxW = std::sqrt(
- fe_values.JxW(q_point) /
- ((jacobians[q_point][0]
- [global_face_coordinate_directions[face][0]] *
- jacobians[q_point][0]
- [global_face_coordinate_directions[face][0]] +
- jacobians[q_point][1]
- [global_face_coordinate_directions[face][0]] *
- jacobians[q_point][1]
- [global_face_coordinate_directions[face][0]] +
- jacobians[q_point][2]
- [global_face_coordinate_directions[face][0]] *
- jacobians[q_point][2]
- [global_face_coordinate_directions[face][0]]) *
- (jacobians[q_point][0]
- [global_face_coordinate_directions[face][1]] *
- jacobians[q_point][0]
- [global_face_coordinate_directions[face][1]] +
- jacobians[q_point][1]
- [global_face_coordinate_directions[face][1]] *
- jacobians[q_point][1]
- [global_face_coordinate_directions[face][1]] +
- jacobians[q_point][2]
- [global_face_coordinate_directions[face][1]] *
- jacobians[q_point][2]
- [global_face_coordinate_directions[face]
- [1]])));
-
- for (unsigned int d = 0; d < dim; ++d)
- assembling_vector(dim * q_point + d) = JxW * tmp[d];
-
- unsigned int index = 0;
-
- for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i)
- if (((dynamic_cast *>(&fe) !=
- nullptr) &&
- (fe.system_to_base_index(
- fe.face_to_cell_index(i, face))
- .first == base_indices) &&
- (fe.base_element(base_indices.first)
- .face_to_cell_index(
- (degree + GeometryInfo::lines_per_face) *
- fe.degree,
- face) <=
- fe.system_to_base_index(
- fe.face_to_cell_index(i, face))
- .second)) ||
- ((dynamic_cast *>(&fe) !=
- nullptr) &&
- ((degree + GeometryInfo::lines_per_face) *
- fe.degree <=
- i)))
- {
- const Tensor<1, dim> shape_value =
- JxW *
- fe_values[vec].value(fe.face_to_cell_index(i, face),
- q_point);
-
- for (unsigned int d = 0; d < dim; ++d)
- assembling_matrix(index, dim * q_point + d) =
- shape_value[d];
-
- ++index;
- }
- }
-
- assembling_matrix.mTmult(cell_matrix, assembling_matrix);
- cell_matrix_inv.invert(cell_matrix);
- assembling_matrix.vmult(cell_rhs, assembling_vector);
- cell_matrix_inv.vmult(solution, cell_rhs);
-
- unsigned int index = 0;
-
- for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i)
- if (((dynamic_cast *>(&fe) != nullptr) &&
- (fe.system_to_base_index(fe.face_to_cell_index(i, face))
- .first == base_indices) &&
- (fe.base_element(base_indices.first)
- .face_to_cell_index(
- (degree + GeometryInfo::lines_per_face) *
- fe.degree,
- face) <=
- fe.system_to_base_index(fe.face_to_cell_index(i, face))
- .second)) ||
- ((dynamic_cast *>(&fe) != nullptr) &&
- ((degree + GeometryInfo::lines_per_face) *
- fe.degree <=
- i)))
- {
- dof_values[i] = solution(index);
- dofs_processed[i] = true;
- ++index;
- }
-
- break;
- }
-
- default:
- Assert(false, ExcNotImplemented());
- }
- }
- } // namespace internals
-
-
-
- template
- void
- project_boundary_values_curl_conforming(
- const DoFHandler & dof_handler,
- const unsigned int first_vector_component,
- const Function & boundary_function,
- const types::boundary_id boundary_component,
- AffineConstraints &constraints,
- const Mapping & mapping)
- {
- // Projection-based interpolation is performed in two (in 2D) respectively
- // three (in 3D) steps. First the tangential component of the function is
- // interpolated on each edge. This gives the values for the degrees of
- // freedom corresponding to the edge shape functions. Now we are done for
- // 2D, but in 3D we possibly have also degrees of freedom, which are
- // located in the interior of the faces. Therefore we compute the residual
- // of the function describing the boundary values and the interpolated
- // part, which we have computed in the last step. On the faces there are
- // two kinds of shape functions, the horizontal and the vertical
- // ones. Thus we have to solve two linear systems of equations of size
- // degree * (degree + 1) to obtain the values for the
- // corresponding degrees of freedom.
- const unsigned int superdegree = dof_handler.get_fe().degree;
- const QGauss reference_face_quadrature(2 * superdegree);
-
- // TODO: the implementation makes the assumption that all faces have the
- // same number of dofs
- AssertDimension(dof_handler.get_fe().n_unique_faces(), 1);
- const unsigned int dofs_per_face = dof_handler.get_fe().n_dofs_per_face(0);
-
- const hp::FECollection &fe_collection(dof_handler.get_fe_collection());
- const hp::MappingCollection mapping_collection(mapping);
- hp::QCollection face_quadrature_collection;
-
- for (unsigned int face : GeometryInfo::face_indices())
- face_quadrature_collection.push_back(
- QProjector::project_to_face(reference_face_quadrature, face));
-
- hp::FEValues fe_face_values(mapping_collection,
- fe_collection,
- face_quadrature_collection,
- update_jacobians | update_JxW_values |
- update_quadrature_points |
- update_values);
-
- std::vector dofs_processed(dofs_per_face);
- std::vector dof_values(dofs_per_face);
- std::vector face_dof_indices(dofs_per_face);
- typename DoFHandler::active_cell_iterator cell =
- dof_handler.begin_active();
-
- switch (dim)
- {
- case 2:
- {
- for (; cell != dof_handler.end(); ++cell)
- if (cell->at_boundary() && cell->is_locally_owned())
- for (const unsigned int face : cell->face_indices())
- if (cell->face(face)->boundary_id() == boundary_component)
- {
- // if the FE is a
- // FE_Nothing object
- // there is no work to
- // do
- if (dynamic_cast *>(
- &cell->get_fe()) != nullptr)
- return;
-
- // This is only
- // implemented, if the
- // FE is a Nedelec
- // element. If the FE
- // is a FESystem, we
- // cannot check this.
- if (dynamic_cast *>(
- &cell->get_fe()) == nullptr)
- {
- AssertThrow(
- dynamic_cast *>(
- &cell->get_fe()) != nullptr,
- (typename FiniteElement<
- dim>::ExcInterpolationNotImplemented()));
- }
-
- for (unsigned int dof = 0; dof < dofs_per_face; ++dof)
- {
- dof_values[dof] = 0.0;
- dofs_processed[dof] = false;
- }
-
- // Compute the
- // projection of the
- // boundary function on
- // the edge.
- internals::compute_face_projection_curl_conforming(
- cell,
- face,
- fe_face_values,
- boundary_function,
- first_vector_component,
- dof_values,
- dofs_processed);
- cell->face(face)->get_dof_indices(
- face_dof_indices, cell->active_fe_index());
-
- // Add the computed constraints to the constraints
- // object, if the degree of freedom is not already
- // constrained.
- for (unsigned int dof = 0; dof < dofs_per_face; ++dof)
- if (dofs_processed[dof] &&
- constraints.can_store_line(face_dof_indices[dof]) &&
- !(constraints.is_constrained(
- face_dof_indices[dof])))
- {
- constraints.add_line(face_dof_indices[dof]);
-
- if (std::abs(dof_values[dof]) > 1e-13)
- constraints.set_inhomogeneity(
- face_dof_indices[dof], dof_values[dof]);
- }
- }
-
- break;
- }
-
- case 3:
- {
- const QGauss reference_edge_quadrature(2 * superdegree);
- const unsigned int degree = superdegree - 1;
- hp::QCollection edge_quadrature_collection;
-
- for (const unsigned int face : GeometryInfo::face_indices())
- for (unsigned int line = 0;
- line < GeometryInfo::lines_per_face;
- ++line)
- edge_quadrature_collection.push_back(
- QProjector::project_to_face(
- QProjector::project_to_face(
- reference_edge_quadrature, line),
- face));
-
- hp::FEValues fe_edge_values(mapping_collection,
- fe_collection,
- edge_quadrature_collection,
- update_jacobians |
- update_JxW_values |
- update_quadrature_points |
- update_values);
-
- for (; cell != dof_handler.end(); ++cell)
- if (cell->at_boundary() && cell->is_locally_owned())
- for (const unsigned int face : cell->face_indices())
- if (cell->face(face)->boundary_id() == boundary_component)
- {
- // if the FE is a
- // FE_Nothing object
- // there is no work to
- // do
- if (dynamic_cast *>(
- &cell->get_fe()) != nullptr)
- return;
-
- // This is only
- // implemented, if the
- // FE is a Nedelec
- // element. If the FE is
- // a FESystem we cannot
- // check this.
- if (dynamic_cast *>(
- &cell->get_fe()) == nullptr)
- {
- AssertThrow(dynamic_cast *>(
- &cell->get_fe()) != nullptr,
- typename FiniteElement<
- dim>::ExcInterpolationNotImplemented());
- }
-
- for (unsigned int dof = 0; dof < dofs_per_face; ++dof)
- {
- dof_values[dof] = 0.0;
- dofs_processed[dof] = false;
- }
-
- // First we compute the
- // projection on the
- // edges.
- for (unsigned int line = 0;
- line < GeometryInfo<3>::lines_per_face;
- ++line)
- internals::compute_edge_projection(
- cell,
- face,
- line,
- fe_edge_values,
- boundary_function,
- first_vector_component,
- dof_values,
- dofs_processed);
-
- // If there are higher
- // order shape
- // functions, there is
- // still some work
- // left.
- if (degree > 0)
- internals::compute_face_projection_curl_conforming(
- cell,
- face,
- fe_face_values,
- boundary_function,
- first_vector_component,
- dof_values,
- dofs_processed);
-
- // Store the computed
- // values in the global
- // vector.
- cell->face(face)->get_dof_indices(
- face_dof_indices, cell->active_fe_index());
-
- for (unsigned int dof = 0; dof < dofs_per_face; ++dof)
- if (dofs_processed[dof] &&
- constraints.can_store_line(face_dof_indices[dof]) &&
- !(constraints.is_constrained(
- face_dof_indices[dof])))
- {
- constraints.add_line(face_dof_indices[dof]);
-
- if (std::abs(dof_values[dof]) > 1e-13)
- constraints.set_inhomogeneity(
- face_dof_indices[dof], dof_values[dof]);
- }
- }
-
- break;
- }
-
- default:
- Assert(false, ExcNotImplemented());
- }
- }
-
-
-
- template
- void
-
- project_boundary_values_curl_conforming(
- const DoFHandler & dof_handler,
- const unsigned int first_vector_component,
- const Function & boundary_function,
- const types::boundary_id boundary_component,
- AffineConstraints & constraints,
- const hp::MappingCollection &mapping_collection)
- {
- const hp::FECollection &fe_collection(dof_handler.get_fe_collection());
- hp::QCollection face_quadrature_collection;
-
- for (unsigned int i = 0; i < fe_collection.size(); ++i)
- {
- const QGauss reference_face_quadrature(
- 2 * fe_collection[i].degree);
-
- for (unsigned int face : GeometryInfo::face_indices())
- face_quadrature_collection.push_back(
- QProjector::project_to_face(reference_face_quadrature, face));
- }
-
- hp::FEValues fe_face_values(mapping_collection,
- fe_collection,
- face_quadrature_collection,
- update_jacobians | update_JxW_values |
- update_quadrature_points |
- update_values);
- std::vector dofs_processed;
- std::vector dof_values;
- std::vector face_dof_indices;
- typename DoFHandler::active_cell_iterator cell =
- dof_handler.begin_active();
-
- switch (dim)
- {
- case 2:
- {
- for (; cell != dof_handler.end(); ++cell)
- if (cell->at_boundary() && cell->is_locally_owned())
- for (const unsigned int face : cell->face_indices())
- if (cell->face(face)->boundary_id() == boundary_component)
- {
- // if the FE is a FE_Nothing object there is no work to do
- if (dynamic_cast *>(
- &cell->get_fe()) != nullptr)
- return;
-
- // This is only implemented, if the FE is a Nedelec
- // element. If the FE is a FESystem we cannot check this.
- if (dynamic_cast *>(
- &cell->get_fe()) == nullptr)
- {
- AssertThrow(dynamic_cast *>(
- &cell->get_fe()) != nullptr,
- typename FiniteElement<
- dim>::ExcInterpolationNotImplemented());
- }
-
- const unsigned int dofs_per_face =
- cell->get_fe().n_dofs_per_face(face);
-
- dofs_processed.resize(dofs_per_face);
- dof_values.resize(dofs_per_face);
-
- for (unsigned int dof = 0; dof < dofs_per_face; ++dof)
- {
- dof_values[dof] = 0.0;
- dofs_processed[dof] = false;
- }
-
- internals::compute_face_projection_curl_conforming(
- cell,
- face,
- fe_face_values,
- boundary_function,
- first_vector_component,
- dof_values,
- dofs_processed);
- face_dof_indices.resize(dofs_per_face);
- cell->face(face)->get_dof_indices(
- face_dof_indices, cell->active_fe_index());
-
- for (unsigned int dof = 0; dof < dofs_per_face; ++dof)
- if (dofs_processed[dof] &&
- constraints.can_store_line(face_dof_indices[dof]) &&
- !(constraints.is_constrained(
- face_dof_indices[dof])))
- {
- constraints.add_line(face_dof_indices[dof]);
-
- if (std::abs(dof_values[dof]) > 1e-13)
- constraints.set_inhomogeneity(
- face_dof_indices[dof], dof_values[dof]);
- }
- }
-
- break;
- }
-
- case 3:
- {
- hp::QCollection edge_quadrature_collection;
-
- for (unsigned int i = 0; i < fe_collection.size(); ++i)
- {
- const QGauss reference_edge_quadrature(
- 2 * fe_collection[i].degree);
-
- for (const unsigned int face :
- GeometryInfo::face_indices())
- for (unsigned int line = 0;
- line < GeometryInfo::lines_per_face;
- ++line)
- edge_quadrature_collection.push_back(
- QProjector::project_to_face(
- QProjector::project_to_face(
- reference_edge_quadrature, line),
- face));
- }
-
- hp::FEValues fe_edge_values(mapping_collection,
- fe_collection,
- edge_quadrature_collection,
- update_jacobians |
- update_JxW_values |
- update_quadrature_points |
- update_values);
-
- for (; cell != dof_handler.end(); ++cell)
- if (cell->at_boundary() && cell->is_locally_owned())
- for (const unsigned int face : cell->face_indices())
- if (cell->face(face)->boundary_id() == boundary_component)
- {
- // if the FE is a FE_Nothing object there is no work to do
- if (dynamic_cast *>(
- &cell->get_fe()) != nullptr)
- return;
-
- // This is only implemented, if the FE is a Nedelec
- // element. If the FE is a FESystem we cannot check this.
- if (dynamic_cast *>(
- &cell->get_fe()) == nullptr)
- {
- AssertThrow(dynamic_cast *>(
- &cell->get_fe()) != nullptr,
- typename FiniteElement<
- dim>::ExcInterpolationNotImplemented());
- }
-
- const unsigned int superdegree = cell->get_fe().degree;
- const unsigned int degree = superdegree - 1;
- const unsigned int dofs_per_face =
- cell->get_fe().n_dofs_per_face(face);
-
- dofs_processed.resize(dofs_per_face);
- dof_values.resize(dofs_per_face);
-
- for (unsigned int dof = 0; dof < dofs_per_face; ++dof)
- {
- dof_values[dof] = 0.0;
- dofs_processed[dof] = false;
- }
-
- for (unsigned int line = 0;
- line < GeometryInfo::lines_per_face;
- ++line)
- internals::compute_edge_projection(
- cell,
- face,
- line,
- fe_edge_values,
- boundary_function,
- first_vector_component,
- dof_values,
- dofs_processed);
-
- // If there are higher order shape functions, there is
- // still some work left.
- if (degree > 0)
- internals::compute_face_projection_curl_conforming(
- cell,
- face,
- fe_face_values,
- boundary_function,
- first_vector_component,
- dof_values,
- dofs_processed);
-
-
- face_dof_indices.resize(dofs_per_face);
- cell->face(face)->get_dof_indices(
- face_dof_indices, cell->active_fe_index());
-
- for (unsigned int dof = 0; dof < dofs_per_face; ++dof)
- if (dofs_processed[dof] &&
- constraints.can_store_line(face_dof_indices[dof]) &&
- !(constraints.is_constrained(
- face_dof_indices[dof])))
- {
- constraints.add_line(face_dof_indices[dof]);
-
- if (std::abs(dof_values[dof]) > 1e-13)
- constraints.set_inhomogeneity(
- face_dof_indices[dof], dof_values[dof]);
- }
- }
-
- break;
- }
-
- default:
- Assert(false, ExcNotImplemented());
- }
- }
-
-
namespace internals
{
template
diff --git a/include/deal.II/numerics/vector_tools_integrate_difference.h b/include/deal.II/numerics/vector_tools_integrate_difference.h
index 09b07f2ba0da..542051608382 100644
--- a/include/deal.II/numerics/vector_tools_integrate_difference.h
+++ b/include/deal.II/numerics/vector_tools_integrate_difference.h
@@ -197,87 +197,6 @@ namespace VectorTools
const Function * weight = nullptr,
const double exponent = 2.);
- /**
- * Compute the cellwise error of the finite element solution. Integrate the
- * difference between a reference function which is given as a continuous
- * function object, and a finite element function. The result of this
- * function is the vector @p difference that contains one value per active
- * cell $K$ of the triangulation. Each of the values of this vector $d$
- * equals
- * @f{align*}{
- * d_K = \| u-u_h \|_X
- * @f}
- * where $X$ denotes the norm chosen and $u$ represents the exact solution.
- *
- * @deprecated Use integrate_difference(const Mapping &, const DoFHandler &, const InVector &, const Function &, OutVector &, const Quadrature &, const NormType &, const Function *, const double) instead.
- */
- template
- DEAL_II_DEPRECATED typename std::enable_if<
- !std::is_same::value>::type
- integrate_difference(const Mapping & mapping,
- const DoFHandler & dof,
- const InVector & fe_function,
- const Function &exact_solution,
- OutVector & difference,
- const Quadrature & q,
- const NormType & norm,
- const Function *weight = nullptr,
- const double exponent = 2.);
-
- /**
- * Call the integrate_difference() function, see above, with
- * mapping=MappingQGeneric@(1).
- *
- * @deprecated Use integrate_difference(const DoFHandler &, const InVector &, const Function &exact_solution, OutVector &, const Quadrature &, const NormType &, const Function *, const double) instead.
- */
- template
- DEAL_II_DEPRECATED typename std::enable_if<
- !std::is_same::value>::type
- integrate_difference(const DoFHandler & dof,
- const InVector & fe_function,
- const Function &exact_solution,
- OutVector & difference,
- const Quadrature & q,
- const NormType & norm,
- const Function *weight = nullptr,
- const double exponent = 2.);
-
- /**
- * Same as above for hp.
- *
- * @deprecated Use integrate_difference(const hp::MappingCollection &, const DoFHandler &, const InVector &, const Function &, OutVector &, const hp::QCollection &, const NormType &, const Function *, const double) instead.
- */
- template
- DEAL_II_DEPRECATED typename std::enable_if<
- !std::is_same::value>::type
- integrate_difference(const hp::MappingCollection &mapping,
- const DoFHandler & dof,
- const InVector & fe_function,
- const Function &exact_solution,
- OutVector & difference,
- const hp::QCollection & q,
- const NormType & norm,
- const Function *weight = nullptr,
- const double exponent = 2.);
-
- /**
- * Call the integrate_difference() function, see above, with
- * mapping=MappingQGeneric@(1).
- *
- * @deprecated Use integrate_difference(const DoFHandler &, const InVector &, const Function &, OutVector &, const hp::QCollection &, const NormType &, const Function *, const double) instead.
- */
- template
- DEAL_II_DEPRECATED typename std::enable_if<
- !std::is_same::value>::type
- integrate_difference(const DoFHandler & dof,
- const InVector & fe_function,
- const Function &exact_solution,
- OutVector & difference,
- const hp::QCollection & q,
- const NormType & norm,
- const Function *weight = nullptr,
- const double exponent = 2.);
-
/**
* Take a Vector @p cellwise_error of errors on each cell with
* tria.n_active_cells() entries and return the global
diff --git a/include/deal.II/numerics/vector_tools_integrate_difference.templates.h b/include/deal.II/numerics/vector_tools_integrate_difference.templates.h
index 19bf68dace23..067208a9a6ec 100644
--- a/include/deal.II/numerics/vector_tools_integrate_difference.templates.h
+++ b/include/deal.II/numerics/vector_tools_integrate_difference.templates.h
@@ -113,82 +113,6 @@ namespace VectorTools
n_q_points, std::vector>(n_components));
}
- template
- struct DEAL_II_DEPRECATED DeprecatedIDScratchData
- {
- DeprecatedIDScratchData(
- const dealii::hp::MappingCollection &mapping,
- const dealii::hp::FECollection & fe,
- const dealii::hp::QCollection & q,
- const UpdateFlags update_flags);
-
- DeprecatedIDScratchData(const DeprecatedIDScratchData &data);
-
- void
- resize_vectors(const unsigned int n_q_points,
- const unsigned int n_components);
-
- std::vector> function_values;
- std::vector>> function_grads;
- std::vector weight_values;
- std::vector> weight_vectors;
-
- std::vector> psi_values;
- std::vector>> psi_grads;
- std::vector psi_scalar;
-
- std::vector tmp_values;
- std::vector> tmp_vector_values;
- std::vector> tmp_gradients;
- std::vector>> tmp_vector_gradients;
-
- dealii::hp::FEValues x_fe_values;
- };
-
-
- template
- DeprecatedIDScratchData::DeprecatedIDScratchData(
- const dealii::hp::MappingCollection &mapping,
- const dealii::hp::FECollection & fe,
- const dealii::hp::QCollection & q,
- const UpdateFlags update_flags)
- : x_fe_values(mapping, fe, q, update_flags)
- {}
-
- template
- DeprecatedIDScratchData::DeprecatedIDScratchData(
- const DeprecatedIDScratchData &data)
- : x_fe_values(data.x_fe_values.get_mapping_collection(),
- data.x_fe_values.get_fe_collection(),
- data.x_fe_values.get_quadrature_collection(),
- data.x_fe_values.get_update_flags())
- {}
-
- template
- void
- DeprecatedIDScratchData::resize_vectors(
- const unsigned int n_q_points,
- const unsigned int n_components)
- {
- function_values.resize(n_q_points, Vector(n_components));
- function_grads.resize(
- n_q_points, std::vector>(n_components));
-
- weight_values.resize(n_q_points);
- weight_vectors.resize(n_q_points, Vector(n_components));
-
- psi_values.resize(n_q_points, Vector(n_components));
- psi_grads.resize(n_q_points,
- std::vector>(n_components));
- psi_scalar.resize(n_q_points);
-
- tmp_values.resize(n_q_points);
- tmp_vector_values.resize(n_q_points, Vector(n_components));
- tmp_gradients.resize(n_q_points);
- tmp_vector_gradients.resize(
- n_q_points, std::vector>(n_components));
- }
-
namespace internal
{
template
@@ -504,296 +428,6 @@ namespace VectorTools
return diff;
}
- template
- DEAL_II_DEPRECATED
- typename std::enable_if::value,
- double>::type
- integrate_difference_inner(
- const Function & exact_solution,
- const NormType & norm,
- const Function * weight,
- const UpdateFlags update_flags,
- const double exponent,
- const unsigned int n_components,
- DeprecatedIDScratchData &data)
- {
- const bool fe_is_system = (n_components != 1);
- const dealii::FEValues &fe_values =
- data.x_fe_values.get_present_fe_values();
- const unsigned int n_q_points = fe_values.n_quadrature_points;
-
- if (weight != nullptr)
- {
- if (weight->n_components > 1)
- weight->vector_value_list(fe_values.get_quadrature_points(),
- data.weight_vectors);
- else
- {
- weight->value_list(fe_values.get_quadrature_points(),
- data.weight_values);
- for (unsigned int k = 0; k < n_q_points; ++k)
- data.weight_vectors[k] = data.weight_values[k];
- }
- }
- else
- {
- for (unsigned int k = 0; k < n_q_points; ++k)
- data.weight_vectors[k] = 1.;
- }
-
-
- if (update_flags & update_values)
- {
- // first compute the exact solution (vectors) at the quadrature
- // points. try to do this as efficient as possible by avoiding a
- // second virtual function call in case the function really has only
- // one component
- //
- // TODO: we have to work a bit here because the Function
- // interface of the argument denoting the exact function only
- // provides us with double/Tensor<1,dim> values, rather than
- // with the correct data type. so evaluate into a temp
- // object, then copy around
- if (fe_is_system)
- {
- exact_solution.vector_value_list(
- fe_values.get_quadrature_points(), data.tmp_vector_values);
- for (unsigned int i = 0; i < n_q_points; ++i)
- data.psi_values[i] = data.tmp_vector_values[i];
- }
- else
- {
- exact_solution.value_list(fe_values.get_quadrature_points(),
- data.tmp_values);
- for (unsigned int i = 0; i < n_q_points; ++i)
- data.psi_values[i](0) = data.tmp_values[i];
- }
-
- // then subtract finite element fe_function
- for (unsigned int q = 0; q < n_q_points; ++q)
- for (unsigned int i = 0; i < data.psi_values[q].size(); ++i)
- data.psi_values[q][i] -= data.function_values[q][i];
- }
-
- // Do the same for gradients, if required
- if (update_flags & update_gradients)
- {
- // try to be a little clever to avoid recursive virtual function
- // calls when calling gradient_list for functions that are really
- // scalar functions
- if (fe_is_system)
- {
- exact_solution.vector_gradient_list(
- fe_values.get_quadrature_points(), data.tmp_vector_gradients);
- for (unsigned int i = 0; i < n_q_points; ++i)
- for (unsigned int comp = 0; comp < data.psi_grads[i].size();
- ++comp)
- data.psi_grads[i][comp] = data.tmp_vector_gradients[i][comp];
- }
- else
- {
- exact_solution.gradient_list(fe_values.get_quadrature_points(),
- data.tmp_gradients);
- for (unsigned int i = 0; i < n_q_points; ++i)
- data.psi_grads[i][0] = data.tmp_gradients[i];
- }
-
- // then subtract finite element function_grads. We need to be
- // careful in the codimension one case, since there we only have
- // tangential gradients in the finite element function, not the full
- // gradient. This is taken care of, by subtracting the normal
- // component of the gradient from the exact function.
- if (update_flags & update_normal_vectors)
- for (unsigned int k = 0; k < n_components; ++k)
- for (unsigned int q = 0; q < n_q_points; ++q)
- {
- // compute (f.n) n
- const typename ProductType::type f_dot_n =
- data.psi_grads[q][k] * fe_values.normal_vector(q);
- const Tensor<1, spacedim, Number> f_dot_n_times_n =
- f_dot_n * fe_values.normal_vector(q);
-
- data.psi_grads[q][k] -=
- (data.function_grads[q][k] + f_dot_n_times_n);
- }
- else
- for (unsigned int k = 0; k < n_components; ++k)
- for (unsigned int q = 0; q < n_q_points; ++q)
- for (unsigned int d = 0; d < spacedim; ++d)
- data.psi_grads[q][k][d] -= data.function_grads[q][k][d];
- }
-
- double diff = 0;
- Number diff_mean = 0;
-
- // First work on function values:
- switch (norm)
- {
- case mean:
- // Compute values in quadrature points and integrate
- for (unsigned int q = 0; q < n_q_points; ++q)
- {
- Number sum = 0;
- for (unsigned int k = 0; k < n_components; ++k)
- if (data.weight_vectors[q](k) != 0)
- sum += data.psi_values[q](k) * data.weight_vectors[q](k);
- diff_mean += sum * fe_values.JxW(q);
- }
- break;
-
- case Lp_norm:
- case L1_norm:
- case W1p_norm:
- // Compute values in quadrature points and integrate
- for (unsigned int q = 0; q < n_q_points; ++q)
- {
- double sum = 0;
- for (unsigned int k = 0; k < n_components; ++k)
- if (data.weight_vectors[q](k) != 0)
- sum += std::pow(static_cast(
- numbers::NumberTraits::abs_square(
- data.psi_values[q](k))),
- exponent / 2.) *
- data.weight_vectors[q](k);
- diff += sum * fe_values.JxW(q);
- }
-
- // Compute the root only if no derivative values are added later
- if (!(update_flags & update_gradients))
- diff = std::pow(diff, 1. / exponent);
- break;
-
- case L2_norm:
- case H1_norm:
- // Compute values in quadrature points and integrate
- for (unsigned int q = 0; q < n_q_points; ++q)
- {
- double sum = 0;
- for (unsigned int k = 0; k < n_components; ++k)
- if (data.weight_vectors[q](k) != 0)
- sum += numbers::NumberTraits::abs_square(
- data.psi_values[q](k)) *
- data.weight_vectors[q](k);
- diff += sum * fe_values.JxW(q);
- }
- // Compute the root only, if no derivative values are added later
- if (norm == L2_norm)
- diff = std::sqrt(diff);
- break;
-
- case Linfty_norm:
- case W1infty_norm:
- for (unsigned int q = 0; q < n_q_points; ++q)
- for (unsigned int k = 0; k < n_components; ++k)
- if (data.weight_vectors[q](k) != 0)
- diff = std::max(diff,
- double(std::abs(data.psi_values[q](k) *
- data.weight_vectors[q](k))));
- break;
-
- case H1_seminorm:
- case Hdiv_seminorm:
- case W1p_seminorm:
- case W1infty_seminorm:
- // function values are not used for these norms
- break;
-
- default:
- Assert(false, ExcNotImplemented());
- break;
- }
-
- // Now compute terms depending on derivatives:
- switch (norm)
- {
- case W1p_seminorm:
- case W1p_norm:
- for (unsigned int q = 0; q < n_q_points; ++q)
- {
- double sum = 0;
- for (unsigned int k = 0; k < n_components; ++k)
- if (data.weight_vectors[q](k) != 0)
- sum += std::pow(data.psi_grads[q][k].norm_square(),
- exponent / 2.) *
- data.weight_vectors[q](k);
- diff += sum * fe_values.JxW(q);
- }
- diff = std::pow(diff, 1. / exponent);
- break;
-
- case H1_seminorm:
- case H1_norm:
- for (unsigned int q = 0; q < n_q_points; ++q)
- {
- double sum = 0;
- for (unsigned int k = 0; k < n_components; ++k)
- if (data.weight_vectors[q](k) != 0)
- sum += data.psi_grads[q][k].norm_square() *
- data.weight_vectors[q](k);
- diff += sum * fe_values.JxW(q);
- }
- diff = std::sqrt(diff);
- break;
-
- case Hdiv_seminorm:
- for (unsigned int q = 0; q < n_q_points; ++q)
- {
- unsigned int idx = 0;
- if (weight != nullptr)
- for (; idx < n_components; ++idx)
- if (data.weight_vectors[0](idx) > 0)
- break;
-
- Assert(
- n_components >= idx + dim,
- ExcMessage(
- "You can only ask for the Hdiv norm for a finite element "
- "with at least 'dim' components. In that case, this function "
- "will find the index of the first non-zero weight and take "
- "the divergence of the 'dim' components that follow it."));
-
- Number sum = 0;
- // take the trace of the derivatives scaled by the weight and
- // square it
- for (unsigned int k = idx; k < idx + dim; ++k)
- if (data.weight_vectors[q](k) != 0)
- sum += data.psi_grads[q][k][k - idx] *
- std::sqrt(data.weight_vectors[q](k));
- diff += numbers::NumberTraits::abs_square(sum) *
- fe_values.JxW(q);
- }
- diff = std::sqrt(diff);
- break;
-
- case W1infty_seminorm:
- case W1infty_norm:
- {
- double t = 0;
- for (unsigned int q = 0; q < n_q_points; ++q)
- for (unsigned int k = 0; k < n_components; ++k)
- if (data.weight_vectors[q](k) != 0)
- for (unsigned int d = 0; d < dim; ++d)
- t = std::max(t,
- double(std::abs(data.psi_grads[q][k][d]) *
- data.weight_vectors[q](k)));
-
- // then add seminorm to norm if that had previously been
- // computed
- diff += t;
- }
- break;
- default:
- break;
- }
-
- if (norm == mean)
- diff = internal::mean_to_double(diff_mean);
-
- // append result of this cell to the end of the vector
- AssertIsFinite(diff);
- return diff;
- }
-
template
@@ -913,124 +547,6 @@ namespace VectorTools
difference(cell->active_cell_index()) = 0;
}
- template
- DEAL_II_DEPRECATED static typename std::enable_if<
- !std::is_same::value>::type
- do_integrate_difference(
- const dealii::hp::MappingCollection &mapping,
- const DoFHandler & dof,
- const InVector & fe_function,
- const Function & exact_solution,
- OutVector & difference,
- const dealii::hp::QCollection & q,
- const NormType & norm,
- const Function * weight,
- const double exponent_1)
- {
- using Number = typename InVector::value_type;
- // we mark the "exponent" parameter to this function "const" since it is
- // strictly incoming, but we need to set it to something different later
- // on, if necessary, so have a read-write version of it:
- double exponent = exponent_1;
-
- const unsigned int n_components = dof.get_fe(0).n_components();
-
- Assert(exact_solution.n_components == n_components,
- ExcDimensionMismatch(exact_solution.n_components, n_components));
-
- if (weight != nullptr)
- {
- Assert((weight->n_components == 1) ||
- (weight->n_components == n_components),
- ExcDimensionMismatch(weight->n_components, n_components));
- }
-
- difference.reinit(dof.get_triangulation().n_active_cells());
-
- switch (norm)
- {
- case L2_norm:
- case H1_seminorm:
- case H1_norm:
- case Hdiv_seminorm:
- exponent = 2.;
- break;
-
- case L1_norm:
- exponent = 1.;
- break;
-
- default:
- break;
- }
-
- UpdateFlags update_flags =
- UpdateFlags(update_quadrature_points | update_JxW_values);
- switch (norm)
- {
- case H1_seminorm:
- case Hdiv_seminorm:
- case W1p_seminorm:
- case W1infty_seminorm:
- update_flags |= UpdateFlags(update_gradients);
- if (spacedim == dim + 1)
- update_flags |= UpdateFlags(update_normal_vectors);
-
- break;
-
- case H1_norm:
- case W1p_norm:
- case W1infty_norm:
- update_flags |= UpdateFlags(update_gradients);
- if (spacedim == dim + 1)
- update_flags |= UpdateFlags(update_normal_vectors);
- DEAL_II_FALLTHROUGH;
-
- default:
- update_flags |= UpdateFlags(update_values);
- break;
- }
-
- const dealii::hp::FECollection &fe_collection =
- dof.get_fe_collection();
- DeprecatedIDScratchData data(mapping,
- fe_collection,
- q,
- update_flags);
-
- // loop over all cells
- for (const auto &cell : dof.active_cell_iterators())
- if (cell->is_locally_owned())
- {
- // initialize for this cell
- data.x_fe_values.reinit(cell);
-
- const dealii::FEValues &fe_values =
- data.x_fe_values.get_present_fe_values();
- const unsigned int n_q_points = fe_values.n_quadrature_points;
- data.resize_vectors(n_q_points, n_components);
-
- if (update_flags & update_values)
- fe_values.get_function_values(fe_function, data.function_values);
- if (update_flags & update_gradients)
- fe_values.get_function_gradients(fe_function,
- data.function_grads);
-
- difference(cell->active_cell_index()) =
- integrate_difference_inner(exact_solution,
- norm,
- weight,
- update_flags,
- exponent,
- n_components,
- data);
- }
- else
- // the cell is a ghost cell or is artificial. write a zero into the
- // corresponding value of the returned vector
- difference(cell->active_cell_index()) = 0;
- }
-
} // namespace internal
template
@@ -1058,31 +574,6 @@ namespace VectorTools
exponent);
}
- template
- DEAL_II_DEPRECATED typename std::enable_if<
- !std::is_same::value>::type
- integrate_difference(const Mapping