Removed blobs, added .gitignore for this directory.

examples/makefile

@@ -22,7 +22,7 @@ MFEM_LIB_FILE = mfem_is_not_built
 -include $(CONFIG_MK)
 
 SEQ_EXAMPLES = ex1 ex2 ex3 ex4 ex5 ex6 ex7 ex8 ex9 ex10 ex14 ex15 ex16 ex17\
-  ex18 ex19 ex20 ex21
+  ex18 ex19 ex20 ex21 drl4amr
 PAR_EXAMPLES = ex1p ex2p ex3p ex4p ex5p ex6p ex7p ex8p ex9p ex10p ex11p ex12p\
  ex13p ex14p ex15p ex16p ex17p ex18p ex19p ex20p ex21p
 
@@ -71,6 +71,10 @@ ifeq ($(MFEM_USE_MPI),YES)
 ex18p: $(SRC)ex18.hpp
 endif
 
+drl4amr:
+	python drl4amr.py build
+	g++ -pthread -shared -Wl,-z,relro build/temp.linux-x86_64-2.7/drl4amr.o -L/usr/lib64 -lpython2.7 -o build/lib.linux-x86_64-2.7/drl4amr.so -L.. -lmfem
+
 MFEM_TESTS = EXAMPLES
 include $(MFEM_TEST_MK)
 test: $(SUBDIRS_TEST)

miniapps/drl4amr/.gitignore

@@ -0,0 +1,3 @@
+drl4amr
+*.so 
+image.jpg

miniapps/drl4amr/drl4amr.cpp

@@ -0,0 +1,392 @@
+//                                MFEM DRL4AMR
+
+#define protected public // FIXME!
+
+#include "mfem.hpp"
+#include "drl4amr.hpp"
+
+#include "linalg/dtensor.hpp"
+
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+#define dbg(...) \
+   { printf("\n\033[33m"); printf(__VA_ARGS__); printf("\033[m"); fflush(0); }
+
+const char *device_config = "cpu";
+
+static int discs;
+static double theta;
+static Array<double> offsets;
+constexpr int nb_discs_max = 6;
+constexpr double sharpness = 100.0;
+
+const bool visualization = true;
+const char *vishost = "localhost";
+const int visport = 19916;
+const int visw = 480;
+const int vish = 480;
+
+// *****************************************************************************
+double x0(const Vector &x)
+{
+   double result = 0.0;
+   const double t = x[0] + tan(theta)*x[1];
+   for (double o : offsets)
+   {
+      result += 1.0 + tanh(sharpness*(o - t));
+   }
+   return result / (discs << 1); // should be in [0,1]
+}
+
+// *****************************************************************************
+Drl4Amr::Drl4Amr(int order):
+   mesh(nx, ny, elem_type, true, sx, sy, false),
+   order(order),
+   dim(mesh.Dimension()),
+   device(device_config),
+   fec(order, dim/*, BasisType::Positive*/),
+   fespace(&mesh, &fec),
+   one(1.0),
+   zero(0.0),
+   integ(new DiffusionIntegrator(one)),
+   xcoeff(x0),
+   x(&fespace),
+   iteration(0),
+   flux_fespace(&mesh, &fec),
+   estimator(*integ, x, flux_fespace),
+   refiner(estimator)
+{
+   //dbg("Drl4Amr order:%d",o);
+   device.Print();
+
+   mesh.EnsureNCMesh();
+   mesh.PrintCharacteristics();
+
+   fespace.Update();
+
+   // Connect to GLVis.
+   if (visualization)
+   {
+      vis[0].open(vishost, visport);
+//      vis[1].open(vishost, visport);
+//      vis[2].open(vishost, visport);
+   }
+
+   // Initialize theta, offsets and x from x0_coeff
+   srand48(4);//time(NULL));
+   theta = M_PI*drand48()/2.0;
+   discs = static_cast<int>(1 + nb_discs_max*drand48());
+   offsets.SetSize(discs);
+   for (int i=0; i < discs; i++)
+   {
+      offsets[i] = drand48();
+   }
+   offsets.Sort();
+   printf("\ntheta = %f, discontinuities:%d", theta, discs);
+   for (double offset: offsets)
+   {
+      printf("\n%f ", offset);
+   }
+   x.ProjectCoefficient(xcoeff); // TODO: call Compute?
+
+/*   if (visualization && vis[0].good())
+   {
+      vis[0].precision(8);
+      vis[0] << "solution" << endl << mesh << x << flush;
+      vis[0] << "window_title '" << "Solution" << "'" << endl
+             << "window_geometry " << 0 << " " << 0 << " " << visw << " " << vish << endl
+             << "keys mgA" << endl;
+
+      vis[1].precision(8);
+      vis[1] << "mesh" << endl << mesh << flush;
+      vis[1] << "window_title '" << "Mesh" << "'" << endl
+             << "window_geometry "
+             << (vish + 10) << " " << 0
+             << " " << visw << " " << vish  << endl
+             << "keys mgA" << endl;
+
+      vis[2].precision(8);
+      vis[2] << "solution" << endl << mesh << x << flush;
+      vis[2] << "window_title '" << "Image" << "'" << endl
+             << "window_geometry "
+             <<  (2 * vish + 10) << " " << 0
+             << " " << visw << " " << vish << endl
+             << "keys RjgA" << endl; // mn
+   }*/
+
+   // 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.
+   //estimator.SetAnisotropic();
+
+   // 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.
+   refiner.SetTotalErrorFraction(0.7);
+}
+
+
+// *****************************************************************************
+int Drl4Amr::Compute()
+{
+   iteration ++;
+   const int cdofs = fespace.GetTrueVSize();
+   cout << "\nAMR iteration " << iteration << endl;
+   cout << "Number of unknowns: " << cdofs << endl;
+
+   // TODO: it would be more proper to actually solve here
+   x.ProjectCoefficient(xcoeff);
+
+   // constrain slave nodes
+   if (fespace.GetProlongationMatrix())
+   {
+      Vector y(fespace.GetTrueVSize());
+      fespace.GetRestrictionMatrix()->Mult(x, y);
+      fespace.GetProlongationMatrix()->Mult(y, x);
+   }
+
+   // Send solution by socket to the GLVis server.
+   if (visualization && vis[0].good())
+   {
+      vis[0] << "solution\n" << mesh << x << flush;
+      fflush(0);
+   }
+   if (cdofs > max_dofs)
+   {
+      cout << "Reached the maximum number of dofs. Stop." << endl;
+      exit(0);
+   }
+   return 0;
+}
+
+
+// *****************************************************************************
+int Drl4Amr::Refine(int el_to_refine)
+{
+   if (el_to_refine >= 0)
+   {
+      //mesh.PrintCharacteristics();
+      mesh.EnsureNCMesh();
+      const int depth = mesh.ncmesh->GetElementDepth(el_to_refine);
+      if (depth == max_depth)
+      {
+         return 0;
+      }
+      MFEM_VERIFY(depth <= max_depth, "max_amr_depth error");
+      //dbg("Refine el:%d, depth:%d", el_to_refine, depth);
+      Array<Refinement> refinements(1);
+      refinements[0] = Refinement(el_to_refine);
+      mesh.GeneralRefinement(refinements, 1, 1);
+      // Send solution by socket to the GLVis server.
+      if (visualization && vis[1].good())
+      {
+         vis[1] << "mesh\n" << mesh << flush;
+         fflush(0);
+      }
+   }
+   else
+   {
+      //dbg("Refine with refiner");
+      // 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;
+         return 1;
+      }
+   }
+
+   // Update the space to reflect the new state of the mesh.
+   fespace.Update();
+   x.Update();
+   return 0;
+}
+
+
+void Drl4Amr::RandomRefine(double p)
+{
+   mesh.RandomRefinement(p);
+   fespace.Update();
+   x.Update();
+}
+
+
+// *****************************************************************************
+double Drl4Amr::GetNorm()
+{
+   //dbg("GetNorm");
+   // Setup all integration rules for any element type
+   const int order_quad = max(2, 2*order+1);
+   const IntegrationRule *irs[Geometry::NumGeom];
+   for (int i=0; i < Geometry::NumGeom; ++i)
+   {
+      irs[i] = &(IntRules.Get(i, order_quad));
+   }
+   const double err_x  = x.ComputeL2Error(xcoeff, irs);
+   const double norm_x = ComputeLpNorm(2., xcoeff, mesh, irs);
+   return err_x / norm_x;
+}
+
+
+// *****************************************************************************
+double *Drl4Amr::GetFullImage()
+{
+   Array<int> vert;
+   Vector sln;
+
+   int width = GetFullWidth();
+   int height = GetFullHeight();
+
+   image.SetSize(width * height);
+
+   // rasterize each element into the image
+   for (int k = 0; k < mesh.GetNE(); k++)
+   {
+      int subdiv = (1 << (max_depth - mesh.ncmesh->GetElementDepth(k))) * order;
+
+      const IntegrationRule &ir =
+         GlobGeometryRefiner.Refine(Geometry::SQUARE, subdiv)->RefPts;
+
+      x.GetValues(k, ir, sln);
+
+      mesh.GetElementVertices(k, vert);
+      const double *v = mesh.GetVertex(vert[0]);
+
+      int ox = int(v[0] * nx*(1 << max_depth) * order);
+      int oy = int(v[1] * ny*(1 << max_depth) * order);
+
+      for (int i = 0; i <= subdiv; i++)
+      for (int j = 0; j <= subdiv; j++)
+      {
+         int n = i*(subdiv+1) + j;
+         int m = (oy + i)*width + (ox + j);
+
+         image(m) = sln(n);
+      }
+   }
+
+   if (visualization) { ShowFullImage(); }
+
+   return image.GetData();
+}
+
+// *****************************************************************************
+double* Drl4Amr::GetLocalImage(int element)
+{
+   NCMesh *ncmesh = mesh.ncmesh;
+   MFEM_ASSERT(ncmesh, "");
+
+   int width = GetLocalWidth();
+   int height = GetLocalHeight();
+   int subdiv = oversample*order;
+
+   local_image.SetSize(width*height*4);
+   local_image = 0.0;
+
+   auto im = Reshape(local_image.Write(false), width, height, 4);
+
+   // find neighbors of the current element
+   // TODO: fix the protected access hack
+   Array<int> neighbors;
+   ncmesh->FindNeighbors(ncmesh->leaf_elements[element], neighbors);
+   for (int i = 0; i < neighbors.Size(); i++)
+   {
+      neighbors[i] = ncmesh->elements[neighbors[i]].index;
+   }
+
+   // add the element itself at the end of the list
+   neighbors.Append(element);
+
+   Array<IsoparametricTransformation*> elemT(neighbors.Size());
+   Array<InverseElementTransformation*> invT(neighbors.Size());
+
+   // prepare forward and inverse transforms for all the elements
+   for (int i = 0; i < neighbors.Size(); i++)
+   {
+      elemT[i] = new IsoparametricTransformation;
+      mesh.GetElementTransformation(neighbors[i], elemT[i]);
+      invT[i] = new InverseElementTransformation(elemT[i]);
+   }
+
+   Vector phys, grad;
+
+   // sample solution at a regular grid of points centered on 'element'
+   for (int y = 0; y < height; y++)
+   for (int x = 0; x < width; x++)
+   {
+      IntegrationPoint ip;
+      ip.Set2((x - context + 0.5) / subdiv,
+              (y - context + 0.5) / subdiv);
+
+      // use the local element to transform forward
+      elemT.Last()->Transform(ip, phys);
+
+      // see what neighbor, if any, the physical point lies in
+      for (int k = 0; k < neighbors.Size(); k++)
+      {
+         if (invT[k]->Transform(phys, ip)
+             == InverseElementTransformation::Inside)
+         {
+            im(x, y, 0) = this->x.GetValue(neighbors[k], ip);
+
+            elemT[k]->SetIntPoint(&ip);
+            this->x.GetGradient(*elemT[k], grad);
+
+            im(x, y, 1) = grad(0);
+            im(x, y, 2) = grad(1);
+            im(x, y, 3) = ncmesh->GetElementDepth(neighbors[k])
+                          - ncmesh->GetElementDepth(element);
+            break;
+         }
+      }
+   }
+
+   for (int i = 0; i < neighbors.Size(); i++)
+   {
+      delete elemT[i];
+      delete invT[i];
+   }
+
+   return local_image.GetData();
+}
+
+
+// *****************************************************************************
+void Drl4Amr::ShowFullImage()
+{
+   if (!vis[2].good()) { return; }
+
+   Mesh imesh(GetFullWidth(), GetFullHeight(), elem_type, false, sx, sy, false);
+
+   L2_FECollection fec(0, imesh.Dimension());
+   FiniteElementSpace fes(&imesh, &fec);
+   GridFunction gridfn(&fes, image.GetData());
+
+   vis[2] << "solution" << endl << imesh << gridfn << flush;
+}
+
+// *****************************************************************************
+extern "C" {
+
+Drl4Amr* Ctrl(int order) { return new Drl4Amr(order); }
+int Compute(Drl4Amr *ctrl) { return ctrl->Compute(); }
+int Refine(Drl4Amr *ctrl, int el) { return ctrl->Refine(el); }
+int GetNDofs(Drl4Amr *ctrl) { return ctrl->GetNDofs(); }
+int GetNE(Drl4Amr *ctrl) { return ctrl->GetNE(); }
+double GetNorm(Drl4Amr *ctrl) { return ctrl->GetNorm(); }
+
+double *GetFullImage(Drl4Amr *ctrl) { return ctrl->GetFullImage(); }
+int GetFullWidth(Drl4Amr *ctrl) { return ctrl->GetFullWidth(); }
+int GetFullHeight(Drl4Amr *ctrl) { return ctrl->GetFullHeight(); }
+
+} // extern "C"

miniapps/drl4amr/drl4amr.hpp

@@ -0,0 +1,70 @@
+#ifndef DRL4AMR_HPP
+#define DRL4AMR_HPP
+
+#include "mfem.hpp"
+
+using namespace std;
+using namespace mfem;
+
+class Drl4Amr
+{
+private:
+   const int nx = 8;
+   const int ny = 8;
+   const int max_dofs = 50000;
+   const int max_depth = 2;
+   const Element::Type elem_type = Element::QUADRILATERAL;
+   const double sx = 1.0;
+   const double sy = 1.0;
+
+   const int oversample = 8;
+   const int context = 4;
+
+   Mesh mesh;
+
+   int order;
+   int dim;
+
+   Device device;
+   socketstream vis[3];
+
+   Vector image;
+   Vector local_image;
+
+   H1_FECollection fec;
+   FiniteElementSpace fespace;
+
+   ConstantCoefficient one;
+   ConstantCoefficient zero;
+   BilinearFormIntegrator *integ;
+   FunctionCoefficient xcoeff;
+
+   GridFunction x;
+   int iteration;
+
+   FiniteElementSpace flux_fespace;
+   ZienkiewiczZhuEstimator estimator;
+   ThresholdRefiner refiner;
+
+public:
+   Drl4Amr(int order);
+
+   int Compute();
+   int Refine(int el =-1);
+   int GetNDofs() { return fespace.GetNDofs();}
+   int GetNE() { return fespace.GetNE(); }
+   double GetNorm();
+
+   double *GetFullImage();
+   int GetFullWidth() const { return 1 + order * (nx << max_depth); }
+   int GetFullHeight() const { return GetFullWidth(); }
+   void ShowFullImage();
+
+   double* GetLocalImage(int element);
+   int GetLocalWidth() const { return oversample*order + 2*context; }
+   int GetLocalHeight() const { return GetLocalWidth(); }
+
+   void RandomRefine(double p = 0.5);
+};
+
+#endif // DRL4AMR_HPP

miniapps/drl4amr/drl4amr.py

@@ -0,0 +1,76 @@
+import numpy as np
+from ctypes import *
+from random import *
+from PIL import Image
+
+c_double_p = POINTER(c_double)
+
+MFEM = cdll.LoadLibrary('./libdrl4amr.so')
+
+def sign(f, args_t, res_t):
+    f.restype = res_t
+    f.argtypes = args_t
+
+class Ctrl(object):
+    def __init__(self, order): self.obj = MFEM.Ctrl(order)
+    sign(MFEM.Ctrl, [c_int], c_void_p)
+
+    def Compute(self): return MFEM.Compute(self.obj)
+    sign(MFEM.Compute, [c_void_p], c_int)
+
+    def Refine(self, el_to_refine): return MFEM.Refine(self.obj, el_to_refine)
+    sign(MFEM.Refine, [c_void_p, c_int], c_int)
+
+    def GetNE(self): return MFEM.GetNE(self.obj)
+    sign(MFEM.GetNE, [c_void_p], c_int)
+
+    def GetNorm(self): return MFEM.GetNorm(self.obj)
+    sign(MFEM.GetNorm, [c_void_p], c_double)
+
+    def GetNDofs(self): return MFEM.GetNDofs(self.obj)
+    sign(MFEM.GetNDofs, [c_void_p], c_int)
+
+    def GetImage(self): return MFEM.GetImage(self.obj)
+    sign(MFEM.GetImage, [c_void_p], c_double_p)
+
+    def GetImageSize(self): return MFEM.GetImageSize(self.obj)
+    sign(MFEM.GetImageSize, [c_void_p], c_int)
+
+    def GetImageX(self): return MFEM.GetImageX(self.obj)
+    sign(MFEM.GetImageX, [c_void_p], c_int)
+
+    def GetImageY(self): return MFEM.GetImageY(self.obj)
+    sign(MFEM.GetImageY, [c_void_p], c_int)
+
+
+order = 2
+
+sim = Ctrl(order)
+
+NE = sim.GetNE()
+NX = sim.GetImageX()
+NY = sim.GetImageY()
+#print(NX,NY)
+
+while sim.GetNorm() > 0.01:
+    sim.Compute()
+
+    #sim.Refine(-1); # Will refine using the internal refiner
+    sim.Refine(int(NE*random()))
+
+    image_d = sim.GetImage()
+    #image_s = sim.GetImageSize()
+    #print("image_s: " + str(image_s))
+
+    # Get the data back, two ways:
+    #data = np.fromiter(image_d, dtype=np.double, count=NX*NY) # copy
+    # or:
+    data = np.frombuffer((c_double*NX*NY).from_address(addressof(image_d.contents))) # address
+
+    # Scale and convert the image
+    image_f = (data * 255 / np.max(data)).astype('uint8')
+    image = Image.fromarray(image_f.reshape(NX,NY),'L')
+    image.save('image.jpg')
+    print("Norm: "+str(sim.GetNorm()))
+
+print("Done, final norm: "+str(sim.GetNorm()))

miniapps/drl4amr/main.cpp

@@ -0,0 +1,66 @@
+#include "mfem.hpp"
+
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+#include "drl4amr.hpp"
+
+
+void ShowImg(socketstream &sock, double *img, int width, int height,
+             const char *title = NULL)
+{
+   Mesh mesh(width, height, Element::QUADRILATERAL, true, 1.0, 1.0, false);
+   L2_FECollection fec(0, 2);
+   FiniteElementSpace fes(&mesh, &fec);
+   GridFunction gf(&fes, img);
+   sock << "solution\n" << mesh << gf;
+   sock << "keys Rjlm\n";
+   if (title)
+   {
+      sock << "window_title '" << title << "'\n";
+   }
+   sock << flush;
+}
+
+
+int main(int argc, char *argv[])
+{
+   const int order = 2;
+   Drl4Amr sim(order);
+
+   socketstream vis1("localhost", 19916);
+   socketstream vis2("localhost", 19916);
+   socketstream vis3("localhost", 19916);
+   socketstream vis4("localhost", 19916);
+
+   //sim.Refine(9);
+   sim.RandomRefine();
+   sim.RandomRefine();
+   sim.RandomRefine();
+
+   sim.Compute();
+
+   while (sim.GetNorm() > 0.01)
+   {
+#if 0
+      const int el = static_cast<int>(drand48()*sim.GetNE());
+      sim.Compute();
+      sim.Refine(el);
+      sim.GetFullImage();
+      sim.GetFullWidth();
+#else
+      double *img = sim.GetLocalImage(142);
+      int w = sim.GetLocalWidth();
+      int h = sim.GetLocalWidth();
+      ShowImg(vis1, img, w, h, "Solution");
+      ShowImg(vis2, img + w*h, w, h, "dx");
+      ShowImg(vis3, img + 2*w*h, w, h, "dy");
+      ShowImg(vis4, img + 3*w*h, w, h, "Depth");
+      break;
+#endif
+   }
+   return 0;
+}

miniapps/drl4amr/makefile

@@ -0,0 +1,69 @@
+# Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at the
+# Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights reserved.
+# See file COPYRIGHT for details.
+#
+# This file is part of the MFEM library. For more information and source code
+# availability see http://mfem.org.
+#
+# MFEM is free software; you can redistribute it and/or modify it under the
+# terms of the GNU Lesser General Public License (as published by the Free
+# Software Foundation) version 2.1 dated February 1999.
+
+# Use the MFEM build directory
+MFEM_DIR ?= ../..
+MFEM_BUILD_DIR ?= ../..
+SRC = $(if $(MFEM_DIR:../..=),$(MFEM_DIR)/miniapps/drl4amr/,)
+CONFIG_MK = $(MFEM_BUILD_DIR)/config/config.mk
+# Use the MFEM install directory
+# MFEM_INSTALL_DIR = ../../mfem
+# CONFIG_MK = $(MFEM_INSTALL_DIR)/share/mfem/config.mk
+
+MFEM_LIB_FILE = mfem_is_not_built
+-include $(CONFIG_MK)
+
+# Some choices below are based on the OS type:
+NOTMAC := $(subst Darwin,,$(shell uname -s))
+ifneq ($(NOTMAC),)
+   SOPREFIX = so
+else
+   SOPREFIX = install_
+endif
+
+SEQ_MINIOBJS = drl4amr.o
+PAR_MINIOBJS =
+ifeq ($(MFEM_USE_MPI),NO)
+   MINIOBJS = $(SEQ_MINIOBJS)
+else
+   MINIOBJS = $(PAR_MINIOBJS) $(SEQ_MINIOBJS)
+endif
+
+.SUFFIXES:
+.SUFFIXES: .o .cpp .mk
+.PHONY: all cln clean clean-build clean-exec
+
+# Remove built-in rule
+%.o: %.cpp
+
+# Replace the default implicit rule for *.cpp files
+%.o: $(SRC)%.cpp $(SRC)%.hpp $(CONFIG_MK)
+	$(MFEM_CXX) -c -fPIC $(MFEM_FLAGS) $< -o $@
+
+all: libdrl4amr.so drl4amr
+
+go:;./drl4amr
+py python:;python drl4amr.py
+tst test:;$(MAKE) go && $(MAKE) python
+
+drl4amr: main.cpp libdrl4amr.so
+	$(MFEM_CXX) $(MFEM_FLAGS) main.cpp -o drl4amr -L. -Wl,-rpath,. -ldrl4amr
+
+libdrl4amr.so: $(MINIOBJS)
+	$(MFEM_CXX) -shared -Wl,-$(SOPREFIX)name,libdrl4amr.so -o libdrl4amr.so $(MINIOBJS) $(MFEM_LIBS)
+
+# Generate an error message if the MFEM library is not built and exit
+$(MFEM_LIB_FILE):
+	$(error The MFEM library is not built)
+
+cln clean:
+	rm -f *.o *~
+	rm -rf *.dSYM *.TVD.*breakpoints $(SEQ_MINIOBJS) drl4amr libdrl4amr.so