Merge pull request #133 from ToFuProject/Issue131_OperatorD0N1

Issue131 operator `D0N1` fixed + removed all `order='F'` from `operators` and `cropbs_flat`

Author: Didou09

bsplines2d/_class01_select.py

@@ -763,6 +763,7 @@ def _mesh2DRect_bsplines_knotscents(
             knots_per_bs_Z = knots_per_bs_Z[:, ind[1]]
 
         nknots = knots_per_bs_R.shape[0]
+        # TBC: reverse repeat and tile ?
         knots_per_bs_R = np.tile(knots_per_bs_R, (nknots, 1))
         knots_per_bs_Z = np.repeat(knots_per_bs_Z, nknots, axis=0)
 
@@ -779,6 +780,7 @@ def _mesh2DRect_bsplines_knotscents(
             cents_per_bs_Z = cents_per_bs_Z[:, ind[1]]
 
         ncents = cents_per_bs_R.shape[0]
+        # TBC: reverse repeat and tile ?
         cents_per_bs_R = np.tile(cents_per_bs_R, (ncents, 1))
         cents_per_bs_Z = np.repeat(cents_per_bs_Z, ncents, axis=0)
 
@@ -787,10 +789,11 @@ def _mesh2DRect_bsplines_knotscents(
 
     if return_knots is True and return_cents is True:
         out = (
-            (knots_per_bs_R, knots_per_bs_Z), (cents_per_bs_R, cents_per_bs_Z)
+            (knots_per_bs_R, knots_per_bs_Z),
+            (cents_per_bs_R, cents_per_bs_Z)
         )
     elif return_knots is True:
         out = (knots_per_bs_R, knots_per_bs_Z)
     else:
         out = (cents_per_bs_R, cents_per_bs_Z)
-    return out
+    return out

bsplines2d/_class02_bsplines_operators_rect.py

@@ -7,7 +7,6 @@
 # Common
 import numpy as np
 import scipy.sparse as scpsp
-import datastock as ds
 
 
 from . import _utils_bsplines_operators as _operators
@@ -67,8 +66,8 @@ def get_mesh2dRect_operators(
     # prepare
 
     nx, ny = knotsx_per_bs.shape[1], knotsy_per_bs.shape[1]
-    kR = np.tile(knotsx_per_bs, ny)
-    kZ = np.repeat(knotsy_per_bs, nx, axis=1)
+    kR = np.repeat(knotsx_per_bs, ny, axis=1)
+    kZ = np.tile(knotsy_per_bs, nx)
     nbs = nx*ny
 
     if cropbs_flat is None:
@@ -219,7 +218,8 @@ def get_mesh2dRect_operators(
         for ir in range(nx):
             for iz in range(ny):
 
-                iflat = ir + iz*nx
+                # iflat = ir + iz*nx
+                iflat = iz + ir*ny
                 if cropbs_flat is not False and not cropbs_flat[iflat]:
                     continue
 
@@ -238,8 +238,8 @@ def get_mesh2dRect_operators(
                     if cropbs_flat is not False and not cropbs_flat[jflat]:
                         continue
 
-                    jr = jflat % nx
-                    jz = jflat // nx
+                    jr = jflat // ny
+                    jz = jflat % ny
 
                     # store (i, j) and (j, i) (symmetric matrix)
                     if jr >= ir:
@@ -277,17 +277,12 @@ def get_mesh2dRect_operators(
             )
 
             # surface elements
-            dZ = np.repeat(knotsy_mult[1:] - knotsy_mult[:-1], nx)
+            dZ = np.tile(knotsy_mult[1:] - knotsy_mult[:-1], nx)
             if geometry == 'linear':
-                dR = np.tile(
-                    knotsx_mult[1:] - knotsx_mult[:-1],
-                    ny,
-                )
+                dR = knotsx_mult[1:] - knotsx_mult[:-1]
             else:
-                dR = np.tile(
-                    0.5*(knotsx_mult[1:]**2 - knotsx_mult[:-1]**2),
-                    ny,
-                )
+                dR = 0.5*(knotsx_mult[1:]**2 - knotsx_mult[:-1]**2)
+            dR = np.repeat(dR, ny)
 
             dS = dR*dZ
             if cropbs_flat is not False:
@@ -328,7 +323,7 @@ def get_mesh2dRect_operators(
             for ir in range(nx):
                 for iz in range(ny):
 
-                    iflat = ir + iz*nx
+                    iflat = iz + ir*ny
                     if cropbs_flat is not False and not cropbs_flat[iflat]:
                         continue
 
@@ -348,8 +343,8 @@ def get_mesh2dRect_operators(
                         if cropbs_flat is not False and not cropbs_flat[jflat]:
                             continue
 
-                        jr = jflat % nx
-                        jz = jflat // nx
+                        jr = jflat // ny
+                        jz = jflat % ny
 
                         # store (i, j) and (j, i) (symmetric matrix)
                         if jr >= ir:
@@ -401,7 +396,7 @@ def get_mesh2dRect_operators(
         for ir in range(nx):
             for iz in range(ny):
 
-                iflat = ir + iz*nx
+                iflat = iz + ir*ny
                 if cropbs_flat is not False and not cropbs_flat[iflat]:
                     continue
 
@@ -422,8 +417,8 @@ def get_mesh2dRect_operators(
                     if cropbs_flat is not False and not cropbs_flat[jflat]:
                         continue
 
-                    jr = jflat % nx
-                    jz = jflat // nx
+                    jr = jflat // ny
+                    jz = jflat % ny
 
                     # store (i, j) and (j, i) (symmetric matrix)
                     if jr >= ir:
@@ -463,4 +458,4 @@ def get_mesh2dRect_operators(
 
         raise NotImplementedError("Integral D3N2 not implemented for deg=3!")
 
-    return opmat, operator, geometry
+    return opmat, operator, geometry

bsplines2d/_class02_bsplines_rect.py

@@ -394,20 +394,20 @@ def _get_overlap(
 ):
     # nb of overlapping, inc. itself in 1d
     nbs0, nbs1 = shapebs
-    ind00 = np.tile(np.arange(0, nbs0), nbs1)
-    ind10 = np.repeat(np.arange(0, nbs1), nbs0)
+    ind00 = np.repeat(np.arange(0, nbs0), nbs1)
+    ind10 = np.tile(np.arange(0, nbs1), nbs0)
 
     # complete
     ntot = 2*deg + 1
 
-    add0= np.tile(np.arange(-deg, deg+1), ntot)
-    add1 = np.repeat(np.arange(-deg, deg+1), ntot)
+    add0 = np.repeat(np.arange(-deg, deg+1), ntot)
+    add1 = np.tile(np.arange(-deg, deg+1), ntot)
 
     inter0 = ind00[None, :] + add0[:, None]
     inter1 = ind10[None, :] + add1[:, None]
 
     # purge
-    inter = inter0 + inter1*nbs0
+    inter = inter1 + inter0*nbs1
     indneg = (
         (inter0 < 0) | (inter0 >= nbs0) | (inter1 < 0) | (inter1 >= nbs1)
     )
@@ -469,4 +469,4 @@ def get_bs_class(
         knots1=knots1,
         deg=deg,
         shapebs=shapebs,
-    )
+    )