From: Ross Kynch Date: Tue, 14 Oct 2014 13:54:57 +0000 (+0100) Subject: Added project_boundary_values_curl_conforming_l2() to VectorTools for X-Git-Tag: v8.3.0-rc1~421^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F197%2Fhead;p=dealii.git Added project_boundary_values_curl_conforming_l2() to VectorTools for use with FE_Nedelec elements in 2D & 3D. This includes a test code included in tests/deal.II to check it works correctly. This new function uses an L2-projection on edges and faces (only edges in 2D, both for 3D) to compute the constaints and factors in the residual from the edge projection when computing the face constraints. The original function project_boundary_values_curl_conforming() has problems with non-rectangular faces, this one does not. The function assumes that the face orientations are standard (e.g. a sphere may cause problems due to mismatched face orientations). The test shows that project_boundary_curl_conforming_l2() works on a distorted cube (non-rectangular faces) for a simple test case (polynomial function). --- diff --git a/include/deal.II/numerics/vector_tools.h b/include/deal.II/numerics/vector_tools.h index 1c86ec796e..0274d29689 100644 --- a/include/deal.II/numerics/vector_tools.h +++ b/include/deal.II/numerics/vector_tools.h @@ -1110,6 +1110,110 @@ namespace VectorTools ConstraintMatrix &constraints, const hp::MappingCollection &mapping_collection = hp::StaticMappingQ1::mapping_collection); + /** + * This function is an updated version of the project_boundary_values_curl_conforming + * function. The intention is to fix a problem when using the previous function in + * conjunction with non-rectangular geometries (i.e. elements with non-rectangular faces). + * The L2-projection method used has been taken from the paper "Electromagnetic scattering + * simulation using an H (curl) conforming hp finite element method in three dimensions" + * by PD Ledger, K Morgan and O Hassan ( Int. J. Num. Meth. Fluids, Volume 53, Issue 8, pages 1267–1296). + * + * This function will compute constraints that correspond to Dirichlet boundary conditions of the form + * $\vec{n}\times\vec{E}=\vec{n}\times\vec{F}$ + * i.e. the tangential components of $\vec{E}$ and $f$ shall coincide. + * + *

Computing constraints

+ * + * To compute the constraints we use a projection method based upon the paper mentioned + * above. In 2D this is done in a single stage for the edge-based shape functions, regardless + * of the order of the finite element. In 3D this is done in two stages, edges first and then + * faces. + * + * For each cell, each edge, $e$, is projected by solving the linear system $Ax=b$ where $x$ + * is the vector of contraints on degrees of freedom on the edge and + * + * $A_{ij} = \int_{e} (\vec{s}_{i}\cdot\vec{t})(\vec{s}_{j}\cdot\vec{t}) dS$ + * + * $b_{i} = \int_{e} (\vec{s}_{i}\cdot\vec{t})(\vec{F}\cdot\vec{t}) dS$ + * + * with $\vec{s}_{i}$ the $i^{th}$ shape function and $\vec{t}$ the tangent vector. + * + * Once all edge constraints, $x$, have been computed, we may compute the face constraints + * in a similar fashion, taking into account the residuals from the edges. + * + * For each face on the cell, $f$, we solve the linear system By=c where $y$ is the vector of + * constraints on degrees of freedom on the face and + * + * B_{ij} = \int_{f} (\vec{n} \cross \vec{s}_{i}) \cdot (\vec{n} \cross \vec{s}_{j}) dS + * + * $c_{i} = \int_{f} (\vec{n} \cross \vec{r}) \cdot (\vec{n} \cross \vec{s}_i} dS + * + * and $\vec{r} = \vec{F} - \sum_{e \in f} \sum{i \in e} \x_{i}\vec{s}_i}$, the edge residual. + * + * The resulting constraints are then given in the solutions $x$ and $y$. + * + * 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. See the discussion on conflicting constraints in the + * module on @ref constraints + * + *

Arguments to this function>

+ * + * This function is explicitly for use with FE_Nedelec elements, or with FESystem + * elements which contain FE_Nedelec elements. It will throw an exception if called + * with any other finite element. The user must ensure that FESystem elements are + * correctly setup when using this function as this check not possible in this case. + * + * The second argument of this function denotes the first vector component of the + * finite element which corresponds to the vector function that you wish to constrain. + * For example, if we are solving Maxwell's equations in 3D and have 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. The @p boundary_function must return 6 components in this example, with the first 3 + * corresponding to $\vec{E}$ and the second 3 corresponding to $\vec{B}$. + * 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. 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. + * + * + * @ingroup constraints + * + * @see @ref GlossBoundaryIndicator "Glossary entry on boundary indicators" + */ + template + void project_boundary_values_curl_conforming_l2 (const DoFHandler &dof_handler, + const unsigned int first_vector_component, + const Function &boundary_function, + const types::boundary_id boundary_component, + ConstraintMatrix &constraints, + const Mapping &mapping = StaticMappingQ1::mapping); + + + /** + * hp-namespace version of project_boundary_values_curl_conforming_l2 (above). + * + * @ingroup constraints + */ + template + void project_boundary_values_curl_conforming_l2 (const hp::DoFHandler &dof_handler, + const unsigned int first_vector_component, + const Function &boundary_function, + const types::boundary_id boundary_component, + ConstraintMatrix &constraints, + const hp::MappingCollection &mapping_collection = hp::StaticMappingQ1::mapping_collection); + /** * Compute constraints that correspond to boundary conditions of the form diff --git a/include/deal.II/numerics/vector_tools.templates.h b/include/deal.II/numerics/vector_tools.templates.h index d5db44f78f..7c5f9ee516 100644 --- a/include/deal.II/numerics/vector_tools.templates.h +++ b/include/deal.II/numerics/vector_tools.templates.h @@ -3801,6 +3801,926 @@ namespace VectorTools } + namespace internals + { + template + void + compute_edge_projection_l2 (const cell_iterator &cell, + const unsigned int face, + const unsigned int line, + hp::FEValues<3> &hp_fe_values, + const Function<3> &boundary_function, + const unsigned int first_vector_component, + std::vector &dof_values, + std::vector &dofs_processed) + { + // This function computes the L2-projection of the given + // boundary function on 3D edges and returns the constraints + // associated with the edge functions for the given cell. + // + // In the context of this function, by associated DoFs we mean: + // the DoFs corresponding to the group of components making up the vector + // with first component first_vector_component (length dim). + const unsigned int dim = 3; + const unsigned int spacedim = 3; + const FiniteElement &fe = cell->get_fe (); + + // reinit for this cell, face and line. + hp_fe_values.reinit + (cell, + (cell->active_fe_index () * GeometryInfo::faces_per_cell + face) + * GeometryInfo::lines_per_face + line); + + // Initialize the required objects. + const FEValues & + fe_values = hp_fe_values.get_present_fe_values (); + // Store degree as fe.degree-1 + // For nedelec elements FE_Nedelec (0) returns fe.degree = 1. + const unsigned int degree = fe.degree - 1; + + const std::vector > & + quadrature_points = fe_values.get_quadrature_points (); + std::vector > values (fe_values.n_quadrature_points, + Vector (fe.n_components ())); + + // Get boundary function values + // at quadrature points. + boundary_function.vector_value_list (quadrature_points, values); + + // Find the group of vector components (dim of them, + // starting at first_vector_component) are within an FESystem. + // + // If not using FESystem then must be using FE_Nedelec, + // which has one base element and one copy of it (with 3 components). + std::pair base_indices (0, 0); + if (dynamic_cast*> (&cell->get_fe ()) != 0) + { + unsigned int fe_index = 0; + unsigned int fe_index_old = 0; + unsigned int i = 0; + + // Find base element: + // base_indices.first + // + // Then select which copy of that base element + // [ each copy is of length fe.base_element(base_indices.first).n_components() ] + // corresponds to first_vector_component: + // base_index.second + 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 (); + } + + // Find DoFs we want to constrain: + // There are fe.dofs_per_line DoFs associated with the + // given line on the given face on the given cell. + // + // Want to know which of these DoFs (there are degree+1 of interest) + // are associated with the components given by first_vector_component. + // Then we can make a map from the associated line DoFs to the face DoFs. + // + // For a single FE_Nedelec<3> element this is simple: + // We know the ordering of local DoFs goes + // lines -> faces -> cells + // + // For a set of FESystem<3> elements we need to pick out the matching base element and + // the index within this ordering. + // + // We call the map associated_edge_dof_to_face_dof + std::vector associated_edge_dof_to_face_dof (degree + 1); + + // Lowest DoF in the base element allowed for this edge: + const unsigned int lower_bound = + fe.base_element(base_indices.first).face_to_cell_index(line * (degree + 1), face); + // Highest DoF in the base element allowed for this edge: + const unsigned int upper_bound = + fe.base_element(base_indices.first).face_to_cell_index((line + 1) * (degree + 1) - 1, face); + + unsigned int associated_edge_dof_index = 0; + // for (unsigned int face_idx = 0; face_idx < fe.dofs_per_face; ++face_idx) + for (unsigned int line_idx = 0; line_idx < fe.dofs_per_line; ++line_idx) + { + // Assuming DoFs on a face are numbered in order by lines then faces. + // i.e. line 0 has degree+1 dofs numbered 0,..,degree + // line 1 has degree+1 dofs numbered (degree+1),..,2*(degree+1) and so on. + const unsigned int face_idx = line*fe.dofs_per_line + line_idx; + // Note, assuming that the edge orientations are "standard" + // i.e. cell->line_orientation(line) = true. + const unsigned int cell_idx = fe.face_to_cell_index (face_idx, face); + + // Check this cell_idx belongs to the correct base_element, component and line: + if (((dynamic_cast*> (&fe) != 0) + && (fe.system_to_base_index (cell_idx).first == base_indices) + && (lower_bound <= fe.system_to_base_index (cell_idx).second) + && (fe.system_to_base_index (cell_idx).second <= upper_bound)) + || ((dynamic_cast*> (&fe) != 0) + && (line * (degree + 1) <= face_idx) + && (face_idx <= (line + 1) * (degree + 1) - 1))) + { + associated_edge_dof_to_face_dof[associated_edge_dof_index] = face_idx; + ++associated_edge_dof_index; + } + } + // Sanity check: + const unsigned int associated_edge_dofs = associated_edge_dof_index; + Assert (associated_edge_dofs == degree + 1, + ExcMessage ("Error: Unexpected number of 3D edge DoFs")); + + // Matrix and RHS vectors to store linear system: + // We have (degree+1) basis functions for an edge + FullMatrix edge_matrix (degree + 1,degree + 1); + FullMatrix edge_matrix_inv (degree + 1,degree + 1); + Vector edge_rhs (degree + 1); + Vector edge_solution (degree + 1); + + const FEValuesExtractors::Vector vec (first_vector_component); + + // coordinate directions of + // the edges of the face. + const unsigned int + edge_coordinate_direction + [GeometryInfo::faces_per_cell] + [GeometryInfo::lines_per_face] + = { { 2, 2, 1, 1 }, + { 2, 2, 1, 1 }, + { 0, 0, 2, 2 }, + { 0, 0, 2, 2 }, + { 1, 1, 0, 0 }, + { 1, 1, 0, 0 } + }; + + const double tol = 0.5 * cell->face (face)->line (line)->diameter () / fe.degree; + const std::vector > & + reference_quadrature_points = fe_values.get_quadrature ().get_points (); + + // Project the boundary function onto the shape functions for this edge + // and set up a linear system of equations to get the values for the DoFs + // associated with this edge. + for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points; ++q_point) + { + // Compute the tangential + // of the edge at + // the quadrature point. + Point shifted_reference_point_1 = reference_quadrature_points[q_point]; + Point shifted_reference_point_2 = reference_quadrature_points[q_point]; + + shifted_reference_point_1 (edge_coordinate_direction[face][line]) += tol; + shifted_reference_point_2 (edge_coordinate_direction[face][line]) -= tol; + Tensor<1,dim> tangential + = (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); + tangential + /= tangential.norm (); + + // Compute the entires of the linear system + // Note the system is symmetric so we could only compute the lower/upper triangle. + // + // The matrix entries are + // \int_{edge} (tangential*edge_shape_function_i)*(tangential*edge_shape_function_j) dS + // + // The RHS entries are: + // \int_{edge} (tangential*boundary_value)*(tangential*edge_shape_function_i) dS. + for (unsigned int j = 0; j < associated_edge_dofs; ++j) + { + const unsigned int j_face_idx = associated_edge_dof_to_face_dof[j]; + const unsigned int j_cell_idx = fe.face_to_cell_index (j_face_idx, face); + for (unsigned int i = 0; i < associated_edge_dofs; ++i) + { + const unsigned int i_face_idx = associated_edge_dof_to_face_dof[i]; + const unsigned int i_cell_idx = fe.face_to_cell_index (i_face_idx, face); + + edge_matrix(i,j) + += fe_values.JxW (q_point) + * (fe_values[vec].value (i_cell_idx, q_point) * tangential) + * (fe_values[vec].value (j_cell_idx, q_point) * tangential); + } + // Compute the RHS entries: + edge_rhs(j) += fe_values.JxW (q_point) + * (values[q_point] (first_vector_component) * tangential [0] + + values[q_point] (first_vector_component + 1) * tangential [1] + + values[q_point] (first_vector_component + 2) * tangential [2]) + * (fe_values[vec].value (j_cell_idx, q_point) * tangential); + } + } + + // Invert linear system + edge_matrix_inv.invert(edge_matrix); + edge_matrix_inv.vmult(edge_solution,edge_rhs); + + // Store computed DoFs + for (unsigned int i = 0; i < associated_edge_dofs; ++i) + { + dof_values[associated_edge_dof_to_face_dof[i]] = edge_solution(i); + dofs_processed[associated_edge_dof_to_face_dof[i]] = true; + } + } + + + template + void + compute_edge_projection_l2 (const cell_iterator &, + const unsigned int, + const unsigned int, + hp::FEValues &, + const Function &, + const unsigned int, + std::vector &, + std::vector &) + { + // dummy implementation of above function + // for all other dimensions + Assert (false, ExcInternalError ()); + } + + template + void + compute_face_projection_curl_conforming_l2 (const cell_iterator &cell, + const unsigned int face, + hp::FEFaceValues &hp_fe_face_values, + const Function &boundary_function, + const unsigned int first_vector_component, + std::vector &dof_values, + std::vector &dofs_processed) + { + // This function computes the L2-projection of the boundary + // function on the interior of faces only. In 3D, this should only be + // called after first calling compute_edge_projection_l2, as it relies on + // edge constraints which are found. + + // In the context of this function, by associated DoFs we mean: + // the DoFs corresponding to the group of components making up the vector + // with first component first_vector_component (with total components dim). + + // Copy to the standard FEFaceValues object: + hp_fe_face_values.reinit (cell, face); + const FEFaceValues &fe_face_values = hp_fe_face_values.get_present_fe_values(); + + // Initialize the required objects. + const FiniteElement &fe = cell->get_fe (); + const std::vector > & + quadrature_points = fe_face_values.get_quadrature_points (); + const unsigned int degree = fe.degree - 1; + + std::vector > + values (fe_face_values.n_quadrature_points, Vector (fe.n_components ())); + + // Get boundary function values at quadrature points. + boundary_function.vector_value_list (quadrature_points, values); + + // Find where the group of vector components (dim of them, + // starting at first_vector_component) are within an FESystem. + // + // If not using FESystem then must be using FE_Nedelec, + // which has one base element and one copy of it (with 3 components). + std::pair base_indices (0, 0); + if (dynamic_cast*> (&cell->get_fe ()) != 0) + { + unsigned int fe_index = 0; + unsigned int fe_index_old = 0; + unsigned int i = 0; + + // Find base element: + // base_indices.first + // + // Then select which copy of that base element + // [ each copy is of length fe.base_element(base_indices.first).n_components() ] + // corresponds to first_vector_component: + // base_index.second + 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 (); + } + + switch (dim) + { + case 2: + // NOTE: This is very similar to compute_edge_projection as used in 3D, + // and contains a lot of overlap with that function. + { + + // Find the DoFs we want to constrain. There are degree+1 in total. + // Create a map from these to the face index + // Note: + // - for a single FE_Nedelec<2> element this is + // simply 0 to fe.dofs_per_face + // - for FESystem<2> this just requires matching the + // base element, fe.system_to_base_index.first.first + // and the copy of the base element we're interested + // in, fe.system_to_base_index.first.second + std::vector associated_edge_dof_to_face_dof (degree + 1); + + unsigned int associated_edge_dof_index = 0; + for (unsigned int face_idx = 0; face_idx < fe.dofs_per_face; ++face_idx) + { + const unsigned int cell_idx = fe.face_to_cell_index (face_idx, face); + if (((dynamic_cast*> (&fe) != 0) + && (fe.system_to_base_index (cell_idx).first == base_indices)) + || (dynamic_cast*> (&fe) != 0)) + { + associated_edge_dof_to_face_dof[associated_edge_dof_index] = face_idx; + ++associated_edge_dof_index; + } + } + // Sanity check: + const unsigned int associated_edge_dofs = associated_edge_dof_index; + Assert (associated_edge_dofs == degree + 1, + ExcMessage ("Error: Unexpected number of 2D edge DoFs")); + + // Matrix and RHS vectors to store: + // We have (degree+1) edge basis functions + FullMatrix edge_matrix (degree + 1,degree + 1); + FullMatrix edge_matrix_inv (degree + 1,degree + 1); + Vector edge_rhs (degree + 1); + Vector edge_solution (degree + 1); + + const FEValuesExtractors::Vector vec (first_vector_component); + + // coordinate directions of the face. + const unsigned int + face_coordinate_direction[GeometryInfo::faces_per_cell] + = { 1, 1, 0, 0 }; + + const double tol = 0.5 * cell->face (face)->diameter () / cell->get_fe ().degree; + const std::vector > & + reference_quadrature_points = fe_face_values.get_quadrature_points (); + + // Project the boundary function onto the shape functions for this edge + // and set up a linear system of equations to get the values for the DoFs + // associated with this edge. + for (unsigned int q_point = 0; + q_point < fe_face_values.n_quadrature_points; ++q_point) + { + // Compute the tangent of the face + // at the quadrature point. + Point shifted_reference_point_1 + = reference_quadrature_points[q_point]; + Point shifted_reference_point_2 + = reference_quadrature_points[q_point]; + + shifted_reference_point_1 (face_coordinate_direction[face]) + += tol; + shifted_reference_point_2 (face_coordinate_direction[face]) + -= tol; + Tensor<1,dim> tangential + = (fe_face_values.get_mapping () + .transform_unit_to_real_cell (cell, + shifted_reference_point_1) + - + fe_face_values.get_mapping () + .transform_unit_to_real_cell (cell, + shifted_reference_point_2)) + / tol; + tangential + /= tangential.norm (); + + // Compute the entires of the linear system + // Note the system is symmetric so we could only compute the lower/upper triangle. + // + // The matrix entries are + // \int_{edge} (tangential * edge_shape_function_i) * (tangential * edge_shape_function_j) dS + // + // The RHS entries are: + // \int_{edge} (tangential* boundary_value) * (tangential * edge_shape_function_i) dS. + for (unsigned int j = 0; j < associated_edge_dofs; ++j) + { + const unsigned int j_face_idx = associated_edge_dof_to_face_dof[j]; + for (unsigned int i = 0; i < associated_edge_dofs; ++i) + { + const unsigned int i_face_idx = associated_edge_dof_to_face_dof[i]; + edge_matrix(i,j) + += fe_face_values.JxW (q_point) + * (fe_face_values[vec].value (fe.face_to_cell_index (i_face_idx, face), q_point) * tangential) + * (fe_face_values[vec].value (fe.face_to_cell_index (j_face_idx, face), q_point) * tangential); + } + edge_rhs(j) + += fe_face_values.JxW (q_point) + * (values[q_point] (first_vector_component) * tangential [0] + + values[q_point] (first_vector_component + 1) * tangential [1]) + * (fe_face_values[vec].value (fe.face_to_cell_index (j_face_idx, face), q_point) * tangential); + } + } + + // Invert linear system + edge_matrix_inv.invert(edge_matrix); + edge_matrix_inv.vmult(edge_solution,edge_rhs); + + // Store computed DoFs + for (unsigned int associated_edge_dof_index = 0; associated_edge_dof_index < associated_edge_dofs; ++associated_edge_dof_index) + { + dof_values[associated_edge_dof_to_face_dof[associated_edge_dof_index]] = edge_solution (associated_edge_dof_index); + dofs_processed[associated_edge_dof_to_face_dof[associated_edge_dof_index]] = true; + } + break; + } + + case 3: + { + const FEValuesExtractors::Vector vec (first_vector_component); + + // First group DoFs associated with edges which we already know. + // Sort these into groups of dofs (0 -> degree+1 of them) by each edge. + // This will help when computing the residual for the face projections. + // + // This matches with the search done in compute_edge_projection. + const unsigned int lines_per_face = GeometryInfo::lines_per_face; + std::vector> + associated_edge_dof_to_face_dof(lines_per_face, std::vector (degree + 1)); + std::vector associated_edge_dofs (lines_per_face); + + for (unsigned int line = 0; line < lines_per_face; ++line) + { + // Lowest DoF in the base element allowed for this edge: + const unsigned int lower_bound = + fe.base_element(base_indices.first).face_to_cell_index(line * (degree + 1), face); + // Highest DoF in the base element allowed for this edge: + const unsigned int upper_bound = + fe.base_element(base_indices.first).face_to_cell_index((line + 1) * (degree + 1) - 1, face); + unsigned int associated_edge_dof_index = 0; + for (unsigned int line_idx = 0; line_idx < fe.dofs_per_line; ++line_idx) + { + const unsigned int face_idx = line*fe.dofs_per_line + line_idx; + const unsigned int cell_idx = fe.face_to_cell_index(face_idx, face); + // Check this cell_idx belongs to the correct base_element, component and line: + if (((dynamic_cast*> (&fe) != 0) + && (fe.system_to_base_index (cell_idx).first == base_indices) + && (lower_bound <= fe.system_to_base_index (cell_idx).second) + && (fe.system_to_base_index (cell_idx).second <= upper_bound)) + || ((dynamic_cast*> (&fe) != 0) + && (line * (degree + 1) <= face_idx) + && (face_idx <= (line + 1) * (degree + 1) - 1))) + { + associated_edge_dof_to_face_dof[line][associated_edge_dof_index] = face_idx; + ++associated_edge_dof_index; + } + } + // Sanity check: + associated_edge_dofs[line] = associated_edge_dof_index; + Assert (associated_edge_dofs[line] == degree + 1, + ExcMessage ("Error: Unexpected number of 3D edge DoFs")); + } + + // Next find the face DoFs associated with the vector components + // we're interested in. There are 2*degree*(degree+1) DoFs associated + // with each face (not including edges!). + // + // Create a map mapping from the consecutively numbered associated_dofs + // to the face DoF (which can be transferred to a local cell index). + // + // For FE_Nedelec<3> we just need to have a face numbering greater than + // the number of edge DoFs (=lines_per_face*(degree+1). + // + // For FE_System<3> we need to base the base_indices (base element and + // copy within that base element) and ensure we're above the number of + // edge DoFs within that base element. + std::vector associated_face_dof_to_face_dof (2*degree*(degree + 1)); + + // Skip the edge DoFs, so we start at lines_per_face*(fe.dofs_per_line). + unsigned int associated_face_dof_index = 0; + for (unsigned int face_idx = lines_per_face*(fe.dofs_per_line); face_idx < fe.dofs_per_face; ++face_idx) + { + const unsigned int cell_idx = fe.face_to_cell_index (face_idx, face); + if (((dynamic_cast*> (&fe) != 0) + && (fe.system_to_base_index (cell_idx).first == base_indices)) + || ((dynamic_cast*> (&fe) != 0))) + { + associated_face_dof_to_face_dof[associated_face_dof_index] = face_idx; + ++associated_face_dof_index; + } + } + // Sanity check: + const unsigned int associated_face_dofs = associated_face_dof_index; + Assert (associated_face_dofs == 2*degree*(degree + 1), + ExcMessage ("Error: Unexpected number of 3D face DoFs")); + + // Storage for the linear system. + // There are 2*degree*(degree+1) DoFs associated with a face in 3D. + // Note this doesn't include the DoFs associated with edges on that face. + FullMatrix face_matrix (2*degree*(degree + 1)); + FullMatrix face_matrix_inv (2*degree*(degree + 1)); + Vector face_rhs (2*degree*(degree + 1)); + Vector face_solution (2*degree*(degree + 1)); + + // Project the boundary function onto the shape functions for this face + // and set up a linear system of equations to get the values for the DoFs + // associated with this face. We also must include the residuals from the + // shape funcations associated with edges. + Tensor<1, dim> tmp; + Tensor<1, dim> normal_vector, + cross_product_i, + cross_product_j, + cross_product_rhs; + + // Store all normal vectors at quad points: + std::vector> normal_vector_list(fe_face_values.get_normal_vectors()); + + // Loop to construct face linear system. + for (unsigned int q_point = 0; + q_point < fe_face_values.n_quadrature_points; ++q_point) + { + // First calculate the residual from the edge functions + // store the result in tmp. + // + // Edge_residual = + // boundary_value - ( + // \sum_(edges on face) + // \sum_(DoFs on edge) edge_dof_value*edge_shape_function + // ) + for (unsigned int d = 0; d < dim; ++d) + { + tmp[d]=0.0; + } + for (unsigned int line = 0; line < lines_per_face; ++line) + { + for (unsigned int associated_edge_dof = 0; associated_edge_dof < associated_edge_dofs[line]; ++associated_edge_dof) + { + const unsigned int face_idx = associated_edge_dof_to_face_dof[line][associated_edge_dof]; + const unsigned int cell_idx = fe.face_to_cell_index(face_idx, face); + tmp -= dof_values[face_idx]*fe_face_values[vec].value(cell_idx, q_point); + } + } + + for (unsigned int d=0; d + void + compute_project_boundary_values_curl_conforming_l2 (const DH &dof_handler, + const unsigned int first_vector_component, + const Function &boundary_function, + const types::boundary_id boundary_component, + ConstraintMatrix &constraints, + const hp::MappingCollection &mapping_collection) + { + // L2-projection based interpolation formed in one (in 2D) or two (in 3D) steps. + // + // In 2D we only need to constrain edge DoFs. + // + // In 3D we need to constrain both edge and face DoFs. This is done in two parts. + // + // For edges, since the face shape functions are zero here ("bubble functions"), + // we project the tangential component of the boundary function and compute + // the L2-projection. This returns the values for the DoFs associated with each edge shape + // function. In 3D, this is computed by internals::compute_edge_projection_l2, in 2D, + // it is handled by compute_face_projection_curl_conforming_l2. + // + // For faces we compute the residual of the boundary function which is satisfied + // by the edge shape functions alone. Which can then be used to calculate the + // remaining face DoF values via a projection which leads to a linear system to + // solve. This is handled by compute_face_projection_curl_conforming_l2 + // + // For details see (for example) section 4.2: + // Electromagnetic scattering simulation using an H (curl) conforming hp finite element + // method in three dimensions, PD Ledger, K Morgan, O Hassan, + // Int. J. Num. Meth. Fluids, Volume 53, Issue 8, pages 1267–1296, 20 March 2007: + // http://onlinelibrary.wiley.com/doi/10.1002/fld.1223/abstract + + // Create hp FEcollection, dof_handler can be either hp or standard type. + // From here on we can treat it like a hp-namespace object. + const hp::FECollection fe_collection (dof_handler.get_fe ()); + + // Create face quadrature collection + hp::QCollection face_quadrature_collection; + for (unsigned int i = 0; i < fe_collection.size (); ++i) + { + const QGauss reference_face_quadrature (2 * fe_collection[i].degree + 1); + face_quadrature_collection.push_back(reference_face_quadrature); + } + + hp::FEFaceValues fe_face_values (mapping_collection, fe_collection, + face_quadrature_collection, + update_values | + update_quadrature_points | + update_normal_vectors | + update_JxW_values); + + // Storage for dof values found and whether they have been processed: + std::vector dofs_processed; + std::vector dof_values; + std::vector face_dof_indices; + typename DH::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 (unsigned int face = 0; face < GeometryInfo::faces_per_cell; ++face) + { + if (cell->face (face)->boundary_indicator () == boundary_component) + { + // If the FE is an FE_Nothing object there is no work to do + if (dynamic_cast*> (&cell->get_fe ()) != 0) + { + return; + } + + // This is only implemented for FE_Nedelec elements. + // If the FE is a FESystem we cannot check this. + if (dynamic_cast*> (&cell->get_fe ()) == 0) + { + typedef FiniteElement FEL; + AssertThrow (dynamic_cast*> (&cell->get_fe ()) != 0, + typename FEL::ExcInterpolationNotImplemented ()); + + } + + const unsigned int dofs_per_face = cell->get_fe ().dofs_per_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; + } + + // Compute the projection of the boundary function on the edge. + // In 2D this is all that's required. + compute_face_projection_curl_conforming_l2 (cell, face, fe_face_values, + boundary_function, + first_vector_component, + dof_values, dofs_processed); + + // store the local->global map: + 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, + // assuming 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: + { + hp::QCollection edge_quadrature_collection; + + // Create equivalent of FEEdgeValues: + for (unsigned int i = 0; i < fe_collection.size (); ++i) + { + const QGauss reference_edge_quadrature (2 * fe_collection[i].degree + 1); + for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; ++face) + { + for (unsigned int line = 0; line < GeometryInfo::lines_per_face; ++line) + { + edge_quadrature_collection.push_back + (QProjector::project_to_face + (QProjector::project_to_face + (reference_edge_quadrature, line), face)); + } + } + } + + hp::FEValues fe_edge_values (mapping_collection, fe_collection, + edge_quadrature_collection, + update_jacobians | + update_JxW_values | + update_quadrature_points | + update_values); + + for (; cell != dof_handler.end (); ++cell) + { + if (cell->at_boundary () && cell->is_locally_owned ()) + { + for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; ++face) + { + if (cell->face (face)->boundary_indicator () == boundary_component) + { + // If the FE is an FE_Nothing object there is no work to do + if (dynamic_cast*> (&cell->get_fe ()) != 0) + { + return; + } + + // This is only implemented for FE_Nedelec elements. + // If the FE is a FESystem we cannot check this. + if (dynamic_cast*> (&cell->get_fe ()) == 0) + { + typedef FiniteElement FEL; + + AssertThrow (dynamic_cast*> (&cell->get_fe ()) != 0, + typename FEL::ExcInterpolationNotImplemented ()); + } + + const unsigned int superdegree = cell->get_fe ().degree; + const unsigned int degree = superdegree - 1; + const unsigned int dofs_per_face = cell->get_fe ().dofs_per_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; + } + + // First compute the projection on the edges. + for (unsigned int line = 0; line < GeometryInfo<3>::lines_per_face; ++line) + { + compute_edge_projection_l2 (cell, face, line, fe_edge_values, + boundary_function, + first_vector_component, + dof_values, dofs_processed); + } + + // If there are higher order shape functions, then we + // still need to compute the face projection + if (degree > 0) + { + compute_face_projection_curl_conforming_l2 (cell, face, fe_face_values, + boundary_function, + first_vector_component, + dof_values, + dofs_processed); + } + + // Store the computed values in the global vector. + 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 ()); + } + } + + } + + + template + void + project_boundary_values_curl_conforming_l2 (const DoFHandler &dof_handler, + const unsigned int first_vector_component, + const Function &boundary_function, + const types::boundary_id boundary_component, + ConstraintMatrix &constraints, + const Mapping &mapping) + { + // non-hp version - calls the internal + // compute_project_boundary_values_curl_conforming_l2() function + // above after recasting the mapping. + + hp::MappingCollection mapping_collection (mapping); + internals:: + compute_project_boundary_values_curl_conforming_l2(dof_handler, + first_vector_component, + boundary_function, + boundary_component, + constraints, + mapping_collection); + } + + template + void + project_boundary_values_curl_conforming_l2 (const hp::DoFHandler &dof_handler, + const unsigned int first_vector_component, + const Function &boundary_function, + const types::boundary_id boundary_component, + ConstraintMatrix &constraints, + const hp::MappingCollection &mapping_collection) + { + // hp version - calls the internal + // compute_project_boundary_values_curl_conforming_l2() function above. + internals:: + compute_project_boundary_values_curl_conforming_l2(dof_handler, + first_vector_component, + boundary_function, + boundary_component, + constraints, + mapping_collection); + } + + + namespace internals { // This function computes the projection of the boundary function on the diff --git a/source/numerics/vector_tools_boundary.inst.in b/source/numerics/vector_tools_boundary.inst.in index 2b006e4727..dec32b7e2a 100644 --- a/source/numerics/vector_tools_boundary.inst.in +++ b/source/numerics/vector_tools_boundary.inst.in @@ -141,6 +141,22 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) const types::boundary_id, ConstraintMatrix&, const hp::MappingCollection&); + template + void project_boundary_values_curl_conforming_l2 + (const DoFHandler&, + const unsigned int, + const Function&, + const types::boundary_id, + ConstraintMatrix&, + const Mapping&); + template + void project_boundary_values_curl_conforming_l2 + (const hp::DoFHandler&, + const unsigned int, + const Function&, + const types::boundary_id, + ConstraintMatrix&, + const hp::MappingCollection&); template void project_boundary_values_div_conforming (const DoFHandler&, diff --git a/tests/deal.II/nedelec_non_rect_face.cc b/tests/deal.II/nedelec_non_rect_face.cc new file mode 100644 index 0000000000..283ad77f4f --- /dev/null +++ b/tests/deal.II/nedelec_non_rect_face.cc @@ -0,0 +1,472 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 1998 - 2014 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +// +// By Ross Kynch +// +// Test to confirm that FE_Nedelec works on non-rectangular faces in 3D. +// The previous project_boundary_values_curl_conforming function only +// produced the correct constraints on rectangular faces. +// +// The new function project_boundary_values_curl_conforming_l2 fixes this +// and should work on any quad face, assuming the underlying mesh conforms to the +// standard orientation (i.e. spheres may still cause problems). +// +// This test solves the real valued curl-curl equation in 3D: +// +// curl(curl(E)) + E = Js +// +// where the solution is: +// E = (x^2 , y^2 + , z^2). +// +// so: +// Js = E. +// +// The domain is a distorted cube which will contain non-rectangular faces. + +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include + +#include +#include +#include + +namespace Maxwell +{ + using namespace dealii; + + // Dirichlet BCs / exact solution:. + template + class ExactSolution : public Function + { + public: + ExactSolution(); + + virtual void vector_value_list (const std::vector > &points, + std::vector > &values) const; + + + void curl_value_list(const std::vector > &points, + std::vector > &value_list); + }; + + template + ExactSolution::ExactSolution() + : + Function (dim) + {} + template + void ExactSolution::vector_value_list (const std::vector > &points, + std::vector > &value_list) const + { + Assert(value_list.size() == points.size(), ExcDimensionMismatch(value_list.size(), points.size())); + const unsigned int n_points = points.size(); + + for (unsigned int i=0; i &p = points[i]; + + /* quadratic: */ + value_list[i](0) = p(0)*p(0); + value_list[i](1) = p(1)*p(1); + value_list[i](2) = p(2)*p(2); + + } + } + // Additional functions to create Neumann conditions, zero in this case. + template + void ExactSolution::curl_value_list(const std::vector > &points, + std::vector > &value_list) + { + Assert(value_list.size() == points.size(), ExcDimensionMismatch(value_list.size(), points.size())); + const unsigned int n_points = points.size(); + + + double exponent; + for (unsigned int i=0; i &p = points[i]; + // Real: + value_list[i](0) = 0.0; + value_list[i](1) = 0.0; + value_list[i](2) = 0.0; + } + } + // END Dirichlet BC + + // MAIN MAXWELL CLASS + template + class MaxwellProblem + { + public: + MaxwellProblem (const unsigned int order); + ~MaxwellProblem (); + void run (); + private: + void setup_system (); + void assemble_system (); + void solve (); + void process_solution(const unsigned int cycle); + + double calcErrorHcurlNorm(); + + Triangulation triangulation; + MappingQ mapping; + DoFHandler dof_handler; + FE_Nedelec fe; + ConstraintMatrix constraints; + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + Vector solution; + Vector system_rhs; + + ConvergenceTable convergence_table; + + ExactSolution exact_solution; + + unsigned int p_order; + unsigned int quad_order; + }; + + template + MaxwellProblem::MaxwellProblem (const unsigned int order) + : + mapping (1), + dof_handler (triangulation), + fe(order), + exact_solution() + { + p_order = order; + quad_order = 2*(p_order+1); + } + + template + MaxwellProblem::~MaxwellProblem () + { + dof_handler.clear (); + } + template + double MaxwellProblem::calcErrorHcurlNorm() + { + QGauss quadrature_formula(quad_order); + const unsigned int n_q_points = quadrature_formula.size(); + + FEValues fe_values (mapping, fe, quadrature_formula, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); + + // Extractor + const FEValuesExtractors::Vector E_re(0); + + const unsigned int dofs_per_cell = fe.dofs_per_cell; + + std::vector local_dof_indices (dofs_per_cell); + + // storage for exact sol: + std::vector > exactsol_list(n_q_points, Vector(fe.n_components())); + std::vector > exactcurlsol_list(n_q_points, Vector(fe.n_components())); + Tensor<1,dim> exactsol; + Tensor<1,dim> exactcurlsol; + + // storage for computed sol: + std::vector > sol(n_q_points); + Tensor<1,dim> curlsol(n_q_points); + + double h_curl_norm=0.0; + + unsigned int block_index_i; + + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell) + { + fe_values.reinit (cell); + + // Store exact values of E and curlE: + exact_solution.vector_value_list(fe_values.get_quadrature_points(), exactsol_list); + exact_solution.curl_value_list(fe_values.get_quadrature_points(), exactcurlsol_list); + + + // Store computed values at quad points: + fe_values[E_re].get_function_values(solution, sol); + + // Calc values of curlE from fe solution: + cell->get_dof_indices (local_dof_indices); + // Loop over quad points to calculate solution: + for (unsigned int q_point=0; q_point + void MaxwellProblem::setup_system () + { + dof_handler.distribute_dofs (fe); + DoFRenumbering::block_wise (dof_handler); + solution.reinit (dof_handler.n_dofs()); + system_rhs.reinit (dof_handler.n_dofs()); + constraints.clear (); + + + DoFTools::make_hanging_node_constraints (dof_handler, constraints); + + // FE_Nedelec boundary condition. + VectorTools::project_boundary_values_curl_conforming_l2 (dof_handler, 0, exact_solution, 0, constraints); + + + constraints.close (); + + CompressedSparsityPattern c_sparsity(dof_handler.n_dofs()); + DoFTools::make_sparsity_pattern(dof_handler, + c_sparsity, + constraints,false); + + sparsity_pattern.copy_from(c_sparsity); + system_matrix.reinit (sparsity_pattern); + + } + template + void MaxwellProblem::assemble_system () + { + QGauss quadrature_formula(quad_order); + QGauss face_quadrature_formula(quad_order); + + const unsigned int n_q_points = quadrature_formula.size(); + const unsigned int n_face_q_points = face_quadrature_formula.size(); + + const unsigned int dofs_per_cell = fe.dofs_per_cell; + + FEValues fe_values (mapping, fe, quadrature_formula, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); + + FEFaceValues fe_face_values(mapping, fe, face_quadrature_formula, + update_values | update_quadrature_points | + update_normal_vectors | update_JxW_values); + + // Extractor + const FEValuesExtractors::Vector E_re(0); + + // Local cell storage: + FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); + Vector cell_rhs (dofs_per_cell); + std::vector local_dof_indices (dofs_per_cell); + + //RHS storage: + std::vector > rhs_value_list(n_q_points, Vector(fe.n_components())); + Tensor<1,dim> rhs_value_vector(dim); + + // Neumann storage + std::vector > neumann_value_list(n_face_q_points, Vector(fe.n_components())); + std::vector> normal_vector_list(fe_face_values.get_normal_vectors()); + Tensor<1,dim> neumann_value_vector(dim); + Tensor<1,dim> neumann_value(dim); + Tensor<1,dim> normal_vector; + + // loop over all cells: + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell) + { + fe_values.reinit(cell); + cell_matrix = 0; + cell_rhs = 0; + + // Calc RHS values: + exact_solution.vector_value_list(fe_values.get_quadrature_points(), rhs_value_list); + + // Loop over all element quad points: + for (unsigned int q_point=0; q_pointget_dof_indices (local_dof_indices); + constraints.distribute_local_to_global(cell_matrix, + cell_rhs, + local_dof_indices, + system_matrix, system_rhs); + } + } + template + void MaxwellProblem::solve () + { + /* Direct */ + SparseDirectUMFPACK A_direct; + A_direct.initialize(system_matrix); + + A_direct.vmult (solution, system_rhs); + constraints.distribute (solution); + + } + template + void MaxwellProblem::process_solution(const unsigned int cycle) + { + + Vector diff_per_cell(triangulation.n_active_cells()); + + + VectorTools::integrate_difference(mapping, dof_handler, solution, exact_solution, + diff_per_cell, QGauss(quad_order+2), + VectorTools::L2_norm); + + + const double L2_error = diff_per_cell.l2_norm(); + + double Hcurl_error = calcErrorHcurlNorm(); + + convergence_table.add_value("cycle", cycle); + convergence_table.add_value("cells", triangulation.n_active_cells()); + convergence_table.add_value("dofs", dof_handler.n_dofs()); + convergence_table.add_value("L2 Error", L2_error); + convergence_table.add_value("H(curl) Error", Hcurl_error); + + } + + template + void MaxwellProblem::run () + { + + unsigned int numcycles=3; + + for (unsigned int cycle=0; cycleset_refine_flag(); // create single refined cell + * triangulation.execute_coarsening_and_refinement(); + */ + + } + else + { + triangulation.refine_global (1); + } + + setup_system (); + assemble_system (); + solve (); + process_solution (cycle); + } + + deallog << "p = " << p_order << std::endl; + + //convergence_table.set_precision("L2 Error",4); + convergence_table.set_scientific("L2 Error",true); + + //convergence_table.set_precision("H(curl) Error",4); + convergence_table.set_scientific("H(curl) Error",true); + + convergence_table.write_text(deallog.get_file_stream()); + + } +} +// END MAXWELL CLASS +int main () +{ + using namespace Maxwell; + + std::ofstream logfile("output"); + deallog.attach(logfile); + deallog.depth_console(0); + + for (unsigned int p=0; p<4; ++p) + { + MaxwellProblem<3> (p).run(); + } + + return 0; +} diff --git a/tests/deal.II/nedelec_non_rect_face.output b/tests/deal.II/nedelec_non_rect_face.output new file mode 100644 index 0000000000..aee54e35d4 --- /dev/null +++ b/tests/deal.II/nedelec_non_rect_face.output @@ -0,0 +1,21 @@ + +DEAL::p = 0 +cycle cells dofs L2 Error H(curl) Error +0 1 12 6.2187e-01 6.0931e-01 +1 8 54 3.3123e-01 3.2983e-01 +2 64 300 1.6777e-01 1.6760e-01 +DEAL::p = 1 +cycle cells dofs L2 Error H(curl) Error +0 1 54 1.5549e-01 1.6230e-01 +1 8 300 3.8832e-02 4.1960e-02 +2 64 1944 9.7185e-03 1.1143e-02 +DEAL::p = 2 +cycle cells dofs L2 Error H(curl) Error +0 1 144 6.0548e-16 1.2581e-15 +1 8 882 4.7208e-15 5.6992e-15 +2 64 6084 3.7516e-14 4.1547e-14 +DEAL::p = 3 +cycle cells dofs L2 Error H(curl) Error +0 1 300 9.4037e-16 1.8758e-15 +1 8 1944 1.5731e-14 1.6809e-14 +2 64 13872 8.7253e-14 9.1275e-14