Line data    Source code

       1             : // The libMesh Finite Element Library.
       2             : // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
       3             : 
       4             : // This library is free software; you can redistribute it and/or
       5             : // modify it under the terms of the GNU Lesser General Public
       6             : // License as published by the Free Software Foundation; either
       7             : // version 2.1 of the License, or (at your option) any later version.
       8             : 
       9             : // This library is distributed in the hope that it will be useful,
      10             : // but WITHOUT ANY WARRANTY; without even the implied warranty of
      11             : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
      12             : // Lesser General Public License for more details.
      13             : 
      14             : // You should have received a copy of the GNU Lesser General Public
      15             : // License along with this library; if not, write to the Free Software
      16             : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
      17             : 
      18             : 
      19             : #include "libmesh/libmesh_config.h"
      20             : #if defined(LIBMESH_HAVE_TRIANGLE) && defined(LIBMESH_HAVE_TETGEN)
      21             : 
      22             : #include "libmesh/fe.h"
      23             : #include "libmesh/quadrature_gauss.h"
      24             : #include "libmesh/quadrature_trap.h"
      25             : #include "libmesh/quadrature_simpson.h"
      26             : #include "libmesh/quadrature_composite.h"
      27             : #include "libmesh/elem.h"
      28             : #include "libmesh/enum_quadrature_type.h"
      29             : 
      30             : 
      31             : 
      32             : namespace libMesh
      33             : {
      34             : 
      35             : 
      36             : template <class QSubCell>
      37           0 : QuadratureType QComposite<QSubCell>::type() const
      38             : {
      39           0 :   return QCOMPOSITE;
      40             : }
      41             : 
      42             : 
      43             : 
      44             : template <class QSubCell>
      45           8 : QComposite<QSubCell>::QComposite(unsigned int d,
      46             :                                  Order o) :
      47             :   QSubCell(d,o), // explicitly call base class constructor
      48             :   _q_subcell(d,o),
      49           8 :   _lagrange_fe(FEBase::build (d, FEType (FIRST, LAGRANGE)))
      50             : {
      51             :   // explicitly call the init function in 1D since the
      52             :   // other tensor-product rules require this one.
      53             :   // note that EDGE will not be used internally, however
      54             :   // if we called the function with INVALID_ELEM it would try to
      55             :   // be smart and return, thinking it had already done the work.
      56           8 :   if (_dim == 1)
      57           0 :     QSubCell::init(EDGE2);
      58             : 
      59           4 :   libmesh_assert (_lagrange_fe.get() != nullptr);
      60             : 
      61           8 :   _lagrange_fe->attach_quadrature_rule (&_q_subcell);
      62           8 : }
      63             : 
      64             : 
      65             : 
      66             : template <class QSubCell>
      67       37888 : void QComposite<QSubCell>::init (const Elem & elem,
      68             :                                  const std::vector<Real> & vertex_distance_func,
      69             :                                  unsigned int p_level)
      70             : {
      71       18944 :   libmesh_assert_equal_to (vertex_distance_func.size(), elem.n_vertices());
      72       18944 :   libmesh_assert_equal_to (_dim, elem.dim());
      73             : 
      74             :   // We already initialized a Lagrange map; we're not supporting
      75             :   // others yet
      76       18944 :   libmesh_assert_equal_to (elem.mapping_type(), LAGRANGE_MAP);
      77             : 
      78             :   // if we are not cut, revert to simple base class init() method.
      79       37888 :   if (!_elem_cutter.is_cut (elem, vertex_distance_func))
      80             :     {
      81       37118 :       _q_subcell.init (elem, p_level);
      82       37118 :       _points  = _q_subcell.get_points();
      83       37118 :       _weights = _q_subcell.get_weights();
      84             : 
      85             :       //this->print_info();
      86       37118 :       return;
      87             :     }
      88             : 
      89             :   // Get a pointer to the element's reference element.  We want to
      90             :   // perform cutting on the reference element such that the quadrature
      91             :   // point locations of the subelements live in the reference
      92             :   // coordinate system, thereby eliminating the need for inverse
      93             :   // mapping.
      94         770 :   const Elem * reference_elem = elem.reference_elem();
      95             : 
      96         385 :   libmesh_assert (reference_elem != nullptr);
      97             : 
      98         770 :   _elem_cutter(*reference_elem, vertex_distance_func);
      99             :   //_elem_cutter(elem, vertex_distance_func);
     100             : 
     101             :   // clear our state & accumulate points from subelements
     102         385 :   _points.clear();
     103         385 :   _weights.clear();
     104             : 
     105             :   // inside subelem
     106             :   {
     107         385 :     const std::vector<Elem const *> & inside_elem (_elem_cutter.inside_elements());
     108             :     // std::cout << inside_elem.size() << " elements inside\n";
     109             : 
     110         770 :     this->add_subelem_values(inside_elem);
     111             :   }
     112             : 
     113             :   // outside subelem
     114             :   {
     115         385 :     const std::vector<Elem const *> & outside_elem (_elem_cutter.outside_elements());
     116             :     // std::cout << outside_elem.size() << " elements outside\n";
     117             : 
     118         770 :     this->add_subelem_values(outside_elem);
     119             :   }
     120             : 
     121             :   // this->print_info();
     122             : }
     123             : 
     124             : 
     125             : 
     126             : template <class QSubCell>
     127        1540 : void QComposite<QSubCell>::add_subelem_values (const std::vector<Elem const *> & subelem)
     128             : 
     129             : {
     130        1540 :   const std::vector<Real>  & subelem_weights = _lagrange_fe->get_JxW();
     131        1540 :   const std::vector<Point> & subelem_points  = _lagrange_fe->get_xyz();
     132             : 
     133        5570 :   for (const auto & elem : subelem)
     134             :     {
     135             :       // tetgen seems to create 0-volume cells on occasion, but we *should*
     136             :       // be catching that appropriately now inside the ElemCutter class.
     137             :       // Just in case trap here, describe the error, and abort.
     138             :       libmesh_try
     139             :         {
     140        4030 :           _lagrange_fe->reinit(elem);
     141        6045 :           _weights.insert(_weights.end(),
     142             :                           subelem_weights.begin(), subelem_weights.end());
     143             : 
     144        6045 :           _points.insert(_points.end(),
     145             :                          subelem_points.begin(), subelem_points.end());
     146             :         }
     147           0 :       libmesh_catch (...)
     148             :         {
     149           0 :           libMesh::err << "ERROR: found a bad cut cell!\n";
     150             : 
     151           0 :           for (auto n : elem->node_index_range())
     152           0 :             libMesh::err << elem->point(n) << std::endl;
     153             : 
     154           0 :           libmesh_error_msg("Tetgen may have created a 0-volume cell during Cutcell integration.");
     155             :         }
     156             :     }
     157        1540 : }
     158             : 
     159             : 
     160             : 
     161             : //--------------------------------------------------------------
     162             : // Explicit instantiations
     163             : template class LIBMESH_EXPORT QComposite<QGauss>;
     164             : template class LIBMESH_EXPORT QComposite<QTrap>;
     165             : template class LIBMESH_EXPORT QComposite<QSimpson>;
     166             : 
     167             : } // namespace libMesh
     168             : 
     169             : #endif // LIBMESH_HAVE_TRIANGLE && LIBMESH_HAVE_TETGEN