From 6f6fc20a22b044916986881f4a08db7f33eb0b8d Mon Sep 17 00:00:00 2001 From: Stephan Hilb <stephan@ecshi.net> Date: Mon, 5 Jul 2021 18:32:56 +0200 Subject: [PATCH] export time series data and refactor --- .gitignore | 1 + src/function.jl | 19 +++++++++++-------- src/mesh.jl | 18 +++++++++++------- src/run.jl | 48 ++++++++++++++++++++++++++++++++++-------------- 4 files changed, 57 insertions(+), 29 deletions(-) diff --git a/.gitignore b/.gitignore index 8d69778..1fcac93 100644 --- a/.gitignore +++ b/.gitignore @@ -1 +1,2 @@ *.vtu +*.pvd diff --git a/src/function.jl b/src/function.jl index 9cbd6e6..e3e3cf0 100644 --- a/src/function.jl +++ b/src/function.jl @@ -183,19 +183,22 @@ end -append_data!(vtkfile, f::FeFunction) = append_data!(vtkfile, f, f.space.element) +vtk_append!(vtkfile, f::FeFunction, fs::FeFunction...) = + (vtk_append!(vtkfile, f); vtk_append!(vtkfile, fs...)) -function append_data!(vtkfile, f::FeFunction, ::P1) - # separate out data per cell +vtk_append!(vtkfile, f::FeFunction) = + vtk_append!(vtkfile, f, f.space.element) + +function vtk_append!(vtkfile, f::FeFunction, ::P1) + # separate out data per cell since different cells have distinct vertices + # in the output vtk fd = vec(f.data[f.space.dofmap]) vtk_point_data(vtkfile, fd, f.name) end -function append_data!(vtkfile, f::FeFunction, ::DP0) +function vtk_append!(vtkfile, f::FeFunction, ::DP0) vtk_cell_data(vtkfile, f.data, f.name) end -append_data!(vtkfile, f::FeFunction, ::DP1) = - append_data!(vtkfile, f, P1()) - - +vtk_append!(vtkfile, f::FeFunction, ::DP1) = + vtk_append!(vtkfile, f, P1()) diff --git a/src/mesh.jl b/src/mesh.jl index 3835b17..02b2d8f 100644 --- a/src/mesh.jl +++ b/src/mesh.jl @@ -43,17 +43,21 @@ init_grid(img::Array{<:Any, 2}, type=:vertex) = # return all(λ .>= 0) && sum(λ) <= 1 #end -save(filename, f::FeFunction, fs::FeFunction...) = - save(filename, f.space.mesh, f, fs...) +function vtk_mesh(filename, mesh::Mesh) + 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, :) + + return vtk_grid(filename, vertices, cells) +end +"convenience function for saving to vtk" function save(filename, mesh::Mesh, fs...) - 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) + vtk = vtk_mesh(filename, mesh) for f in fs f.space.mesh == mesh || throw(ArgumentError("meshes do not match")) - append_data!(vtkfile, f) + vtk_append!(vtk, f) end - vtk_save(vtkfile) + vtk_save(vtk) end diff --git a/src/run.jl b/src/run.jl index 47a49c7..6b81628 100644 --- a/src/run.jl +++ b/src/run.jl @@ -5,6 +5,7 @@ using LinearAlgebra: norm function myrun() + name = "test" mesh = init_grid(50) d = 2 m = 1 @@ -63,6 +64,21 @@ function myrun() return -p2 + lambda / m2 * (nablau + nabladu) - cond end + function p1_project!() + p1.space.element::DP0 + p1.data .= clamp.(p1.data, -alpha1, alpha1) + end + + function p2_project!() + p2.space.element::DP1 + p2d = reshape(p2.data, m * d * 3, :) # no copy + for i in axes(p2d, 2) + p2in = norm(p2d[:, i]) + if p2in > lambda + p2d[:, i] .*= lambda ./ p2in + end + end + end @inline function du_a(x, du, nabladu, phi, nablaphi; g, u, nablau, p1, p2) m1 = max(gamma1, norm(T(u) - g)) @@ -97,31 +113,35 @@ function myrun() return -aB - p1part - p2part + gpart end + function save_step!(pvd, i) + vtk = vtk_mesh("$(name)_$(lpad(i, 5, '0')).vtu", mesh) + vtk_append!(vtk, g, u, p1, p2) + vtk_save(vtk) + pvd[i] = vtk + end + + pvd = paraview_collection("$(name).pvd") + save_step!(pvd, 0) for i = 1:6 + # solve du A = assemble(Vu, du_a; g, u, nablau, p1, p2) b = assemble_rhs(Vu, du_l; g, u, nablau) - du.data .= A \ b + # solve dp1, dp2 interpolate!(dp1, dp1_update; g, u, p1, du) interpolate!(dp2, dp2_update; u, nablau, p2, du, nabladu) + # newton update u.data .+= du.data p1.data .+= dp1.data p2.data .+= dp2.data - # reproject - p1.space.element::DP0 - p1.data .= clamp.(p1.data, -alpha1, alpha1) - p2.space.element::DP1 - p2d = reshape(p2.data, m * d * 3, :) # no copy - for i in axes(p2d, 2) - p2in = norm(p2d[:, i]) - if p2in > lambda - p2d[:, i] .*= lambda ./ p2in - end - end - end + # reproject p1, p2 + p1_project!() + p2_project!() - save("test.vtu", g, u, p1, p2) + save_step!(pvd, i) + end + vtk_save(pvd) end -- GitLab