"""
A single bar problem

    0----1
         ^ Fy
         |

    ^ y, v
    |
    ----> x, u
    
    bar: stiffness E*A = 1
         length    L   = 1
    Load: Fy = 1
    Constraints(displacements: u,v): u_0 = 0
                                     v_0 = 0
"""
"""
It yields:

KU = F
[ 1.0   0.0   -1.0  0.0  ][  0.0   ]    [  0.0   ]
[ 0.0   0.5   0.0   -0.5 ][  0.0   ]    [  0.0   ]
[ -1.0  0.0   1.0   0.0  ][  0.0   ]    [  0.0   ]
[ 0.0   -0.5  0.0   0.5  ][  2.0   ]    [  1.0   ]

The bar has 1.0 stiffness in x direction.
It is because stiffness E*A = 1.
On the other hand, the bar has 0.5 stiffness in y direction.

Can we know why?
"""

import getfem as gf
import numpy as np

m = gf.Mesh('empty', 2)

# set points
pts = np.array([[0, 1],
                [0, 0]])
pids = m.add_point(pts)

# set convexes
gt = gf.GeoTrans('GT_PK(1, 1)')
cvids = m.add_convex(gt, pts)

# set regions
cvfids = m.faces_from_pid(pids)  # cvfids: [[0, 0], [0, 1]]
DIRICHLET_BOUNDARY = 1
m.set_region(DIRICHLET_BOUNDARY, np.array([[0], [1]]))
NEUMANN_BOUNDARY = 2
m.set_region(NEUMANN_BOUNDARY, np.array([[0], [0]]))

# set fem im 
mfu = gf.MeshFem(m, 2)  # 2: vector on xy-plane
mfu.set_fem(gf.Fem('FEM_PK(1, 1)'))  # 2-dof on one segment
mim = gf.MeshIm(m, gf.Integ('IM_GAUSS1D(1)'))

# handle DIRICHLET_BOUNDARY
mfd = gf.MeshFem(m, 1)  #1: single value
mfd.set_fem(gf.Fem('FEM_PK(1, 0)')) # 1-dof on one segment
Hs = mfd.eval('[[1, 0], [0, 1]]')
Rs = mfd.eval('[0, 0]')
H, R = gf.asm_dirichlet(DIRICHLET_BOUNDARY, mim, mfu, mfd, Hs, Rs)
N, U0 = H.dirichlet_nullspace(R)

# handle NEUMANN_BOUNDARY
Fs = np.repeat([[0.], [1.]], mfd.nbdof(), axis=1)  # Fx=0, Fy=400
F = gf.asm_boundary_source(NEUMANN_BOUNDARY, mim, mfu, mfd, Fs)

# get K: stiffness matrix
E = 1.
Nu = 0.  # no interaction of deformation in x- and y-direction
A = 1.
Lambda = A * E * Nu / ((1. + Nu) * (1. - 2. * Nu))
Mu = A * E / (2. * (1. + Nu))
K = gf.asm_linear_elasticity(mim, mfu, mfd, 
                             np.repeat([Lambda], mfd.nbdof()), 
                             np.repeat([Mu], mfd.nbdof()), -1)

# solve
Nt = gf.Spmat('copy', N)
Nt.transpose()
KK = Nt * K * N
FF = Nt * F
P = gf.Precond('ildlt', KK)
UU = gf.linsolve_cg(KK, FF, P)
U = N * UU + U0

# print KU = F
print('KU = F')
for i in range(len(U)):
    print('[', end='')
    for kj in K.full()[i]:
        print(str(kj).center(6), end='')
    print(']', end='')
    print('[', str(U[i]).center(6), ']    [', str(F[i]).center(6), ']')