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             : #include "ThreadedRadialAverageLoop.h"
      11             : #include "Function.h"
      12             : 
      13             : // Make newer nanoflann API spelling compatible with older nanoflann
      14             : // versions
      15             : #if NANOFLANN_VERSION < 0x150
      16             : namespace nanoflann
      17             : {
      18             : typedef SearchParams SearchParameters;
      19             : 
      20             : template <typename T, typename U>
      21             : using ResultItem = std::pair<T, U>;
      22             : }
      23             : #endif
      24             : 
      25         120 : ThreadedRadialAverageLoop::ThreadedRadialAverageLoop(RadialAverage & green) : _radavg(green) {}
      26             : 
      27             : // Splitting Constructor
      28          12 : ThreadedRadialAverageLoop::ThreadedRadialAverageLoop(const ThreadedRadialAverageLoop & x,
      29          12 :                                                      Threads::split /*split*/)
      30          12 :   : _radavg(x._radavg)
      31             : {
      32          12 : }
      33             : 
      34             : void
      35         132 : ThreadedRadialAverageLoop::operator()(const QPDataRange & qpdata_range)
      36             : {
      37             :   // fetch data from parent
      38         132 :   const auto radius = _radavg._radius;
      39         132 :   const auto & qp_data = _radavg._qp_data;
      40         132 :   const auto & kd_tree = _radavg._kd_tree;
      41         132 :   const auto & weights_type = _radavg._weights_type;
      42             : 
      43             :   // tree search data structures
      44         132 :   std::vector<nanoflann::ResultItem<std::size_t, Real>> ret_matches;
      45         132 :   nanoflann::SearchParameters search_params;
      46             : 
      47             :   // result map entry
      48         132 :   const auto end_it = _radavg._average.end();
      49         132 :   auto it = end_it;
      50             : 
      51             :   // iterate over qp range
      52       75732 :   for (auto && local_qp : qpdata_range)
      53             :   {
      54             :     // Look up result map iterator only if we enter a new element. this saves a bunch
      55             :     // of map lookups because same element entries are consecutive in the qp_data vector.
      56       75600 :     if (it == end_it || it->first != local_qp._elem_id)
      57        8400 :       it = _radavg._average.find(local_qp._elem_id);
      58             : 
      59             :     // initialize result entry
      60             :     mooseAssert(it != end_it, "Current element id not found in result set.");
      61       75600 :     auto & sum = it->second[local_qp._qp];
      62       75600 :     sum = 0.0;
      63             : 
      64       75600 :     ret_matches.clear();
      65             :     std::size_t n_result =
      66       75600 :         kd_tree->radiusSearch(&(local_qp._q_point(0)), radius * radius, ret_matches, search_params);
      67       75600 :     Real total_vol = 0.0;
      68       75600 :     Real weight = 1.0;
      69     4288032 :     for (std::size_t j = 0; j < n_result; ++j)
      70             :     {
      71     4212432 :       const auto & other_qp = qp_data[ret_matches[j].first];
      72     4212432 :       switch (weights_type)
      73             :       {
      74     1404144 :         case RadialAverage::WeightsType::CONSTANT:
      75     1404144 :           break;
      76             : 
      77     1404144 :         case RadialAverage::WeightsType::LINEAR:
      78     1404144 :           weight = radius - std::sqrt(ret_matches[j].second);
      79     1404144 :           break;
      80             : 
      81     1404144 :         case RadialAverage::WeightsType::COSINE:
      82     1404144 :           weight = std::cos(std::sqrt(ret_matches[j].second) / radius * libMesh::pi) + 1.0;
      83     1404144 :           break;
      84             :       }
      85             : 
      86     4212432 :       sum += other_qp._value * other_qp._volume * weight;
      87     4212432 :       total_vol += other_qp._volume * weight;
      88             :     }
      89       75600 :     sum /= total_vol;
      90             :   }
      91         132 : }