Ferrite-FEM · termi-official · May 6, 2024 · May 6, 2024 · Jul 16, 2024
diff --git a/docs/Manifest.toml b/docs/Manifest.toml
diff --git a/docs/Project.toml b/docs/Project.toml
@@ -1,10 +1,12 @@
 [deps]
 BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+Ferrite = "c061ca5d-56c9-439f-9c0e-210fe06d3992"
 FerriteDistributed = "570c3397-5fe4-4792-be0d-48dbf0d14605"
 HYPRE = "b5ffcf37-a2bd-41ab-a3da-4bd9bc8ad771"
 IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
+MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195"
 MPIPreferences = "3da0fdf6-3ccc-4f1b-acd9-58baa6c99267"
 Metis = "2679e427-3c69-5b7f-982b-ece356f1e94b"
 PartitionedArrays = "5a9dfac6-5c52-46f7-8278-5e2210713be9"

diff --git a/docs/make.jl b/docs/make.jl
@@ -25,8 +25,8 @@ makedocs(
     sitename = "FerriteDistributed.jl",
     format = Documenter.HTML(),
     doctest = false,
-    strict = false,
-    draft = true,
+    warnonly = true,
+    draft = liveserver,
     pages = Any[
         "Home" => "index.md",
         "Examples" => [GENERATEDEXAMPLES;],

diff --git a/docs/src/literate/heat_equation_hypre.jl b/docs/src/literate/heat_equation_hypre.jl
@@ -22,24 +22,28 @@ import FerriteDistributed: getglobalgrid, global_comm, local_dof_range #TODO REM
 # Launch MPI and HYPRE
 MPI.Init()
 HYPRE.Init()
+#md nothing # hide
 
 # We start generating a simple grid with 10x10x10 hexahedral elements
 # and distribute it across our processors using `generate_distributed_grid`.
-dgrid = generate_nod_grid(MPI.COMM_WORLD, Hexahedron, (10, 10, 10); partitioning_alg=FerriteDistributed.PartitioningAlgorithm.Metis(:RECURSIVE));
+dgrid = generate_nod_grid(MPI.COMM_WORLD, Hexahedron, (10, 10, 10); partitioning_alg=FerriteDistributed.PartitioningAlgorithm.Metis(:RECURSIVE))
+#md nothing # hide
 
 # ### Trial and test functions
 # Nothing changes here.
 ref = RefHexahedron
 ip = Lagrange{ref, 2}()
 ip_geo = Lagrange{ref, 1}()
 qr = QuadratureRule{ref}(3)
-cellvalues = CellValues(qr, ip, ip_geo);
+cellvalues = CellValues(qr, ip, ip_geo)
+#md nothing # hide
 
 # ### Degrees of freedom
 # Nothing changes here, too. The constructor takes care of creating the correct distributed dof handler.
 dh = DofHandler(dgrid)
 push!(dh, :u, 1, ip)
-close!(dh);
+close!(dh)
+#md nothing # hide
 
 # ### Boundary conditions
 # Nothing has to be changed here either.
@@ -48,65 +52,67 @@ ch = ConstraintHandler(dh);
 dbc = Dirichlet(:u, ∂Ω, (x, t) -> 0)
 dbc_val = 0                                 #src
 dbc = Dirichlet(:u, ∂Ω, (x, t) -> dbc_val)  #src
-add!(ch, dbc);
+add!(ch, dbc)
 close!(ch)
+#md nothing # hide
 
 # ### Assembling the linear system
-# Assembling the system works also mostly analogue.
-function doassemble(cellvalues::CellValues, dh::FerriteDistributed.NODDofHandler{dim}, ch::ConstraintHandler) where {dim}
+# Assembling the system works also mostly analogue. We use the same function to assemble the element as in the serial version.
+function assemble_element!(Ke, fe, cellvalues, cell_coords::AbstractVector{<:Vec{dim}}) where dim
+    fill!(Ke, 0)
+    fill!(fe, 0)
+
     n_basefuncs = getnbasefunctions(cellvalues)
-    Ke = zeros(n_basefuncs, n_basefuncs)
-    fe = zeros(n_basefuncs)
-
-    # --------------------- Distributed assembly --------------------
-    # The synchronization with the global sparse matrix is handled by
-    # an assembler again. You can choose from different backends, which
-    # are described in the docs and will be expaned over time. This call
-    # may trigger a large amount of communication.
+    for q_point in 1:getnquadpoints(cellvalues)
+        dΩ = getdetJdV(cellvalues, q_point)
+
+        for i in 1:n_basefuncs
+            v  = shape_value(cellvalues, q_point, i)
+            ∇v = shape_gradient(cellvalues, q_point, i)
+            ## Manufactured solution of Π cos(xᵢπ)
+            x = spatial_coordinate(cellvalues, q_point, cell_coords)
+            fe[i] += (π/2)^2 * dim * prod(cos, x*π/2) * v * dΩ
+
+            for j in 1:n_basefuncs
+                ∇u = shape_gradient(cellvalues, q_point, j)
+                Ke[i, j] += (∇v ⋅ ∇u) * dΩ
+            end
+        end
+    end
+end
+#md nothing # hide
 
-    # TODO how to put this into an interface?
+# The synchronization with the global sparse matrix is handled by
+# an assembler again. You can choose from different backends, which
+# are described in the docs and will be expaned over time. This call
+# may trigger a large amount of communication.
+function doassemble(cellvalues::CellValues, dh::FerriteDistributed.NODDofHandler, ch::ConstraintHandler)
+    ## TODO how to put this into an interface?
     dgrid = getglobalgrid(dh)
     comm = global_comm(dgrid)
     ldofrange = local_dof_range(dh)
     K = HYPREMatrix(comm, first(ldofrange), last(ldofrange))
     f = HYPREVector(comm, first(ldofrange), last(ldofrange))
-
     assembler = start_assemble(K, f)
 
-    # For the local assembly nothing changes
-    for cell in CellIterator(dh)
-        fill!(Ke, 0)
-        fill!(fe, 0)
+    n_basefuncs = getnbasefunctions(cellvalues)
+    Ke = zeros(n_basefuncs, n_basefuncs)
+    fe = zeros(n_basefuncs)
 
+    ## For the local assembly nothing changes
+    for cell in CellIterator(dh)
         reinit!(cellvalues, cell)
-        coords = getcoordinates(cell)
-
-        for q_point in 1:getnquadpoints(cellvalues)
-            dΩ = getdetJdV(cellvalues, q_point)
-
-            for i in 1:n_basefuncs
-                v  = shape_value(cellvalues, q_point, i)
-                ∇v = shape_gradient(cellvalues, q_point, i)
-                # Manufactured solution of Π cos(xᵢπ)
-                x = spatial_coordinate(cellvalues, q_point, coords)
-                fe[i] += (π/2)^2 * dim * prod(cos, x*π/2) * v * dΩ
-
-                for j in 1:n_basefuncs
-                    ∇u = shape_gradient(cellvalues, q_point, j)
-                    Ke[i, j] += (∇v ⋅ ∇u) * dΩ
-                end
-            end
-        end
-
+        assemble_element!(Ke, fe, cellvalues, getcoordinates(cell))
+        ## Local elimination of boundary conditions, because global 
+        ## elimination is not implemented for the HYPRE extension.
         apply_local!(Ke, fe, celldofs(cell), ch)
-
-        # TODO how to put this into an interface.
+        ## TODO how to put this into an interface?
         assemble!(assembler, dh.ldof_to_gdof[celldofs(cell)], fe, Ke)
     end
 
-    # Finally, for the `HYPREAssembler` we have to call
-    # `end_assemble` to construct the global sparse matrix and the global
-    # right hand side vector.
+    ## Finally, for the `HYPREAssembler` we have to call
+    ## `end_assemble` to construct the global sparse matrix and the global
+    ## right hand side vector.
     end_assemble(assembler)
 
     return K, f
@@ -115,30 +121,34 @@ end
 
 # ### Solution of the system
 # Again, we assemble our problem. Note that we applied the constraints locally.
-K, f = doassemble(cellvalues, dh, ch);
+K, f = doassemble(cellvalues, dh, ch)
+#md nothing # hide
 
 # We use CG with AMG preconditioner to solve the system.
 precond = HYPRE.BoomerAMG()
 solver = HYPRE.PCG(; Precond = precond)
 uₕ = HYPRE.solve(solver, K, f)
+#md nothing # hide
 
 # And convert the solution from HYPRE to Ferrite
 u_local = Vector{Float64}(undef, FerriteDistributed.num_local_dofs(dh))
 FerriteDistributed.extract_local_part!(u_local, uₕ, dh)
+#md nothing # hide
 
 # ### Exporting via PVTK
 # To visualize the result we export the grid and our field `u`
 # to a VTK-file, which can be viewed in e.g. [ParaView](https://www.paraview.org/).
+# For debugging purposes it can be helpful to enrich
+# the visualization with some meta  information about
+# the grid and its partitioning.
 vtk_grid("heat_equation_distributed", dh) do vtk
     vtk_point_data(vtk, dh, u_local)
-    # For debugging purposes it can be helpful to enrich
-    # the visualization with some meta  information about
-    # the grid and its partitioning
     vtk_shared_vertices(vtk, dgrid)
     vtk_shared_faces(vtk, dgrid)
     vtk_shared_edges(vtk, dgrid) #src
     vtk_partitioning(vtk, dgrid)
 end
+#md nothing # hide
 
 ## Test the result against the manufactured solution                    #src
 using Test                                                              #src

diff --git a/docs/src/literate/heat_equation_pa.jl b/docs/src/literate/heat_equation_pa.jl
@@ -20,24 +20,28 @@ using IterativeSolvers
 
 # Launch MPI
 MPI.Init()
+#md nothing # hide
 
 # We start generating a simple grid with 20x20 quadrilateral elements
 # and distribute it across our processors using `generate_distributed_grid`. 
-dgrid = generate_nod_grid(MPI.COMM_WORLD, Hexahedron, (10, 10, 10); partitioning_alg=FerriteDistributed.PartitioningAlgorithm.Metis(:RECURSIVE));
+dgrid = generate_nod_grid(MPI.COMM_WORLD, Hexahedron, (10, 10, 10); partitioning_alg=FerriteDistributed.PartitioningAlgorithm.Metis(:RECURSIVE))
+#md nothing # hide
 
 # ### Trial and test functions
 # Nothing changes here.
 ref = RefHexahedron
 ip = Lagrange{ref, 2}()
 ip_geo = Lagrange{ref, 1}()
 qr = QuadratureRule{ref}(3)
-cellvalues = CellValues(qr, ip, ip_geo);
+cellvalues = CellValues(qr, ip, ip_geo)
+#md nothing # hide
 
 # ### Degrees of freedom
 # Nothing changes here, too. The constructor takes care of creating the correct distributed dof handler.
 dh = DofHandler(dgrid)
 add!(dh, :u, ip)
-close!(dh);
+close!(dh)
+#md nothing # hide
 
 # ### Boundary conditions
 # Nothing has to be changed here either.
@@ -46,94 +50,101 @@ ch = ConstraintHandler(dh);
 dbc = Dirichlet(:u, ∂Ω, (x, t) -> 0)
 dbc_val = 0                                 #src
 dbc = Dirichlet(:u, ∂Ω, (x, t) -> dbc_val)  #src
-add!(ch, dbc);
+add!(ch, dbc)
 close!(ch)
-update!(ch, 0.0);
+#md nothing # hide
 
 # ### Assembling the linear system
 # Assembling the system works also mostly analogue. Note that the dof handler type changed.
-function doassemble(cellvalues::CellValues, dh::FerriteDistributed.NODDofHandler{dim}) where {dim}
+# We use the same function to assemble the element as in the serial version.
+function assemble_element!(Ke, fe, cellvalues, cell_coords::AbstractVector{<:Vec{dim}}) where dim
+    fill!(Ke, 0)
+    fill!(fe, 0)
+
     n_basefuncs = getnbasefunctions(cellvalues)
-    Ke = zeros(n_basefuncs, n_basefuncs)
-    fe = zeros(n_basefuncs)
+    for q_point in 1:getnquadpoints(cellvalues)
+        dΩ = getdetJdV(cellvalues, q_point)
+
+        for i in 1:n_basefuncs
+            v  = shape_value(cellvalues, q_point, i)
+            ∇v = shape_gradient(cellvalues, q_point, i)
+            ## Manufactured solution of Π cos(xᵢπ)
+            x = spatial_coordinate(cellvalues, q_point, cell_coords)
+            fe[i] += (π/2)^2 * dim * prod(cos, x*π/2) * v * dΩ
+
+            for j in 1:n_basefuncs
+                ∇u = shape_gradient(cellvalues, q_point, j)
+                Ke[i, j] += (∇v ⋅ ∇u) * dΩ
+            end
+        end
+    end
+end
+#md nothing # hide
 
-    # --------------------- Distributed assembly --------------------
-    # The synchronization with the global sparse matrix is handled by 
-    # an assembler again. You can choose from different backends, which
-    # are described in the docs and will be expaned over time. This call
-    # may trigger a large amount of communication.
-    # NOTE: At the time of writing the only backend available is a COO 
-    #       assembly via PartitionedArrays.jl .
+# The synchronization with the global sparse matrix is handled by 
+# an assembler again. You can choose from different backends, which
+# are described in the docs and will be expaned over time. This call
+# may trigger a large amount of communication.
+# NOTE: At the time of writing the only backend available is a COO 
+#       assembly via PartitionedArrays.jl.
+function doassemble(cellvalues::CellValues, dh::FerriteDistributed.NODDofHandler{dim}) where {dim}
     assembler = start_assemble(dh, distribute_with_mpi(LinearIndices((MPI.Comm_size(MPI.COMM_WORLD),))))
 
-    # For the local assembly nothing changes
+    ## For the local assembly nothing changes
+    n_basefuncs = getnbasefunctions(cellvalues)
+    Ke = zeros(n_basefuncs, n_basefuncs)
+    fe = zeros(n_basefuncs)
     for cell in CellIterator(dh)
         fill!(Ke, 0)
         fill!(fe, 0)
-
         reinit!(cellvalues, cell)
-        coords = getcoordinates(cell)
-
-        for q_point in 1:getnquadpoints(cellvalues)
-            dΩ = getdetJdV(cellvalues, q_point)
-
-            for i in 1:n_basefuncs
-                v  = shape_value(cellvalues, q_point, i)
-                ∇v = shape_gradient(cellvalues, q_point, i)
-                # Manufactured solution of Π cos(xᵢπ)
-                x = spatial_coordinate(cellvalues, q_point, coords)
-                fe[i] += (π/2)^2 * dim * prod(cos, x*π/2) * v * dΩ
-
-                for j in 1:n_basefuncs
-                    ∇u = shape_gradient(cellvalues, q_point, j)
-                    Ke[i, j] += (∇v ⋅ ∇u) * dΩ
-                end
-            end
-        end
-
-        # Note that this call should be communication-free!
+        assemble_element!(Ke, fe, cellvalues, getcoordinates(cell))
+        ## Note that this call should be communication-free!
         Ferrite.assemble!(assembler, celldofs(cell), fe, Ke)
     end
 
-    # Finally, for the `PartitionedArraysCOOAssembler` we have to call
-    # `end_assemble` to construct the global sparse matrix and the global
-    # right hand side vector.
+    ## Finally, for the `PartitionedArraysCOOAssembler` we have to call
+    ## `end_assemble` to construct the global sparse matrix and the global
+    ## right hand side vector.
     return end_assemble(assembler)
 end
 #md nothing # hide
 
 # ### Solution of the system
 # Again, we assemble our problem and apply the constraints as needed.
-K, f = doassemble(cellvalues, dh);
+K, f = doassemble(cellvalues, dh)
 apply!(K, f, ch)
+#md nothing # hide
 
 # To compute the solution we utilize conjugate gradients because at the time of writing 
 # this is the only available scalable working solver.
 # Additional note: At the moment of writing this we have no good preconditioners for PSparseMatrix in Julia, 
 # partly due to unimplemented multiplication operators for the matrix data type.
 u = cg(K, f)
+#md nothing # hide
 
 # ### Exporting via PVTK
 # To visualize the result we export the grid and our field `u`
 # to a VTK-file, which can be viewed in e.g. [ParaView](https://www.paraview.org/).
+# For debugging purposes it can be helpful to enrich 
+# the visualization with some meta  information about 
+# the grid and its partitioning.
 vtk_grid("heat_equation_distributed", dh) do vtk
     vtk_point_data(vtk, dh, u)
-    # For debugging purposes it can be helpful to enrich 
-    # the visualization with some meta  information about 
-    # the grid and its partitioning
     vtk_shared_vertices(vtk, dgrid)
     vtk_shared_faces(vtk, dgrid)
     vtk_shared_edges(vtk, dgrid) #src
     vtk_partitioning(vtk, dgrid)
 end
+#md nothing # hide
 
 ## Test the result against the manufactured solution                    #src
 using Test                                                              #src
 for cell in CellIterator(dh)                                            #src
     reinit!(cellvalues, cell)                                           #src
     n_basefuncs = getnbasefunctions(cellvalues)                         #src
     coords = getcoordinates(cell)                                       #src
-    map(local_values(u)) do u_local                         #src
+    map(local_values(u)) do u_local                                     #src
         uₑ = u_local[celldofs(cell)]                                    #src
         for q_point in 1:getnquadpoints(cellvalues)                     #src
             x = spatial_coordinate(cellvalues, q_point, coords)         #src

diff --git a/ext/FerriteDistributedHYPREAssembly.jl b/ext/FerriteDistributedHYPREAssembly.jl
@@ -5,7 +5,7 @@ module FerriteDistributedHYPREAssembly
 
 using FerriteDistributed
 # TODO remove me. These are merely hotfixes to split the extensions trasiently via an internal API.
-import FerriteDistributed: getglobalgrid, num_local_true_dofs, num_local_dofs, global_comm, interface_comm, global_rank, compute_owner, remote_entities, num_fields
+import FerriteDistributed: getglobalgrid, num_local_true_dofs, num_local_dofs, global_comm, interface_comm, global_rank, compute_owner, remote_entities, num_fields, local_entities
 using MPI
 using HYPRE
 using Base: @propagate_inbounds