What is the proper way to extract a linear data set from mfem?

[mclarsen]

Here is the scenario. I am exporting an mfem mesh and grid functions to another visualization tool that only knows about linear elements. I would like to refine the elements some number of times and extract a linear data set. The methods in VisIt (blueprint reader and mfem reader) just iterate of the geometry and add cells and fields. This is sub-optimal since the shared connectivity between nodes is lost. What is the right way to do this? Perhaps use a linear finite element collection?

[tzanio]

Hi Matt,

If I understand correctly, you asking how to extract the "low-order refined" mesh with vertices at the high-order nodes (as opposed to the low-order mesh with vertices at the high-order vertices) and methods such as Mesh::PrintVTK(std::ostream &out, int ref, int field_data) are not sufficient because they print a disjoint mesh. Correct?

One way to achieve that is to use this constructor: https://github.com/mfem/mfem/blob/master/mesh/mesh.hpp#L576-L587, e.g. like this:

Mesh *mesh_lor = new Mesh(mesh, order, BasisType::GaussLobatto);

Hope this helps, Tzanio

[tzanio]

Here is a version of Example 1 that further clarifies the high-order (HO) -to- low-order-refined (LOR) construction: ex1-lor.cpp.zip.

The codes starts with a HO mesh and a HO H1 finite element space on it (for simplicity of same order), defines the correspondinG LOR mesh and a 1st order space on it, and uses both discretizations to solve the same PDE (simple Poisson).

The HO and LOR meshes are synchronized so the corresponding finite element spaces have the same degrees of freedom in this case, which could be transferred simply with the identity operator (also illustrated).

Let me know if you have any further questions.

@cyrush -- this will fix the level lines in VisIt.

//                                MFEM Example 1
//                             LOR version for Matt
//
// Compile with: make ex1-lor
//
// Sample runs:  ex1-lor -m ../data/square-disc.mesh
//               ex1-lor -m ../data/star.mesh
//               ex1-lor -m ../data/escher.mesh
//               ex1-lor -m ../data/fichera.mesh
//               ex1-lor -m ../data/square-disc-p2.vtk -o 2
//               ex1-lor -m ../data/square-disc-p3.mesh -o 3
//               ex1-lor -m ../data/square-disc-nurbs.mesh -o -1
//               ex1-lor -m ../data/disc-nurbs.mesh -o -1
//               ex1-lor -m ../data/pipe-nurbs.mesh -o -1
//               ex1-lor -m ../data/star-surf.mesh
//               ex1-lor -m ../data/square-disc-surf.mesh
//               ex1-lor -m ../data/inline-segment.mesh
//               ex1-lor -m ../data/amr-quad.mesh
//               ex1-lor -m ../data/amr-hex.mesh
//               ex1-lor -m ../data/fichera-amr.mesh
//               ex1-lor -m ../data/mobius-strip.mesh
//               ex1-lor -m ../data/mobius-strip.mesh -o -1 -sc
//
// Description:  This example code demonstrates the use of MFEM to define a
//               simple finite element discretization of the Laplace problem
//               -Delta u = 1 with homogeneous Dirichlet boundary conditions.
//               Specifically, we discretize using a FE space of the specified
//               order, or if order < 1 using an isoparametric/isogeometric
//               space (i.e. quadratic for quadratic curvilinear mesh, NURBS for
//               NURBS mesh, etc.)
//
//               The example highlights the use of mesh refinement, finite
//               element grid functions, as well as linear and bilinear forms
//               corresponding to the left-hand side and right-hand side of the
//               discrete linear system. We also cover the explicit elimination
//               of essential boundary conditions, static condensation, and the
//               optional connection to the GLVis tool for visualization.

#include "mfem.hpp"
#include <fstream>
#include <iostream>

using namespace std;
using namespace mfem;

int main(int argc, char *argv[])
{
   // 1. Parse command-line options.
   const char *mesh_file = "../data/star-q3.mesh";
   int order = 1;
   bool static_cond = false;
   bool visualization = 1;

   OptionsParser args(argc, argv);
   args.AddOption(&mesh_file, "-m", "--mesh",
                  "Mesh file to use.");
   args.AddOption(&order, "-o", "--order",
                  "Finite element order (polynomial degree) or -1 for"
                  " isoparametric space.");
   args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
                  "--no-static-condensation", "Enable static condensation.");
   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
                  "--no-visualization",
                  "Enable or disable GLVis visualization.");
   args.Parse();
   if (!args.Good())
   {
      args.PrintUsage(cout);
      return 1;
   }
   args.PrintOptions(cout);

   // 2. Read the mesh from the given mesh file. We can handle triangular,
   //    quadrilateral, tetrahedral, hexahedral, surface and volume meshes with
   //    the same code.
   Mesh *mesh = new Mesh(mesh_file, 1, 1);
   int dim = mesh->Dimension();

   // 3. Refine the mesh to increase the resolution. In this example we do
   //    'ref_levels' of uniform refinement. We choose 'ref_levels' to be the
   //    largest number that gives a final mesh with no more than 50,000
   //    elements.
   {
      int ref_levels = 0;
      for (int l = 0; l < ref_levels; l++)
      {
         mesh->UniformRefinement();
      }
   }

   // For simplicity, enforce the same mesh & solution order. This is not a
   // requirement, any GridFunction could be transferred between the HO and LOR
   // mesh, see FiniteElementSpace::GetTransferOperator.
   int mesh_order =  mesh->GetNodalFESpace()->GetOrder(0);
   if (mesh_order > order)
   {
      order = mesh_order;
   }
   if (mesh_order < order)
   {
      mesh_order = order;
      mesh->SetCurvature(mesh_order);
   }

   // 4. Define a finite element space on the mesh. Here we use continuous
   //    Lagrange finite elements of the specified order. If order < 1, we
   //    instead use an isoparametric/isogeometric space.
   FiniteElementCollection *fec;
   if (order > 0)
   {
      fec = new H1_FECollection(order, dim);
   }
   else if (mesh->GetNodes())
   {
      fec = mesh->GetNodes()->OwnFEC();
      cout << "Using isoparametric FEs: " << fec->Name() << endl;
   }
   else
   {
      fec = new H1_FECollection(order = 1, dim);
   }
   FiniteElementSpace *fespace = new FiniteElementSpace(mesh, fec);
   cout << "Number of HO finite element unknowns: "
        << fespace->GetTrueVSize() << endl;

   // Define the LOR mesh and finite element space. In the case of equal orders
   // and default basis choices, the HO and LOR spaces will have the same number
   // of dofs and "solutions" could be transfered between them with the identity
   // operator (see below).
   Mesh *mesh_lor = new Mesh(mesh, mesh_order, BasisType::GaussLobatto);
   FiniteElementCollection *fec_lor = new H1_FECollection(1, dim);
   FiniteElementSpace *fespace_lor = new FiniteElementSpace(mesh_lor, fec_lor);
   cout << "Number of LOR finite element unknowns: "
        << fespace_lor->GetTrueVSize() << endl;

   // HO solve -- the original Example 1 code
   {
      // 5. Determine the list of true (i.e. conforming) essential boundary dofs.
      //    In this example, the boundary conditions are defined by marking all
      //    the boundary attributes from the mesh as essential (Dirichlet) and
      //    converting them to a list of true dofs.
      Array<int> ess_tdof_list;
      if (mesh->bdr_attributes.Size())
      {
         Array<int> ess_bdr(mesh->bdr_attributes.Max());
         ess_bdr = 1;
         fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
      }

      // 6. Set up the linear form b(.) which corresponds to the right-hand side of
      //    the FEM linear system, which in this case is (1,phi_i) where phi_i are
      //    the basis functions in the finite element fespace.
      LinearForm *b = new LinearForm(fespace);
      ConstantCoefficient one(1.0);
      b->AddDomainIntegrator(new DomainLFIntegrator(one));
      b->Assemble();

      // 7. Define the solution vector x as a finite element grid function
      //    corresponding to fespace. Initialize x with initial guess of zero,
      //    which satisfies the boundary conditions.
      GridFunction x(fespace);
      x = 0.0;

      // 8. Set up the bilinear form a(.,.) on the finite element space
      //    corresponding to the Laplacian operator -Delta, by adding the Diffusion
      //    domain integrator.
      BilinearForm *a = new BilinearForm(fespace);
      a->AddDomainIntegrator(new DiffusionIntegrator(one));

      // 9. Assemble the bilinear form and the corresponding linear system,
      //    applying any necessary transformations such as: eliminating boundary
      //    conditions, applying conforming constraints for non-conforming AMR,
      //    static condensation, etc.
      if (static_cond) { a->EnableStaticCondensation(); }
      a->Assemble();

      SparseMatrix A;
      Vector B, X;
      a->FormLinearSystem(ess_tdof_list, x, *b, A, X, B);

      cout << "Size of linear system: " << A.Height() << endl;

#ifndef MFEM_USE_SUITESPARSE
      // 10. Define a simple symmetric Gauss-Seidel preconditioner and use it to
      //     solve the system A X = B with PCG.
      GSSmoother M(A);
      PCG(A, M, B, X, 1, 200, 1e-12, 0.0);
#else
      // 10. If MFEM was compiled with SuiteSparse, use UMFPACK to solve the system.
      UMFPackSolver umf_solver;
      umf_solver.Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
      umf_solver.SetOperator(A);
      umf_solver.Mult(B, X);
#endif

      // 11. Recover the solution as a finite element grid function.
      a->RecoverFEMSolution(X, *b, x);

      // 12. Save the refined mesh and the solution. This output can be viewed later
      //     using GLVis: "glvis -m refined.mesh -g sol.gf".
      ofstream mesh_ofs("refined.mesh");
      mesh_ofs.precision(8);
      mesh->Print(mesh_ofs);
      ofstream sol_ofs("sol.gf");
      sol_ofs.precision(8);
      x.Save(sol_ofs);

      // 13. Send the solution by socket to a GLVis server.
      if (visualization)
      {
         char vishost[] = "localhost";
         int  visport   = 19916;
         socketstream sol_sock(vishost, visport);
         sol_sock.precision(8);
         sol_sock << "solution\n" << *mesh << x << flush
                  << "plot_caption 'HO Mesh + HO Solution'" << endl;
      }

      // Use the HO solution as an approximation on the LOR mesh
      if (visualization)
      {
         char vishost[] = "localhost";
         int  visport   = 19916;
         socketstream sol_sock(vishost, visport);
         sol_sock.precision(8);

         // Transfer HO solution as a GridFunction on the LOR space
         GridFunction x_lor(fespace_lor);
         x_lor = x;

         sol_sock << "solution\n" << *mesh_lor << x_lor << flush
                  << "plot_caption 'LOR Mesh + HO Solution'" << endl;
      }

      // 14. Free the used memory.
      delete a;
      delete b;
   }

   // LOR solve -- Repeat the solve with the LOR mesh and space
   {
      // 5. Determine the list of true (i.e. conforming) essential boundary dofs.
      //    In this example, the boundary conditions are defined by marking all
      //    the boundary attributes from the mesh as essential (Dirichlet) and
      //    converting them to a list of true dofs.
      Array<int> ess_tdof_list;
      if (mesh_lor->bdr_attributes.Size())
      {
         Array<int> ess_bdr(mesh_lor->bdr_attributes.Max());
         ess_bdr = 1;
         fespace_lor->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
      }

      // 6. Set up the linear form b(.) which corresponds to the right-hand side of
      //    the FEM linear system, which in this case is (1,phi_i) where phi_i are
      //    the basis functions in the finite element fespace_lor.
      LinearForm *b = new LinearForm(fespace_lor);
      ConstantCoefficient one(1.0);
      b->AddDomainIntegrator(new DomainLFIntegrator(one));
      b->Assemble();

      // 7. Define the solution vector x as a finite element grid function
      //    corresponding to fespace_lor. Initialize x with initial guess of zero,
      //    which satisfies the boundary conditions.
      GridFunction x(fespace_lor);
      x = 0.0;

      // 8. Set up the bilinear form a(.,.) on the finite element space
      //    corresponding to the Laplacian operator -Delta, by adding the Diffusion
      //    domain integrator.
      BilinearForm *a = new BilinearForm(fespace_lor);
      a->AddDomainIntegrator(new DiffusionIntegrator(one));

      // 9. Assemble the bilinear form and the corresponding linear system,
      //    applying any necessary transformations such as: eliminating boundary
      //    conditions, applying conforming constraints for non-conforming AMR,
      //    static condensation, etc.
      if (static_cond) { a->EnableStaticCondensation(); }
      a->Assemble();

      SparseMatrix A;
      Vector B, X;
      a->FormLinearSystem(ess_tdof_list, x, *b, A, X, B);

      cout << "Size of linear system: " << A.Height() << endl;

#ifndef MFEM_USE_SUITESPARSE
      // 10. Define a simple symmetric Gauss-Seidel preconditioner and use it to
      //     solve the system A X = B with PCG.
      GSSmoother M(A);
      PCG(A, M, B, X, 1, 200, 1e-12, 0.0);
#else
      // 10. If MFEM was compiled with SuiteSparse, use UMFPACK to solve the system.
      UMFPackSolver umf_solver;
      umf_solver.Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
      umf_solver.SetOperator(A);
      umf_solver.Mult(B, X);
#endif

      // 11. Recover the solution as a finite element grid function.
      a->RecoverFEMSolution(X, *b, x);

      // 12. Save the refined mesh and the solution. This output can be viewed later
      //     using GLVis: "glvis -m refined.mesh -g sol.gf".
      ofstream mesh_lor_ofs("refined_lor.mesh");
      mesh_lor_ofs.precision(8);
      mesh_lor->Print(mesh_lor_ofs);
      ofstream sol_ofs("sol_lor.gf");
      sol_ofs.precision(8);
      x.Save(sol_ofs);

      // 13. Send the solution by socket to a GLVis server.
      if (visualization)
      {
         char vishost[] = "localhost";
         int  visport   = 19916;
         socketstream sol_sock(vishost, visport);
         sol_sock.precision(8);
         sol_sock << "solution\n" << *mesh_lor << x << flush
                  << "plot_caption 'LOR Mesh + LOR Solution'" << endl;
      }

      // Use the LOR solution as an approximation on the HO mesh
      if (visualization)
      {
         char vishost[] = "localhost";
         int  visport   = 19916;
         socketstream sol_sock(vishost, visport);
         sol_sock.precision(8);

         // Transfer LOR solution as a GridFunction on the HO space
         GridFunction x_ho(fespace);
         x_ho = x;

         sol_sock << "solution\n" << *mesh << x_ho << flush
                  << "plot_caption 'HO Mesh + LOR Solution'" << endl;
      }

      // 14. Free the used memory.
      delete a;
      delete b;
   }

   delete fespace;
   if (order > 0) { delete fec; }
   delete mesh;

   delete fespace_lor;
   if (order > 0) { delete fec_lor; }
   delete mesh_lor;

   return 0;
}

[mclarsen]

I had found the constructor the creates the refined linear mesh. What was unclear to me was how to map the un-refined higher order grid functions (fields) back onto the refined linear mesh. I am new to this space so there are undoubtably concepts that I do not understand. I was pointed to some code in MARBL that does what I am after so I am going to see if I can get that to work.

To be a little clearer, I want to run example one as high order, but save out a low order solution.

[tzanio]

This is the LOR mesh + HO solution above.