/*@{*/
/**
- * 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 <i>some level</i> of continuity.
*
* <h3>Dice Rule</h3>
* Since any function in the P1 nonconforming space is piecewise linear on each element,
* 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.
*
* <h3>Shape functions</h3>
- * 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
* @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:
*
* <ul>
* <li> shape function $\phi_0$:
*
* </ul>
*
- * 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.
*
- * <h3>DOFs</h3>
- * 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
+ * <h3>Degrees of freedom</h3>
+ * We next have to consider the <i>global</i> 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.
*
* <h3>Unit support points</h3>
* 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.
*
* <h3>References</h3>
- * 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},
* }
* @endcode
*
+ * @author Jaeryun Yim, 2015, 2016.
*/
class FE_P1NC : public FiniteElement<2,2>
{
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));
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>());
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>());
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>());
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
const unsigned int n_q_points = mapping_data.quadrature_points.size();
- std::vector<double> values(flags & update_values ? this->dofs_per_cell : 0);
- std::vector<Tensor<1,2> > grads(flags & update_gradients ? this->dofs_per_cell : 0);
-
// linear shape functions
std_cxx11::array<std_cxx11::array<double,3>,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; i<n_q_points; ++i)
- {
- for (unsigned int k=0; k<this->dofs_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; k<this->dofs_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; i<n_q_points; ++i)
+ for (unsigned int k=0; k<this->dofs_per_cell; ++k)
+ output_data.shape_gradients[k][i] = Point<2>(coeffs[k][0],
+ coeffs[k][1]);
}
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<double> values(flags & update_values ? this->dofs_per_cell : 0);
- std::vector<Tensor<1,2> > grads(flags & update_gradients ? this->dofs_per_cell : 0);
-
// linear shape functions
- std_cxx11::array<std_cxx11::array<double,3>,4> coeffs = get_linear_shape_coefficients (cell);
+ const std_cxx11::array<std_cxx11::array<double,3>,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.size(); ++i)
- {
+ const Quadrature<2> quadrature_on_face = QProjector<2>::project_to_face(quadrature, face_no);
+
+ if (flags & update_values)
+ for (unsigned int i=0; i<quadrature_on_face.size(); ++i)
for (unsigned int k=0; k<this->dofs_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; i<quadrature_on_face.size(); ++i)
+ for (unsigned int k=0; k<this->dofs_per_cell; ++k)
+ output_data.shape_gradients[k][i] = Point<2>(coeffs[k][0],
+ coeffs[k][1]);
}
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<double> values(flags & update_values ? this->dofs_per_cell : 0);
- std::vector<Tensor<1,2> > grads(flags & update_gradients ? this->dofs_per_cell : 0);
-
// linear shape functions
- std_cxx11::array<std_cxx11::array<double,3>,4> coeffs = get_linear_shape_coefficients (cell);
+ const std_cxx11::array<std_cxx11::array<double,3>,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.size(); ++i)
- {
+ const Quadrature<2> quadrature_on_subface = QProjector<2>::project_to_subface(quadrature, face_no, sub_no);
+
+ if (flags & update_values)
+ for (unsigned int i=0; i<quadrature_on_subface.size(); ++i)
+ {
+ for (unsigned int k=0; k<this->dofs_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; i<quadrature_on_subface.size(); ++i)
for (unsigned int k=0; k<this->dofs_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<Point<1> > 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());