'''''
PyMFEM
built on mfem 4.0 (commit SHA=4d900b0c5fd6352c92173e74678bcbeeb11c8691)
'''''

PyMFEM is a python wrapper for MFEM, ligith-weight FEM (finite element
method) library developed by LLNL (http://mfem.org).
PyMFEM is tested with python = 2.7, 3.7

This wrapper is meant for a rapid-prototyping of FEM program, and
is built using SWIG 3.0.12
With PyMFEM, a user can create c++ MFEM objects and call their
method from python. We strongly recommend to visit the MFEM web site
to find more detail of the MFEM libirary.

-- Module structure --
  <root>/mfem/ser : wrapper for serial MFEM
             /par : wrapper for parallel MFEM
	     /common : helper modules 
             /Makefile_tempates: example template
	     /examples : example scripts

-- INSTALLATION --

  Please see INSTALL.txt

-- Features --
 Following features are realized in PyMFEM using SWIG interface
 in order to use MFEM more conveientely from python. 
 Some of them may be e integrated to MFEM itself in future.
 
  4-1) mfem::Vector
     Constructor using python list/numpy array:
     
        v = mfem.Vector([1,2,3.])
        v = mfem.Vector(np.array([1,2,3,4.])

        Note that when list is passed. Contents of list is copied and data is
	owned by mfem::Vector. When numpy arra is passed. Data is not copied.
	mfem::Vector holds the passed numpy array as _link_to_data attribute,
	so that python does not garbage collect the numpy array untill mfem::Vector
	is deleted.

        Constructor also get SWIG object to <double *>. However, memory
	management for this pointer object is left to a user.
	
     Array element access:
         reading element: print v[0]
	 element assinment v[0] = 3

         negative index is also possible; v[-1] is equivalent to v[size-1]
	 
     Assign method: operator= is replaced by Assine method 
        v.Assign(1.) (equivalent to v=1 in c++)
        v.Assign(<some mfem::Vector>) (copy data from mfem vector)
        v.Assign(np.array([1., 2., 3.])  (copy data from np.array)
	
     Data Access as Numpy:  v.GetDataArray()

     In c++, a new partial view to existing vector is acquired as
	Vector v(vec.GetData() + sc, sc);
     In Python, one can call either
        v = mfem.Vector(vec, offset, size)
     or use slice
        vec[offset:offset+size]  # more python-ish...

     Note that this slice returns mere new nmemory view, sharing the same vector
     object. However, be carefull that, when slice requires non contigeous momory
     view, PyMFEM returns a new vector object NOT sharing the memory.
     If one wants to get the access to a memory view with non-trivial striding,
     use GetDataArray discussed below to get a numpy array pointing the
     moemory allocated for Vector, and use slice of the newly created numpy object.
     
     Vector::GetDataArray returns numpy object viewing the data stored in C++ Vector
     object. Note that the returned numpy array does not own the data.

  4-2) mfem::Coefficient
       PyCoefficient is derived from FunctionCoefficient
       PyCoefficientT is time-dependent version
       VectorPyCoefficient  is derveid from VectorFunctionCoefficient
       VectorPyCoefficientT is time-dependent version

       When using these class, inherit one of them and define
       cls::EvalValue method as follows

       class Jt(mfem.VectorPyCoefficient):
           def EvalValue(self, x):
               return [0.,0.,  np.cos(np.abs(x[2]-0.03)/0.03*np.pi/2)]
       class Sigma(mfem.PyCoefficient):
           def EvalValue(self, x):
               return 0.01   // need to return 

       then use it

       Cof_Jt = Jt(3) // VectorPyCoefficient constructor need sdim 
       dd = mfem.VectorFEBoundaryTangentLFIntegrator(Cof_Jt);

       If one wants to call Eval with a transformation object and
       an integration point, and one can inherit PyCoefficientBase
       (or their Vector/Matrix version) and overwrite Eval methods.
       These classes are defined as a director classes in c++.
       Therefore, if one overwrite Eval method, SWIG reroute the
       call to the Python side.

       V/M ConstantCoefficeint supports natrual python object as input
          MatrixConstantCoefficient([[1,2], [2,3]])
          VecotrConstantCoefficient([1,2,3])

          Note: thie argument is first tested if it can be converted using
	        np.array(x,  dtype=float)

  4-3) mfem::GridFunction
       GetNodalValues(i) will perform GetNodalValue(Vector(), i) and
       return numpy array of Vector()
  
       = operator is renamed as Assign method.
         (python)
	     g = mfem.GridFunction(fespace)
	     g.Assign(0)
         (c++)
	     GridFunction *g = new GridFunction(&fespace);
	     g = 0;
	     
      GridFunction(FiniteElementSpace *f, double *data) is extended to
      support offsets
          v = mfem.GridFunction(fec, <pointer to double>, offset)
      so that follwoing can be wrappeds
	  GridFucntion(fes, vec.GetData() + sc)
      the same is used for ParGridFunction

      GridFunction::Save(std::ostream &out) is wrapped as Save(filename, precision=8)
      QuadratureFunction::Save(std::ostream &out) is wrapped as Save(filename, precision=8)      
	     
  4-4) mfem::Mesh
         Wrapper of following methods are customized to
	 return either python list object or tuple of lists
	  
            GetBdrElementVertices(i)
            GetElementVertices(i)
            GetElementVEdges(i)
            GetBdrElementEdges(i)
            GetFaceEdges(i)
            GetEdgeVertices(i)
            GetFaceVertices(i)
            GetElementFaces(i)
	    GetBdrElementAdjacentElement()

         These methods returns an IsoparametricTransformation object
	 
            GetElementTransformation(i)
            GetBdrElementTransformation(i)
            GetFaceTransformation(i)
            GetEdgeTransformation(i)
	    
	 pointer passing in following methos are returned
	 as tuple
            elem1, elem2 = mesh.GetFaceElements(i)
	    
         Additional constructors:
	    # this one takes file name as string
	    Mesh(const char *mesh_file, int generate_edges, int refine,
                 bool fix_orientation = True)
		 
	    # this one takes element type as string
            Mesh(int nx, int ny, int nz, const char *type, int generate_edges = 0,
              double sx = 1.0, double sy = 1.0, double sz = 1.0)
            Mesh(int nx, int ny, const char *type, int generate_edges = 0,
              double sx = 1.0, double sy = 1.0)
	      
	    (note) an issue is Element::Type (c++ enum) is treated as int in
	    python, and therefore, the wrapper can not disingish it from the
	    following constructor.
              Mesh(int _Dim, int NVert, int NElem, int NBdrElem = 0,
    	           int _spaceDim= -1)
         
         GetVertexArray can be used to obtain Vertex point data
	 as numpy object

         GetBdrElementFace(i) returns tuple
         GetBdrArray(int idx) returns array of boundary element whose BdrAttribute is idx
         GetBdrAttributeArray() returns array of boundary attribute
         GetAttributeArray() returns array of aqttribute
         GetDomainArray(int idx) returns array of element whose Attribute is idx
	 
         SwapNodes returns nodes and onws_nodes as follows.
	    nodes, owns_nodes = mesh.SwapNodes(nodes, owns_nodes)
	    nodes, owns_nodes = pmesh.SwapNodes(nodes, owns_nodes)

         Following method filename
	    virtual void PrintXG(std::ostream &out = mfem::out) const;
            virtual void Print(std::ostream &out = mfem::out) const { Printer(out); }
            void PrintVTK(std::ostream &out);	    
         Following method accept empty argument (=std:cout) and filename
            virtual void PrintInfo(std::ostream &out = mfem::out)


  4-5) sparsemat
         RAP has two different implementations. One accept three references.
	 The other accept a pointer as a third argument, instead.
	 These two are not distinguishable from python. So, RAP_P and
	 RAP_R are added. P indicates the third argment is pointer.
	 
         Add functions are accessed as add_sparse.

         GetIArray, GetJArray, and GetDataArray. These methos gives numpy
	 array of CSR matrix data. It borrows data from SparseMatrix.
	 Therefore, be careful not to access after the matrix is freeed.

         Following method accept empty argument (=std:cout) and filename
           void Print(std::ostream &out = mfem::out, int width_ = 4) const;
           void PrintMatlab(std::ostream &out = mfem::out) const;
           void PrintMM(std::ostream &out = mfem::out) const;
           void PrintCSR(std::ostream &out) const;
           void PrintCSR2(std::ostream &out) const;
           void PrintInfo(std::ostream &out) const;

        Constructor accept numpy array
           indptr = np.array([0, 2, 3, 6], dtype=np.int32)
           indices = np.array([0, 2, 2, 0, 1, 2], dtype=np.int32)
           data = np.array([1, 2, 3, 4, 5, 6], dtype=float)
           mat = mfem.SparseMatrix([indptr, indices, data, 3,3])

  4-6) densemat
         Add functions are accessed as add_dense.
	 
	 elements are accessed similar to a regular python array
	 Assign method replaces = operator
         GetDataArray returns a memory view to the internal data as numpy array.
	 
         (example)
              x = DenseMatrix(3, 3)	 
	      x.Assign(0)
              v = Vector([1,2,3,4,5,6,7,8,9])
  	      m = DenseMatrix(v.GetData(), 3, 3)
	      x.Assign(m)

              x[i,j] = xxx
	      print x[i,j]

              # Assigning (copy) from numpy array
	      m = DenseMatrix(3,3)
	      m.Assign(np.arange(9.).reshape(3,3))
	      
              x = DenseTensor(2, 3, 5)
	      x.Assign(1) # set all 1
	      x[5] # returns DenseMatrix with col=2 and row=3
	      print x[1,2,3]
	      x[1,2,3]=3

              # NOTE the index order is different in the numpy array returned
	      # by DenseTensor::GetDataArray
              x[1 2,3] = x.GetDataArray()[3,1,2]
	      
	 Note that internally MFEM densmatrix is column major, while numpy
	 is row major. Therefore,
            v = Vector([1,2,3,4,5,6,7,8,9])
	    m = DenseMatrix(v.GetData(), 3, 3)
         will make
	    [1 4,7]
	    [2,5,8]
	    [3,6,9],
	 whereas
	    numpy.arange(9).reshape(3,3)
         would give you
	    [1 2,3]
	    [4,5,6]
	    [7,8,9].


	  
  4-7) estimators
        ZienkiewiczZhuEstimator and L2ZienkiewZhuEstimator
	have two versions of constructor. flux_fes
	can be eithor passed as pointer or reference. In python class, you
	need to pass 5th argument if you want to pass flux_fes as
	reference. By default,  own_flux_fes is  False, so a user
	can skip passing this keyword when passing by reference
	
	pass by reference 
          ZienkiewiczZhuEstimator(integ, sol, flux_fes, own_flux_fes = False)	  
        pass by pointer
          ZienkiewiczZhuEstimator(integ, sol, flux_fes, own_flux_fes = True)
	pass by reference 
          L2ZienkiewiczZhuEstimator(integ, sol, flux_fes, own_flux_fes = False)	  
        pass by pointer
          L2ZienkiewiczZhuEstimator(integ, sol, flux_fes, own_flux_fes = True)	  
	  

  4-8) operator
        PyTimeDependentOperatorBase and PyOperatorBase is added to
	allow for implementing those operator in python. A user
	can inherit these class and overwrite Mult in python. See
	ex9.py for example.
	
  4-9) ode
        Step uses reference passing for t and dt. To get the result
	from this method, use
	   t, dt = ode_solver.Step(u, t, dt)

  4-10) hypre
        Following methods are added trhough SWIG interface
        HypreParVector::GetPartitioningArray  returns Partitioning as numpy array
        HypreParMatrix::GetColPartArray  returns ColPart as numpy array
        HypreParMatrix::GetRowPartArray  returns RowPart as numpy array	
        HypreParMatrix::GetCooDataArray  returns data to construct local coo matrix

        mfem.common.chypre
        This module defines CHypreVec and CHhypreMat classes. Those are classes
	supports...
	    1) create HypreParCSR and HypreParVector from scipy.sparse.csr_matrix
	    and numpy.ndarray, respectively
	    2) convert HypreMatrix/Vector to numpy array or scipy.sparse matrix
	    3) handle complex number using a pair of Hypre object (real and imag)
        Class methods of these classes are named in a similar way to numpy array.

   4-11) socketstream

         socketstream in PyMFEM has send_text and send_solution. These
	 method also send endl and flush. One can use them as follows
	 
         (example)
           sock = mfem.socketstream("localhost", 19916)
           sock.precision(8)
           sock.send_text("parallel " + str(num_procs) +  " " + str(myid))    
           sock.send_solution(pmesh, x)
	 
         socketstream object support << operator partially, but not works
	 perfectly. This was done by adding __lshift__ and other methos to
	 socketstream class in .i file. We don't want to wrap the whole
	 iostream. Instead, we are adding methods used in examples such as
	 precision. Also, note that sock << flush in c++ should be rephrased
	 sock.flush()
	   
   4-12) std::ostream&

         methods which takes std::ostream& are wrapped so that it acccept
	 either empty argument (stdout) or filename.
	 Not all methos are wrapped yet. Look for OSTREAM_ADD_DEFAULT_FILE
	 (or OSTREAM_ADD_DEFAULT_STDOUT_FILE) macros in the SWIG interface
	 file to find which method is wrapped.

         (example)
	    sparsmat.PrintInfo()
       
   4-13) mfem::table
           GetRowList returns a version of GetRow which returns Python
	   list instead of int pointer.

   4-14) mfem::ParMesh
   
         ParMesh(comm, filename)
            load parallel file
	 
         ParPrintToFile(filename)
            write parallel file
	    
   4-14) mfem::ParGridFunction

         ParGridFunction(pmesh, filename) : read parallel GridFunction

         example to read ParMesh and ParGridFunction:
   	    smyid = '{:0>6d}'.format(myid)  # myid = MPI.COMM_WORLD.rank
            mesh1 = mfem.ParMesh(comm, 'wg.mesh.'+smyid)
            mesh1.ReorientTetMesh()
            x = mfem.ParGridFunction(mesh1, 'wg.gf.'+smyid)

   4-15) mfem:BlockNonlinearForm, mfem:ParBlockNonlinearForm
        
         Constructor accept tuple/list of (Par)FiniteElementSpace
	     R_space = mfem.FiniteElementSpace(...)
             W_space = mfem.FiniteElementSpace(...)
             spaces = [R_space, W_space]
     	     mfem.BlockNonlinearForm(spaces)

             note that user must pay attention that Python does not
	     garbage collect spaces. (this will be fix in future)

         Similarly, mfem::(Par)BlockNonlinearForm::SetEssentialBC()
	 accept tuple/list of intArray (= wrapped Array<int>) and
	 mfem::Vector

   4-16) mfem:ComplexOperator
         ownReal, onwImag for Constructor is optional.
	 hermitan keyword is used to specify the type of matrix

             op_complex = ComplexOperator(M1, M2, ownReal=False, ownImag=False,
	                                  hermitan=True)
             hermitan = False results in BLOCK_SYMMETRIC format

         Returned ComplexOperator hold the passed real and imag operators
	 so that Python does not garbage collect them.
	 If thisown of real/imag operatros are 1, a user has an option to
	 set ownReal, ownImag.

    4-17) strumpack
         if you have strumpack linked with MFEM, you can enable wrapper by
	 ENABLE_STRUMPACK option in make

	    make ENABLE_STRUMPACK=1
	    make ENABLE_STRUMPACK=1 STRUMPACK_INCLUDE=${HOME}/sandbox/include


         Here is sample script to use strumpack from PyMFEM (only in mfem.par)
	 
	    import mfem.par.strumpack as strmpk
            Arow = strmpk.STRUMPACKRowLocMatrix(A)
            args = ["--sp_hss_min_sep_size", "128", "--sp_enable_hss"]
            strumpack = strmpk.STRUMPACKSolver(args, MPI.COMM_WORLD)
            strumpack.SetPrintFactorStatistics(True)
            strumpack.SetPrintSolveStatistics(False)
            strumpack.SetKrylovSolver(strmpk.KrylovSolver_DIRECT);
            strumpack.SetReorderingStrategy(strmpk.ReorderingStrategy_METIS)
            strumpack.SetMC64Job(strmpk.MC64Job_NONE)
            # strumpack.SetSymmetricPattern(True)
            strumpack.SetOperator(Arow)
            strumpack.SetFromCommandLine()
            strumpack.Mult(B, X);
	     
-- Known issues --

  PyMFEM is an on-going effort. Not all MFEM functioality is propary
  wrapped yet. Presnetly all serial/parallel examples avaiable in
  the top-level mfem/example directory are rewritten in python.
  This would be sufficient for certain projects, but for other projects,
  you may find various features are missing. Please inform us if you
  notice something missing.

  1) t*.hpp files are not wrapped
  
  2) Mesh::SwapNodes is not working perfectly.
  
  3) Memory management.
      python uses garbage collection. Since wrapper does not know
      exactly how objects are used in c++ layer, MFEM object may be
      collected before being ready to be freed. We are addressing
      this by adjusting thisown flag or %newobject directive. This
      needs to be done method by method, and not yet completed.

  4) Visualization through visit is not included

  5) SuperLU and Petsc are not supported.
  
*This work is supported by US DoE contract DE-FC02-99ER54512