Line data    Source code

       1             : /**********************************************************************/
       2             : /*                     DO NOT MODIFY THIS HEADER                      */
       3             : /* MAGPIE - Mesoscale Atomistic Glue Program for Integrated Execution */
       4             : /*                                                                    */
       5             : /*            Copyright 2017 Battelle Energy Alliance, LLC            */
       6             : /*                        ALL RIGHTS RESERVED                         */
       7             : /**********************************************************************/
       8             : 
       9             : #pragma once
      10             : 
      11             : #include "ElementUserObject.h"
      12             : #include "DataIO.h"
      13             : 
      14             : #include "libmesh/data_type.h"
      15             : 
      16             : #include <set>
      17             : #include <map>
      18             : #include <vector>
      19             : #include <array>
      20             : #include <memory>
      21             : #include <tuple>
      22             : 
      23             : using namespace TIMPI;
      24             : 
      25             : class ThreadedRadialGreensConvolutionLoop;
      26             : 
      27             : /**
      28             :  * Gather and communicate a full list of all quadrature points and the values of
      29             :  * a selected variable at each point. Use a KD-Tree to integrate the weighted
      30             :  * neighborhood of each QP to obtain the convolution.
      31             :  */
      32             : class RadialGreensConvolution : public ElementUserObject
      33             : {
      34             : public:
      35             :   static InputParameters validParams();
      36             : 
      37             :   RadialGreensConvolution(const InputParameters & parameters);
      38             : 
      39             :   virtual void initialSetup() override;
      40             : 
      41             :   virtual void initialize() override;
      42             :   virtual void execute() override;
      43             :   virtual void finalize() override;
      44             :   virtual void threadJoin(const UserObject & y) override;
      45             :   virtual void meshChanged() override;
      46             : 
      47             :   /// quaddrature point data
      48             :   struct QPData
      49             :   {
      50             :     /// physical coordinates of the quadrature point
      51             :     Point _q_point;
      52             :     /// element id
      53             :     dof_id_type _elem_id;
      54             :     /// index of the quadrature point
      55             :     short _qp;
      56             :     /// current value * _JxW
      57             :     Real _volume;
      58             :     /// variable value
      59             :     Real _value;
      60             : 
      61       15202 :     QPData() : _q_point(), _elem_id(libMesh::invalid_uint), _qp(0), _volume(0.0), _value(0.0) {}
      62             :     QPData(const Point & q_point, dof_id_type elem_id, short qp, Real volume, Real value)
      63      142920 :       : _q_point(q_point), _elem_id(elem_id), _qp(qp), _volume(volume), _value(value)
      64             :     {
      65             :     }
      66             :   };
      67             : 
      68             :   using Result = std::map<dof_id_type, std::vector<Real>>;
      69          40 :   const Result & getConvolution() const { return _convolution; }
      70             : 
      71             : protected:
      72             :   Real attenuationIntegral(Real h1, Real h2, Real r, unsigned int dim) const;
      73             :   void insertNotLocalPointNeighbors(dof_id_type);
      74             :   void insertNotLocalPeriodicPointNeighbors(dof_id_type, const Node *);
      75             :   void findNotLocalPeriodicPointNeighbors(const Node *);
      76             :   void updateCommunicationLists();
      77             : 
      78             :   /// variable field to be gathered
      79             :   const VariableValue & _v;
      80             : 
      81             :   /// index of field variable
      82             :   unsigned int _v_var;
      83             : 
      84             :   /// Green's function
      85             :   const Function & _function;
      86             : 
      87             :   /// Green's function cut-off radius
      88             :   const Real _r_cut;
      89             : 
      90             :   /// Normalize the Green's function to one to make the integral of the convolution
      91             :   /// the same as the integral of the original data.
      92             :   const bool _normalize;
      93             : 
      94             :   /// mesh dimension
      95             :   unsigned int _dim;
      96             : 
      97             :   /// greens function integral table for correction of lower dimensional convolutions
      98             :   std::vector<Real> _correction_integral;
      99             : 
     100             :   /// gathered data
     101             :   std::vector<QPData> _qp_data;
     102             : 
     103             :   /// convolution result
     104             :   Result _convolution;
     105             : 
     106             :   /// is the mesh translated periodic in a given cardinal direction
     107             :   std::array<bool, LIBMESH_DIM> _periodic;
     108             : 
     109             :   ///@{ periodic size per component
     110             :   std::array<Real, LIBMESH_DIM> _periodic_min;
     111             :   std::array<Real, LIBMESH_DIM> _periodic_max;
     112             :   std::array<Point, LIBMESH_DIM> _periodic_vector;
     113             :   ///@}
     114             : 
     115             :   using KDTreeType = nanoflann::KDTreeSingleIndexAdaptor<
     116             :       nanoflann::L2_Simple_Adaptor<Real, PointListAdaptor<QPData>>,
     117             :       PointListAdaptor<QPData>,
     118             :       LIBMESH_DIM,
     119             :       std::size_t>;
     120             : 
     121             :   /// spatial index (nanoflann guarantees this to be threadsafe under read-only operations)
     122             :   std::unique_ptr<KDTreeType> _kd_tree;
     123             : 
     124             :   Real _zero_dh;
     125             : 
     126             :   /// DOF map
     127             :   const DofMap & _dof_map;
     128             : 
     129             :   /// A pointer to the periodic boundary constraints object
     130             :   PeriodicBoundaries * _pbs;
     131             : 
     132             :   /// PointLocator for finding topological neighbors
     133             :   std::unique_ptr<PointLocatorBase> _point_locator;
     134             : 
     135             :   /**
     136             :    * The data structure which is a list of nodes that are constrained to other nodes
     137             :    * based on the imposed periodic boundary conditions.
     138             :    */
     139             :   std::multimap<dof_id_type, dof_id_type> _periodic_node_map;
     140             : 
     141             :   /// The data structure used to find neighboring elements give a node ID
     142             :   std::unordered_map<dof_id_type, std::vector<const Elem *>> _nodes_to_elem_map;
     143             : 
     144             :   // list of direct point neighbor elements of the current processor domain
     145             :   std::set<const Elem *> _point_neighbors;
     146             : 
     147             :   // list of periodic point neighbor elements of the current processor domain
     148             :   std::set<std::tuple<const Elem *, const Node *, const Node *>> _periodic_point_neighbors;
     149             : 
     150             :   /// QPData indices to send to the various processors
     151             :   std::vector<std::set<std::size_t>> _communication_lists;
     152             :   bool _update_communication_lists;
     153             : 
     154             :   processor_id_type _my_pid;
     155             : 
     156             :   //@{ PerfGraph identifiers
     157             :   PerfID _perf_meshchanged;
     158             :   PerfID _perf_updatelists;
     159             :   PerfID _perf_finalize;
     160             :   //@}
     161             : 
     162             :   friend class ThreadedRadialGreensConvolutionLoop;
     163             : };
     164             : 
     165             : template <>
     166             : const Point & PointListAdaptor<RadialGreensConvolution::QPData>::getPoint(
     167             :     const RadialGreensConvolution::QPData & item) const;
     168             : 
     169             : namespace TIMPI
     170             : {
     171             : 
     172             : template <>
     173             : class StandardType<RadialGreensConvolution::QPData> : public DataType
     174             : {
     175             : public:
     176             :   explicit StandardType(const RadialGreensConvolution::QPData * example = nullptr);
     177             :   StandardType(const StandardType<RadialGreensConvolution::QPData> & t);
     178          36 :   ~StandardType() { this->free(); }
     179             : };
     180             : } // namespace TIMPI