This page was generated from wta/elasticity.ipynb.

Elasticity

import netgen.geom2d as geom2d
from ngsolve import *
from ngsolve.webgui import Draw

a rectangle with refinement at corners:

geo = geom2d.SplineGeometry()
pnums = [geo.AddPoint (x,y,maxh=0.01) for x,y in [(0,0), (1,0), (1,0.1), (0,0.1)] ]
for p1,p2,bc in [(0,1,"bot"), (1,2,"right"), (2,3,"top"), (3,0,"left")]:
     geo.Append(["line", pnums[p1], pnums[p2]], bc=bc)
mesh = Mesh(geo.GenerateMesh(maxh=0.05))

 Generate Mesh from spline geometry
 Boundary mesh done, np = 68
 CalcLocalH: 68 Points 0 Elements 0 Surface Elements
 Meshing domain 1 / 1
 load internal triangle rules
 Surface meshing done
 Edgeswapping, topological
 Smoothing
 Split improve
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Split improve
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Split improve
 Combine improve
 Smoothing
 Update mesh topology
 Update clusters

Cauchy-Green tensor and hyperelastic energy density:

E, nu = 210, 0.2
mu  = E / 2 / (1+nu)
lam = E * nu / ((1+nu)*(1-2*nu))

def C(u):
    F = Id(2) + Grad(u)
    return F.trans * F

def NeoHooke (C):
    return 0.5*mu*(Trace(C-Id(2)) + 2*mu/lam*Det(C)**(-lam/2/mu)-1)

stationary point of total energy:

\[\delta \int W(C(u)) - f u = 0\]

factor = Parameter(0)
force = CoefficientFunction( (0,factor) )

fes = H1(mesh, order=4, dirichlet="left", dim=mesh.dim)
u  = fes.TrialFunction()

a = BilinearForm(fes, symmetric=True)
a += Variation(NeoHooke(C(u)).Compile()*dx)
a += Variation((-InnerProduct(force,u)).Compile()*dx)

u = GridFunction(fes)
u.vec[:] = 0

Compiled CF:
N5ngfem27IdentityCoefficientFunctionE
N5ngfem13ProxyFunctionE
N5ngfem13cl_BinaryOpCFINS_11GenericPlusEEE
N5ngfem28TransposeCoefficientFunctionE
N5ngfem29MultMatMatCoefficientFunctionE
N5ngfem27IdentityCoefficientFunctionE
N5ngfem13cl_BinaryOpCFINS_12GenericMinusEEE
N5ngfem24TraceCoefficientFunctionE
N5ngfem30DeterminantCoefficientFunctionILi2EEE
N5ngfem27ConstantCoefficientFunctionE
N5ngfem13cl_BinaryOpCFI10GenericPowEE
N5ngfem24ScaleCoefficientFunctionE
N5ngfem13cl_BinaryOpCFINS_11GenericPlusEEE
N5ngfem27ConstantCoefficientFunctionE
N5ngfem13cl_BinaryOpCFINS_12GenericMinusEEE
N5ngfem24ScaleCoefficientFunctionE
inputs =
0:
1:
2: 0 1
3: 2
4: 3 2
5:
6: 4 5
7: 6
8: 4
9:
10: 8 9
11: 10
12: 7 11
13:
14: 12 13
15: 14

Compiled CF:
N5ngfem27ConstantCoefficientFunctionE
N5ngfem28ParameterCoefficientFunctionIdEE
N5ngfem28VectorialCoefficientFunctionE
N5ngfem12cl_UnaryOpCFI15GenericIdentityEE
N5ngfem13ProxyFunctionE
N5ngfem31T_MultVecVecCoefficientFunctionILi2EEE
N5ngfem24ScaleCoefficientFunctionE
inputs =
0:
1:
2: 0 1
3: 2
4:
5: 3 4
6: 5

a simple Newton solver, using automatic differentiation for residual and tangential stiffness:

def SolveNewton():
    res = u.vec.CreateVector()

    for it in range(10):
        print ("it", it, "energy = ", a.Energy(u.vec))
        a.Apply(u.vec, res)
        a.AssembleLinearization(u.vec)
        inv = a.mat.Inverse(fes.FreeDofs() )
        u.vec.data -= inv*res

factor.Set(0.4)
SolveNewton()
scene = Draw (C(u)[0,0], mesh, deformation=u)

it 0 energy =  8.749999999999972
Assemble linearization
assemble VOL element 170/170
call pardiso ... done
it 1 energy =  8.811178274176932
Assemble linearization
assemble VOL element 170/170
call pardiso ... done
it 2 energy =  8.74811675282931
Assemble linearization
assemble VOL element 170/170
call pardiso ... done
it 3 energy =  8.747829195439456
Assemble linearization
assemble VOL element 170/170
call pardiso ... done
it 4 energy =  8.747829114564311
Assemble linearization
assemble VOL element 170/170
call pardiso ... done
it 5 energy =  8.747829114549656
Assemble linearization
assemble VOL element 170/170
call pardiso ... done
it 6 energy =  8.74782911454966
Assemble linearization
assemble VOL element 170/170
call pardiso ... done
it 7 energy =  8.74782911454966
Assemble linearization
assemble VOL element 170/170
call pardiso ... done
it 8 energy =  8.74782911454966
Assemble linearization
assemble VOL element 170/170
call pardiso ... done
it 9 energy =  8.74782911454966
Assemble linearization
assemble VOL element 170/170
call pardiso ... done

factor.Set(factor.Get()+0.4)
SolveNewton()
scene.Redraw()

it 0 energy =  8.743591354967782
Assemble linearization
assemble VOL element 170/170
call pardiso ... done
it 1 energy =  8.781886826184756
Assemble linearization
assemble VOL element 170/170
call pardiso ... done
it 2 energy =  8.741979230980744
Assemble linearization
assemble VOL element 170/170
call pardiso ... done
it 3 energy =  8.74188819248704
Assemble linearization
assemble VOL element 170/170
call pardiso ... done
it 4 energy =  8.741861474021357
Assemble linearization
assemble VOL element 170/170
call pardiso ... done
it 5 energy =  8.74185078504465
Assemble linearization
assemble VOL element 170/170
call pardiso ... done
it 6 energy =  8.741850620429458
Assemble linearization
assemble VOL element 170/170
call pardiso ... done
it 7 energy =  8.741850620044117
Assemble linearization
assemble VOL element 170/170
call pardiso ... done
it 8 energy =  8.74185062004412
Assemble linearization
assemble VOL element 170/170
call pardiso ... done
it 9 energy =  8.74185062004412
Assemble linearization
assemble VOL element 170/170
call pardiso ... done

Compute \(2^{nd}\) Piola Kirchhoff stress tensor by symbolic differentiation:

C_=C(u)
sigma = NeoHooke(C_).Diff(C_)

Draw (sigma[0,0], mesh, "Sxx", deformation=u, min=-10.001, max=10.001)

BaseWebGuiScene