From: Daniel Arndt Date: Fri, 4 Jun 2021 14:58:24 +0000 (-0400) Subject: Remove deprecated VectorTools X-Git-Tag: v9.4.0-rc1~1278^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F12389%2Fhead;p=dealii.git Remove deprecated VectorTools --- diff --git a/doc/news/changes/incompatibilities/20210604DanielArndt b/doc/news/changes/incompatibilities/20210604DanielArndt new file mode 100644 index 0000000000..9e5d6d39e9 --- /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 273a1d0a2b..056d53c915 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 6fb5ec1c58..e7e456aedd 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 09b07f2ba0..5420516083 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 19bf68dace..067208a9a6 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 & mapping, - const DoFHandler &dof, - const InVector & fe_function, - const Function & exact_solution, - OutVector & difference, - const Quadrature & q, - const NormType & norm, - const Function * weight, - const double exponent) - { - internal::do_integrate_difference(hp::MappingCollection( - mapping), - dof, - fe_function, - exact_solution, - difference, - hp::QCollection(q), - norm, - weight, - exponent); - } - template void @@ -1109,32 +600,6 @@ namespace VectorTools } - 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, - const double exponent) - { - internal::do_integrate_difference( - hp::StaticMappingQ1::mapping_collection, - dof, - fe_function, - exact_solution, - difference, - hp::QCollection(q), - norm, - weight, - exponent); - } - - - template void integrate_difference( @@ -1160,33 +625,6 @@ namespace VectorTools exponent); } - template - DEAL_II_DEPRECATED typename std::enable_if< - !std::is_same::value>::type - integrate_difference( - const dealii::hp::MappingCollection &mapping, - const dealii::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) - { - internal::do_integrate_difference(hp::MappingCollection( - mapping), - dof, - fe_function, - exact_solution, - difference, - q, - norm, - weight, - exponent); - } - - template void integrate_difference( @@ -1211,30 +649,6 @@ namespace VectorTools exponent); } - template - DEAL_II_DEPRECATED typename std::enable_if< - !std::is_same::value>::type - integrate_difference(const dealii::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) - { - internal::do_integrate_difference( - hp::StaticMappingQ1::mapping_collection, - dof, - fe_function, - exact_solution, - difference, - q, - norm, - weight, - exponent); - } - template double compute_global_error(const Triangulation &tria, diff --git a/source/numerics/vector_tools_boundary.inst.in b/source/numerics/vector_tools_boundary.inst.in index f2e24d69f3..1a7cfba47a 100644 --- a/source/numerics/vector_tools_boundary.inst.in +++ b/source/numerics/vector_tools_boundary.inst.in @@ -233,22 +233,6 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) # if deal_II_dimension != 1 template void - project_boundary_values_curl_conforming( - const DoFHandler &, - const unsigned int, - const Function &, - const types::boundary_id, - AffineConstraints &, - const Mapping &); - template void - project_boundary_values_curl_conforming( - const DoFHandler &, - const unsigned int, - const Function &, - const types::boundary_id, - AffineConstraints &, - const hp::MappingCollection &); - template void project_boundary_values_curl_conforming_l2( const DoFHandler &, const unsigned int, diff --git a/source/numerics/vector_tools_integrate_difference.inst.in b/source/numerics/vector_tools_integrate_difference.inst.in index 826fd7cdd9..1a5ac1ce5b 100644 --- a/source/numerics/vector_tools_integrate_difference.inst.in +++ b/source/numerics/vector_tools_integrate_difference.inst.in @@ -143,653 +143,6 @@ for (VEC : VECTOR_TYPES; deal_II_dimension : DIMENSIONS; #endif } -// All the functions in this for loop are deprecated -for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) - { -#if deal_II_dimension <= deal_II_space_dimension - namespace VectorTools - \{ - - template void - integrate_difference, - Vector, - deal_II_space_dimension>( - const Mapping &, - const DoFHandler &, - const Vector &, - const Function &, - Vector &, - const Quadrature &, - const NormType &, - const Function *, - const double); - - template void - integrate_difference, - Vector, - deal_II_space_dimension>( - const Mapping &, - const DoFHandler &, - const BlockVector &, - const Function &, - Vector &, - const Quadrature &, - const NormType &, - const Function *, - const double); - - template void - integrate_difference, - Vector, - deal_II_space_dimension>( - const Mapping &, - const DoFHandler &, - const LinearAlgebra::Vector &, - const Function &, - Vector &, - const Quadrature &, - const NormType &, - const Function *, - const double); - - template void - integrate_difference, - Vector, - deal_II_space_dimension>( - const Mapping &, - const DoFHandler &, - const LinearAlgebra::distributed::Vector &, - const Function &, - Vector &, - const Quadrature &, - const NormType &, - const Function *, - const double); - - template void - integrate_difference, - Vector, - deal_II_space_dimension>( - const Mapping &, - const DoFHandler &, - const LinearAlgebra::distributed::BlockVector &, - const Function &, - Vector &, - const Quadrature &, - const NormType &, - const Function *, - const double); - -# ifdef DEAL_II_TRILINOS_WITH_TPETRA - template void - integrate_difference, - Vector, - deal_II_space_dimension>( - const Mapping &, - const DoFHandler &, - const LinearAlgebra::TpetraWrappers::Vector &, - const Function &, - Vector &, - const Quadrature &, - const NormType &, - const Function *, - const double); -# endif - - template void - integrate_difference, - Vector, - deal_II_space_dimension>( - const DoFHandler &, - const Vector &, - const Function &, - Vector &, - const Quadrature &, - const NormType &, - const Function *, - const double); - - template void - integrate_difference, - Vector, - deal_II_space_dimension>( - const DoFHandler &, - const BlockVector &, - const Function &, - Vector &, - const Quadrature &, - const NormType &, - const Function *, - const double); - - template void - integrate_difference, - Vector, - deal_II_space_dimension>( - const DoFHandler &, - const LinearAlgebra::Vector &, - const Function &, - Vector &, - const Quadrature &, - const NormType &, - const Function *, - const double); - - template void - integrate_difference, - Vector, - deal_II_space_dimension>( - const DoFHandler &, - const LinearAlgebra::distributed::Vector &, - const Function &, - Vector &, - const Quadrature &, - const NormType &, - const Function *, - const double); - - template void - integrate_difference, - Vector, - deal_II_space_dimension>( - const DoFHandler &, - const LinearAlgebra::distributed::BlockVector &, - const Function &, - Vector &, - const Quadrature &, - const NormType &, - const Function *, - const double); - -# ifdef DEAL_II_TRILINOS_WITH_TPETRA - template void - integrate_difference, - Vector, - deal_II_space_dimension>( - const DoFHandler &, - const LinearAlgebra::TpetraWrappers::Vector &, - const Function &, - Vector &, - const Quadrature &, - const NormType &, - const Function *, - const double); -# endif - - template void - integrate_difference, - Vector, - deal_II_space_dimension>( - const Mapping &, - const DoFHandler &, - const BlockVector &, - const Function &, - Vector &, - const Quadrature &, - const NormType &, - const Function *, - const double); - - template void - integrate_difference, - Vector, - deal_II_space_dimension>( - const Mapping &, - const DoFHandler &, - const LinearAlgebra::Vector &, - const Function &, - Vector &, - const Quadrature &, - const NormType &, - const Function *, - const double); - - - template void - integrate_difference, - Vector, - deal_II_space_dimension>( - const Mapping &, - const DoFHandler &, - const LinearAlgebra::distributed::Vector &, - const Function &, - Vector &, - const Quadrature &, - const NormType &, - const Function *, - const double); - - template void - integrate_difference, - Vector, - deal_II_space_dimension>( - const Mapping &, - const DoFHandler &, - const LinearAlgebra::distributed::BlockVector &, - const Function &, - Vector &, - const Quadrature &, - const NormType &, - const Function *, - const double); - -# ifdef DEAL_II_TRILINOS_WITH_TPETRA - template void - integrate_difference, - Vector, - deal_II_space_dimension>( - const Mapping &, - const DoFHandler &, - const LinearAlgebra::TpetraWrappers::Vector &, - const Function &, - Vector &, - const Quadrature &, - const NormType &, - const Function *, - const double); -# endif - - template void - integrate_difference, - Vector, - deal_II_space_dimension>( - const DoFHandler &, - const BlockVector &, - const Function &, - Vector &, - const Quadrature &, - const NormType &, - const Function *, - const double); - - template void - integrate_difference, - Vector, - deal_II_space_dimension>( - const DoFHandler &, - const LinearAlgebra::Vector &, - const Function &, - Vector &, - const Quadrature &, - const NormType &, - const Function *, - const double); - - template void - integrate_difference, - Vector, - deal_II_space_dimension>( - const DoFHandler &, - const LinearAlgebra::distributed::Vector &, - const Function &, - Vector &, - const Quadrature &, - const NormType &, - const Function *, - const double); - - template void - integrate_difference, - Vector, - deal_II_space_dimension>( - const DoFHandler &, - const LinearAlgebra::distributed::BlockVector &, - const Function &, - Vector &, - const Quadrature &, - const NormType &, - const Function *, - const double); - -# ifdef DEAL_II_TRILINOS_WITH_TPETRA - template void - integrate_difference, - Vector, - deal_II_space_dimension>( - const DoFHandler &, - const LinearAlgebra::TpetraWrappers::Vector &, - const Function &, - Vector &, - const Quadrature &, - const NormType &, - const Function *, - const double); -# endif - - template void - integrate_difference, - Vector, - deal_II_space_dimension>( - const hp::MappingCollection - &, - const DoFHandler &, - const BlockVector &, - const Function &, - Vector &, - const hp::QCollection &, - const NormType &, - const Function *, - const double); - - template void - integrate_difference, - Vector, - deal_II_space_dimension>( - const hp::MappingCollection - &, - const DoFHandler &, - const LinearAlgebra::Vector &, - const Function &, - Vector &, - const hp::QCollection &, - const NormType &, - const Function *, - const double); - - template void - integrate_difference, - Vector, - deal_II_space_dimension>( - const hp::MappingCollection - &, - const DoFHandler &, - const LinearAlgebra::distributed::Vector &, - const Function &, - Vector &, - const hp::QCollection &, - const NormType &, - const Function *, - const double); - - template void - integrate_difference, - Vector, - deal_II_space_dimension>( - const hp::MappingCollection - &, - const DoFHandler &, - const LinearAlgebra::distributed::BlockVector &, - const Function &, - Vector &, - const hp::QCollection &, - const NormType &, - const Function *, - const double); - -# ifdef DEAL_II_TRILINOS_WITH_TPETRA - template void - integrate_difference, - Vector, - deal_II_space_dimension>( - const hp::MappingCollection - &, - const DoFHandler &, - const LinearAlgebra::TpetraWrappers::Vector &, - const Function &, - Vector &, - const hp::QCollection &, - const NormType &, - const Function *, - const double); -# endif - - template void - integrate_difference, - Vector, - deal_II_space_dimension>( - const DoFHandler &, - const BlockVector &, - const Function &, - Vector &, - const hp::QCollection &, - const NormType &, - const Function *, - const double); - - template void - integrate_difference, - Vector, - deal_II_space_dimension>( - const DoFHandler &, - const LinearAlgebra::Vector &, - const Function &, - Vector &, - const hp::QCollection &, - const NormType &, - const Function *, - const double); - - template void - integrate_difference, - Vector, - deal_II_space_dimension>( - const DoFHandler &, - const LinearAlgebra::distributed::Vector &, - const Function &, - Vector &, - const hp::QCollection &, - const NormType &, - const Function *, - const double); - - template void - integrate_difference, - Vector, - deal_II_space_dimension>( - const DoFHandler &, - const LinearAlgebra::distributed::BlockVector &, - const Function &, - Vector &, - const hp::QCollection &, - const NormType &, - const Function *, - const double); - -# ifdef DEAL_II_TRILINOS_WITH_TPETRA - template void - integrate_difference, - Vector, - deal_II_space_dimension>( - const DoFHandler &, - const LinearAlgebra::TpetraWrappers::Vector &, - const Function &, - Vector &, - const hp::QCollection &, - const NormType &, - const Function *, - const double); -# endif - - template void - integrate_difference, - Vector, - deal_II_space_dimension>( - const hp::MappingCollection - &, - const DoFHandler &, - const BlockVector &, - const Function &, - Vector &, - const hp::QCollection &, - const NormType &, - const Function *, - const double); - - template void - integrate_difference, - Vector, - deal_II_space_dimension>( - const hp::MappingCollection - &, - const DoFHandler &, - const LinearAlgebra::Vector &, - const Function &, - Vector &, - const hp::QCollection &, - const NormType &, - const Function *, - const double); - - template void - integrate_difference, - Vector, - deal_II_space_dimension>( - const hp::MappingCollection - &, - const DoFHandler &, - const LinearAlgebra::distributed::Vector &, - const Function &, - Vector &, - const hp::QCollection &, - const NormType &, - const Function *, - const double); - - template void - integrate_difference, - Vector, - deal_II_space_dimension>( - const hp::MappingCollection - &, - const DoFHandler &, - const LinearAlgebra::distributed::BlockVector &, - const Function &, - Vector &, - const hp::QCollection &, - const NormType &, - const Function *, - const double); - -# ifdef DEAL_II_TRILINOS_WITH_TPETRA - template void - integrate_difference, - Vector, - deal_II_space_dimension>( - const hp::MappingCollection - &, - const DoFHandler &, - const LinearAlgebra::TpetraWrappers::Vector &, - const Function &, - Vector &, - const hp::QCollection &, - const NormType &, - const Function *, - const double); -# endif - - template void - integrate_difference, - Vector, - deal_II_space_dimension>( - const DoFHandler &, - const BlockVector &, - const Function &, - Vector &, - const hp::QCollection &, - const NormType &, - const Function *, - const double); - - template void - integrate_difference, - Vector, - deal_II_space_dimension>( - const DoFHandler &, - const LinearAlgebra::Vector &, - const Function &, - Vector &, - const hp::QCollection &, - const NormType &, - const Function *, - const double); - - template void - integrate_difference, - Vector, - deal_II_space_dimension>( - const DoFHandler &, - const LinearAlgebra::distributed::Vector &, - const Function &, - Vector &, - const hp::QCollection &, - const NormType &, - const Function *, - const double); - - template void - integrate_difference, - Vector, - deal_II_space_dimension>( - const DoFHandler &, - const LinearAlgebra::distributed::BlockVector &, - const Function &, - Vector &, - const hp::QCollection &, - const NormType &, - const Function *, - const double); - -# ifdef DEAL_II_TRILINOS_WITH_TPETRA - template void - integrate_difference, - Vector, - deal_II_space_dimension>( - const DoFHandler &, - const LinearAlgebra::TpetraWrappers::Vector &, - const Function &, - Vector &, - const hp::QCollection &, - const NormType &, - const Function *, - const double); -# endif - - \} -#endif - } - for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) { #if deal_II_dimension <= deal_II_space_dimension