Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://mooseframework.inl.gov
       3             : //*
       4             : //* All rights reserved, see COPYRIGHT for full restrictions
       5             : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
       6             : //*
       7             : //* Licensed under LGPL 2.1, please see LICENSE for details
       8             : //* https://www.gnu.org/licenses/lgpl-2.1.html
       9             : 
      10             : #pragma once
      11             : 
      12             : // For MooseIndex
      13             : #include "MooseTypes.h"
      14             : 
      15             : #include "libmesh/libmesh_config.h"
      16             : #include "libmesh/libmesh_common.h"
      17             : #include "libmesh/mesh_base.h"
      18             : #include "libmesh/point.h"
      19             : #include "libmesh/elem.h"
      20             : #ifdef LIBMESH_HAVE_NANOFLANN
      21             : #include "libmesh/nanoflann.hpp"
      22             : #else
      23             : SORRY THIS APPLICATION REQUIRES NANOFLANN
      24             : #endif
      25             : 
      26             : // Using statements. These are used so that the code below does not
      27             : // need to be littered with libMesh:: disambiguations.
      28             : using libMesh::Point;
      29             : using libMesh::MeshBase;
      30             : using libMesh::Real;
      31             : using libMesh::subdomain_id_type;
      32             : using libMesh::dof_id_type;
      33             : using libMesh::Elem;
      34             : 
      35             : /**
      36             :  * This allows us to adapt the MeshBase class for use with nanoflann.
      37             :  *
      38             :  * Taken from detect_slit.cc
      39             :  */
      40             : template <unsigned int Dim>
      41             : class NanoflannMeshAdaptor
      42             : {
      43             : private:
      44             :   // Constant reference to the Mesh we are adapting for use in Nanoflann
      45             :   const MeshBase & _mesh;
      46             : 
      47             : public:
      48             :   NanoflannMeshAdaptor(const MeshBase & mesh) : _mesh(mesh) {}
      49             : 
      50             :   /**
      51             :    * libMesh \p Point coordinate type
      52             :    */
      53             :   typedef Real coord_t;
      54             : 
      55             :   /**
      56             :    * Must return the number of data points
      57             :    */
      58             :   inline size_t kdtree_get_point_count() const { return _mesh.n_nodes(); }
      59             : 
      60             :   /**
      61             :    * Returns the distance between the vector "p1[0:size-1]"
      62             :    * and the data point with index "idx_p2" stored in _mesh
      63             :    */
      64             :   inline coord_t kdtree_distance(const coord_t * p1, const size_t idx_p2, size_t size) const
      65             :   {
      66             :     libmesh_assert_equal_to(size, Dim);
      67             : 
      68             :     // Construct a libmesh Point object from the input coord_t.  This
      69             :     // assumes LIBMESH_DIM==3.
      70             :     Point point1(p1[0], size > 1 ? p1[1] : 0., size > 2 ? p1[2] : 0.);
      71             : 
      72             :     // Get the referred-to point from the Mesh
      73             :     const Point & point2 = _mesh.point(idx_p2);
      74             : 
      75             :     // Compute Euclidean distance, squared
      76             :     return (point1 - point2).norm_sq();
      77             :   }
      78             : 
      79             :   /**
      80             :    * Returns the dim'th component of the idx'th point in the class.
      81             :    */
      82             :   inline coord_t kdtree_get_pt(const size_t idx, int dim) const
      83             :   {
      84             :     libmesh_assert_less(dim, (int)Dim);
      85             :     libmesh_assert_less(idx, _mesh.n_nodes());
      86             :     libmesh_assert_less(dim, 3);
      87             : 
      88             :     return _mesh.point(idx)(dim);
      89             :   }
      90             : 
      91             :   /**
      92             :    * Optional bounding-box computation: return false to default to a
      93             :    * standard bbox computation loop.
      94             :    */
      95             :   template <class BBOX>
      96             :   bool kdtree_get_bbox(BBOX & /* bb */) const
      97             :   {
      98             :     return false;
      99             :   }
     100             : };
     101             : 
     102             : // Useful typedefs for working with NanoflannMeshAdaptors.
     103             : 
     104             : // Declare a type templated on NanoflannMeshAdaptor
     105             : typedef nanoflann::L2_Simple_Adaptor<Real, NanoflannMeshAdaptor<3>> adatper_t;
     106             : 
     107             : // Declare a KDTree type based on NanoflannMeshAdaptor
     108             : typedef nanoflann::KDTreeSingleIndexAdaptor<adatper_t, NanoflannMeshAdaptor<3>, 3> kd_tree_t;
     109             : 
     110             : /**
     111             :  * Special adaptor that works with subdomains of the Mesh.  When
     112             :  * nanoflann asks for the distance between points p1 and p2, and p2 is
     113             :  * not a Point of an element in the required subdomain, a "large"
     114             :  * distance is returned. Likewise, when the coordinates of a point
     115             :  * not in the subdomain are requested, a point at "infinity" is returned.
     116             :  */
     117             : template <unsigned int Dim>
     118             : class NanoflannMeshSubdomainAdaptor
     119             : {
     120             : private:
     121             :   // Constant reference to the Mesh we are adapting for use in Nanoflann
     122             :   const MeshBase & _mesh;
     123             : 
     124             :   // This could be generalized to a std::set of subodmain ids.
     125             :   subdomain_id_type _sid;
     126             : 
     127             :   // Indices of points that are attached to elements in the requested subdomain.
     128             :   std::set<dof_id_type> _legal_point_indices;
     129             : 
     130             : public:
     131        8737 :   NanoflannMeshSubdomainAdaptor(const MeshBase & mesh, subdomain_id_type s) : _mesh(mesh), _sid(s)
     132             :   {
     133             :     // Loop over the elements of the Mesh, for those in the requested
     134             :     // subdomain, add its node ids to the _legal_point_indices set.
     135     1747388 :     for (const auto & elem : _mesh.active_element_ptr_range())
     136      864957 :       if (elem->subdomain_id() == _sid)
     137      143322 :         for (MooseIndex(elem->n_vertices()) n = 0; n < elem->n_vertices(); ++n)
     138       99294 :           _legal_point_indices.insert(elem->node_id(n));
     139        8737 :   }
     140             : 
     141             :   /**
     142             :    * libMesh \p Point coordinate type
     143             :    */
     144             :   typedef Real coord_t;
     145             : 
     146             :   /**
     147             :    * Must return the number of data points
     148             :    */
     149       61159 :   inline size_t kdtree_get_point_count() const { return _mesh.n_nodes(); }
     150             : 
     151             :   /**
     152             :    * Returns the distance between the vector "p1[0:size-1]"
     153             :    * and the data point with index "idx_p2" stored in _mesh
     154             :    */
     155             :   inline coord_t kdtree_distance(const coord_t * p1, const size_t idx_p2, size_t size) const
     156             :   {
     157             :     libmesh_assert_equal_to(size, Dim);
     158             : 
     159             :     // If this is not a valid point, then return a "large" distance.
     160             :     if (!_legal_point_indices.count(static_cast<dof_id_type>(idx_p2)))
     161             :       return std::numeric_limits<coord_t>::max();
     162             : 
     163             :     // Construct a libmesh Point object from the input coord_t.  This
     164             :     // assumes LIBMESH_DIM==3.
     165             :     Point point1(p1[0], size > 1 ? p1[1] : 0., size > 2 ? p1[2] : 0.);
     166             : 
     167             :     // Get the referred-to point from the Mesh
     168             :     const Point & point2 = _mesh.point(idx_p2);
     169             : 
     170             :     // Compute Euclidean distance, squared
     171             :     return (point1 - point2).norm_sq();
     172             :   }
     173             : 
     174             :   /**
     175             :    * Returns the dim'th component of the idx'th point in the class.
     176             :    */
     177    74441270 :   inline coord_t kdtree_get_pt(const size_t idx, int dim) const
     178             :   {
     179             :     libmesh_assert_less(dim, (int)Dim);
     180             :     libmesh_assert_less(idx, _mesh.n_nodes());
     181             :     libmesh_assert_less(dim, 3);
     182             : 
     183             :     // If this is not a valid point, then return a "large" distance.
     184    74441270 :     if (!_legal_point_indices.count(static_cast<dof_id_type>(idx)))
     185    71057122 :       return std::numeric_limits<coord_t>::max();
     186             : 
     187     3384148 :     return _mesh.point(idx)(dim);
     188             :   }
     189             : 
     190             :   /**
     191             :    * Optional bounding-box computation: return false to default to a
     192             :    * standard bbox computation loop.
     193             :    */
     194             :   template <class BBOX>
     195       17474 :   bool kdtree_get_bbox(BBOX & /* bb */) const
     196             :   {
     197       17474 :     return false;
     198             :   }
     199             : };
     200             : 
     201             : // Useful typedefs for working with NanoflannMeshAdaptors.
     202             : 
     203             : // Declare a type templated on NanoflannMeshAdaptor
     204             : typedef nanoflann::L2_Simple_Adaptor<Real, NanoflannMeshSubdomainAdaptor<3>> subdomain_adatper_t;
     205             : 
     206             : // Declare a KDTree type based on NanoflannMeshAdaptor
     207             : typedef nanoflann::
     208             :     KDTreeSingleIndexAdaptor<subdomain_adatper_t, NanoflannMeshSubdomainAdaptor<3>, 3>
     209             :         subdomain_kd_tree_t;