a lot more stuff for juliac. Got a full app working!

examples/app/Project.toml

@@ -3,5 +3,9 @@ name = "MyApp"
 [deps]
 Exodus = "f57ae99e-f805-4780-bdca-96e224be1e5a"
 FiniteElementContainers = "d08262e4-672f-4e7f-a976-f2cea5767631"
+Krylov = "ba0b0d4f-ebba-5204-a429-3ac8c609bfb7"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 ReferenceFiniteElements = "6dc62d09-f8eb-43fd-9672-074e490a997f"
+Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
+TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"

examples/app/input-file.toml

@@ -1,26 +1,26 @@
 [device]
 backend = "cpu"
 
-[functions.zero]
+[functions.zero_ic]
 type = "constant"
 expression = "0.0"
 variables = ["x", "y"]
 
-[functions.one]
+[functions.zero_bc]
 type = "constant"
-expression = "1.0"
+expression = "0.0"
 variables = ["x", "y", "t"]
 
 [mesh]
 file_path = "poisson.g"
 file_type = "exodus"
 
 [[initial_conditions]]
-function = "zero"
-block = "block_1"
-variable = "u"
+function = "zero_ic"
+blocks = ["block_1"]
+variables = ["u"]
 
 [[boundary_conditions.dirichlet]]
-function = "one"
-side_set = "sset_1"
-variable = "u"
+function = "zero_bc"
+side_sets = ["sset_1", "sset_2", "sset_3", "sset_4"]
+variables = ["u"]

examples/app/build.jl → examples/app/make.jl

@@ -10,7 +10,7 @@ img = ImageRecipe(
     file           = "$src_path/src/MyApp.jl",
     trim_mode      = "safe",
     add_ccallables = false,
-    verbose        = true,
+    verbose        = false,
 )
 
 link = LinkRecipe(

examples/app/src/MyApp.jl

@@ -2,46 +2,92 @@ import FiniteElementContainers as FEC
 import FiniteElementContainers.AppTools as AT
 using Exodus
 using FiniteElementContainers
+using Krylov
+using LinearAlgebra
+using TimerOutputs
 
 # include files
 include("Physics.jl")
 
-f(X, _) = 2. * π^2 * sin(π * X[1]) * sin(π * X[2])
-bc_zero_func(_, _) = 0.0
+f(X, _) = 2. * π^2 * sin(2π * X[1]) * sin(2π * X[2])
 
 function app_main(ARGS::Vector{String})
     app = AT.App("MyApp")
     AT.add_cli_arg!(app, "--backend"; default = "cpu")
     sim = AT.setup(app, ARGS)
 
+    #####################################
+    # need to define some types
+    #####################################
+    FT = AT.ExpressionFunction{Float64}
+    SPT = FEC.CSCMatrix()
+
     #####################################
     # setup function space
     #####################################
+    # block_ids = FEC._setup_juliac_safe_block_to_ref_fe_id(sim.mesh)
+    # for n in 1:16
+    #     println(Core.stdout, "Block $n element type $(block_ids[n])")
+    # end
     V = FunctionSpace{true}(sim.mesh, H1Field, Lagrange)
     u = ScalarFunction(V, "u")
     dof = DofManager{false}(u)
-    sp_type = FiniteElementContainers.CSCMatrix
-    asm = SparseMatrixAssembler{sp_type, false, false}(dof)
+    asm = SparseMatrixAssembler{SPT, false, false}(dof)
 
-    physics = Poisson(f)
-    props = create_properties(physics)
+    # setup physics and properties
+    physics = Poisson{typeof(f)}[]
+    props = Vector{Float64}[]
+    for _ in 1:length(sim.mesh.element_conns)
+        temp_physics = Poisson(f)
+        push!(physics, temp_physics)
+        push!(props, create_properties(temp_physics))
+    end
 
-    # U = create_unknowns(asm)
-    U = create_field(asm)
+    times = TimeStepper(0.0, 0.0, 1)
+    p = FEC.TypeStableParameters{FT}(sim.mesh, asm, physics, props, sim.ics, sim.dbcs, times)
+    # U = create_field(asm)
+    # Uu = create_unknowns(asm)
 
-    # setting up ics
-    FT = AT.ExpressionFunction{Float64}
-    ics = InitialConditions{FT}(sim.mesh, dof, sim.ics)
-    update_ic_values!(ics, sim.mesh.nodal_coords)
-    update_field_ics!(U, ics)
+    # # setting up ics
+    # # ics = InitialConditions{FT}(sim.mesh, dof, sim.ics)
+    # X = sim.mesh.nodal_coords
+    # # println(Core.stdout, "X size 1 = $(size(X, 1))")
+    # # println(Core.stdout, "X size 2 = $(size(X, 2))")
+    # # for block_id in 1:1
+    # #     ref_fe = FEC.block_reference_element(V, block_id)
+    # #     println(Core.stdout, ref_fe)
+    # # end
+
+    FEC.initialize!(p)
+    # FEC._update_for_assembly!(p, dof, Uu)
 
-    # setting up dbcs
-    dbcs = DirichletBCs{FT}(sim.mesh, dof, sim.dbcs)
-    update_bc_values!(dbcs, sim.mesh.nodal_coords, 1.0)
-    update_field_dirichlet_bcs!(U, dbcs)
+    # # lsolver = DirectLinearSolver(asm)
+    # # lsolver = x -> IterativeLinearSolver(x, :cg)
+    preconditioner = I
+    timer = TimerOutput()
+    workspace = CgWorkspace(stiffness(asm), residual(asm))
+    ΔUu = create_unknowns(asm)
+    lsolver = IterativeLinearSolver(asm, preconditioner, workspace, timer, ΔUu)
+    nlsolver = NewtonSolver(lsolver)
 
-    # # p = create_parameters(mesh, asm, physics, props)
+    # println(Core.stdout, "Number of blocks = $(length(sim.mesh.element_types))")
+    # # N = length(sim.mesh.element_types)
+    # # # N = Val(N)
+    # # N = Val{N}()
+    # # ref_fe = _get_ref_fe(V, N)
+    # # indices = ntuple(x -> x, N)
+    # # indices = ntuple(x -> x, Val(N))
+    # # s = MyStruct{N}(N)::MyStruct{N}
+    integrator = QuasiStaticIntegrator(nlsolver)
     println(sim.log_file.io, "Setup complete")
+    println(sim.log_file.io, "Solving...")
+    evolve!(integrator, p)
+    println(sim.log_file.io, "Solve complete")
+    println(Core.stdout, maximum(p.field.data))
+    println(Core.stdout, minimum(p.field.data))
+
+    # pp = PostProcessor(sim.mesh, "juliac_output.e", u)
+
 end
 
 function @main(ARGS::Vector{String})

src/AppTools.jl

@@ -14,6 +14,9 @@ import ..FileMesh
 import ..FunctionSpace
 import ..H1Field
 import ..InitialCondition
+import ..NeumannBC
+import ..RobinBC
+import ..Source
 import ..UnstructuredMesh
 import ..element_blocks
 import ..element_ids
@@ -284,6 +287,9 @@ print_dict(l::LogFile, d::Dict{String, String}) = print_dict(l, d)
 #######################################################
 struct BCSettings{T <: Number}
     dirichlet::Vector{DirichletBC{ExpressionFunction{T}}}
+    neumann::Vector{NeumannBC{ExpressionFunction{T}}}
+    robin::Vector{RobinBC{ExpressionFunction{T}}}
+    source::Vector{Source{ExpressionFunction{T}}}
 end
 
 struct FunctionSettings{T <: Number}
@@ -318,35 +324,102 @@ end
 
 function _parse_boundary_condition_settings(log_file, data, functions::FunctionSettings{T}) where T <: Number
     print_banner(log_file, "Boundary conditions")
-    # if length(data["boundary_conditions"])
     if haskey(data, "boundary_conditions")
         bc_settings = data["boundary_conditions"]::Dict{String, Any}
     else
         bc_settings = Dict{String, Any}()
     end
-    dbc_settings = bc_settings["dirichlet"]::Vector{Any}
     dbcs = DirichletBC{ExpressionFunction{T}}[]
-    for bc in dbc_settings
-        # has_key_check(log_file, bc, "function", "key \"function\" not found for DirichletBC")
-        # has_key_check(log_file, bc, "variable", "key \"variable\" not found for DirichletBC")
-        temp = bc::Dict{String, Any}
-        func = temp["function"]::String
-        func = functions.expr_funcs[func]
-        var = temp["variable"]::String
-        if haskey(temp, "block")
-            block = temp["block"]::String
-            push!(dbcs, DirichletBC(var, func; block_name = block))
-        elseif haskey(temp, "node_set")
-            node_set = temp["node_set"]::String
-            push!(dbcs, DirichletBC(var, func; nodeset_name = node_set))
-        elseif haskey(temp, "side_set")
-            side_set = temp["side_set"]::String
-            push!(dbcs, DirichletBC(var, func; sideset_name = side_set))
-        else
-            error_message(log_file, "Requires block, node_set, or side_set")
+    nbcs = NeumannBC{ExpressionFunction{T}}[]
+    rbcs = RobinBC{ExpressionFunction{T}}[]
+    srcs = Source{ExpressionFunction{T}}[]
+    if haskey(bc_settings, "dirichlet")
+        dbc_settings = bc_settings["dirichlet"]::Vector{Any}
+        for bc in dbc_settings
+            # has_key_check(log_file, bc, "function", "key \"function\" not found for DirichletBC")
+            # has_key_check(log_file, bc, "variable", "key \"variable\" not found for DirichletBC")
+            temp = bc::Dict{String, Any}
+            func = temp["function"]::String
+            func = functions.expr_funcs[func]
+            vars = temp["variables"]::Vector{String}
+            if haskey(temp, "blocks")
+                blocks = temp["blocks"]::Vector{String}
+                for block in blocks
+                    for var in vars
+                        push!(dbcs, DirichletBC(var, func; block_name = block))
+                    end
+                end
+            elseif haskey(temp, "node_sets")
+                node_sets = temp["node_sets"]::Vector{String}
+                for node_set in node_sets
+                    for var in vars
+                        push!(dbcs, DirichletBC(var, func; nodeset_name = node_set))
+                    end
+                end
+            elseif haskey(temp, "side_sets")
+                side_sets = temp["side_sets"]::Vector{String}
+                for side_set in side_sets
+                    for var in vars
+                        push!(dbcs, DirichletBC(var, func; sideset_name = side_set))
+                    end
+                end
+            else
+                error_message(log_file, "Requires block, node_set, or side_set")
+            end
+        end
+    end
+
+    if haskey(bc_settings, "neumann")
+        nbc_settings = bc_settings["neumann"]::Vector{Any}
+        for bc in nbc_settings
+            temp = bc::Dict{String, Any}
+            func = temp["function"]::String
+            func = functions.expr_funcs[func]
+            sidesets = temp["side_sets"]::Vector{String}
+            vars = temp["variables"]::Vector{String}
+            for side_set in sidesets
+                for var in vars
+                    push!(nbcs, NeumannBC(var, func, side_set))
+                end
+            end
         end
     end
-    return BCSettings(dbcs)
+
+    if haskey(bc_settings, "robin")
+        rbc_settings = bc_settings["robin"]::Vector{Any}
+        for bc in rbc_settings
+            temp = bc::Dict{String, Any}
+            func = temp["function"]::String
+            func = functions.expr_funcs[func]
+            sidesets = temp["side_sets"]::Vector{String}
+            vars = temp["variables"]::Vector{String}
+            for side_set in sidesets
+                for var in vars
+                    push!(nbcs, RobinBC(var, func, side_set))
+                end
+            end
+        end
+    end
+
+    if haskey(bc_settings, "source")
+        src_settings = bc_settings["source"]::Vector{Any}
+        for src in src_settings
+            temp = src::Dict{String, Any}
+            blocks = temp["blocks"]::Vector{String}
+            func = temp["function"]::String
+            func = functions.expr_funcs[func]
+            sidesets = temp["side_sets"]::Vector{String}
+            vars = temp["variables"]::Vector{String}
+            for block in blocks
+                for var in vars
+                    push!(nbcs, Source(var, func, block))
+                end
+            end
+        end
+    end
+
+    # nbc_settings = 
+    return BCSettings(dbcs, nbcs, rbcs, srcs)
 end
 
 function _parse_function_space_settings(log_file, data)
@@ -395,11 +468,15 @@ function _parse_initial_condition_settings(log_file, data, functions)
     ics = InitialCondition{ExpressionFunction{Float64}}[]
     for ic in ic_settings
         temp = ic::Dict{String, Any}
-        block = temp["block"]::String
+        blocks = temp["blocks"]::Vector{String}
         func = temp["function"]::String
-        var = temp["variable"]::String
         func = functions.expr_funcs[func]
-        push!(ics, InitialCondition(var, func, block))
+        vars = temp["variables"]::Vector{String}
+        for block in blocks
+            for var in vars
+                push!(ics, InitialCondition(var, func, block))
+            end
+        end
     end
     return ics
 end
@@ -585,6 +662,9 @@ struct Simulation{T <: Number, IO, Mesh}
     ics::Vector{InitialCondition{ExpressionFunction{T}}}
     log_file::LogFile{IO}
     mesh::Mesh
+    nbcs::Vector{NeumannBC{ExpressionFunction{T}}}
+    rbcs::Vector{RobinBC{ExpressionFunction{T}}}
+    srcs::Vector{Source{ExpressionFunction{T}}}
 
     function Simulation(settings::InputSettings, log_file::LogFile{IO}) where IO
         print_banner(log_file, "Mesh")
@@ -603,7 +683,10 @@ struct Simulation{T <: Number, IO, Mesh}
         else
             T = Float64
         end
-        new{T, IO, typeof(mesh)}(settings.bcs.dirichlet, settings.ics, log_file, mesh)
+        new{T, IO, typeof(mesh)}(
+            settings.bcs.dirichlet, settings.ics, log_file, mesh,
+            settings.bcs.neumann, settings.bcs.robin, settings.bcs.source
+        )
     end
 end
 
@@ -699,7 +782,6 @@ function generate_app(
 
         function app_main(ARGS::Vector{String})
             app = AT.App(\"$(name)\")
-            AT.add_cli_arg(app, \"--backend\"; default = \"cpu\")
             sim = AT.setup(app, ARGS)
             println(sim.log_file.io, "Setup complete")
         end
@@ -725,14 +807,14 @@ function generate_app(
         src_path = joinpath(@__DIR__)
         @show build_path
         @show src_path
-        rm(build_path; force=true, recursive=true)
+        rm(build_path; force = true, recursive = true)
 
         img = ImageRecipe(
             output_type    = "--output-exe",
             file           = "\$src_path/src/$name.jl",
             trim_mode      = "$trim_str",
             add_ccallables = false,
-            verbose        = true,
+            verbose        = false,
         )
 
         link = LinkRecipe(
@@ -742,9 +824,9 @@ function generate_app(
 
         bun = BundleRecipe(
             link_recipe = link,
-            output_dir  = build_path # or `nothing` to skip bundling
+            #output_dir  = build_path # or `nothing` to skip bundling
+            output_dir  = nothing
         )
-
         compile_products(img)
         link_products(link)
         bundle_products(bun)