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 "PointSamplerBase.h"
      11             : 
      12             : // MOOSE includes
      13             : #include "MooseMesh.h"
      14             : #include "Assembly.h"
      15             : 
      16             : #include "libmesh/mesh_tools.h"
      17             : 
      18             : InputParameters
      19       14289 : PointSamplerBase::validParams()
      20             : {
      21       14289 :   InputParameters params = GeneralVectorPostprocessor::validParams();
      22             : 
      23       14289 :   params += SamplerBase::validParams();
      24             : 
      25       42867 :   params.addParam<PostprocessorName>(
      26       28578 :       "scaling", 1.0, "The postprocessor that the variables are multiplied with");
      27       42867 :   params.addParam<bool>(
      28             :       "warn_discontinuous_face_values",
      29       28578 :       true,
      30             :       "Whether to return a warning if a discontinuous variable is sampled on a face");
      31             : 
      32       14289 :   return params;
      33           0 : }
      34             : 
      35         731 : PointSamplerBase::PointSamplerBase(const InputParameters & parameters)
      36             :   : GeneralVectorPostprocessor(parameters),
      37             :     SamplerBase(parameters, this, _communicator),
      38        1462 :     _mesh(_subproblem.mesh()),
      39         731 :     _pp_value(getPostprocessorValue("scaling")),
      40         731 :     _warn_discontinuous_face_values(getParam<bool>("warn_discontinuous_face_values")),
      41         731 :     _discontinuous_at_faces(false)
      42             : {
      43         731 : }
      44             : 
      45             : void
      46        3866 : PointSamplerBase::initialize()
      47             : {
      48             :   // Generate new Ids if the point vector has grown (non-negative counting numbers)
      49        3866 :   if (_points.size() > _ids.size())
      50             :   {
      51          11 :     auto old_size = _ids.size();
      52          11 :     _ids.resize(_points.size());
      53          11 :     std::iota(_ids.begin() + old_size, _ids.end(), old_size);
      54             :   }
      55             :   // Otherwise sync the ids array to be smaller if the point vector has been shrunk
      56        3855 :   else if (_points.size() < _ids.size())
      57           0 :     _ids.resize(_points.size());
      58             : 
      59        3866 :   SamplerBase::initialize();
      60             : 
      61             :   // We do this here just in case it's been destroyed and recreated because of mesh adaptivity.
      62        3866 :   _pl = _mesh.getPointLocator();
      63             : 
      64             :   // We may not find a requested point on a distributed mesh, and
      65             :   // that's okay.
      66        3866 :   _pl->enable_out_of_mesh_mode();
      67             : 
      68             :   // Reset the point arrays
      69        3866 :   _found_points.assign(_points.size(), false);
      70             : 
      71        3866 :   _point_values.resize(_points.size());
      72        3866 :   std::for_each(
      73       57666 :       _point_values.begin(), _point_values.end(), [](std::vector<Real> & vec) { vec.clear(); });
      74        3866 : }
      75             : 
      76             : void
      77        3866 : PointSamplerBase::finalize()
      78             : {
      79             :   // Save off for speed
      80        3866 :   const auto pid = processor_id();
      81             : 
      82             :   // Consolidate _found_points across processes to know which points were found
      83        3866 :   auto _global_found_points = _found_points;
      84        3866 :   _comm.sum(_global_found_points);
      85             : 
      86             :   // Keep track of maximum process ids for each point to only add it once
      87        3866 :   std::vector<unsigned int> max_pid(_found_points.size());
      88        3866 :   _comm.maxloc(_found_points, max_pid);
      89             : 
      90       61528 :   for (MooseIndex(_found_points) i = 0; i < _found_points.size(); ++i)
      91             :   {
      92             :     // _global_found_points should contain all 1's at this point (ie every point was found by a
      93             :     // proc)
      94       57666 :     if (pid == 0 && !_global_found_points[i])
      95           4 :       mooseError("In ", name(), ", sample point not found: ", _points[i]);
      96             : 
      97             :     // only process that found the point has the value, and only process with max id should add
      98       57662 :     if (pid == max_pid[i] && _found_points[i])
      99       40345 :       SamplerBase::addSample(_points[i], _ids[i], _point_values[i]);
     100             :   }
     101             : 
     102        3862 :   SamplerBase::finalize();
     103        3862 : }
     104             : 
     105             : const Elem *
     106       42895 : PointSamplerBase::getLocalElemContainingPoint(const Point & p)
     107             : {
     108       42895 :   const Elem * elem = nullptr;
     109       42895 :   if (_discontinuous_at_faces)
     110             :   {
     111             :     libmesh_parallel_only(comm());
     112             : 
     113             :     // Get all possible elements the point may be in
     114        4770 :     std::set<const Elem *> candidate_elements;
     115        4770 :     (*_pl)(p, candidate_elements);
     116             : 
     117             :     // Look at all the element IDs
     118        4770 :     std::set<dof_id_type> candidate_ids;
     119        9477 :     for (auto candidate : candidate_elements)
     120        4707 :       candidate_ids.insert(candidate->id());
     121             : 
     122        4770 :     comm().set_union(candidate_ids);
     123             : 
     124             :     // Domains without candidate elements will not own the lowest ID one
     125        4770 :     if (candidate_elements.size())
     126             :     {
     127             :       // If we know of the minimum, this will be valid. Otherwise, it'll be
     128             :       // nullptr which is fine
     129        4533 :       elem = _mesh.queryElemPtr(*(candidate_ids.begin()));
     130             : 
     131             :       // Print a warning if it's on a face and a variable is discontinuous
     132        4533 :       if (_warn_discontinuous_face_values && candidate_ids.size() > 1)
     133          18 :         mooseDoOnce(mooseWarning("A discontinuous variable is sampled on a face, at ", p));
     134             :     }
     135        4770 :   }
     136             :   else // continuous variables
     137             :   {
     138             :     // Get all possible elements the point may be in
     139             :     // We cant just locate in one element because at the edge between two process domains, it could
     140             :     // be that both domains find the element that is not within their domain
     141       38125 :     std::set<const Elem *> candidate_elements;
     142       38125 :     (*_pl)(p, candidate_elements);
     143             : 
     144             :     // Only keep the one that may be local
     145      103253 :     for (auto candidate : candidate_elements)
     146       65128 :       if (candidate->processor_id() == processor_id())
     147       61630 :         elem = candidate;
     148       38125 :   }
     149             : 
     150       42895 :   if (elem && elem->processor_id() == processor_id())
     151       40836 :     return elem;
     152             : 
     153        2059 :   return nullptr;
     154             : }