quadrature_composite.C

Go to the documentation of this file.

// The libMesh Finite Element Library.
// Copyright (C) 2002-2019 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
 
 
#include "libmesh/libmesh_config.h"
#if defined(LIBMESH_HAVE_TRIANGLE) && defined(LIBMESH_HAVE_TETGEN)
 
#include "libmesh/fe.h"
#include "libmesh/quadrature_gauss.h"
#include "libmesh/quadrature_trap.h"
#include "libmesh/quadrature_simpson.h"
#include "libmesh/quadrature_composite.h"
#include "libmesh/elem.h"
#include "libmesh/enum_quadrature_type.h"
 
 
 
namespace libMesh
{
 
 
template <class QSubCell>
QuadratureType QComposite<QSubCell>::type() const
{
  return QCOMPOSITE;
}
 
 
 
template <class QSubCell>
QComposite<QSubCell>::QComposite(unsigned int d,
                                 Order o) :
  QSubCell(d,o), // explicitly call base class constructor
  _q_subcell(d,o),
  _lagrange_fe(FEBase::build (d, FEType (FIRST, LAGRANGE)))
{
  // explicitly call the init function in 1D since the
  // other tensor-product rules require this one.
  // note that EDGE will not be used internally, however
  // if we called the function with INVALID_ELEM it would try to
  // be smart and return, thinking it had already done the work.
  if (_dim == 1)
    QSubCell::init(EDGE2);
 
  libmesh_assert (_lagrange_fe.get() != nullptr);
 
  _lagrange_fe->attach_quadrature_rule (&_q_subcell);
}
 
 
 
template <class QSubCell>
void QComposite<QSubCell>::init (const Elem & elem,
                                 const std::vector<Real> & vertex_distance_func,
                                 unsigned int p_level)
{
  libmesh_assert_equal_to (vertex_distance_func.size(), elem.n_vertices());
  libmesh_assert_equal_to (_dim, elem.dim());
 
  // We already initialized a Lagrange map; we're not supporting
  // others yet
  libmesh_assert_equal_to (elem.mapping_type(), LAGRANGE_MAP);
 
  // if we are not cut, revert to simple base class init() method.
  if (!_elem_cutter.is_cut (elem, vertex_distance_func))
    {
      _q_subcell.init (elem.type(), p_level);
      _points  = _q_subcell.get_points();
      _weights = _q_subcell.get_weights();
 
      //this->print_info();
      return;
    }
 
  // Get a pointer to the element's reference element.  We want to
  // perform cutting on the reference element such that the quadrature
  // point locations of the subelements live in the reference
  // coordinate system, thereby eliminating the need for inverse
  // mapping.
  const Elem * reference_elem = elem.reference_elem();
 
  libmesh_assert (reference_elem != nullptr);
 
  _elem_cutter(*reference_elem, vertex_distance_func);
  //_elem_cutter(elem, vertex_distance_func);
 
  // clear our state & accumulate points from subelements
  _points.clear();
  _weights.clear();
 
  // inside subelem
  {
    const std::vector<Elem const *> & inside_elem (_elem_cutter.inside_elements());
    std::cout << inside_elem.size() << " elements inside\n";
 
    this->add_subelem_values(inside_elem);
  }
 
  // outside subelem
  {
    const std::vector<Elem const *> & outside_elem (_elem_cutter.outside_elements());
    std::cout << outside_elem.size() << " elements outside\n";
 
    this->add_subelem_values(outside_elem);
  }
 
  this->print_info();
}
 
 
 
template <class QSubCell>
void QComposite<QSubCell>::add_subelem_values (const std::vector<Elem const *> & subelem)
 
{
  const std::vector<Real>  & subelem_weights = _lagrange_fe->get_JxW();
  const std::vector<Point> & subelem_points  = _lagrange_fe->get_xyz();
 
  for (const auto & elem : subelem)
    {
      // tetgen seems to create 0-volume cells on occasion, but we *should*
      // be catching that appropriately now inside the ElemCutter class.
      // Just in case trap here, describe the error, and abort.
      libmesh_try
        {
          _lagrange_fe->reinit(elem);
          _weights.insert(_weights.end(),
                          subelem_weights.begin(), subelem_weights.end());
 
          _points.insert(_points.end(),
                         subelem_points.begin(), subelem_points.end());
        }
      libmesh_catch (...)
        {
          libMesh::err << "ERROR: found a bad cut cell!\n";
 
          for (auto n : elem->node_index_range())
            libMesh::err << elem->point(n) << std::endl;
 
          libmesh_error_msg("Tetgen may have created a 0-volume cell during Cutcell integration.");
        }
    }
}
 
 
 
//--------------------------------------------------------------
// Explicit instantiations
template class QComposite<QGauss>;
template class QComposite<QTrap>;
template class QComposite<QSimpson>;
 
} // namespace libMesh
 
#endif // LIBMESH_HAVE_TRIANGLE && LIBMESH_HAVE_TETGEN