only use rhsidx pour csr

python/CMakeLists.txt

@@ -31,7 +31,7 @@ set(SRC
   petsc4pylsys.py
   scipylsys.py
   VTK.py
-  dof_numbering.py
+  # dof_numbering.py
 )
 
 if(BUILD_MUMPS)

python/csr.py

@@ -2,11 +2,11 @@ import numpy as np
 import time
 class CSR :
 
-    def __init__(self, idx, rhsidx,constraints):
+    def __init__(self, rhsidx,constraints):
         shift = 2**32
-        num = np.max(idx)
+        num = np.max(rhsidx)
         self.ndof = num+1
-        idx = rhsidx.reshape([idx.shape[0],-1]).astype('uint64')
+        idx = rhsidx.astype('uint64')
         pairs = (idx[:,:,None]*shift+idx[:,None,:]).flatten()
         allpairs = [pairs]
         self.constraints = constraints
@@ -29,7 +29,7 @@ class CSR :
                               dtype=np.int32)])
         self.col = (pairs%shift).astype(np.int32)
         self.size = self.row.size-1
-        self.rhsidx = rhsidx
+        self.rhsidx = rhsidx.reshape([-1])
 
     def assemble_rhs(self, localv, u, constraints_value) :
         rhs = np.bincount(self.rhsidx,localv,self.size)

python/petsc4pylsys.py

@@ -5,36 +5,12 @@ from petsc4py import PETSc
 import numpy as np
 from ._tools import timeit
 from .csr import CSR
-# try :
-    # from .dof_numbering import MeshCl, LagrangeSpace
-# except :
-    # print("dof_numbering not available !!!")
 class LinearSystemAIJ:
 
     def __init__(self, idx, 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.idx = idx
         self.localsize = idx.shape[1]
         self.constraints = list(c[:,0]*n_fields + c[:,1] for c in constraints)
-        elements = np.asarray(idx[:,:self.localsize//n_fields]/n_fields, np.int32)
-        idx_matrix = ((elements*n_fields)[:,:,None]+np.arange(n_fields)[None,None,:]).reshape([elements.shape[0],-1])
-        self.idx = (elements[:,None,:]*n_fields+np.arange(n_fields)[None,:,None]).reshape([-1])
-        
-        # self.idx = (elements[:,None,:]*n_fields+np.arange(n_fields)[None,:,None]).reshape([-1])
-        # self.space_v = LagrangeSpace(MeshCl(elements), order=1).element_nodes
-        # self.space_p = LagrangeSpace(MeshCl(elements), order=1).element_nodes
-        # dimension = n_fields-1
-        # nmax_v = self.space_v.max()+1
-        # self.localsize = dimension*self.space_v.shape[1] + self.space_p.shape[1]
-        # self.globalsize = nmax_v*dimension + self.space_p.max()+1
-        # format u0,u1,...,v0,v1,...,p0,p1,...
-        # self.idx = np.tile(self.space_v, dimension) \
-            # + (np.repeat(np.arange(dimension), self.space_v.shape[1])*nmax_v)[None,:]
-        # pmap = dimension*nmax_v + self.space_p
-        # self.idx = np.concatenate([self.idx, pmap], axis=1)
-        
-        
+        self.idx = idx
         PETSc.Options().prefixPush("fluid_")
         PETSc.Options().insertString("-pc_type ilu -pc_factor_levels 5 -pc_factor_shift_type INBLOCKS -pc_factor_shift_amount 1e-16")
         PETSc.Options().insertString(options)
@@ -43,7 +19,7 @@ class LinearSystemAIJ:
         self.ksp = PETSc.KSP().create()
         self.ksp.setOptionsPrefix(b"fluid_")
         self.ksp.setFromOptions()
-        self.csr = CSR(idx_matrix,self.idx,self.constraints)
+        self.csr = CSR(self.idx,self.constraints)
         self.val = np.zeros(self.csr.col.size)
         self.mat = PETSc.Mat().createAIJWithArrays(self.csr.size,(self.csr.row,self.csr.col,self.val),bsize=1)