Example312_PeriodicElasticity3D.jl

#=

# 312 : Periodic Boundary 3D
([source code](@__SOURCE_URL__))

This is a simple demonstration and validation of the generic periodic boundary operator.

We construct an unstructured periodic 3D grid and solve a simple linear elastic problem
with periodic coupling along the x-axis.

![](example312.png)
=#

module Example312_PeriodicElasticity3D

using ExtendableFEM
using ExtendableGrids
using SimplexGridFactory
using GridVisualize
using TetGen
using UnicodePlots
using StaticArrays
using LinearAlgebra
using Test #hide

const reg_left = 1
const reg_right = 2
const reg_dirichlet = 3
const reg_default = 4


# define the Hooke tensor for AlN (from DOI 10.1063/1.1368156)
function material_tensor()
    c11 = 396.0
    c12 = 137.0
    c13 = 108.0
    c33 = 373.0
    c44 = 116.0

    return @SArray [
        c11 c12 c13 0   0   0
        c12 c11 c13 0   0   0
        c13 c13 c33 0   0   0
        0   0   0   c44 0   0
        0   0   0   0   c44 0
        0   0   0   0   0   c44
    ]
end

## generate the kernels for the linear problem
## 𝐂: Hooke tensor, 𝑓: body force
function make_kernels(𝐂, 𝑓)

    ## linear stress-strain mapping
    bilinear_kernel!(σ, εv, qpinfo) = mul!(σ, 𝐂, εv)

    ## body force
    linear_kernel!(result, qpinfo) = (result .= 𝑓)

    return bilinear_kernel!, linear_kernel!
end


"""
    create 3D grid with Dirichlet boundary region at the bottom center
"""
function create_grid(; h, height, width, depth)
    builder = SimplexGridBuilder(; Generator = TetGen)

    ## bottom points
    b01 = point!(builder, 0, 0, 0)
    b02 = point!(builder, 0.45 * width, 0, 0)
    b03 = point!(builder, 0.55 * width, 0, 0)
    b04 = point!(builder, width, 0, 0)

    b11 = point!(builder, 0, depth, 0)
    b12 = point!(builder, 0.45 * width, depth, 0)
    b13 = point!(builder, 0.55 * width, depth, 0)
    b14 = point!(builder, width, depth, 0)

    ## top points
    t01 = point!(builder, 0, 0, height)
    t02 = point!(builder, width, 0, height)

    t11 = point!(builder, 0, depth, height)
    t12 = point!(builder, width, depth, height)

    ## center points
    c01 = point!(builder, 0.5 * width, 0, 0)
    c02 = point!(builder, 0.5 * width, 0, height)
    c11 = point!(builder, 0.5 * width, depth, 0)
    c12 = point!(builder, 0.5 * width, depth, height)

    ## default faces
    facetregion!(builder, reg_default)
    facet!(builder, b01, b02, b12, b11)
    facet!(builder, b03, b04, b14, b13)
    facet!(builder, [t01, c02, c12, t11])
    facet!(builder, [c02, t02, t12, c12])
    facet!(builder, [b01, b02, c01, c02, t01])
    facet!(builder, [c01, b03, b04, t02, c02])
    facet!(builder, [b11, b12, c11, c12, t11])
    facet!(builder, [c11, b13, b14, t12, c12])
    facet!(builder, c01, c02, c12, c11)

    ## left face
    facetregion!(builder, reg_left)
    facet!(builder, b01, t01, t11, b11)

    ## right face
    facetregion!(builder, reg_right)
    facet!(builder, b04, t02, t12, b14)

    ## Dirichlet face
    facetregion!(builder, reg_dirichlet)
    facet!(builder, [b02, c01, c11, b12])
    facet!(builder, [c01, b03, b13, c11])

    cellregion!(builder, 1)
    maxvolume!(builder, h)
    regionpoint!(builder, width / 3, depth / 2, height / 2)
    cellregion!(builder, 2)
    ## finer grid on the right half to make the periodic coupling non-trivial
    maxvolume!(builder, 0.3 * h)
    regionpoint!(builder, 2 * width / 3, depth / 2, height / 2)

    return simplexgrid(builder)
end

function main(;
        order = 1,
        periodic_coupling = :high_level_restriction, # :restriction, :operator, :high_level_restriction
        Plotter = nothing,
        force = 1.0,
        h = 1.0e-4,
        width = 6.0,
        height = 0.2,
        depth = 1,
        threads = 1,
        kwargs...
    )
    ## print options for better logs
    @info "selected options" periodic_coupling threads

    xgrid = create_grid(; h, width, height, depth)

    ## create finite element space and solution vector
    if order == 1
        FES = FESpace{H1P1{3}}(xgrid)
    elseif order == 2
        FES = FESpace{H1P2{3, 3}}(xgrid)
    end

    ## problem description
    PD = ProblemDescription()
    u = Unknown("u"; name = "displacement")
    assign_unknown!(PD, u)

    𝐂 = material_tensor()
    𝑓 = force * [0, 0, 1]

    bilinear_kernel!, linear_kernel! = make_kernels(𝐂, 𝑓)
    assign_operator!(PD, BilinearOperator(bilinear_kernel!, [εV(u, 1.0)]; kwargs...))
    assign_operator!(PD, LinearOperator(linear_kernel!, [id(u)]; kwargs...))

    assign_operator!(PD, HomogeneousBoundaryData(u; regions = [reg_dirichlet]))

    if periodic_coupling == :high_level_restriction

        # new high-level call
        assign_restriction!(PD, CoupledDofsRestriction(u, reg_left, reg_right; parallel = threads > 1, threads))

    elseif periodic_coupling != :none

        function give_opposite!(y, x)
            y .= x
            y[1] = width - x[1]
            return nothing
        end
        @showtime coupling_matrix = get_periodic_coupling_matrix(FES, reg_left, reg_right, give_opposite!; parallel = threads > 1, threads)
        if periodic_coupling == :restriction
            assign_restriction!(PD, CoupledDofsRestriction(coupling_matrix))
        else # :operator
            assign_operator!(PD, CombineDofs(u, u, coupling_matrix; kwargs...))
        end
    end

    ## solve
    sol, SC = solve(PD, FES; return_config = true, kwargs...)
    residual(SC) < 1.0e-10 || error("Residual is not zero!")

    @info "Lagrange residuals" SC.statistics[:restriction_residuals]

    displace_mesh!(xgrid, sol[u])
    plt = plot([grid(u)], sol; Plotter, do_vector_plots = false, width = 1200, height = 800, title = "displaced mesh", scene3d = :LScene)

    return sol, plt

end

generateplots = ExtendableFEM.default_generateplots(Example312_PeriodicElasticity3D, "example312.png") #hide
function runtests()                                                                                  #hide
    sol1, _ = main(periodic_coupling = :operator, threads = 1)                                       #hide
    @test abs(maximum(view(sol1[1])) - 1.8004602502175202) < 2.0e-3                                  #hide

    sol2, _ = main(periodic_coupling = :operator, threads = 4)                                       #hide
    @test sol1.entries ≈ sol2.entries                                                                #hide

    sol3, _ = main(periodic_coupling = :restriction, threads = 4)                                    #hide
    @test sol1.entries ≈ sol3.entries                                                                #hide

    sol4, _ = main(periodic_coupling = :high_level_restriction, threads = 4)                         #hide
    @test sol1.entries ≈ sol4.entries                                                                #hide

    return nothing                                                                                   #hide
end                                                                                                  #hide

end # module