}
+ namespace internals
+ {
+ template<typename cell_iterator>
+ 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<double> &dof_values,
+ std::vector<bool> &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<dim> &fe = cell->get_fe ();
+
+ // reinit for this cell, face and line.
+ hp_fe_values.reinit
+ (cell,
+ (cell->active_fe_index () * GeometryInfo<dim>::faces_per_cell + face)
+ * GeometryInfo<dim>::lines_per_face + line);
+
+ // Initialize the required objects.
+ const FEValues<dim> &
+ fe_values = hp_fe_values.get_present_fe_values ();
+ // Store degree as fe.degree-1
+ // For nedelec elements FE_Nedelec<dim> (0) returns fe.degree = 1.
+ const unsigned int degree = fe.degree - 1;
+
+ const std::vector<Point<dim> > &
+ quadrature_points = fe_values.get_quadrature_points ();
+ std::vector<Vector<double> > values (fe_values.n_quadrature_points,
+ Vector<double> (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<unsigned int, unsigned int> base_indices (0, 0);
+ if (dynamic_cast<const FESystem<dim>*> (&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<unsigned int> 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<const FESystem<dim>*> (&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<const FE_Nedelec<dim>*> (&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<double> edge_matrix (degree + 1,degree + 1);
+ FullMatrix<double> edge_matrix_inv (degree + 1,degree + 1);
+ Vector<double> edge_rhs (degree + 1);
+ Vector<double> 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<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 double tol = 0.5 * cell->face (face)->line (line)->diameter () / fe.degree;
+ const std::vector<Point<dim> > &
+ 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<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;
+ 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<int dim, typename cell_iterator>
+ void
+ compute_edge_projection_l2 (const cell_iterator &,
+ const unsigned int,
+ const unsigned int,
+ hp::FEValues<dim> &,
+ const Function<dim> &,
+ const unsigned int,
+ std::vector<double> &,
+ std::vector<bool> &)
+ {
+ // dummy implementation of above function
+ // for all other dimensions
+ Assert (false, ExcInternalError ());
+ }
+
+ template<int dim, typename cell_iterator>
+ void
+ compute_face_projection_curl_conforming_l2 (const cell_iterator &cell,
+ const unsigned int face,
+ hp::FEFaceValues<dim> &hp_fe_face_values,
+ const Function<dim> &boundary_function,
+ const unsigned int first_vector_component,
+ std::vector<double> &dof_values,
+ std::vector<bool> &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<dim> &fe_face_values = hp_fe_face_values.get_present_fe_values();
+
+ // Initialize the required objects.
+ const FiniteElement<dim> &fe = cell->get_fe ();
+ const std::vector<Point<dim> > &
+ quadrature_points = fe_face_values.get_quadrature_points ();
+ const unsigned int degree = fe.degree - 1;
+
+ std::vector<Vector<double> >
+ values (fe_face_values.n_quadrature_points, Vector<double> (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<unsigned int, unsigned int> base_indices (0, 0);
+ if (dynamic_cast<const FESystem<dim>*> (&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<unsigned int> 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<const FESystem<dim>*> (&fe) != 0)
+ && (fe.system_to_base_index (cell_idx).first == base_indices))
+ || (dynamic_cast<const FE_Nedelec<dim>*> (&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<double> edge_matrix (degree + 1,degree + 1);
+ FullMatrix<double> edge_matrix_inv (degree + 1,degree + 1);
+ Vector<double> edge_rhs (degree + 1);
+ Vector<double> edge_solution (degree + 1);
+
+ const FEValuesExtractors::Vector vec (first_vector_component);
+
+ // coordinate directions of the face.
+ const unsigned int
+ face_coordinate_direction[GeometryInfo<dim>::faces_per_cell]
+ = { 1, 1, 0, 0 };
+
+ const double tol = 0.5 * cell->face (face)->diameter () / cell->get_fe ().degree;
+ const std::vector<Point<dim> > &
+ 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<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;
+ 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<dim>::lines_per_face;
+ std::vector<std::vector<unsigned int>>
+ associated_edge_dof_to_face_dof(lines_per_face, std::vector<unsigned int> (degree + 1));
+ std::vector<unsigned int> 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<const FESystem<dim>*> (&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<const FE_Nedelec<dim>*> (&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<unsigned int> 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<const FESystem<dim>*> (&fe) != 0)
+ && (fe.system_to_base_index (cell_idx).first == base_indices))
+ || ((dynamic_cast<const FE_Nedelec<dim>*> (&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<double> face_matrix (2*degree*(degree + 1));
+ FullMatrix<double> face_matrix_inv (2*degree*(degree + 1));
+ Vector<double> face_rhs (2*degree*(degree + 1));
+ Vector<double> 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<Point<dim>> 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<dim; ++d)
+ {
+ tmp[d] += values[q_point] (first_vector_component + d);
+ }
+
+ // Tensor of normal vector on the face at q_point;
+ for (unsigned int d = 0; d < dim; ++d)
+ {
+ normal_vector[d] = normal_vector_list[q_point](d);
+ }
+ // Now compute the linear system:
+ // On a face:
+ // The matrix entries are:
+ // \int_{face} (n x face_shape_function_i) \cdot ( n x face_shape_function_j) dS
+ //
+ // The RHS entries are:
+ // \int_{face} (n x (Edge_residual) \cdot (n x face_shape_function_i) dS
+
+ for (unsigned int j = 0; j < associated_face_dofs; ++j)
+ {
+ const unsigned int j_face_idx = associated_face_dof_to_face_dof[j];
+ const unsigned int cell_j = fe.face_to_cell_index (j_face_idx, face);
+
+ cross_product(cross_product_j,
+ normal_vector,
+ fe_face_values[vec].value(cell_j, q_point));
+
+ for (unsigned int i = 0; i < associated_face_dofs; ++i)
+ {
+ const unsigned int i_face_idx = associated_face_dof_to_face_dof[i];
+ const unsigned int cell_i = fe.face_to_cell_index (i_face_idx, face);
+ cross_product(cross_product_i,
+ normal_vector,
+ fe_face_values[vec].value(cell_i, q_point));
+
+ face_matrix(i,j)
+ += fe_face_values.JxW (q_point)
+ *cross_product_i
+ *cross_product_j;
+
+ }
+ // compute rhs
+ cross_product(cross_product_rhs,
+ normal_vector,
+ tmp);
+ face_rhs(j)
+ += fe_face_values.JxW (q_point)
+ *cross_product_rhs
+ *cross_product_j;
+
+ }
+ }
+
+ // Solve lienar system:
+ face_matrix_inv.invert(face_matrix);
+ face_matrix_inv.vmult(face_solution, face_rhs);
+
+
+ // Store computed DoFs:
+ for (unsigned int associated_face_dof = 0; associated_face_dof < associated_face_dofs; ++associated_face_dof)
+ {
+ dof_values[associated_face_dof_to_face_dof[associated_face_dof]] = face_solution(associated_face_dof);
+ dofs_processed[associated_face_dof_to_face_dof[associated_face_dof]] = true;
+ }
+ break;
+ }
+ default:
+ Assert (false, ExcNotImplemented ());
+ }
+ }
+
+ template <int dim, class DH>
+ void
+ compute_project_boundary_values_curl_conforming_l2 (const DH &dof_handler,
+ const unsigned int first_vector_component,
+ const Function<dim> &boundary_function,
+ const types::boundary_id boundary_component,
+ ConstraintMatrix &constraints,
+ const hp::MappingCollection<dim, dim> &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<dim> fe_collection (dof_handler.get_fe ());
+
+ // Create face quadrature collection
+ hp::QCollection<dim - 1> 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 + 1);
+ face_quadrature_collection.push_back(reference_face_quadrature);
+ }
+
+ hp::FEFaceValues<dim> 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<bool> dofs_processed;
+ std::vector<double> dof_values;
+ std::vector<types::global_dof_index> 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<dim>::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<const FE_Nothing<dim>*> (&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<const FESystem<dim>*> (&cell->get_fe ()) == 0)
+ {
+ typedef FiniteElement<dim> FEL;
+ AssertThrow (dynamic_cast<const FE_Nedelec<dim>*> (&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<dim> edge_quadrature_collection;
+
+ // Create equivalent of FEEdgeValues:
+ for (unsigned int i = 0; i < fe_collection.size (); ++i)
+ {
+ const QGauss<dim-2> reference_edge_quadrature (2 * fe_collection[i].degree + 1);
+ for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
+ {
+ 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 (unsigned int face = 0; face < GeometryInfo<dim>::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<const FE_Nothing<dim>*> (&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<const FESystem<dim>*> (&cell->get_fe ()) == 0)
+ {
+ typedef FiniteElement<dim> FEL;
+
+ AssertThrow (dynamic_cast<const FE_Nedelec<dim>*> (&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 <int dim>
+ void
+ project_boundary_values_curl_conforming_l2 (const DoFHandler<dim> &dof_handler,
+ const unsigned int first_vector_component,
+ const Function<dim> &boundary_function,
+ const types::boundary_id boundary_component,
+ ConstraintMatrix &constraints,
+ const Mapping<dim> &mapping)
+ {
+ // non-hp version - calls the internal
+ // compute_project_boundary_values_curl_conforming_l2() function
+ // above after recasting the mapping.
+
+ hp::MappingCollection<dim> mapping_collection (mapping);
+ internals::
+ compute_project_boundary_values_curl_conforming_l2(dof_handler,
+ first_vector_component,
+ boundary_function,
+ boundary_component,
+ constraints,
+ mapping_collection);
+ }
+
+ template <int dim>
+ void
+ project_boundary_values_curl_conforming_l2 (const hp::DoFHandler<dim> &dof_handler,
+ const unsigned int first_vector_component,
+ const Function<dim> &boundary_function,
+ const types::boundary_id boundary_component,
+ ConstraintMatrix &constraints,
+ const hp::MappingCollection<dim, dim> &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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/convergence_table.h>
+
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/sparse_direct.h>
+#include <deal.II/lac/compressed_sparsity_pattern.h>
+#include <deal.II/lac/constraint_matrix.h>
+
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_in.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/dofs/dof_renumbering.h>
+
+#include <deal.II/fe/fe_nedelec.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/fe_nothing.h>
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/fe/mapping_q1.h>
+
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/numerics/matrix_tools.h>
+
+#include <fstream>
+#include <iostream>
+#include <sstream>
+
+namespace Maxwell
+{
+ using namespace dealii;
+
+ // Dirichlet BCs / exact solution:.
+ template<int dim>
+ class ExactSolution : public Function<dim>
+ {
+ public:
+ ExactSolution();
+
+ virtual void vector_value_list (const std::vector<Point<dim> > &points,
+ std::vector<Vector<double> > &values) const;
+
+
+ void curl_value_list(const std::vector<Point<dim> > &points,
+ std::vector<Vector<double> > &value_list);
+ };
+
+ template<int dim>
+ ExactSolution<dim>::ExactSolution()
+ :
+ Function<dim> (dim)
+ {}
+ template <int dim>
+ void ExactSolution<dim>::vector_value_list (const std::vector<Point<dim> > &points,
+ std::vector<Vector<double> > &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<n_points; ++i)
+ {
+ const Point<dim> &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 <int dim>
+ void ExactSolution<dim>::curl_value_list(const std::vector<Point<dim> > &points,
+ std::vector<Vector<double> > &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<n_points; ++i)
+ {
+ const Point<dim> &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 <int dim>
+ 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<dim> triangulation;
+ MappingQ<dim> mapping;
+ DoFHandler<dim> dof_handler;
+ FE_Nedelec<dim> fe;
+ ConstraintMatrix constraints;
+ SparsityPattern sparsity_pattern;
+ SparseMatrix<double> system_matrix;
+ Vector<double> solution;
+ Vector<double> system_rhs;
+
+ ConvergenceTable convergence_table;
+
+ ExactSolution<dim> exact_solution;
+
+ unsigned int p_order;
+ unsigned int quad_order;
+ };
+
+ template <int dim>
+ MaxwellProblem<dim>::MaxwellProblem (const unsigned int order)
+ :
+ mapping (1),
+ dof_handler (triangulation),
+ fe(order),
+ exact_solution()
+ {
+ p_order = order;
+ quad_order = 2*(p_order+1);
+ }
+
+ template <int dim>
+ MaxwellProblem<dim>::~MaxwellProblem ()
+ {
+ dof_handler.clear ();
+ }
+ template<int dim>
+ double MaxwellProblem<dim>::calcErrorHcurlNorm()
+ {
+ QGauss<dim> quadrature_formula(quad_order);
+ const unsigned int n_q_points = quadrature_formula.size();
+
+ FEValues<dim> 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<types::global_dof_index> local_dof_indices (dofs_per_cell);
+
+ // storage for exact sol:
+ std::vector<Vector<double> > exactsol_list(n_q_points, Vector<double>(fe.n_components()));
+ std::vector<Vector<double> > exactcurlsol_list(n_q_points, Vector<double>(fe.n_components()));
+ Tensor<1,dim> exactsol;
+ Tensor<1,dim> exactcurlsol;
+
+ // storage for computed sol:
+ std::vector<Tensor<1,dim> > sol(n_q_points);
+ Tensor<1,dim> curlsol(n_q_points);
+
+ double h_curl_norm=0.0;
+
+ unsigned int block_index_i;
+
+ typename DoFHandler<dim>::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<n_q_points; ++q_point)
+ {
+ // Split exact solution into real/imaginary parts:
+ for (unsigned int component=0; component<dim; component++)
+ {
+ exactsol[component] = exactsol_list[q_point][component];
+ exactcurlsol[component] = exactcurlsol_list[q_point][component];
+ }
+ // Loop over DoFs to calculate curl of solution @ quad point
+ curlsol=0.0;
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ {
+ // Construct local curl value @ quad point
+ curlsol += solution(local_dof_indices[i])*fe_values[E_re].curl(i,q_point);
+ }
+ // Integrate difference at each point:
+ h_curl_norm += ( (exactsol-sol[q_point])*(exactsol-sol[q_point])
+ + (exactcurlsol-curlsol)*(exactcurlsol-curlsol)
+ )*fe_values.JxW(q_point);
+ }
+ }
+ return sqrt(h_curl_norm);
+ }
+
+ template <int dim>
+ void MaxwellProblem<dim>::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 <int dim>
+ void MaxwellProblem<dim>::assemble_system ()
+ {
+ QGauss<dim> quadrature_formula(quad_order);
+ QGauss<dim-1> 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<dim> fe_values (mapping, fe, quadrature_formula,
+ update_values | update_gradients |
+ update_quadrature_points | update_JxW_values);
+
+ FEFaceValues<dim> 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<double> cell_matrix (dofs_per_cell, dofs_per_cell);
+ Vector<double> cell_rhs (dofs_per_cell);
+ std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
+
+ //RHS storage:
+ std::vector<Vector<double> > rhs_value_list(n_q_points, Vector<double>(fe.n_components()));
+ Tensor<1,dim> rhs_value_vector(dim);
+
+ // Neumann storage
+ std::vector<Vector<double> > neumann_value_list(n_face_q_points, Vector<double>(fe.n_components()));
+ std::vector<Point<dim>> 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<dim>::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_point<n_q_points; ++q_point)
+ {
+ // store rhs value at this q point & turn into tensor
+ for (unsigned int component=0; component<dim; component++)
+ {
+ rhs_value_vector[component] = rhs_value_list[q_point](component);
+ }
+
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ {
+ // Construct local matrix:
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ {
+ cell_matrix(i,j) += ( (fe_values[E_re].curl(i,q_point)*fe_values[E_re].curl(j,q_point))
+ + fe_values[E_re].value(i,q_point)*fe_values[E_re].value(j,q_point)
+ )*fe_values.JxW(q_point);
+
+ }
+ // construct local RHS:
+ cell_rhs(i) += rhs_value_vector*fe_values[E_re].value(i,q_point)*fe_values.JxW(q_point);
+
+ }
+ }
+
+ cell->get_dof_indices (local_dof_indices);
+ constraints.distribute_local_to_global(cell_matrix,
+ cell_rhs,
+ local_dof_indices,
+ system_matrix, system_rhs);
+ }
+ }
+ template <int dim>
+ void MaxwellProblem<dim>::solve ()
+ {
+ /* Direct */
+ SparseDirectUMFPACK A_direct;
+ A_direct.initialize(system_matrix);
+
+ A_direct.vmult (solution, system_rhs);
+ constraints.distribute (solution);
+
+ }
+ template<int dim>
+ void MaxwellProblem<dim>::process_solution(const unsigned int cycle)
+ {
+
+ Vector<double> diff_per_cell(triangulation.n_active_cells());
+
+
+ VectorTools::integrate_difference(mapping, dof_handler, solution, exact_solution,
+ diff_per_cell, QGauss<dim>(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 <int dim>
+ void MaxwellProblem<dim>::run ()
+ {
+
+ unsigned int numcycles=3;
+
+ for (unsigned int cycle=0; cycle<numcycles; ++cycle)
+ {
+ if (cycle == 0)
+ {
+ /* Cube mesh */
+ GridGenerator::hyper_cube (triangulation, 0, 1);
+ GridTools::distort_random (0.2, triangulation, false);
+
+
+ /* Testing hanging nodes
+ * triangulation.refine_global(1);
+ * cell = triangulation.begin_active();
+ * cell->set_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;
+}
projects / dealii.git / commitdiff