fe_p1nc.cc

Go to the documentation of this file.

// -----------------------------------------------------------------------------
//
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
// Copyright (C) 2016 - 2026 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Detailed license information governing the source code and contributions
// can be found in LICENSE.md and CONTRIBUTING.md at the top level directory.
//
// -----------------------------------------------------------------------------
 
 
#include <deal.II/base/config.h>
 
#include <deal.II/base/exception_macros.h>
#include <deal.II/base/geometry_info.h>
#include <deal.II/base/point.h>
#include <deal.II/base/qprojector.h>
#include <deal.II/base/quadrature.h>
#include <deal.II/base/table.h>
#include <deal.II/base/tensor.h>
#include <deal.II/base/types.h>
 
#include <deal.II/fe/component_mask.h>
#include <deal.II/fe/fe.h>
#include <deal.II/fe/fe_data.h>
#include <deal.II/fe/fe_p1nc.h>
#include <deal.II/fe/mapping_related_data.h>
 
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
 
#include <deal.II/hp/collection.h>
#include <deal.II/hp/q_collection.h>
 
#include <deal.II/lac/full_matrix.h>
 
#include <Kokkos_Macros.hpp>
 
#include <array>
#include <memory>
#include <ostream>
 
DEAL_II_NAMESPACE_OPEN
 
 
FE_P1NC::FE_P1NC()
  : FiniteElement<2, 2>(
      FiniteElementData<2>(get_dpo_vector(), 1, 1, FiniteElementData<2>::L2),
      std::vector<bool>(1, false),
      std::vector<ComponentMask>(4, ComponentMask(1, true)))
{
  // face support points: 2 end vertices
  unit_face_support_points[0].resize(2);
  unit_face_support_points[0][0][0] = 0.0;
  unit_face_support_points[0][1][0] = 1.0;
 
  // initialize constraints matrix
  initialize_constraints();
}
 
 
 
std::string
FE_P1NC::get_name() const
{
  return "FE_P1NC";
}
 
 
 
UpdateFlags
FE_P1NC::requires_update_flags(const UpdateFlags flags) const
{
  UpdateFlags out = update_default;
 
  if (flags & update_values)
    out |= update_values | update_quadrature_points;
  if (flags & update_gradients)
    out |= update_gradients;
  if (flags & update_normal_vectors)
    out |= update_normal_vectors | update_JxW_values;
  if (flags & update_hessians)
    out |= update_hessians;
 
  return out;
}
 
 
 
std::unique_ptr<FiniteElement<2, 2>>
FE_P1NC::clone() const
{
  return std::make_unique<FE_P1NC>(*this);
}
 
 
 
std::vector<unsigned int>
FE_P1NC::get_dpo_vector()
{
  std::vector<unsigned int> dpo(3);
  dpo[0] = 1; // dofs per object: vertex
  dpo[1] = 0; // line
  dpo[2] = 0; // quad
  return dpo;
}
 
 
 
ndarray<double, 4, 3>
FE_P1NC::get_linear_shape_coefficients(
  const Triangulation<2, 2>::cell_iterator &cell)
{
  // edge midpoints
  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
  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]);
 
  ndarray<double, 4, 3> coeffs;
  coeffs[0][0] =
    ((mpt[2][1] - mpt[3][1]) * (0.5) - (mpt[0][1] - mpt[1][1]) * (0.5)) / det;
  coeffs[1][0] =
    ((mpt[2][1] - mpt[3][1]) * (-0.5) - (mpt[0][1] - mpt[1][1]) * (0.5)) / det;
  coeffs[2][0] =
    ((mpt[2][1] - mpt[3][1]) * (0.5) - (mpt[0][1] - mpt[1][1]) * (-0.5)) / det;
  coeffs[3][0] =
    ((mpt[2][1] - mpt[3][1]) * (-0.5) - (mpt[0][1] - mpt[1][1]) * (-0.5)) / det;
 
  coeffs[0][1] =
    (-(mpt[2][0] - mpt[3][0]) * (0.5) + (mpt[0][0] - mpt[1][0]) * (0.5)) / det;
  coeffs[1][1] =
    (-(mpt[2][0] - mpt[3][0]) * (-0.5) + (mpt[0][0] - mpt[1][0]) * (0.5)) / det;
  coeffs[2][1] =
    (-(mpt[2][0] - mpt[3][0]) * (0.5) + (mpt[0][0] - mpt[1][0]) * (-0.5)) / det;
  coeffs[3][1] =
    (-(mpt[2][0] - mpt[3][0]) * (-0.5) + (mpt[0][0] - mpt[1][0]) * (-0.5)) /
    det;
 
  coeffs[0][2] = 0.25 - cpt[0] * coeffs[0][0] - cpt[1] * coeffs[0][1];
  coeffs[1][2] = 0.25 - cpt[0] * coeffs[1][0] - cpt[1] * coeffs[1][1];
  coeffs[2][2] = 0.25 - cpt[0] * coeffs[2][0] - cpt[1] * coeffs[2][1];
  coeffs[3][2] = 0.25 - cpt[0] * coeffs[3][0] - cpt[1] * coeffs[3][1];
 
  return coeffs;
}
 
 
 
std::unique_ptr<FiniteElement<2, 2>::InternalDataBase>
FE_P1NC::get_data(
  const UpdateFlags update_flags,
  const Mapping<2, 2> &,
  const Quadrature<2> &quadrature,
  ::internal::FEValuesImplementation::FiniteElementRelatedData<2, 2>
    &output_data) const
{
  auto data_ptr = std::make_unique<FiniteElement<2, 2>::InternalDataBase>();
 
  data_ptr->update_each = requires_update_flags(update_flags);
 
  const unsigned int n_q_points = quadrature.size();
  output_data.initialize(n_q_points, FE_P1NC(), data_ptr->update_each);
 
  // this is a linear element, so its second derivatives are zero
  if (data_ptr->update_each & update_hessians)
    output_data.shape_hessians.fill(Tensor<2, 2>());
 
  return data_ptr;
}
 
 
 
std::unique_ptr<FiniteElement<2, 2>::InternalDataBase>
FE_P1NC::get_face_data(
  const UpdateFlags update_flags,
  const Mapping<2, 2> &,
  const hp::QCollection<1> &quadrature,
  ::internal::FEValuesImplementation::FiniteElementRelatedData<2, 2>
    &output_data) const
{
  AssertDimension(quadrature.size(), 1);
 
  auto data_ptr = std::make_unique<FiniteElement<2, 2>::InternalDataBase>();
 
  data_ptr->update_each = requires_update_flags(update_flags);
 
  const unsigned int n_q_points = quadrature[0].size();
  output_data.initialize(n_q_points, FE_P1NC(), data_ptr->update_each);
 
  // this is a linear element, so its second derivatives are zero
  if (data_ptr->update_each & update_hessians)
    output_data.shape_hessians.fill(Tensor<2, 2>());
 
  return data_ptr;
}
 
 
 
std::unique_ptr<FiniteElement<2, 2>::InternalDataBase>
FE_P1NC::get_subface_data(
  const UpdateFlags update_flags,
  const Mapping<2, 2> &,
  const Quadrature<1> &quadrature,
  ::internal::FEValuesImplementation::FiniteElementRelatedData<2, 2>
    &output_data) const
{
  auto data_ptr = std::make_unique<FiniteElement<2, 2>::InternalDataBase>();
 
  data_ptr->update_each = requires_update_flags(update_flags);
 
  const unsigned int n_q_points = quadrature.size();
  output_data.initialize(n_q_points, FE_P1NC(), data_ptr->update_each);
 
  // this is a linear element, so its second derivatives are zero
  if (data_ptr->update_each & update_hessians)
    output_data.shape_hessians.fill(Tensor<2, 2>());
 
  return data_ptr;
}
 
 
 
void
FE_P1NC::fill_fe_values(
  const Triangulation<2, 2>::cell_iterator &cell,
  const CellSimilarity::Similarity,
  const Quadrature<2> &,
  const Mapping<2, 2> &,
  const Mapping<2, 2>::InternalDataBase &,
  const internal::FEValuesImplementation::MappingRelatedData<2, 2>
                                              &mapping_data,
  const FiniteElement<2, 2>::InternalDataBase &fe_internal,
  internal::FEValuesImplementation::FiniteElementRelatedData<2, 2> &output_data)
  const
{
  const UpdateFlags flags(fe_internal.update_each);
 
  const unsigned int n_q_points = mapping_data.quadrature_points.size();
 
  // linear shape functions
  ndarray<double, 4, 3> coeffs = get_linear_shape_coefficients(cell);
 
  // compute on the cell
  if (flags & update_values)
    for (unsigned int i = 0; i < n_q_points; ++i)
      for (unsigned int k = 0; k < this->n_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->n_dofs_per_cell(); ++k)
        output_data.shape_gradients[k][i] =
          Point<2>(coeffs[k][0], coeffs[k][1]);
}
 
 
 
void
FE_P1NC::fill_fe_face_values(
  const Triangulation<2, 2>::cell_iterator &cell,
  const unsigned int                        face_no,
  const hp::QCollection<1>                 &quadrature,
  const Mapping<2, 2>                      &mapping,
  const Mapping<2, 2>::InternalDataBase &,
  const internal::FEValuesImplementation::MappingRelatedData<2, 2> &,
  const InternalDataBase &fe_internal,
  ::internal::FEValuesImplementation::FiniteElementRelatedData<2, 2>
    &output_data) const
{
  AssertDimension(quadrature.size(), 1);
 
  const UpdateFlags flags(fe_internal.update_each);
 
  // linear shape functions
  const ndarray<double, 4, 3> coeffs = get_linear_shape_coefficients(cell);
 
  // compute on the face
  const Quadrature<2> quadrature_on_face =
    QProjector<2>::project_to_face(this->reference_cell(),
                                   quadrature[0],
                                   face_no,
                                   numbers::default_geometric_orientation);
 
  if (flags & update_values)
    for (unsigned int i = 0; i < quadrature_on_face.size(); ++i)
      for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++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->n_dofs_per_cell(); ++k)
        output_data.shape_gradients[k][i] =
          Point<2>(coeffs[k][0], coeffs[k][1]);
}
 
 
 
void
FE_P1NC::fill_fe_subface_values(
  const Triangulation<2, 2>::cell_iterator &cell,
  const unsigned int                        face_no,
  const unsigned int                        sub_no,
  const Quadrature<1>                      &quadrature,
  const Mapping<2, 2>                      &mapping,
  const Mapping<2, 2>::InternalDataBase &,
  const internal::FEValuesImplementation::MappingRelatedData<2, 2> &,
  const InternalDataBase &fe_internal,
  ::internal::FEValuesImplementation::FiniteElementRelatedData<2, 2>
    &output_data) const
{
  const UpdateFlags flags(fe_internal.update_each);
 
  // linear shape functions
  const ndarray<double, 4, 3> coeffs = get_linear_shape_coefficients(cell);
 
  // compute on the subface
  const Quadrature<2> quadrature_on_subface =
    QProjector<2>::project_to_subface(this->reference_cell(),
                                      quadrature,
                                      face_no,
                                      sub_no,
                                      numbers::default_geometric_orientation,
                                      RefinementCase<1>::isotropic_refinement);
 
  if (flags & update_values)
    for (unsigned int i = 0; i < quadrature_on_subface.size(); ++i)
      {
        for (unsigned int k = 0; k < this->n_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->n_dofs_per_cell(); ++k)
        output_data.shape_gradients[k][i] =
          Point<2>(coeffs[k][0], coeffs[k][1]);
}
 
 
 
void
FE_P1NC::initialize_constraints()
{
  // coefficient relation between children and mother
  interface_constraints.TableBase<2, double>::reinit(
    interface_constraints_size());
 
  interface_constraints(0, 0) = 0.5;
  interface_constraints(0, 1) = 0.5;
}
 
 
DEAL_II_NAMESPACE_CLOSE