Merge pull request #161 from GalerkinToolkit/assembly_matrix

matrix assembly update

GalerkinToolkitExamples/Project.toml

@@ -4,6 +4,8 @@ authors = ["Francesc Verdugo <[email protected]>"]
 version = "0.1.0"
 
 [deps]
+AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 GalerkinToolkit = "5e3ba9c4-dd81-444d-b69a-0e7bd7bf60a4"
 IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153"
@@ -20,14 +22,12 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192"
-BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+
+[compat]
+PartitionedSolvers = "0.3"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
 test = ["Test"]
-
-[compat]
-PartitionedSolvers = "0.3"
-

GalerkinToolkitExamples/src/poisson.jl

@@ -13,6 +13,8 @@ using SparseArrays
 using PartitionedArrays
 using TimerOutputs
 using WriteVTK
+using AbstractTrees
+
 
 function main(params_in)
     params_default = default_params()
@@ -26,18 +28,15 @@ function main(params_in)
     end
 end
 
-gradient(u) = x->ForwardDiff.gradient(u,x)
-jacobian(u) = x->ForwardDiff.jacobian(u,x)
 laplacian(u,x) = tr(ForwardDiff.jacobian(y->ForwardDiff.gradient(u,y),x))
-laplacian(u) = x->laplacian(u,x)
+# TODO hide this
+laplacian(u::GT.AbstractQuantity,x::GT.AbstractQuantity) = GT.call(laplacian,u,x)
 
-Δ(u) = GT.call(laplacian,u)
-∇(u) = GT.call(gradient,u)
-Δ(u,x) = Δ(u)(x)
-∇(u,x) = ∇(u)(x)
+const ∇ = ForwardDiff.gradient
+const Δ = laplacian
 
-mean(u,x) = 0.5*(u(x)[+]+u(x)[-])
-jump(u,n,x) = u(x)[+]*n[+](x) + u(x)[-]*n[-](x)
+mean(u) = 0.5*(u[1]+u[2])
+jump(u,n) = u[2]*n[2] + u[1]*n[1]
 
 function main_automatic(params)
     timer = params[:timer]
@@ -54,7 +53,7 @@ function main_automatic(params)
 
     u = GT.analytical_field(params[:u],Ω)
     f(x) = -Δ(u,x)
-    n_Γn = GT.unit_normal(Γn,Ω)
+    n_Γn = GT.unit_normal(mesh,D-1)
     g(x) = n_Γn(x)⋅∇(u,x)
 
     interpolation_degree = params[:interpolation_degree]
@@ -68,7 +67,7 @@ function main_automatic(params)
         GT.label_interior_faces!(mesh;physical_name="__INTERIOR_FACES__")
         Λ = GT.skeleton(mesh;physical_names=["__INTERIOR_FACES__"])
         dΛ = GT.measure(Λ,integration_degree)
-        n_Λ = GT.unit_normal(Λ,Ω)
+        n_Λ = GT.unit_normal(mesh,D-1)
         h_Λ = GT.face_diameter_field(Λ)
     else
         conformity = :default
@@ -79,7 +78,7 @@ function main_automatic(params)
         uhd = GT.dirichlet_field(Float64,V)
         GT.interpolate_dirichlet!(u,uhd)
     else
-        n_Γd = GT.unit_normal(Γd,Ω)
+        n_Γd = GT.unit_normal(mesh,D-1)
         h_Γd = GT.face_diameter_field(Γd)
         dΓd = GT.measure(Γd,integration_degree)
         V = GT.lagrange_space(Ω,interpolation_degree;conformity)
@@ -98,7 +97,7 @@ function main_automatic(params)
         end
         if params[:discretization_method] === :interior_penalty
             r += ∫( x->
-                   (γ/h_Λ(x))*jump(v,n_Λ,x)⋅jump(u,n_Λ,x)-jump(v,n_Λ,x)⋅mean(∇(u),x)-mean(∇(v),x)⋅jump(u,n_Λ,x), dΛ)
+                   (γ/h_Λ(x))*jump(v(x),n_Λ(x))⋅jump(u(x),n_Λ(x))-jump(v(x),n_Λ(x))⋅mean(∇(u,x))-mean(∇(v,x))⋅jump(u(x),n_Λ(x)), dΛ)
         end
         r
     end
@@ -202,7 +201,7 @@ function main_hand_written(params)
     dΩ = GT.measure(Ω,integration_degree)
 
     u = params[:u]
-    f(x) = -laplacian(u,x)
+    f(x) = -Δ(u,x)
 
     # Set interpolation space -> Lagrange polynomial basis function, the element type and geometry
     V = GT.lagrange_space(Ω,interpolation_degree;dirichlet_boundary=Γd)
@@ -564,8 +563,8 @@ function default_params()
     params[:dirichlet_tags] = ["boundary"]
     params[:neumann_tags] = String[]
     params[:dirichlet_method] = :strong
-    params[:integration_degree] = 1
-    params[:interpolation_degree] = 2*params[:integration_degree]
+    params[:interpolation_degree] = 1
+    params[:integration_degree] = 2*params[:interpolation_degree]
     params[:solver] = PS.LinearAlgebra_lu
     params[:example_path] = joinpath(mkpath(joinpath(@__DIR__,"..","output")),"example_001")
     params[:export_vtu] = true

GalerkinToolkitExamples/test/poisson_tests.jl

@@ -19,7 +19,7 @@ for discretization_method in (:continuous_galerkin,)
         params[:dirichlet_tags] = ["1-face-1","1-face-2","1-face-3","1-face-4"]
         params[:discretization_method] = discretization_method
         params[:dirichlet_method] = dirichlet_method
-        params[:integration_degree] = 2
+        params[:interpolation_degree] = 2
         results = Poisson.main(params)
         @test results[:error_l2_norm] < tol
         @test results[:error_h1_norm] < tol
@@ -35,7 +35,7 @@ for discretization_method in (:interior_penalty,:continuous_galerkin)
         params[:neumann_tags] = ["1-face-2"]
         params[:discretization_method] = discretization_method
         params[:dirichlet_method] = dirichlet_method
-        params[:integration_degree] = 2
+        params[:interpolation_degree] = 2
         results = Poisson.main(params)
         @test results[:error_l2_norm] < tol
         @test results[:error_h1_norm] < tol

docs/src/examples/methods.jl

@@ -11,21 +11,22 @@ import FileIO # hide
 domain = (0,1,0,1)
 cells = (10,10)
 mesh = GT.cartesian_mesh(domain,cells)
+D = GT.num_dims(mesh)
 GT.label_boundary_faces!(mesh;physical_name="Γd")
 Ω = GT.interior(mesh)
 Γd = GT.boundary(mesh;physical_names=["Γd"])
-n_Γd = GT.unit_normal(Γd,Ω)
+n_Γd = GT.unit_normal(mesh,D-1)
 u = GT.analytical_field(x->sum(x),Ω)
 #TODO these two definitions should be enough
-#∇ = ForwardDiff.gradient
-#Δ = (u,x) -> tr(ForwardDiff.jacobian(y->∇(u,y),x))
-gradient(u) = x->ForwardDiff.gradient(u,x)
-jacobian(u) = x->ForwardDiff.jacobian(u,x)
-laplacian(u) = x-> tr(ForwardDiff.jacobian(y->ForwardDiff.gradient(u,y),x))
-Δ(u) = GT.call(laplacian,u)
-∇(u) = GT.call(gradient,u)
-Δ(u,x) = Δ(u)(x)
-∇(u,x) = ∇(u)(x)
+∇ = ForwardDiff.gradient
+Δ = (u,x) -> tr(ForwardDiff.jacobian(y->∇(u,y),x))
+#gradient(u) = x->ForwardDiff.gradient(u,x)
+#jacobian(u) = x->ForwardDiff.jacobian(u,x)
+#laplacian(u) = x-> tr(ForwardDiff.jacobian(y->ForwardDiff.gradient(u,y),x))
+#Δ(u) = GT.call(laplacian,u)
+#∇(u) = GT.call(gradient,u)
+#Δ(u,x) = Δ(u)(x)
+#∇(u,x) = ∇(u)(x)
 tol = 1e-10
 order = 2
 γ = order*(order+1)/10
@@ -54,7 +55,7 @@ el2 = sqrt(sum(GT.∫( x->abs2(eh(x)), dΩ)))
 # ## Nitche Method
 
 V = GT.lagrange_space(Ω,order)
-n_Γd = GT.unit_normal(Γd,Ω)
+n_Γd = GT.unit_normal(mesh,D-1)
 h_Γd = GT.face_diameter_field(Γd)
 dΓd = GT.measure(Γd,2*order)
 a_Γd = (u,v) -> GT.∫(
@@ -90,18 +91,18 @@ el2 = sqrt(sum(GT.∫( x->abs2(eh(x)), dΩ)))
 
 # ## Interior Penalty
 
-mean(u,x) = 0.5*(u(x)[+]+u(x)[-])
-jump(u,n,x) = u(x)[+]*n[+](x) + u(x)[-]*n[-](x)
+mean(u) = 0.5*(u[1]+u[2])
+jump(u,n) = u[2]*n[2] + u[1]*n[1]
 GT.label_interior_faces!(mesh;physical_name="Λ")
 Λ = GT.skeleton(mesh;physical_names=["Λ"])
 dΛ = GT.measure(Λ,2*order)
-n_Λ = GT.unit_normal(Λ,Ω)
+n_Λ = GT.unit_normal(mesh,D-1)
 h_Λ = GT.face_diameter_field(Λ)
 V = GT.lagrange_space(Ω,order;conformity=:L2)
 a_Λ = (u,v) -> GT.∫(
-    x-> (γ/h_Λ(x))*jump(v,n_Λ,x)⋅jump(u,n_Λ,x)-
-    jump(v,n_Λ,x)⋅mean(∇(u),x)-
-    mean(∇(v),x)⋅jump(u,n_Λ,x), dΛ)
+                    x-> (γ/h_Λ(x))*jump(v(x),n_Λ(x))⋅jump(u(x),n_Λ(x))-
+                    jump(v(x),n_Λ(x))⋅mean(∇(u,x))-
+                    mean(∇(v,x))⋅jump(u(x),n_Λ(x)), dΛ)
 a_ip = (u,v) -> a_Ω(u,v) + a_Γd(u,v) + a_Λ(u,v)
 l_ip = l_ni
 p = GT.linear_problem(Float64,V,a_ip,l_ip)