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             : #include "ElementUserObject.h"
      13             : #include "DataIO.h"
      14             : #include "PointListAdaptor.h"
      15             : #include "Function.h"
      16             : 
      17             : #include "libmesh/data_type.h"
      18             : 
      19             : #include <set>
      20             : #include <map>
      21             : #include <vector>
      22             : #include <memory>
      23             : 
      24             : class ThreadedRadialAverageLoop;
      25             : 
      26             : /**
      27             :  * Gather and communicate a full list of all quadrature points and the values of
      28             :  * a selected material property at each point. Use a KD-Tree to get the weighted spatial
      29             :  * average of a material property. This code borrows heavily from
      30             :  * RadialGreensConvolution in MAGPIE. This code does not include support for
      31             :  * periodic BCs, but RadialGreensConvolution shows how that can be supported.
      32             :  */
      33             : class RadialAverage : public ElementUserObject
      34             : {
      35             : public:
      36             :   static InputParameters validParams();
      37             : 
      38             :   RadialAverage(const InputParameters & parameters);
      39             : 
      40             :   virtual void initialSetup() override;
      41             : 
      42             :   virtual void initialize() override;
      43             :   virtual void execute() override;
      44             :   virtual void finalize() override;
      45             :   virtual void threadJoin(const UserObject & y) override;
      46             :   virtual void meshChanged() override;
      47             : 
      48             :   /// quadrature point data
      49             :   struct QPData
      50             :   {
      51             :     /// physical coordinates of the quadrature point
      52             :     Point _q_point;
      53             :     /// element id
      54             :     dof_id_type _elem_id;
      55             :     /// index of the quadrature point
      56             :     short _qp;
      57             :     /// current value * _JxW
      58             :     Real _volume;
      59             :     /// variable value
      60             :     Real _value;
      61             : 
      62        9714 :     QPData() : _q_point(), _elem_id(libMesh::invalid_uint), _qp(0), _volume(0.0), _value(0.0) {}
      63       75600 :     QPData(const Point & q_point, dof_id_type elem_id, short qp, Real volume, Real value)
      64       75600 :       : _q_point(q_point), _elem_id(elem_id), _qp(qp), _volume(volume), _value(value)
      65             :     {
      66       75600 :     }
      67             :   };
      68             : 
      69             :   using Result = std::map<dof_id_type, std::vector<Real>>;
      70          33 :   const Result & getAverage() const { return _average; }
      71             : 
      72             : protected:
      73             :   void updateCommunicationLists();
      74             : 
      75             :   /// distance based weight function
      76             :   enum class WeightsType
      77             :   {
      78             :     CONSTANT,
      79             :     LINEAR,
      80             :     COSINE
      81             :   } _weights_type;
      82             : 
      83             :   /// material to be averaged
      84             :   const MaterialProperty<Real> & _prop;
      85             : 
      86             :   /// cut-off radius
      87             :   const Real _radius;
      88             : 
      89             :   /// communication padding distance
      90             :   const Real _padding;
      91             : 
      92             :   /// gathered data
      93             :   std::vector<QPData> _qp_data;
      94             : 
      95             :   /// average result
      96             :   Result _average;
      97             : 
      98             :   using KDTreeType = nanoflann::KDTreeSingleIndexAdaptor<
      99             :       nanoflann::L2_Simple_Adaptor<Real, PointListAdaptor<QPData>>,
     100             :       PointListAdaptor<QPData>,
     101             :       LIBMESH_DIM,
     102             :       std::size_t>;
     103             : 
     104             :   /// spatial index (nanoflann guarantees this to be threadsafe under read-only operations)
     105             :   std::unique_ptr<KDTreeType> _kd_tree;
     106             : 
     107             :   /// The data structure used to find neighboring elements give a node ID
     108             :   std::vector<std::vector<const Elem *>> _nodes_to_elem_map;
     109             : 
     110             :   /// set of nodes on the boundary of the current processor domain
     111             :   std::set<Point> _boundary_nodes;
     112             : 
     113             :   /// set of all _qp_data indices that are within _radius of any _boundary_nodes
     114             :   std::set<std::size_t> _boundary_data_indices;
     115             : 
     116             :   /// QPData indices to send to the various processors
     117             :   std::vector<std::set<std::size_t>> _communication_lists;
     118             :   bool _update_communication_lists;
     119             : 
     120             :   /// processors to send (potentially empty) data to
     121             :   std::vector<processor_id_type> _candidate_procs;
     122             : 
     123             :   processor_id_type _my_pid;
     124             : 
     125             :   //@{ PerfGraph identifiers
     126             :   PerfID _perf_meshchanged;
     127             :   PerfID _perf_updatelists;
     128             :   PerfID _perf_finalize;
     129             :   //@}
     130             : 
     131             :   friend class ThreadedRadialAverageLoop;
     132             : };
     133             : 
     134             : namespace TIMPI
     135             : {
     136             : 
     137             : template <>
     138             : class StandardType<RadialAverage::QPData> : public DataType
     139             : {
     140             : public:
     141             :   explicit StandardType(const RadialAverage::QPData * example = nullptr);
     142             :   StandardType(const StandardType<RadialAverage::QPData> & t);
     143         216 :   ~StandardType() { this->free(); }
     144             : };
     145             : 
     146             : } // namespace TIMPI