khameelbm / pymfem Goto Github PK

'''''
PyMFEM
'''''
PyMFEM is a python wrapper for MFEM, ligith-weight FEM (finite element
method) library developed by LLNL (http://mfem.org).
This wrapper is meant for a rapid-prototyping of FEM program, and
is built using SWIG 3.0.
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 --

0: 3rd parth library
  0-1) mpi4py
         needed either you use parallel MFEM
  0-2) MFEM
	 serial version and parallel version need to be
	 compiled in separate directory.

         By default, Makefile assumes that MFEM is installed
	 under /usr/local/mfem-3.2 and serial version of MFEM
	 is installed under /usr/local/mfem-3.2ser.
	 A user can change it by editing Makefile.local

         < note about METIS version>
         MUMPS uses METIS5. Therefore, if you are going to use
	 this with MUMPS in PyMFEM_pi,
	 parallel MFEM is required to build MFEM with METIS5.
	 Open config.mk in mfem/config (see INSTALL in MFEM for
	 detail) change 
   	    MFEM_USE_MPI         = YES 
            MFEM_USE_METIS_5     = YES
	 You also need to specify METIS_DIR, METIS_OPT, METIS_LIB
	 Then, try.
	    make config
	    mamke -j	 

  0-2) HYPRE
	 for parallel MFEM is used

1:  recompile MFEM with -fPIC option
    This may not be necessary. (I didn't need it on MacOSX, 
    I don't know why)     

2: Prepare Makefile.local
   copy template from Makefile_template/* to repository root.
   rename it as Makefile.local
   edit it for your enviroment

   Leaving MFEMSER in Makefile.local  skips building serial version.
   Leaving MFEM in Makefile.local  skips building parallel version.   
   
2-1 (Optional): Build SWIG cxx files
This is needed only when you edit SWIG interface files.
   * make cxx

   or
   
   * make sercxx (sereal MFEM wrapepr)
   * make parcxx (parallel MFEM wrapper)
	
3: Compile Extensions
   * make
   or 
   * make ser (sereal MFEM wrapepr)
   * make par (parallel MFEM wrapper)

4) run example
   Add PyMFEM to your PYTHONPATH. In bash,
      > export PYTHONPATH=<hogehogehoge>:$PYTHONPATH
   Then, launch python and see if it import mfem
      > import mfem
      > import mfem.ser (serial)
      or
      > import mfem      
      > import mfem.par (parallel)
   
   Many of example c++ programs in MFEM are converted in python and
   found in example directory. Once PYTHONPATH is set, you
   can run them as follows.
     
     > python ex1.py  (serial)
     > mpiexec -n 2 python ex1p.py (parallel)

   Note that python version  examples are simplfied. It does not
   parse arguments nor visulaization via socket connection.

-- 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:  print v[0] or v[0] = 3
     Assign method:         v.Assign(1) (equivalent to v=1 in c++)
     Data Access as Numpy:  v.GetDataArray()


     In c++, a new partial view to existing vector is acquired
     as
	Vector v(vec.GetData() + sc, sc);
     , where as in python, In python, we need to
        v = mfem.Vector(vec, offset, size)

  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);

  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;
	     
  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)
	 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 whose BdrAttribute is idx
         GetBdrAttributeArray(int idx) returns array of boundary attribute

  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.	 
	 
  4-6) densemat
         Add functions are accessed as add_dense.
	 
  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 usins 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 object support << operator partially. This is 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()

         (example)
           sol_sock = mfem.socketstream("localhost", 19916)
           sol_sock.precision(8)
           sol_sock << "solution\n" << mesh << x
           sol_sock.flush()
	
-- Known issues --

  PyMFEM is a on-going effort. Not all MFEM functioality is propary
  wrapped yet. Presnetly most of serial example (except for
  ex15.cpp) and 5 parallel examples 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, but the approach may not 
      work in a different scenario.

  4) Visualization through visit is not included. Visulaization
     using socketsteram is on-going work.
  
*This work is supported by US DoE contract DE-FC02-99ER54512