fe_xyz_boundary.C

Go to the documentation of this file.

// The libMesh Finite Element Library.
// Copyright (C) 2002-2026 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
 
// This library is free software; you can 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.
 
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.
 
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 
 
 
// C++ includes
#include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
#include <cmath> // for std::sqrt
 
 
// Local includes
#include "libmesh/libmesh_common.h"
#include "libmesh/fe.h"
#include "libmesh/quadrature.h"
#include "libmesh/elem.h"
#include "libmesh/libmesh_logging.h"
 
namespace libMesh
{
 
 
 
 
//-------------------------------------------------------
// Full specialization for 0D when this is a useless method
template <>
void FEXYZ<0>::reinit(const Elem *,
                      const unsigned int,
                      const Real,
                      const std::vector<Point> * const,
                      const std::vector<Real> * const)
{
  libmesh_error_msg("ERROR: This method only makes sense for 2D/3D elements!");
}
 
 
//-------------------------------------------------------
// Full specialization for 1D when this is a useless method
template <>
void FEXYZ<1>::reinit(const Elem *,
                      const unsigned int,
                      const Real,
                      const std::vector<Point> * const,
                      const std::vector<Real> * const)
{
  libmesh_error_msg("ERROR: This method only makes sense for 2D/3D elements!");
}
 
 
//-------------------------------------------------------
// Method for 2D, 3D
template <unsigned int Dim>
void FEXYZ<Dim>::reinit(const Elem * elem,
                        const unsigned int s,
                        const Real,
                        const std::vector<Point> * const pts,
                        const std::vector<Real> * const weights)
{
  libmesh_assert(elem);
  libmesh_assert (this->qrule != nullptr || pts != nullptr);
  // We don't do this for 1D elements!
  libmesh_assert_not_equal_to (Dim, 1);
 
  // Build the side of interest
  const std::unique_ptr<const Elem> side(elem->build_side_ptr(s));
 
  // Find the max p_level to select
  // the right quadrature rule for side integration
  unsigned int side_p_level = elem->p_level();
  if (elem->neighbor_ptr(s) != nullptr)
    side_p_level = std::max(side_p_level, elem->neighbor_ptr(s)->p_level());
 
  // Initialize the shape functions at the user-specified
  // points
  if (pts != nullptr)
    {
      // We can't get away without recomputing shape functions next
      // time
      this->shapes_on_quadrature = false;
 
      // Set the element type
      this->_elem = elem;
      this->_elem_type = elem->type();
 
      // Initialize the face shape functions
      this->_fe_map->template init_face_shape_functions<Dim>(*pts,  side.get());
      if (weights != nullptr)
        {
          this->compute_face_values (elem, side.get(), *weights);
        }
      else
        {
          std::vector<Real> dummy_weights (pts->size(), 1.);
          // Compute data on the face for integration
          this->compute_face_values (elem, side.get(), dummy_weights);
        }
    }
  else
    {
      // initialize quadrature rule
      this->qrule->init(*side, side_p_level);
 
      {
        // Set the element
        this->_elem = elem;
        this->_elem_type = elem->type();
 
        // Initialize the face shape functions
        this->_fe_map->template init_face_shape_functions<Dim>(this->qrule->get_points(),  side.get());
      }
      // We can't get away without recomputing shape functions next
      // time
      this->shapes_on_quadrature = false;
      // Compute data on the face for integration
      this->compute_face_values (elem, side.get(), this->qrule->get_weights());
    }
}
 
 
 
template <unsigned int Dim>
void FEXYZ<Dim>::compute_face_values(const Elem * elem,
                                     const Elem * side,
                                     const std::vector<Real> & qw)
{
  libmesh_assert(elem);
  libmesh_assert(side);
 
  LOG_SCOPE("compute_face_values()", "FEXYZ");
 
  // The number of quadrature points.
  const std::size_t n_qp = qw.size();
 
  // Number of shape functions in the finite element approximation
  // space.
  const unsigned int n_approx_shape_functions =
    this->n_dofs(elem, this->get_order());
 
  // Resize the shape functions and their gradients
  this->phi.resize    (n_approx_shape_functions);
  this->dphi.resize   (n_approx_shape_functions);
  this->dphidx.resize (n_approx_shape_functions);
  this->dphidy.resize (n_approx_shape_functions);
  this->dphidz.resize (n_approx_shape_functions);
 
  for (unsigned int i=0; i<n_approx_shape_functions; i++)
    {
      this->phi[i].resize    (n_qp);
      this->dphi[i].resize   (n_qp);
      this->dphidx[i].resize (n_qp);
      this->dphidy[i].resize (n_qp);
      this->dphidz[i].resize (n_qp);
    }
 
  this->_fe_map->compute_face_map(this->dim, qw, side);
 
  const std::vector<libMesh::Point> & xyz = this->_fe_map->get_xyz();
 
  switch (this->dim)
    {
      // A 2D finite element living in either 2D or 3D space.
      // This means the boundary is a 1D finite element, i.e.
      // and EDGE2 or EDGE3.
    case 2:
      {
        // compute the shape function values & gradients
        for (unsigned int i=0; i<n_approx_shape_functions; i++)
          for (std::size_t p=0; p<n_qp; p++)
            {
              this->phi[i][p] = FE<Dim,XYZ>::shape (elem, this->fe_type.order, i, xyz[p]);
 
              this->dphi[i][p](0) =
                this->dphidx[i][p] = FE<Dim,XYZ>::shape_deriv (elem, this->fe_type.order, i, 0, xyz[p]);
 
              this->dphi[i][p](1) =
                this->dphidy[i][p] = FE<Dim,XYZ>::shape_deriv (elem, this->fe_type.order, i, 1, xyz[p]);
 
#if LIBMESH_DIM == 3
              this->dphi[i][p](2) = // can only assign to the Z component if LIBMESH_DIM==3
#endif
                this->dphidz[i][p] = 0.;
            }
 
        // done computing face values
        break;
      }
 
      // A 3D finite element living in 3D space.
    case 3:
      {
        // compute the shape function values & gradients
        for (unsigned int i=0; i<n_approx_shape_functions; i++)
          for (std::size_t p=0; p<n_qp; p++)
            {
              this->phi[i][p] = FE<Dim,XYZ>::shape (elem, this->fe_type.order, i, xyz[p]);
 
              this->dphi[i][p](0) =
                this->dphidx[i][p] = FE<Dim,XYZ>::shape_deriv (elem, this->fe_type.order, i, 0, xyz[p]);
 
              this->dphi[i][p](1) =
                this->dphidy[i][p] = FE<Dim,XYZ>::shape_deriv (elem, this->fe_type.order, i, 1, xyz[p]);
 
              this->dphi[i][p](2) =
                this->dphidz[i][p] = FE<Dim,XYZ>::shape_deriv (elem, this->fe_type.order, i, 2, xyz[p]);
            }
 
        // done computing face values
        break;
      }
 
    default:
      libmesh_error_msg("Invalid dim " << this->dim);
    }
}
 
 
 
 
//--------------------------------------------------------------
// Explicit instantiations (doesn't make sense in 1D!) using fe_macro.h's macro
 
template class LIBMESH_EXPORT FEXYZ<2>;
template class LIBMESH_EXPORT FEXYZ<3>;
 
 
} // namespace libMesh