From: Markus Buerg Date: Fri, 13 Aug 2010 08:47:34 +0000 (+0000) Subject: Functions for computing Dirichlet boundary conditions for Nedelec elements. X-Git-Tag: v8.0.0~5734 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=1ee70ff616e70e71e896925553c9e6ea5c8c8fc8;p=dealii.git Functions for computing Dirichlet boundary conditions for Nedelec elements. git-svn-id: https://svn.dealii.org/trunk@21650 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/numerics/vectors.h b/deal.II/deal.II/include/numerics/vectors.h index 0acec10fbb..ea3ae3c338 100644 --- a/deal.II/deal.II/include/numerics/vectors.h +++ b/deal.II/deal.II/include/numerics/vectors.h @@ -18,7 +18,8 @@ #include #include #include -#include +#include +#include #include #include @@ -944,6 +945,99 @@ class VectorTools std::vector component_mapping = std::vector()); + /** + * 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 ConstraintMatrix @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. + * + * This function is explecitly 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_indicator of the face. 255 + * 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 \v{S}olin, Segeth and Dole\v{z}el + * (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. + */ + template + static void project_boundary_values_curl_conforming (const DoFHandler& dof_handler, + const unsigned int first_vector_component, + const Function& boundary_function, + const unsigned char boundary_component, + ConstraintMatrix& constraints, + const Mapping& mapping = StaticMappingQ1::mapping); + + /** + * Same as above for the hp-namespace. + */ + template + static void project_boundary_values_curl_conforming (const hp::DoFHandler& dof_handler, + const unsigned int first_vector_component, + const Function& boundary_function, + const unsigned char boundary_component, + ConstraintMatrix& constraints, + const hp::MappingCollection& mapping_collection = hp::StaticMappingQ1::mapping_collection); + /** * Compute the constraints that * correspond to boundary conditions of diff --git a/deal.II/deal.II/include/numerics/vectors.templates.h b/deal.II/deal.II/include/numerics/vectors.templates.h index 132b0fd96f..0260ff79de 100644 --- a/deal.II/deal.II/include/numerics/vectors.templates.h +++ b/deal.II/deal.II/include/numerics/vectors.templates.h @@ -32,14 +32,19 @@ #include #include #include +#include #include #include #include #include +#include #include #include #include #include +#include +#include +#include #include #include @@ -2913,10 +2918,811 @@ namespace internal Assert (false, ExcNotImplemented()); } } + } +} + +namespace internals { + namespace VectorTools { + + // 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, FEValues& fe_values, const Quadrature& quadrature, + const Function& boundary_function, const unsigned int first_vector_component, + std::vector& dof_values) { + fe_values.reinit (cell); + + // Initialize the required objects. + std::vector > jacobians = fe_values.get_jacobians (); + std::vector > tangentials (fe_values.n_quadrature_points); + std::vector > quadrature_points = fe_values.get_quadrature_points (); + std::vector > values (fe_values.n_quadrature_points, + Vector (dim)); + + // Get boundary function values at quadrature points. + boundary_function.vector_value_list (quadrature_points, values); + quadrature_points = quadrature.get_points (); + + const unsigned int superdegree = cell->get_fe ().degree; + const unsigned int degree = superdegree - 1; + Point shifted_reference_point_1; + Point shifted_reference_point_2; + unsigned int edge_coordinate_direction[4]; + + // Get coordinate directions of the edges of the face. + switch (face) { + case 0: case 1: { + edge_coordinate_direction[0] = 2; + edge_coordinate_direction[1] = 2; + edge_coordinate_direction[2] = 1; + edge_coordinate_direction[3] = 1; + break; + } + + case 2: case 3: { + edge_coordinate_direction[0] = 0; + edge_coordinate_direction[1] = 0; + edge_coordinate_direction[2] = 2; + edge_coordinate_direction[3] = 2; + break; + } + + default: { + edge_coordinate_direction[0] = 1; + edge_coordinate_direction[1] = 1; + edge_coordinate_direction[2] = 0; + edge_coordinate_direction[3] = 0; + } + } + + // 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. + for (unsigned int d = 0; d < dim; ++d) { + shifted_reference_point_1 (d) = quadrature_points[q_point] (d); + shifted_reference_point_2 (d) = quadrature_points[q_point] (d); + } + + shifted_reference_point_1 (edge_coordinate_direction[line]) += 1e-13; + shifted_reference_point_2 (edge_coordinate_direction[line]) -= 1e-13; + tangentials[q_point] = 2e13 * + (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)); + tangentials[q_point] /= std::sqrt (tangentials[q_point].square ()); + // Compute the mean value. + dof_values[line * superdegree] += fe_values.JxW (q_point) + * (values[q_point] (0) * tangentials[q_point] (0) + values[q_point] (1) + * tangentials[q_point] (1) + values[q_point] (2) * tangentials[q_point] (2)) + / (jacobians[q_point][0][edge_coordinate_direction[line]] + * jacobians[q_point][0][edge_coordinate_direction[line]] + + jacobians[q_point][1][edge_coordinate_direction[line]] + * jacobians[q_point][1][edge_coordinate_direction[line]] + + jacobians[q_point][2][edge_coordinate_direction[line]] + * jacobians[q_point][2][edge_coordinate_direction[line]]); + } + + // If there are also higher order shape functions we have + // still some work left. + if (degree > 0) { + const FEValuesExtractors::Vector vec (first_vector_component); + FullMatrix assembling_matrix (degree, fe_values.n_quadrature_points); + Tensor<1, dim> shape_value; + Tensor<1, dim> tmp; + Vector assembling_vector (fe_values.n_quadrature_points); + + // We set up a linear system of equations to get the values + // for the remaining degrees of freedom associated with + // the edge. + 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 + // tangential components of the residual of the boundary + // function and the interpolated part above. + tmp = std::sqrt (fe_values.JxW (q_point) + / (jacobians[q_point][0][edge_coordinate_direction[line]] + * jacobians[q_point][0][edge_coordinate_direction[line]] + + jacobians[q_point][1][edge_coordinate_direction[line]] + * jacobians[q_point][1][edge_coordinate_direction[line]] + + jacobians[q_point][2][edge_coordinate_direction[line]] + * jacobians[q_point][2][edge_coordinate_direction[line]])) + * tangentials[q_point]; + shape_value + = fe_values[vec].value (cell->get_fe ().face_to_cell_index (line * superdegree, face), + q_point); + // In the weak form the right hand side function is multiplicated + // by the higher order shape functions. + assembling_vector (q_point) = (values[q_point] (0) + - dof_values[line * superdegree] * shape_value[0]) * tmp[0] + (values[q_point] (1) + - dof_values[line * superdegree] * shape_value[1]) * tmp[1] + (values[q_point] (2) + - dof_values[line * superdegree] * shape_value[2]) * tmp[2]; + + for (unsigned int i = 0; i < degree; ++i) + assembling_matrix (i, q_point) + = fe_values[vec].value (cell->get_fe ().face_to_cell_index (i + line * superdegree + 1, face), + q_point) * tmp; + } + + FullMatrix cell_matrix (degree, degree); + + // Create the system matrix by multiplying the assembling + // matrix with its transposed. + assembling_matrix.mTmult (cell_matrix, assembling_matrix); + + Vector cell_rhs (degree); + + // Create the system right hand side vector by multiplying + // the assembling matrix with the assembling vector. + assembling_matrix.vmult (cell_rhs, assembling_vector); + + PreconditionJacobi > precondition; + + // Use Jacobi preconditioner with the PCG method to solve the + // problem. + precondition.initialize (cell_matrix); + + SolverControl solver_control (degree, 1e-15, false, false); + SolverCG<> cg (solver_control); + Vector solution (degree); + + cg.solve (cell_matrix, solution, cell_rhs, precondition); + + // Store the computed values. + for (unsigned int i = 0; i < degree; ++i) + dof_values[i + line * superdegree + 1] = solution (i); + } + } + + // This function computes the projection of the + // boundary function on the interior of faces in + // 3D. + template + void + compute_face_projection (const cell_iterator& cell, const unsigned int face, + FEValues& fe_values, const Function& boundary_function, const unsigned + int first_vector_component, + std::vector& dof_values) { + fe_values.reinit (cell); + + // Initialize the required objects. + std::vector > jacobians = fe_values.get_jacobians (); + std::vector > quadrature_points = fe_values.get_quadrature_points (); + std::vector > values (fe_values.n_quadrature_points, + Vector (dim)); + + // Get boundary function values at quadrature points. + boundary_function.vector_value_list (quadrature_points, values); + + const FEValuesExtractors::Vector vec (first_vector_component); + const unsigned int superdegree = cell->get_fe ().degree; + const unsigned int degree = superdegree - 1; + double JxW; + FullMatrix assembling_matrix (degree * superdegree, + 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 ()); + Vector solution (cell_matrix.m ()); + SolverControl solver_control (cell_matrix.m (), 1e-15, false, false); + SolverCG<> cg (solver_control); + PreconditionJacobi > precondition; + Tensor<1, dim> tmp; + Tensor<1, dim> shape_value; + unsigned int global_face_coordinate_directions[2]; + unsigned int local_face_coordinate_directions[2]; + + // Get coordinate directions of the face. + switch (face) { + case 0: case 1: { + global_face_coordinate_directions[0] = 1; + global_face_coordinate_directions[1] = 2; + local_face_coordinate_directions[0] = 1; + local_face_coordinate_directions[1] = 0; + break; + } + + case 2: case 3: { + global_face_coordinate_directions[0] = 0; + global_face_coordinate_directions[1] = 2; + local_face_coordinate_directions[0] = 0; + local_face_coordinate_directions[1] = 1; + break; + } + + default: { + global_face_coordinate_directions[0] = 0; + global_face_coordinate_directions[1] = 1; + local_face_coordinate_directions[0] = 1; + local_face_coordinate_directions[1] = 0; + } + } + + // The projection is divided into two steps. In the first step we + // project the boundary function on the horizontal shape functions. + // Then the bounary 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. + for (unsigned int d = 0; d < dim; ++d) + tmp[d] = values[q_point] (d); + + for (unsigned int i = 0; i < 2; ++i) + for (unsigned int j = 0; j <= degree; ++j) + tmp -= dof_values[(i + 2 * local_face_coordinate_directions[0]) * superdegree + j] + * fe_values[vec].value (cell->get_fe ().face_to_cell_index + ((i + 2 * local_face_coordinate_directions[0]) * superdegree + j, + face), q_point); + + JxW = std::sqrt (fe_values.JxW (q_point) + / ((jacobians[q_point][0][global_face_coordinate_directions[0]] + * jacobians[q_point][0][global_face_coordinate_directions[0]] + + jacobians[q_point][1][global_face_coordinate_directions[0]] + * jacobians[q_point][1][global_face_coordinate_directions[0]] + + jacobians[q_point][2][global_face_coordinate_directions[0]] + * jacobians[q_point][2][global_face_coordinate_directions[0]]) + * (jacobians[q_point][0][global_face_coordinate_directions[1]] + * jacobians[q_point][0][global_face_coordinate_directions[1]] + + jacobians[q_point][1][global_face_coordinate_directions[1]] + * jacobians[q_point][1][global_face_coordinate_directions[1]] + + jacobians[q_point][2][global_face_coordinate_directions[1]] + * jacobians[q_point][2][global_face_coordinate_directions[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]; + + for (unsigned int i = 0; i <= degree; ++i) + for (unsigned int j = 0; j < degree; ++j) { + shape_value = JxW + * fe_values[vec].value (cell->get_fe ().face_to_cell_index + ((i + GeometryInfo::lines_per_face) * degree + j + + GeometryInfo::lines_per_face, face), q_point); + + for (unsigned int d = 0; d < dim; ++d) + assembling_matrix (i * degree + j, dim * q_point + d) = shape_value[d]; + } + } + + // Create the system matrix by multiplying the assembling + // matrix with its transposed and the right hand side vector + // by mutliplying the assembling matrix with the assembling + // vector. The problem is solved by the PCG method. + assembling_matrix.mTmult (cell_matrix, assembling_matrix); + assembling_matrix.vmult (cell_rhs, assembling_vector); + precondition.initialize (cell_matrix); + cg.solve (cell_matrix, solution, cell_rhs, precondition); + + // Store the computed values. + for (unsigned int i = 0; i <= degree; ++i) + for (unsigned int j = 0; j < degree; ++j) + dof_values[(i + GeometryInfo::lines_per_face) * degree + j + + GeometryInfo::lines_per_face] = solution (i * degree + j); + + // 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) { + for (unsigned int d = 0; d < dim; ++d) + tmp[d] = values[q_point] (d); + + for (unsigned int i = 0; i < 2; ++i) + for (unsigned int j = 0; j <= degree; ++j) + tmp + -= dof_values[(i + 2 * local_face_coordinate_directions[1]) * superdegree + j] + * fe_values[vec].value (cell->get_fe ().face_to_cell_index + ((i + 2 * local_face_coordinate_directions[1]) * superdegree + j, + face), q_point); + + JxW = std::sqrt (fe_values.JxW (q_point) + / ((jacobians[q_point][0][global_face_coordinate_directions[0]] + * jacobians[q_point][0][global_face_coordinate_directions[0]] + + jacobians[q_point][1][global_face_coordinate_directions[0]] + * jacobians[q_point][1][global_face_coordinate_directions[0]] + + jacobians[q_point][2][global_face_coordinate_directions[0]] + * jacobians[q_point][2][global_face_coordinate_directions[0]]) + * (jacobians[q_point][0][global_face_coordinate_directions[1]] + * jacobians[q_point][0][global_face_coordinate_directions[1]] + + jacobians[q_point][1][global_face_coordinate_directions[1]] + * jacobians[q_point][1][global_face_coordinate_directions[1]] + + jacobians[q_point][2][global_face_coordinate_directions[1]] + * jacobians[q_point][2][global_face_coordinate_directions[1]]))); + + for (unsigned int d = 0; d < dim; ++d) + assembling_vector (dim * q_point + d) = JxW * tmp[d]; + + for (unsigned int i = 0; i < degree; ++i) + for (unsigned int j = 0; j <= degree; ++j) { + shape_value = JxW + * fe_values[vec].value (cell->get_fe ().face_to_cell_index + ((i + degree + GeometryInfo::lines_per_face) * superdegree + j, + face), q_point); + + for (unsigned int d = 0; d < dim; ++d) + assembling_matrix (i * superdegree + j, dim * q_point + d) = shape_value[d]; + } + } + + assembling_matrix.mTmult (cell_matrix, assembling_matrix); + assembling_matrix.vmult (cell_rhs, assembling_vector); + precondition.initialize (cell_matrix); + cg.solve (cell_matrix, solution, cell_rhs, precondition); + + for (unsigned int i = 0; i < degree; ++i) + for (unsigned int j = 0; j <= degree; ++j) + dof_values[(i + degree + GeometryInfo::lines_per_face) * superdegree + j] + = solution (i * superdegree + j); + } + + + // This function computes the projection of the + // boundary function on the faces in 2D. + template + void + compute_face_projection (const cell_iterator& cell, const unsigned int face, + FEValues& fe_values, const Quadrature& quadrature, const Function& + boundary_function, const unsigned int first_vector_component, std::vector& + dof_values) { + fe_values.reinit (cell); + + // Initialize the required objects. + std::vector > tangentials (fe_values.n_quadrature_points); + std::vector > jacobians = fe_values.get_jacobians (); + std::vector > quadrature_points = fe_values.get_quadrature_points (); + std::vector > values (fe_values.n_quadrature_points, + Vector (dim)); + + // Get boundary function values at quadrature points. + boundary_function.vector_value_list (quadrature_points, values); + quadrature_points = quadrature.get_points (); + + const unsigned int degree = cell->get_fe ().degree - 1; + Point shifted_reference_point_1; + Point shifted_reference_point_2; + Tensor<1, dim> tmp; + Tensor<1, dim> shape_value; + unsigned int face_coordinate_direction; + + // Get coordinate directions of the face. + switch (face) { + case 0: case 1: { + face_coordinate_direction = 1; + break; + } + + default: + face_coordinate_direction = 0; + } + + // 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. + for (unsigned int d = 0; d < dim; ++d) { + shifted_reference_point_1 (d) = quadrature_points[q_point] (d); + shifted_reference_point_2 (d) = quadrature_points[q_point] (d); + } + + shifted_reference_point_1 (face_coordinate_direction) += 1e-13; + shifted_reference_point_2 (face_coordinate_direction) -= 1e-13; + tangentials[q_point] = 2e13 + * (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)); + tangentials[q_point] /= std::sqrt (tangentials[q_point].square ()); + // Compute the mean value. + dof_values[0] += fe_values.JxW (q_point) * (values[q_point] (0) + * tangentials[q_point] (0) + values[q_point] (1) * tangentials[q_point] (1)) + / (jacobians[q_point][0][face_coordinate_direction] + * jacobians[q_point][0][face_coordinate_direction] + + jacobians[q_point][1][face_coordinate_direction] + * jacobians[q_point][1][face_coordinate_direction]); + } + + // If there are also higher order shape functions we have + // still some work left. + if (degree > 0) { + const FEValuesExtractors::Vector vec (first_vector_component); + FullMatrix assembling_matrix (degree, fe_values.n_quadrature_points); + Vector assembling_vector (fe_values.n_quadrature_points); + + // We set up a linear system of equations to get the values + // for the remaining degrees of freedom associated with + // 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 + // tangential components of the residual of the boundary + // function and the interpolated part above. + tmp = std::sqrt (fe_values.JxW (q_point) + / std::sqrt (jacobians[q_point][0][face_coordinate_direction] + * jacobians[q_point][0][face_coordinate_direction] + + jacobians[q_point][1][face_coordinate_direction] + * jacobians[q_point][1][face_coordinate_direction])) * tangentials[q_point]; + shape_value + = fe_values[vec].value (cell->get_fe ().face_to_cell_index (0, face), q_point); + assembling_vector (q_point) = (values[q_point] (0) + - dof_values[0] * shape_value[0]) * tmp[0] + (values[q_point] (1) + - dof_values[1] * shape_value[1]) * tmp[1]; + + // In the weak form the right hand side function is multiplicated + // by the higher order shape functions. + for (unsigned int i = 0; i < degree; ++i) + assembling_matrix (i, q_point) + = fe_values[vec].value (cell->get_fe ().face_to_cell_index (i + 1, face), + q_point) * tmp; + } + + FullMatrix cell_matrix (degree, degree); + + // Create the system matrix by multiplying the assembling + // matrix with its transposed. + assembling_matrix.mTmult (cell_matrix, assembling_matrix); + + Vector cell_rhs (degree); + + // Create the system right hand side vector by multiplying + // the assembling matrix with the assembling vector. + assembling_matrix.vmult (cell_rhs, assembling_vector); + + PreconditionJacobi > precondition; + + // Use Jacobi preconditioner with the PCG method to solve the + // problem. + precondition.initialize (cell_matrix); + + SolverControl solver_control (degree, 1e-15, false, false); + SolverCG<> cg (solver_control); + Vector solution (degree); + + cg.solve (cell_matrix, solution, cell_rhs, precondition); + + // Store the computed values. + for (unsigned int i = 0; i < degree; ++i) + dof_values[i + 1] = solution (i); + } + } } } + + + // 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 lowest order edge shape functions. Then the interpolated + // part of the function is subtracted and we project the tangential + // component of the residual onto the space of the remaining + // (higher order) edge shape functions. This is done by building + // a linear system of equations of dimension degree. The + // solution gives us the values for the degrees of freedom + // corresponding to the remaining 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 two steps. On the faces there are two kinds of shape + // functions, the horizontal and the vertical ones. Thus we have + // two solve two linear systems of equations of size + // degree * (degree + 1) to obtain the values for the + // corresponding degrees of freedom. + +template +void VectorTools::project_boundary_values_curl_conforming (const DoFHandler& dof_handler, + const unsigned int first_vector_component, + const Function& boundary_function, + const unsigned char boundary_component, + ConstraintMatrix& constraints, + const Mapping& mapping) +{ + std::vector dof_values; + std::vector face_dof_indices; + typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active (); + unsigned int dofs_per_face; + unsigned int superdegree; + + switch (dim) { + case 2: { + for (; cell != dof_handler.end (); ++cell) + if (cell->at_boundary ()) + for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; ++face) + if (cell->face (face)->boundary_indicator () == boundary_component) { + // this is only implemented, if the FE is a Nedelec element + typedef FiniteElement FEL; + + AssertThrow ((cell->get_fe ().get_name ().find ("FE_Nedelec<") == 0), + typename FEL::ExcInterpolationNotImplemented ()); + + dofs_per_face = cell->get_fe ().dofs_per_face; + dof_values.resize (dofs_per_face); + + for (unsigned int dof = 0; dof < dofs_per_face; ++dof) + dof_values[dof] = 0.0; + + superdegree = cell->get_fe ().degree; + + QGauss reference_face_quadrature (2 * superdegree); + Quadrature face_quadrature + = QProjector::project_to_face (reference_face_quadrature, face); + FEValues fe_face_values (mapping, cell->get_fe (), face_quadrature, + update_jacobians | update_JxW_values | update_quadrature_points | update_values); + + // Compute the projection of the boundary function on the edge. + internals::VectorTools::compute_face_projection (cell, face, fe_face_values, + face_quadrature, boundary_function, first_vector_component, dof_values); + face_dof_indices.resize (dofs_per_face); + cell->face (face)->get_dof_indices (face_dof_indices, cell->active_fe_index ()); + + // Add the computed constraints to the constraint matrix. + for (unsigned int dof = 0; dof < dofs_per_face; ++dof) { + constraints.add_line (face_dof_indices[dof]); + + if (std::abs (dof_values[dof]) > 1e-14) + constraints.set_inhomogeneity (face_dof_indices[dof], dof_values[dof]); + } + } + + break; + } + + case 3: { + const unsigned int n_dofs = dof_handler.n_dofs (); + std::vector computed_constraints (n_dofs); + std::vector projected_dofs (n_dofs); + unsigned int degree; + unsigned int superdegree; + + for (unsigned int dof = 0; dof < n_dofs; ++dof) + projected_dofs[dof] = -1; + + for (; cell != dof_handler.end (); ++cell) + if (cell->at_boundary ()) + for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; ++face) + if (cell->face (face)->boundary_indicator () == boundary_component) { + // this is only implemented, if the FE is a Nedelec element + typedef FiniteElement FEL; + + AssertThrow ((cell->get_fe ().get_name ().find ("FE_Nedelec<") == 0), + typename FEL::ExcInterpolationNotImplemented ()); + + superdegree = cell->get_fe ().degree; + degree = superdegree - 1; + + QGauss reference_edge_quadrature (2 * superdegree); + + dofs_per_face = cell->get_fe ().dofs_per_face; + dof_values.resize (dofs_per_face); + + for (unsigned int dof = 0; dof < dofs_per_face; ++dof) + dof_values[dof] = 0.0; + + face_dof_indices.resize (dofs_per_face); + cell->face (face)->get_dof_indices (face_dof_indices, cell->active_fe_index ()); + + // First we compute the projection on the edges. + for (unsigned int line = 0; line < GeometryInfo<3>::lines_per_face; ++line) { + // If we have reached this edge through another cell before, we do + // not do here anything unless we have a good reason, i.e. a higher + // polynomial degree. + if (projected_dofs[face_dof_indices[line * superdegree]] < (int) degree) { + Quadrature edge_quadrature + = QProjector::project_to_face (QProjector::project_to_face + (reference_edge_quadrature, line), face); + FEValues fe_edge_values (mapping, cell->get_fe (), edge_quadrature, + update_JxW_values | update_jacobians | update_quadrature_points | update_values); + // Compute the projection of the boundary function on the edge. + internals::VectorTools::compute_edge_projection (cell, face, line, + fe_edge_values, edge_quadrature, boundary_function, first_vector_component, + dof_values); + // Mark the projected degrees of freedom. + for (unsigned int dof = line * superdegree; dof < (line + 1) * superdegree; + ++dof) + projected_dofs[face_dof_indices[dof]] = degree; + } + + // If we have computed the values in a previous step of the loop, + // we just copy the values in the local vector. + else + for (unsigned int dof = line * superdegree; dof < (line + 1) * superdegree; + ++dof) + dof_values[dof] = computed_constraints[face_dof_indices[dof]]; + } + + // If there are higher order shape functions, there is still some + // work left. + if (degree > 0) { + QGauss reference_face_quadrature (2 * superdegree); + Quadrature face_quadrature + = QProjector::project_to_face (reference_face_quadrature, face); + FEValues fe_face_values (mapping, cell->get_fe (), face_quadrature, + update_JxW_values | update_jacobians | update_quadrature_points | update_values); + + // Compute the projection of the boundary function on the interior + // of the face. + internals::VectorTools::compute_face_projection (cell, face, fe_face_values, + boundary_function, first_vector_component, dof_values); + + // Mark the projected degrees of freedom. + for (unsigned int dof = GeometryInfo::lines_per_face * superdegree; + dof < dofs_per_face; ++dof) + projected_dofs[face_dof_indices[dof]] = degree; + } + + // Store the computed values in the global vector. + for (unsigned int dof = 0; dof < dofs_per_face; ++dof) + if (std::abs (dof_values[dof]) > 1e-14) + computed_constraints[face_dof_indices[dof]] = dof_values[dof]; + } + + // Add the computed constraints to the constraint matrix. + for (unsigned int dof = 0; dof < n_dofs; ++dof) + if (projected_dofs[dof] != -1) { + constraints.add_line (dof); + constraints.set_inhomogeneity (dof, computed_constraints[dof]); + } + } + } +} + +template +void VectorTools::project_boundary_values_curl_conforming (const hp::DoFHandler& dof_handler, + const unsigned int first_vector_component, + const Function& boundary_function, + const unsigned char boundary_component, + ConstraintMatrix& constraints, + const hp::MappingCollection& mapping_collection) +{ + std::vector dof_values; + std::vector face_dof_indices; + typename hp::DoFHandler::active_cell_iterator cell = dof_handler.begin_active (); + unsigned int dofs_per_face; + + switch (dim) { + case 2: { + for (; cell != dof_handler.end (); ++cell) + if (cell->at_boundary ()) + for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; ++face) + if (cell->face (face)->boundary_indicator () == boundary_component) { + // this is only implemented, if the FE is a Nédélec element + typedef FiniteElement FEL; + + AssertThrow ((cell->get_fe ().get_name ().find ("FE_Nedelec<") == 0), + typename FEL::ExcInterpolationNotImplemented ()); + + dofs_per_face = cell->get_fe ().dofs_per_face; + dof_values.resize (dofs_per_face); + + for (unsigned int dof = 0; dof < dofs_per_face; ++dof) + dof_values[dof] = 0.0; + + QGauss reference_face_quadrature (2 * (cell->get_fe ().degree)); + Quadrature face_quadrature + = QProjector::project_to_face (reference_face_quadrature, face); + FEValues fe_face_values (mapping_collection[cell->active_fe_index ()], + cell->get_fe (), face_quadrature, update_jacobians | update_JxW_values + | update_quadrature_points | update_values); + + internals::VectorTools::compute_face_projection (cell, face, fe_face_values, + face_quadrature, boundary_function, first_vector_component, dof_values); + 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) { + constraints.add_line (face_dof_indices[dof]); + + if (std::abs (dof_values[dof]) > 1e-14) + constraints.set_inhomogeneity (face_dof_indices[dof], dof_values[dof]); + } + } + + break; + } + + case 3: { + const unsigned int n_dofs = dof_handler.n_dofs (); + std::vector computed_constraints (n_dofs); + std::vector projected_dofs (n_dofs); + unsigned int degree; + unsigned int superdegree; + + for (unsigned int dof = 0; dof < n_dofs; ++dof) + projected_dofs[dof] = -1; + + for (; cell != dof_handler.end (); ++cell) + if (cell->at_boundary ()) + for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; ++face) + if (cell->face (face)->boundary_indicator () == boundary_component) { + // this is only implemented, if the FE is a Nédélec element + typedef FiniteElement FEL; + + AssertThrow ((cell->get_fe ().get_name ().find ("FE_Nedelec<") == 0), + typename FEL::ExcInterpolationNotImplemented ()); + + superdegree = cell->get_fe ().degree; + degree = superdegree - 1; + + QGauss reference_edge_quadrature (2 * superdegree); + + dofs_per_face = cell->get_fe ().dofs_per_face; + dof_values.resize (dofs_per_face); + + for (unsigned int dof = 0; dof < dofs_per_face; ++dof) + dof_values[dof] = 0.0; + + face_dof_indices.resize (dofs_per_face); + cell->face (face)->get_dof_indices (face_dof_indices, cell->active_fe_index ()); + + for (unsigned int line = 0; line < GeometryInfo::lines_per_face; ++line) { + if (projected_dofs[face_dof_indices[line * superdegree]] < (int) degree) { + Quadrature edge_quadrature = QProjector::project_to_face + (QProjector::project_to_face (reference_edge_quadrature, line), + face); + FEValues fe_edge_values (mapping_collection[cell->active_fe_index ()], + cell->get_fe (), edge_quadrature, update_JxW_values | update_jacobians + | update_quadrature_points | update_values); + + internals::VectorTools::compute_edge_projection (cell, face, line, + fe_edge_values, edge_quadrature, boundary_function, first_vector_component, + dof_values); + + for (unsigned int dof = line * superdegree; dof < (line + 1) * superdegree; + ++dof) + projected_dofs[face_dof_indices[dof]] = degree; + } + + else + for (unsigned int dof = line * superdegree; dof < (line + 1) * superdegree; + ++dof) + dof_values[dof] = computed_constraints[face_dof_indices[dof]]; + } + + if (degree > 0) { + QGauss reference_face_quadrature (2 * superdegree); + Quadrature face_quadrature + = QProjector::project_to_face (reference_face_quadrature, face); + FEValues fe_face_values (mapping_collection[cell->active_fe_index ()], + cell->get_fe (), face_quadrature, update_JxW_values | update_jacobians + | update_quadrature_points | update_values); + + internals::VectorTools::compute_face_projection (cell, face, fe_face_values, + boundary_function, first_vector_component, dof_values); + + for (unsigned int dof = GeometryInfo::lines_per_face * superdegree; + dof < dofs_per_face; ++dof) + projected_dofs[face_dof_indices[dof]] = degree; + } + + for (unsigned int dof = 0; dof < dofs_per_face; ++dof) + if (std::abs (dof_values[dof]) > 1e-14) + computed_constraints[face_dof_indices[dof]] = dof_values[dof]; + } + + for (unsigned int dof = 0; dof < n_dofs; ++dof) + if (projected_dofs[dof] != -1) { + constraints.add_line (dof); + constraints.set_inhomogeneity (dof, computed_constraints[dof]); + } + } + } +} template class DH, int spacedim> diff --git a/deal.II/deal.II/source/fe/fe_values.cc b/deal.II/deal.II/source/fe/fe_values.cc index 98a9980d01..3c62232cce 100644 --- a/deal.II/deal.II/source/fe/fe_values.cc +++ b/deal.II/deal.II/source/fe/fe_values.cc @@ -862,7 +862,6 @@ namespace FEValuesViews case 0: { for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points; ++q_point) { - curls[q_point][0] = 0.0; curls[q_point][1] += value * (*shape_gradient_ptr)[2]; curls[q_point][2] -= value * (*shape_gradient_ptr++)[1]; } @@ -874,7 +873,6 @@ namespace FEValuesViews for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points; ++q_point) { curls[q_point][0] -= value * (*shape_gradient_ptr)[2]; - curls[q_point][1] = 0.0; curls[q_point][2] += value * (*shape_gradient_ptr++)[0]; } @@ -886,7 +884,6 @@ namespace FEValuesViews q_point < fe_values.n_quadrature_points; ++q_point) { curls[q_point][0] += value * (*shape_gradient_ptr)[1]; curls[q_point][1] -= value * (*shape_gradient_ptr++)[0]; - curls[q_point][2] = 0.0; } } } diff --git a/deal.II/deal.II/source/numerics/vectors.cc b/deal.II/deal.II/source/numerics/vectors.cc index 4b6bc3e594..92bddb10dc 100644 --- a/deal.II/deal.II/source/numerics/vectors.cc +++ b/deal.II/deal.II/source/numerics/vectors.cc @@ -219,6 +219,22 @@ void VectorTools::project_boundary_values #if deal_II_dimension != 1 template +void VectorTools::project_boundary_values_curl_conforming +(const DoFHandler&, + const unsigned int, + const Function&, + const unsigned char, + ConstraintMatrix&, + const Mapping&); +template +void VectorTools::project_boundary_values_curl_conforming +(const hp::DoFHandler&, + const unsigned int, + const Function&, + const unsigned char, + ConstraintMatrix&, + const hp::MappingCollection&); +template void VectorTools::compute_no_normal_flux_constraints (const DoFHandler &dof_handler, const unsigned int first_vector_component, @@ -228,6 +244,62 @@ VectorTools::compute_no_normal_flux_constraints (const DoFHandler::cell_iterator&, + const unsigned int, + FEValues&, + const Quadrature&, + const Function&, + const unsigned int, + std::vector&); +template +void +internals::VectorTools::compute_face_projection (const hp::DoFHandler::cell_iterator&, + const unsigned int, + FEValues&, + const Quadrature&, + const Function&, + const unsigned int, + std::vector&); +template +void +internals::VectorTools::compute_edge_projection (const DoFHandler::cell_iterator&, + const unsigned int, + const unsigned int, + FEValues&, + const Quadrature&, + const Function&, + const unsigned int, + std::vector&); +template +void +internals::VectorTools::compute_edge_projection (const hp::DoFHandler::cell_iterator&, + const unsigned int, + const unsigned int, + FEValues&, + const Quadrature&, + const Function&, + const unsigned int, + std::vector&); +template +void +internals::VectorTools::compute_face_projection (const DoFHandler::cell_iterator&, + const unsigned int, + FEValues&, + const Function&, + const unsigned int, + std::vector&); +template +void +internals::VectorTools::compute_face_projection (const hp::DoFHandler::cell_iterator&, + const unsigned int, + FEValues&, + const Function&, + const unsigned int, + std::vector&); + + // // Due to introducing the DoFHandler as a template parameter, // // the following instantiations are required in 1d // #if deal_II_dimension == 1 diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index e068cbfd19..06531d6f73 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -218,6 +218,24 @@ inconvenience this causes.
    +
  1. +

    + New: The functions VectorTools::project_boundary_values_curl_conforming + are added. They can compute Dirichlet boundary conditions for Nedelec + elements. +
    + (Markus Buerg 2010/08/13) +

    + +
  2. +

    + Fixed: The function FEValuesViews::Vector::get_function_curls produced + wrong results in some cases, because it erased the given vector first. + This is now fixed. +
    + (Markus Buerg 2010/08/13) +

    +
  3. New: Ability to project second-order SymmetricTensor and first-order Tensor objects from the quadrature points to the support points of the cell using FETools::compute_projection_from_quadrature_points