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 "ProjectionAux.h"
      11             : #include "SystemBase.h"
      12             : #include "libmesh/system.h"
      13             : 
      14             : registerMooseObjectRenamed("MooseApp", SelfAux, "01/30/2024 24:00", ProjectionAux);
      15             : registerMooseObject("MooseApp", ProjectionAux);
      16             : 
      17             : InputParameters
      18       37187 : ProjectionAux::validParams()
      19             : {
      20       37187 :   InputParameters params = AuxKernel::validParams();
      21       37187 :   params.addClassDescription(
      22             :       "Returns the specified variable as an auxiliary variable with a projection of the source "
      23             :       "variable. If they are the same type, this amounts to a simple copy.");
      24       37187 :   params.addRequiredCoupledVar("v", "Variable to take the value of.");
      25             : 
      26      111561 :   params.addParam<bool>("use_block_restriction_for_source",
      27       74374 :                         false,
      28             :                         "Whether to use the auxkernel block restriction to also restrict the "
      29             :                         "locations selected for source variable values");
      30             : 
      31             :   // Technically possible to project from nodal to elemental and back
      32       37187 :   params.set<bool>("_allow_nodal_to_elemental_coupling") = true;
      33             : 
      34             :   // We need some ghosting for all elemental to nodal projections
      35       37187 :   params.addParam<unsigned short>("ghost_layers", 1, "The number of layers of elements to ghost.");
      36       37187 :   params.addRelationshipManager("ElementPointNeighborLayers",
      37             :                                 Moose::RelationshipManagerType::ALGEBRAIC,
      38           0 :                                 [](const InputParameters & obj_params, InputParameters & rm_params)
      39             :                                 {
      40       12444 :                                   rm_params.set<unsigned short>("layers") =
      41       12444 :                                       obj_params.get<unsigned short>("ghost_layers");
      42       12444 :                                   rm_params.set<bool>("use_displaced_mesh") =
      43       12444 :                                       obj_params.get<bool>("use_displaced_mesh");
      44       12444 :                                 });
      45       37187 :   return params;
      46           0 : }
      47             : 
      48        4493 : ProjectionAux::ProjectionAux(const InputParameters & parameters)
      49             :   : AuxKernel(parameters),
      50        4493 :     _v(coupledValue("v")),
      51        4493 :     _source_variable(*getFieldVar("v", 0)),
      52        4493 :     _source_sys(_c_fe_problem.getSystem(coupledName("v"))),
      53        8986 :     _use_block_restriction_for_source(getParam<bool>("use_block_restriction_for_source"))
      54             : {
      55             :   // Output some messages to user
      56        4493 :   if (_source_variable.order() > _var.order())
      57        1641 :     mooseInfo("Projection lowers order, please expect a loss of accuracy");
      58             :   // Check the dimension of the block restriction
      59        8981 :   for (const auto & sub_id : blockIDs())
      60        4488 :     if (_mesh.isLowerD(sub_id))
      61           0 :       paramError("block",
      62             :                  "ProjectionAux's block restriction must not include lower dimensional blocks");
      63        4493 : }
      64             : 
      65             : Real
      66    14210332 : ProjectionAux::computeValue()
      67             : {
      68    14210332 :   if (!isNodal() || (_source_variable.isNodal() && _source_variable.order() >= _var.order()))
      69    14037604 :     return _v[_qp];
      70             :   // projecting continuous elemental variable onto a nodal one
      71             :   // AND nodal low order -> nodal higher order
      72      215163 :   else if (isNodal() && _source_variable.getContinuity() != DISCONTINUOUS &&
      73       42435 :            _source_variable.getContinuity() != SIDE_DISCONTINUOUS)
      74             :   {
      75       84870 :     return _source_sys.point_value(
      76       84870 :         _source_variable.number(), *_current_node, elemOnNodeVariableIsDefinedOn());
      77             :   }
      78             :   // Handle discontinuous elemental variable projection into a nodal variable
      79             :   else
      80             :   {
      81             :     // Custom projection rule : use nodal values, if discontinuous the one coming from each element,
      82             :     // weighted by element volumes
      83             :     // First, find all the elements that this node is part of
      84      130293 :     auto elem_ids = _mesh.nodeToElemMap().find(_current_node->id());
      85             :     mooseAssert(elem_ids != _mesh.nodeToElemMap().end(),
      86             :                 "Should have found an element around node " + std::to_string(_current_node->id()));
      87             : 
      88             :     // Get the neighbor element centroid values & element volumes
      89      130293 :     Real sum_weighted_values = 0;
      90      130293 :     Real sum_volumes = 0;
      91      579231 :     for (auto & id : elem_ids->second)
      92             :     {
      93      448938 :       const auto & elem = _mesh.elemPtr(id);
      94      448938 :       const auto block_id = elem->subdomain_id();
      95             :       // Only use higher D elements
      96             :       // We allow full-dimensional elements in a higher dimension mesh
      97      897660 :       if (_source_variable.hasBlocks(block_id) && (!_mesh.isLowerD(block_id)) &&
      98      448722 :           (!_use_block_restriction_for_source || hasBlocks(block_id)))
      99             :       {
     100      448722 :         const auto elem_volume = elem->volume();
     101      448722 :         sum_weighted_values +=
     102      448722 :             _source_sys.point_value(_source_variable.number(), *_current_node, elem) * elem_volume;
     103      448722 :         sum_volumes += elem_volume;
     104             :       }
     105             :     }
     106      130293 :     if (sum_volumes == 0)
     107           0 :       mooseError("Did not find a valid source variable value for node: ", *_current_node);
     108      130293 :     return sum_weighted_values / sum_volumes;
     109             :   }
     110             : }
     111             : 
     112             : const Elem *
     113       42435 : ProjectionAux::elemOnNodeVariableIsDefinedOn() const
     114             : {
     115       42435 :   for (const auto & elem_id : _mesh.nodeToElemMap().find(_current_node->id())->second)
     116             :   {
     117       42435 :     const auto & elem = _mesh.elemPtr(elem_id);
     118       42435 :     const auto block_id = elem->subdomain_id();
     119       42435 :     if (_source_variable.hasBlocks(block_id) &&
     120       84870 :         (!_mesh.isLowerD(block_id) && elem->dim() == _mesh.dimension()) &&
     121       42435 :         (!_use_block_restriction_for_source || hasBlocks(block_id)))
     122       42435 :       return elem;
     123             :   }
     124           0 :   mooseError("Source variable is not defined everywhere the target variable is");
     125             : }