[Getfem-commits] [getfem-commits] branch master updated: Fix empty response error in exported vtu files

[Konstantinos Poulios]

This is an automated email from the git hooks/post-receive script.

logari81 pushed a commit to branch master
in repository getfem.

The following commit(s) were added to refs/heads/master by this push:
     new 6af44d9  Fix empty response error in exported vtu files
6af44d9 is described below

commit 6af44d98ceb7ca69edf7ecb430ca84fd2d8f4e9a
Author: Tetsuo Koyama <tkoyama010@gmail.com>
AuthorDate: Fri Jan 8 12:17:28 2021 +0100

    Fix empty response error in exported vtu files
---
 ...ty_finite_strain_linear_hardening_tension_3D.py |   5 +-
 ...strain_linear_hardening_tension_axisymmetric.py |   5 +-
 ...strain_linear_hardening_tension_plane_strain.py |   5 +-
 interface/tests/python/check_export_vtu.py         | 120 +++++++++++++++++----
 src/getfem/getfem_export.h                         |   4 +-
 src/getfem_export.cc                               | 109 ++++++++++---------
 6 files changed, 165 insertions(+), 83 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..431131e 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)
@@ -166,7 +166,6 @@ md.add_initialized_data("disp", [0.])
 
 md.add_initialized_data("K", E/(3.*(1.-2.*nu))) # Bulk modulus
 md.add_initialized_data("mu", E/(2*(1+nu)))     # Shear modulus
-md.add_initialized_data("Y", np.sqrt(2./3.)*pl_sigma_0) # Initial yield limit
 md.add_macro("F", "Id(3)+Grad_u")
 md.add_macro("J", "Det(F)")
 md.add_macro("tauH", "K*log(J)")
@@ -238,7 +237,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..7969b23 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)
@@ -156,7 +156,6 @@ md.add_initialized_data("disp", [0.])
 
 md.add_initialized_data("K", E/(3.*(1.-2.*nu))) # Bulk modulus
 md.add_initialized_data("mu", E/(2*(1+nu)))     # Shear modulus
-md.add_initialized_data("Y", np.sqrt(2./3.)*pl_sigma_0) # Initial yield limit
 md.add_macro("F", "Id(2)+Grad_u")
 md.add_macro("F3d", 
"Id(3)+[0,0,0;0,0,0;0,0,1/X(2)]*u(2)+[1,0;0,1;0,0]*Grad_u*[1,0,0;0,1,0]")
 md.add_macro("J", "Det(F)*(1+u(2)/X(2))")
@@ -227,7 +226,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..a80d325 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)
@@ -155,7 +155,6 @@ md.add_initialized_data("disp", [0.])
 
 md.add_initialized_data("K", E/(3.*(1.-2.*nu))) # Bulk modulus
 md.add_initialized_data("mu", E/(2*(1+nu)))     # Shear modulus
-md.add_initialized_data("Y", np.sqrt(2./3.)*pl_sigma_0) # Initial yield limit
 md.add_macro("F", "Id(2)+Grad_u")
 md.add_macro("F3d", "Id(3)+[1,0;0,1;0,0]*Grad_u*[1,0,0;0,1,0]")
 md.add_macro("J", "Det(F)")
@@ -225,7 +224,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..0176288 100644
--- a/src/getfem_export.cc
+++ b/src/getfem_export.cc
@@ -433,8 +433,8 @@ namespace getfem
     /* count total number of simplexes, and total number of entries */
     size_type cells_cnt = 0, splx_cnt = 0;
     for (size_type ic=0; ic < psl->nb_convex(); ++ic) {
-      for (size_type i=0; i < psl->simplexes(ic).size(); ++i)
-       cells_cnt += psl->simplexes(ic)[i].dim() + 2;
+      for (const slice_simplex &s : psl->simplexes(ic))
+       cells_cnt += s.dim() + 2;
       splx_cnt += psl->simplexes(ic).size();
     }
     if (vtk) {
@@ -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,15 +474,21 @@ 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) {
+          for (const slice_simplex &s : psl->simplexes(ic))
+            size += int((s.dim()+1)*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 (const auto &val : s) {
-        if (vtk) {
-          write_val(int(val.dim()+1));
-        }
-        for (size_type j=0; j < val.dim()+1; ++j)
-          write_val(int(val.inodes[j] + nodes_cnt));
+      for (const slice_simplex &s : psl->simplexes(ic)) {
+        if (vtk)
+          write_val(int(s.dim()+1));
+        for (size_type j=0; j < s.dim()+1; ++j)
+          write_val(int(s.inodes[j] + nodes_cnt));
         write_separ();
       }
       nodes_cnt += psl->nodes(ic).size();
@@ -494,19 +501,24 @@ 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)
+          size += int(psl->simplexes(ic).size()*sizeof(int));
+        write_val(size);
+      }
     }
     int cnt = 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 (const slice_simplex &s : psl->simplexes(ic))
         if (vtk) {
-          write_val(int(vtk_simplex_code[val.dim()]));
+          write_val(int(vtk_simplex_code[s.dim()]));
         } else {
-          cnt += int(val.dim()+1);
+          cnt += int(s.dim()+1);
           write_val(cnt);
         }
       write_separ();
-      splx_cnt -= s.size();
+      splx_cnt -= psl->simplexes(ic).size();
     }
     write_vals();
     assert(splx_cnt == 0); // sanity check
@@ -514,14 +526,17 @@ namespace getfem
       os << (ascii ? "" : "\n") << "</DataArray>\n";
       os << "<DataArray type=\"Int32\" Name=\"types\" ";
       os << (ascii ? "format=\"ascii\">\n" : "format=\"binary\">\n");
-      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) {
-          write_val(int(vtk_simplex_code[s[i].dim()]));
-        }
+      if (!ascii) {
+        int size = 0;
+        for (size_type ic=0; ic < psl->nb_convex(); ++ic)
+          size += int(psl->simplexes(ic).size()*sizeof(int));
+        write_val(size);
       }
+      for (size_type ic=0; ic < psl->nb_convex(); ++ic)
+        for (const slice_simplex &s : psl->simplexes(ic))
+          write_val(int(vtk_simplex_code[s.dim()]));
       write_vals();
-      os << (ascii ? "" : "\n") << "</DataArray>\n";
+      os << "\n" << "</DataArray>\n";
       os << "</Cells>\n";
     }
     state = STRUCTURE_WRITTEN;
@@ -542,11 +557,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 +582,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,12 +608,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;
-        for (dal::bv_visitor cv(pmf->convex_index()); !cv.finished(); ++cv)
-          size += int(sizeof(int));
-        write_val(size);
-      }
+      if (!ascii)
+        write_val(int(pmf->convex_index().card()*sizeof(int)));
       cnt = 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]);
@@ -610,12 +620,8 @@ namespace getfem
       os << "\n" << "</DataArray>\n";
       os << "<DataArray type=\"Int32\" Name=\"types\" ";
       os << (ascii ? "format=\"ascii\">\n" : "format=\"binary\">\n");
-      if (!ascii) {
-        size = 0;
-        for (dal::bv_visitor cv(pmf->convex_index()); !cv.finished(); ++cv)
-          size += int(sizeof(int));
-        write_val(size);
-      }
+      if (!ascii)
+        write_val(int(pmf->convex_index().card()*sizeof(int)));
     }
     for (dal::bv_visitor cv(pmf->convex_index()); !cv.finished(); ++cv) {
       write_val(int(select_vtk_type(pmf_mapping_type[cv])));
@@ -964,16 +970,15 @@ namespace getfem
 
     size_type nodes_cnt = 0; /* <- a virer , global_index le remplace */
     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) {
-        if (s[i].dim() == connections_dim) {
-          for (size_type j=0; j < s[i].dim()+1; ++j) {
-           size_type k;
-           if (psl_use_merged)
-             k = psl->merged_index(ic, s[i].inodes[j]);
-           else k = psl->global_index(ic, s[i].inodes[j]);
-            write_val(int(k));
-         }
+      for (const slice_simplex &s : psl->simplexes(ic)) {
+        if (s.dim() == connections_dim) {
+          for (size_type j=0; j < s.dim()+1; ++j) {
+            size_type k;
+            if (psl_use_merged)
+              k = psl->merged_index(ic, s.inodes[j]);
+            else k = psl->global_index(ic, s.inodes[j]);
+              write_val(int(k));
+          }
           write_separ();
         }
       }
@@ -1202,10 +1207,9 @@ namespace getfem
                 << int(dim) << "D slice (not supported)");
 
     for (size_type ic=0, pcnt=0; ic < psl->nb_convex(); ++ic) {
-      for (getfem::mesh_slicer::cs_simplexes_ct::const_iterator 
it=psl->simplexes(ic).begin();
-           it != psl->simplexes(ic).end(); ++it) {
+      for (const slice_simplex &s : psl->simplexes(ic)) {
         int t = -1;
-        switch ((*it).dim()){
+        switch (s.dim()){
           case 0: t = POS_PT; break;
           case 1: t = POS_LN; break;
           case 2: t = POS_TR; break;
@@ -1220,15 +1224,14 @@ namespace getfem
         std::vector<unsigned> cell_dof;
         cell_dof.resize(dmap.size(),unsigned(-1));
         for (size_type i=0; i < dmap.size(); ++i)
-          cell_dof[i] = unsigned(it->inodes[dmap[i]] + pcnt);
+          cell_dof[i] = unsigned(s.inodes[dmap[i]] + pcnt);
         pos_cell_dof.push_back(cell_dof);
       }
-      for(getfem::mesh_slicer::cs_nodes_ct::const_iterator 
it=psl->nodes(ic).begin();
-          it != psl->nodes(ic).end(); ++it) {
+      for (const slice_node &n : psl->nodes(ic)) {
         std::vector<float> pt;
         pt.resize(dim,float(0));
         for (size_type i=0; i<dim; ++i)
-          pt[i] = float(it->pt[i]);
+          pt[i] = float(n.pt[i]);
         pos_pts.push_back(pt);
       }
       pcnt += psl->nodes(ic).size();