Remove white spaces

fem/fespace.cpp

@@ -4051,8 +4051,11 @@ void FiniteElementSpace::GetTransferOperator(
    const FiniteElementSpace &coarse_fes, OperatorHandle &T) const
 {
    // Assumptions: see the declaration of the method.
-
-   if (T.Type() == Operator::MFEM_SPARSEMAT)
+   if (coarse_fes.FEColl()->Name() != fec->Name())
+   {
+      T.Reset(new GenericTransferOperator(coarse_fes, *this));
+   }
+   else if (T.Type() == Operator::MFEM_SPARSEMAT)
    {
       if (!IsVariableOrder())
       {

fem/transfer.cpp

@@ -138,6 +138,29 @@ const Operator &GridTransfer::MakeTrueOperator(
    return *t_oper.Ptr();
 }
 
+const Operator &GenericGridTransfer::ForwardOperator()
+{
+   if (F.Ptr())
+   {
+      return *F.Ptr();
+   }
+
+   F.Reset(new GenericTransferOperator(ran_fes, dom_fes));
+   return *F.Ptr();
+}
+
+const Operator &GenericGridTransfer::BackwardOperator()
+{
+   if (B.Ptr())
+   {
+      return *B.Ptr();
+   }
+
+   B.Reset(new GenericTransferOperator(dom_fes, ran_fes));
+   return *B.Ptr();
+}
+
+
 
 InterpolationGridTransfer::~InterpolationGridTransfer()
 {
@@ -2032,6 +2055,75 @@ bool L2ProjectionGridTransfer::SupportsBackwardsOperator() const
    return ran_fes.GetTrueVSize() >= dom_fes.GetTrueVSize();
 }
 
+GenericTransferOperator::GenericTransferOperator(FiniteElementSpace& dom_fes,
+                                                 FiniteElementSpace& ran_fes)
+   : Operator(ran_fes.GetVSize(), dom_fes.GetVSize()),
+     dom_gf(new GridFunction(&dom_fes)),
+     ran_gf(new GridFunction(&ran_fes))
+{
+   MFEM_VERIFY(dom_fes.GetVectorDim() == ran_fes.GetVectorDim(),
+               "GenericTransferOperator: domainn and range VectorDim do not match");
+
+   if (dom_fes.GetVectorDim() == 1)
+   {
+      dom_cf = new GridFunctionCoefficient(dom_gf);
+      ran_cf = new GridFunctionCoefficient(ran_gf);
+   }
+   else
+   {
+      dom_vcf = new VectorGridFunctionCoefficient(dom_gf);
+      ran_vcf = new VectorGridFunctionCoefficient(ran_gf);
+   }
+}
+
+GenericTransferOperator::~GenericTransferOperator()
+{
+   delete dom_gf;
+   delete ran_gf;
+   if (dom_cf) { delete dom_cf; }
+   if (ran_cf) { delete ran_cf; }
+   if (dom_vcf) { delete dom_vcf; }
+   if (ran_vcf) { delete ran_vcf; }
+}
+
+void GenericTransferOperator::Mult(const Vector& x, Vector& y) const
+{
+   dom_gf->SetFromTrueDofs(x);
+   if (dom_cf)
+   {
+      ran_gf->ProjectCoefficient(*dom_cf);
+   }
+   else if (dom_vcf)
+   {
+      ran_gf->ProjectCoefficient(*dom_vcf);
+   }
+   else
+   {
+      mfem_error("GenericTransferOperator::Mult\n"
+                 " coefficient not defined");
+   }
+   ran_gf->GetTrueDofs(y);
+}
+
+void GenericTransferOperator::MultTranspose(const Vector& x, Vector& y) const
+{
+   ran_gf->SetFromTrueDofs(x);
+   if (ran_cf)
+   {
+      dom_gf->ProjectCoefficient(*ran_cf);
+   }
+   else if (ran_vcf)
+   {
+      dom_gf->ProjectCoefficient(*ran_vcf);
+   }
+   else
+   {
+      mfem_error("GenericTransferOperator::MultTranspose\n"
+                 " coefficient not defined");
+   }
+   dom_gf->GetTrueDofs(y);
+}
+
 
 TransferOperator::TransferOperator(const FiniteElementSpace& lFESpace_,
                                    const FiniteElementSpace& hFESpace_)

fem/transfer.hpp

@@ -116,6 +116,29 @@ public:
 };
 
 
+/** @brief Transfer data between two FiniteElementSpace% based on embedded
+    refined meshes but using arbitrary FiniteElementCollections. */
+class GenericGridTransfer : public GridTransfer
+{
+protected:
+   OperatorHandle F; ///< Forward, coarse-to-fine, operator
+   OperatorHandle B; ///< Backward, fine-to-coarse, operator
+
+public:
+   GenericGridTransfer(FiniteElementSpace &dom_fes,
+                       FiniteElementSpace &ran_fes)
+      : GridTransfer(dom_fes, ran_fes)
+   { }
+
+   virtual ~GenericGridTransfer() {}
+
+   const Operator &ForwardOperator() override;
+
+   const Operator &BackwardOperator() override;
+};
+
+
+
 /** @brief Transfer data between a coarse mesh and an embedded refined mesh
     using interpolation. */
 /** The forward, coarse-to-fine, transfer uses nodal interpolation. The
@@ -531,6 +554,47 @@ private:
    void BuildF();
 };
 
+
+/// Matrix-free transfer operator between finite element spaces
+class GenericTransferOperator : public Operator
+{
+private:
+   GridFunction*  dom_gf = nullptr;
+   GridFunction*  ran_gf = nullptr;
+
+   Coefficient*  dom_cf = nullptr;
+   Coefficient*  ran_cf = nullptr;
+
+   VectorCoefficient*  dom_vcf = nullptr;
+   VectorCoefficient*  ran_vcf = nullptr;
+public:
+   /// Constructs a transfer operator from \p dom_fes to \p ran_fes.
+   /** No matrices are assembled, only the action to a vector is being computed.
+       The assumption is that grid%s are related. Meaning they are either equal
+       or refined. This class leverages GridFunctionCoefficient or
+       GridFunctionCoefficient. Both use RefinedToCoarse to establish a
+       connection between the meshes.*/
+   GenericTransferOperator(FiniteElementSpace& dom_fes,
+                           FiniteElementSpace& ran_fes);
+
+   GenericTransferOperator(const FiniteElementSpace& dom_fes,
+                           const FiniteElementSpace& ran_fes)
+      :GenericTransferOperator(const_cast<FiniteElementSpace&>(dom_fes),
+                               const_cast<FiniteElementSpace&>(ran_fes)) {};
+
+   /// Destructor
+   virtual ~GenericTransferOperator();
+
+   /// @brief Interpolation or prolongation of a vector \p x corresponding to
+   /// the coarse space to the vector \p y corresponding to the fine space.
+   void Mult(const Vector& x, Vector& y) const override;
+
+   /// Restriction by applying the transpose of the Mult method.
+   /** The vector \p x corresponding to the fine space is restricted to the
+       vector \p y corresponding to the coarse space. */
+   void MultTranspose(const Vector& x, Vector& y) const override;
+};
+
 /// Matrix-free transfer operator between finite element spaces
 class TransferOperator : public Operator
 {

tests/unit/fem/test_transfer.cpp

@@ -460,6 +460,141 @@ TEST_CASE("Restriction Transpose Operator")
    REQUIRE(y3.Normlinf() == MFEM_Approx(0.0));
 }
 
+void TestGenericTransfer(Mesh *mesh, int order, int lor)
+{
+   // Define NURBS gridfunction
+   NURBSFECollection nurbs_coll(order);
+   FiniteElementSpace nurbs_fes(mesh, new NURBSExtension(mesh->NURBSext, order),
+                                &nurbs_coll);
+   GridFunction nurbs_gf(&nurbs_fes);
+
+   // Define H1 gridfunction on refined mesh
+   Mesh h1_mesh = Mesh::MakeRefined(*mesh, lor, BasisType::GaussLobatto);
+   H1_FECollection h1_coll(order, mesh->Dimension());
+   FiniteElementSpace h1_fes(mesh, &h1_coll);
+   GridFunction h1_gf(&h1_fes);
+
+   // Get Transfer Operator
+   OperatorHandle nurbs_to_h1;
+   h1_fes.GetTransferOperator(nurbs_fes, nurbs_to_h1);
+
+   // Project coefficient of NURBS gridfunction
+   CartesianXCoefficient xcf;
+   CartesianYCoefficient ycf;
+   ProductCoefficient  cf(xcf, ycf);
+   nurbs_gf.ProjectCoefficient(cf);
+
+   // Transfer NURBS gridfunction to H1 gridfunction
+   nurbs_to_h1.Ptr()->Mult(nurbs_gf, h1_gf);
+   REQUIRE(h1_gf.ComputeL2Error(cf) == MFEM_Approx(0.0));
+
+   // Transfer H1 gridfunction to NURBS gridfunction
+   nurbs_to_h1.Ptr()->MultTranspose(h1_gf, nurbs_gf);
+   REQUIRE(nurbs_gf.ComputeL2Error(cf) == MFEM_Approx(0.0));
+}
+
+TEST_CASE("Generic Transfer Operator", "[Dimension][Order][LOR]")
+{
+   Mesh mesh2d("../../data/square-nurbs.mesh", 1, 1);
+   mesh2d.UniformRefinement();
+
+   Mesh mesh3d("../../data/cube-nurbs.mesh", 1, 1);
+   mesh3d.UniformRefinement();
+
+   dimension = GENERATE(2, 3);
+   int order = GENERATE(2, 3, 4);
+   int lor = GENERATE(2, 3, 4, 5);
+
+   if (dimension == 2)
+   {
+      TestGenericTransfer(&mesh2d, order, lor);
+   }
+   else
+   {
+      TestGenericTransfer(&mesh3d, order, lor);
+   }
+}
+
+/* This test requires PR #4326
+TEST_CASE("Generic Transfer Operator -- Vector", "[Dimension][Order][LOR]")
+{
+   dimension = GENERATE(2, 3);
+   int order = GENERATE(2, 3, 4);
+   int lor = GENERATE(2, 3, 4, 5);
+   auto vectorspace = GENERATE(VecSpace::VectorH1, VecSpace::ND,
+                               VecSpace::RT);
+
+   Mesh mesh;
+   if (dimension == 2)
+   {
+      mesh = Mesh::LoadFromFile("../../data/square-nurbs.mesh");
+   }
+   else
+   {
+      mesh = Mesh::LoadFromFile("../../data/cube-nurbs.mesh");
+   }
+   for (int r = 0; r < 2; r++)
+   {
+      mesh.UniformRefinement();
+   }
+
+   FiniteElementCollection* h1_fec = new H1_FECollection(order, dimension);
+   FiniteElementCollection* nurbs_fec = nullptr;
+   int vdim = 1;
+   switch (vectorspace)
+   {
+      case VecSpace::VectorH1:
+         nurbs_fec = new NURBSFECollection(order);
+         vdim = dimension;
+         break;
+      case VecSpace::ND:
+         nurbs_fec = new NURBS_HCurlFECollection(order);
+         break;
+      case VecSpace::RT:
+         nurbs_fec = new NURBS_HDivFECollection(order);
+         break;
+      case VecSpace::H1:
+         mfem_error("Only for the vector case");
+   }
+
+   // Define NURBS gridfunction
+   FiniteElementSpace nurbs_fes(&mesh, new NURBSExtension(mesh.NURBSext, order),
+                                nurbs_fec, vdim);
+   GridFunction nurbs_gf(&nurbs_fes);
+
+   // Define H1 gridfunction on refined mesh
+   Mesh h1_mesh = Mesh::MakeRefined(mesh, lor, BasisType::GaussLobatto);
+   H1_FECollection h1_coll(order, mesh.Dimension());
+   FiniteElementSpace h1_fes(&mesh, h1_fec, dimension);
+   GridFunction h1_gf(&h1_fes);
+
+   // Project coefficient of NURBS gridfunction
+   CartesianXCoefficient xcf;
+   CartesianYCoefficient ycf;
+   ProductCoefficient  pcf(xcf, ycf);
+   VectorArrayCoefficient cf(dimension);
+   cf.Set(0, &pcf, false);
+   cf.Set(1, &pcf, false);
+   if ( dimension == 3 ) { cf.Set(2, &pcf, false); }
+   nurbs_gf.ProjectCoefficient(cf);
+
+   // Transfer NURBS gridfunction to H1 gridfunction
+   OperatorHandle nurbs_to_h1;
+   h1_fes.GetTransferOperator(nurbs_fes, nurbs_to_h1);
+   nurbs_to_h1.Ptr()->Mult(nurbs_gf, h1_gf);
+
+   // Check
+   REQUIRE(h1_gf.ComputeL2Error(cf) == MFEM_Approx(0.0));
+
+   // Transfer NURBS gridfunction to H1 gridfunction
+   nurbs_gf = 0.0;
+   nurbs_to_h1.Ptr()->MultTranspose(h1_gf, nurbs_gf);
+
+   // Check
+   REQUIRE(nurbs_gf).ComputeL2Error(cf) == MFEM_Approx(0.0));
+}*/
+
+
 #ifdef MFEM_USE_MPI
 
 TEST_CASE("Parallel Transfer", "[Transfer][Parallel]")