From: Wolfgang Bangerth Date: Fri, 30 Sep 2016 22:25:44 +0000 (-0600) Subject: Update the P1NC element implementation. X-Git-Tag: v8.5.0-rc1~611^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=aceae3e58f0eec9332b098c021f6d37aea0f7548;p=dealii.git Update the P1NC element implementation. Specifically, this updates the class documentation for language and clarity. It also tightens up the code in a number of places. --- diff --git a/include/deal.II/fe/fe_p1nc.h b/include/deal.II/fe/fe_p1nc.h index a920ae647e..b6fa0008c2 100644 --- a/include/deal.II/fe/fe_p1nc.h +++ b/include/deal.II/fe/fe_p1nc.h @@ -29,25 +29,27 @@ DEAL_II_NAMESPACE_OPEN /*@{*/ /** - * Implementation for the scalar version of the P1 nonconforming finite + * Implementation of the scalar version of the P1 nonconforming finite * element, a piecewise linear element on quadrilaterals in 2D. - * This implementation is only for 2D and codimension = 0. + * This implementation is only for 2D cells in a 2D space (i.e., codimension 0). * - * Unlike any continuous conforming finite element which belongs to $H^1_0$, - * the P1 nonconforming element does not enforce the continuity across edges. - * But it requires the continuity just in integral sense: + * Unlike the usual continuous, $H^1$ conforming finite elements, + * the P1 nonconforming element does not enforce continuity across edges. + * However, it requires the continuity in an integral sense: * any function in the space should have the same integral values * on two sides of the common edge shared by two adjacent elements. * - * Thus each function in the nonconforming element space can be discontinuous, - * not included in $H^1_0$, as functions in Discontinuous Galerkin (DG) finite - * element spaces. - * Although any function in DG space also has nonconformity, - * it is completely discontinuous across edges without any relation. - * This is a reason why usual weak formulations for DG schemes contain - * additional penalty terms for jump across edges to control discontinuity. - * However nonconforming elements usually do not need additional terms - * in their weak formulations due to the continuity in integral on edges. + * Thus, each function in the nonconforming element space can be + * discontinuous, and consequently not included in $H^1_0$, just like + * the basis functions in Discontinuous Galerkin (DG) finite element + * spaces. On the other hand, basis functions in DG spaces are + * completely discontinuous across edges without any relation between + * the values from both sides. This is a reason why usual weak + * formulations for DG schemes contain additional penalty terms for + * jump across edges to control discontinuity. However, nonconforming + * elements usually do not need additional terms in their weak + * formulations because their integrals along edges are the same from + * both sides, i.e., there is some level of continuity. * *

Dice Rule

* Since any function in the P1 nonconforming space is piecewise linear on each element, @@ -56,35 +58,35 @@ DEAL_II_NAMESPACE_OPEN * the continuity of the midpoint value of each edge in this case. * * Thus for the P1 nonconforming element, the function values at midpoints on edges of a cell are important. - * The first attempt to define (local) degrees of freedom (DOFs) on a quadrilateral + * The first attempt to define (local) degrees of freedom (DoFs) on a quadrilateral * is by using midpoint values of a function. * * However, these 4 functionals are not linearly independent * because a linear function on 2D is uniquely determined by only 3 independent values. - * A simple observation reads that any linear function on a quadrilateral should satisfies the 'dice rule': - * the sum of two function values at two midpoints of the edge pair on opposite - * position is equal to the sum of those of the other edge pair. + * A simple observation reads that any linear function on a quadrilateral should satisfy the 'dice rule': + * the sum of two function values at the midpoints of the edge pair on opposite + * sides of a cell is equal to the sum of those at the midpoints of the other edge pair. * This is called the 'dice rule' because the number of points on opposite sides of a dice always * adds up to the same number as well (in the case of dice, to seven). * * In formulas, the dice rule is written as $\phi(m_0) + \phi(m_1) = \phi(m_2) + \phi(m_3)$ * for all $\phi$ in the function space where $m_j$ is the midpoint of the edge $e_j$. * Here, we assume the standard numbering convention for edges used in deal.II - * and described for class GeometryInfo. + * and described in class GeometryInfo. * - * Conversely if 4 values at midpoints satisfying the dice rule are just given, + * Conversely if 4 values at midpoints satisfying the dice rule are given, * then there always exists the unique linear function which coincides with 4 midpoints values. * * Due to the dice rule, three values at any three midpoints can determine * the last value at the last midpoint. * It means that the number of independent local functionals on a cell is 3, - * and it is same as the dimension of the linear polynomial space on a cell in 2D. + * and this is also the dimension of the linear polynomial space on a cell in 2D. * *

Shape functions

- * Before introduction of the DOFs, we present 4 local shape functions on a cell. + * Before introducing the degrees of freedom, we present 4 local shape functions on a cell. * Due to the dice rule, we need a special construction for shape functions. * Although the following 4 shape functions are not linearly independent within a cell, - * they are helpful to define the global basis functions which are linearly independent on whole domain. + * they are helpful to define the global basis functions which are linearly independent on the whole domain. * Again, we assume the standard numbering for vertices used in deal.II. * * @verbatim @@ -102,13 +104,17 @@ DEAL_II_NAMESPACE_OPEN * @endverbatim * * For each vertex $v_j$ of given cell, there are two edges of which $v_j$ is one of end points. - * Consider a linear function such that 0.5 value at two midpoints of such edges, - * and 0.0 at two midpoints of other edges. + * Consider a linear function such that it has value 0.5 at the midpoints of two adjacent edges, + * and 0.0 at the two midpoints of the other edges. * Note that the set of these values satisfies the dice rule which is described above. * We denote such a function associated with vertex $v_j$ by $\phi_j$. * Then the set of 4 shape functions is a partition of unity on a cell: $\sum_{j=0}^{3} \phi_j = 1$. + * (This is easy to see: at each edge midpoint, the sum of the four function adds up to one + * because two functions have value 0.5 and the other value 0.0. Because the function is globally + * linear, the only function that can have value 1 at four points must also be globally + * equal to one.) * - * The following figures represent $\phi_j$ for $j=0,\cdots,3$ with its midpoint values. + * The following figures represent $\phi_j$ for $j=0,\cdots,3$ with their midpoint values: * * * - * The local DOFs are defined by the coefficients of the shape functions associated vertices, respectively. - * Although these 4 local DOFs are not linearly independent within a single cell as well, - * this definition is a good start point for the definition of the global DOFs. + * The local DoFs are defined by the coefficients of the shape functions associated with vertices, respectively. + * Although these 4 local DoFs are not linearly independent within a single cell, + * this definition is a good start point for the definition of the global DoFs. * * We want to emphasize that the shape functions are constructed on each cell, not on the reference cell only. * Usual finite elements are defined based on a 'parametric' concept. - * It means that a function space for a finite element is defined on one reference cell, and it is transfomed + * It means that a function space for a finite element is defined on one reference cell, and it is transformed * into each cell via a mapping from the reference cell. * However the P1 nonconforming element does not follow such concept. It defines a function space with * linear shape functions on each cell without any help of a function space on the reference cell. * In other words, the element is defined in real space, not via a mapping from a reference cell. + * In this, it is similar to the FE_DGPNonparametric element. * * Thus this implementation does not have to compute shape values on the reference cell. - * Rather than, the shape values are computed by construction of the shape functions + * Rather, the shape values are computed by construction of the shape functions * on each cell independently. * - *

DOFs

- * We have to consider the basis function for the element space in global domain - * because the system of equations which we have to solve at last is for a global system, not local. - * The global basis function associated with a node is defined by a cell-wise composition of + *

Degrees of freedom

+ * We next have to consider the global basis functions for the element + * because the system of equations which we ultimately have to solve is for a global system, not local. + * The global basis functions associated with a node are defined by a cell-wise composition of * local shape functions associated with the node on each element. - * And we define a global DOF associated with a node by a coefficient of the basis function associated with that node. * * There is a theoretical result about the linear independency of the global basis functions * depending on the type of the boundary condition we consider. * - * When the homogeneous Dirichlet boundary condition is given, + * When homogeneous Dirichlet boundary conditions are given, * the global basis functions associated with interior nodes are linearly independent. - * And the number of DOFs is equal to the number of interior nodes, - * same as the number of DOFs for the standard bilinear finite element @p Q_1. + * Then, the number of DoFs is equal to the number of interior nodes, + * and consequently the same as the number of DoFs for the standard bilinear $Q_1$ finite element. * - * When the Neumann boundary condition is given, + * When Neumann boundary conditions are given, * the global basis functions associated with all nodes (including boundary nodes) * are actually not linearly independent. There exists one redundancy. - * Thus in this case, the number of DOFs is equal to the number of all nodes minus 1. + * Thus in this case, the number of DoFs is equal to the number of all nodes minus 1. This is, again + * as for the regular $Q_1$ element. * *

Unit support points

* For a smooth function, we construct a piecewise linear function which belongs to the element space by - * using its nodal values as DOF values. + * using its nodal values as DoF values. * * Note that for the P1 nonconforming element two nodal values of a smooth function and its interpolant do not - * coincide in general, contrast with ordinary Lagrange finite elements. + * coincide in general, in contrast with ordinary Lagrange finite elements. * Of course, it is meaningless to refer 'nodal value' because the element space has nonconformity. - * But it is also true even though the single global basis function associated a node is considered - * with the unique 'nodal value' at the node. + * But it is also true even though the single global basis function associated with a node is considered + * the unique 'nodal value' at the node. * For instance, consider the basis function associated with a node. * Consider two lines representing the level sets for value 0.5 and 0, respectively, by connecting two midpoints. * Then we cut the quad into two sub-triangles by the diagonal which is placed along those two lines. * It gives another level set for value 0.25 which coincides with the cutting diagonal. * Therefore these three level sets are all parallel in the quad and it gives the value 0.75 at the base node, not value 1. - * Even though a general quad is given, this is also true. + * This is true whether the quadrilateral is a rectangle, parallelogram, or any other shape. * *

References

- * You can find the original paper for the P1 nonconforming element which is - * accessible at http://epubs.siam.org/doi/abs/10.1137/S0036142902404923. + * The original paper for the P1 nonconforming element is + * accessible at http://epubs.siam.org/doi/abs/10.1137/S0036142902404923 + * and has the following complete reference: * @code(.bib) * @article{park2003p, * title = {P 1-nonconforming quadrilateral finite element methods for second-order elliptic problems}, @@ -241,6 +249,7 @@ DEAL_II_NAMESPACE_OPEN * } * @endcode * + * @author Jaeryun Yim, 2015, 2016. */ class FE_P1NC : public FiniteElement<2,2> { diff --git a/source/fe/fe_p1nc.cc b/source/fe/fe_p1nc.cc index 873cf159b4..43753aaf14 100644 --- a/source/fe/fe_p1nc.cc +++ b/source/fe/fe_p1nc.cc @@ -91,24 +91,17 @@ std_cxx11::array,4> FE_P1NC::get_linear_shape_coefficients (const Triangulation<2,2>::cell_iterator &cell) { // edge midpoints - Point<2> mpt[4]; - - mpt[0](0) = (cell->vertex(0)(0) + cell->vertex(2)(0))*0.5; - mpt[0](1) = (cell->vertex(0)(1) + cell->vertex(2)(1))*0.5; - - mpt[1](0) = (cell->vertex(1)(0) + cell->vertex(3)(0))*0.5; - mpt[1](1) = (cell->vertex(1)(1) + cell->vertex(3)(1))*0.5; - - mpt[2](0) = (cell->vertex(0)(0) + cell->vertex(1)(0))*0.5; - mpt[2](1) = (cell->vertex(0)(1) + cell->vertex(1)(1))*0.5; - - mpt[3](0) = (cell->vertex(2)(0) + cell->vertex(3)(0))*0.5; - mpt[3](1) = (cell->vertex(2)(1) + cell->vertex(3)(1))*0.5; + const Point<2> mpt[4] = { (cell->vertex(0) + cell->vertex(2)) / 2, + (cell->vertex(1) + cell->vertex(3)) / 2, + (cell->vertex(0) + cell->vertex(1)) / 2, + (cell->vertex(2) + cell->vertex(3)) / 2 + }; // center point - Point<2> cpt; - cpt(0) = (mpt[0](0) + mpt[1](0) + mpt[2](0) + mpt[3](0))*0.25; - cpt(1) = (mpt[0](1) + mpt[1](1) + mpt[2](1) + mpt[3](1))*0.25; + const Point<2> cpt = (cell->vertex(0) + + cell->vertex(1) + + cell->vertex(2) + + cell->vertex(3)) / 4; const double det = (mpt[0](0)-mpt[1](0))*(mpt[2](1)-mpt[3](1)) - (mpt[2](0)-mpt[3](0))*(mpt[0](1)-mpt[1](1)); @@ -143,9 +136,10 @@ FE_P1NC::get_data (const UpdateFlags update_flags, data->update_each = requires_update_flags(update_flags); - // hessian const unsigned int n_q_points = quadrature.size(); output_data.initialize (n_q_points, FE_P1NC(), data->update_each); + + // this is a linear element, so its second derivatives are zero if (data->update_each & update_hessians) output_data.shape_hessians.fill(Tensor<2,2>()); @@ -164,9 +158,10 @@ FE_P1NC::get_face_data (const UpdateFlags update_flags, data->update_each = requires_update_flags(update_flags); - // hessian const unsigned int n_q_points = quadrature.size(); output_data.initialize (n_q_points, FE_P1NC(), data->update_each); + + // this is a linear element, so its second derivatives are zero if (data->update_each & update_hessians) output_data.shape_hessians.fill(Tensor<2,2>()); @@ -185,9 +180,10 @@ FE_P1NC::get_subface_data (const UpdateFlags update_flags, data->update_each = requires_update_flags(update_flags); - // hessian const unsigned int n_q_points = quadrature.size(); output_data.initialize (n_q_points, FE_P1NC(), data->update_each); + + // this is a linear element, so its second derivatives are zero if (data->update_each & update_hessians) output_data.shape_hessians.fill(Tensor<2,2>()); @@ -198,10 +194,10 @@ FE_P1NC::get_subface_data (const UpdateFlags update_flags, void FE_P1NC::fill_fe_values (const Triangulation<2,2>::cell_iterator &cell, - const CellSimilarity::Similarity cell_similarity, - const Quadrature<2> &quadrature, - const Mapping<2,2> &mapping, - const Mapping<2,2>::InternalDataBase &mapping_internal, + const CellSimilarity::Similarity , + const Quadrature<2> &, + const Mapping<2,2> &, + const Mapping<2,2>::InternalDataBase &, const internal::FEValues::MappingRelatedData<2,2> &mapping_data, const FiniteElement<2,2>::InternalDataBase &fe_internal, internal::FEValues::FiniteElementRelatedData<2,2> &output_data) const @@ -210,32 +206,22 @@ FE_P1NC::fill_fe_values (const Triangulation<2,2>::cell_iterator &cell const unsigned int n_q_points = mapping_data.quadrature_points.size(); - std::vector values(flags & update_values ? this->dofs_per_cell : 0); - std::vector > grads(flags & update_gradients ? this->dofs_per_cell : 0); - // linear shape functions std_cxx11::array,4> coeffs = get_linear_shape_coefficients (cell); // compute on the cell - if (flags & (update_values | update_gradients)) + if (flags & update_values) for (unsigned int i=0; idofs_per_cell; ++k) - { - if (flags & update_values) - { - values[k] = coeffs[k][0]*mapping_data.quadrature_points[i](0) + coeffs[k][1]*mapping_data.quadrature_points[i](1) + coeffs[k][2]; - output_data.shape_values[k][i] = values[k]; - } - - if (flags & update_gradients) - { - grads[k][0] = coeffs[k][0]; - grads[k][1] = coeffs[k][1]; - output_data.shape_gradients[k][i] = grads[k]; - } - } - } + for (unsigned int k=0; kdofs_per_cell; ++k) + output_data.shape_values[k][i] = (coeffs[k][0]*mapping_data.quadrature_points[i](0) + + coeffs[k][1]*mapping_data.quadrature_points[i](1) + + coeffs[k][2]); + + if (flags & update_gradients) + for (unsigned int i=0; idofs_per_cell; ++k) + output_data.shape_gradients[k][i] = Point<2>(coeffs[k][0], + coeffs[k][1]); } @@ -245,43 +231,37 @@ FE_P1NC::fill_fe_face_values (const Triangulation<2,2>::cell_iterator const unsigned int face_no, const Quadrature<1> &quadrature, const Mapping<2,2> &mapping, - const Mapping<2,2>::InternalDataBase &mapping_internal, - const dealii::internal::FEValues::MappingRelatedData<2,2> &mapping_data, + const Mapping<2,2>::InternalDataBase &, + const dealii::internal::FEValues::MappingRelatedData<2,2> &, const InternalDataBase &fe_internal, dealii::internal::FEValues::FiniteElementRelatedData<2,2> &output_data) const { const UpdateFlags flags(fe_internal.update_each); - const unsigned int n_q_points = mapping_data.quadrature_points.size(); - - std::vector values(flags & update_values ? this->dofs_per_cell : 0); - std::vector > grads(flags & update_gradients ? this->dofs_per_cell : 0); - // linear shape functions - std_cxx11::array,4> coeffs = get_linear_shape_coefficients (cell); + const std_cxx11::array,4> coeffs + = get_linear_shape_coefficients (cell); // compute on the face - Quadrature<2> quadrature_on_face = QProjector<2>::project_to_face(quadrature, face_no); - Point<2> quadrature_point; - for (unsigned int i=0; i quadrature_on_face = QProjector<2>::project_to_face(quadrature, face_no); + + if (flags & update_values) + for (unsigned int i=0; idofs_per_cell; ++k) { - if (flags & update_values) - { - quadrature_point = mapping.transform_unit_to_real_cell(cell, quadrature_on_face.point(i)); - values[k] = coeffs[k][0]*quadrature_point(0) + coeffs[k][1]*quadrature_point(1) + coeffs[k][2]; - output_data.shape_values[k][i] = values[k]; - } - - if (flags & update_gradients) - { - grads[k][0] = coeffs[k][0]; - grads[k][1] = coeffs[k][1]; - output_data.shape_gradients[k][i] = grads[k]; - } + const Point<2> quadrature_point + = mapping.transform_unit_to_real_cell(cell, quadrature_on_face.point(i)); + + output_data.shape_values[k][i] = (coeffs[k][0]*quadrature_point(0) + + coeffs[k][1]*quadrature_point(1) + + coeffs[k][2]); } - } + + if (flags & update_gradients) + for (unsigned int i=0; idofs_per_cell; ++k) + output_data.shape_gradients[k][i] = Point<2>(coeffs[k][0], + coeffs[k][1]); } @@ -292,54 +272,45 @@ FE_P1NC::fill_fe_subface_values (const Triangulation<2,2>::cell_iterator const unsigned int sub_no, const Quadrature<1> &quadrature, const Mapping<2,2> &mapping, - const Mapping<2,2>::InternalDataBase &mapping_internal, - const dealii::internal::FEValues::MappingRelatedData<2,2> &mapping_data, + const Mapping<2,2>::InternalDataBase &, + const dealii::internal::FEValues::MappingRelatedData<2,2> &, const InternalDataBase &fe_internal, dealii::internal::FEValues::FiniteElementRelatedData<2,2> &output_data) const { const UpdateFlags flags(fe_internal.update_each); - const unsigned int n_q_points = mapping_data.quadrature_points.size(); - - std::vector values(flags & update_values ? this->dofs_per_cell : 0); - std::vector > grads(flags & update_gradients ? this->dofs_per_cell : 0); - // linear shape functions - std_cxx11::array,4> coeffs = get_linear_shape_coefficients (cell); + const std_cxx11::array,4> coeffs + = get_linear_shape_coefficients (cell); // compute on the subface - Quadrature<2> quadrature_on_subface = QProjector<2>::project_to_subface(quadrature, face_no, sub_no); - Point<2> quadrature_point; - for (unsigned int i=0; i quadrature_on_subface = QProjector<2>::project_to_subface(quadrature, face_no, sub_no); + + if (flags & update_values) + for (unsigned int i=0; idofs_per_cell; ++k) + { + const Point<2> quadrature_point + = mapping.transform_unit_to_real_cell(cell, quadrature_on_subface.point(i)); + + output_data.shape_values[k][i] = (coeffs[k][0]*quadrature_point(0) + + coeffs[k][1]*quadrature_point(1) + + coeffs[k][2]); + } + } + + if (flags & update_gradients) + for (unsigned int i=0; idofs_per_cell; ++k) - { - if (flags & update_values) - { - quadrature_point = mapping.transform_unit_to_real_cell(cell, quadrature_on_subface.point(i)); - values[k] = coeffs[k][0]*quadrature_point(0) + coeffs[k][1]*quadrature_point(1) + coeffs[k][2]; - output_data.shape_values[k][i] = values[k]; - } - - if (flags & update_gradients) - { - grads[k][0] = coeffs[k][0]; - grads[k][1] = coeffs[k][1]; - output_data.shape_gradients[k][i] = grads[k]; - } - } - } + output_data.shape_gradients[k][i] = Point<2>(coeffs[k][0], + coeffs[k][1]); } void FE_P1NC::initialize_constraints () { - std::vector > constraint_points; - - // Add midpoint - constraint_points.push_back (Point<1> (0.5)); - // coefficient relation between children and mother interface_constraints .TableBase<2,double>::reinit (interface_constraints_size());