Indicial notation for constitutive laws in finite elasticity

[bhaveshshrimali]

Hi all,

Finally started experimenting with Gridap after a long overdue :)

I had a very basic question What is the best way in Gridap to include indicial notation as in continuum mechanics. For instance, in the following I would like to calculate the tangent modulus L = dS/dF which could be given component wise as mu * one[i,k]*one[j,l] + ...

# deformation gradient
F(∇u) = one(∇u) + ∇u
J = det(F)

# constitutive relation
@law function S(∇u)
    Js = J(F(∇u))
    μ * (F - (inv(F))') + λ * Js * (Js-1) * (inv(F(∇u)))' 
end

@law function L(∇u)

end

Albeit for this simple case it may not be needed, the idea is to do this for more general constitutive relations, which could be cumbersom to write abstractly

[santiagobadia]

Probably, the easiest way to implement this is defining the law after applying the direction, e.g.,

@law function L(∇u,∇du)
  SymmetricTensorValue(your components)
end

You need a SymmetricTensorValue, because you will use L(∇u,∇du)⊙ ε(dv) in your bilinear form.

Alternatively, you could use a SymmetricFourthOrderTensor for L(∇u) but probably more complicated.

[santiagobadia]

By the way, even though not in the tests yet, you can also use autodiff. Just define the residual and the tangent operator is automatically calculated.

Take a look here:

https://github.com/gridap/Gridap.jl/blob/ce0e1f1dd24aca6fac88fa699047393faac753ba/test/GridapTests/PLaplacianWithAutodiffTests.jl#L34

You just do FETerm(res,trian,quad) and the Jacobian is computer internally

[bhaveshshrimali]

You need a SymmetricTensorValue, because you will use L(∇u,∇du)⊙ ε(dv) in your bilinear form.

Thanks @santiagobadia !! Right, I was thinking more like L(∇u,∇du)⊙ ∇v but in any case the tangent (right now) will have major symmetry so either should work.

You just do FETerm(res,trian,quad) and the Jacobian is computer internally

This is fantastic. I will definitely check it out.

On a slightly related note, I am trying to implement a simple 3D hyperelasticity example for understanding Gridap usage. In fact it is the FEniCS hyperelasticity example. I wanted to check the numbering of the facet IDs to apply Dirichlet boundary conditions, but somehow I am not able to write the model to a vtk file.

using Gridap
domain = (0, 1, 0, 1, 0, 1)
partition = (24, 16, 16)
model = CartesianDiscreteModel(domain, partition)
writevtk(model, "model")

fails with

ERROR: LoadError: MethodError: Cannot `convert` an object of type
  WriteVTK.MeshCell{WriteVTK.VTKCellTypes.VTKCellType,Array{Int64,1}} to an object of type
  WriteVTK.MeshCell{Array{Int64,1},V} where V<:(Union{AbstractArray{T,1}, Tuple{Vararg{T,N}} where N} where T<:Integer)
Closest candidates are:
  convert(::Type{T}, !Matched::T) where T at essentials.jl:171

Any clues as to what could be going wrong ?

[fverdugo]

Hi @bhaveshshrimali

This looks like an old error (issue #286). It should be fixed. Which version are you using?

[bhaveshshrimali]

This looks like an old error (issue #286). It should be fixed. Which version are you using?

Ah. I am using 0.9.2

[fverdugo]

This is way too old. Use the latest one v0.14.0

[bhaveshshrimali]

This is way too old. Use the latest one v0.14.0

In trying to update Gridap to the master branch

(@v1.4) pkg>add Gridap#master

it seems to be incompatible with GridapPardiso and GridapODEs for some reasons. Removing those in the meantime still gave errors with DifferentialEquations.jl. So I added the v0.13.4 version. Note that doing

(@v1.4) pkg> add Gridap@0.14.0

threw a bunch of compatibility errors

ERROR: Unsatisfiable requirements detected for package Parameters [d96e819e]:
 Parameters [d96e819e] log:
 ├─possible versions are: [0.9.1-0.9.2, 0.10.0-0.10.3, 0.11.0, 0.12.0-0.12.1] or uninstalled
 ├─restricted by compatibility requirements with LineSearches [d3d80556] to versions: [0.9.1-0.9.2, 0.10.0-0.10.3, 0.11.0, 0.12.0-0.12.1]
 │ └─LineSearches [d3d80556] log:
 │   ├─possible versions are: [7.0.0-7.0.1, 7.1.0] or uninstalled
 │   ├─restricted to versions * by an explicit requirement, leaving only versions [7.0.0-7.0.1, 7.1.0]
 │   └─restricted by compatibility requirements with Gridap [56d4f2e9] to versions: [7.0.1, 7.1.0]

Edit: I need to clean up my environment it seems.

[bhaveshshrimali]

Making some progress... I am trying to enforce non-homogeneous Dirichlet data on boundary. I essentially followed the syntax from the Poisson tutorial

using Gridap
domain = (0, 1, 0, 1, 0, 1)
partition = (5, 5, 5)
model = CartesianDiscreteModel(domain, partition)
labels = get_face_labeling(model)

# create meshes in Gmsh and use GridapGmsh, but this should do meanwhile?

add_tag_from_tags!(labels,"rightFace",[26])
add_tag_from_tags!(labels,"leftFace",[25])

V = TestFESpace(
  model=model, valuetype=VectorValue{2,Float64},
  reffe=:Lagrangian,conformity=:H1,order=1,
  dirichlet_tags = ["leftFace", "rightFace"])

g0(x) = VectorValue(0.0,0.0,0.0)
g1(x) = VectorValue(
    0.0,
    (0.5 + (x[2]-0.5) * cos(π/3) - (x[3] - 0.5) * sin(π/3) - x[2])/2,
    (0.5 + (x[2]-0.5) * sin(π/3) + (x[3] - 0.5) * cos(π/3) - x[3])/2
)

U = TrialFESpace(V, [g0, g1])

But this throws up BoundsError

ERROR: LoadError: BoundsError: attempt to access (1, 9)
  at index [3]

when trying to create the TrialFESpace. I can post the entire traceback if you want. Sorry for bugging with trivial queries. I'll pick things up in time.

Thanks so much for your help in advance!

Edit: I pasted things from the hyperelasticity tutorial and forgot to change the definition of V to

V = TestFESpace(
  model=model,valuetype=VectorValue{3,Float64},
  reffe=:Lagrangian,conformity=:H1,order=1,
  dirichlet_tags = ["leftFace", "rightFace"])

[bhaveshshrimali]

So after a couple of more tries, I managed to get a working code without writing the jacobian explicitly, and letting it calculate the jacobian using AD, as Santiago suggested above:

using Gridap
using LinearAlgebra: inv, det
using LineSearches: BackTracking

E = 10.
ν = 0.3
const μ = E/2/(1+ν)
const λ = 2*μ*ν/(1-2*ν)

# deformation gradient

F(∇u) = one(∇u) + ∇u
J(F) = det(F)
C(F) = (F') * F

# body force
B(x) = VectorValue(
    0.0, -0.5, 0.0
)

# surfaceTraction
T(x) = VectorValue(
    0.1, 0.0, 0.0
)

# constitutive relation
@law function S(∇u)
    Js = J(F(∇u))
    μ * (F(∇u) - (inv(F(∇u)))') + λ * Js * (Js - 1) * (inv(F(∇u)))' 
end

function res(u, v)
    S(∇(u)) ⊙ ∇(v) - B ⊙ v
end

domain = (0, 1, 0, 1, 0, 1)
partition = (24, 16, 16)
model = CartesianDiscreteModel(domain, partition)
labels = get_face_labeling(model)
add_tag_from_tags!(labels,"rightFace",[26])
add_tag_from_tags!(labels,"leftFace",[25])
add_tag_from_tags!(labels,"topFace",[24])

V = TestFESpace(
  model=model,valuetype=VectorValue{3,Float64},
  reffe=:Lagrangian,conformity=:H1,order=1,
  dirichlet_tags = ["leftFace", "rightFace"])

trian = Triangulation(model)
degree=2
quad=CellQuadrature(trian, degree)

neumannTags = ["topFace"]
btrian = BoundaryTriangulation(model, neumannTags)
bquad = CellQuadrature(btrian, degree)

a_Ω = FETerm(res, trian, quad)

b_Γ(v) = T ⊙ v
t_Γ = FESource(b_Γ, btrian, bquad)

nls = NLSolver(
  show_trace=true,
  method=:newton,
  linesearch=BackTracking()
)

g0(x) = VectorValue(0.0,0.0,0.0)
g1(x) = VectorValue(
    0.0,
    (0.5 + (x[2]-0.5) * cos(π/3) - (x[3] - 0.5) * sin(π/3) - x[2])/2,
    (0.5 + (x[2]-0.5) * sin(π/3) + (x[3] - 0.5) * cos(π/3) - x[3])/2
)

U = TrialFESpace(V, [g0, g1])

#FE problem
op = FEOperator(U, V, a_Ω, t_Γ)
solver = FESolver(nls)

x0 = zeros(Float64, num_free_dofs(V))
uh = FEFunction(U, x0)
uh, = solve!(uh, solver, op)

writevtk(trian,"results_hyper",cellfields=["uh"=>uh,"sigma"=>S(∇(uh))])

however the results seem to be a bit off (I will check the Dirichlet boundary conditions). I am guessing either the facets aren't properly identified (which of course you could tell better) or there is something else that I am doing incorrectly.

Result

Convergence

Albeit poor in this case, I would be happy to know on how to tune knobs on the solver to accelerate this. :)

Iter     f(x) inf-norm    Step 2-norm
------   --------------   --------------
     0     3.112098e-01              NaN
     1     3.438451e-01     9.883077e-01
     2     3.101497e-01     8.094973e-01
     3     3.081619e-01     8.202590e-03
     4     2.584847e-01     7.362241e-01
     5     2.561357e-01     7.525048e-03
     6     2.443304e-01     2.779464e-01
     7     2.333284e-01     1.297588e-01
     8     2.304624e-01     1.591403e-02
     9     2.112249e-01     4.549931e-01
    10     2.015671e-01     8.803574e-02
    11     1.847135e-01     2.841682e-01
    12     1.760663e-01     6.510608e-02
    13     1.574417e-01     3.385638e-01
    14     1.411452e-01     1.721633e+01
    15     6.264860e-02     1.482406e-02
    16     6.257339e-02     1.374721e-04
    17     5.961920e-02     2.429494e-02
    18     5.428662e-02     1.533825e-02
    19     5.395815e-02     9.677739e-04
    20     4.525693e-02     3.074569e-02
    21     4.462859e-02     1.753714e-02
    22     4.016006e-02     3.056113e-03
    23     3.290873e-02     7.791026e-03
    24     2.964262e-02     1.851017e-03
    25     1.798754e-02     1.967587e-02
    26     1.453088e-02     2.984918e-03
    27     1.556418e-02     6.476555e-02
    28     1.401099e-02     3.537729e-05
    29     1.116159e-02     1.143109e-04
    30     6.233009e-03     1.672470e-03
    31     2.843089e-03     5.852692e-04
    32     3.018364e-04     2.283751e-05
    33     4.181298e-05     5.906123e-06
    34     7.248264e-08     1.036177e-08
    35     6.395863e-12     8.371601e-13

Thanks again!

[bhaveshshrimali]

Just for completeness, when I try to define the tangent as

function Lv(∇u, ∇du)
    Js = J(F(∇u))
    Finv = inv(F(∇u))
    μ * (∇du + Finv' ⋅ ∇du' ⋅ Finv') + λ * (2 * Js - 1) * Js * ((Finv' ⊙ ∇du) * Finv') - λ * Js * (Js - 1) * (Finv' ⋅ ∇du' ⋅ Finv') 
end

jac(u, du, v) = Lv(∇(u), ∇(du)) ⊙ ∇(v)

and change the FETerm to

a_Ω = FETerm(res, jac, trian, quad)

it throws a method error on one which seems really weird. Am I missing anything obvious ?

ERROR: LoadError: MethodError: no method matching one(::Gridap.Geometry.GenericCellField{true})
Closest candidates are:
  one(!Matched::Type{Missing}) at missing.jl:103
  one(!Matched::BitArray{2}) at bitarray.jl:422
  one(!Matched::Missing) at missing.jl:100

Edit:

@law function Lv(∇du, ∇u)
    Js = J(F(∇u))
    Finv = inv(F(∇u))
    μ * (∇du + Finv' ⋅ ∇du' ⋅ Finv') + λ * (2 * Js - 1) * Js * ((Finv' ⊙ ∇du) * Finv') - λ * Js * (Js - 1) * (Finv' ⋅ ∇du' ⋅ Finv') 
end

fixes it, but I might have to check the expressions once again. If you happen to notice something in the application of Dirichlet boundary conditions, I'd be glad to know. No hurries :)

[fverdugo]

@bhaveshshrimali you are imposing Dirichlet BCs only at the interior of the faces, you also need to impose on the edges and vertices of their boundary. Just add the corresponding ids when defining "rightFace" etc.

[bhaveshshrimali]

Thanks a lot @fverdugo ! That worked perfectly! Is there a way to look for topological entities/subdomains in the mesh from inside Gridap? For instance, generating the tags using a bool function?

[fverdugo]

Is there a way to look for topological entities/subdomains in the mesh from inside Gridap? For instance, generating the tags using a bool function?

Issue https://github.com/gridap/Gridap.jl/issues/367

[bhaveshshrimali]

Here is the complete code in case if someone ends up wanting to take a look in the future. The timings are really good (16s) even when compared to FEniCS (12s) on my laptop, which wasn't the point of this exercise but that it ends up being as such is a bonus

Edit: I made a mistake in noting the timings. Gridap is about twice as fast (not exactly the same mesh) 12.916s (Gridap) vs 24.9 s (FEniCS)

Code:

using Gridap
using LinearAlgebra: inv, det
using LineSearches: BackTracking, StrongWolfe, HagerZhang
using BenchmarkTools

E = 10.
ν = 0.3
const μ = E/2/(1+ν)
const λ = 2*μ*ν/(1-2*ν)

function run()

    # deformation gradient

    F(∇u) = one(∇u) + ∇u
    J(F) = det(F)
    C(F) = (F') * F

    # body force
    B(x) = VectorValue(
        0.0, -0.5, 0.0
    )

    # surfaceTraction
    T(x) = VectorValue(
        0.1, 0.0, 0.0
    )

    # constitutive relation
    @law function S(∇u)
        Js = J(F(∇u))
        μ * (F(∇u) - (inv(F(∇u)))') + λ * Js * (Js - 1) * (inv(F(∇u)))' 
    end

    function res(u, v)
        S(∇(u)) ⊙ ∇(v) - B ⊙ v
    end

    @law function Lv(∇du, ∇u)
        Js = J(F(∇u))
        Finv = inv(F(∇u))
        μ * (∇du + Finv' ⋅ ∇du' ⋅ Finv') + λ * (2. * Js - 1.) * Js * ((Finv' ⊙ ∇du) * Finv') - λ * Js * (Js - 1) * (Finv' ⋅ ∇du' ⋅ Finv') 
    end

    jac(u, du, v) = Lv(∇(du), ∇(u)) ⊙ ∇(v)

    domain = (0, 1, 0, 1, 0, 1)
    partition = (24, 16, 16)
    model = CartesianDiscreteModel(domain, partition)
    labels = get_face_labeling(model)
    add_tag_from_tags!(labels,"rightFace",[2, 4, 6, 8, 14, 16, 18, 20, 26])
    add_tag_from_tags!(labels,"leftFace",[1, 3, 5, 7, 13, 15, 17, 19, 25])
    # add_tag_from_tags!(labels,"topFace",[4, 24])

    V = TestFESpace(
    model=model,valuetype=VectorValue{3,Float64},
    reffe=:Lagrangian,conformity=:H1,order=1,
    dirichlet_tags = ["leftFace", "rightFace"])

    trian = Triangulation(model)
    degree=2
    quad=CellQuadrature(trian, degree)

    neumannTags = ["boundary"]
    btrian = BoundaryTriangulation(model, neumannTags)
    bquad = CellQuadrature(btrian, degree)

    a_Ω = FETerm(res, trian, quad)

    b_Γ(v) = T ⊙ v
    t_Γ = FESource(b_Γ, btrian, bquad)

    nls = NLSolver(
    show_trace=true,
    method=:newton,
    linesearch=BackTracking()
    )

    g0(x) = VectorValue(0.0,0.0,0.0)
    g1(x) = VectorValue(
        0.0,
        (0.5 + (x[2]-0.5) * cos(π/3) - (x[3] - 0.5) * sin(π/3) - x[2])/2,
        (0.5 + (x[2]-0.5) * sin(π/3) + (x[3] - 0.5) * cos(π/3) - x[3])/2
    )

    U = TrialFESpace(V, [g0, g1])

    #FE problem
    op = FEOperator(U, V, a_Ω, t_Γ)
    solver = FESolver(nls)

    x0 = zeros(Float64, num_free_dofs(V))
    uh = FEFunction(U, x0)
    uh, = solve!(uh, solver, op)

    writevtk(trian,"results_hyper",cellfields=["uh"=>uh,"sigma"=>S(∇(uh))])
end
@btime run()

Newton Iterations

Iter     f(x) inf-norm    Step 2-norm
------   --------------   --------------
     0     1.078283e-01              NaN
     1     2.437673e-02     8.468532e+01
     2     4.816547e-03     1.918918e+00
     3     7.730276e-04     5.282747e-02
     4     1.071959e-05     3.422126e-05
     5     9.712950e-09     1.616911e-08

Output

[bhaveshshrimali]

Just one last question before I close this: how to create a mesh with tets instead of hexs ? I am planning to take a look here but that would require some more reading and understanding on my part.

Also for I/O is XDMF a format that you might be considering in the future especially with GridapDistributed ?

[santiagobadia]

Great!

Just simplexify your current model:

model = CartesianDiscreteModel(domain,partition)
tmodel = simplexify(model)

You could also use unstructed meshes, take a look at the corresponding tutorial.

The timing results are amazing. We have always tried to implement things the right way, type-stability, zero allocations in numerically intensive parts, etc. But we have not optimised the code yet. So, as a starting point, it is VERY sound.

[bhaveshshrimali]

Thanks @santiagobadia @fverdugo This was a nice introductory exercise. Will continue to read and explore Gridap in the weeks to come. My long term goal is to also use Gridap for solving larger problems, and if possible on a cluster.

[bhaveshshrimali]

With the newer syntax via 0.15.0 (slightly different timings, but mostly the same)

Iter     f(x) inf-norm    Step 2-norm 
------   --------------   --------------
     0     1.078283e-01              NaN
     1     2.166898e-02     8.536010e+01
     2     6.889270e-03     1.942855e+00
     3     8.970176e-04     6.448036e-02
     4     1.966333e-05     5.180504e-05
     5     1.936679e-08     2.823595e-08
     6     1.600196e-14     1.315987e-14
  20.919 s (2131237 allocations: 3.54 GiB)

Code

using Gridap
using LineSearches: BackTracking, StrongWolfe, HagerZhang
using BenchmarkTools

E = 10.
ν = 0.3
const μ = E/2/(1+ν)
const λ = 2*μ*ν/(1-2*ν)

function run()

    # deformation gradient

    F(∇u) = one(∇u) + ∇u
    J(F) = det(F)
    C(F) = (F') * F

    # body force
    B(x) = VectorValue(
        0.0, -0.5, 0.0
    )

    # surfaceTraction
    T(x) = VectorValue(
        -0.1, 0.0, 0.0
    )

    # constitutive relation
    function S(∇u)
        Js = J(F(∇u))
        μ * (F(∇u) - (inv(F(∇u)))') + λ * Js * (Js - 1) * (inv(F(∇u)))' 
    end

    function res(u, v)
        ∫( (S∘∇(u) ⊙ ∇(v)) - B ⊙ v )*dΩ - ∫( T ⊙ v )*dΓ_t
    end

    function Lv(∇du, ∇u)
        Js = J(F(∇u))
        Finv = inv(F(∇u))
        μ * (∇du + Finv' ⋅ ∇du' ⋅ Finv') + λ * (2. * Js - 1.) * Js * ((Finv' ⊙ ∇du) * Finv') - λ * Js * (Js - 1) * (Finv' ⋅ ∇du' ⋅ Finv') 
    end

    jac(u, du, v) = ∫( Lv∘(∇(du), ∇(u)) ⊙ ∇(v) )*dΩ

    domain = (0.0, 1.0, 0.0, 1.0, 0.0, 1.0)
    partition = (24, 16, 16)
    model = CartesianDiscreteModel(domain, partition)

    labels = get_face_labeling(model)
    add_tag_from_tags!(labels, "rightFace", [2, 4, 6, 8, 14, 16, 18, 20, 26])
    add_tag_from_tags!(labels, "leftFace", [1, 3, 5, 7, 13, 15, 17, 19, 25])

    reffe = ReferenceFE(:Lagrangian, VectorValue{3, Float64}, 1)
    V = TestFESpace(model, reffe, conformity=:H1, dirichlet_tags = ["leftFace", "rightFace"])

    Ω = Triangulation(model)
    degree = 2
    dΩ = LebesgueMeasure(Ω, degree)

    neumannTags = ["boundary"]
    Γ_t = BoundaryTriangulation(model, tags=neumannTags)
    dΓ_t = LebesgueMeasure(Γ_t, degree)

    nls = NLSolver(
    show_trace=true,
    method=:newton,
    linesearch=BackTracking()
    )

    g0(x) = VectorValue(0.0, 0.0, 0.0)
    g1(x) = VectorValue(
        0.0,
        (0.5 + (x[2]-0.5) * cos(π/3) - (x[3] - 0.5) * sin(π/3) - x[2])/2,
        (0.5 + (x[2]-0.5) * sin(π/3) + (x[3] - 0.5) * cos(π/3) - x[3])/2
    )

    U = TrialFESpace(V, [g0, g1])

    #FE problem
    op = FEOperator(res, jac, U, V)
    solver = FESolver(nls)

    x0 = zeros(Float64, num_free_dofs(V))
    uh = FEFunction(U, x0)
    uh, = solve!(uh, solver, op)

    writevtk(Ω, "resultsHyperNew", cellfields=["uh"=>uh, "sigma"=>S∘∇(uh)])
end
run()

Result