branch: devel-tetsuo-fix-export-vtu commit 2d3812300b9be7d0587dcac30d7feb73b187a126 Author: Tetsuo Koyama <tkoyama...@gmail.com> AuthorDate: Sun Dec 27 20:02:44 2020 +0900
Fix empty response error of export vtu --- ...ty_finite_strain_linear_hardening_tension_3D.py | 4 +- ...strain_linear_hardening_tension_axisymmetric.py | 4 +- ...strain_linear_hardening_tension_plane_strain.py | 4 +- interface/tests/python/check_export_vtu.py | 120 +++++++++++++++++---- src/getfem/getfem_export.h | 4 +- src/getfem_export.cc | 44 ++++++-- 6 files changed, 145 insertions(+), 35 deletions(-) diff --git a/contrib/continuum_mechanics/plasticity_finite_strain_linear_hardening_tension_3D.py b/contrib/continuum_mechanics/plasticity_finite_strain_linear_hardening_tension_3D.py index f041500..bd5eca9 100644 --- a/contrib/continuum_mechanics/plasticity_finite_strain_linear_hardening_tension_3D.py +++ b/contrib/continuum_mechanics/plasticity_finite_strain_linear_hardening_tension_3D.py @@ -128,7 +128,7 @@ if dH > 0: pts[2,i] -= (z*dH)/(2*H) * (1 + np.cos(2.*np.pi*x/L)) mesh.set_pts(pts) -mesh.export_to_vtk("%s/mesh.vtk" % resultspath) +mesh.export_to_vtu("%s/mesh.vtu" % resultspath) # FEM mfN = gf.MeshFem(mesh, N) @@ -238,7 +238,7 @@ with open("%s/tension_3D.dat" % resultspath, "w") as f1: mfu, md.variable("u"), "Displacements", mfout2, md.interpolation("dirmult", mfout2, XP_RG), "Nominal reaction traction", mfout2, md.local_projection(mim, "gamma", mfout2), "plastic strain") - mfout2.export_to_vtk("%s/tension_3D_%i.vtk" % (resultspath, step), *output) + mfout2.export_to_vtu("%s/tension_3D_%i.vtu" % (resultspath, step), *output) md.set_variable("gamma0", md.interpolation("gamma", mimd1, -1)) md.set_variable("invCp0vec", diff --git a/contrib/continuum_mechanics/plasticity_finite_strain_linear_hardening_tension_axisymmetric.py b/contrib/continuum_mechanics/plasticity_finite_strain_linear_hardening_tension_axisymmetric.py index 8fc92ae..5ee76d8 100644 --- a/contrib/continuum_mechanics/plasticity_finite_strain_linear_hardening_tension_axisymmetric.py +++ b/contrib/continuum_mechanics/plasticity_finite_strain_linear_hardening_tension_axisymmetric.py @@ -119,7 +119,7 @@ if dH > 0: pts[1,i] -= (y*dH)/(2*H) * (1 + np.cos(2.*np.pi*x/L)) mesh.set_pts(pts) -mesh.export_to_vtk("%s/mesh.vtk" % resultspath) +mesh.export_to_vtu("%s/mesh.vtu" % resultspath) # FEM mfN = gf.MeshFem(mesh, N) @@ -227,7 +227,7 @@ with open("%s/tension_axisymmetric.dat" % resultspath, "w") as f1: mfu, md.variable("u"), "Displacements", mfout2, md.interpolation("dirmult", mfout2, XP_RG), "Nominal reaction traction", mfout2, md.local_projection(mim, "gamma", mfout2), "plastic strain") - mfout2.export_to_vtk("%s/tension_axisymmetric_%i.vtk" % (resultspath, step), *output) + mfout2.export_to_vtu("%s/tension_axisymmetric_%i.vtu" % (resultspath, step), *output) md.set_variable("gamma0", md.interpolation("gamma", mimd1, -1)) md.set_variable("invCp0vec", diff --git a/contrib/continuum_mechanics/plasticity_finite_strain_linear_hardening_tension_plane_strain.py b/contrib/continuum_mechanics/plasticity_finite_strain_linear_hardening_tension_plane_strain.py index 422efe1..af489fd 100644 --- a/contrib/continuum_mechanics/plasticity_finite_strain_linear_hardening_tension_plane_strain.py +++ b/contrib/continuum_mechanics/plasticity_finite_strain_linear_hardening_tension_plane_strain.py @@ -118,7 +118,7 @@ if dH > 0: pts[1,i] -= (y*dH)/(2*H) * (1 + np.cos(2.*np.pi*x/L)) mesh.set_pts(pts) -mesh.export_to_vtk("%s/mesh.vtk" % resultspath) +mesh.export_to_vtu("%s/mesh.vtu" % resultspath) # FEM mfN = gf.MeshFem(mesh, N) @@ -225,7 +225,7 @@ with open("%s/tension_plane_strain.dat" % resultspath, "w") as f1: mfu, md.variable("u"), "Displacements", mfout2, md.interpolation("dirmult", mfout2, XP_RG), "Nominal reaction traction", mfout2, md.local_projection(mim, "gamma", mfout2), "plastic strain") - mfout2.export_to_vtk("%s/tension_plane_strain_%i.vtk" % (resultspath, step), *output) + mfout2.export_to_vtu("%s/tension_plane_strain_%i.vtu" % (resultspath, step), *output) md.set_variable("gamma0", md.interpolation("gamma", mimd1, -1)) md.set_variable("invCp0vec", diff --git a/interface/tests/python/check_export_vtu.py b/interface/tests/python/check_export_vtu.py index d6c45f9..0f20332 100644 --- a/interface/tests/python/check_export_vtu.py +++ b/interface/tests/python/check_export_vtu.py @@ -2,7 +2,7 @@ # -*- coding: utf-8 -*- # Python GetFEM interface # -# Copyright (C) 2004-2020 Yves Renard, Julien Pommier. +# Copyright (C) 2020-2020 Tetsuo Koyama. # # This file is a part of GetFEM # @@ -29,6 +29,7 @@ import getfem as gf import numpy as np import sys + try: import pyvista as pv except: @@ -41,26 +42,105 @@ mesh = gf.Mesh("cartesian", [0.0, 1.0, 2.0]) pts = mesh.pts()[0] -filenames = [ - "check_mesh_ascii.vtk", - "check_mesh_binary.vtk", - "check_mesh_ascii.vtu", - "check_mesh_binary.vtu", -] +file_name = "check_mesh_ascii.vtk" +mesh.export_to_vtk(file_name, "ascii") +unstructured_grid = pv.read(file_name) +expected = pts +actual = unstructured_grid.points[:, 0] +np.testing.assert_equal(expected, actual, "export of mesh pts is not correct.") +expected = convex_connectivity +actual = unstructured_grid.cell_connectivity +np.testing.assert_equal(expected, actual, "export of mesh convex is not correct.") + +file_name = "check_mesh_binary.vtk" +mesh.export_to_vtk(file_name) +unstructured_grid = pv.read(file_name) +expected = pts +actual = unstructured_grid.points[:, 0] +np.testing.assert_equal(expected, actual, "export of mesh pts is not correct.") +expected = convex_connectivity +actual = unstructured_grid.cell_connectivity +np.testing.assert_equal(expected, actual, "export of mesh convex is not correct.") + +file_name = "check_mesh_ascii.vtu" +mesh.export_to_vtu(file_name, "ascii") +unstructured_grid = pv.read(file_name) +expected = pts +actual = unstructured_grid.points[:, 0] +np.testing.assert_equal(expected, actual, "export of mesh pts is not correct.") +expected = convex_connectivity +actual = unstructured_grid.cell_connectivity +np.testing.assert_equal(expected, actual, "export of mesh convex is not correct.") + +file_name = "check_mesh_binary.vtu" +mesh.export_to_vtu(file_name) +unstructured_grid = pv.read(file_name) +expected = pts +actual = unstructured_grid.points[:, 0] +np.testing.assert_equal(expected, actual, "export of mesh pts is not correct.") +expected = convex_connectivity +actual = unstructured_grid.cell_connectivity +np.testing.assert_equal(expected, actual, "export of mesh convex is not correct.") + +mfu = gf.MeshFem(mesh, 1) +mfu.set_classical_fem(1) +U1 = np.array([2.0, 1.0, 0.0]) + +file_name = "check_meshfem_ascii.vtk" +mfu.export_to_vtk(file_name, "ascii", U1, "U1") +unstructured_grid = pv.read(file_name) +expected = U1 +actual = unstructured_grid.point_arrays["U1"] +np.testing.assert_equal(expected, actual, "export of U1 is not correct.") + +file_name = "check_meshfem_binary.vtk" +mfu.export_to_vtk(file_name, U1, "U1") +unstructured_grid = pv.read(file_name) +expected = U1 +actual = unstructured_grid.point_arrays["U1"] +np.testing.assert_equal(expected, actual, "export of U1 is not correct.") + +file_name = "check_meshfem_ascii.vtu" +mfu.export_to_vtu(file_name, "ascii", U1, "U1") +unstructured_grid = pv.read(file_name) +expected = U1 +actual = unstructured_grid.point_arrays["U1"] +np.testing.assert_equal(expected, actual, "export of U1 is not correct.") + +file_name = "check_meshfem_binary.vtu" +mfu.export_to_vtu(file_name, U1, "U1") +unstructured_grid = pv.read(file_name) +expected = U1 +actual = unstructured_grid.point_arrays["U1"] +np.testing.assert_equal(expected, actual, "export of U1 is not correct.") + +sl = gf.Slice(("boundary",), mesh, 1) +U2 = np.array([3.0, 2.0, 1.0, 0.0]) -mesh.export_to_vtk(filenames[0], "ascii") -mesh.export_to_vtk(filenames[1]) -mesh.export_to_vtu(filenames[2], "ascii") -mesh.export_to_vtu(filenames[3]) +file_name = "check_slice_ascii.vtk" +sl.export_to_vtk(file_name, "ascii", U2, "U2") +unstructured_grid = pv.read(file_name) +expected = U2 +actual = unstructured_grid.point_arrays["U2"] +np.testing.assert_equal(expected, actual, "export of U2 is not correct.") -for filename in filenames: - print(filename) - unstructured_grid = pv.read(filename) +file_name = "check_slice_binary.vtk" +sl.export_to_vtk(file_name, U2, "U2") +unstructured_grid = pv.read(file_name) +expected = U2 +actual = unstructured_grid.point_arrays["U2"] +np.testing.assert_equal(expected, actual, "export of U2 is not correct.") - expected = pts - actual = unstructured_grid.points[:, 0] - np.testing.assert_equal(expected, actual, "export of mesh pts is not correct.") +file_name = "check_slice_ascii.vtu" +sl.export_to_vtu(file_name, "ascii", U2, "U2") +unstructured_grid = pv.read(file_name) +expected = U2 +actual = unstructured_grid.point_arrays["U2"] +np.testing.assert_equal(expected, actual, "export of U2 is not correct.") - expected = convex_connectivity - actual = unstructured_grid.cell_connectivity - np.testing.assert_equal(expected, actual, "export of mesh convex is not correct.") +file_name = "check_slice_binary.vtu" +sl.export_to_vtu(file_name, U2, "U2") +unstructured_grid = pv.read(file_name) +expected = U2 +actual = unstructured_grid.point_arrays["U2"] +np.testing.assert_equal(expected, actual, "export of U2 is not correct.") diff --git a/src/getfem/getfem_export.h b/src/getfem/getfem_export.h index d54c0e5..5d03795 100644 --- a/src/getfem/getfem_export.h +++ b/src/getfem/getfem_export.h @@ -255,6 +255,7 @@ namespace getfem { "inconsistency in the size of the dataset: " << gmm::vect_size(U) << " != " << nb_val << "*" << Q); if (vtk) write_separ(); + if (!vtk && !ascii) write_val(float(gmm::vect_size(U))); if (Q == 1) { if (vtk) os << "SCALARS " << remove_spaces(name) << " float 1\n" @@ -288,8 +289,9 @@ namespace getfem { } else GMM_ASSERT1(false, std::string(vtk ? "vtk" : "vtu") + " does not accept vectors of dimension > 3"); + write_vals(); if (vtk) write_separ(); - if (!vtk) os << (ascii ? "" : "\n") << "</DataArray>\n"; + if (!vtk) os << "\n" << "</DataArray>\n"; } diff --git a/src/getfem_export.cc b/src/getfem_export.cc index 0d3fe84..ab331de 100644 --- a/src/getfem_export.cc +++ b/src/getfem_export.cc @@ -448,6 +448,7 @@ namespace getfem os << "<DataArray type=\"Float32\" Name=\"Points\" "; os << "NumberOfComponents=\"3\" "; os << (ascii ? "format=\"ascii\">\n" : "format=\"binary\">\n"); + if (!ascii) write_val(int(sizeof(float)*psl->nb_points()*3)); } /* points are not merge, vtk is mostly fine with that (except for @@ -473,6 +474,16 @@ namespace getfem os << "<Cells>\n"; os << "<DataArray type=\"Int32\" Name=\"connectivity\" "; os << (ascii ? "format=\"ascii\">\n" : "format=\"binary\">\n"); + if (!ascii) { + int size = 0; + for (size_type ic=0; ic < psl->nb_convex(); ++ic) { + const getfem::mesh_slicer::cs_simplexes_ct& s = psl->simplexes(ic); + for (const auto &val : s) + for (size_type j=0; j < val.dim()+1; ++j) + size += int(sizeof(int)); + } + write_val(size); + } } for (size_type ic=0; ic < psl->nb_convex(); ++ic) { const getfem::mesh_slicer::cs_simplexes_ct& s = psl->simplexes(ic); @@ -494,6 +505,15 @@ namespace getfem os << (ascii ? "" : "\n") << "</DataArray>\n"; os << "<DataArray type=\"Int32\" Name=\"offsets\" "; os << (ascii ? "format=\"ascii\">\n" : "format=\"binary\">\n"); + if (!ascii) { + int size = 0; + for (size_type ic=0; ic < psl->nb_convex(); ++ic) { + const getfem::mesh_slicer::cs_simplexes_ct& s = psl->simplexes(ic); + for (const auto &val : s) + size += int(sizeof(int)); + } + write_val(size); + } } int cnt = 0; for (size_type ic=0; ic < psl->nb_convex(); ++ic) { @@ -514,6 +534,15 @@ namespace getfem os << (ascii ? "" : "\n") << "</DataArray>\n"; os << "<DataArray type=\"Int32\" Name=\"types\" "; os << (ascii ? "format=\"ascii\">\n" : "format=\"binary\">\n"); + if (!ascii) { + int size = 0; + for (size_type ic=0; ic < psl->nb_convex(); ++ic) { + const getfem::mesh_slicer::cs_simplexes_ct& s = psl->simplexes(ic); + for (size_type i=0; i < s.size(); ++i) + size += int(sizeof(int)); + } + write_val(size); + } for (size_type ic=0; ic < psl->nb_convex(); ++ic) { const getfem::mesh_slicer::cs_simplexes_ct& s = psl->simplexes(ic); for (size_type i=0; i < s.size(); ++i) { @@ -521,7 +550,7 @@ namespace getfem } } write_vals(); - os << (ascii ? "" : "\n") << "</DataArray>\n"; + os << "\n" << "</DataArray>\n"; os << "</Cells>\n"; } state = STRUCTURE_WRITTEN; @@ -542,11 +571,10 @@ namespace getfem os << "<DataArray type=\"Float32\" Name=\"Points\" "; os << "NumberOfComponents=\"3\" "; os << (ascii ? "format=\"ascii\">\n" : "format=\"binary\">\n"); + if (!ascii) write_val(int(sizeof(float)*pmf_dof_used.card()*3)); } std::vector<int> dofmap(pmf->nb_dof()); int cnt = 0; - int size = int(sizeof(float)*pmf_dof_used.card()*3); - if (!vtk && !ascii) write_val(size); for (dal::bv_visitor d(pmf_dof_used); !d.finished(); ++d) { dofmap[d] = cnt++; base_node P = pmf->point_of_basic_dof(d); @@ -568,8 +596,8 @@ namespace getfem os << "<Cells>\n"; os << "<DataArray type=\"Int32\" Name=\"connectivity\" "; os << (ascii ? "format=\"ascii\">\n" : "format=\"binary\">\n"); - if (!vtk && !ascii) { - size = 0; + if (!ascii) { + int size = 0; for (dal::bv_visitor cv(pmf->convex_index()); !cv.finished(); ++cv) { const std::vector<unsigned> &dmap = select_vtk_dof_mapping(pmf_mapping_type[cv]); size += int(sizeof(int)*dmap.size()); @@ -594,8 +622,8 @@ namespace getfem os << (ascii ? "" : "\n") << "</DataArray>\n"; os << "<DataArray type=\"Int32\" Name=\"offsets\" "; os << (ascii ? "format=\"ascii\">\n" : "format=\"binary\">\n"); - if (!vtk && !ascii) { - size = 0; + if (!ascii) { + int size = 0; for (dal::bv_visitor cv(pmf->convex_index()); !cv.finished(); ++cv) size += int(sizeof(int)); write_val(size); @@ -611,7 +639,7 @@ namespace getfem os << "<DataArray type=\"Int32\" Name=\"types\" "; os << (ascii ? "format=\"ascii\">\n" : "format=\"binary\">\n"); if (!ascii) { - size = 0; + int size = 0; for (dal::bv_visitor cv(pmf->convex_index()); !cv.finished(); ++cv) size += int(sizeof(int)); write_val(size);