Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://www.mooseframework.org
       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 "AverageSectionValueSampler.h"
      11             : #include "MooseMesh.h"
      12             : #include "SystemBase.h"
      13             : #include <limits>
      14             : 
      15             : registerMooseObject("TensorMechanicsApp", AverageSectionValueSampler);
      16             : 
      17             : InputParameters
      18          18 : AverageSectionValueSampler::validParams()
      19             : {
      20          18 :   InputParameters params = GeneralVectorPostprocessor::validParams();
      21          18 :   params.addClassDescription("Compute the section's variable average in three-dimensions "
      22             :                              "given a user-defined definition of the cross section.");
      23          36 :   params.addParam<std::vector<SubdomainName>>(
      24             :       "block",
      25             :       "The list of blocks in which to search for cross sectional nodes to compute the variable "
      26             :       "average.");
      27          36 :   params.addRequiredParam<Point>("axis_direction", "Direction of the structural component's axis");
      28          36 :   params.addRequiredParam<Point>("reference_point",
      29             :                                  "Structural component reference starting point from which the "
      30             :                                  "input parameter 'lengths' applies.");
      31          36 :   params.addParam<Real>("cross_section_maximum_radius",
      32          36 :                         std::numeric_limits<double>::max(),
      33             :                         "Radial distance with respect to the body axis within which nodes are "
      34             :                         "considered to belong to this "
      35             :                         "structural component. Used to disambiguate multiple components that share "
      36             :                         "the same mesh block.");
      37          36 :   params.addRequiredParam<std::vector<VariableName>>(
      38             :       "variables", "Variables for the cross section output. These variables must be nodal.");
      39          36 :   params.addRequiredParam<std::vector<Real>>(
      40             :       "lengths",
      41             :       "Distance(s) along axis of from reference_point at which to compute average values.");
      42          36 :   params.addParam<Real>("tolerance",
      43          36 :                         1.0e-6,
      44             :                         "Maximum axial distance of nodes from the specified axial lengths to "
      45             :                         "consider them in the cross-section average");
      46          18 :   return params;
      47           0 : }
      48             : 
      49           9 : AverageSectionValueSampler::AverageSectionValueSampler(const InputParameters & parameters)
      50             :   : GeneralVectorPostprocessor(parameters),
      51           9 :     _displaced_mesh(_app.actionWarehouse().displacedMesh()),
      52           9 :     _mesh(_app.actionWarehouse().mesh()),
      53          18 :     _variables(getParam<std::vector<VariableName>>("variables")),
      54          18 :     _direction(getParam<Point>("axis_direction")),
      55          18 :     _reference_point(getParam<Point>("reference_point")),
      56          18 :     _lengths(getParam<std::vector<Real>>("lengths")),
      57          18 :     _tolerance(getParam<Real>("tolerance")),
      58           9 :     _number_of_nodes(_lengths.size()),
      59          36 :     _cross_section_maximum_radius(getParam<Real>("cross_section_maximum_radius"))
      60             : {
      61           9 :   if (_mesh->dimension() != 3)
      62           0 :     mooseError("The AverageSectionValueSampler postprocessor can only be used with three "
      63             :                "dimensional meshes.");
      64             : 
      65           9 :   if (!MooseUtils::absoluteFuzzyEqual(_direction.norm_sq(), 1.0))
      66           0 :     paramError("axis_direction", "Axis direction must be a unit vector.");
      67             : 
      68           9 :   _output_vector.resize(_variables.size());
      69          36 :   for (const auto j : make_range(_variables.size()))
      70          27 :     _output_vector[j] = &declareVector(_variables[j]);
      71             : 
      72          36 :   for (const auto j : make_range(_variables.size()))
      73             :   {
      74          27 :     const MooseVariable & variable = _sys.getFieldVariable<Real>(_tid, _variables[j]);
      75          27 :     if (!variable.isNodal())
      76           0 :       paramError(
      77             :           "variables",
      78             :           "The variables provided to this vector postprocessor must be defined at the nodes.");
      79             :   }
      80           9 : }
      81             : 
      82             : void
      83           9 : AverageSectionValueSampler::initialize()
      84             : {
      85          36 :   for (const auto j : make_range(_variables.size()))
      86             :   {
      87          27 :     _output_vector[j]->clear();
      88          27 :     _output_vector[j]->resize(_lengths.size());
      89             :   }
      90           9 :   _number_of_nodes.clear();
      91           9 :   _number_of_nodes.resize(_lengths.size());
      92           9 : }
      93             : 
      94             : void
      95           9 : AverageSectionValueSampler::finalize()
      96             : {
      97          36 :   for (const auto j : make_range(_variables.size()))
      98          27 :     _communicator.sum(*(_output_vector[j]));
      99             : 
     100           9 :   _communicator.sum(_number_of_nodes);
     101             : 
     102          27 :   for (const auto li : index_range(_lengths))
     103          18 :     if (_number_of_nodes[li] < 1)
     104           0 :       mooseError("No nodes were found in AverageSectionValueSampler postprocessor.");
     105             : 
     106          27 :   for (const auto li : index_range(_lengths))
     107          72 :     for (const auto j : make_range(_variables.size()))
     108          54 :       (*_output_vector[j])[li] /= _number_of_nodes[li];
     109           9 : }
     110             : 
     111             : void
     112           9 : AverageSectionValueSampler::execute()
     113             : {
     114          18 :   _block_ids = _displaced_mesh->getSubdomainIDs(getParam<std::vector<SubdomainName>>("block"));
     115             : 
     116           9 :   auto * active_nodes = _mesh->getActiveSemiLocalNodeRange();
     117             : 
     118             :   std::vector<std::vector<Real>> output_vector_partial(_variables.size(),
     119           9 :                                                        std::vector<Real>(_lengths.size()));
     120             : 
     121           9 :   const NumericVector<Number> & _sol = *_sys.currentSolution();
     122             : 
     123        5333 :   for (const auto & node : *active_nodes)
     124             :   {
     125        5324 :     const std::set<SubdomainID> & node_blk_ids = _displaced_mesh->getNodeBlockIds(*node);
     126             : 
     127       15972 :     for (const auto li : index_range(_lengths))
     128       21296 :       for (const auto i : _block_ids)
     129       10648 :         if (node_blk_ids.find(i) != node_blk_ids.end())
     130             :         {
     131             :           // Check if node is close enough to user-prescribed plane
     132       10648 :           if (distancePointToPlane(*node, _reference_point, _lengths[li]) > _tolerance)
     133       10064 :             continue;
     134         584 :           if ((*node).processor_id() != processor_id())
     135           8 :             continue;
     136             : 
     137             :           // Retrieve nodal variables
     138        2304 :           for (const auto j : make_range(_variables.size()))
     139             :           {
     140        1728 :             const MooseVariable & variable = _sys.getFieldVariable<Real>(_tid, _variables[j]);
     141             :             const auto var_num = variable.number();
     142        1728 :             output_vector_partial[j][li] += _sol(node->dof_number(_sys.number(), var_num, 0));
     143             :           }
     144             : 
     145         576 :           _number_of_nodes[li]++;
     146             :         }
     147             :   }
     148             : 
     149          36 :   for (const auto j : make_range(_variables.size()))
     150          27 :     *(_output_vector[j]) = output_vector_partial[j];
     151           9 : }
     152             : 
     153             : Real
     154       10648 : AverageSectionValueSampler::distancePointToPlane(const Node & node,
     155             :                                                  const Point & reference_point,
     156             :                                                  const Real length) const
     157             : {
     158             :   // Compute node location w.r.t. structural component length
     159             :   const Point relative_distance{
     160       10648 :       node(0) - reference_point(0), node(1) - reference_point(1), node(2) - reference_point(2)};
     161             : 
     162             :   const Real axial_distance = _direction * relative_distance;
     163             :   const Real in_plane_distance =
     164       10648 :       (relative_distance - relative_distance * _direction * _direction).norm();
     165             : 
     166             :   // If the in-plane distance is greater than the specified cross-section radius, the point is not
     167             :   // in this component
     168       10648 :   if (in_plane_distance > _cross_section_maximum_radius)
     169             :     return std::numeric_limits<double>::max();
     170             : 
     171        6424 :   return std::abs(axial_distance - length);
     172             : }