implement basic operator assembly

Manifest.toml

@@ -37,22 +37,6 @@ version = "3.31.0"
 deps = ["Artifacts", "Libdl"]
 uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
 
-[[Crayons]]
-git-tree-sha1 = "3f71217b538d7aaee0b69ab47d9b7724ca8afa0d"
-uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
-version = "4.0.4"
-
-[[DataAPI]]
-git-tree-sha1 = "ee400abb2298bd13bfc3df1c412ed228061a2385"
-uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
-version = "1.7.0"
-
-[[DataStructures]]
-deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
-git-tree-sha1 = "4437b64df1e0adccc3e5d1adbc3ac741095e4677"
-uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
-version = "0.18.9"
-
 [[Dates]]
 deps = ["Printf"]
 uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"
@@ -87,24 +71,6 @@ version = "0.8.5"
 deps = ["ArgTools", "LibCURL", "NetworkOptions"]
 uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
 
-[[FastGaussQuadrature]]
-deps = ["LinearAlgebra", "SpecialFunctions", "StaticArrays"]
-git-tree-sha1 = "5829b25887e53fb6730a9df2ff89ed24baa6abf6"
-uuid = "442a2c76-b920-505d-bb47-c5924d526838"
-version = "0.4.7"
-
-[[Ferrite]]
-deps = ["LinearAlgebra", "Reexport", "SparseArrays", "Tensors", "WriteVTK"]
-path = "/home/stev47/.julia/dev/Ferrite"
-uuid = "c061ca5d-56c9-439f-9c0e-210fe06d3992"
-version = "0.3.0"
-
-[[Fewaku]]
-deps = ["FastGaussQuadrature", "ForwardDiff", "LinearAlgebra", "LinearMaps", "StaticArrays", "WriteVTK"]
-path = "/home/stev47/stuff/Fewaku.jl"
-uuid = "eaeb555d-1c71-5d84-9937-0fdf4ae0c3db"
-version = "0.1.0"
-
 [[FillArrays]]
 deps = ["LinearAlgebra", "Random", "SparseArrays"]
 git-tree-sha1 = "31939159aeb8ffad1d4d8ee44d07f8558273120a"
@@ -162,12 +128,6 @@ version = "0.9.0"
 deps = ["Libdl"]
 uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 
-[[LinearMaps]]
-deps = ["LinearAlgebra", "SparseArrays"]
-git-tree-sha1 = "ebec45cc8cfe3c947b01c2c491e7e1eefba73c4f"
-uuid = "7a12625a-238d-50fd-b39a-03d52299707e"
-version = "3.3.0"
-
 [[LogExpFunctions]]
 deps = ["DocStringExtensions", "LinearAlgebra"]
 git-tree-sha1 = "1ba664552f1ef15325e68dc4c05c3ef8c2d5d885"
@@ -191,12 +151,6 @@ uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
 deps = ["Artifacts", "Libdl"]
 uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"
 
-[[Missings]]
-deps = ["DataAPI"]
-git-tree-sha1 = "4ea90bd5d3985ae1f9a908bd4500ae88921c5ce7"
-uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
-version = "1.0.0"
-
 [[Mmap]]
 uuid = "a63ad114-7e13-5084-954f-fe012c677804"
 
@@ -217,11 +171,6 @@ git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1"
 uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
 version = "0.5.5+0"
 
-[[OrderedCollections]]
-git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
-uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
-version = "1.4.1"
-
 [[Pkg]]
 deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
 uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
@@ -244,19 +193,9 @@ uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
 deps = ["Serialization"]
 uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 
-[[Reexport]]
-git-tree-sha1 = "5f6c21241f0f655da3952fd60aa18477cf96c220"
-uuid = "189a3867-3050-52da-a836-e630ba90ab69"
-version = "1.1.0"
-
 [[SHA]]
 uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
 
-[[SIMD]]
-git-tree-sha1 = "9ba33637b24341aba594a2783a502760aa0bff04"
-uuid = "fdea26ae-647d-5447-a871-4b548cad5224"
-version = "3.3.1"
-
 [[Serialization]]
 uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
 
@@ -267,12 +206,6 @@ uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"
 [[Sockets]]
 uuid = "6462fe0b-24de-5631-8697-dd941f90decc"
 
-[[SortingAlgorithms]]
-deps = ["DataStructures"]
-git-tree-sha1 = "2ec1962eba973f383239da22e75218565c390a96"
-uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
-version = "1.0.0"
-
 [[SparseArrays]]
 deps = ["LinearAlgebra", "Random"]
 uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
@@ -293,17 +226,6 @@ version = "1.2.4"
 deps = ["LinearAlgebra", "SparseArrays"]
 uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
-[[StatsAPI]]
-git-tree-sha1 = "1958272568dc176a1d881acb797beb909c785510"
-uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
-version = "1.0.0"
-
-[[StatsBase]]
-deps = ["DataAPI", "DataStructures", "LinearAlgebra", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "StatsAPI"]
-git-tree-sha1 = "2f6792d523d7448bbe2fec99eca9218f06cc746d"
-uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
-version = "0.33.8"
-
 [[TOML]]
 deps = ["Dates"]
 uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
@@ -312,12 +234,6 @@ uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
 deps = ["ArgTools", "SHA"]
 uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"
 
-[[Tensors]]
-deps = ["ForwardDiff", "LinearAlgebra", "SIMD", "Statistics"]
-git-tree-sha1 = "09a1b69613b0b5014a02a9ab619fa722417e5c3a"
-uuid = "48a634ad-e948-5137-8d70-aa71f2a747f4"
-version = "1.5.0"
-
 [[Test]]
 deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
 uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
@@ -335,12 +251,6 @@ uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
 [[Unicode]]
 uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
 
-[[UnicodePlots]]
-deps = ["Crayons", "Dates", "SparseArrays", "StatsBase"]
-git-tree-sha1 = "1a63e6eea76b291378ff9f95801f8b6d96213208"
-uuid = "b8865327-cd53-5732-bb35-84acbb429228"
-version = "1.3.0"
-
 [[WriteVTK]]
 deps = ["Base64", "CodecZlib", "FillArrays", "LightXML", "Random", "TranscodingStreams"]
 git-tree-sha1 = "2d559f58ff021c0f0f3d2394d2688a84a84257aa"

Project.toml

@@ -4,9 +4,8 @@ authors = ["Stephan Hilb <stephan@ecshi.net>"]
 version = "0.1.0"
 
 [deps]
-Ferrite = "c061ca5d-56c9-439f-9c0e-210fe06d3992"
-Fewaku = "eaeb555d-1c71-5d84-9937-0fdf4ae0c3db"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
-UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228"
 WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192"

src/SemiSmoothNewton.jl

@@ -3,6 +3,7 @@ module SemiSmoothNewton
 include("image.jl")
 include("mesh.jl")
 include("function.jl")
+include("operator.jl")
 
 function clip_line(r0, r1, l0, l1)
     d_l = -(l1[1] - l0[1])

src/function.jl

@@ -6,35 +6,21 @@ struct P1 end
 struct DP0 end
 struct DP1 end
 
+
 # Fe: finite element type
-struct Mapper{Fe}
-    mesh
+struct Mapper{M,Fe}
+    mesh::M
     ndofs::Int
     map::Array{Int, 2} # (ldof, cid) -> gdof
 end
 
-# dof ordering for vector valued functions:
-#   (ldof, fdims...)
-
-struct FeFunction
-    mapper::Mapper
-    size::NTuple
-    data::Array{Float64, 2} # (rdim, gdof) -> data
-    name::String
-end
-
-function FeFunction(mapper, size=(1,); name=string(gensym("f")))
-    data = Array{Float64, 2}(undef, prod(size), mapper.ndofs)
-    return FeFunction(mapper, size, data, name)
-end
-
 # P1 Elements (1st order Lagrange)
 
 function Mapper{P1}(mesh)
     # special case for P1: dofs correspond to vertices
     ndofs = size(mesh.vertices, 2)
     map = copy(mesh.cells)
-    return Mapper{P1}(mesh, ndofs, map)
+    return Mapper{typeof(mesh),P1}(mesh, ndofs, map)
 end
 
 # DP0 Elements (0th order discontinuous Lagrange)
@@ -43,12 +29,46 @@ function Mapper{DP0}(mesh)
     # special case for DP0: dofs correspond to cells
     ndofs = size(mesh.cells, 2)
     map = collect(reshape(axes(mesh.cells, 2), 1, :))
-    return Mapper{DP0}(mesh, ndofs, map)
+    return Mapper{typeof(mesh),DP0}(mesh, ndofs, map)
+end
+
+# DP1 Elements (1st order discontinuous Lagrange)
+
+function Mapper{DP1}(mesh)
+    # special case for DP0: dofs correspond to cells
+    ndofs = size(mesh.cells, 2)
+    map = collect(reshape(axes(mesh.cells, 2), 1, :))
+    return Mapper{typeof(mesh),DP0}(mesh, ndofs, map)
 end
 
+Base.show(io::IO, ::MIME"text/plain", m::Mapper) =
+    print("$(nameof(typeof(m))), $(element(m)) elements, $(size(m.map, 1)) local dofs each")
+
+element(::Mapper{<:Any,Fe}) where Fe = Fe
+
+
+
+# dof ordering for vector valued functions:
+#   (ldof, fdims...)
+
+struct FeFunction
+    mapper::Mapper
+    size::NTuple
+    data::Array{Float64, 2} # (rdim, gdof) -> data
+    name::String
+end
+
+function FeFunction(mapper, size=(1,); name=string(gensym("f")))
+    data = Array{Float64, 2}(undef, prod(size), mapper.ndofs)
+    return FeFunction(mapper, size, data, name)
+end
+
+Base.show(io::IO, ::MIME"text/plain", f::FeFunction) =
+    print("$(nameof(typeof(f))), size $(f.size) with $(length(f.data)) dofs")
+
 interpolate!(dst::FeFunction, src::Function) = interpolate!(dst, dst.mapper, src)
 
-function interpolate!(dst::FeFunction, mapper::Mapper{P1}, src::Function)
+function interpolate!(dst::FeFunction, mapper::Mapper{<:Any,P1}, src::Function)
     mesh = mapper.mesh
     for cid in axes(mesh.cells, 2)
 	for ldof in axes(mapper.map, 1)
@@ -62,7 +82,7 @@ function interpolate!(dst::FeFunction, mapper::Mapper{P1}, src::Function)
     end
 end
 
-function interpolate!(dst::FeFunction, mapper::Mapper{DP0}, src::Function)
+function interpolate!(dst::FeFunction, mapper::Mapper{<:Any,DP0}, src::Function)
     mesh = mapper.mesh
     for cid in axes(mesh.cells, 2)
 	vertices = mesh.vertices[:, mesh.cells[:, cid]]
@@ -75,10 +95,12 @@ end
 
 append_data!(vtkfile, f::FeFunction) = append_data!(vtkfile, f, f.mapper)
 
-function append_data!(vtkfile, f::FeFunction, ::Mapper{P1})
+function append_data!(vtkfile, f::FeFunction, ::Mapper{<:Any, P1})
     vtk_point_data(vtkfile, f.data, f.name)
 end
 
-function append_data!(vtkfile, f::FeFunction, ::Mapper{DP0})
+function append_data!(vtkfile, f::FeFunction, ::Mapper{<:Any, DP0})
     vtk_cell_data(vtkfile, f.data, f.name)
 end
+
+

src/mesh.jl

@@ -8,6 +8,9 @@ struct Mesh
     cells::Array{Int, 2}
 end
 
+Base.show(io::IO, ::MIME"text/plain", f::Mesh) =
+    print("$(nameof(typeof(f))), $(size(f.cells, 2)) cells")
+
 function init_grid(n)
     r = LinRange(0., 1., n + 1)
     coords = collect(Iterators.product(r, r))

src/operator.jl

+using SparseArrays: sparse
+using LinearAlgebra: det, dot
+using ForwardDiff: jacobian
+
+export Poisson, L2Projection, init_point!, assemble
+
+abstract type Operator end
+
+struct L2Projection <: Operator
+    f::FeFunction # target
+    # space info
+    mapper::Mapper
+    size::NTuple
+end
+
+L2Projection(f, gspace) = L2Projection(f, g.mapper, g.size)
+
+a(op::L2Projection, xloc, u, du, v, dv) = dot(u, v)
+l(op::L2Projection, xloc, v, dv) = dot(op.f, v)
+
+
+struct Poisson{M,S} <: Operator
+    mapper::M
+    size::S
+end
+
+a(op::Poisson, xloc, u, du, v, dv) = dot(du, dv)
+
+
+
+# degree 2 quadrature on triangle
+quadrature() = [1/6, 1/6, 1/6], [1/6 4/6 1/6; 1/6 1/6 4/6]
+
+# evaluate all (1-dim) local basis functions against x
+# linear 1st order lagrange
+localbasis(x) = [1 - x[1] - x[2], x[1], x[2]]
+
+elmap(mesh, cell, x) =
+    mesh.vertices[:, mesh.cells[:, cell]] * [1 - x[1] - x[2], x[1], x[2]]
+
+
+function assemble(op::Operator)
+    mapper = op.mapper
+    mesh = mapper.mesh
+
+    # precompute basis at quadrature points
+
+    qw, qx = quadrature()
+
+    d = size(qx, 1) # domain dimension
+    nrdims = prod(op.size)
+    nldofs = size(mapper.map, 1) # number of local dofs (not counting range dimensions)
+    nqpts = length(qw) # number of quadrature points
+
+    qphi = zeros(nrdims, nqpts, nrdims, nldofs)
+    dqphi = zeros(nrdims, d, nqpts, nrdims, nldofs)
+    for r in 1:nrdims
+	for k in axes(qx, 2)
+	    qphi[r, k, r, :] .= localbasis(qx[:, k])
+	    dqphi[r, :, k, r, :] .= transpose(jacobian(localbasis, qx[:, k])::Array{Float64,2})
+	end
+    end
+
+    I = Float64[]
+    J = Float64[]
+    V = Float64[]
+    gdof = LinearIndices((nrdims, mapper.ndofs))
+    for cell in axes(mesh.cells, 2)
+	delmap = jacobian(x -> elmap(mesh, cell, x), [0., 0.])::Array{Float64,2} # constant on element
+	delmapinv = inv(delmap) # constant on element
+	intel = abs(det(delmap))
+
+	# local dofs
+	for idim in 1:nrdims, ldofi in 1:nldofs
+	    mdofi = mapper.map[ldofi, cell]
+	    gdofi = gdof[idim, mdofi]
+	    for jdim in 1:nrdims, ldofj in 1:nldofs
+		mdofj = mapper.map[ldofj, cell]
+		gdofj = gdof[jdim, mdofj]
+
+		# quadrature points
+		for k in axes(qx, 2)
+		    xhat = qx[:, k]
+		    phii = reshape(qphi[:, k, idim, ldofi], op.size)
+		    dphii = reshape(dqphi[:, :, k, idim, ldofi] * delmapinv, (op.size..., :))
+		    phij = reshape(qphi[:, k, jdim, ldofj], op.size)
+		    dphij = reshape(dqphi[:, :, k, jdim, ldofj] * delmapinv, (op.size..., :))
+
+		    gdofv = qw[k] * a(op, xhat, phii, dphii, phij, dphij) * intel
+
+		    push!(I, gdofi)
+		    push!(J, gdofj)
+		    push!(V, gdofv)
+
+		    #display(phii)
+		    #display(dphii)
+		    #display(phij)
+		    #display(dphij)
+
+		    #println("$cell, $gdofi, $gdofj, $gdofv")
+		end
+	    end
+	end
+    end
+
+    ngdofs = nrdims * mapper.ndofs
+    A = sparse(I, J, V, ngdofs, ngdofs)
+    return A
+end
+
+