implement basic vtk export

.gitignore

+*.vtu

src/SemiSmoothNewton.jl

@@ -2,11 +2,7 @@ module SemiSmoothNewton
 
 include("image.jl")
 include("mesh.jl")
-
-function run()
-
-end
-
+include("function.jl")
 
 function clip_line(r0, r1, l0, l1)
     d_l = -(l1[1] - l0[1])
@@ -40,15 +36,6 @@ function clip_line(r0, r1, l0, l1)
     return c0, c1
 end
 
-
-
-# meshfunction
-#   mesh
-
-# mesh
-#   vertices: idx -> vertex
-#   cells: idx -> cell
-
 "project array data onto a grid function"
 function project_array!(u, img::Array)
     vertices = u.mesh.vertices

src/function.jl

+struct P1Mapper
+    mesh
+    ndofs::Int
+    map::Array{Int, 2}
+end
+
+function P1Mapper(mesh)
+    # special case for P1: dofs correspond to vertices
+    ndofs = size(mesh.vertices, 2)
+    map = copy(mesh.cells)
+    return P1Mapper(mesh, ndofs, map)
+end
+
+struct P1Function
+    mapper::P1Mapper
+    data::Vector{Float64}
+end
+
+function P1Function(mapper)
+    data = Vector{Float64}(undef, mapper.ndofs)
+    return P1Function(mapper, data)
+end
+
+function interpolate!(dst::P1Function, src::Function)
+    mapper = dst.mapper
+    mesh = mapper.mesh
+    for cid in axes(mesh.cells, 2)
+	for ldof in axes(mapper.map, 1)
+	    # TODO: how do ldofs generically map to position?
+	    xid = mesh.cells[ldof, cid]
+	    x = mesh.vertices[:, xid]
+
+	    # currently scalar
+	    fx = src(x)
+
+	    gdof = mapper.map[ldof, cid]
+	    dst.data[gdof] = fx
+	end
+    end
+end

src/mesh.jl

@@ -3,24 +3,40 @@ using WriteVTK
 # 2d, simplex grid
 struct Mesh
     vertices::Array{Float64, 2}
-    cells::Array{NTuple{3, Int64}}
+    cells::Array{Int, 2}
 end
 
+function init_grid(n)
+    r = LinRange(0., 1., n + 1)
+    coords = collect(Iterators.product(r, r))
 
+    vertices = [x[i] for i in 1:2, x in vec(coords)]
+    cells = Array{Int, 2}(undef, 3, 2 * n^2)
 
+    vidx = LinearIndices((n+1, n+1))
+    e1 = CartesianIndex(1, 0)
+    e2 = CartesianIndex(0, 1)
+    k = 0
+    for I in CartesianIndices((n, n))
+	cells[:, k += 1] .= (vidx[I], vidx[I + e1], vidx[I + e1 + e2])
+	cells[:, k += 1] .= (vidx[I], vidx[I + e1 + e2], vidx[I + e2])
+    end
+
+    return Mesh(vertices, cells)
+end
 
-#function append_data!(vtkfile, f::GlobalMeshFunction{<:TriangleGrid, n}) where n
-#    fgrid = mesh(f)
-#    fdata = Array{Float64, 2}(undef, n, size(fgrid.vertices, 2))
-#    for i in axes(fdata, 2)
-#        fdata[:, i] .= f.f(fgrid.vertices[:, i])
-#    end
-#    vtk_point_data(vtkfile, fdata, "f")
-#end
-#
-#cells = [MeshCell(VTKCellTypes.VTK_TRIANGLE, grid.cells[:, i]) for i in axes(grid.cells, 2)]
-#vtkfile = vtk_grid(filename, grid.vertices, cells)
-#for (i, f) in enumerate(fs)
-#    append_data!(vtkfile, f)
-#end
-#vtk_save(vtkfile)
+function append_data!(vtkfile, f::P1Function)
+    # FIXME: vector-valued data
+    fdata = reshape(f.data, 1, :)
+    vtk_point_data(vtkfile, fdata, "f")
+end
+
+function save(filename, mesh::Mesh, fs...)
+    cells = [MeshCell(VTKCellTypes.VTK_TRIANGLE, mesh.cells[:, i]) for i in axes(mesh.cells, 2)]
+    #vertices = [mesh.vertices[i, j] for i in axes(mesh.vertices, 1), j in axes(mesh.vertices, 2)]
+    vtkfile = vtk_grid(filename, mesh.vertices, cells)
+    for f in fs
+        append_data!(vtkfile, f)
+    end
+    vtk_save(vtkfile)
+end