- Thank you received: 0
Forum Message
We have moved the forum to https://forum.ngsolve.org . This is an archived version of the topics until 05/05/23. All the topics were moved to the new forum and conversations can be continued there. This forum is just kept as legacy to not invalidate old links. If you want to continue a conversation just look for the topic in the new forum.
Notice
The forum is in read only mode.
Solving Hyperleasticity: Correctly defining boundaries
- bhaveshshrimali
- Topic Author
- Offline
- Junior Member
-
Less
More
3 years 5 months ago #3324
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.
gives
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
-
Less
More
- Thank you received: 101
3 years 5 months ago - 3 years 5 months 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:
Best
Christopher
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"))
Christopher
Last edit: 3 years 5 months ago by christopher.
- bhaveshshrimali
- Topic Author
- Offline
- Junior Member
-
Less
More
- Thank you received: 0
3 years 5 months 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...
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
-
Less
More
- Thank you received: 0
3 years 5 months 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.
I will take a look at implementing non-homogenous BC's more carefully.
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
-
Less
More
- Thank you received: 0
3 years 5 months 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:
gives
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.
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
-
Less
More
- Thank you received: 59
3 years 5 months 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:
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:
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
which will provide a better initial guess helping Newton to converge.
Best,
Michael
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)
Best,
Michael
The following user(s) said Thank You: bhaveshshrimali
Time to create page: 0.113 seconds
