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             : #ifdef MOOSE_MFEM_ENABLED
      11             : 
      12             : #include "MFEMValueSamplerBase.h"
      13             : 
      14             : #include "MFEMProblem.h"
      15             : #include "MFEMVectorUtils.h"
      16             : #include "MooseError.h"
      17             : 
      18             : #include "mfem/fem/fespace.hpp"
      19             : 
      20             : namespace
      21             : {
      22             : /** Enum for values returned by gslib for point location relative to mesh
      23             :  */
      24             : enum class GSLibLocationCode : unsigned int
      25             : {
      26             :   INTERNAL = 0,
      27             :   BORDER = 1,
      28             :   NOT_FOUND = 2,
      29             : };
      30             : 
      31             : void
      32         972 : MFEMVectorToPostprocessorPoints(
      33             :     const mfem::Vector & mfem_points,
      34             :     std::vector<std::reference_wrapper<VectorPostprocessorValue>> & points,
      35             :     const unsigned int num_dims,
      36             :     const mfem::Ordering::Type ordering)
      37             : {
      38         972 :   const unsigned int num_points = mfem_points.Size() / num_dims;
      39        5116 :   for (unsigned int i_point = 0; i_point < num_points; i_point++)
      40             :   {
      41       12930 :     for (unsigned int i_dim = 0; i_dim < num_dims; i_dim++)
      42             :     {
      43        8786 :       const size_t idx = Moose::MFEM::MFEMIndex(i_dim, i_point, num_dims, num_points, ordering);
      44             : 
      45        8786 :       points[i_dim].get()[i_point] = mfem_points(idx);
      46             :     }
      47             :   }
      48         972 : }
      49             : 
      50             : mfem::FindPointsGSLIB::AvgType
      51         366 : getAvgType(const MooseEnum & avg_type)
      52             : {
      53         366 :   if (avg_type == "NONE")
      54          13 :     return mfem::FindPointsGSLIB::AvgType::NONE;
      55         353 :   else if (avg_type == "ARITHMETIC")
      56         340 :     return mfem::FindPointsGSLIB::AvgType::ARITHMETIC;
      57          13 :   else if (avg_type == "HARMONIC")
      58          13 :     return mfem::FindPointsGSLIB::AvgType::HARMONIC;
      59             :   else
      60           0 :     mooseError(
      61             :         "Unknown average type: ",
      62             :         avg_type,
      63             :         ", this is an internal MOOSE error, new enum variants must be handled by this function.");
      64             : }
      65             : }
      66             : 
      67             : InputParameters
      68        4978 : MFEMValueSamplerBase::validParams()
      69             : {
      70        4978 :   InputParameters params = MFEMVectorPostprocessor::validParams();
      71             : 
      72       19912 :   MFEMExecutedObject::addRequiredDependencyParam<VariableName>(
      73             :       params, "variable", "The names of the variables that this VectorPostprocessor operates on");
      74       19912 :   MooseEnum ordering("NODES VDIM", "VDIM", false);
      75       19912 :   params.addParam<MooseEnum>(
      76             :       "point_ordering", ordering, "Ordering style to use for point vector DoFs.");
      77       19912 :   MooseEnum avg_type("NONE ARITHMETIC HARMONIC", "ARITHMETIC", false);
      78       19912 :   params.addParam<MooseEnum>("average_type",
      79             :                              avg_type,
      80             :                              "Average type used when sampling L2 functions at element boundaries.");
      81        9956 :   params.addParam<double>("mesh_boundary_tolerance",
      82        9956 :                           1e-8,
      83             :                           "Distance from point to mesh boundary below which the point is "
      84             :                           "considered to be on the boundary rather than outside the mesh.");
      85             : 
      86        9956 :   return params;
      87        4978 : }
      88             : 
      89         366 : MFEMValueSamplerBase::MFEMValueSamplerBase(const InputParameters & parameters,
      90         366 :                                            const std::vector<Point> & points)
      91             :   : MFEMVectorPostprocessor(parameters),
      92         366 :     _var_name(getParam<VariableName>("variable")),
      93         366 :     _var(*getMFEMProblem().getGridFunction(_var_name)),
      94         366 :     _mesh(const_cast<mfem::ParMesh &>(getMFEMProblem().getMFEMVariableMesh(_var_name))),
      95         366 :     _finder(this->comm().get()),
      96         732 :     _points_ordering(getParam<MooseEnum>("point_ordering") == "NODES" ? mfem::Ordering::byNODES
      97             :                                                                       : mfem::Ordering::byVDIM),
      98         366 :     _points(
      99         366 :         Moose::MFEM::libMeshPointsToMFEMVector(points, _mesh.SpaceDimension(), _points_ordering)),
     100         732 :     _interp_vals(points.size())
     101             : {
     102         366 :   if (getMFEMProblem().mesh().shouldDisplace())
     103           0 :     mooseError("MFEMValueSamplerBase does not yet support problems with displacement.");
     104             : 
     105         732 :   _finder.SetL2AvgType(getAvgType(getParam<MooseEnum>("average_type")));
     106         732 :   _finder.SetDistanceToleranceForPointsFoundOnBoundary(getParam<double>("mesh_boundary_tolerance"));
     107             : 
     108             :   // set up points vector
     109         366 :   _mesh.EnsureNodes();
     110         366 :   _finder.Setup(_mesh);
     111         366 :   _finder.FindPoints(_points, _points_ordering);
     112             : 
     113             :   bool fe_boundary_discontinuous =
     114         366 :       _var.FESpace()->FEColl()->GetContType() == mfem::FiniteElementCollection::DISCONTINUOUS;
     115             : 
     116             :   // check all points were found
     117         366 :   mfem::Array<unsigned int> point_codes = _finder.GetCode();
     118        3850 :   for (size_t i = 0; i < points.size(); i++)
     119             :   {
     120        3490 :     switch (GSLibLocationCode(point_codes[i]))
     121             :     {
     122        2711 :       case GSLibLocationCode::INTERNAL:
     123        2711 :         break;
     124         775 :       case GSLibLocationCode::BORDER:
     125         775 :         if (fe_boundary_discontinuous)
     126         233 :           mooseWarning("MFEMValueSamplerBase found a point on an element boundary but "
     127             :                        "the FE space is discontinuous at boundaries: ",
     128         233 :                        points[i],
     129             :                        ".");
     130         773 :         break;
     131           4 :       default:
     132           4 :         mooseError("MFEMValueSamplerBase could not find point at ", points[i], ".");
     133             :         break;
     134             :     }
     135             :   }
     136             : 
     137             :   // declare points vectors for outputting
     138         360 :   const auto mesh_dim = _mesh.SpaceDimension();
     139        1196 :   for (int i = 0; i < mesh_dim; i++)
     140             :   {
     141             :     std::reference_wrapper<VectorPostprocessorValue> declared_dim =
     142         836 :         this->declareVector("x_" + std::to_string(i));
     143         836 :     declared_dim.get().resize(points.size());
     144         836 :     _declared_points.push_back(declared_dim);
     145             :   }
     146             : 
     147             :   // declare value vectors for outputting
     148         360 :   const auto val_dim = _var.VectorDim();
     149         784 :   for (int i = 0; i < val_dim; i++)
     150             :   {
     151             :     std::reference_wrapper<VectorPostprocessorValue> declared_dim =
     152         424 :         this->declareVector(_var_name + "_" + std::to_string(i));
     153         424 :     declared_dim.get().resize(points.size());
     154         424 :     _declared_vals.push_back(declared_dim);
     155             :   }
     156         360 : }
     157             : 
     158             : void
     159         972 : MFEMValueSamplerBase::execute()
     160             : {
     161         972 :   _finder.Interpolate(_var, _interp_vals);
     162         972 : }
     163             : 
     164             : void
     165         972 : MFEMValueSamplerBase::finalize()
     166             : {
     167         972 :   _interp_vals.HostReadWrite();
     168         972 :   _points.HostReadWrite();
     169             : 
     170         972 :   const auto mesh_dim = _mesh.SpaceDimension();
     171         972 :   MFEMVectorToPostprocessorPoints(_points, _declared_points, mesh_dim, _points_ordering);
     172         972 :   const auto val_dims = _var.VectorDim();
     173         972 :   const auto num_points = _declared_points[0].get().size();
     174         972 :   const auto val_fespace_ordering = _var.FESpace()->GetOrdering();
     175        2008 :   for (int i_dim = 0; i_dim < val_dims; i_dim++)
     176             :   {
     177        5638 :     for (size_t i_point = 0; i_point < num_points; i_point++)
     178             :     {
     179             :       const auto mfem_idx =
     180        4602 :           Moose::MFEM::MFEMIndex(i_dim, i_point, val_dims, num_points, val_fespace_ordering);
     181        4602 :       _declared_vals[i_dim].get()[i_point] = _interp_vals[mfem_idx];
     182             :     }
     183             :   }
     184         972 : }
     185             : 
     186             : #endif // MOOSE_MFEM_ENABLED