From c4ccd9e88fbc00c50e93d2a8dfa63440035a5e2f Mon Sep 17 00:00:00 2001
From: Stephan Hilb <stephan@ecshi.net>
Date: Mon, 5 Jul 2021 18:06:18 +0200
Subject: [PATCH] add vtk export support for DP1
---
src/function.jl | 7 ++++++-
src/mesh.jl | 6 +++---
src/run.jl | 2 +-
3 files changed, 10 insertions(+), 5 deletions(-)
diff --git a/src/function.jl b/src/function.jl
index 0e547e0..9cbd6e6 100644
--- a/src/function.jl
+++ b/src/function.jl
@@ -186,11 +186,16 @@ end
append_data!(vtkfile, f::FeFunction) = append_data!(vtkfile, f, f.space.element)
function append_data!(vtkfile, f::FeFunction, ::P1)
- vtk_point_data(vtkfile, f.data, f.name)
+ # separate out data per cell
+ fd = vec(f.data[f.space.dofmap])
+ vtk_point_data(vtkfile, fd, f.name)
end
function append_data!(vtkfile, f::FeFunction, ::DP0)
vtk_cell_data(vtkfile, f.data, f.name)
end
+append_data!(vtkfile, f::FeFunction, ::DP1) =
+ append_data!(vtkfile, f, P1())
+
diff --git a/src/mesh.jl b/src/mesh.jl
index fb54af3..3835b17 100644
--- a/src/mesh.jl
+++ b/src/mesh.jl
@@ -47,9 +47,9 @@ save(filename, f::FeFunction, fs::FeFunction...) =
save(filename, f.space.mesh, f, fs...)
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)
+ cells = [MeshCell(VTKCellTypes.VTK_TRIANGLE, 3*(i-1)+1:3*(i-1)+3) for i in axes(mesh.cells, 2)]
+ vertices = reshape(mesh.vertices[:, mesh.cells], 2, :)
+ vtkfile = vtk_grid(filename, vertices, cells)
for f in fs
f.space.mesh == mesh ||
throw(ArgumentError("meshes do not match"))
diff --git a/src/run.jl b/src/run.jl
index 37a20ab..47a49c7 100644
--- a/src/run.jl
+++ b/src/run.jl
@@ -123,5 +123,5 @@ function myrun()
end
end
- save("test.vtu", g, u, p1)
+ save("test.vtu", g, u, p1, p2)
end
--
GitLab