Ritz-Method

[jon21as]

Hello,

is there a way to implement a functional and solve it with the ritz method? I have an energy functional that represents elasto-plastic material behaviour at small strains and I would like to solve and implement the problem directly with the functional and not with the weak form. The specific functional is on the picture. The solution is u and pi, where the functional is minimized. Can you link a code example?

Thanks

[JordiManyer]

Hi @jon21as ,

I don't have a code example in mind, but you can just use any julia optimiser + gridap's automatic differentiation to do this.

Code your energy as a function of your variable, using Gridap's API. You can then take the gradient of it like we do in most non-linear tutorials.

If you work with regular optimiser libraries, you may need to wrap these such that your energy and gradient functions take in vectors instead of FEFunctions. You can use uh = FEFunction(V,x) and x = get_free_dof_values(uh) to map between vector and FEFunction representations of your unknowns.

There are more complex things being done in the context of Topology Optimisation in one of the tutorials.

Good luck!

[jon21as]

Hello,

I tried a lot of things and I wanted to share my results. It was not possible to use a regular optimiser library like Optim with the method NelderMead because no solution was found. That´s why I implemented a solution method by myself, which is not optimal, because I used the Newton method but my functional has non-differentiable points. But my code worked and my prof liked it. Thank you!

using Gridap
using Gridap.TensorValues

#tensor norm für elasticity region: Gestaltänderungshypothese
function GEH(x)
    v10=VectorValue([1.0,0.0])
    v01=VectorValue([0.0,1.0])
    return sqrt((((x)⋅v10)⋅v10)*(((x)⋅v10)⋅v10)+(((x)⋅v01)⋅v01)*(((x)⋅v01)⋅v01)-(((x)⋅v10)⋅v10)*(((x)⋅v01)⋅v01)+3*(((x)⋅v01)⋅v10)*(((x)⋅v01)⋅v10)+0.0001)
end

#helperfunction
function get_index(i,j,k,l)
    index=SymFourthOrderTensorValue((1,2,3,4,5,6,7,8,9))
    return index[i,j,k,l]
end

#Kronecker delta
function δ(i,j)
    if i==j
        return 1.0
    else
        return 0.0
    end
end

#Elastizitätstensor
function material_law(λ,μ)
    data=[0.0;0.0;0.0;0.0;0.0;0.0;0.0;0.0;0.0]
    d=2
    for i=1:d
        for k=1:d
            for j=i:d 
                for l=k:d
                    data[get_index(i,j,k,l)]=(λ)*δ(i,j)*δ(k,l)+(μ)*(δ(i,k)*δ(j,l)+δ(i,l)*δ(j,k))
                end
            end
        end
    end
    return SymFourthOrderTensorValue((data[1],data[2],data[3],data[4],data[5],data[6],data[7],data[8],data[9]))
end

#solver function
function find_minimizer(J::Function, x::Vector{Float64}, V::FESpace; tol=1e-7,maxn=20)
    error=1
    println("")
    println("start calculation with tol=$tol and maxn=$maxn")
    println("")
    println("     Step | |J(x(i-1))-J(x(i))| | J(x(i))")
    println("-------------------------------------------------")
    println("=====  0  |       ----        | $(sum(J(FEFunction(V,x))))")

    
    for i in 1:maxn+1
        if error<tol
            println("")
            println("calculation finished after $(i-1) steps with J(x(i-1))-J(x(i))=$error")
            println("")
            println("")
            break
        end
        if i==maxn
            @error "calculation did not converge after $maxn steps, calculation stopped"
        end
        x0=x
        gJ=∇(J,FEFunction(V,x))
        b=assemble_vector(gJ,V)
        print("=")
        H=hessian(J,FEFunction(V,x))
        A=assemble_matrix(H,V,V)
        print("=")
        A=Matrix(A)
        print("=")
        A=inv(A)
        print("=")
        x=x0.-A*b
        print("=")
        error=abs(sum(J(FEFunction(V,x0))-J(FEFunction(V,x))))
        println("  $i  | $error | $(sum(J(FEFunction(V,x))))")
    end
    return FEFunction(V,x)
end

#main computation with setting up the model
function main()
    #create model: 2D rectangle with length laenge and width breite, divided into d1*d2 rectangles
    laenge=50.0
    breite=10.0
    d1=50
    d2=10
    order=2
    model=CartesianDiscreteModel((0,laenge,0,breite),(d1,d2))

    #labels anpassen
    labels=get_face_labeling(model)
    add_tag_from_tags!(labels,"left",[1;3;7])
    add_tag_from_tags!(labels,"right",[2;4;8])
    add_tag_from_tags!(labels,"bottom",[1;2;5])
    add_tag_from_tags!(labels,"top",[3;4;6])
    add_tag_from_tags!(labels,"point",[1])

    #define reference finite element and vector function space for displacement
    reffedisp=ReferenceFE(lagrangian,VectorValue{2,Float64},order)
    V=TestFESpace(model,reffedisp;conformity=:H1,labels=labels,dirichlet_tags=["left","right","point"],
    dirichlet_masks=[(true,false), (true,false), (true,true)])

    #define function to implement inhomogenous Dirichlet boundary conditions
    g(x,t)=VectorValue([0.2*(x[2]-5.0)*t,0.0])

    # define reference finite element and tensor function space for plasticity
    reffeplast=ReferenceFE(lagrangian,SymTracelessTensorValue{2,Float64,3},order)
    Q=TestFESpace(model,reffeplast;conformity=:H1)

    #create trial space for plasticity
    P=TrialFESpace(Q)

    #define start conditions for plasticity
    p_old=FEFunction(Q,zeros(num_free_dofs(P)))
    p0=FEFunction(Q,0.001.*ones(num_free_dofs(P)))

    #define material parameters
    κ=0.001
    E=200000
    Et=1450
    ν=0.3#nu
    λ(E)=(E*ν)/((1+ν)*(1-2*ν))
    μ(E)=E/2/(1+ν)
    σy=250

    C=material_law(λ(E),μ(E))
    H=material_law(λ(Et),μ(Et))


    #define triangulation and measure
    Ω=Triangulation(model)
    dΩ=Measure(Ω,4)

    #define boundary conditions for displacement
    g1=VectorValue([0.0,0.0])
    g2(x)=g(x,0.0)
    g3=VectorValue([0.0,0.0])

    #assigne boundary conditions and create trial space for displacement
    U=TrialFESpace(V,[g1;g2;g3])

    #define start conditions for displacement
    u0=FEFunction(U,0.01.*ones(num_free_dofs(U)))
    x0=[get_free_dof_values(u0);get_free_dof_values(p0)]

    k=1
    for t=0.1:0.1:1.0

        #define boundary conditions for displacement
        g1=VectorValue([0.0,0.0])
        g2new(x)=g(x,t)
        g3=VectorValue([0.0,0.0])

        #assigne boundary conditions and create trial space for displacement
        U=TrialFESpace(V,[g1;g2new;g3])

        #create multi physics space
        UxP=MultiFieldFESpace([U,P])

        println("=====")
        println("t=$t")
        println("=====")

        #functional for displacement and plasticity
        J((u,p))=∫(1/2*(C⊙(ε(u)-p))⊙(ε(u)-p)+1/2*(H⊙p)⊙p+κ/2*∇(p)⊙∇(p)+σy*(GEH∘(p-p_old)))*dΩ
        
        #solve
        res=find_minimizer(J,x0,UxP,tol=1e-2,maxn=20)
        x0=get_free_dof_values(res)
        p0=res.single_fe_functions[2]
        u=res.single_fe_functions[1]
        p_old=p0

        writevtk(Ω,"results$k.vtu",cellfields=["u"=>u,"π"=>p0,"e"=>ε(u),"s"=>(C)⊙(ε(u)-p0),"el"=>ε(u)-p0,"p_GEH"=>GEH∘(p0),"s_GEH"=>GEH∘((C)⊙(ε(u)-p0)),"el_GEH"=>GEH∘(ε(u)-p0)])#,"d"=>ζ_old])

        k=k+1
        println("")
        println("")
    end 
end

main()

[JordiManyer]

I'm glad it worked out!