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 "DiracKernelInfo.h"
      11             : #include "MooseMesh.h"
      12             : #include "MooseEnum.h"
      13             : #include "DiracKernelBase.h"
      14             : #include "MooseBase.h"
      15             : 
      16             : // LibMesh
      17             : #include "libmesh/point_locator_base.h"
      18             : #include "libmesh/elem.h"
      19             : #include "libmesh/enum_point_locator_type.h"
      20             : #include "libmesh/point.h"
      21             : 
      22             : using namespace libMesh;
      23             : 
      24       64619 : DiracKernelInfo::DiracKernelInfo()
      25       64619 :   : _point_locator(), _point_equal_distance_sq(libMesh::TOLERANCE * libMesh::TOLERANCE)
      26             : {
      27       64619 : }
      28             : 
      29       61479 : DiracKernelInfo::~DiracKernelInfo() {}
      30             : 
      31             : void
      32      711980 : DiracKernelInfo::addPoint(const Elem * elem, const Point & p, const Real & value)
      33             : {
      34      711980 :   _elements.insert(elem);
      35             : 
      36      711980 :   std::pair<std::vector<Point>, std::vector<Real>> & multi_point_list = _points[elem];
      37             : 
      38      711980 :   const unsigned int npoint = multi_point_list.first.size();
      39             :   mooseAssert(npoint == multi_point_list.second.size(),
      40             :               "Different sizes for location and point value data");
      41             : 
      42      731540 :   for (unsigned int i = 0; i < npoint; ++i)
      43       63828 :     if (pointsFuzzyEqual(multi_point_list.first[i], p))
      44             :     {
      45             :       // A point at the same (within a tolerance) location as p exists, accumulate its value.
      46       44268 :       multi_point_list.second[i] += value;
      47       44268 :       return;
      48             :     }
      49             : 
      50             :   // No prior point found at this location, add it with its initial value.
      51      667712 :   multi_point_list.first.push_back(p);
      52      667712 :   multi_point_list.second.push_back(value);
      53             : }
      54             : 
      55             : void
      56     3703355 : DiracKernelInfo::clearPoints()
      57             : {
      58     3703355 :   _elements.clear();
      59     3703355 :   _points.clear();
      60     3703355 : }
      61             : 
      62             : bool
      63      315876 : DiracKernelInfo::hasPoint(const Elem * elem, const Point & p)
      64             : {
      65      315876 :   std::vector<Point> & point_list = _points[elem].first;
      66             : 
      67      324576 :   for (const auto & pt : point_list)
      68      324576 :     if (pointsFuzzyEqual(pt, p))
      69      315876 :       return true;
      70             : 
      71             :   // If we haven't found it, we don't have it.
      72           0 :   return false;
      73             : }
      74             : 
      75             : void
      76      155211 : DiracKernelInfo::updatePointLocator(const MooseMesh & mesh)
      77             : {
      78             :   // Note: we could update the PointLocator *every* time we call this
      79             :   // function, but that may introduce an unacceptable overhead in
      80             :   // problems which don't need a PointLocator at all.  This issue will
      81             :   // most likely become a moot point when we eventually add a shared
      82             :   // "CachingPointLocator" to MOOSE.
      83             :   // _point_locator = PointLocatorBase::build(TREE_LOCAL_ELEMENTS, mesh);
      84             : 
      85             :   // Construct the PointLocator object if *any* processors have Dirac
      86             :   // points.  Note: building a PointLocator object is a parallel_only()
      87             :   // function, so this is an all-or-nothing thing.
      88      155211 :   unsigned pl_needs_rebuild = _elements.size();
      89      155211 :   mesh.comm().max(pl_needs_rebuild);
      90             : 
      91      155211 :   if (pl_needs_rebuild)
      92             :   {
      93             :     // PointLocatorBase::build() is a parallel_only function!  So we
      94             :     // can't skip building it just becuase our local _elements is
      95             :     // empty, it might be non-empty on some other processor!
      96        2423 :     _point_locator = PointLocatorBase::build(TREE_LOCAL_ELEMENTS, mesh);
      97             : 
      98             :     // We may be querying for points which are not in the semilocal
      99             :     // part of a distributed mesh.
     100        2423 :     _point_locator->enable_out_of_mesh_mode();
     101             :   }
     102             :   else
     103             :   {
     104             :     // There are no elements with Dirac points, but we have been
     105             :     // requested to update the PointLocator so we have to assume the
     106             :     // old one is invalid.  Therefore we reset it to NULL... however
     107             :     // adding this line causes the code to hang because it triggers
     108             :     // the PointLocator to be rebuilt in a non-parallel-only segment
     109             :     // of the code later... so it's commented out for now even though
     110             :     // it's probably the right behavior.
     111             :     // _point_locator.reset(NULL);
     112             :   }
     113      155211 : }
     114             : 
     115             : const Elem *
     116       39344 : DiracKernelInfo::findPoint(const Point & p,
     117             :                            const MooseMesh & mesh,
     118             :                            const std::set<SubdomainID> & blocks,
     119             :                            const PointNotFoundBehavior point_not_found_behavior,
     120             :                            const MooseBase & consumer)
     121             : {
     122             :   // If the PointLocator has never been created, do so now.  NOTE - WE
     123             :   // CAN'T DO THIS if findPoint() is only called on some processors,
     124             :   // PointLocatorBase::build() is a 'parallel_only' method!
     125       39344 :   if (_point_locator.get() == NULL)
     126             :   {
     127         567 :     _point_locator = PointLocatorBase::build(TREE_LOCAL_ELEMENTS, mesh);
     128         567 :     _point_locator->enable_out_of_mesh_mode();
     129             :   }
     130             : 
     131             :   // Check that the PointLocator is ready to start locating points.
     132             :   // So far I do not have any tests that trip this...
     133       39344 :   if (_point_locator->initialized() == false)
     134           0 :     consumer.mooseError("PointLocator is not initialized!");
     135             : 
     136             :   // Note: The PointLocator object returns NULL when the Point is not
     137             :   // found within the Mesh.  This is not considered to be an error as
     138             :   // far as the DiracKernels are concerned: sometimes the Mesh moves
     139             :   // out from the Dirac point entirely and in that case the Point just
     140             :   // gets "deactivated".
     141       39344 :   const Elem * elem = (*_point_locator)(p, &blocks);
     142             : 
     143             :   // The processors may not agree on which Elem the point is in.  This
     144             :   // can happen if a Dirac point lies on the processor boundary, and
     145             :   // two or more neighboring processors think the point is in the Elem
     146             :   // on *their* side.
     147       39344 :   dof_id_type elem_id = elem ? elem->id() : DofObject::invalid_id;
     148             : 
     149             :   // We are going to let the element with the smallest ID "win", all other
     150             :   // procs will return NULL.
     151       39344 :   dof_id_type min_elem_id = elem_id;
     152       39344 :   mesh.comm().min(min_elem_id);
     153             : 
     154       39344 :   if (min_elem_id == DofObject::invalid_id)
     155             :   {
     156        3169 :     std::stringstream msg;
     157       12676 :     msg << "Point " << p << " not found in block(s) " << Moose::stringify(blocks, ", ") << ".\n";
     158        3169 :     switch (point_not_found_behavior)
     159             :     {
     160           3 :       case PointNotFoundBehavior::ERROR:
     161           3 :         consumer.mooseError(msg.str());
     162             :         break;
     163         298 :       case PointNotFoundBehavior::WARNING:
     164         298 :         mooseDoOnce(consumer.mooseWarning(msg.str() + "This message will not be repeated."));
     165         298 :         break;
     166        2868 :       case PointNotFoundBehavior::IGNORE:
     167        2868 :         break;
     168           0 :       default:
     169           0 :         consumer.mooseError("Internal enum error.");
     170             :     }
     171        3166 :   }
     172             : 
     173             :   // But we notably need the processor which owns elem_id to return it!
     174       39341 :   if (min_elem_id != DofObject::invalid_id)
     175       36175 :     if (const auto min_elem = mesh.queryElemPtr(min_elem_id);
     176       36175 :         min_elem && min_elem->processor_id() == mesh.processor_id())
     177       26406 :       return min_elem;
     178       12935 :   return nullptr;
     179             : }
     180             : 
     181             : bool
     182      388404 : DiracKernelInfo::pointsFuzzyEqual(const Point & a, const Point & b)
     183             : {
     184      388404 :   const Real dist_sq = (a - b).norm_sq();
     185      388404 :   return dist_sq < _point_equal_distance_sq;
     186             : }