Solving Hyperleasticity: Correctly defining boundaries

bhaveshshrimali
Topic Author
Offline
Junior Member

4 years 1 week ago #3324 by bhaveshshrimali

Solving Hyperleasticity: Correctly defining boundaries was created by bhaveshshrimali

Hi everyone,

I was trying to translate the hyperelasticity example from FEniCS to ngsolve, following the nonlinear elasticity tutorial, and Navier Stokes tutorial to implement homogenous and non-homogenous dirichlet boundary conditions on different parts of the boundary.

However, it seems that the boundary conditions are not getting assigned properly (looking at the GUI). Would anyone be able to help me pin down the error? You can see from the norm of the error, that it does not converge.

Code:

import netgen.gui
from math import pi as myPi
from ngsolve import *
from netgen.csg import unit_cube

mesh = Mesh(unit_cube.GenerateMesh(maxh=0.1))
V = VectorH1(mesh, order=1, dirichlet="left|right")
u = GridFunction(V, str="u")
u.vec[:] = 0.
du, v = V.TnT()


# material parameters and forces
E, nu = 10.0, 0.3
mu = E / 2./(1+nu)
lmbda = 2 * mu * nu/(1-2*nu)
bodyForce = CoefficientFunction((0, -0.5, 0.))

# definition of the boundaries 
scale = 0.5
yo, zo = 0.5, 0.5
thta = myPi/3.
uLeft = CoefficientFunction((0., 0., 0.))
uRight = CoefficientFunction((
    0.,
    scale * (yo + (y - yo) * cos(thta) - (z - zo) * sin(thta) - y),
    scale * (zo + (y - yo) * sin(thta) + (z-zo) * cos(thta) - z)
))
u.Set(uLeft, definedon=mesh.Boundaries("left"))
u.Set(uRight, definedon=mesh.Boundaries("right"))

# Define strain measures
I = Id(mesh.dim)
F = I + grad(du)
I1 = Trace(F.trans * F)
J = Det(F)
psi = mu * (I1 - 3) - mu * log(J) + lmbda/2. * (J-1)**2

# definition of bilinear and linear forms
a = BilinearForm(V, symmetric=False)
a += Variation(psi.Compile() * dx)
a += Variation((-1 * InnerProduct(bodyForce, du)).Compile() * dx)

# create the residual and iteration vectors
res = u.vec.CreateVector()
w = u.vec.CreateVector()

resNorm = 1.0

while resNorm > 1.e-8:
    a.Apply(u.vec, res)
    a.AssembleLinearization(u.vec)
    inv = a.mat.Inverse(V.FreeDofs())
    w.data = inv * res
    resNorm = InnerProduct(w, res)
    print(f"Norm: {resNorm}")
    u.vec.data -= w

    Draw(u, mesh, "displacement")
    SetVisualization(deformation=True)

gives

Code:

Norm: 6.686688205442815
Norm: 0.49337768606525656
Norm: 0.34295300814873925
Norm: 0.15001895765763176
Norm: 0.06517533253371934
Norm: 0.019404532237654374
Norm: 0.18895838479738908
Norm: 2.024789522475844
Norm: 5815.697445306668
Norm: 49695.17791867306
Norm: 58282.3479021152
Norm: 11584203088583.24
Norm: 2.3448080388890567e+22
Norm: 3.154397165784633e+26
Norm: 8.26906306359288e+25
Norm: 2.167685693089545e+25
Norm: 5.682457983539969e+24

christopher
Offline
Administrator

4 years 1 week ago - 4 years 1 week ago #3325 by christopher

Replied by christopher on topic Solving Hyperleasticity: Correctly defining boundaries

u.Set(cf, definedon=...)
sets u to cf on the definedon part and to 0 everywhere else, set the boundary values like this:

Code:

u.Set(mesh.BoundaryCF({"left" : uLeft, "right" : uRight}),  definedon=mesh.Boundaries("left|right"))

Best
Christopher

Last edit: 4 years 1 week ago by christopher.

bhaveshshrimali
Topic Author
Offline
Junior Member

4 years 1 week ago #3328 by bhaveshshrimali

Replied by bhaveshshrimali on topic Solving Hyperleasticity: Correctly defining boundaries

Thanks Christopher!

Is this a newer feature? NGSolve-6.2.1909 does not seem to have it. I will try updating NGSolve to the latest version and test it again...

Code:

AttributeError: 'ngsolve.comp.Mesh' object has no attribute 'BoundaryCF'

bhaveshshrimali
Topic Author
Offline
Junior Member

4 years 1 week ago #3329 by bhaveshshrimali

Replied by bhaveshshrimali on topic Solving Hyperleasticity: Correctly defining boundaries

Ok this works with the newer version. But even with the newer version (NGSolve-6.2.2009), I get a diverging Newton.

Code:

Norm: 2.6808208867801797
Norm: 577.2816648847577
Norm: 152009.87301735353
Norm: 7366462.124902223
Norm: 65416518034032.02
Norm: 5.554793589592823e+16
Norm: 1.4678499972939826e+16
Norm: 3848066201629095.0
Norm: 1008756719615877.0
Norm: 264477504711945.2

I will take a look at implementing non-homogenous BC's more carefully.

bhaveshshrimali
Topic Author
Offline
Junior Member

4 years 1 week ago #3331 by bhaveshshrimali

Replied by bhaveshshrimali on topic Solving Hyperleasticity: Correctly defining boundaries

After inspecting a bit more the location of the boundaries "left", "right", ... I noticed that the boundaries that I need (x=0 and x=1) correspond to "back" and "front". Indeed after replacing them accordingly and applying the boundary condition in steps, rather than at once, the solution seems to be reasonable till a point, and then it diverges:

Code:

for i in range(1, 11):
    factor.Set(i/10)
    u.Set(mesh.BoundaryCF({"back" : uLeft, "front" : uRight}),  definedon=mesh.Boundaries("back|front"))
    print(i)
    resNorm = 1.0
    iter = 0    
    while resNorm > 1.e-8 and iter < maxiters:
        a.Apply(u.vec, res)
        a.AssembleLinearization(u.vec)
        inv = a.mat.Inverse(V.FreeDofs())
        w.data = inv * res
        resNorm = sqrt(abs(InnerProduct(w, w)))
        print(f"Norm: {resNorm}")
        u.vec.data -= w
        iter += 1
        Draw(u, mesh, "displacement")
        SetVisualization(deformation=True)

gives

Code:

1
Norm: 0.840274935028221
Norm: 0.010631338197297788
Norm: 2.840043781994656e-05
Norm: 5.625360578002826e-10
2
Norm: 1.3145809712371617
Norm: 0.03864670201207834
Norm: 0.0003690373463827511
Norm: 8.934665769419485e-08
Norm: 1.54388135479357e-14
3
Norm: 1.8523594043744418
Norm: 0.08405462602954474
Norm: 0.0018541634265310004
Norm: 1.9279392913312044e-06
Norm: 8.193550778213933e-12
4
Norm: 2.4123885921399872
Norm: 0.14726391288538263
Norm: 0.0060619644879045765
Norm: 1.9854480113520826e-05
Norm: 1.1251011825818185e-09
5
Norm: 2.816514170226768
Norm: 0.4758279175877349
Norm: 0.07470240485416736
Norm: 0.14433503658998675
Norm: 0.2576930420365274
Norm: 0.15947077282736843
Norm: 0.04193568718574943
Norm: 0.04151801466555532
Norm: 0.0044976089152604015
Norm: 5.306089515763066e-05
Norm: 8.324583871537812e-09
6
Norm: 4.623385880239417
Norm: 2.913358091896825
Norm: 2.617097160078447
Norm: 2.9116418166392988
Norm: 6.296120155081865
Norm: 12.355324294951478
Norm: 19.201276563870323
Norm: 58.62691253475249
Norm: 230.766474679662

I think, this is due to inaccurate approximation of the jacobian in automatic differentiation of the residual. I don't see any other reason on why Newton wouldn't converge after a point for this problem. I checked the deformation and it seems reasonable till the 6th loading step.

mneunteufel
Offline
Premium Member

4 years 1 week ago #3332 by mneunteufel

Replied by mneunteufel on topic Solving Hyperleasticity: Correctly defining boundaries

Hi bhaveshshrimali,

are you sure your material law is correct? For zero external forces and zero Dirichlet boundary condition the identity function is not the trivial solution. For the Neo-Hookean material law mu/2 should be the correct scaling for the first term:

Code:

psi = mu/2 * (I1 - 3) - mu * log(J) + lmbda/2. * (J-1)**2

By setting the boundary conditions for every loadstep you reset your solution rather than using the solution of the previous loadstep as starting point for the next one. I would suggest to set the boundary conditions only at the beginning:

Code:

u.Set(mesh.BoundaryCF({"back" : uLeft, "front" : uRight}),  definedon=mesh.Boundaries("back|front"))
for i in range(1, 11):
    factor.Set(i/10)
    ...

The boundary condition at the right boundary is quite large. As on the opposite boundary you have zero Dirichlet condition you can set them by

Code:

u.Set(x*uRight)

which will provide a better initial guess helping Newton to converge.

Best,
Michael

The following user(s) said Thank You: bhaveshshrimali

Time to create page: 0.108 seconds