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             : 
      13             : // LibMesh
      14             : #include "libmesh/point_locator_base.h"
      15             : #include "libmesh/elem.h"
      16             : #include "libmesh/enum_point_locator_type.h"
      17             : #include "libmesh/point.h"
      18             : 
      19             : using namespace libMesh;
      20             : 
      21       60987 : DiracKernelInfo::DiracKernelInfo()
      22       60987 :   : _point_locator(), _point_equal_distance_sq(libMesh::TOLERANCE * libMesh::TOLERANCE)
      23             : {
      24       60987 : }
      25             : 
      26       56616 : DiracKernelInfo::~DiracKernelInfo() {}
      27             : 
      28             : void
      29      710602 : DiracKernelInfo::addPoint(const Elem * elem, const Point & p)
      30             : {
      31      710602 :   _elements.insert(elem);
      32             : 
      33      710602 :   std::pair<std::vector<Point>, std::vector<unsigned int>> & multi_point_list = _points[elem];
      34             : 
      35      710602 :   const unsigned int npoint = multi_point_list.first.size();
      36             :   mooseAssert(npoint == multi_point_list.second.size(),
      37             :               "Different sizes for location and multiplicity data");
      38             : 
      39      730042 :   for (unsigned int i = 0; i < npoint; ++i)
      40       63215 :     if (pointsFuzzyEqual(multi_point_list.first[i], p))
      41             :     {
      42             :       // a point at the same (within a tolerance) location as p exists, increase its multiplicity
      43       43775 :       multi_point_list.second[i]++;
      44       43775 :       return;
      45             :     }
      46             : 
      47             :   // no prior point found at this location, add it with a multiplicity of one
      48      666827 :   multi_point_list.first.push_back(p);
      49      666827 :   multi_point_list.second.push_back(1);
      50             : }
      51             : 
      52             : void
      53     3715306 : DiracKernelInfo::clearPoints()
      54             : {
      55     3715306 :   _elements.clear();
      56     3715306 :   _points.clear();
      57     3715306 : }
      58             : 
      59             : bool
      60      315754 : DiracKernelInfo::hasPoint(const Elem * elem, const Point & p)
      61             : {
      62      315754 :   std::vector<Point> & point_list = _points[elem].first;
      63             : 
      64      324394 :   for (const auto & pt : point_list)
      65      324394 :     if (pointsFuzzyEqual(pt, p))
      66      315754 :       return true;
      67             : 
      68             :   // If we haven't found it, we don't have it.
      69           0 :   return false;
      70             : }
      71             : 
      72             : void
      73      156844 : DiracKernelInfo::updatePointLocator(const MooseMesh & mesh)
      74             : {
      75             :   // Note: we could update the PointLocator *every* time we call this
      76             :   // function, but that may introduce an unacceptable overhead in
      77             :   // problems which don't need a PointLocator at all.  This issue will
      78             :   // most likely become a moot point when we eventually add a shared
      79             :   // "CachingPointLocator" to MOOSE.
      80             :   // _point_locator = PointLocatorBase::build(TREE_LOCAL_ELEMENTS, mesh);
      81             : 
      82             :   // Construct the PointLocator object if *any* processors have Dirac
      83             :   // points.  Note: building a PointLocator object is a parallel_only()
      84             :   // function, so this is an all-or-nothing thing.
      85      156844 :   unsigned pl_needs_rebuild = _elements.size();
      86      156844 :   mesh.comm().max(pl_needs_rebuild);
      87             : 
      88      156844 :   if (pl_needs_rebuild)
      89             :   {
      90             :     // PointLocatorBase::build() is a parallel_only function!  So we
      91             :     // can't skip building it just becuase our local _elements is
      92             :     // empty, it might be non-empty on some other processor!
      93        2379 :     _point_locator = PointLocatorBase::build(TREE_LOCAL_ELEMENTS, mesh);
      94             : 
      95             :     // We may be querying for points which are not in the semilocal
      96             :     // part of a distributed mesh.
      97        2379 :     _point_locator->enable_out_of_mesh_mode();
      98             :   }
      99             :   else
     100             :   {
     101             :     // There are no elements with Dirac points, but we have been
     102             :     // requested to update the PointLocator so we have to assume the
     103             :     // old one is invalid.  Therefore we reset it to NULL... however
     104             :     // adding this line causes the code to hang because it triggers
     105             :     // the PointLocator to be rebuilt in a non-parallel-only segment
     106             :     // of the code later... so it's commented out for now even though
     107             :     // it's probably the right behavior.
     108             :     // _point_locator.reset(NULL);
     109             :   }
     110      156844 : }
     111             : 
     112             : const Elem *
     113       41247 : DiracKernelInfo::findPoint(const Point & p,
     114             :                            const MooseMesh & mesh,
     115             :                            const std::set<SubdomainID> & blocks)
     116             : {
     117             :   // If the PointLocator has never been created, do so now.  NOTE - WE
     118             :   // CAN'T DO THIS if findPoint() is only called on some processors,
     119             :   // PointLocatorBase::build() is a 'parallel_only' method!
     120       41247 :   if (_point_locator.get() == NULL)
     121             :   {
     122         543 :     _point_locator = PointLocatorBase::build(TREE_LOCAL_ELEMENTS, mesh);
     123         543 :     _point_locator->enable_out_of_mesh_mode();
     124             :   }
     125             : 
     126             :   // Check that the PointLocator is ready to start locating points.
     127             :   // So far I do not have any tests that trip this...
     128       41247 :   if (_point_locator->initialized() == false)
     129           0 :     mooseError("Error, PointLocator is not initialized!");
     130             : 
     131             :   // Note: The PointLocator object returns NULL when the Point is not
     132             :   // found within the Mesh.  This is not considered to be an error as
     133             :   // far as the DiracKernels are concerned: sometimes the Mesh moves
     134             :   // out from the Dirac point entirely and in that case the Point just
     135             :   // gets "deactivated".
     136       41247 :   const Elem * elem = (*_point_locator)(p, &blocks);
     137             : 
     138             :   // The processors may not agree on which Elem the point is in.  This
     139             :   // can happen if a Dirac point lies on the processor boundary, and
     140             :   // two or more neighboring processors think the point is in the Elem
     141             :   // on *their* side.
     142       41247 :   dof_id_type elem_id = elem ? elem->id() : DofObject::invalid_id;
     143             : 
     144             :   // We are going to let the element with the smallest ID "win", all other
     145             :   // procs will return NULL.
     146       41247 :   dof_id_type min_elem_id = elem_id;
     147       41247 :   mesh.comm().min(min_elem_id);
     148             : 
     149       41247 :   return min_elem_id == elem_id ? elem : NULL;
     150             : }
     151             : 
     152             : bool
     153      387609 : DiracKernelInfo::pointsFuzzyEqual(const Point & a, const Point & b)
     154             : {
     155      387609 :   const Real dist_sq = (a - b).norm_sq();
     156      387609 :   return dist_sq < _point_equal_distance_sq;
     157             : }