do not use tkinter in show_problem

paraview_utils/show_problem.py

 from os.path import isfile
-import tkinter as tk
-from tkinter.filedialog import askdirectory
-
-# Create the box to select the output directory
-boxroot = tk.Tk()
-boxroot.withdraw()
-dirname = askdirectory(parent=boxroot)
-boxroot.destroy()
-
+import os.path
 from paraview.simple import *
+try :
+    dirname = os.path.dirname(GetActiveSource().FileName)
+except:
+    raise ValueError("Cannot get the directory of the active source.")
+Delete()
 # Show the fluid velocity
 def showFluid(fluid, renderView, animationScene):
     fluidDisplay = Show(fluid, renderView)
@@ -65,7 +62,7 @@ def showContacts(particleProblem, renderView):
         RenameSource('Contacts', contacts)
         return contacts
     except NameError:
-        print("MigFlow extract contacts Filter has not been charged !\nPlease see the wiki to add it in your Paraview filter.") 
+        print("MigFlow extract contacts Filter has not been charged !\nPlease see the wiki to add it in your Paraview filter.")
         return None
 
 def showFluidVelocity(fluid, renderView):
@@ -85,7 +82,7 @@ def showFluidVelocity(fluid, renderView):
 animationScene = GetAnimationScene()
 animationScene.PlayMode = 'Snap To TimeSteps'
 renderView = GetActiveViewOrCreate('RenderView')
-renderView.InteractionMode = '3D'
+renderView.InteractionMode = '2D'
 
 # Read the pvd files
 particleProblemFile = str(dirname)+"/particle_problem.pvd"
@@ -114,4 +111,4 @@ if isfile(fluidFile):
 # Show the first time step
 animationScene.GoToFirst()
 renderView.Update()
-renderView.ResetCamera()
\ No newline at end of file
+renderView.ResetCamera()

python/petsc4pylsys.py

+import petsc4py
+import sys
+petsc4py.init(sys.argv)
+from petsc4py import PETSc
+import numpy as np
+import atexit
+from ._tools import timeit
+
+def gen_csr(idx) :
+    pairs = np.ndarray([idx.shape[0],idx.shape[1],idx.shape[1]],dtype=([('i0',np.int32),('i1',np.int32)]))
+    pairs['i0'][:,:,:] = idx[:,:,None]
+    pairs['i1'][:,:,:] = idx[:,None,:]
+    pairs, csrmap = np.unique(pairs,return_inverse=True)
+    csr_rowidx = np.hstack([np.array([0],dtype=np.int32),
+                            np.cumsum(np.bincount(pairs["i0"]),
+                           dtype=np.int32)])
+    csr_col = pairs['i1'].reshape(-1)
+    return csr_rowidx, csr_col, csrmap
+
+class LinearSystemAIJ:
+
+    def __init__(self, elements, n_fields, options, constraints = []) :
+        idx = ((elements*n_fields)[:,:,None]+np.arange(n_fields)[None,None,:]).reshape([elements.shape[0],-1])
+        self.localsize = elements.shape[1]*n_fields
+        self.constraints = list(c[:,0]*n_fields + c[:,1] for c in constraints)
+        self.idx = (elements[:,None,:]*n_fields+np.arange(n_fields)[None,:,None]).reshape([-1,9]).astype(np.int32)
+        self.idx2 = (elements[:,:,None]*n_fields+np.arange(n_fields)[None,None,:]).reshape([-1,9]).astype(np.int32)
+        options="-ksp_monitor -pc_type lu -info"
+        PETSc.Options().prefixPush("fluid_")
+        PETSc.Options().insertString(options)
+        PETSc.Options().prefixPop()
+        self.ksp = PETSc.KSP().create()
+        self.ksp.setOptionsPrefix(b"fluid_")
+        self.ksp.setFromOptions()
+        self.csr_rowidx, self.csr_col, self.csrmap = gen_csr(idx)
+        self.val = np.zeros(self.csr_col.size)
+        self.mat = PETSc.Mat()
+        self.mat.createAIJWithArrays(self.csr_rowidx.shape[0]-1,(self.csr_rowidx,self.csr_col,self.val),bsize=3)
+
+    #@timeit
+    #def local_to_global(self,localv,localm, u, constraints_value = []):
+    #    msize = self.csr_rowidx.shape[0]-1
+    #    rhs = np.bincount(self.idx,localv,msize)
+    #    self.val[:] = np.bincount(self.csrmap,localm,self.csr_col.size)
+    #    return rhs
+
+    @timeit
+    def local_to_global2(self,localv,localm, u, constraints_value = []):
+        self.mat.zeroEntries()
+        for e,m in zip(self.idx2,localm.reshape([-1,self.localsize**2])) :
+            self.mat.setValues(e,e,m,PETSc.InsertMode.ADD)
+        self.mat.assemble()
+
+    def local_to_global(self,localv,localm, u, constraints_value = []):
+        print(self.idx.shape)
+        msize = self.csr_rowidx.shape[0]-1
+        rhs = np.zeros((msize,))
+        #rhs = np.zeros((self.ndof + len(constraints_value,)))
+        np.add.at(rhs.reshape([-1]),self.idx.reshape([-1]),localv)
+        self.local_to_global2(localv,localm,u,constraints_value)
+        return rhs
+
+    @timeit
+    def solve(self,rhs) :
+        x = np.ndarray(rhs.shape)
+        prhs = PETSc.Vec().createWithArray(rhs.reshape([-1]))
+        px = PETSc.Vec().createWithArray(x.reshape([-1]))
+        #self.ksp.getPC().setOperators(None)
+        #self.ksp.getPC().setOperators(self.mat)
+        #self.ksp.setOperators(None)
+        self.ksp.setOperators(self.mat)
+        self.ksp.solve(prhs,px)
+        return x
+    
+
+class LinearSystemBAIJ :
+    def __init__(self, elements, n_fields, options, constraints = []) :
+        nnodes = np.max(elements)+1
+        nn = elements.shape[1]
+        pairs = np.ndarray((elements.shape[0],nn*(nn-1)),dtype=([('i0',np.int32),('i1',np.int32)]))
+        pairs['i0'] = elements[:,np.repeat(range(nn),nn-1)]
+        idx = np.hstack([[i for i in range(nn) if i!=j] for j in range(nn)]) 
+        pairs['i1'] = elements[:,idx]
+        pairs = np.unique(pairs)
+        nnz = (np.bincount(pairs["i0"])+1).astype(np.int32)
+        PETSc.Options().prefixPush("fluid_")
+        PETSc.Options().insertString(options)
+        PETSc.Options().prefixPop()
+        self.mat = PETSc.Mat().createBAIJ(nnodes*n_fields,n_fields,nnz)
+        self.ksp = PETSc.KSP().create()
+        self.ksp.setOptionsPrefix(b"fluid_")
+        self.ksp.setFromOptions()
+        self.mat.zeroEntries()
+        self.constraints = list(c[:,0]*n_fields + c[:,1] for c in constraints)
+        self.ndof = (np.max(elements) + 1)*n_fields
+        ###
+        self.localsize = elements.shape[1]*n_fields
+        self.elements = elements
+        self.idx = (self.elements[:,None,:]*n_fields+np.arange(n_fields)[None,:,None]).reshape([-1])
+
+    @timeit
+    def local_to_global2(self,localv,localm, u, constraints_value = []):
+        self.mat.zeroEntries()
+        for e,m in zip(self.elements,localm.reshape([-1,self.localsize**2])) :
+            self.mat.setValuesBlocked(e,e,m,PETSc.InsertMode.ADD)
+        self.mat.assemble()
+
+    def local_to_global(self,localv,localm, u, constraints_value = []):
+        rhs = np.zeros((self.ndof + len(constraints_value,)))
+        np.add.at(rhs.reshape([-1]),self.idx,localv)
+        self.local_to_global2(localv,localm,u,constraints_value)
+        return rhs
+
+    @timeit
+    def solve(self,rhs) :
+        viewer = PETSc.Viewer.ASCII("matinfo")
+        viewer.pushFormat(5)
+        self.mat.view(viewer)
+        x = np.ndarray(rhs.shape)
+        prhs = PETSc.Vec().createWithArray(rhs.reshape([-1]))
+        px = PETSc.Vec().createWithArray(x.reshape([-1]))
+        self.ksp.setOperators(self.mat)
+        self.ksp.solve(prhs,px)
+        return x
+
+LinearSystem = LinearSystemAIJ

python/petsclsys.py

-# MigFlow - Copyright (C) <2010-2020>
+# MigFlow - Copyright (C) <2010-2018>
 # <Universite catholique de Louvain (UCL), Belgium
 #  Universite de Montpellier, France>
 # 	
@@ -17,50 +17,214 @@
 # 
 # You should have received a copy of the GNU Lesser General Public License
 # along with this program (see COPYING and COPYING.LESSER files).  If not, 
-# see <http://www.gnu.org/licenses/>.
-import petsc4py
-#petsc4py.init()
-from petsc4py import PETSc
+# see self.mat<http://www.gnu.org/licenses/>.
 import numpy as np
-import time
-class LinearSystem :
-    def __init__(self, elements, n_fields, options) :
+import atexit
+import weakref
+
+from ._tools import timeit
+
+import ctypes
+import ctypes.util
+import os
+
+print(os.environ["LD_LIBRARY_PATH"])
+petscpath = os.environ["PETSC_DIR"]+"/"+os.environ["PETSC_ARCH"]+"/lib"
+os.environ["LD_LIBRARY_PATH"] += ":"+petscpath
+print(os.environ["LD_LIBRARY_PATH"])
+lib = np.ctypeslib.load_library("libpetsc",petscpath)
+
+lib.PetscInitialize(None, None, None)
+COMM_WORLD = ctypes.c_void_p.in_dll(lib,"PETSC_COMM_WORLD")
+#atexit.register(lib.PetscFinalize)
+
+def _np2c(a,dtype=float) :
+    tmp = np.require(a,dtype,"C")
+    r = ctypes.c_void_p(tmp.ctypes.data)
+    r.tmp = tmp
+    return r
+
+
+class PetscObj:
+    
+    def __init__(self, constructor, destructor, *args):
+        self._p = ctypes.c_void_p()
+        self._destructor = destructor
+        lib[constructor](*args,ctypes.byref(self._p))
+
+    def __del__(self):
+        lib[self._destructor](ctypes.byref(self._p))
+
+    @property
+    def _as_parameter_(self) :
+        return self._p
+
+
+class Viewer(PetscObj) :
+    def __init__(self):
+        super().__init__("PetscViewerASCIIOpen","PetscViewerDestroy",COMM_WORLD,"matinfo".encode())
+        lib.PetscViewerPushFormat(self,4)
+
+
+class KSP(PetscObj) :
+
+    def __init__(self, options) :
+        super().__init__("KSPCreate","KSPDestroy",COMM_WORLD)
+        lib.PetscOptionsPrefixPush(None, b"fluid_")
+        lib.PetscOptionsInsertString(None, options.encode())
+        lib.PetscOptionsInsertString(None, b"-pc_type lu -ksp_monitor")
+        lib.PetscOptionsPrefixPop(None)
+        lib.KSPSetOptionsPrefix(self,b"fluid_")
+        lib.KSPSetFromOptions(self)
+
+    def solve(self, mat, rhs):
+        x = np.zeros_like(rhs)
+        #lib.KSPSetReusePreconditioner(self, 0)
+        lib.KSPSetOperators(self, mat, mat)
+        lib.KSPSolve(self, Vec(rhs), Vec(x))
+        return x
+
+
+class MatSeqBAIJ(PetscObj) :
+
+    def __init__(self, ms, bs, csrrow, csrcol) :
+        super().__init__("MatCreate","MatDestroy",COMM_WORLD)
+        lib.MatSetType(self, b"seqbaij")
+        lib.MatSetSizes(self,ctypes.c_int(ms),ctypes.c_int(ms),ctypes.c_int(ms),ctypes.c_int(ms))
+        lib.MatSeqBAIJSetPreallocationCSR(self,
+                ctypes.c_int(bs), _np2c(csrrow, np.int32), _np2c(csrcol, np.int32), None)
+        self.ms = ms
+
+    def set(self, v) :
+        aptr = ctypes.POINTER(ctypes.c_double)()
+        lib.MatSeqBAIJGetArray(self, ctypes.byref(aptr))
+        np.ctypeslib.as_array(aptr,shape=(v.size,))[:] = v
+        lib.MatSeqBAIJRestoreArray(self, ctypes.byref(aptr))
+        lib.MatAssemblyBegin(self,0)
+        lib.MatAssemblyEnd(self,0)
+
+
+class MatSeqAIJ(PetscObj) :
+
+    def __init__(self, csr) :
+        super().__init__("MatCreate","MatDestroy",COMM_WORLD)
+        lib.MatSetType(self, b"seqaij")
+        lib.MatSetSizes(self,
+                ctypes.c_int(csr.size),ctypes.c_int(csr.size),
+                ctypes.c_int(csr.size),ctypes.c_int(csr.size))
+        lib.MatSeqAIJSetPreallocationCSR(self, _np2c(csr.row,np.int32), _np2c(csr.col,np.int32), None)
+
+    def set(self, v) :
+        aptr = ctypes.POINTER(ctypes.c_double)()
+        lib.MatSeqAIJGetArray(self, ctypes.byref(aptr))
+        np.ctypeslib.as_array(aptr,shape=(v.size,))[:] = v
+        lib.MatSeqAIJRestoreArray(self, ctypes.byref(aptr))
+        lib.MatAssemblyBegin(self,0)
+        lib.MatAssemblyEnd(self,0)
+
+
+class Vec(PetscObj) :
+
+    def __init__(self, v) :
+        super().__init__("VecCreateSeqWithArray","VecDestroy",COMM_WORLD,
+                ctypes.c_int(v.shape[-1]), ctypes.c_int(v.size),ctypes.c_void_p(v.ctypes.data))
+
+
+class CSR :
+
+    def __init__(self, idx, rhsidx,constraints):
+        pairs = np.ndarray([idx.shape[0],idx.shape[1],idx.shape[1]],dtype=([('i0',np.int32),('i1',np.int32)]))
+        pairs['i0'][:,:,:] = idx[:,:,None]
+        pairs['i1'][:,:,:] = idx[:,None,:]
+        pairs = pairs.reshape([-1])
+        allpairs = [pairs.reshape([-1])]
+        num = np.max(idx)
+        self.constraints = constraints
+        for c in constraints :
+            num += 1
+            pairs = np.ndarray([c.size*2],dtype=([('i0',np.int32),('i1',np.int32)]))
+            pairs['i0'][:c.size] = c
+            pairs['i1'][:c.size] = num
+            pairs['i0'][c.size:] = num
+            pairs['i1'][c.size:] = c
+            allpairs.append(pairs)
+        pairs = np.vstack(allpairs)
+        pairs, pmap = np.unique(pairs,return_inverse=True)
+        self.map = []
+        count = 0
+        for p in allpairs :
+            self.map.append(pmap[count:count+p.size])
+            count += p.size
+        self.row = np.hstack([np.array([0],dtype=np.int32),
+                              np.cumsum(np.bincount(pairs["i0"]),
+                              dtype=np.int32)])
+        self.col = pairs['i1']
+        self.size = self.row.size-1
+        self.rhsidx = rhsidx
+
+    def assemble_rhs(self, localv, u, constraints_value) :
+        rhs = np.bincount(self.rhsidx,localv,self.size)
+        for i, (c, cv) in enumerate(zip(self.constraints, constraints_value)):
+            localm = np.hstack((localm, cv[0], cv[0]))
+            rhs[i+self.ndof] = cv[1] + np.sum(u[c]*cv[0])
+        return rhs
+
+    def assemble_mat(self, localm, constraints_value) :
+        m = np.bincount(self.map[0],localm,self.col.size)
+        for cmap, cv in zip (self.map[1:], constraints_value) :
+            m += np.bincount(np.hstack(cmap, [cv[0],cv[0]], self.col.size))
+        return m
+
+
+
+class LinearSystemAIJ:
+
+    def __init__(self, elements, n_fields, options, constraints = []) :
+        idx = ((elements*n_fields)[:,:,None]+np.arange(n_fields)[None,None,:]).reshape([elements.shape[0],-1])
+        self.localsize = elements.shape[1]*n_fields
+        self.constraints = list(c[:,0]*n_fields + c[:,1] for c in constraints)
+        self.idx = (elements[:,None,:]*n_fields+np.arange(n_fields)[None,:,None]).reshape([-1])
+        self.ksp = KSP(options)
+        self.csr = CSR(idx, self.idx,[])
+        self.mat = MatSeqAIJ(self.csr)
+
+    @timeit
+    def local_to_global(self,localv,localm, u, constraints_value = []):
+        self.mat.set(self.csr.assemble_mat(localm, constraints_value))
+        return self.csr.assemble_rhs(localv, u , constraints_value)
+
+    @timeit
+    def solve(self,rhs) :
+        return self.ksp.solve(self.mat, rhs.reshape([-1])).reshape(rhs.shape)
+
+
+class LinearSystemBAIJ :
+
+    def __init__(self, elements, n_fields, options, constraints = []) :
+        if len(constraints) != 0 :
+            raise ValueError("petsc BAIJ not implemented with constraints")
         nnodes = np.max(elements)+1
         nn = elements.shape[1]
-        pairs = np.ndarray((elements.shape[0],nn*(nn-1)),dtype=([('i0',np.int32),('i1',np.int32)]))
-        pairs['i0'] = elements[:,np.repeat(range(nn),nn-1)]
-        idx = np.hstack([[i for i in range(nn) if i!=j] for j in range(nn)]) 
-        pairs['i1'] = elements[:,idx]
-        pairs = np.unique(pairs)
-        nnz = (np.bincount(pairs["i0"])+1).astype(np.int32)
-        PETSc.Options().prefixPush("fluid_")
-        PETSc.Options().insertString("-pc_type lu")
-        PETSc.Options().insertString(options)
-        PETSc.Options().prefixPop()
-        self.mat = PETSc.Mat().createBAIJ(nnodes*n_fields,n_fields,nnz)
-        self.ksp = PETSc.KSP().create()
-        self.ksp.setOptionsPrefix(b"fluid_")
-        self.ksp.setFromOptions()
-        self.mat.zeroEntries()
-        ###
+        self.ksp = KSP(options)
+        self.csr = CSR(elements,[])
+        self.mat = MatSeqBAIJ(self.csr.size*n_fields, n_fields, self.csr.row, self.csr.col)
         self.localsize = elements.shape[1]*n_fields
         self.elements = elements
         self.idx = (self.elements[:,None,:]*n_fields+np.arange(n_fields)[None,:,None]).reshape([-1])
-        self.globalsize = nnodes*n_fields;
+        self.size = nnodes*n_fields
+        self.n_fields = n_fields
+        csrmap = ((self.csr.map[0]*9)[:,None,None]+np.arange(0,9,3)[None,None,:]+np.arange(0,3)[None,:,None])
+        self.csrmap = np.copy(np.swapaxes(csrmap.reshape([elements.shape[0],nn,nn, n_fields,n_fields]),2,3).reshape([-1]))
 
-    def local_to_global(self,localv,localm,rhs):
-        self.mat.zeroEntries()
-        np.add.at(rhs.reshape([-1]),self.idx,localv)
-        for e,m in zip(self.elements,localm.reshape([-1,self.localsize**2])) :
-            self.mat.setValuesBlocked(e,e,m,PETSc.InsertMode.ADD)
-        self.mat.assemble()
+    @timeit
+    def local_to_global(self,localv,localm, u, constraints_value = []):
+        rhs = np.bincount(self.idx,localv,self.csr.size*self.n_fields)
+        self.mat.set(np.bincount(self.csrmap,localm,self.csr.col.size*self.n_fields**2))
+        return rhs
 
+    @timeit
     def solve(self,rhs) :
-        tic = time.time()
-        x = np.ndarray(rhs.shape)
-        prhs = PETSc.Vec().createWithArray(rhs.reshape([-1]))
-        px = PETSc.Vec().createWithArray(x.reshape([-1]))
-        self.ksp.setOperators(self.mat)
-        self.ksp.solve(prhs,px)
-        print("solve time : ",time.time()-tic)
-        return x
+        return self.ksp.solve(self.mat, rhs.reshape([-1,self.n_fields])).reshape(rhs.shape)
+
+LinearSystem = LinearSystemAIJ
+