add lshape example

examples/drl_lshape.cpp

@@ -0,0 +1,340 @@
+// Compile with: make drl_shock_wave
+//
+// drl_shock_wave -o 2 -m ../data/inline-quad.mesh
+// for multi agent local, set the mesh to use 20x20 grid because that is what 
+// was used for training.
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+#include "multi_agent_local_refiner.hpp"
+
+#define MFEM_USE_RLLIB
+#ifdef MFEM_USE_RLLIB
+#include <Python.h>
+#include "numpy/arrayobject.h"
+#endif
+
+using namespace std;
+using namespace mfem;
+
+#define alpha 2.0/3.0
+
+double exact_vel(const Vector &x)
+{
+   double xv = x(0), yv = x(1);
+   double rv = xv*xv + yv*yv;
+   if (rv > 0) { rv = pow(rv, 0.5); };
+   double theta = atan2(yv, xv);
+   if (theta < 0.0) { theta += 2*M_PI; }
+   return pow(rv, alpha)*sin(alpha*theta);
+}
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command-line options.
+   const char *mesh_file = "lshape.mesh";
+   int order = 2;
+   bool pa = false;
+   const char *device_config = "cpu";
+   bool visualization = true;
+   int jobid = 0;
+   double error_threshold = 0.10;
+   double max_elem_error = 5.0e-3;
+   int refinement_levels = 2;
+
+
+#ifdef MFEM_USE_RLLIB
+   Py_Initialize();
+   import_array(); // numpy init
+#endif
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree).");
+   args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
+                  "--no-partial-assembly", "Enable Partial Assembly.");
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&jobid, "-j", "--jobid",
+                  "slurb_jobid.");
+   args.AddOption(&error_threshold, "-err", "--err",
+                  "Total error fraction for zz or max_elem_error for policy.");
+   args.AddOption(&refinement_levels, "-r", "--ref",
+                 "Refinement levels");
+ 
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   args.PrintOptions(cout);
+
+   // 2. Enable hardware devices such as GPUs, and programming models such as
+   //    CUDA, OCCA, RAJA and OpenMP based on command line options.
+   Device device(device_config);
+   device.Print();
+
+   // 3. 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(mesh_file, 1, 1);
+   int dim = mesh.Dimension();
+   int sdim = mesh.SpaceDimension();
+
+   mesh.SetCurvature(2);
+
+   // 4. Since a NURBS mesh can currently only be refined uniformly, we need to
+   //    convert it to a piecewise-polynomial curved mesh. First we refine the
+   //    NURBS mesh a bit more and then project the curvature to quadratic Nodes.
+   if (mesh.NURBSext)
+   {
+      for (int i = 0; i < 2; i++)
+      {
+         mesh.UniformRefinement();
+      }
+      mesh.SetCurvature(2);
+   }
+   else {
+      // mesh.UniformRefinement();
+      //mesh.UniformRefinement();
+      for (int i = 0; i < refinement_levels; i++)
+      {
+         mesh.UniformRefinement();
+      }
+      mesh.EnsureNCMesh();
+   }
+
+   // 5. Define a finite element space on the mesh. The polynomial order is
+   //    one (linear) by default, but this can be changed on the command line.
+   H1_FECollection fec(order, dim);
+   FiniteElementSpace fespace(&mesh, &fec);
+
+   // Create 0-order L2 gridfunction to hold errors
+   L2_FECollection fec0(0, dim);
+   FiniteElementSpace fes0(&mesh, &fec0);
+   GridFunction err(&fes0);
+
+   // 6. As in Example 1, we set up bilinear and linear forms corresponding to
+   //    the Laplace problem -\Delta u = 1. We don't assemble the discrete
+   //    problem yet, this will be done in the main loop.
+   BilinearForm a(&fespace);
+   if (pa)
+   {
+      a.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+      a.SetDiagonalPolicy(Operator::DIAG_ONE);
+   }
+   LinearForm b(&fespace);
+
+   ConstantCoefficient rhs(0.0);
+   ConstantCoefficient one(1.0);
+   FunctionCoefficient exact(exact_vel);
+
+   BilinearFormIntegrator *integ = new DiffusionIntegrator(one);
+   a.AddDomainIntegrator(integ);
+   int int_order = 8;
+   int geom_type = mesh.GetElementBaseGeometry(0);
+   DomainLFIntegrator* dlfi = new DomainLFIntegrator(rhs);
+   dlfi->SetIntRule(&IntRules.Get(geom_type, int_order));
+   b.AddDomainIntegrator(dlfi);
+
+   // 7. The solution vector x and the associated finite element grid function
+   //    will be maintained over the AMR iterations. We initialize it to zero.
+   GridFunction x(&fespace);
+   x = 0.0;
+
+   // 8. All boundary attributes will be used for essential (Dirichlet) BC.
+   MFEM_VERIFY(mesh.bdr_attributes.Size() > 0,
+               "Boundary attributes required in the mesh.");
+   Array<int> ess_bdr(mesh.bdr_attributes.Max());
+   ess_bdr = 1;
+
+   // 9. Connect to GLVis.
+   char vishost[] = "localhost";
+   int  visport   = 19916;
+   socketstream sol_sock;
+   socketstream err_sock;
+   if (visualization)
+   {
+      sol_sock.open(vishost, visport);
+      err_sock.open(vishost, visport);
+   }
+
+   // 10. Set up an error estimator. Here we use the Zienkiewicz-Zhu estimator
+   //     that uses the ComputeElementFlux method of the DiffusionIntegrator to
+   //     recover a smoothed flux (gradient) that is subtracted from the element
+   //     flux to get an error indicator. We need to supply the space for the
+   //     smoothed flux: an (H1)^sdim (i.e., vector-valued) space is used here.
+   FiniteElementSpace flux_fespace(&mesh, &fec, sdim);
+   ZienkiewiczZhuEstimator estimator(*integ, x, flux_fespace);
+   //KellyErrorEstimator estimator2(*integ, x, flux_fespace);
+   //estimator.SetAnisotropic();
+
+   // 11. A refiner selects and refines elements based on a refinement strategy.
+   //     The strategy here is to refine elements with errors larger than a
+   //     fraction of the maximum element error. Other strategies are possible.
+   //     The refiner will call the given error estimator.
+
+   bool zz = false;
+#if 0
+   ThresholdRefiner refiner(estimator);
+   refiner.SetTotalErrorFraction(error_threshold);
+   zz = true;
+#else
+   MAL_DRLRefiner refiner(x, error_threshold);
+#endif
+
+   bool derefine = false;
+   ThresholdDerefiner derefiner(estimator);
+   derefiner.SetThreshold(0.05);
+   derefiner.SetNCLimit(0);
+
+   string errorfilename;
+   errorfilename = to_string(jobid) + "_lshape_error.txt";
+   ofstream myfile;
+   myfile.open(errorfilename, ofstream::in | ofstream::out | ofstream::app);
+
+
+   // 12. The main AMR loop. In each iteration we solve the problem on the
+   //     current mesh, visualize the solution, and refine the mesh.
+   const int max_dofs = 200000;
+   for (int it = 0; it < 6; it++)
+   {
+      int cdofs = fespace.GetTrueVSize();
+      cout << "\nAMR iteration " << it << endl;
+      cout << "Number of unknowns: " << cdofs << endl;
+
+      // 13. Assemble the right-hand side.
+      b.Assemble();
+
+      // 14. Set Dirichlet boundary values in the GridFunction x.
+      //     Determine the list of Dirichlet true DOFs in the linear system.
+      Array<int> ess_tdof_list;
+      x.ProjectBdrCoefficient(exact, ess_bdr);
+      fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+
+      // 15. Assemble the stiffness matrix.
+      a.Assemble();
+      
+      // 16. Create the linear system: eliminate boundary conditions, constrain
+      //     hanging nodes and possibly apply other transformations. The system
+      //     will be solved for true (unconstrained) DOFs only.
+      OperatorPtr A;
+      Vector B, X;
+
+      const int copy_interior = 1;
+      a.FormLinearSystem(ess_tdof_list, x, b, A, X, B, copy_interior);
+
+      // 17. Solve the linear system A X = B.
+      if (!pa)
+      {
+#ifndef MFEM_USE_SUITESPARSE
+         // Use a simple symmetric Gauss-Seidel preconditioner with PCG.
+         GSSmoother M((SparseMatrix&)(*A));
+         PCG(*A, M, B, X, 3, 200, 1e-12, 0.0);
+#else
+         // 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
+      }
+      else // Diagonal preconditioning in partial assembly mode.
+      {
+         OperatorJacobiSmoother M(a, ess_tdof_list);
+         PCG(*A, M, B, X, 3, 2000, 1e-12, 0.0);
+      }
+
+      // 18. After solving the linear system, reconstruct the solution as a
+      //     finite element GridFunction. Constrained nodes are interpolated
+      //     from true DOFs (it may therefore happen that x.Size() >= X.Size()).
+      a.RecoverFEMSolution(X, b, x);
+
+      // Compute error against exact solution
+
+      x.ComputeElementL2Errors(exact, err);
+      int int_order = std::max(20 - it, 2*order+1);
+      double error;
+      error = err.Norml2();
+
+      if (derefine) {
+         myfile << error_threshold << " " << cdofs << " " << error << endl;
+      }
+      else {
+         myfile << -error_threshold << " " << cdofs << " " << error << endl;
+      }
+
+      // 19. Send solution by socket to the GLVis server.
+      if (visualization && sol_sock.good())
+      {
+         sol_sock.precision(8);
+         sol_sock << "solution\n" << mesh << x << flush;
+      }
+      if (visualization && err_sock.good())
+      {
+         err_sock.precision(8);
+         err_sock << "solution\n" << mesh << err << flush;
+      }
+
+      if (cdofs > max_dofs)
+      {
+         cout << "Reached the maximum number of dofs. Stop." << endl;
+         break;
+      }
+
+      // 20. Call the refiner to modify the mesh. The refiner calls the error
+      //     estimator to obtain element errors, then it selects elements to be
+      //     refined and finally it modifies the mesh. The Stop() method can be
+      //     used to determine if a stopping criterion was met.
+
+      refiner.Apply(mesh);
+      if (refiner.Stop())
+      {
+         cout << "Stopping criterion satisfied. Stop." << endl;
+         break;
+      }
+      fespace.Update();fes0.Update();x.Update();err.Update();
+      a.Update();b.Update();
+
+      if (derefine) {
+         derefiner.Apply(mesh);
+         fespace.Update();fes0.Update();x.Update();err.Update();
+         a.Update();b.Update();
+      }
+e
+      {
+         string solname = to_string(jobid) + "_lshape_amr" + to_string(it) + ".gf";
+         ofstream sol_ofs(solname);
+         x.Save(sol_ofs);
+      }
+
+      {
+         string meshname = to_string(jobid) + "_lshape_amr" + to_string(it) + ".mesh";
+         ofstream mesh_ofs(meshname);
+         mesh_ofs.precision(14);
+         mesh.Print(mesh_ofs);
+      }
+   }
+   myfile.close();
+
+
+   {
+      ofstream sol_ofs("lshape_amr.gf");
+      x.Save(sol_ofs);
+   }
+
+   {
+      ofstream mesh_ofs("lshape_amr.mesh");
+      mesh_ofs.precision(14);
+      mesh.Print(mesh_ofs);
+   }
+
+   return 0;
+}

examples/drl_shock_wave.cpp

@@ -0,0 +1,330 @@
+// Compile with: make drl_shock_wave
+//
+// drl_shock_wave -o 2 -m ../data/inline-quad.mesh
+// for multi agent local, set the mesh to use 20x20 grid because that is what 
+// was used for training.
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+#include "drl_shock_wave.hpp"
+#include "multi_agent_local_refiner.hpp"
+
+#define MFEM_USE_RLLIB
+#ifdef MFEM_USE_RLLIB
+#include <Python.h>
+#include "numpy/arrayobject.h"
+#endif
+
+using namespace std;
+using namespace mfem;
+
+double bdry_function(const Vector& x)
+{
+   return 0.0; // default
+}
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command-line options.
+   const char *mesh_file = "../data/inline-quad.mesh";
+   int order = 2;
+   bool pa = false;
+   const char *device_config = "cpu";
+   bool visualization = true;
+   int jobid = 0;
+
+#ifdef MFEM_USE_RLLIB
+   Py_Initialize();
+   import_array(); // numpy init
+#endif
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree).");
+   args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
+                  "--no-partial-assembly", "Enable Partial Assembly.");
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&jobid, "-j", "--jobid",
+                  "slurb_jobid.");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   args.PrintOptions(cout);
+
+   // 2. Enable hardware devices such as GPUs, and programming models such as
+   //    CUDA, OCCA, RAJA and OpenMP based on command line options.
+   Device device(device_config);
+   device.Print();
+
+   // 3. 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(mesh_file, 1, 1);
+   int dim = mesh.Dimension();
+   int sdim = mesh.SpaceDimension();
+
+   mesh.SetCurvature(2);
+
+   // 4. Since a NURBS mesh can currently only be refined uniformly, we need to
+   //    convert it to a piecewise-polynomial curved mesh. First we refine the
+   //    NURBS mesh a bit more and then project the curvature to quadratic Nodes.
+   if (mesh.NURBSext)
+   {
+      for (int i = 0; i < 2; i++)
+      {
+         mesh.UniformRefinement();
+      }
+      mesh.SetCurvature(2);
+   }
+   else {
+      // mesh.UniformRefinement();
+      //mesh.UniformRefinement();
+      mesh.EnsureNCMesh();
+   }
+
+   // 5. Define a finite element space on the mesh. The polynomial order is
+   //    one (linear) by default, but this can be changed on the command line.
+   H1_FECollection fec(order, dim);
+   FiniteElementSpace fespace(&mesh, &fec);
+
+   // Create 0-order L2 gridfunction to hold errors
+   L2_FECollection fec0(0, dim);
+   FiniteElementSpace fes0(&mesh, &fec0);
+   GridFunction err(&fes0);
+
+   // 6. As in Example 1, we set up bilinear and linear forms corresponding to
+   //    the Laplace problem -\Delta u = 1. We don't assemble the discrete
+   //    problem yet, this will be done in the main loop.
+   BilinearForm a(&fespace);
+   if (pa)
+   {
+      a.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+      a.SetDiagonalPolicy(Operator::DIAG_ONE);
+   }
+   LinearForm b(&fespace);
+
+   FunctionCoefficient bdry(bdry_function);
+   FunctionCoefficient rhs(layer2_laplace);
+   FunctionCoefficient exact(layer2_exsol);
+   ConstantCoefficient one(1.0);
+
+   BilinearFormIntegrator *integ = new DiffusionIntegrator(one);
+   a.AddDomainIntegrator(integ);
+   int int_order = 8;
+   int geom_type = mesh.GetElementBaseGeometry(0);
+   DomainLFIntegrator* dlfi = new DomainLFIntegrator(rhs);
+   dlfi->SetIntRule(&IntRules.Get(geom_type, int_order));
+   b.AddDomainIntegrator(dlfi);
+
+   // 7. The solution vector x and the associated finite element grid function
+   //    will be maintained over the AMR iterations. We initialize it to zero.
+   GridFunction x(&fespace);
+   x = 0.0;
+
+   // 8. All boundary attributes will be used for essential (Dirichlet) BC.
+   MFEM_VERIFY(mesh.bdr_attributes.Size() > 0,
+               "Boundary attributes required in the mesh.");
+   Array<int> ess_bdr(mesh.bdr_attributes.Max());
+   ess_bdr = 1;
+
+   // 9. Connect to GLVis.
+   char vishost[] = "localhost";
+   int  visport   = 19916;
+   socketstream sol_sock;
+   socketstream err_sock;
+   if (visualization)
+   {
+      sol_sock.open(vishost, visport);
+      err_sock.open(vishost, visport);
+   }
+
+   // 10. Set up an error estimator. Here we use the Zienkiewicz-Zhu estimator
+   //     that uses the ComputeElementFlux method of the DiffusionIntegrator to
+   //     recover a smoothed flux (gradient) that is subtracted from the element
+   //     flux to get an error indicator. We need to supply the space for the
+   //     smoothed flux: an (H1)^sdim (i.e., vector-valued) space is used here.
+   FiniteElementSpace flux_fespace(&mesh, &fec, sdim);
+   ZienkiewiczZhuEstimator estimator(*integ, x, flux_fespace);
+   //KellyErrorEstimator estimator2(*integ, x, flux_fespace);
+   //estimator.SetAnisotropic();
+
+   // 11. A refiner selects and refines elements based on a refinement strategy.
+   //     The strategy here is to refine elements with errors larger than a
+   //     fraction of the maximum element error. Other strategies are possible.
+   //     The refiner will call the given error estimator.
+
+#if 0
+   ThresholdRefiner refiner(estimator);
+   refiner.SetTotalErrorFraction(0.10);
+#else
+   MAL_DRLRefiner refiner(x);
+   //DRLRefiner refiner(x);
+#endif
+  
+#if 0
+   x.ProjectCoefficient(exact);
+   refiner.Apply(mesh);
+   fespace.Update();
+   fes0.Update();
+   x.Update();
+   x.ProjectCoefficient(exact);
+   refiner.Apply(mesh);
+   fespace.Update();
+   fes0.Update();
+   x.Update();
+   {
+      ofstream sol_ofs("sol.gf");
+      x.Save(sol_ofs);
+   }
+   {
+      ofstream mesh_ofs("amr.mesh");
+      mesh_ofs.precision(14);
+      mesh.Print(mesh_ofs);
+   }
+   MFEM_ABORT(" ");
+#endif
+
+   // 12. The main AMR loop. In each iteration we solve the problem on the
+   //     current mesh, visualize the solution, and refine the mesh.
+   const int max_dofs = 20000;
+   for (int it = 0; it < 4; it++)
+   {
+      int cdofs = fespace.GetTrueVSize();
+      cout << "\nAMR iteration " << it << endl;
+      cout << "Number of unknowns: " << cdofs << endl;
+
+      // 13. Assemble the right-hand side.
+      b.Assemble();
+
+      // 14. Set Dirichlet boundary values in the GridFunction x.
+      //     Determine the list of Dirichlet true DOFs in the linear system.
+      Array<int> ess_tdof_list;
+      x.ProjectBdrCoefficient(exact, ess_bdr);
+      fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+
+      // 15. Assemble the stiffness matrix.
+      a.Assemble();
+      
+      // 16. Create the linear system: eliminate boundary conditions, constrain
+      //     hanging nodes and possibly apply other transformations. The system
+      //     will be solved for true (unconstrained) DOFs only.
+      OperatorPtr A;
+      Vector B, X;
+
+      const int copy_interior = 1;
+      a.FormLinearSystem(ess_tdof_list, x, b, A, X, B, copy_interior);
+
+      // 17. Solve the linear system A X = B.
+      if (!pa)
+      {
+#ifndef MFEM_USE_SUITESPARSE
+         // Use a simple symmetric Gauss-Seidel preconditioner with PCG.
+         GSSmoother M((SparseMatrix&)(*A));
+         PCG(*A, M, B, X, 3, 200, 1e-12, 0.0);
+#else
+         // 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
+      }
+      else // Diagonal preconditioning in partial assembly mode.
+      {
+         OperatorJacobiSmoother M(a, ess_tdof_list);
+         PCG(*A, M, B, X, 3, 2000, 1e-12, 0.0);
+      }
+
+      // 18. After solving the linear system, reconstruct the solution as a
+      //     finite element GridFunction. Constrained nodes are interpolated
+      //     from true DOFs (it may therefore happen that x.Size() >= X.Size()).
+      a.RecoverFEMSolution(X, b, x);
+
+      // Compute error against exact solution
+
+      x.ComputeElementL2Errors(exact, err);
+
+      // 19. Send solution by socket to the GLVis server.
+      if (visualization && sol_sock.good())
+      {
+         sol_sock.precision(8);
+         sol_sock << "solution\n" << mesh << x << flush;
+      }
+      if (visualization && err_sock.good())
+      {
+         err_sock.precision(8);
+         err_sock << "solution\n" << mesh << err << flush;
+      }
+
+      if (cdofs > max_dofs || it == 3)
+      {
+         cout << "Reached the maximum number of dofs. Stop." << endl;
+         break;
+      }
+
+      // 20. Call the refiner to modify the mesh. The refiner calls the error
+      //     estimator to obtain element errors, then it selects elements to be
+      //     refined and finally it modifies the mesh. The Stop() method can be
+      //     used to determine if a stopping criterion was met.
+
+      refiner.Apply(mesh);
+      if (refiner.Stop())
+      {
+         cout << "Stopping criterion satisfied. Stop." << endl;
+         break;
+      }
+
+      // 21. Update the space to reflect the new state of the mesh. Also,
+      //     interpolate the solution x so that it lies in the new space but
+      //     represents the same function. This saves solver iterations later
+      //     since we'll have a good initial guess of x in the next step.
+      //     Internally, FiniteElementSpace::Update() calculates an
+      //     interpolation matrix which is then used by GridFunction::Update().
+      fespace.Update();
+      fes0.Update();
+      x.Update();
+      err.Update();
+
+      {
+         string solname = to_string(jobid) + "_amrsol" + to_string(it) + ".gf";
+         ofstream sol_ofs(solname);
+         x.Save(sol_ofs);
+      }
+
+      {
+         string meshname = to_string(jobid) + "_amr" + to_string(it) + ".mesh";
+         ofstream mesh_ofs(meshname);
+         mesh_ofs.precision(14);
+         mesh.Print(mesh_ofs);
+      }
+
+      // 22. Inform also the bilinear and linear forms that the space has
+      //     changed.
+      a.Update();
+      b.Update();
+   }
+
+
+   {
+      ofstream sol_ofs("amrsol.gf");
+      x.Save(sol_ofs);
+   }
+
+   {
+      ofstream mesh_ofs("amr.mesh");
+      mesh_ofs.precision(14);
+      mesh.Print(mesh_ofs);
+   }
+
+   return 0;
+}

examples/drl_shock_wave.hpp

@@ -0,0 +1,136 @@
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace mfem;
+using namespace std;
+
+const double alpha =  200.0; // standard params
+const double center = -0.05;
+const double radius = 0.7;
+template<typename T> T sqr(T x) { return x*x; }
+double layer2_exsol(Vector &p)
+{
+   double x = p(0), y = p(1);
+   double r = sqrt(sqr(x - center) + sqr(y - center));
+   return atan(alpha * (r - radius));
+}
+
+void layer2_exgrad(const Vector &p, Vector &grad)
+{
+   double x = p(0), y = p(1);
+   double r = sqrt(sqr(x - center) + sqr(y - center));
+   double u = r * (sqr(alpha) * sqr(r - radius) + 1);
+   grad(0) = alpha * (x - center) / u;
+   grad(1) = alpha * (y - center) / u;
+}
+
+double layer2_laplace(Vector &p)
+{
+   double x = p(0), y = p(1);
+   double r = sqr(y - center) + sqr(x - center);
+   double u = sqr(alpha) * sqr(sqrt(r) - radius) + 1;
+
+   return 2 * pow(alpha,3) * (sqrt(r) - radius) * sqr(y - center) / (r * sqr(u))
+          + alpha * sqr(y - center) / (pow(r, 1.5) * u)
+          - 2 * alpha / (sqrt(r) * u)
+          + 2 * pow(alpha,3) * (sqrt(r) - radius) * sqr(x - center) / (r * sqr(u))
+          + alpha * sqr(x - center) / (pow(r, 1.5) * u);
+}
+
+static double safeSqrt(double x)
+{
+   if (x < 0.0)
+      return -sqrt(-x);
+   else
+      return sqrt(x);
+}
+
+double CalculateH10Error(GridFunction *sol, VectorCoefficient *exgrad,
+                         Array<double> *elemError, Array<int> *elemRef,
+                         int intOrder)
+{
+   const FiniteElementSpace *fes = sol->FESpace();
+   Mesh* mesh = fes->GetMesh();
+
+   Vector e_grad, a_grad, el_dofs, q_grad;
+   DenseMatrix dshape, dshapet, Jinv;
+   Array<int> vdofs;
+   const FiniteElement *fe;
+   ElementTransformation *transf;
+
+   int dim = mesh->Dimension();
+   e_grad.SetSize(dim);
+   a_grad.SetSize(dim);
+   q_grad.SetSize(dim);
+   Jinv.SetSize(dim);
+
+   double error = 0.0;
+   if (elemError) elemError->SetSize(mesh->GetNE());
+   if (elemRef)   elemRef->SetSize(mesh->GetNE());
+
+   for (int i = 0; i < mesh->GetNE(); i++)
+   {
+      fe = fes->GetFE(i);
+      int fdof = fe->GetDof();
+      transf = mesh->GetElementTransformation(i);
+      el_dofs.SetSize(fdof);
+      dshape.SetSize(fdof, dim);
+      dshapet.SetSize(fdof, dim);
+
+      fes->GetElementVDofs(i, vdofs);
+      for (int k = 0; k < fdof; k++)
+         if (vdofs[k] >= 0)
+            el_dofs(k) =  (*sol)(vdofs[k]);
+         else
+            el_dofs(k) = -(*sol)(-1-vdofs[k]);
+
+      const IntegrationRule &ir = IntRules.Get(fe->GetGeomType(), intOrder);
+
+      double el_err = 0.0, a_dxyz[3] = { 0, 0, 0 };
+      for (int j = 0; j < ir.GetNPoints(); j++)
+      {
+         const IntegrationPoint &ip = ir.IntPoint(j);
+
+         transf->SetIntPoint(&ip);
+         CalcInverse(transf->Jacobian(), Jinv);
+         double w = ip.weight * transf->Weight();
+
+         exgrad->Eval(e_grad, *transf, ip);
+
+         fe->CalcDShape(ip, dshape);
+         Mult(dshape, Jinv, dshapet);
+         dshapet.MultTranspose(el_dofs, a_grad);
+
+         e_grad -= a_grad;
+         el_err += w * (e_grad * e_grad);
+
+         transf->Jacobian().MultTranspose(e_grad, q_grad);
+         for (int k = 0; k < dim; k++)
+         {
+            a_dxyz[k] += w * (q_grad[k] * q_grad[k]);
+         }
+      }
+
+      error += el_err;
+      if (elemError)
+         (*elemError)[i] = sqrt(fabs(el_err));
+
+      if (elemRef)
+      {
+         double sum = 0;
+         for (int k = 0; k < dim; k++)
+            sum += a_dxyz[k];
+
+         const double thresh = 0.2 * 3/dim;
+         int ref = 0;
+         for (int k = 0; k < dim; k++)
+            if (a_dxyz[k] / sum > thresh)
+               ref |= (1 << k);
+
+         (*elemRef)[i] = ref;
+      }
+   }
+
+   return safeSqrt(error);
+}

examples/ex6.cpp

@@ -43,21 +43,159 @@
 #include <fstream>
 #include <iostream>
 
+#define MFEM_USE_RLLIB
+#ifdef MFEM_USE_RLLIB
+#include <Python.h>
+#include "numpy/arrayobject.h"
+#endif
+
 using namespace std;
 using namespace mfem;
 
+int problem; // problem number, controls source term and bdry condition
+
+// Returns either rhs function f, or exact solution u.
+
+double rhs_function(const Vector &x, bool exact = false)
+{
+   int dim = x.Size();
+
+   if (0 == problem) {
+
+      // NIST: "Peak 2D" problem
+      // u = e^(-alpha*(x^2+y^2)), alpha = 1e3 or alpha = 1e5
+
+      double xc = 0.5;
+      double yc = 0.5;
+
+      double x0 = x(0)-xc;
+      double y0 = x(1)-yc;
+
+      double alpha = 1000;
+
+      double xx = x0*x0;
+      double yy = y0*y0;
+
+      if (exact) {
+
+         double u = exp(-alpha*(xx+yy));
+
+         return u;
+      }
+      else {
+
+         double f =
+            alpha*(4.*alpha*xx -2.0)*exp(-alpha*(xx+yy)) +
+            alpha*(4.*alpha*yy -2.0)*exp(-alpha*(xx+yy));
+
+         return -f;
+      }
+   }
+   if (1 == problem) {
+
+      // NIST "arctan circular wavefront" problem, w/ minor
+      // modifications.
+
+      double r0 = 0.25;
+      double a = 100.0;
+      double h = 1./3.;
+      double c = 0.5;
+
+      double x0 = 0.5;
+      double y0 = 0.5;
+
+      double dx = x(0)-x0;
+      double dy = x(1)-y0;
+      double dxdx = dx*dx;
+      double dydy = dy*dy;
+      double r = sqrt(dxdx+dydy);
+
+      if (exact) {
+         double u = c +h*atan(a*(r-r0));
+         return u;
+      }
+      else {
+         double aa = a*a;
+         double dr = r-r0;
+         double drdr = dr*dr;
+         double rr = r*r;
+         double t = 1+aa*drdr;
+         double tt = t*t;
+
+         double fx =
+            a/(t*r) -
+            a*dxdx/(t*rr*r) -
+            2*aa*a*dxdx*dr/(tt*rr);
+
+         double fy =
+            a/(t*r) -
+            a*dydy/(t*rr*r) -
+            2*aa*a*dydy*dr/(tt*rr);
+
+         return -h*(fx+fy);
+      }
+   }
+   if (2 == problem) {
+
+      // cross-shaped source
+
+      double x0 = x(0)-0.5;
+      double y0 = x(1)-0.5;
+
+      double w1 = 0.04;
+      double w2 = 0.20;
+
+      if (x(0) > 0.5-w1 && x(0) < 0.5+w1 &&
+          x(1) > 0.5-w2 && x(1) < 0.5+w2 ) {
+         return 100.0;
+      }
+      if (x(1) > 0.5-w1 && x(1) < 0.5+w1 &&
+          x(0) > 0.5-w2 && x(0) < 0.5+w2 ) {
+         return 100.0;
+      }
+      return 0.0;
+   }
+
+   return 1.0; // default
+}
+
+double exact_soln(const Vector& x)
+{
+   return rhs_function(x, true);
+}
+
+double bdry_function(const Vector& x)
+{
+   if (problem == 0) {
+      return 0.0;
+   }
+   if (problem == 1) {
+      return 1.0;
+   }
+
+   return 0.0; // default
+}
+
 int main(int argc, char *argv[])
 {
    // 1. Parse command-line options.
+   problem = 0;
    const char *mesh_file = "../data/star.mesh";
    int order = 1;
    bool pa = false;
    const char *device_config = "cpu";
    bool visualization = true;
 
+#ifdef MFEM_USE_RLLIB
+   Py_Initialize();
+   import_array(); // numpy init
+#endif
+
    OptionsParser args(argc, argv);
    args.AddOption(&mesh_file, "-m", "--mesh",
                   "Mesh file to use.");
+   args.AddOption(&problem, "-p", "--problem",
+                  "Problem setup to use. See options in rhs_function().");
    args.AddOption(&order, "-o", "--order",
                   "Finite element order (polynomial degree).");
    args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
@@ -98,12 +236,22 @@ int main(int argc, char *argv[])
       }
       mesh.SetCurvature(2);
    }
+   else {
+      // mesh.UniformRefinement();
+      //mesh.UniformRefinement();
+      mesh.EnsureNCMesh();
+   }
 
    // 5. Define a finite element space on the mesh. The polynomial order is
    //    one (linear) by default, but this can be changed on the command line.
    H1_FECollection fec(order, dim);
    FiniteElementSpace fespace(&mesh, &fec);
 
+   // Create 0-order L2 gridfunction to hold errors
+   L2_FECollection fec0(0, dim);
+   FiniteElementSpace fes0(&mesh, &fec0);
+   GridFunction err(&fes0);
+
    // 6. As in Example 1, we set up bilinear and linear forms corresponding to
    //    the Laplace problem -\Delta u = 1. We don't assemble the discrete
    //    problem yet, this will be done in the main loop.
@@ -115,12 +263,14 @@ int main(int argc, char *argv[])
    }
    LinearForm b(&fespace);
 
+   FunctionCoefficient rhs(rhs_function);
+   FunctionCoefficient exact(exact_soln);
+   FunctionCoefficient bdry(bdry_function);
    ConstantCoefficient one(1.0);
-   ConstantCoefficient zero(0.0);
 
    BilinearFormIntegrator *integ = new DiffusionIntegrator(one);
    a.AddDomainIntegrator(integ);
-   b.AddDomainIntegrator(new DomainLFIntegrator(one));
+   b.AddDomainIntegrator(new DomainLFIntegrator(rhs));
 
    // 7. The solution vector x and the associated finite element grid function
    //    will be maintained over the AMR iterations. We initialize it to zero.
@@ -137,9 +287,11 @@ int main(int argc, char *argv[])
    char vishost[] = "localhost";
    int  visport   = 19916;
    socketstream sol_sock;
+   socketstream err_sock;
    if (visualization)
    {
       sol_sock.open(vishost, visport);
+      err_sock.open(vishost, visport);
    }
 
    // 10. Set up an error estimator. Here we use the Zienkiewicz-Zhu estimator
@@ -149,18 +301,24 @@ int main(int argc, char *argv[])
    //     smoothed flux: an (H1)^sdim (i.e., vector-valued) space is used here.
    FiniteElementSpace flux_fespace(&mesh, &fec, sdim);
    ZienkiewiczZhuEstimator estimator(*integ, x, flux_fespace);
-   estimator.SetAnisotropic();
+   //KellyErrorEstimator estimator2(*integ, x, flux_fespace);
+   //estimator.SetAnisotropic();
 
    // 11. A refiner selects and refines elements based on a refinement strategy.
    //     The strategy here is to refine elements with errors larger than a
    //     fraction of the maximum element error. Other strategies are possible.
    //     The refiner will call the given error estimator.
+
+#if 1
    ThresholdRefiner refiner(estimator);
-   refiner.SetTotalErrorFraction(0.7);
+   refiner.SetTotalErrorFraction(0.10);
+#else
+   DRLRefiner refiner(x);
+#endif
 
    // 12. The main AMR loop. In each iteration we solve the problem on the
    //     current mesh, visualize the solution, and refine the mesh.
-   const int max_dofs = 50000;
+   const int max_dofs = 1000;
    for (int it = 0; ; it++)
    {
       int cdofs = fespace.GetTrueVSize();
@@ -173,7 +331,7 @@ int main(int argc, char *argv[])
       // 14. Set Dirichlet boundary values in the GridFunction x.
       //     Determine the list of Dirichlet true DOFs in the linear system.
       Array<int> ess_tdof_list;
-      x.ProjectBdrCoefficient(zero, ess_bdr);
+      x.ProjectBdrCoefficient(bdry, ess_bdr);
       fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
 
       // 15. Assemble the stiffness matrix.
@@ -214,14 +372,23 @@ int main(int argc, char *argv[])
       //     from true DOFs (it may therefore happen that x.Size() >= X.Size()).
       a.RecoverFEMSolution(X, b, x);
 
+      // Compute error against exact solution
+
+      x.ComputeElementL2Errors(exact, err);
+
       // 19. Send solution by socket to the GLVis server.
       if (visualization && sol_sock.good())
       {
          sol_sock.precision(8);
          sol_sock << "solution\n" << mesh << x << flush;
       }
+      if (visualization && err_sock.good())
+      {
+         err_sock.precision(8);
+         err_sock << "solution\n" << mesh << err << flush;
+      }
 
-      if (cdofs > max_dofs)
+      if (cdofs > max_dofs || it == 3)
       {
          cout << "Reached the maximum number of dofs. Stop." << endl;
          break;
@@ -231,6 +398,7 @@ int main(int argc, char *argv[])
       //     estimator to obtain element errors, then it selects elements to be
       //     refined and finally it modifies the mesh. The Stop() method can be
       //     used to determine if a stopping criterion was met.
+
       refiner.Apply(mesh);
       if (refiner.Stop())
       {
@@ -245,7 +413,9 @@ int main(int argc, char *argv[])
       //     Internally, FiniteElementSpace::Update() calculates an
       //     interpolation matrix which is then used by GridFunction::Update().
       fespace.Update();
+      fes0.Update();
       x.Update();
+      err.Update();
 
       // 22. Inform also the bilinear and linear forms that the space has
       //     changed.

examples/lshape_amr.mesh

@@ -0,0 +1,640 @@
+MFEM mesh v1.1
+
+#
+# MFEM Geometry Types (see mesh/geom.hpp):
+#
+# POINT       = 0
+# SEGMENT     = 1
+# TRIANGLE    = 2
+# SQUARE      = 3
+# TETRAHEDRON = 4
+# CUBE        = 5
+# PRISM       = 6
+#
+
+dimension
+2
+
+elements
+84
+1 3 0 21 53 24
+1 3 21 8 22 53
+1 3 53 22 18 23
+1 3 24 53 23 11
+1 3 8 25 54 22
+1 3 25 1 26 54
+1 3 54 26 9 27
+1 3 22 54 27 18
+1 3 18 65 92 68
+1 3 68 92 67 30
+1 3 92 66 55 67
+1 3 65 27 66 92
+1 3 27 9 28 55
+1 3 55 69 93 72
+1 3 69 28 70 93
+1 3 93 70 4 71
+1 3 72 93 71 29
+1 3 95 72 29 73
+1 3 74 95 73 10
+1 3 30 67 95 74
+1 3 67 55 72 95
+1 3 11 23 56 32
+1 3 23 18 30 56
+1 3 56 30 10 31
+1 3 32 56 31 3
+1 3 3 31 57 35
+1 3 31 10 33 57
+1 3 57 33 19 34
+1 3 35 57 34 14
+1 3 10 73 96 75
+1 3 75 96 77 33
+1 3 96 76 58 77
+1 3 73 29 76 96
+1 3 29 71 97 76
+1 3 71 94 108 98
+1 3 94 110 114 109
+1 3 110 4 111 114
+1 3 114 111 99 112
+1 3 109 114 112 108
+1 3 108 99 78 100
+1 3 98 108 100 97
+1 3 97 78 36 79
+1 3 76 97 79 58
+1 3 58 79 101 82
+1 3 79 36 80 101
+1 3 101 80 12 81
+1 3 82 101 81 37
+1 3 33 58 37 19
+1 3 19 37 59 40
+1 3 37 12 38 59
+1 3 59 38 7 39
+1 3 40 59 39 13
+1 3 14 34 60 42
+1 3 34 19 40 60
+1 3 60 40 13 41
+1 3 42 60 41 6
+1 3 4 103 113 99
+1 3 103 83 104 113
+1 3 113 104 102 105
+1 3 99 113 105 78
+1 3 78 102 85 36
+1 3 102 84 61 85
+1 3 83 43 84 102
+1 3 43 86 106 84
+1 3 86 15 87 106
+1 3 106 87 44 88
+1 3 84 106 88 61
+1 3 61 44 20 45
+1 3 36 61 45 12
+1 3 15 89 107 87
+1 3 87 107 91 44
+1 3 107 90 62 91
+1 3 89 46 90 107
+1 3 46 5 47 62
+1 3 62 47 16 48
+1 3 44 62 48 20
+1 3 20 48 63 51
+1 3 48 16 49 63
+1 3 63 49 2 50
+1 3 51 63 50 17
+1 3 12 45 64 38
+1 3 45 20 51 64
+1 3 64 51 17 52
+1 3 38 64 52 7
+
+boundary
+37
+1 1 0 21
+1 1 24 0
+1 1 21 8
+1 1 11 24
+1 1 8 25
+1 1 25 1
+1 1 1 26
+1 1 26 9
+1 1 9 28
+1 1 28 70
+1 1 70 4
+1 1 32 11
+1 1 3 32
+1 1 35 3
+1 1 14 35
+1 1 7 39
+1 1 39 13
+1 1 42 14
+1 1 13 41
+1 1 41 6
+1 1 6 42
+1 1 4 103
+1 1 103 83
+1 1 83 43
+1 1 43 86
+1 1 86 15
+1 1 15 89
+1 1 89 46
+1 1 46 5
+1 1 5 47
+1 1 47 16
+1 1 16 49
+1 1 49 2
+1 1 2 50
+1 1 50 17
+1 1 17 52
+1 1 52 7
+
+vertex_parents
+50
+65 18 27
+66 27 55
+67 30 55
+68 18 30
+69 28 55
+70 4 28
+71 4 29
+72 29 55
+73 10 29
+74 10 30
+75 10 33
+76 29 58
+77 33 58
+78 4 36
+79 36 58
+80 12 36
+81 12 37
+82 37 58
+83 4 43
+84 43 61
+85 36 61
+86 15 43
+87 15 44
+88 44 61
+89 15 46
+90 46 62
+91 44 62
+92 65 67
+93 69 71
+94 4 71
+95 67 73
+96 73 77
+97 71 79
+98 71 97
+99 4 78
+100 78 97
+101 79 81
+102 83 85
+103 4 83
+104 83 102
+105 78 102
+106 86 88
+107 89 91
+108 94 100
+109 94 108
+110 4 94
+111 4 99
+112 99 108
+113 103 105
+114 110 112
+
+coarse_elements
+12
+3 8 11 10 9
+3 13 14 15 16
+3 19 20 17 18
+3 29 32 31 30
+3 35 36 37 38
+3 34 88 39 40
+3 33 89 41 42
+3 43 44 45 46
+3 56 57 58 59
+3 92 62 61 60
+3 63 64 65 66
+3 69 72 71 70
+
+vertices
+115
+
+nodes
+FiniteElementSpace
+FiniteElementCollection: H1_2D_P2
+VDim: 2
+Ordering: 1
+
+-1 -1
+0 -1
+1 1
+-1 0
+0 0
+1 0
+-1 1
+0 1
+-0.5 -1
+0 -0.5
+-0.5 0
+-1 -0.5
+0 0.5
+-0.5 1
+-1 0.5
+0.5 0
+1 0.5
+0.5 1
+-0.5 -0.5
+-0.5 0.5
+0.5 0.5
+-0.75 -1
+-0.5 -0.75
+-0.75 -0.5
+-1 -0.75
+-0.25 -1
+0 -0.75
+-0.25 -0.5
+0 -0.25
+-0.25 0
+-0.5 -0.25
+-0.75 0
+-1 -0.25
+-0.5 0.25
+-0.75 0.5
+-1 0.25
+0 0.25
+-0.25 0.5
+0 0.75
+-0.25 1
+-0.5 0.75
+-0.75 1
+-1 0.75
+0.25 0
+0.5 0.25
+0.25 0.5
+0.75 0
+1 0.25
+0.75 0.5
+1 0.75
+0.75 1
+0.5 0.75
+0.25 1
+-0.75 -0.75
+-0.25 -0.75
+-0.25 -0.25
+-0.75 -0.25
+-0.75 0.25
+-0.25 0.25
+-0.25 0.75
+-0.75 0.75
+0.25 0.25
+0.75 0.25
+0.75 0.75
+0.25 0.75
+-0.375 -0.5
+-0.25 -0.375
+-0.375 -0.25
+-0.5 -0.375
+-0.125 -0.25
+0 -0.125
+-0.125 0
+-0.25 -0.125
+-0.375 0
+-0.5 -0.125
+-0.5 0.125
+-0.25 0.125
+-0.375 0.25
+0 0.125
+-0.125 0.25
+0 0.375
+-0.125 0.5
+-0.25 0.375
+0.125 0
+0.25 0.125
+0.125 0.25
+0.375 0
+0.5 0.125
+0.375 0.25
+0.625 0
+0.75 0.125
+0.625 0.25
+-0.375 -0.375
+-0.125 -0.125
+-0.0625 0
+-0.375 -0.125
+-0.375 0.125
+-0.125 0.125
+-0.125 0.0625
+0 0.0625
+-0.0625 0.125
+-0.125 0.375
+0.125 0.125
+0.0625 0
+0.125 0.0625
+0.0625 0.125
+0.375 0.125
+0.625 0.125
+-0.0625 0.0625
+-0.0625 0.03125
+-0.03125 0
+0 0.03125
+-0.03125 0.0625
+0.0625 0.0625
+-0.03125 0.03125
+-0.875 -1
+-0.75 -0.875
+-0.875 -0.75
+-1 -0.875
+-0.625 -1
+-0.5 -0.875
+-0.625 -0.75
+-0.5 -0.625
+-0.625 -0.5
+-0.75 -0.625
+-0.875 -0.5
+-1 -0.625
+-0.375 -1
+-0.25 -0.875
+-0.375 -0.75
+-0.125 -1
+0 -0.875
+-0.125 -0.75
+0 -0.625
+-0.125 -0.5
+-0.25 -0.625
+-0.375 -0.5
+-0.4375 -0.5
+-0.375 -0.4375
+-0.4375 -0.375
+-0.5 -0.4375
+-0.375 -0.3125
+-0.4375 -0.25
+-0.5 -0.3125
+-0.3125 -0.375
+-0.25 -0.3125
+-0.3125 -0.25
+-0.3125 -0.5
+-0.25 -0.4375
+0 -0.375
+-0.125 -0.25
+-0.25 -0.375
+-0.1875 -0.25
+-0.125 -0.1875
+-0.1875 -0.125
+-0.25 -0.1875
+-0.0625 -0.25
+0 -0.1875
+-0.0625 -0.125
+0 -0.0625
+-0.0625 0
+-0.125 -0.0625
+-0.1875 0
+-0.25 -0.0625
+-0.3125 -0.125
+-0.3125 0
+-0.375 -0.0625
+-0.4375 -0.125
+-0.4375 0
+-0.5 -0.0625
+-0.375 -0.1875
+-0.5 -0.1875
+-0.75 -0.375
+-0.875 -0.25
+-1 -0.375
+-0.5 -0.375
+-0.625 -0.25
+-0.5 -0.125
+-0.625 0
+-0.75 -0.125
+-0.875 0
+-1 -0.125
+-0.75 0.125
+-0.875 0.25
+-1 0.125
+-0.5 0.125
+-0.625 0.25
+-0.5 0.375
+-0.625 0.5
+-0.75 0.375
+-0.875 0.5
+-1 0.375
+-0.375 0.0625
+-0.4375 0.125
+-0.5 0.0625
+-0.375 0.1875
+-0.4375 0.25
+-0.5 0.1875
+-0.3125 0.125
+-0.25 0.1875
+-0.3125 0.25
+-0.25 0.0625
+-0.125 0.0625
+-0.1875 0.125
+-0.09375 0
+-0.0625 0.03125
+-0.09375 0.0625
+-0.125 0.03125
+-0.046875 0
+-0.03125 0.015625
+-0.046875 0.03125
+-0.0625 0.015625
+-0.015625 0
+0 0.015625
+-0.015625 0.03125
+0 0.046875
+-0.015625 0.0625
+-0.03125 0.046875
+-0.046875 0.0625
+-0.0625 0.046875
+-0.03125 0.0625
+0 0.09375
+-0.03125 0.125
+-0.0625 0.09375
+-0.09375 0.125
+-0.125 0.09375
+-0.0625 0.125
+0 0.1875
+-0.0625 0.25
+-0.125 0.1875
+-0.1875 0.25
+-0.125 0.3125
+-0.1875 0.375
+-0.25 0.3125
+0 0.3125
+-0.0625 0.375
+0 0.4375
+-0.0625 0.5
+-0.125 0.4375
+-0.1875 0.5
+-0.25 0.4375
+-0.375 0.25
+-0.25 0.375
+-0.375 0.5
+-0.25 0.625
+-0.375 0.75
+-0.5 0.625
+-0.125 0.5
+0 0.625
+-0.125 0.75
+0 0.875
+-0.125 1
+-0.25 0.875
+-0.375 1
+-0.5 0.875
+-0.75 0.625
+-0.875 0.75
+-1 0.625
+-0.625 0.75
+-0.625 1
+-0.75 0.875
+-0.875 1
+-1 0.875
+0.03125 0
+0.0625 0.03125
+0.03125 0.0625
+0 0.03125
+0.09375 0
+0.125 0.03125
+0.09375 0.0625
+0.125 0.09375
+0.09375 0.125
+0.0625 0.09375
+0.03125 0.125
+0.0625 0.125
+0.125 0.1875
+0.0625 0.25
+0.1875 0.125
+0.25 0.1875
+0.1875 0.25
+0.1875 0
+0.25 0.0625
+0.125 0.0625
+0.3125 0
+0.375 0.0625
+0.3125 0.125
+0.4375 0
+0.5 0.0625
+0.4375 0.125
+0.5 0.1875
+0.4375 0.25
+0.375 0.1875
+0.3125 0.25
+0.375 0.25
+0.5 0.375
+0.375 0.5
+0.25 0.375
+0.125 0.25
+0.125 0.5
+0 0.375
+0.5625 0
+0.625 0.0625
+0.5625 0.125
+0.625 0.1875
+0.5625 0.25
+0.6875 0.125
+0.75 0.1875
+0.6875 0.25
+0.6875 0
+0.75 0.0625
+0.875 0
+1 0.125
+0.875 0.25
+0.75 0.125
+1 0.375
+0.875 0.5
+0.75 0.375
+0.625 0.25
+0.625 0.5
+0.75 0.625
+0.625 0.75
+0.5 0.625
+1 0.625
+0.875 0.75
+1 0.875
+0.875 1
+0.75 0.875
+0.625 1
+0.5 0.875
+0.25 0.625
+0.125 0.75
+0.375 0.75
+0.375 1
+0.25 0.875
+0.125 1
+-0.875 -0.875
+-0.625 -0.875
+-0.625 -0.625
+-0.875 -0.625
+-0.375 -0.875
+-0.125 -0.875
+-0.125 -0.625
+-0.375 -0.625
+-0.4375 -0.4375
+-0.4375 -0.3125
+-0.3125 -0.3125
+-0.3125 -0.4375
+-0.125 -0.375
+-0.1875 -0.1875
+-0.0625 -0.1875
+-0.0625 -0.0625
+-0.1875 -0.0625
+-0.3125 -0.0625
+-0.4375 -0.0625
+-0.4375 -0.1875
+-0.3125 -0.1875
+-0.875 -0.375
+-0.625 -0.375
+-0.625 -0.125
+-0.875 -0.125
+-0.875 0.125
+-0.625 0.125
+-0.625 0.375
+-0.875 0.375
+-0.4375 0.0625
+-0.4375 0.1875
+-0.3125 0.1875
+-0.3125 0.0625
+-0.1875 0.0625
+-0.09375 0.03125
+-0.046875 0.015625
+-0.015625 0.015625
+-0.015625 0.046875
+-0.046875 0.046875
+-0.03125 0.09375
+-0.09375 0.09375
+-0.0625 0.1875
+-0.1875 0.1875
+-0.1875 0.3125
+-0.0625 0.3125
+-0.0625 0.4375
+-0.1875 0.4375
+-0.375 0.375
+-0.375 0.625
+-0.125 0.625
+-0.125 0.875
+-0.375 0.875
+-0.875 0.625
+-0.625 0.625
+-0.625 0.875
+-0.875 0.875
+0.03125 0.03125
+0.09375 0.03125
+0.09375 0.09375
+0.03125 0.09375
+0.0625 0.1875
+0.1875 0.1875
+0.1875 0.0625
+0.3125 0.0625
+0.4375 0.0625
+0.4375 0.1875
+0.3125 0.1875
+0.375 0.375
+0.125 0.375
+0.5625 0.0625
+0.5625 0.1875
+0.6875 0.1875
+0.6875 0.0625
+0.875 0.125
+0.875 0.375
+0.625 0.375
+0.625 0.625
+0.875 0.625
+0.875 0.875
+0.625 0.875
+0.125 0.625
+0.375 0.625
+0.375 0.875
+0.125 0.875

examples/mal_rllib_eval.py

@@ -0,0 +1,188 @@
+
+import commentjson
+from ray.rllib.agents.registry import get_agent_class
+import amr_env
+import gym
+from gym import spaces
+import ray
+import ray.rllib.agents.ppo as ppo
+import tensorflow as tf
+
+import numpy as np
+
+
+import argparse
+import random
+from amr import models
+import pytest
+import os
+import copy
+from math import sqrt, nan, inf, isnan
+from ray.rllib.models import ModelCatalog
+from amr.models.cnn import CNNSmall
+
+# rllib requires you to give the policy an env with the same action
+# and observation spaces as used in training. The rest of it can be
+# "fake" if you provide your own observation data some other way.
+
+# USER INPUT - solution, 20x20 mesh, sine,tanh,steps,steps2, norm-diff reward with random threshold
+#866764
+#checkpoint_folder = "/p/lustre1/mittal3/local_deref/PPO/PPO_LocalAMR-v0_5b99d_00000_0_2021-06-07_12-12-52/"
+#checkpoint_number = 900
+
+# USER INPUT - solution, 20x20 mesh, sine,tanh,steps,steps2, binary reward with random threshold
+#866762
+# checkpoint_folder = "/p/lustre1/mittal3/local_deref/PPO/PPO_LocalAMR-v0_8d234_00000_0_2021-06-07_12-07-06/"
+# checkpoint_number = 900
+
+# solution, 20x20, steps2,sine,tanh,bumps,  binary with random
+# 880258 - fixed threshold 1.e-5
+checkpoint_folder = "/p/lustre1/mittal3/local_deref/PPO/PPO_LocalAMR-v0_cc8cb_00000_0_2021-06-16_13-38-25"
+checkpoint_number = 1400
+
+# 880259 - random threshold [1.e-2, 1.e-6]
+checkpoint_folder = "/p/lustre1/mittal3/local_deref/PPO/PPO_LocalAMR-v0_e1561_00000_0_2021-06-16_13-31-50"
+checkpoint_number = 900
+# with 1.e-3 - the second and third refinements are really good. still more than needed in first
+# with 1.e-2 - picks the right amount of elements.
+
+#880328 - fixed threshold 1.e-2
+#checkpoint_folder = "/p/lustre1/mittal3/local_deref/PPO/PPO_LocalAMR-v0_3a8a0_00000_0_2021-06-16_16-18-58"
+#checkpoint_number = 300
+#refines 15 elements at first iteration.. not good with 300.
+
+# solution, 10x10, steps2,sine,tanh,bumps,  binary with random [1.e-2, 1.e-6]
+# slurm-880260.out
+#checkpoint_folder = "/p/lustre1/mittal3/local_deref/PPO/PPO_LocalAMR-v0_cc8cb_00000_0_2021-06-16_13-38-25/"
+#checkpoint_number = 1500
+#1.e-3 -> right region but first iteration has too many
+#1.e-2 is almost perfect
+
+
+
+
+#this has the policy without error threshold
+#checkpoint_folder = "/p/lustre1/mittal3/local_deref/PPO/PPO_LocalAMR-v0_c6ca5_00000_0_2021-05-24_16-07-13/"
+
+# END OF USER INPUT
+
+# Read info from json file used for training.
+full_checkpoint_path = checkpoint_folder + '/checkpoint_' + str(checkpoint_number)+ '/checkpoint-' + str(checkpoint_number)
+path_env_config_file = checkpoint_folder + '/params.json'
+with open(path_env_config_file) as json_file:
+    env_trainer_config = commentjson.load(json_file)
+
+trainer_config = env_trainer_config
+trainer_config['env_config']['mesh_params']['nx'] = 1
+trainer_config['env_config']['mesh_params']['ny'] = 1
+local_sample = env_trainer_config['env_config']['local_sample']
+local_context = env_trainer_config['env_config']['local_context']
+reward_params = env_trainer_config['env_config']['reward_function_params']
+
+#set some default params
+observe_error = False
+observe_values = True
+observe_grads = False
+
+#get observing quantities
+observe_values = env_trainer_config['env_config']['observe_values']
+observe_depth = env_trainer_config['env_config']['observe_depth']
+observe_jacobian = env_trainer_config['env_config']['observe_jacobian']
+observe_ar = env_trainer_config['env_config']['observe_ar']
+observe_grads = env_trainer_config['env_config']['observe_grads']
+normalization = env_trainer_config['env_config']['normalization']
+
+#get reward_params
+reward_params = env_trainer_config['env_config']['reward_function_params']
+reward_params_name = reward_params['name']
+if reward_params_name == "random_penalized_norm_diff":
+    observe_error = True
+if reward_params_name == "random_binary":
+    observe_error = True
+
+class DummyEnv(gym.Env):
+
+    def __init__(self, config): # the config param is required by rllib
+
+        # image size is 42x42 (a size which uses CNN by default in rllib)
+        self.obsx = local_sample+2*local_context
+        self.obsy = self.obsx
+
+        # Either do nothing (0) or refine (1)
+        self.action_space = spaces.Discrete(2)
+
+        n_channels = 1
+        n_channels = observe_values + observe_depth + observe_grads
+
+        low = -np.inf
+        
+        high = np.inf
+        self.observation_space = spaces.Dict({
+                "scalar_info": spaces.Box(low=low, high=high, shape=(1 + observe_jacobian + observe_error, ), dtype=np.float32),
+                "obs_data":    spaces.Box(low=low, high=high,
+                shape=(self.obsx,
+                       self.obsy,
+                       n_channels), dtype=np.float32)
+                })
+
+        self.state = None
+        
+    def step(self, action):
+        pass
+    def reset(self):
+        pass
+    def render(self):
+        pass
+
+class Evaluator():
+
+    def __init__(self):
+
+        print("starting ray...")
+        ray.shutdown()
+        ray.init()
+        print("ray up...")
+
+        ModelCatalog.register_custom_model("cnn_small", CNNSmall)
+        trainer = ppo.PPOTrainer(config=trainer_config,env=DummyEnv)
+        trainer.restore(full_checkpoint_path)
+
+        trainer_config['evaluation_num_workers'] = 1
+        trainer_config['evaluation_interval'] = 0
+        trainer_config['num_workers'] = 0
+        trainer_config['num_envs_per_worker'] = 1
+
+        self.agent = ppo.PPOTrainer(config=trainer_config,env=DummyEnv)
+        self.agent.restore(full_checkpoint_path)
+        
+        self.env = DummyEnv({})
+
+    def eval(self,obso,scalar):
+        if observe_values and normalization:
+            obso -= np.mean(obso)
+        
+        obs = {
+            "obs_data" : obso,
+            "scalar_info" : scalar
+            }
+
+        pick = self.agent.compute_action(obs, explore=False)
+        return pick
+
+    def get_local_sample(self):
+        return local_sample
+
+    def get_local_context(self):
+        return local_context
+
+    def get_observe_error(self):
+        return observe_error
+
+    def get_observe_jacobian(self):
+        return observe_jacobian
+
+    def get_observe_values(self):
+        return observe_values
+
+    def get_observe_gradient(self):
+        return observe_grads

examples/multi_agent_local_refiner.hpp

@@ -0,0 +1,419 @@
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace mfem;
+using namespace std;
+
+class MAL_DRLRefiner : public MeshOperator
+{
+protected:
+   GridFunction& u;
+
+   Array<Refinement> marked_elements;
+   long current_sequence;
+
+   int imgsz;
+   int local_sample;
+   int local_context;
+   bool observe_jacobian;
+   bool observe_error;
+   bool observe_values;
+   bool observe_gradient;
+
+   int nc_limit;
+
+   PyObject* eval_method;
+   PyObject* get_local_sample_method;
+   PyObject* get_local_context_method;
+   PyObject* get_observe_error_method;
+   PyObject* get_observe_jacobian_method;
+   PyObject* get_observe_values_method;
+   PyObject* get_observe_gradient_method;
+
+   FindPointsGSLIB *gslib;
+
+   /** @brief Apply the operator to the mesh.
+       @return STOP if a stopping criterion is satisfied or no elements were
+       marked for refinement; REFINED + CONTINUE otherwise. */
+   virtual int ApplyImpl(Mesh &mesh);
+
+public:
+
+   /// Construct a MAL_DRLRefiner that will operate on u.
+   MAL_DRLRefiner(GridFunction &u);
+
+   // default destructor (virtual)
+
+   /** @brief Set the maximum ratio of refinement levels of adjacent elements
+       (0 = unlimited). */
+   void SetNCLimit(int nc_limit)
+   {
+      MFEM_ASSERT(nc_limit >= 0, "Invalid NC limit");
+      this->nc_limit = nc_limit;
+   }
+
+   virtual void Reset();
+};
+
+MAL_DRLRefiner::MAL_DRLRefiner(GridFunction& u_) : u(u_)
+{
+   int ret = _import_array();
+   if (ret < 0) {
+      printf("problem with import_array\n");
+   }
+
+   PyRun_SimpleString("import sys");
+   PyRun_SimpleString("sys.path.append('.')");
+
+   // This is a workaround for something in tensorflow that dies without it.
+   PyRun_SimpleString("if not hasattr(sys, 'argv'):\n"
+                      "  sys.argv  = ['']");
+
+   PyObject* eval_mod = PyImport_ImportModule("mal_rllib_eval");
+   if (eval_mod == 0) {
+      PyErr_Print();
+      exit(1);
+   }
+
+   PyObject* eval_class = PyObject_GetAttrString(eval_mod, "Evaluator");
+   if (eval_class == 0) {
+      PyErr_Print();
+      exit(1);
+   }
+   Py_DECREF(eval_mod);
+
+   PyObject* args = Py_BuildValue("()");
+   if (args == 0) {
+      PyErr_Print();
+      exit(1);
+   }
+
+   PyObject* eval_obj = PyEval_CallObject(eval_class, args);
+   if (eval_obj == NULL) {
+      PyErr_Print();
+      exit(1);
+   }
+
+   eval_method = PyObject_GetAttrString(eval_obj, "eval");
+   if (eval_method == 0) {
+      PyErr_Print();
+      exit(1);
+   }
+
+   get_local_sample_method = PyObject_GetAttrString(eval_obj, "get_local_sample");
+   get_local_context_method = PyObject_GetAttrString(eval_obj, "get_local_context");
+   get_observe_error_method = PyObject_GetAttrString(eval_obj, "get_observe_error");
+   get_observe_jacobian_method = PyObject_GetAttrString(eval_obj, "get_observe_jacobian");
+   get_observe_values_method = PyObject_GetAttrString(eval_obj, "get_observe_values");
+   get_observe_gradient_method = PyObject_GetAttrString(eval_obj, "get_observe_gradient");
+
+   PyObject* local_sample_p = PyObject_CallFunctionObjArgs(
+                              get_local_sample_method, nullptr);
+   PyObject* local_context_p = PyObject_CallFunctionObjArgs(
+                              get_local_context_method, nullptr);
+   PyObject* observe_jacobian_p = PyObject_CallFunctionObjArgs(
+                              get_observe_jacobian_method, nullptr);
+   PyObject* observe_error_p = PyObject_CallFunctionObjArgs(
+                              get_observe_error_method, nullptr);
+   PyObject* observe_values_p = PyObject_CallFunctionObjArgs(
+                                get_observe_values_method, nullptr);
+   PyObject* observe_gradient_p = PyObject_CallFunctionObjArgs(
+                                  get_observe_gradient_method, nullptr);
+
+   PyArg_Parse(local_sample_p, "i", &local_sample);
+   PyArg_Parse(local_context_p, "i", &local_context);
+   int observe_jacobian_i, observe_error_i, observe_values_i, observe_gradient_i;
+   PyArg_Parse(observe_jacobian_p, "i", &observe_jacobian_i);
+   PyArg_Parse(observe_error_p, "i", &observe_error_i);
+   PyArg_Parse(observe_values_p, "i", &observe_values_i);
+   PyArg_Parse(observe_gradient_p, "i", &observe_gradient_i);
+
+   observe_jacobian = bool(observe_jacobian_i);
+   observe_error = bool(observe_error_i);
+   observe_values = bool(observe_values_i);
+   observe_gradient = bool(observe_gradient_i);
+
+#ifdef MFEM_USE_GSLIB
+   // setup gslib
+   gslib = new FindPointsGSLIB();
+   std::cout << " initialize findpts\n";
+#endif
+}
+
+int MAL_DRLRefiner::ApplyImpl(Mesh &mesh)
+{
+#ifdef MFEM_USE_GSLIB
+   // setup gslib
+   gslib->FreeData();
+   gslib->Setup(mesh);
+#endif
+   double u_min = u.Min();
+   double u_max = u.Max();
+   // u -= u_min;
+   // u /= (u_max-u_min); 
+
+   marked_elements.SetSize(0);
+   imgsz = local_sample + 2 * local_context;
+
+   GridFunction ugrad(u.FESpace());
+   GridFunction ugradmag(u.FESpace());
+
+   if (observe_gradient) {
+      const int s = ugrad.Size();
+
+      for (int d = 0; d < 2; d++)
+      {
+         u.GetDerivative(1, d, ugrad);
+         for (int i = 0; i < s; i++)
+         {
+            ugradmag(i) += pow(ugrad(i), 2.0);
+         }
+      }
+      for (int i = 0; i < s; i++)
+      {
+         ugradmag(i) = sqrt(ugradmag(i) + 1e-12);
+      }
+   }
+      
+   for (int k = 0; k < mesh.GetNE(); k++) {
+
+      // get scalar info
+      const int scalar_size = 1 + observe_error + observe_jacobian;
+      double* scalar_obs1 = new double[scalar_size];
+      scalar_obs1[0] = 1;
+      bool boundary = false;
+      Array<int> fcs, cor;
+      int e1, e2, inf1, inf2, ncf;
+      mesh.GetElementEdges(k, fcs, cor);
+      for (int f = 0; f < fcs.Size() && boundary == false; f++) {
+         mesh.GetFaceElements(fcs[f], &e1, &e2);
+         mesh.GetFaceInfos(fcs[f], &inf1, &inf2, &ncf);
+         if (e2 < 0 && inf2 < 0 && ncf == -1) {
+            boundary = true;
+         }
+      }
+      if (!boundary) {
+         scalar_obs1[0] = 2;
+      }
+
+
+      double error_threshold = 1.0e-2;
+      if (observe_jacobian && observe_error) {
+         scalar_obs1[1] = mesh.GetElementVolume(k);
+         scalar_obs1[2] = error_threshold;
+      }
+      else if (observe_jacobian) {
+         scalar_obs1[1] = mesh.GetElementVolume(k);
+      }
+      else if (observe_error) {
+         scalar_obs1[1] = error_threshold;
+      }
+
+      // assemble matrix of sample points, ref_space -> phys_space
+      ElementTransformation* trk = mesh.GetElementTransformation(k);
+      IntegrationPoint ipk;
+      Vector xk(2);
+      DenseMatrix m(2,imgsz*imgsz);
+      int c = 0;
+      double r_init = 0.001,
+             r_final = 1.0-r_init;
+      for (int j = 0; j < imgsz; ++j) {
+         for (int i = 0; i < imgsz; ++i) {
+            if (!boundary) {
+               ipk.y = (j - local_context + 0.5)/local_sample;
+               ipk.x = (i - local_context + 0.5)/local_sample;
+            }
+            else {
+               ipk.y = (r_init + j*(r_final-r_init))/(imgsz-1);
+               ipk.x = (r_init + i*(r_final-r_init))/(imgsz-1);
+            }
+
+            trk->Transform(ipk, xk);
+            m.SetCol(c++,xk);
+         }
+      } 
+
+      // phys_space -> elements, ips
+      Array<int> elems(imgsz*imgsz);
+      Array<IntegrationPoint> ips(imgsz*imgsz);
+      int n;
+      bool complete;
+      Vector ui;
+      double* obs1 = new double[imgsz*imgsz];
+#ifdef MFEM_USE_GSLIB
+      DenseMatrix mt(imgsz*imgsz, 2);
+      mt.Transpose(m);
+      Vector xyz(mt.GetData(), imgsz*imgsz*2);
+      ui.SetDataAndSize(obs1, imgsz*imgsz);
+      gslib->Interpolate(xyz, u, ui);
+#else
+      n = mesh.FindPoints(m, elems, ips, false);
+
+      // Build observation from GridFunction using elements, ips
+      n = 0;
+      complete = true;
+      for (int j = 0; j < imgsz; ++j) {
+         for (int i = 0; i < imgsz; ++i) {
+            int el = elems[n];
+            if (el == -1) {
+               obs1[i*imgsz+j] = 0.0;
+               complete = false;
+            }
+            else {
+               IntegrationPoint& ip = ips[n];
+               if (observe_gradient) {
+                  //std::cout << i << " " << j << " k10getgradmag\n";
+                  obs1[i*imgsz+j] = ugradmag.GetValue(el, ip);
+                  obs1[i*imgsz+j] = u.GetValue(el, ip);
+               }
+               else {
+                  obs1[i*imgsz+j] = u.GetValue(el, ip);
+               }
+            }
+            n++;
+         }
+      }
+#endif
+
+      // TODO: More efficient way to do the below mirroring
+      // invert i
+      double* obs2 = new double[imgsz*imgsz];
+      n = 0;
+      complete = true;
+      for (int j = 0; j < imgsz; ++j) {
+         for (int i = imgsz-1; i >= 0; --i) {
+            // int el = elems[n];
+            // if (el == -1) {
+            //    obs2[i*imgsz+j] = 0.0;
+            //    complete = false;
+            // }
+            // else {
+            //    IntegrationPoint& ip = ips[n];
+            //    obs2[i*imgsz+j] = u.GetValue(el, ip);
+            // }
+            obs2[i*imgsz+j] = obs1[n];
+            n++;
+         }
+      }
+
+      // invert j
+      double* obs3 = new double[imgsz*imgsz];
+      n = 0;
+      complete = true;
+      for (int j = imgsz-1; j >= 0; --j) {
+         for (int i = 0; i < imgsz; ++i) {
+            // int el = elems[n];
+            // if (el == -1) {
+            //    obs3[i*imgsz+j] = 0.0;
+            //    complete = false;
+            // }
+            // else {
+            //    IntegrationPoint& ip = ips[n];
+            //    obs3[i*imgsz+j] = u.GetValue(el, ip);
+            // }
+            obs3[i*imgsz+j] = obs1[n];
+            n++;
+         }
+      }
+
+      // invert i and j
+      double* obs4 = new double[imgsz*imgsz];
+      n = 0;
+      complete = true;
+      for (int j = imgsz-1; j >= 0; --j) {
+         for (int i = imgsz-1; i >= 0; --i) {
+            // int el = elems[n];
+            // if (el == -1) {
+            //    obs4[i*imgsz+j] = 0.0;
+            //    complete = false;
+            // }
+            // else {
+            //    IntegrationPoint& ip = ips[n];
+            //    obs4[i*imgsz+j] = u.GetValue(el, ip);
+            // }
+            obs4[i*imgsz+j] = obs1[n];
+            n++;
+         }
+      }
+
+
+      // apply policy: state -> action
+      bool refine = false;
+      if (complete) {
+
+         // convert to numpy array
+         npy_intp dims[3];
+         dims[0] = imgsz;
+         dims[1] = imgsz;
+         dims[2] = 1;
+
+         npy_intp scalar_dims[1];
+         scalar_dims[0] = scalar_size;
+
+         PyObject *pArray1 = PyArray_SimpleNewFromData(
+            3, dims, NPY_DOUBLE, reinterpret_cast<void*>(obs1));
+         PyObject *pArray2 = PyArray_SimpleNewFromData(
+            3, dims, NPY_DOUBLE, reinterpret_cast<void*>(obs2));
+         PyObject *pArray3 = PyArray_SimpleNewFromData(
+            3, dims, NPY_DOUBLE, reinterpret_cast<void*>(obs3));
+         PyObject *pArray4 = PyArray_SimpleNewFromData(
+            3, dims, NPY_DOUBLE, reinterpret_cast<void*>(obs4));   
+         PyObject *sArray1 = PyArray_SimpleNewFromData(
+            1, scalar_dims, NPY_DOUBLE, reinterpret_cast<void*>(scalar_obs1));
+         if (pArray1 == NULL) printf("pArray1 NULL!\n");
+
+         PyObject* action1 = PyObject_CallFunctionObjArgs(
+                     eval_method, pArray1, sArray1, nullptr);
+         PyObject* action2 = PyObject_CallFunctionObjArgs(
+                     eval_method, pArray2, sArray1, nullptr);
+         PyObject* action3 = PyObject_CallFunctionObjArgs(
+                     eval_method, pArray3, sArray1, nullptr);
+         PyObject* action4 = PyObject_CallFunctionObjArgs(
+                     eval_method, pArray4, sArray1, nullptr);
+
+         if (action1 == 0 || action2 == 0) {
+            PyErr_Print();
+            exit(1);
+         }
+
+         // parse integer return value
+         int action1_val;
+         PyArg_Parse(action1, "i", &action1_val);
+         int action2_val;
+         PyArg_Parse(action2, "i", &action2_val);
+         int action3_val;
+         PyArg_Parse(action3, "i", &action3_val);
+         int action4_val;
+         PyArg_Parse(action4, "i", &action4_val);
+         refine =
+            bool(action1_val) ||
+            bool(action2_val) ||
+            bool(action3_val) ||
+            bool(action4_val);
+      }
+
+      delete scalar_obs1;
+      delete obs1;
+      delete obs2;
+      delete obs3;
+      delete obs4;
+
+      if (refine) {
+         marked_elements.Append(Refinement(k));
+      }
+   }
+
+   long int num_marked_elements = mesh.ReduceInt(marked_elements.Size());
+   printf("marked %d elements\n",num_marked_elements);
+   if (num_marked_elements == 0) { return STOP; }
+
+   bool nonconforming = true;
+   mesh.GeneralRefinement(marked_elements, nonconforming, nc_limit);
+   return CONTINUE + REFINED;
+}
+
+void MAL_DRLRefiner::Reset()
+{
+   current_sequence = -1;
+}

examples/rllib_eval.py

@@ -0,0 +1,64 @@
+
+
+import gym
+from gym import spaces
+import ray
+import ray.rllib.agents.ppo as ppo
+import tensorflow as tf
+
+import numpy as np
+
+# rllib requires you to give the policy an env with the same action
+# and observation spaces as used in training. The rest of it can be
+# "fake" if you provide your own observation data some other way.
+
+class DummyEnv(gym.Env):
+
+    def __init__(self, config): # the config param is required by rllib
+
+        # image size is 42x42 (a size which uses CNN by default in rllib)
+        self.obsx = 42
+        self.obsy = 42
+
+        # Either do nothing (0) or refine (1)
+        self.action_space = spaces.Discrete(2)
+        self.observation_space = spaces.Box(-1.0, 2.0, shape=(self.obsx,self.obsy,1))
+
+        self.state = None
+        
+    def step(self, action):
+        pass
+    def reset(self):
+        pass
+    def render(self):
+        pass
+
+class Evaluator():
+
+    def __init__(self):
+
+        print("starting ray...")
+        ray.shutdown()
+        ray.init()
+        print("ray up...")
+
+        config = ppo.DEFAULT_CONFIG.copy()
+        config["log_level"] = "WARN"
+
+        # Create agent from checkpoint
+        self.agent = ppo.PPOTrainer(config,env=DummyEnv)
+        self.agent.restore("DRLRefinePolicy/checkpoint_210/checkpoint-210")
+
+        self.env = DummyEnv({})
+
+    def eval(self,obs):
+        pick = self.agent.compute_action(obs, explore=False)
+        return pick
+
+#evaluator = Evaluator()
+#obs = np.ones((42,42,1))
+#ref = evaluator.eval(obs)
+#evaluator.show_logits(obs)
+# evaluator.eval(np.random.rand(8))
+# evaluator.eval(np.random.rand(8))
+

makefile

@@ -234,6 +234,9 @@ ifeq ($(MFEM_USE_CUDA),YES)
    endif
 endif
 
+MFEM_CXX += $(shell /usr/bin/python3-config --includes)
+#ALL_LIBS += $(shell /usr/bin/python3-config --ldflags)
+
 # HIP configuration
 ifeq ($(MFEM_USE_HIP),YES)
    ifeq ($(MFEM_USE_MPI),YES)
@@ -311,6 +314,8 @@ $(foreach dep,$(MFEM_DEPENDENCIES),$(eval $(call mfem_add_dependency,$(dep))))
 $(foreach dep,$(MFEM_LEGACY_DEPENDENCIES),$(eval $(call \
    mfem_add_legacy_dependency,$(dep))))
 
+ALL_LIBS += $(shell /usr/bin/python3-config --ldflags)
+
 # Timer option
 ifeq ($(MFEM_TIMER_TYPE),2)
    ALL_LIBS += $(POSIX_CLOCKS_LIB)

mesh/mesh_operators.cpp

@@ -145,6 +145,242 @@ void ThresholdRefiner::Reset()
    // marked_elements.SetSize(0); // not necessary
 }
 
+DRLRefiner::DRLRefiner(GridFunction& u_) : u(u_)
+{
+   int ret = _import_array();
+   if (ret < 0) {
+      printf("problem with import_array\n");
+   }
+
+   obs_x = 42;
+   obs_y = 42;
+
+   PyRun_SimpleString("import sys");
+   PyRun_SimpleString("sys.path.append('.')");
+
+   // This is a workaround for something in tensorflow that dies without it.
+   PyRun_SimpleString("if not hasattr(sys, 'argv'):\n"
+                      "  sys.argv  = ['']");
+
+   PyObject* eval_mod = PyImport_ImportModule("rllib_eval");
+   if (eval_mod == 0) {
+      PyErr_Print();
+      exit(1);
+   }
+
+   PyObject* eval_class = PyObject_GetAttrString(eval_mod, "Evaluator");
+   if (eval_class == 0) {
+      PyErr_Print();
+      exit(1);
+   }
+   Py_DECREF(eval_mod);
+
+   PyObject* args = Py_BuildValue("()");
+   if (args == 0) {
+      PyErr_Print();
+      exit(1);
+   }
+
+   PyObject* eval_obj = PyEval_CallObject(eval_class, args);
+   if (eval_obj == NULL) {
+      PyErr_Print();
+      exit(1);
+   }
+
+   eval_method = PyObject_GetAttrString(eval_obj, "eval");
+   if (eval_method == 0) {
+      PyErr_Print();
+      exit(1);
+   }
+}
+
+int DRLRefiner::ApplyImpl(Mesh &mesh)
+{
+   marked_elements.SetSize(0);
+
+   const int imgsz = 42;
+
+   double ref_w = 2.0;
+   assert(obs_x == obs_y);
+   int ref_n = obs_x;
+   double ref_lo = 0.0 -ref_w;
+   double ref_hi = 1.0 +ref_w;
+   double ref_dx = (ref_hi -ref_lo)/ref_n;
+
+   for (int k = 0; k < mesh.GetNE(); k++) {
+
+      // assemble matrix of sample points, ref_space -> phys_space
+      ElementTransformation* trk = mesh.GetElementTransformation(k);
+      IntegrationPoint ipk;
+      Vector xk(2);
+      DenseMatrix m(2,imgsz*imgsz);
+      int c = 0;
+      for (int j = 0; j < obs_y; ++j) {
+         ipk.y = ref_lo +(j+0.5)*ref_dx;
+         for (int i = 0; i < obs_x; ++i) {
+            ipk.x = ref_lo +(i+0.5)*ref_dx;
+            trk->Transform(ipk, xk);
+            m.SetCol(c++,xk);
+         }
+      }
+
+      // phys_space -> elements, ips
+      Array<int> elems;
+      Array<IntegrationPoint> ips;
+      int n = mesh.FindPoints(m, elems, ips, false);
+
+      // Build observation from GridFunction using elements, ips
+      double* obs1 = new double[imgsz*imgsz];
+      n = 0;
+      bool complete = true;
+      for (int j = 0; j < obs_y; ++j) {
+         for (int i = 0; i < obs_x; ++i) {
+            int el = elems[n];
+            if (el == -1) {
+               obs1[i*imgsz+j] = 0.0;
+               complete = false;
+            }
+            else {
+               IntegrationPoint& ip = ips[n];
+               obs1[i*imgsz+j] = u.GetValue(el, ip);
+            }
+            n++;
+         }
+      }
+
+      // TODO: More efficient way to do the below mirroring
+
+      // invert i
+      double* obs2 = new double[imgsz*imgsz];
+      n = 0;
+      complete = true;
+      for (int j = 0; j < obs_y; ++j) {
+         for (int i = obs_x-1; i >= 0; --i) {
+            int el = elems[n];
+            if (el == -1) {
+               obs2[i*imgsz+j] = 0.0;
+               complete = false;
+            }
+            else {
+               IntegrationPoint& ip = ips[n];
+               obs2[i*imgsz+j] = u.GetValue(el, ip);
+            }
+            n++;
+         }
+      }
+
+      // invert j
+      double* obs3 = new double[imgsz*imgsz];
+      n = 0;
+      complete = true;
+      for (int j = obs_y-1; j >= 0; --j) {
+         for (int i = 0; i < obs_x; ++i) {
+            int el = elems[n];
+            if (el == -1) {
+               obs3[i*imgsz+j] = 0.0;
+               complete = false;
+            }
+            else {
+               IntegrationPoint& ip = ips[n];
+               obs3[i*imgsz+j] = u.GetValue(el, ip);
+            }
+            n++;
+         }
+      }
+
+      // invert i and j
+      double* obs4 = new double[imgsz*imgsz];
+      n = 0;
+      complete = true;
+      for (int j = obs_y-1; j >= 0; --j) {
+         for (int i = obs_x-1; i >= 0; --i) {
+            int el = elems[n];
+            if (el == -1) {
+               obs4[i*imgsz+j] = 0.0;
+               complete = false;
+            }
+            else {
+               IntegrationPoint& ip = ips[n];
+               obs4[i*imgsz+j] = u.GetValue(el, ip);
+            }
+            n++;
+         }
+      }
+
+      // apply policy: state -> action
+      bool refine = false;
+      if (complete) {
+
+         // convert to numpy array
+         npy_intp dims[3];
+         dims[0] = imgsz;
+         dims[1] = imgsz;
+         dims[2] = 1;
+
+         PyObject *pArray1 = PyArray_SimpleNewFromData(
+            3, dims, NPY_DOUBLE, reinterpret_cast<void*>(obs1));
+         if (pArray1 == NULL) printf("pArray1 NULL!\n");
+
+         PyObject *pArray2 = PyArray_SimpleNewFromData(
+            3, dims, NPY_DOUBLE, reinterpret_cast<void*>(obs2));
+         if (pArray2 == NULL) printf("pArray2 NULL!\n");
+
+         PyObject *pArray3 = PyArray_SimpleNewFromData(
+            3, dims, NPY_DOUBLE, reinterpret_cast<void*>(obs3));
+         if (pArray3 == NULL) printf("pArray3 NULL!\n");
+
+         PyObject *pArray4 = PyArray_SimpleNewFromData(
+            3, dims, NPY_DOUBLE, reinterpret_cast<void*>(obs4));
+         if (pArray4 == NULL) printf("pArray4 NULL!\n");
+
+         PyObject* action1 = PyObject_CallFunctionObjArgs(
+            eval_method, pArray1, nullptr);
+         PyObject* action2 = PyObject_CallFunctionObjArgs(
+            eval_method, pArray2, nullptr);
+         PyObject* action3 = PyObject_CallFunctionObjArgs(
+            eval_method, pArray3, nullptr);
+         PyObject* action4 = PyObject_CallFunctionObjArgs(
+            eval_method, pArray4, nullptr);
+
+         if (action1 == 0 || action2 == 0) {
+            PyErr_Print();
+            exit(1);
+         }
+
+         // parse integer return value
+         int action1_val;
+         PyArg_Parse(action1, "i", &action1_val);
+         int action2_val;
+         PyArg_Parse(action2, "i", &action2_val);
+         int action3_val;
+         PyArg_Parse(action3, "i", &action3_val);
+         int action4_val;
+         PyArg_Parse(action4, "i", &action4_val);
+         refine =
+            bool(action1_val) ||
+            bool(action2_val) ||
+            bool(action3_val) ||
+            bool(action4_val);
+      }
+
+      if (refine) {
+         marked_elements.Append(Refinement(k));
+      }
+   }
+
+   long int num_marked_elements = mesh.ReduceInt(marked_elements.Size());
+   printf("marked %d elements\n",num_marked_elements);
+   if (num_marked_elements == 0) { return STOP; }
+
+   bool nonconforming = true;
+   mesh.GeneralRefinement(marked_elements, nonconforming, nc_limit);
+   return CONTINUE + REFINED;
+}
+
+void DRLRefiner::Reset()
+{
+   current_sequence = -1;
+}
 
 int ThresholdDerefiner::ApplyImpl(Mesh &mesh)
 {

mesh/mesh_operators.hpp

@@ -19,6 +19,12 @@
 
 #include <limits>
 
+#define MFEM_USE_RLLIB
+#ifdef MFEM_USE_RLLIB
+#include <Python.h>
+#include "numpy/arrayobject.h"
+#endif
+
 namespace mfem
 {
 
@@ -254,6 +260,48 @@ public:
    virtual void Reset();
 };
 
+/** @brief Mesh refinement operator using a DRL policy from rllib.
+
+*/
+
+class DRLRefiner : public MeshOperator
+{
+protected:
+   GridFunction& u;
+
+   Array<Refinement> marked_elements;
+   long current_sequence;
+
+   int obs_x;
+   int obs_y;
+
+   int nc_limit;
+
+   PyObject* eval_method;
+
+   /** @brief Apply the operator to the mesh.
+       @return STOP if a stopping criterion is satisfied or no elements were
+       marked for refinement; REFINED + CONTINUE otherwise. */
+   virtual int ApplyImpl(Mesh &mesh);
+
+public:
+
+   /// Construct a DRLRefiner that will operate on u.
+   DRLRefiner(GridFunction &u);
+
+   // default destructor (virtual)
+
+   /** @brief Set the maximum ratio of refinement levels of adjacent elements
+       (0 = unlimited). */
+   void SetNCLimit(int nc_limit)
+   {
+      MFEM_ASSERT(nc_limit >= 0, "Invalid NC limit");
+      this->nc_limit = nc_limit;
+   }
+
+   virtual void Reset();
+};
+
 // TODO: BulkRefiner to refine a portion of the global error