```python
# Header starts here.
from sympy.physics.units import *
from sympy import *

# Rounding:
import decimal
from decimal import Decimal as DX
from copy import deepcopy
def iso_round(obj, pv, rounding=decimal.ROUND_HALF_EVEN):
    import sympy
    """
    Rounding acc. to DIN EN ISO 80000-1:2013-08
    place value = Rundestellenwert
    """
    assert pv in set([
        # place value   #  round to:
        1,              #  1
        0.1,            #  1st digit after decimal
        0.01,           #  2nd
        0.001,          #  3rd
        0.0001,         #  4th
        0.00001,        #  5th
        0.000001,       #  6th
        0.0000001,      #  7th
        0.00000001,     #  8th
        0.000000001,    #  9th
        0.0000000001,   # 10th
        ])
    objc = deepcopy(obj)
    try:
        tmp = DX(str(float(objc)))
        objc = tmp.quantize(DX(str(pv)), rounding=rounding)
    except:
        for i in range(len(objc)):
            tmp = DX(str(float(objc[i])))
            objc[i] = tmp.quantize(DX(str(pv)), rounding=rounding)
    return objc

# LateX:
kwargs = {}
kwargs["mat_str"] = "bmatrix"
kwargs["mat_delim"] = ""

# kwargs["symbol_names"] = {FB: "F^{\mathsf B}", }

# Units:
(k, M, G ) = ( 10**3, 10**6, 10**9 )
(mm, cm) = ( m/1000, m/100 )
Newton = kg*m/s**2
Pa     = Newton/m**2
MPa    = M*Pa
GPa    = G*Pa
kN     = k*Newton
deg    = pi/180

half = S(1)/2

# Header ends here.
#
a, F, EA = var("a, F, EA")
(l1, l2, l3) = (a, sqrt(2)*a, a)

S1 = - F
S2 = sqrt(2)*F
S3 = - F

(dl1, dl2, dl3) = ( l1/EA*S1, l2/EA*S2, l3/EA*S3 )

s22 = sqrt(2)/2
e1 = Matrix([1, 0])
e2 = Matrix([-s22, s22])
e3 = Matrix([0, 1])

u2x, u2y, u3y = var("u2x, u2y, u3y")
u1 = Matrix([0,0])
u2 = Matrix([u2x,u2y])
u3 = Matrix([0,u3y])

eq1 = Eq(dl1, e1.dot(u2-u1))
eq2 = Eq(dl2, e2.dot(u3-u2))
eq3 = Eq(dl3, e3.dot(u3-u1))

sol = solve([eq1, eq2, eq3],[u2x, u2y, u3y])

(u2x, u2y, u3y) = (sol[u2x], sol[u2y], sol[u3y])

pprint("\nClassical Solution:")
pprint("u2x, u2y, u3y / ( Fa / EA ):")
for u in (u2x, u2y, u3y):
    u*=EA/F/a
    pprint(u)

# FEM-Solution:
EA, a = var("EA, a")
# psi:
p = sqrt(2)/4

u2x, u2y, u3y = var("u2x, u2y, u3y")

K = Matrix((
    [p+1, -p,  p ],
    [ -p,  p, -p ],
    [  p, -p, p+1]
    ))
K*=EA/a
u = Matrix([u2x, u2y, u3y])
f = Matrix([0,-F,0])
eq = Eq( K*u, f )
sol = solve(eq,[u2x, u2y, u3y])
(u2x, u2y, u3y) = (sol[u2x], sol[u2y], sol[u3y])

pprint("\nFEM Solution:")
pprint("u2x, u2y, u3y / (Fa/EA):")
for u in (u2x, u2y, u3y):
    u*=EA/F/a
    pprint(u)

F1x, F1y, F3x = var("F1x, F1y, F3x")

f = EA/a
eq1 = Eq( f * (-u2x), F1x)
eq2 = Eq( f * (-u3y), F1y)
eq3 = Eq( f * p * ( -u2x + u2y - u3y ), F3x)

sol = solve([eq1, eq2, eq3], [F1x, F1y, F3x])

pprint("\n")
pprint(sol)

pprint("\nS1, S2, S3:")
f = sqrt(2)/2
dl1 = u2x
dl2 = f * (u2x + u3y - u2y)
dl3 = u3y
S1 = EA *dl1/l1
S2 = EA *dl2/l2
S3 = EA *dl3/l3
for S in [S1, S2, S3]:
    S = S.simplify()
    pprint(S)

# Classical Solution:
# u2x, u2y, u3y / ( Fa / EA ):
# -1
# -2⋅√2 - 2
# -1
#
# FEM Solution:
# u2x, u2y, u3y / (Fa/EA):
# -1
# -2⋅√2 - 2
# -1
#
# {F1x: F, F1y: F, F3x: -F}
#
# S1, S2, S3:
# -F
# √2⋅F
# -F

```

1905a5605f82a0cd79572597f578394d4a62ff6a
ipynb
Jupyter Notebook
ipynb/TM_2/1_ZD/2.1.N-FEM/2.1.N-FEM_cc.ipynb
kassbohm/tm-snippets
5e0621ba2470116e54643b740d1b68b9f28bff12
MIT
ipynb/TM_2/1_ZD/2.1.N-FEM/2.1.N-FEM_cc.ipynb
kassbohm/tm-snippets
5e0621ba2470116e54643b740d1b68b9f28bff12
MIT
ipynb/TM_2/1_ZD/2.1.N-FEM/2.1.N-FEM_cc.ipynb
kassbohm/tm-snippets
5e0621ba2470116e54643b740d1b68b9f28bff12
MIT
Qwen/Qwen-72B
1. YES
2. YES