CALFEM / calfem-python

CALFEM for Python is the Python port of the CALFEM finite element toolkit. It also implements meshing function based on GMSH and triangle. Visualisation routines are implemented using visvis and matplotlib.
MIT License
115 stars 189 forks source link

LinAlgError: Singular matrix by `a, r = cfc.solveq(K,f,bc)` #58

Closed Martin15135215 closed 10 months ago

Martin15135215 commented 11 months ago

This is my example truss:

image

This is the code:

import numpy as np
import calfem.core as cfc

# Topology matrix Edof
Edof = np.array([
    [1,2,3,4,5,6], # Element 1
    [1,2,3,7,8,9], # Element 2
    [1,2,3,10,11,12], # Element 3
])

K = np.matrix(np.zeros((12,12)))
f = np.matrix(np.zeros((12,1)))

# Area and E-Modulus
E = 1000 # N/mm²
A = 1000 # mm²
ep = [E,A]

#Element coordinates
ex1 = np.array([0., 0.])  # Element 1
ex2 = np.array([0., -4.]) # Element 2
ex3 = np.array([0., -4.]) # Element 3

ey1 = np.array([0., 0.]) # Element 1
ey2 = np.array([0., 0.]) # Element 2
ey3 = np.array([0., 0.]) # Element 3

ez1 = np.array([3., 0.])  # Element 1
ez2 = np.array([3., 3.])  # Element 2
ez3 = np.array([3., 6.])  # Element 3

# Element stiffness matrices
Ke1 = cfc.bar3e(ex1,ey1,ez1,ep)
Ke2 = cfc.bar3e(ex2,ey2,ez2,ep)
Ke3 = cfc.bar3e(ex3,ey3,ez3,ep)

# Assemble Ke into K
cfc.assem(Edof[0,:],K,Ke1)
cfc.assem(Edof[1,:],K,Ke2)
cfc.assem(Edof[2,:],K,Ke3)

# prescribed displacements in bc
bc = np.array([4,5,6,7,8,9,10,11,12])

f[1] = 100000 # N
f[2] = -100000 # N

a, r = cfc.solveq(K,f,bc)

When I run this, I get this error message:

LinAlgError                               Traceback (most recent call last)
Cell In[14], line 49
     46 f[1] = 100000 # N
     47 f[2] = -100000 # N
---> 49 a, r = cfc.solveq(K,f,bc)

File ~/.local/share/python3.11/site-packages/calfem/core.py:5317, in solveq(K, f, bcPrescr, bcVal)
   5313 bcDofs = bcDofs[bc]
   5315 fsys = f[bcDofs]-K[np.ix_((bcDofs), (bcPrescr-1))] * \
   5316     np.asmatrix(bcVal).reshape(nPdofs, 1)
-> 5317 asys = np.linalg.solve(K[np.ix_((bcDofs), (bcDofs))], fsys)
   5319 a = np.zeros([nDofs, 1])
   5320 a[np.ix_(bcPrescr-1)] = np.asmatrix(bcVal).reshape(nPdofs, 1)

File <__array_function__ internals>:200, in solve(*args, **kwargs)

File ~/.local/share/python3.11/site-packages/numpy/linalg/linalg.py:386, in solve(a, b)
    384 signature = 'DD->D' if isComplexType(t) else 'dd->d'
    385 extobj = get_linalg_error_extobj(_raise_linalgerror_singular)
--> 386 r = gufunc(a, b, signature=signature, extobj=extobj)
    388 return wrap(r.astype(result_t, copy=False))

File ~/.local/share/python3.11/site-packages/numpy/linalg/linalg.py:89, in _raise_linalgerror_singular(err, flag)
     88 def _raise_linalgerror_singular(err, flag):
---> 89     raise LinAlgError("Singular matrix")

LinAlgError: Singular matrix

My example truss is 2-dimensional, but I wanted to test bar3e whether it also works in 3 dimensions, and accordingly set all coordinates to zero on one axis.

Martin15135215 commented 10 months ago

I forgot to add a support:

Now, the code is working without error.

import numpy as np
import calfem.core as cfc

# Topology matrix Edof
Edof = np.array([
    [1,2,3,4,5,6], # Element 1
    [1,2,3,7,8,9], # Element 2
    [1,2,3,10,11,12], # Element 3
])

K = np.matrix(np.zeros((12,12)))
f = np.matrix(np.zeros((12,1)))

# Area and E-Modulus
E =  1000000.0 # E-Modulus kN/m²
A = 0.001 # Area m²
ep = [E,A]

#Element coordinates
ex1 = np.array([0., 0.])  # Element 1
ex2 = np.array([0., -4.]) # Element 2
ex3 = np.array([0., -4.]) # Element 3

ey1 = np.array([0., 0.]) # Element 1
ey2 = np.array([0., 0.]) # Element 2
ey3 = np.array([0., 0.]) # Element 3

ez1 = np.array([3., 0.])  # Element 1
ez2 = np.array([3., 3.])  # Element 2
ez3 = np.array([3., 6.])  # Element 3

# Element stiffness matrices
Ke1 = cfc.bar3e(ex1,ey1,ez1,ep)
Ke2 = cfc.bar3e(ex2,ey2,ez2,ep)
Ke3 = cfc.bar3e(ex3,ey3,ez3,ep)

# Assemble Ke into K
cfc.assem(Edof[0,:],K,Ke1)
cfc.assem(Edof[1,:],K,Ke2)
cfc.assem(Edof[2,:],K,Ke3)

# prescribed displacements in bc
bc = np.array([2, # has also a support
    4,5,6,7,8,9,10,11,12])

f[1] = 100 # kN
f[2] = -100 # kN

a, r = cfc.solveq(K,f,bc)