--- /dev/null
+Removed: The deprecated specializations of VectorTools::integrate_difference()
+and VectorTools::project_boundary_values_curl_conforming() have been removed.
+<br>
+(Daniel Arndt, 2021/06/04)
* 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.
AffineConstraints<number> &constraints,
std::vector<unsigned int> 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 <code>dim</code> 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.
- *
- * <h4>Computing constraints</h4>
- *
- * 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 <int dim>
- DEAL_II_DEPRECATED void
- project_boundary_values_curl_conforming(
- const DoFHandler<dim, dim> & dof_handler,
- const unsigned int first_vector_component,
- const Function<dim, double> &boundary_function,
- const types::boundary_id boundary_component,
- AffineConstraints<double> & constraints,
- const Mapping<dim> & 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 <int dim>
- DEAL_II_DEPRECATED void
- project_boundary_values_curl_conforming(
- const DoFHandler<dim, dim> & dof_handler,
- const unsigned int first_vector_component,
- const Function<dim, double> & boundary_function,
- const types::boundary_id boundary_component,
- AffineConstraints<double> & constraints,
- const hp::MappingCollection<dim, dim> &mapping_collection =
- hp::StaticMappingQ1<dim>::mapping_collection);
-
/**
* This function is an updated version of the
* project_boundary_values_curl_conforming function. The intention is to fix
}
- namespace internals
- {
- // This function computes the
- // projection of the boundary
- // function on edges for 3D.
- template <typename cell_iterator>
- 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<double> &dof_values,
- std::vector<bool> & 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<dim>::lines_per_face +
- line);
-
- // Initialize the required
- // objects.
- const FEValues<dim> &fe_values = hp_fe_values.get_present_fe_values();
- const FiniteElement<dim> & fe = cell->get_fe();
- const std::vector<DerivativeForm<1, dim, spacedim>> &jacobians =
- fe_values.get_jacobians();
- const std::vector<Point<dim>> &quadrature_points =
- fe_values.get_quadrature_points();
-
- std::vector<Tensor<1, dim>> tangentials(fe_values.n_quadrature_points);
- std::vector<Vector<double>> values(fe_values.n_quadrature_points,
- Vector<double>(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<Point<dim>> &reference_quadrature_points =
- fe_values.get_quadrature().get_points();
- std::pair<unsigned int, unsigned int> base_indices(0, 0);
-
- if (dynamic_cast<const FESystem<dim> *>(&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<dim>::faces_per_cell]
- [GeometryInfo<dim>::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<dim> shifted_reference_point_1 =
- reference_quadrature_points[q_point];
- Point<dim> 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<const FESystem<dim> *>(&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<const FE_Nedelec<dim> *>(&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 <int dim, typename cell_iterator>
- void
- compute_edge_projection(const cell_iterator &,
- const unsigned int,
- const unsigned int,
- hp::FEValues<dim> &,
- const Function<dim> &,
- const unsigned int,
- std::vector<double> &,
- std::vector<bool> &)
- {
- Assert(false, ExcInternalError());
- }
-
- // This function computes the
- // projection of the boundary
- // function on the interior of
- // faces.
- template <int dim, typename cell_iterator, typename number>
- void
- compute_face_projection_curl_conforming(
- const cell_iterator & cell,
- const unsigned int face,
- hp::FEValues<dim> & hp_fe_values,
- const Function<dim, number> &boundary_function,
- const unsigned int first_vector_component,
- std::vector<double> & dof_values,
- std::vector<bool> & 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<dim> &fe_values = hp_fe_values.get_present_fe_values();
- const FiniteElement<dim> & fe = cell->get_fe();
- const std::vector<DerivativeForm<1, dim, spacedim>> &jacobians =
- fe_values.get_jacobians();
- const std::vector<Point<dim>> &quadrature_points =
- fe_values.get_quadrature_points();
- const unsigned int degree = fe.degree - 1;
- std::pair<unsigned int, unsigned int> base_indices(0, 0);
-
- if (dynamic_cast<const FESystem<dim> *>(&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<Vector<double>> values(fe_values.n_quadrature_points,
- Vector<double>(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<Tensor<1, dim>> tangentials(
- fe_values.n_quadrature_points);
-
- const std::vector<Point<dim>> &reference_quadrature_points =
- fe_values.get_quadrature().get_points();
-
- // coordinate directions
- // of the face.
- const unsigned int
- face_coordinate_direction[GeometryInfo<dim>::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<dim> shifted_reference_point_1 =
- reference_quadrature_points[q_point];
- Point<dim> 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<const FESystem<dim> *>(&fe) !=
- nullptr) &&
- (fe.system_to_base_index(
- fe.face_to_cell_index(i, face))
- .first == base_indices)) ||
- (dynamic_cast<const FE_Nedelec<dim> *>(&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<double> assembling_matrix(
- degree * fe.degree, dim * fe_values.n_quadrature_points);
- Vector<double> assembling_vector(assembling_matrix.n());
- Vector<double> cell_rhs(assembling_matrix.m());
- FullMatrix<double> cell_matrix(assembling_matrix.m(),
- assembling_matrix.m());
- FullMatrix<double> cell_matrix_inv(assembling_matrix.m(),
- assembling_matrix.m());
- Vector<double> 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<const FESystem<dim> *>(&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<const FE_Nedelec<dim> *>(&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<const FESystem<dim> *>(&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<dim>::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<dim>::lines_per_face) *
- fe.degree,
- face))) ||
- ((dynamic_cast<const FE_Nedelec<dim> *>(&fe) !=
- nullptr) &&
- (GeometryInfo<dim>::lines_per_face * fe.degree <= i) &&
- (i < (degree + GeometryInfo<dim>::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<const FESystem<dim> *>(&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<dim>::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<dim>::lines_per_face) *
- fe.degree,
- face))) ||
- ((dynamic_cast<const FE_Nedelec<dim> *>(&fe) !=
- nullptr) &&
- (GeometryInfo<dim>::lines_per_face * fe.degree <= i) &&
- (i < (degree + GeometryInfo<dim>::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<const FESystem<dim> *>(&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<const FE_Nedelec<dim> *>(&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<const FESystem<dim> *>(&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<dim>::lines_per_face) *
- fe.degree,
- face) <=
- fe.system_to_base_index(
- fe.face_to_cell_index(i, face))
- .second)) ||
- ((dynamic_cast<const FE_Nedelec<dim> *>(&fe) !=
- nullptr) &&
- ((degree + GeometryInfo<dim>::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<const FESystem<dim> *>(&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<dim>::lines_per_face) *
- fe.degree,
- face) <=
- fe.system_to_base_index(fe.face_to_cell_index(i, face))
- .second)) ||
- ((dynamic_cast<const FE_Nedelec<dim> *>(&fe) != nullptr) &&
- ((degree + GeometryInfo<dim>::lines_per_face) *
- fe.degree <=
- i)))
- {
- dof_values[i] = solution(index);
- dofs_processed[i] = true;
- ++index;
- }
-
- break;
- }
-
- default:
- Assert(false, ExcNotImplemented());
- }
- }
- } // namespace internals
-
-
-
- template <int dim>
- void
- project_boundary_values_curl_conforming(
- const DoFHandler<dim> & dof_handler,
- const unsigned int first_vector_component,
- const Function<dim> & boundary_function,
- const types::boundary_id boundary_component,
- AffineConstraints<double> &constraints,
- const Mapping<dim> & 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
- // <tt>degree * (degree + 1)<tt> to obtain the values for the
- // corresponding degrees of freedom.
- const unsigned int superdegree = dof_handler.get_fe().degree;
- const QGauss<dim - 1> 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<dim> &fe_collection(dof_handler.get_fe_collection());
- const hp::MappingCollection<dim> mapping_collection(mapping);
- hp::QCollection<dim> face_quadrature_collection;
-
- for (unsigned int face : GeometryInfo<dim>::face_indices())
- face_quadrature_collection.push_back(
- QProjector<dim>::project_to_face(reference_face_quadrature, face));
-
- hp::FEValues<dim> fe_face_values(mapping_collection,
- fe_collection,
- face_quadrature_collection,
- update_jacobians | update_JxW_values |
- update_quadrature_points |
- update_values);
-
- std::vector<bool> dofs_processed(dofs_per_face);
- std::vector<double> dof_values(dofs_per_face);
- std::vector<types::global_dof_index> face_dof_indices(dofs_per_face);
- typename DoFHandler<dim>::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<const FE_Nothing<dim> *>(
- &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<const FESystem<dim> *>(
- &cell->get_fe()) == nullptr)
- {
- AssertThrow(
- dynamic_cast<const FE_Nedelec<dim> *>(
- &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<dim - 2> reference_edge_quadrature(2 * superdegree);
- const unsigned int degree = superdegree - 1;
- hp::QCollection<dim> edge_quadrature_collection;
-
- for (const unsigned int face : GeometryInfo<dim>::face_indices())
- for (unsigned int line = 0;
- line < GeometryInfo<dim>::lines_per_face;
- ++line)
- edge_quadrature_collection.push_back(
- QProjector<dim>::project_to_face(
- QProjector<dim - 1>::project_to_face(
- reference_edge_quadrature, line),
- face));
-
- hp::FEValues<dim> 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<const FE_Nothing<dim> *>(
- &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<const FESystem<dim> *>(
- &cell->get_fe()) == nullptr)
- {
- AssertThrow(dynamic_cast<const FE_Nedelec<dim> *>(
- &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 <int dim>
- void
-
- project_boundary_values_curl_conforming(
- const DoFHandler<dim> & dof_handler,
- const unsigned int first_vector_component,
- const Function<dim> & boundary_function,
- const types::boundary_id boundary_component,
- AffineConstraints<double> & constraints,
- const hp::MappingCollection<dim> &mapping_collection)
- {
- const hp::FECollection<dim> &fe_collection(dof_handler.get_fe_collection());
- hp::QCollection<dim> face_quadrature_collection;
-
- for (unsigned int i = 0; i < fe_collection.size(); ++i)
- {
- const QGauss<dim - 1> reference_face_quadrature(
- 2 * fe_collection[i].degree);
-
- for (unsigned int face : GeometryInfo<dim>::face_indices())
- face_quadrature_collection.push_back(
- QProjector<dim>::project_to_face(reference_face_quadrature, face));
- }
-
- hp::FEValues<dim> fe_face_values(mapping_collection,
- fe_collection,
- face_quadrature_collection,
- update_jacobians | update_JxW_values |
- update_quadrature_points |
- update_values);
- std::vector<bool> dofs_processed;
- std::vector<double> dof_values;
- std::vector<types::global_dof_index> face_dof_indices;
- typename DoFHandler<dim>::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<const FE_Nothing<dim> *>(
- &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<const FESystem<dim> *>(
- &cell->get_fe()) == nullptr)
- {
- AssertThrow(dynamic_cast<const FE_Nedelec<dim> *>(
- &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<dim> edge_quadrature_collection;
-
- for (unsigned int i = 0; i < fe_collection.size(); ++i)
- {
- const QGauss<dim - 2> reference_edge_quadrature(
- 2 * fe_collection[i].degree);
-
- for (const unsigned int face :
- GeometryInfo<dim>::face_indices())
- for (unsigned int line = 0;
- line < GeometryInfo<dim>::lines_per_face;
- ++line)
- edge_quadrature_collection.push_back(
- QProjector<dim>::project_to_face(
- QProjector<dim - 1>::project_to_face(
- reference_edge_quadrature, line),
- face));
- }
-
- hp::FEValues<dim> 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<const FE_Nothing<dim> *>(
- &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<const FESystem<dim> *>(
- &cell->get_fe()) == nullptr)
- {
- AssertThrow(dynamic_cast<const FE_Nedelec<dim> *>(
- &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<dim>::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 <int dim, typename cell_iterator, typename number>
const Function<spacedim, double> * 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<dim, spacedim> &, const DoFHandler<dim, spacedim> &, const InVector &, const Function<spacedim, typename InVector::value_type> &, OutVector &, const Quadrature<dim> &, const NormType &, const Function<spacedim, double> *, const double) instead.
- */
- template <int dim, class InVector, class OutVector, int spacedim>
- DEAL_II_DEPRECATED typename std::enable_if<
- !std::is_same<typename InVector::value_type, double>::value>::type
- integrate_difference(const Mapping<dim, spacedim> & mapping,
- const DoFHandler<dim, spacedim> & dof,
- const InVector & fe_function,
- const Function<spacedim, double> &exact_solution,
- OutVector & difference,
- const Quadrature<dim> & q,
- const NormType & norm,
- const Function<spacedim, double> *weight = nullptr,
- const double exponent = 2.);
-
- /**
- * Call the integrate_difference() function, see above, with
- * <tt>mapping=MappingQGeneric@<dim@>(1)</tt>.
- *
- * @deprecated Use integrate_difference(const DoFHandler<dim, spacedim> &, const InVector &, const Function<spacedim, typename InVector::value_type> &exact_solution, OutVector &, const Quadrature<dim> &, const NormType &, const Function<spacedim, double> *, const double) instead.
- */
- template <int dim, class InVector, class OutVector, int spacedim>
- DEAL_II_DEPRECATED typename std::enable_if<
- !std::is_same<typename InVector::value_type, double>::value>::type
- integrate_difference(const DoFHandler<dim, spacedim> & dof,
- const InVector & fe_function,
- const Function<spacedim, double> &exact_solution,
- OutVector & difference,
- const Quadrature<dim> & q,
- const NormType & norm,
- const Function<spacedim, double> *weight = nullptr,
- const double exponent = 2.);
-
- /**
- * Same as above for hp.
- *
- * @deprecated Use integrate_difference(const hp::MappingCollection<dim, spacedim> &, const DoFHandler<dim, spacedim> &, const InVector &, const Function<spacedim, typename InVector::value_type> &, OutVector &, const hp::QCollection<dim> &, const NormType &, const Function<spacedim, double> *, const double) instead.
- */
- template <int dim, class InVector, class OutVector, int spacedim>
- DEAL_II_DEPRECATED typename std::enable_if<
- !std::is_same<typename InVector::value_type, double>::value>::type
- integrate_difference(const hp::MappingCollection<dim, spacedim> &mapping,
- const DoFHandler<dim, spacedim> & dof,
- const InVector & fe_function,
- const Function<spacedim, double> &exact_solution,
- OutVector & difference,
- const hp::QCollection<dim> & q,
- const NormType & norm,
- const Function<spacedim, double> *weight = nullptr,
- const double exponent = 2.);
-
- /**
- * Call the integrate_difference() function, see above, with
- * <tt>mapping=MappingQGeneric@<dim@>(1)</tt>.
- *
- * @deprecated Use integrate_difference(const DoFHandler<dim, spacedim> &, const InVector &, const Function<spacedim, typename InVector::value_type> &, OutVector &, const hp::QCollection<dim> &, const NormType &, const Function<spacedim, double> *, const double) instead.
- */
- template <int dim, class InVector, class OutVector, int spacedim>
- DEAL_II_DEPRECATED typename std::enable_if<
- !std::is_same<typename InVector::value_type, double>::value>::type
- integrate_difference(const DoFHandler<dim, spacedim> & dof,
- const InVector & fe_function,
- const Function<spacedim, double> &exact_solution,
- OutVector & difference,
- const hp::QCollection<dim> & q,
- const NormType & norm,
- const Function<spacedim, double> *weight = nullptr,
- const double exponent = 2.);
-
/**
* Take a Vector @p cellwise_error of errors on each cell with
* <tt>tria.n_active_cells()</tt> entries and return the global
n_q_points, std::vector<Tensor<1, spacedim, Number>>(n_components));
}
- template <int dim, int spacedim, typename Number>
- struct DEAL_II_DEPRECATED DeprecatedIDScratchData
- {
- DeprecatedIDScratchData(
- const dealii::hp::MappingCollection<dim, spacedim> &mapping,
- const dealii::hp::FECollection<dim, spacedim> & fe,
- const dealii::hp::QCollection<dim> & 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<Vector<Number>> function_values;
- std::vector<std::vector<Tensor<1, spacedim, Number>>> function_grads;
- std::vector<double> weight_values;
- std::vector<Vector<double>> weight_vectors;
-
- std::vector<Vector<Number>> psi_values;
- std::vector<std::vector<Tensor<1, spacedim, Number>>> psi_grads;
- std::vector<Number> psi_scalar;
-
- std::vector<double> tmp_values;
- std::vector<Vector<double>> tmp_vector_values;
- std::vector<Tensor<1, spacedim>> tmp_gradients;
- std::vector<std::vector<Tensor<1, spacedim>>> tmp_vector_gradients;
-
- dealii::hp::FEValues<dim, spacedim> x_fe_values;
- };
-
-
- template <int dim, int spacedim, typename Number>
- DeprecatedIDScratchData<dim, spacedim, Number>::DeprecatedIDScratchData(
- const dealii::hp::MappingCollection<dim, spacedim> &mapping,
- const dealii::hp::FECollection<dim, spacedim> & fe,
- const dealii::hp::QCollection<dim> & q,
- const UpdateFlags update_flags)
- : x_fe_values(mapping, fe, q, update_flags)
- {}
-
- template <int dim, int spacedim, typename Number>
- DeprecatedIDScratchData<dim, spacedim, Number>::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 <int dim, int spacedim, typename Number>
- void
- DeprecatedIDScratchData<dim, spacedim, Number>::resize_vectors(
- const unsigned int n_q_points,
- const unsigned int n_components)
- {
- function_values.resize(n_q_points, Vector<Number>(n_components));
- function_grads.resize(
- n_q_points, std::vector<Tensor<1, spacedim, Number>>(n_components));
-
- weight_values.resize(n_q_points);
- weight_vectors.resize(n_q_points, Vector<double>(n_components));
-
- psi_values.resize(n_q_points, Vector<Number>(n_components));
- psi_grads.resize(n_q_points,
- std::vector<Tensor<1, spacedim, Number>>(n_components));
- psi_scalar.resize(n_q_points);
-
- tmp_values.resize(n_q_points);
- tmp_vector_values.resize(n_q_points, Vector<double>(n_components));
- tmp_gradients.resize(n_q_points);
- tmp_vector_gradients.resize(
- n_q_points, std::vector<Tensor<1, spacedim>>(n_components));
- }
-
namespace internal
{
template <typename number>
return diff;
}
- template <int dim, int spacedim, typename Number>
- DEAL_II_DEPRECATED
- typename std::enable_if<!std::is_same<Number, double>::value,
- double>::type
- integrate_difference_inner(
- const Function<spacedim> & exact_solution,
- const NormType & norm,
- const Function<spacedim> * weight,
- const UpdateFlags update_flags,
- const double exponent,
- const unsigned int n_components,
- DeprecatedIDScratchData<dim, spacedim, Number> &data)
- {
- const bool fe_is_system = (n_components != 1);
- const dealii::FEValues<dim, spacedim> &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<dim,double>
- // 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<Number, double>::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<double>(
- numbers::NumberTraits<Number>::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<Number>::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<Number>::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 <int dim, int spacedim, class InVector, class OutVector>
difference(cell->active_cell_index()) = 0;
}
- template <int dim, int spacedim, class InVector, class OutVector>
- DEAL_II_DEPRECATED static typename std::enable_if<
- !std::is_same<typename InVector::value_type, double>::value>::type
- do_integrate_difference(
- const dealii::hp::MappingCollection<dim, spacedim> &mapping,
- const DoFHandler<dim, spacedim> & dof,
- const InVector & fe_function,
- const Function<spacedim> & exact_solution,
- OutVector & difference,
- const dealii::hp::QCollection<dim> & q,
- const NormType & norm,
- const Function<spacedim> * 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<dim, spacedim> &fe_collection =
- dof.get_fe_collection();
- DeprecatedIDScratchData<dim, spacedim, Number> 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<dim, spacedim> &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<dim, spacedim, Number>(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 <int dim, class InVector, class OutVector, int spacedim>
exponent);
}
- template <int dim, class InVector, class OutVector, int spacedim>
- DEAL_II_DEPRECATED typename std::enable_if<
- !std::is_same<typename InVector::value_type, double>::value>::type
- integrate_difference(const Mapping<dim, spacedim> & mapping,
- const DoFHandler<dim, spacedim> &dof,
- const InVector & fe_function,
- const Function<spacedim> & exact_solution,
- OutVector & difference,
- const Quadrature<dim> & q,
- const NormType & norm,
- const Function<spacedim> * weight,
- const double exponent)
- {
- internal::do_integrate_difference(hp::MappingCollection<dim, spacedim>(
- mapping),
- dof,
- fe_function,
- exact_solution,
- difference,
- hp::QCollection<dim>(q),
- norm,
- weight,
- exponent);
- }
-
template <int dim, class InVector, class OutVector, int spacedim>
void
}
- template <int dim, class InVector, class OutVector, int spacedim>
- DEAL_II_DEPRECATED typename std::enable_if<
- !std::is_same<typename InVector::value_type, double>::value>::type
- integrate_difference(const DoFHandler<dim, spacedim> &dof,
- const InVector & fe_function,
- const Function<spacedim> & exact_solution,
- OutVector & difference,
- const Quadrature<dim> & q,
- const NormType & norm,
- const Function<spacedim> * weight,
- const double exponent)
- {
- internal::do_integrate_difference(
- hp::StaticMappingQ1<dim, spacedim>::mapping_collection,
- dof,
- fe_function,
- exact_solution,
- difference,
- hp::QCollection<dim>(q),
- norm,
- weight,
- exponent);
- }
-
-
-
template <int dim, class InVector, class OutVector, int spacedim>
void
integrate_difference(
exponent);
}
- template <int dim, class InVector, class OutVector, int spacedim>
- DEAL_II_DEPRECATED typename std::enable_if<
- !std::is_same<typename InVector::value_type, double>::value>::type
- integrate_difference(
- const dealii::hp::MappingCollection<dim, spacedim> &mapping,
- const dealii::DoFHandler<dim, spacedim> & dof,
- const InVector & fe_function,
- const Function<spacedim> & exact_solution,
- OutVector & difference,
- const dealii::hp::QCollection<dim> & q,
- const NormType & norm,
- const Function<spacedim> * weight,
- const double exponent)
- {
- internal::do_integrate_difference(hp::MappingCollection<dim, spacedim>(
- mapping),
- dof,
- fe_function,
- exact_solution,
- difference,
- q,
- norm,
- weight,
- exponent);
- }
-
-
template <int dim, class InVector, class OutVector, int spacedim>
void
integrate_difference(
exponent);
}
- template <int dim, class InVector, class OutVector, int spacedim>
- DEAL_II_DEPRECATED typename std::enable_if<
- !std::is_same<typename InVector::value_type, double>::value>::type
- integrate_difference(const dealii::DoFHandler<dim, spacedim> &dof,
- const InVector & fe_function,
- const Function<spacedim> & exact_solution,
- OutVector & difference,
- const dealii::hp::QCollection<dim> & q,
- const NormType & norm,
- const Function<spacedim> * weight,
- const double exponent)
- {
- internal::do_integrate_difference(
- hp::StaticMappingQ1<dim, spacedim>::mapping_collection,
- dof,
- fe_function,
- exact_solution,
- difference,
- q,
- norm,
- weight,
- exponent);
- }
-
template <int dim, int spacedim, class InVector>
double
compute_global_error(const Triangulation<dim, spacedim> &tria,
# if deal_II_dimension != 1
template void
- project_boundary_values_curl_conforming<deal_II_dimension>(
- const DoFHandler<deal_II_dimension> &,
- const unsigned int,
- const Function<deal_II_dimension> &,
- const types::boundary_id,
- AffineConstraints<double> &,
- const Mapping<deal_II_dimension> &);
- template void
- project_boundary_values_curl_conforming<deal_II_dimension>(
- const DoFHandler<deal_II_dimension> &,
- const unsigned int,
- const Function<deal_II_dimension> &,
- const types::boundary_id,
- AffineConstraints<double> &,
- const hp::MappingCollection<deal_II_dimension> &);
- template void
project_boundary_values_curl_conforming_l2(
const DoFHandler<deal_II_dimension> &,
const unsigned int,
#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<deal_II_dimension,
- Vector<float>,
- Vector<float>,
- deal_II_space_dimension>(
- const Mapping<deal_II_dimension, deal_II_space_dimension> &,
- const DoFHandler<deal_II_dimension, deal_II_space_dimension> &,
- const Vector<float> &,
- const Function<deal_II_space_dimension> &,
- Vector<float> &,
- const Quadrature<deal_II_dimension> &,
- const NormType &,
- const Function<deal_II_space_dimension> *,
- const double);
-
- template void
- integrate_difference<deal_II_dimension,
- BlockVector<float>,
- Vector<float>,
- deal_II_space_dimension>(
- const Mapping<deal_II_dimension, deal_II_space_dimension> &,
- const DoFHandler<deal_II_dimension, deal_II_space_dimension> &,
- const BlockVector<float> &,
- const Function<deal_II_space_dimension> &,
- Vector<float> &,
- const Quadrature<deal_II_dimension> &,
- const NormType &,
- const Function<deal_II_space_dimension> *,
- const double);
-
- template void
- integrate_difference<deal_II_dimension,
- LinearAlgebra::Vector<float>,
- Vector<float>,
- deal_II_space_dimension>(
- const Mapping<deal_II_dimension, deal_II_space_dimension> &,
- const DoFHandler<deal_II_dimension, deal_II_space_dimension> &,
- const LinearAlgebra::Vector<float> &,
- const Function<deal_II_space_dimension> &,
- Vector<float> &,
- const Quadrature<deal_II_dimension> &,
- const NormType &,
- const Function<deal_II_space_dimension> *,
- const double);
-
- template void
- integrate_difference<deal_II_dimension,
- LinearAlgebra::distributed::Vector<float>,
- Vector<float>,
- deal_II_space_dimension>(
- const Mapping<deal_II_dimension, deal_II_space_dimension> &,
- const DoFHandler<deal_II_dimension, deal_II_space_dimension> &,
- const LinearAlgebra::distributed::Vector<float> &,
- const Function<deal_II_space_dimension> &,
- Vector<float> &,
- const Quadrature<deal_II_dimension> &,
- const NormType &,
- const Function<deal_II_space_dimension> *,
- const double);
-
- template void
- integrate_difference<deal_II_dimension,
- LinearAlgebra::distributed::BlockVector<float>,
- Vector<float>,
- deal_II_space_dimension>(
- const Mapping<deal_II_dimension, deal_II_space_dimension> &,
- const DoFHandler<deal_II_dimension, deal_II_space_dimension> &,
- const LinearAlgebra::distributed::BlockVector<float> &,
- const Function<deal_II_space_dimension> &,
- Vector<float> &,
- const Quadrature<deal_II_dimension> &,
- const NormType &,
- const Function<deal_II_space_dimension> *,
- const double);
-
-# ifdef DEAL_II_TRILINOS_WITH_TPETRA
- template void
- integrate_difference<deal_II_dimension,
- LinearAlgebra::TpetraWrappers::Vector<float>,
- Vector<float>,
- deal_II_space_dimension>(
- const Mapping<deal_II_dimension, deal_II_space_dimension> &,
- const DoFHandler<deal_II_dimension, deal_II_space_dimension> &,
- const LinearAlgebra::TpetraWrappers::Vector<float> &,
- const Function<deal_II_space_dimension> &,
- Vector<float> &,
- const Quadrature<deal_II_dimension> &,
- const NormType &,
- const Function<deal_II_space_dimension> *,
- const double);
-# endif
-
- template void
- integrate_difference<deal_II_dimension,
- Vector<float>,
- Vector<float>,
- deal_II_space_dimension>(
- const DoFHandler<deal_II_dimension, deal_II_space_dimension> &,
- const Vector<float> &,
- const Function<deal_II_space_dimension> &,
- Vector<float> &,
- const Quadrature<deal_II_dimension> &,
- const NormType &,
- const Function<deal_II_space_dimension> *,
- const double);
-
- template void
- integrate_difference<deal_II_dimension,
- BlockVector<float>,
- Vector<float>,
- deal_II_space_dimension>(
- const DoFHandler<deal_II_dimension, deal_II_space_dimension> &,
- const BlockVector<float> &,
- const Function<deal_II_space_dimension> &,
- Vector<float> &,
- const Quadrature<deal_II_dimension> &,
- const NormType &,
- const Function<deal_II_space_dimension> *,
- const double);
-
- template void
- integrate_difference<deal_II_dimension,
- LinearAlgebra::Vector<float>,
- Vector<float>,
- deal_II_space_dimension>(
- const DoFHandler<deal_II_dimension, deal_II_space_dimension> &,
- const LinearAlgebra::Vector<float> &,
- const Function<deal_II_space_dimension> &,
- Vector<float> &,
- const Quadrature<deal_II_dimension> &,
- const NormType &,
- const Function<deal_II_space_dimension> *,
- const double);
-
- template void
- integrate_difference<deal_II_dimension,
- LinearAlgebra::distributed::Vector<float>,
- Vector<float>,
- deal_II_space_dimension>(
- const DoFHandler<deal_II_dimension, deal_II_space_dimension> &,
- const LinearAlgebra::distributed::Vector<float> &,
- const Function<deal_II_space_dimension> &,
- Vector<float> &,
- const Quadrature<deal_II_dimension> &,
- const NormType &,
- const Function<deal_II_space_dimension> *,
- const double);
-
- template void
- integrate_difference<deal_II_dimension,
- LinearAlgebra::distributed::BlockVector<float>,
- Vector<float>,
- deal_II_space_dimension>(
- const DoFHandler<deal_II_dimension, deal_II_space_dimension> &,
- const LinearAlgebra::distributed::BlockVector<float> &,
- const Function<deal_II_space_dimension> &,
- Vector<float> &,
- const Quadrature<deal_II_dimension> &,
- const NormType &,
- const Function<deal_II_space_dimension> *,
- const double);
-
-# ifdef DEAL_II_TRILINOS_WITH_TPETRA
- template void
- integrate_difference<deal_II_dimension,
- LinearAlgebra::TpetraWrappers::Vector<float>,
- Vector<float>,
- deal_II_space_dimension>(
- const DoFHandler<deal_II_dimension, deal_II_space_dimension> &,
- const LinearAlgebra::TpetraWrappers::Vector<float> &,
- const Function<deal_II_space_dimension> &,
- Vector<float> &,
- const Quadrature<deal_II_dimension> &,
- const NormType &,
- const Function<deal_II_space_dimension> *,
- const double);
-# endif
-
- template void
- integrate_difference<deal_II_dimension,
- BlockVector<float>,
- Vector<double>,
- deal_II_space_dimension>(
- const Mapping<deal_II_dimension, deal_II_space_dimension> &,
- const DoFHandler<deal_II_dimension, deal_II_space_dimension> &,
- const BlockVector<float> &,
- const Function<deal_II_space_dimension> &,
- Vector<double> &,
- const Quadrature<deal_II_dimension> &,
- const NormType &,
- const Function<deal_II_space_dimension> *,
- const double);
-
- template void
- integrate_difference<deal_II_dimension,
- LinearAlgebra::Vector<float>,
- Vector<double>,
- deal_II_space_dimension>(
- const Mapping<deal_II_dimension, deal_II_space_dimension> &,
- const DoFHandler<deal_II_dimension, deal_II_space_dimension> &,
- const LinearAlgebra::Vector<float> &,
- const Function<deal_II_space_dimension> &,
- Vector<double> &,
- const Quadrature<deal_II_dimension> &,
- const NormType &,
- const Function<deal_II_space_dimension> *,
- const double);
-
-
- template void
- integrate_difference<deal_II_dimension,
- LinearAlgebra::distributed::Vector<float>,
- Vector<double>,
- deal_II_space_dimension>(
- const Mapping<deal_II_dimension, deal_II_space_dimension> &,
- const DoFHandler<deal_II_dimension, deal_II_space_dimension> &,
- const LinearAlgebra::distributed::Vector<float> &,
- const Function<deal_II_space_dimension> &,
- Vector<double> &,
- const Quadrature<deal_II_dimension> &,
- const NormType &,
- const Function<deal_II_space_dimension> *,
- const double);
-
- template void
- integrate_difference<deal_II_dimension,
- LinearAlgebra::distributed::BlockVector<float>,
- Vector<double>,
- deal_II_space_dimension>(
- const Mapping<deal_II_dimension, deal_II_space_dimension> &,
- const DoFHandler<deal_II_dimension, deal_II_space_dimension> &,
- const LinearAlgebra::distributed::BlockVector<float> &,
- const Function<deal_II_space_dimension> &,
- Vector<double> &,
- const Quadrature<deal_II_dimension> &,
- const NormType &,
- const Function<deal_II_space_dimension> *,
- const double);
-
-# ifdef DEAL_II_TRILINOS_WITH_TPETRA
- template void
- integrate_difference<deal_II_dimension,
- LinearAlgebra::TpetraWrappers::Vector<float>,
- Vector<double>,
- deal_II_space_dimension>(
- const Mapping<deal_II_dimension, deal_II_space_dimension> &,
- const DoFHandler<deal_II_dimension, deal_II_space_dimension> &,
- const LinearAlgebra::TpetraWrappers::Vector<float> &,
- const Function<deal_II_space_dimension> &,
- Vector<double> &,
- const Quadrature<deal_II_dimension> &,
- const NormType &,
- const Function<deal_II_space_dimension> *,
- const double);
-# endif
-
- template void
- integrate_difference<deal_II_dimension,
- BlockVector<float>,
- Vector<double>,
- deal_II_space_dimension>(
- const DoFHandler<deal_II_dimension, deal_II_space_dimension> &,
- const BlockVector<float> &,
- const Function<deal_II_space_dimension> &,
- Vector<double> &,
- const Quadrature<deal_II_dimension> &,
- const NormType &,
- const Function<deal_II_space_dimension> *,
- const double);
-
- template void
- integrate_difference<deal_II_dimension,
- LinearAlgebra::Vector<float>,
- Vector<double>,
- deal_II_space_dimension>(
- const DoFHandler<deal_II_dimension, deal_II_space_dimension> &,
- const LinearAlgebra::Vector<float> &,
- const Function<deal_II_space_dimension> &,
- Vector<double> &,
- const Quadrature<deal_II_dimension> &,
- const NormType &,
- const Function<deal_II_space_dimension> *,
- const double);
-
- template void
- integrate_difference<deal_II_dimension,
- LinearAlgebra::distributed::Vector<float>,
- Vector<double>,
- deal_II_space_dimension>(
- const DoFHandler<deal_II_dimension, deal_II_space_dimension> &,
- const LinearAlgebra::distributed::Vector<float> &,
- const Function<deal_II_space_dimension> &,
- Vector<double> &,
- const Quadrature<deal_II_dimension> &,
- const NormType &,
- const Function<deal_II_space_dimension> *,
- const double);
-
- template void
- integrate_difference<deal_II_dimension,
- LinearAlgebra::distributed::BlockVector<float>,
- Vector<double>,
- deal_II_space_dimension>(
- const DoFHandler<deal_II_dimension, deal_II_space_dimension> &,
- const LinearAlgebra::distributed::BlockVector<float> &,
- const Function<deal_II_space_dimension> &,
- Vector<double> &,
- const Quadrature<deal_II_dimension> &,
- const NormType &,
- const Function<deal_II_space_dimension> *,
- const double);
-
-# ifdef DEAL_II_TRILINOS_WITH_TPETRA
- template void
- integrate_difference<deal_II_dimension,
- LinearAlgebra::TpetraWrappers::Vector<float>,
- Vector<double>,
- deal_II_space_dimension>(
- const DoFHandler<deal_II_dimension, deal_II_space_dimension> &,
- const LinearAlgebra::TpetraWrappers::Vector<float> &,
- const Function<deal_II_space_dimension> &,
- Vector<double> &,
- const Quadrature<deal_II_dimension> &,
- const NormType &,
- const Function<deal_II_space_dimension> *,
- const double);
-# endif
-
- template void
- integrate_difference<deal_II_dimension,
- BlockVector<float>,
- Vector<double>,
- deal_II_space_dimension>(
- const hp::MappingCollection<deal_II_dimension, deal_II_space_dimension>
- &,
- const DoFHandler<deal_II_dimension, deal_II_space_dimension> &,
- const BlockVector<float> &,
- const Function<deal_II_space_dimension> &,
- Vector<double> &,
- const hp::QCollection<deal_II_dimension> &,
- const NormType &,
- const Function<deal_II_space_dimension> *,
- const double);
-
- template void
- integrate_difference<deal_II_dimension,
- LinearAlgebra::Vector<float>,
- Vector<double>,
- deal_II_space_dimension>(
- const hp::MappingCollection<deal_II_dimension, deal_II_space_dimension>
- &,
- const DoFHandler<deal_II_dimension, deal_II_space_dimension> &,
- const LinearAlgebra::Vector<float> &,
- const Function<deal_II_space_dimension> &,
- Vector<double> &,
- const hp::QCollection<deal_II_dimension> &,
- const NormType &,
- const Function<deal_II_space_dimension> *,
- const double);
-
- template void
- integrate_difference<deal_II_dimension,
- LinearAlgebra::distributed::Vector<float>,
- Vector<double>,
- deal_II_space_dimension>(
- const hp::MappingCollection<deal_II_dimension, deal_II_space_dimension>
- &,
- const DoFHandler<deal_II_dimension, deal_II_space_dimension> &,
- const LinearAlgebra::distributed::Vector<float> &,
- const Function<deal_II_space_dimension> &,
- Vector<double> &,
- const hp::QCollection<deal_II_dimension> &,
- const NormType &,
- const Function<deal_II_space_dimension> *,
- const double);
-
- template void
- integrate_difference<deal_II_dimension,
- LinearAlgebra::distributed::BlockVector<float>,
- Vector<double>,
- deal_II_space_dimension>(
- const hp::MappingCollection<deal_II_dimension, deal_II_space_dimension>
- &,
- const DoFHandler<deal_II_dimension, deal_II_space_dimension> &,
- const LinearAlgebra::distributed::BlockVector<float> &,
- const Function<deal_II_space_dimension> &,
- Vector<double> &,
- const hp::QCollection<deal_II_dimension> &,
- const NormType &,
- const Function<deal_II_space_dimension> *,
- const double);
-
-# ifdef DEAL_II_TRILINOS_WITH_TPETRA
- template void
- integrate_difference<deal_II_dimension,
- LinearAlgebra::TpetraWrappers::Vector<float>,
- Vector<double>,
- deal_II_space_dimension>(
- const hp::MappingCollection<deal_II_dimension, deal_II_space_dimension>
- &,
- const DoFHandler<deal_II_dimension, deal_II_space_dimension> &,
- const LinearAlgebra::TpetraWrappers::Vector<float> &,
- const Function<deal_II_space_dimension> &,
- Vector<double> &,
- const hp::QCollection<deal_II_dimension> &,
- const NormType &,
- const Function<deal_II_space_dimension> *,
- const double);
-# endif
-
- template void
- integrate_difference<deal_II_dimension,
- BlockVector<float>,
- Vector<double>,
- deal_II_space_dimension>(
- const DoFHandler<deal_II_dimension, deal_II_space_dimension> &,
- const BlockVector<float> &,
- const Function<deal_II_space_dimension> &,
- Vector<double> &,
- const hp::QCollection<deal_II_dimension> &,
- const NormType &,
- const Function<deal_II_space_dimension> *,
- const double);
-
- template void
- integrate_difference<deal_II_dimension,
- LinearAlgebra::Vector<float>,
- Vector<double>,
- deal_II_space_dimension>(
- const DoFHandler<deal_II_dimension, deal_II_space_dimension> &,
- const LinearAlgebra::Vector<float> &,
- const Function<deal_II_space_dimension> &,
- Vector<double> &,
- const hp::QCollection<deal_II_dimension> &,
- const NormType &,
- const Function<deal_II_space_dimension> *,
- const double);
-
- template void
- integrate_difference<deal_II_dimension,
- LinearAlgebra::distributed::Vector<float>,
- Vector<double>,
- deal_II_space_dimension>(
- const DoFHandler<deal_II_dimension, deal_II_space_dimension> &,
- const LinearAlgebra::distributed::Vector<float> &,
- const Function<deal_II_space_dimension> &,
- Vector<double> &,
- const hp::QCollection<deal_II_dimension> &,
- const NormType &,
- const Function<deal_II_space_dimension> *,
- const double);
-
- template void
- integrate_difference<deal_II_dimension,
- LinearAlgebra::distributed::BlockVector<float>,
- Vector<double>,
- deal_II_space_dimension>(
- const DoFHandler<deal_II_dimension, deal_II_space_dimension> &,
- const LinearAlgebra::distributed::BlockVector<float> &,
- const Function<deal_II_space_dimension> &,
- Vector<double> &,
- const hp::QCollection<deal_II_dimension> &,
- const NormType &,
- const Function<deal_II_space_dimension> *,
- const double);
-
-# ifdef DEAL_II_TRILINOS_WITH_TPETRA
- template void
- integrate_difference<deal_II_dimension,
- LinearAlgebra::TpetraWrappers::Vector<float>,
- Vector<double>,
- deal_II_space_dimension>(
- const DoFHandler<deal_II_dimension, deal_II_space_dimension> &,
- const LinearAlgebra::TpetraWrappers::Vector<float> &,
- const Function<deal_II_space_dimension> &,
- Vector<double> &,
- const hp::QCollection<deal_II_dimension> &,
- const NormType &,
- const Function<deal_II_space_dimension> *,
- const double);
-# endif
-
- template void
- integrate_difference<deal_II_dimension,
- BlockVector<float>,
- Vector<float>,
- deal_II_space_dimension>(
- const hp::MappingCollection<deal_II_dimension, deal_II_space_dimension>
- &,
- const DoFHandler<deal_II_dimension, deal_II_space_dimension> &,
- const BlockVector<float> &,
- const Function<deal_II_space_dimension> &,
- Vector<float> &,
- const hp::QCollection<deal_II_dimension> &,
- const NormType &,
- const Function<deal_II_space_dimension> *,
- const double);
-
- template void
- integrate_difference<deal_II_dimension,
- LinearAlgebra::Vector<float>,
- Vector<float>,
- deal_II_space_dimension>(
- const hp::MappingCollection<deal_II_dimension, deal_II_space_dimension>
- &,
- const DoFHandler<deal_II_dimension, deal_II_space_dimension> &,
- const LinearAlgebra::Vector<float> &,
- const Function<deal_II_space_dimension> &,
- Vector<float> &,
- const hp::QCollection<deal_II_dimension> &,
- const NormType &,
- const Function<deal_II_space_dimension> *,
- const double);
-
- template void
- integrate_difference<deal_II_dimension,
- LinearAlgebra::distributed::Vector<float>,
- Vector<float>,
- deal_II_space_dimension>(
- const hp::MappingCollection<deal_II_dimension, deal_II_space_dimension>
- &,
- const DoFHandler<deal_II_dimension, deal_II_space_dimension> &,
- const LinearAlgebra::distributed::Vector<float> &,
- const Function<deal_II_space_dimension> &,
- Vector<float> &,
- const hp::QCollection<deal_II_dimension> &,
- const NormType &,
- const Function<deal_II_space_dimension> *,
- const double);
-
- template void
- integrate_difference<deal_II_dimension,
- LinearAlgebra::distributed::BlockVector<float>,
- Vector<float>,
- deal_II_space_dimension>(
- const hp::MappingCollection<deal_II_dimension, deal_II_space_dimension>
- &,
- const DoFHandler<deal_II_dimension, deal_II_space_dimension> &,
- const LinearAlgebra::distributed::BlockVector<float> &,
- const Function<deal_II_space_dimension> &,
- Vector<float> &,
- const hp::QCollection<deal_II_dimension> &,
- const NormType &,
- const Function<deal_II_space_dimension> *,
- const double);
-
-# ifdef DEAL_II_TRILINOS_WITH_TPETRA
- template void
- integrate_difference<deal_II_dimension,
- LinearAlgebra::TpetraWrappers::Vector<float>,
- Vector<float>,
- deal_II_space_dimension>(
- const hp::MappingCollection<deal_II_dimension, deal_II_space_dimension>
- &,
- const DoFHandler<deal_II_dimension, deal_II_space_dimension> &,
- const LinearAlgebra::TpetraWrappers::Vector<float> &,
- const Function<deal_II_space_dimension> &,
- Vector<float> &,
- const hp::QCollection<deal_II_dimension> &,
- const NormType &,
- const Function<deal_II_space_dimension> *,
- const double);
-# endif
-
- template void
- integrate_difference<deal_II_dimension,
- BlockVector<float>,
- Vector<float>,
- deal_II_space_dimension>(
- const DoFHandler<deal_II_dimension, deal_II_space_dimension> &,
- const BlockVector<float> &,
- const Function<deal_II_space_dimension> &,
- Vector<float> &,
- const hp::QCollection<deal_II_dimension> &,
- const NormType &,
- const Function<deal_II_space_dimension> *,
- const double);
-
- template void
- integrate_difference<deal_II_dimension,
- LinearAlgebra::Vector<float>,
- Vector<float>,
- deal_II_space_dimension>(
- const DoFHandler<deal_II_dimension, deal_II_space_dimension> &,
- const LinearAlgebra::Vector<float> &,
- const Function<deal_II_space_dimension> &,
- Vector<float> &,
- const hp::QCollection<deal_II_dimension> &,
- const NormType &,
- const Function<deal_II_space_dimension> *,
- const double);
-
- template void
- integrate_difference<deal_II_dimension,
- LinearAlgebra::distributed::Vector<float>,
- Vector<float>,
- deal_II_space_dimension>(
- const DoFHandler<deal_II_dimension, deal_II_space_dimension> &,
- const LinearAlgebra::distributed::Vector<float> &,
- const Function<deal_II_space_dimension> &,
- Vector<float> &,
- const hp::QCollection<deal_II_dimension> &,
- const NormType &,
- const Function<deal_II_space_dimension> *,
- const double);
-
- template void
- integrate_difference<deal_II_dimension,
- LinearAlgebra::distributed::BlockVector<float>,
- Vector<float>,
- deal_II_space_dimension>(
- const DoFHandler<deal_II_dimension, deal_II_space_dimension> &,
- const LinearAlgebra::distributed::BlockVector<float> &,
- const Function<deal_II_space_dimension> &,
- Vector<float> &,
- const hp::QCollection<deal_II_dimension> &,
- const NormType &,
- const Function<deal_II_space_dimension> *,
- const double);
-
-# ifdef DEAL_II_TRILINOS_WITH_TPETRA
- template void
- integrate_difference<deal_II_dimension,
- LinearAlgebra::TpetraWrappers::Vector<float>,
- Vector<float>,
- deal_II_space_dimension>(
- const DoFHandler<deal_II_dimension, deal_II_space_dimension> &,
- const LinearAlgebra::TpetraWrappers::Vector<float> &,
- const Function<deal_II_space_dimension> &,
- Vector<float> &,
- const hp::QCollection<deal_II_dimension> &,
- const NormType &,
- const Function<deal_II_space_dimension> *,
- const double);
-# endif
-
- \}
-#endif
- }
-
for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS)
{
#if deal_II_dimension <= deal_II_space_dimension
projects / dealii.git / commitdiff