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 "MFEMNodalProjector.h"
      13             : #include "MFEMVectorUtils.h"
      14             : #include "MooseError.h"
      15             : 
      16             : void
      17         542 : MFEMNodalProjector::extractNodePositions(const mfem::ParFiniteElementSpace & fespace,
      18             :                                          mfem::Vector & node_positions,
      19             :                                          mfem::Ordering::Type & node_ordering)
      20             : {
      21         542 :   const int nelem = fespace.GetParMesh()->GetNE();
      22         542 :   const int dim = fespace.GetParMesh()->Dimension();
      23             : 
      24             :   // Find total number of local nodes/interpolation points
      25         542 :   int nnodes = 0;
      26       56710 :   for (const auto i : make_range(nelem))
      27       56168 :     nnodes += fespace.GetFE(i)->GetNodes().GetNPoints();
      28             : 
      29         542 :   node_positions.SetSize(nnodes * dim);
      30             : 
      31         542 :   int nodal_offset = 0;
      32       56710 :   for (const auto i : make_range(nelem))
      33             :   {
      34       56168 :     const mfem::IntegrationRule ir = fespace.GetFE(i)->GetNodes();
      35       56168 :     const int nqpt = ir.GetNPoints();
      36             : 
      37       56168 :     mfem::DenseMatrix pos;
      38       56168 :     fespace.GetElementTransformation(i)->Transform(ir, pos);
      39       56168 :     mfem::Vector row;
      40             : 
      41      198076 :     for (const auto d : make_range(dim))
      42             :     {
      43      141908 :       row.SetDataAndSize(node_positions.GetData() + nodal_offset + d * nnodes, nqpt);
      44      141908 :       pos.GetRow(d, row);
      45             :     }
      46             : 
      47       56168 :     nodal_offset += nqpt;
      48       56168 :   }
      49         542 :   node_ordering = mfem::Ordering::Type::byNODES;
      50         542 : }
      51             : 
      52             : void
      53         549 : MFEMNodalProjector::projectNodalValues(const mfem::Vector & nodal_vals,
      54             :                                        const mfem::Ordering::Type & nodal_val_ordering,
      55             :                                        mfem::ParGridFunction & gridfunction)
      56             : {
      57         549 :   mfem::ParFiniteElementSpace & fespace = *gridfunction.ParFESpace();
      58         549 :   const int nelem = fespace.GetParMesh()->GetNE();
      59         549 :   const int gf_ncomp = gridfunction.VectorDim();
      60             : 
      61         549 :   fespace.GetParMesh()->EnsureNodes();
      62             : 
      63             :   // Check FESpaces can be transferred
      64         549 :   const auto map = fespace.FEColl()->GetMapType(fespace.GetParMesh()->Dimension());
      65         549 :   const bool H1L2 = map == mfem::FiniteElement::VALUE || map == mfem::FiniteElement::INTEGRAL;
      66         549 :   const bool RTND = map == mfem::FiniteElement::H_DIV || map == mfem::FiniteElement::H_CURL;
      67         549 :   if (!H1L2 && !RTND)
      68           0 :     mooseError("FESpace type not supported yet in transfers.");
      69             : 
      70             :   // Project the interpolated values to the target FiniteElementSpace.
      71         549 :   mfem::Ordering::Type dof_ordering = H1L2 ? mfem::Ordering::byNODES : mfem::Ordering::byVDIM;
      72         549 :   mfem::Array<int> vdofs;
      73         549 :   mfem::Vector vals;
      74         549 :   int nodal_offset = 0;
      75       57417 :   for (const auto el : make_range(nelem))
      76             :   {
      77       56868 :     const int nqpt = fespace.GetFE(el)->GetNodes().GetNPoints();
      78       56868 :     mfem::Vector dof_vals(nqpt * gf_ncomp);
      79       56868 :     fespace.GetElementVDofs(el, vdofs); // Returned vdofs always indexed with ordering byNODES
      80      340640 :     for (const auto qp : make_range(nqpt))
      81      707512 :       for (const auto d : make_range(gf_ncomp))
      82      423740 :         dof_vals(Moose::MFEM::MFEMIndex(d, qp, gf_ncomp, nqpt, dof_ordering)) = nodal_vals(
      83      423740 :             Moose::MFEM::MFEMIndex(d, qp, gf_ncomp, nqpt, nodal_val_ordering) + nodal_offset);
      84             : 
      85       56868 :     if (H1L2)
      86       49380 :       gridfunction.SetSubVector(vdofs, dof_vals);
      87        7488 :     else if (RTND)
      88             :     {
      89        7488 :       vals.SetSize(vdofs.Size());
      90        7488 :       fespace.GetFE(el)->ProjectFromNodes(dof_vals, *fespace.GetElementTransformation(el), vals);
      91        7488 :       gridfunction.SetSubVector(vdofs, vals);
      92             :     }
      93             : 
      94       56868 :     nodal_offset += nqpt * gf_ncomp;
      95       56868 :   }
      96         549 :   gridfunction.SetTrueVector();
      97         549 : }
      98             : 
      99             : #endif