2-D FEM coordinate and connectivity matrices

import numpy as np
from tabulate import tabulate

Let’s first create the global coordinate vector. The columns refer to the nodes, and the rows to the x,y coordinates.

One issue is the index counting in python. Remember that the first element in an array is stored at index 0. We have two options: if we want to start counting our elements/nodes at 1, we always have to substrat 1 to get the storage index; or we use the less intuitive way of counting and have the first node/element be 0. Let’s use the second option.

Another is the direction in which we count. We can either count in the row direction, or in the collumn direction. These two ways are called column major or row major. Related to this is the way that data is stored in the memory. If the data is stored F contiguous, the data is stored in the column direction; if it is C contiguous then in the row direction. Matlab is F contiguous, while python is C contiguous (by default, it can handle both).

You can find more info here: https://stackoverflow.com/questions/26998223/what-is-the-difference-between-contiguous-and-non-contiguous-arrays

Bottom line is, in python we should count in the row direction.

Let’s make some simple examples:

A = np.arange(1,13)
print(A)
print(A.flags)

B = np.reshape(A,(3,4))
print(B)
print(B.flags)
print(np.take(B, 3)) # linear index

[ 1  2  3  4  5  6  7  8  9 10 11 12]
  C_CONTIGUOUS : True
  F_CONTIGUOUS : True
  OWNDATA : True
  WRITEABLE : True
  ALIGNED : True
  WRITEBACKIFCOPY : False
  UPDATEIFCOPY : False

[[ 1  2  3  4]
 [ 5  6  7  8]
 [ 9 10 11 12]]
  C_CONTIGUOUS : True
  F_CONTIGUOUS : False
  OWNDATA : False
  WRITEABLE : True
  ALIGNED : True
  WRITEBACKIFCOPY : False
  UPDATEIFCOPY : False

4

Array B is C_CONTIGUOUS and the natural way to store the data and to not make mistakes during reshapes is to count in the row direction.

Here are some useful links on indexing: https://kanoki.org/2020/07/05/numpy-index-array-with-another-array/ https://realpython.com/numpy-array-programming/

Alright, let’s do the coordinate vector. First node should have the index 0 and we count in row direction!

lx          = 1
ly          = 1
nx          = 5
ny          = 4
nnodel      = 4
dx          = lx/(nx-1)
dy          = ly/(ny-1)

nex         = nx-1
ney         = ny-1
nnod        = nx*ny
nel         = nex*ney
GCOORD      = np.zeros((nnod,2))

id = 0

#for i in ??:
#    for j in ??:
#        GCOORD[id,0] = ??
#        GCOORD[id,1] = ??
#        id          = id + 1

#print(tabulate(GCOORD, headers=['Node #', 'X', 'Y'], showindex=True))

# hide: the code in this cell is hidden by the author

  Node #     X         Y
--------  ----  --------
       0  0     0
       1  0.25  0
       2  0.5   0
       3  0.75  0
       4  1     0
       5  0     0.333333
       6  0.25  0.333333
       7  0.5   0.333333
       8  0.75  0.333333
       9  1     0.333333
      10  0     0.666667
      11  0.25  0.666667
      12  0.5   0.666667
      13  0.75  0.666667
      14  1     0.666667
      15  0     1
      16  0.25  1
      17  0.5   1
      18  0.75  1
      19  1     1

Now we do the connectivity. Here again the first element will have index 0 and we count in the row direction!

#E2N   = np.zeros((nel,nnodel), dtype=int) # connectivity matrix

#for iel in range(0,nel):
#    row        = ??? np.ceil()  / np.floor
#    ind        = ???
#    E2N[iel,:] = [ind, ind+1, ?, ?]

#print(tabulate(E2N, headers=['Element #', 'Node 1', 'Node 2', 'Node 3', 'Node 4'], showindex=True))

#hide: the code in this cell is hidden by the author

  Element #    Node 1    Node 2    Node 3    Node 4
-----------  --------  --------  --------  --------
          0         0         1         6         5
          1         1         2         7         6
          2         2         3         8         7
          3         3         4         9         8
          4         5         6        11        10
          5         6         7        12        11
          6         7         8        13        12
          7         8         9        14        13
          8        10        11        16        15
          9        11        12        17        16
         10        12        13        18        17
         11        13        14        19        18

If we now want to get the coordinates of a specific element, we can do the following

iel = 5 # current element

print(np.take(GCOORD[:,0], E2N[iel,:])) # x coordinates
print(np.take(GCOORD[:,1], E2N[iel,:])) # y coordinates

[0.25 0.5  0.5  0.25]
[0.33333333 0.33333333 0.66666667 0.66666667]