DiracKernelInfo.C

Go to the documentation of this file.

//* This file is part of the MOOSE framework
//* https://mooseframework.inl.gov
//*
//* All rights reserved, see COPYRIGHT for full restrictions
//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
//*
//* Licensed under LGPL 2.1, please see LICENSE for details
//* https://www.gnu.org/licenses/lgpl-2.1.html

#include "DiracKernelInfo.h"
#include "MooseMesh.h"
#include "MooseEnum.h"
#include "DiracKernelBase.h"
#include "MooseBase.h"

// LibMesh
#include "libmesh/point_locator_base.h"
#include "libmesh/elem.h"
#include "libmesh/enum_point_locator_type.h"
#include "libmesh/point.h"

using namespace libMesh;

DiracKernelInfo::DiracKernelInfo()
  : _point_locator(), _point_equal_distance_sq(libMesh::TOLERANCE * libMesh::TOLERANCE)
{
}

DiracKernelInfo::~DiracKernelInfo() {}

void
DiracKernelInfo::addPoint(const Elem * elem, const Point & p, const Real & value)
{
  _elements.insert(elem);

  std::pair<std::vector<Point>, std::vector<Real>> & multi_point_list = _points[elem];

  const unsigned int npoint = multi_point_list.first.size();
  mooseAssert(npoint == multi_point_list.second.size(),
              "Different sizes for location and point value data");

  for (unsigned int i = 0; i < npoint; ++i)
    if (pointsFuzzyEqual(multi_point_list.first[i], p))
    {
      // A point at the same (within a tolerance) location as p exists, accumulate its value.
      multi_point_list.second[i] += value;
      return;
    }

  // No prior point found at this location, add it with its initial value.
  multi_point_list.first.push_back(p);
  multi_point_list.second.push_back(value);
}

void
DiracKernelInfo::clearPoints()
{
  _elements.clear();
  _points.clear();
}

bool
DiracKernelInfo::hasPoint(const Elem * elem, const Point & p)
{
  std::vector<Point> & point_list = _points[elem].first;

  for (const auto & pt : point_list)
    if (pointsFuzzyEqual(pt, p))
      return true;

  // If we haven't found it, we don't have it.
  return false;
}

void
DiracKernelInfo::updatePointLocator(const MooseMesh & mesh)
{
  // Note: we could update the PointLocator *every* time we call this
  // function, but that may introduce an unacceptable overhead in
  // problems which don't need a PointLocator at all.  This issue will
  // most likely become a moot point when we eventually add a shared
  // "CachingPointLocator" to MOOSE.
  // _point_locator = PointLocatorBase::build(TREE_LOCAL_ELEMENTS, mesh);

  // Construct the PointLocator object if *any* processors have Dirac
  // points.  Note: building a PointLocator object is a parallel_only()
  // function, so this is an all-or-nothing thing.
  unsigned pl_needs_rebuild = _elements.size();
  mesh.comm().max(pl_needs_rebuild);

  if (pl_needs_rebuild)
  {
    // PointLocatorBase::build() is a parallel_only function!  So we
    // can't skip building it just becuase our local _elements is
    // empty, it might be non-empty on some other processor!
    _point_locator = PointLocatorBase::build(TREE_LOCAL_ELEMENTS, mesh);

    // We may be querying for points which are not in the semilocal
    // part of a distributed mesh.
    _point_locator->enable_out_of_mesh_mode();
  }
  else
  {
    // There are no elements with Dirac points, but we have been
    // requested to update the PointLocator so we have to assume the
    // old one is invalid.  Therefore we reset it to NULL... however
    // adding this line causes the code to hang because it triggers
    // the PointLocator to be rebuilt in a non-parallel-only segment
    // of the code later... so it's commented out for now even though
    // it's probably the right behavior.
    // _point_locator.reset(NULL);
  }
}

const Elem *
DiracKernelInfo::findPoint(const Point & p,
                           const MooseMesh & mesh,
                           const std::set<SubdomainID> & blocks,
                           const PointNotFoundBehavior point_not_found_behavior,
                           const MooseBase & consumer)
{
  // If the PointLocator has never been created, do so now.  NOTE - WE
  // CAN'T DO THIS if findPoint() is only called on some processors,
  // PointLocatorBase::build() is a 'parallel_only' method!
  if (_point_locator.get() == NULL)
  {
    _point_locator = PointLocatorBase::build(TREE_LOCAL_ELEMENTS, mesh);
    _point_locator->enable_out_of_mesh_mode();
  }

  // Check that the PointLocator is ready to start locating points.
  // So far I do not have any tests that trip this...
  if (_point_locator->initialized() == false)
    consumer.mooseError("PointLocator is not initialized!");

  // Note: The PointLocator object returns NULL when the Point is not
  // found within the Mesh.  This is not considered to be an error as
  // far as the DiracKernels are concerned: sometimes the Mesh moves
  // out from the Dirac point entirely and in that case the Point just
  // gets "deactivated".
  const Elem * elem = (*_point_locator)(p, &blocks);

  // The processors may not agree on which Elem the point is in.  This
  // can happen if a Dirac point lies on the processor boundary, and
  // two or more neighboring processors think the point is in the Elem
  // on *their* side.
  dof_id_type elem_id = elem ? elem->id() : DofObject::invalid_id;

  // We are going to let the element with the smallest ID "win", all other
  // procs will return NULL.
  dof_id_type min_elem_id = elem_id;
  mesh.comm().min(min_elem_id);

  if (min_elem_id == DofObject::invalid_id)
  {
    std::stringstream msg;
    msg << "Point " << p << " not found in block(s) " << Moose::stringify(blocks, ", ") << ".\n";
    switch (point_not_found_behavior)
    {
      case PointNotFoundBehavior::ERROR:
        consumer.mooseError(msg.str());
        break;
      case PointNotFoundBehavior::WARNING:
        mooseDoOnce(consumer.mooseWarning(msg.str() + "This message will not be repeated."));
        break;
      case PointNotFoundBehavior::IGNORE:
        break;
      default:
        consumer.mooseError("Internal enum error.");
    }
  }

  // But we notably need the processor which owns elem_id to return it!
  if (min_elem_id != DofObject::invalid_id)
    if (const auto min_elem = mesh.queryElemPtr(min_elem_id);
        min_elem && min_elem->processor_id() == mesh.processor_id())
      return min_elem;
  return nullptr;
}

bool
DiracKernelInfo::pointsFuzzyEqual(const Point & a, const Point & b)
{
  const Real dist_sq = (a - b).norm_sq();
  return dist_sq < _point_equal_distance_sq;
}