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             : #pragma once
      11             : 
      12             : // MOOSE includes
      13             : #include "GeneralVectorPostprocessor.h"
      14             : #include "SamplerBase.h"
      15             : #include "BlockRestrictable.h"
      16             : #include "LineSegment.h"
      17             : #include "RayTracing.h" // Moose::elementsIntersectedByLine()
      18             : #include "Assembly.h"   // Assembly::qRule()
      19             : #include "MooseMesh.h"  // MooseMesh::getMesh()
      20             : #include "SwapBackSentinel.h"
      21             : #include "FEProblem.h"
      22             : 
      23             : #include "libmesh/quadrature.h" // _qrule->n_points()
      24             : 
      25             : // Forward Declarations
      26             : class MooseMesh;
      27             : 
      28             : /**
      29             :  * This is a base class for sampling material properties for the
      30             :  * integration points in all elements that are intersected by a
      31             :  * user-defined line.  The positions of those points are output in x,
      32             :  * y, z coordinates, as well as in terms of the projected positions of
      33             :  * those points along the line.  Derived classes can be created to
      34             :  * sample arbitrary types of material properties.
      35             :  */
      36             : template <typename T>
      37             : class LineMaterialSamplerBase : public GeneralVectorPostprocessor,
      38             :                                 public SamplerBase,
      39             :                                 public BlockRestrictable
      40             : {
      41             : public:
      42             :   static InputParameters validParams();
      43             : 
      44             :   /**
      45             :    * Class constructor
      46             :    * Sets up variables for output based on the properties to be output
      47             :    * @param parameters The input parameters
      48             :    */
      49             :   LineMaterialSamplerBase(const InputParameters & parameters);
      50             : 
      51             :   /**
      52             :    * Initialize
      53             :    * Calls through to base class's initialize()
      54             :    */
      55             :   virtual void initialize() override;
      56             : 
      57             :   /**
      58             :    * Finds all elements along the user-defined line, loops through them, and samples their
      59             :    * material properties.
      60             :    */
      61             :   virtual void execute() override;
      62             : 
      63             :   /**
      64             :    * Finalize
      65             :    * Calls through to base class's finalize()
      66             :    */
      67             :   virtual void finalize() override;
      68             : 
      69             :   /**
      70             :    * Reduce the material property to a scalar for output
      71             :    * @param property The material property
      72             :    * @param curr_point The point corresponding to this material property
      73             :    * @return A scalar value from this material property to be output
      74             :    */
      75             :   virtual Real getScalarFromProperty(const T & property, const Point & curr_point) = 0;
      76             : 
      77             : protected:
      78             :   /// The beginning of the line
      79             :   Point _start;
      80             : 
      81             :   /// The end of the line
      82             :   Point _end;
      83             : 
      84             :   /// The material properties to be output
      85             :   std::vector<const MaterialProperty<T> *> _material_properties;
      86             : 
      87             :   /// The mesh
      88             :   MooseMesh & _mesh;
      89             : 
      90             :   /// The quadrature rule
      91             :   const QBase * const & _qrule;
      92             : 
      93             :   /// The quadrature points
      94             :   const MooseArray<Point> & _q_point;
      95             : };
      96             : 
      97             : template <typename T>
      98             : InputParameters
      99       14309 : LineMaterialSamplerBase<T>::validParams()
     100             : {
     101       14309 :   InputParameters params = GeneralVectorPostprocessor::validParams();
     102       14309 :   params += SamplerBase::validParams();
     103       14309 :   params += BlockRestrictable::validParams();
     104       14309 :   params.addRequiredParam<Point>("start", "The beginning of the line");
     105       14309 :   params.addRequiredParam<Point>("end", "The end of the line");
     106       14309 :   params.addRequiredParam<std::vector<std::string>>(
     107             :       "property", "Name of the material property to be output along a line");
     108             : 
     109       14309 :   return params;
     110           0 : }
     111             : 
     112             : template <typename T>
     113          22 : LineMaterialSamplerBase<T>::LineMaterialSamplerBase(const InputParameters & parameters)
     114             :   : GeneralVectorPostprocessor(parameters),
     115             :     SamplerBase(parameters, this, _communicator),
     116             :     BlockRestrictable(this),
     117          22 :     _start(getParam<Point>("start")),
     118          22 :     _end(getParam<Point>("end")),
     119          22 :     _mesh(_subproblem.mesh()),
     120          22 :     _qrule(_subproblem.assembly(_tid, 0).qRule()),
     121          66 :     _q_point(_subproblem.assembly(_tid, 0).qPoints())
     122             : {
     123             :   // See https://github.com/idaholab/moose/issues/21865
     124          22 :   _mesh.errorIfDistributedMesh("LineMaterialSamplerBase");
     125             : 
     126          22 :   std::vector<std::string> material_property_names = getParam<std::vector<std::string>>("property");
     127          44 :   for (unsigned int i = 0; i < material_property_names.size(); ++i)
     128             :   {
     129          22 :     if (!hasMaterialProperty<T>(material_property_names[i]))
     130           0 :       mooseError("In LineMaterialSamplerBase material property: " + material_property_names[i] +
     131             :                  " does not exist.");
     132          22 :     _material_properties.push_back(&getMaterialProperty<T>(material_property_names[i]));
     133             :   }
     134             : 
     135          22 :   SamplerBase::setupVariables(material_property_names);
     136          22 : }
     137             : 
     138             : template <typename T>
     139             : void
     140          20 : LineMaterialSamplerBase<T>::initialize()
     141             : {
     142          20 :   SamplerBase::initialize();
     143          20 : }
     144             : 
     145             : template <typename T>
     146             : void
     147          20 : LineMaterialSamplerBase<T>::execute()
     148             : {
     149          20 :   std::vector<Elem *> intersected_elems;
     150          20 :   std::vector<LineSegment> segments;
     151             : 
     152          20 :   std::unique_ptr<libMesh::PointLocatorBase> pl = _mesh.getPointLocator();
     153          20 :   Moose::elementsIntersectedByLine(_start, _end, _mesh, *pl, intersected_elems, segments);
     154             : 
     155          20 :   const RealVectorValue line_vec = _end - _start;
     156          20 :   const Real line_length(line_vec.norm());
     157          20 :   const RealVectorValue line_unit_vec = line_vec / line_length;
     158          20 :   std::vector<Real> values(_material_properties.size());
     159             : 
     160          20 :   std::unordered_set<unsigned int> needed_mat_props;
     161          20 :   const auto & mp_deps = getMatPropDependencies();
     162          20 :   needed_mat_props.insert(mp_deps.begin(), mp_deps.end());
     163          20 :   _fe_problem.setActiveMaterialProperties(needed_mat_props, _tid);
     164             : 
     165         274 :   for (const auto & elem : intersected_elems)
     166             :   {
     167         146 :     if (elem->processor_id() != processor_id())
     168          38 :       continue;
     169             : 
     170         108 :     if (!hasBlocks(elem->subdomain_id()))
     171           0 :       continue;
     172             : 
     173         108 :     _subproblem.setCurrentSubdomainID(elem, _tid);
     174         108 :     _subproblem.prepare(elem, _tid);
     175         108 :     _subproblem.reinitElem(elem, _tid);
     176             : 
     177             :     // Set up Sentinel class so that, even if reinitMaterials() throws, we
     178             :     // still remember to swap back during stack unwinding.
     179         108 :     SwapBackSentinel sentinel(_fe_problem, &FEProblem::swapBackMaterials, _tid);
     180         108 :     _fe_problem.reinitMaterials(elem->subdomain_id(), _tid);
     181             : 
     182         380 :     for (unsigned int qp = 0; qp < _qrule->n_points(); ++qp)
     183             :     {
     184         272 :       const RealVectorValue qp_pos(_q_point[qp]);
     185             : 
     186         272 :       const RealVectorValue start_to_qp(qp_pos - _start);
     187         272 :       const Real qp_proj_dist_along_line = start_to_qp * line_unit_vec;
     188             : 
     189         272 :       if (qp_proj_dist_along_line < 0 || qp_proj_dist_along_line > line_length)
     190          28 :         continue;
     191             : 
     192         488 :       for (unsigned int j = 0; j < _material_properties.size(); ++j)
     193         244 :         values[j] = getScalarFromProperty((*_material_properties[j])[qp], _q_point[qp]);
     194             : 
     195         244 :       addSample(_q_point[qp], qp_proj_dist_along_line, values);
     196             :     }
     197             :   }
     198          20 :   _fe_problem.clearActiveMaterialProperties(_tid);
     199          20 : }
     200             : 
     201             : template <typename T>
     202             : void
     203          20 : LineMaterialSamplerBase<T>::finalize()
     204             : {
     205          20 :   SamplerBase::finalize();
     206          20 : }