gridap_EM/fea.jl at main · bmdaj/gridap_EM · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
using Gridap
using Gridap.Geometry
using LinearAlgebra

function fea_init(model, order, degree, dirichlet_tags="None", neumann_tags="None", source_tags="None")

    # Define the FEM space

    reffe = ReferenceFE(lagrangian,Float64,order)

    if dirichlet_tags == "None"
        V = TestFESpace(model,reffe, vector_type=Vector{ComplexF64})
    else
        V = TestFESpace(model,reffe, dirichlet_tags=dirichlet_tags, vector_type=Vector{ComplexF64})
    end

    U = V # mathematically equivalent to TrialFESpace(V,0)
    Ω = Triangulation(model)
    dΩ = Measure(Ω,degree)

    if neumann_tags == "None"

        Γ_n = false
        dΓ_n = false

    else

        Γ_n = BoundaryTriangulation(model;tags=neumann_tags)
        dΓ_n = Measure(Γ_n,degree)

    end


    if source_tags == "None"

        Γ_s = false
        dΓ_s = false

    else

        Γ_s = BoundaryTriangulation(model;tags=source_tags)
        dΓ_s = Measure(Γ_s,degree)


    end

    return U, V, Ω, dΩ, Γ_n, dΓ_n, Γ_s, dΓ_s

end

function fea_init_topopt(model, order, degree, dirichlet_tags="None", neumann_tags="None", design_tags="None", source_tags="None")

    U, V, Ω, dΩ, Γ_n, dΓ_n, Γ_s, dΓ_s = fea_init(model, order, degree, dirichlet_tags, neumann_tags, source_tags)

    if design_tags == "None"

        Ω_d = false
        dΩ_d = false

    else

        Ω_d = Triangulation(model, tags="Design")
        dΩ_d = Measure(Ω_d, degree)

    end

    return U, V, Ω, dΩ, Ω_d, dΩ_d, Γ_n, dΓ_n, Γ_s, dΓ_s

end

function set_tags(model, εₛ, tag_list, ε₀, pol="TE")

    labels = get_face_labeling(model)
    dimension = num_cell_dims(model)
    tags = get_face_tag(labels,dimension)
    τ = CellField(tags,Ω)

    function ξ(tag)

        for i in range(1, length(tag_list))

            tag_name = get_tag_from_name(labels,tag_list[i])

            if tag == tag_name
                if pol == "TE"
                    return εₛ[i]
                elseif pol == "TM"
                    return 1/εₛ[i]
                end
            end
        end
        if pol == "TE"
            return ε₀
        elseif pol == "TM"
            return 1/ε₀
        end
    end

    return ξ, τ

end