PrintVTK/SaveVTK

[NikolayYavich]

Dear All,

I carefully reviewed #797 and #558 before posting this.
My question is how to avoid refinement and extra points when exporting a field
modelled with linear elements to a VTK file.
From the mention posts it follows that setting ref = 1 in PrintVTK and SaveVTK does this.
However, this simple code bellow shows that this is not the case.

#include "mfem.hpp"
#include <fstream>
#include <iostream>

using namespace mfem;
using namespace std;
 
 void main()
 {
    Mesh *mesh = new Mesh(5, 5, 5, Element::TETRAHEDRON, true, 1.0, 1.0, 1.0);
   FiniteElementCollection *fec(new H1_FECollection(1, 3));
   FiniteElementSpace *W_space = new FiniteElementSpace(mesh, fec);

   const int n = mesh->GetNV();
   cout << "n = " << n << "\n";
   Vector pres(n);
   pres = 0.0;
   for (int i = 0; i < n; i++)
     pres(i) = i % 2;

   GridFunction pres_field(W_space, pres.GetData());

   fstream vtkFs("saving_test.vtk", ios::out);
   const int ref = 1;
   mesh->PrintVTK( vtkFs, ref);
   pres_field.SaveVTK( vtkFs, "pressure", ref);

   delete mesh;
   delete fec;
   delete W_space;
}

When I check the generated VTK file I see:

$ grep "POINTS"  saving_test.vtk
POINTS 3000 double

while the correct number of vertices in this example is 216.

Will appreciate any hint,
Thanks,
Nikolay.

[pazner]

Hi @NikolayYavich,

When you set ref = 1, MFEM will output four points per tetrahedron. Since your mesh has 750 tetrahedra, this is why the VTK file has 3000 points. The reason why there are more points than vertices in the mesh is because of duplicated points (and not because of extra refinement points).

Hope this helps.

[NikolayYavich]

Hi @pazner ,

Thanks a lot for the answer.
So, how to avoid duplicating points when exporting modelled fields?
This is a big obstacle when dealing with larger problems - files get even larger, Paraview hangs...

[pazner]

There is not currently a way to output a GridFunction in VTK format without duplicated points.

If the file size is the main concern, you can try using the ParaViewDataCollection with binary compression enabled, which should generate much smaller files.

[pazner]

If you're specifically using linear finite elements as in your example, then you can use mesh->PrintVTK(vtkFs) (i.e. without any ref parameter) to output the linear mesh with no duplication of vertices.

The following function will write a linear GridFunction in VTK format defined on such a mesh (not this will not work if you use finite elements with order greater than 1):

void SaveLinearGridFunctionVTK(const GridFunction &gf, std::ostream &out,
                               const std::string &field_name)
{
   const FiniteElementSpace *fes = gf.FESpace();
   Mesh *mesh = fes->GetMesh();
   Vector val;
   DenseMatrix vval, pmat;
   int vec_dim = gf.VectorDim();

   MFEM_VERIFY(gf.Size()/vec_dim == mesh->GetNV(),
               "GridFunction must have as many dofs as mesh vertices");

   out << "POINT_DATA " << fes->GetNDofs() << '\n' << flush;

   if (vec_dim == 1)
   {
      // scalar data
      out << "SCALARS " << field_name << " double 1\n"
          << "LOOKUP_TABLE default\n";
      for (int i=0; i < fes->GetNDofs(); ++i)
      {
         out << gf[i] << '\n';
      }
   }
   else if ((vec_dim == 2 || vec_dim == 3) && mesh->SpaceDimension() > 1)
   {
      // vector data
      out << "VECTORS " << field_name << " double\n";
      Array<int> vdofs(vec_dim);
      for (int i=0; i < fes->GetNDofs(); ++i)
      {
         vdofs.SetSize(1);
         vdofs[0] = i;
         fes->DofsToVDofs(vdofs);
         out << gf[vdofs[0]] << ' ' << gf[vdofs[1]] << ' ';
         if (vec_dim == 2)
         {
            out << 0.0;
         }
         else
         {
            out << gf[vdofs[2]];
         }
         out << '\n';
      }
   }
   else
   {
      // other data: save the components as separate scalars
      for (int vd = 0; vd < vec_dim; vd++)
      {
         out << "SCALARS " << field_name << vd << " double 1\n"
             << "LOOKUP_TABLE default\n";
         Array<int> vdofs(vec_dim);
         for (int i=0; i < fes->GetNDofs(); ++i)
         {
            vdofs.SetSize(1);
            vdofs[0] = i;
            fes->DofsToVDofs(vdofs);
            out << gf[vdofs[vd]] << '\n';
         }
      }
   }
   out.flush();
}

[NikolayYavich]

Hi @pazner

Thanks a lot for the detailed answers, I'll check these options out.