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       36557 : ProjectionAux::validParams()
      19             : {
      20       36557 :   InputParameters params = AuxKernel::validParams();
      21       36557 :   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       36557 :   params.addRequiredCoupledVar("v", "Variable to take the value of.");
      25             : 
      26      109671 :   params.addParam<bool>("use_block_restriction_for_source",
      27       73114 :                         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       36557 :   params.set<bool>("_allow_nodal_to_elemental_coupling") = true;
      33             : 
      34             :   // We need some ghosting for all elemental to nodal projections
      35       36557 :   params.addParam<unsigned short>("ghost_layers", 1, "The number of layers of elements to ghost.");
      36       36557 :   params.addRelationshipManager("ElementPointNeighborLayers",
      37             :                                 Moose::RelationshipManagerType::ALGEBRAIC,
      38           0 :                                 [](const InputParameters & obj_params, InputParameters & rm_params)
      39             :                                 {
      40       11499 :                                   rm_params.set<unsigned short>("layers") =
      41       11499 :                                       obj_params.get<unsigned short>("ghost_layers");
      42       11499 :                                   rm_params.set<bool>("use_displaced_mesh") =
      43       11499 :                                       obj_params.get<bool>("use_displaced_mesh");
      44       11499 :                                 });
      45       36557 :   return params;
      46           0 : }
      47             : 
      48        4178 : ProjectionAux::ProjectionAux(const InputParameters & parameters)
      49             :   : AuxKernel(parameters),
      50        4178 :     _v(coupledValue("v")),
      51        4178 :     _source_variable(*getFieldVar("v", 0)),
      52        4178 :     _source_sys(_c_fe_problem.getSystem(coupledName("v"))),
      53        8356 :     _use_block_restriction_for_source(getParam<bool>("use_block_restriction_for_source"))
      54             : {
      55             :   // Output some messages to user
      56        4178 :   if (_source_variable.order() > _var.order())
      57        1524 :     mooseInfo("Projection lowers order, please expect a loss of accuracy");
      58             :   // Check the dimension of the block restriction
      59        8351 :   for (const auto & sub_id : blockIDs())
      60        4173 :     if (_mesh.isLowerD(sub_id))
      61           0 :       paramError("block",
      62             :                  "ProjectionAux's block restriction must not include lower dimensional blocks");
      63        4178 : }
      64             : 
      65             : Real
      66    12646224 : ProjectionAux::computeValue()
      67             : {
      68    12646224 :   if (!isNodal() || (_source_variable.isNodal() && _source_variable.order() >= _var.order()))
      69    12492688 :     return _v[_qp];
      70             :   // projecting continuous elemental variable onto a nodal one
      71             :   // AND nodal low order -> nodal higher order
      72      191256 :   else if (isNodal() && _source_variable.getContinuity() != DISCONTINUOUS &&
      73       37720 :            _source_variable.getContinuity() != SIDE_DISCONTINUOUS)
      74             :   {
      75       75440 :     return _source_sys.point_value(
      76       75440 :         _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      115816 :     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      115816 :     Real sum_weighted_values = 0;
      90      115816 :     Real sum_volumes = 0;
      91      514872 :     for (auto & id : elem_ids->second)
      92             :     {
      93      399056 :       const auto & elem = _mesh.elemPtr(id);
      94      399056 :       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      797920 :       if (_source_variable.hasBlocks(block_id) && (!_mesh.isLowerD(block_id)) &&
      98      398864 :           (!_use_block_restriction_for_source || hasBlocks(block_id)))
      99             :       {
     100      398864 :         const auto elem_volume = elem->volume();
     101      398864 :         sum_weighted_values +=
     102      398864 :             _source_sys.point_value(_source_variable.number(), *_current_node, elem) * elem_volume;
     103      398864 :         sum_volumes += elem_volume;
     104             :       }
     105             :     }
     106      115816 :     if (sum_volumes == 0)
     107           0 :       mooseError("Did not find a valid source variable value for node: ", *_current_node);
     108      115816 :     return sum_weighted_values / sum_volumes;
     109             :   }
     110             : }
     111             : 
     112             : const Elem *
     113       37720 : ProjectionAux::elemOnNodeVariableIsDefinedOn() const
     114             : {
     115       37720 :   for (const auto & elem_id : _mesh.nodeToElemMap().find(_current_node->id())->second)
     116             :   {
     117       37720 :     const auto & elem = _mesh.elemPtr(elem_id);
     118       37720 :     const auto block_id = elem->subdomain_id();
     119       37720 :     if (_source_variable.hasBlocks(block_id) &&
     120       75440 :         (!_mesh.isLowerD(block_id) && elem->dim() == _mesh.dimension()) &&
     121       37720 :         (!_use_block_restriction_for_source || hasBlocks(block_id)))
     122       37720 :       return elem;
     123             :   }
     124           0 :   mooseError("Source variable is not defined everywhere the target variable is");
     125             : }