Handling 3D surface quad meshes

[patrickb314]

I have a miniapp I'm working on that outputs a 2D (i/j-indexed) surface on which a mesh resides and variables reside on that Mesh. I've struggled to get this working under Silo/Visit, though I don't know if this is due to my own mistakes, Silo limitations, Visit limitations, or some combination of the three.

1. First I tried to use a 3D curvilinear mesh with one node in the Z dimension. Silo seemed to write things fine (though I've had trouble knowing how to verify this from the command line), but I couldn't get Visit to read the resulting curvilinear mesh. I don't know if this is because such a mesh has no (3D) zones and hence isn't legal, or if I was just making mistakes with the Silo calls that wrote it.
2. I've wirten the data successfully as a 2D regular mesh with a 3D variable at each point describing its position, but haven't been able to get this to render as a 2D surface with these points in 3-space.
3. I haven't yet tried stating this as a UCD mesh, though that seems like the next potential direction.

What's the right way to express what I want in Silo and render it in Visit and/or Paraview, or should I be looking at a different I/O library and/or visualization system to do what I need to do?

[markcmiller86]

Are you willing to attach your example Silo files here? If so, I could give very clear guidance on what might be missing from them.

For example, for a 3D curvilinear mesh (e.g. DBPutQuadmesh()) which is really a 2D surface in 3 space, you probably need to set DBOPT_TOPO_DIM DBOptlist option to 2 so that VisIt knows its really a surface in 3 space and not a volumetric mesh.

Your option 2, if I understand it, is a 2D mesh (e.g. DBPutQuadmesh()) but with some variables (e.g. DBPutQuadvar()) that have 3 components...maybe those variables represent the x, y and z cartesian coordinates...but VisIt will never know that because the coordinates of the mesh are those writte in the Silo call that creates the mesh in the file. They are not a wholly new variable on the mesh.

You can do it with a UCD mesh too but that is less interesting and useful if you trully have a structured gridded situation.

Again, having some example files would be great if you are willing to attach 'em here 😉

[markcmiller86]

In looking at the Silo source code, I am getting the funny feeling Silo's DBPutQuadmesh() call doesn't allow us to construct the object we want. I know we can do it with a DBPutUcdmesh() call but we loose the structured grid. I am still investigating.

[patrickb314]

Here's a try that Visit doesn't like.

BeatnikOutput0000000000.pdb.gz
.

[markcmiller86]

Ok, I have something working (attached below) but not without some caveats (below). I modified VisIt's quad_disk.C example to elevate the disk out of the 2D plane.

• I had to patch VisIt's ghost zone communication routines slightly and don't like what it needed to fix this.
• The 3D surface is indeed a 2D curvilinear mesh of quads but it is effectively an elevated surface...not fully arbitrary 3D surface. I will try a fully arbitrary example next.
• I needed to output it to Silo in an unusual way (which may relate and work only for the limited case I identified above).
• I can see some minor fixes to make it fully general.

Silo Code

    double *coordsD[3];
    coordsD[0] = xcoords;
    coordsD[1] = ycoords;
    coordsD[2] = sphElevD;

    dbfile = DBCreate("quad_disk.silo", DB_CLOBBER, DB_LOCAL,
                      "2D logical grid deformed into a disk", driver);

    char *coordnames[3];
    coordnames[0] = "xcoords";
    coordnames[1] = "ycoords";
    coordnames[2] = "elevations";
    int ndims = 3;
    int dims[3];
    dims[0] = N; 
    dims[1] = N;
    dims[2] = 1;

    // Write 3D elevated surface
    DBPutQuadmesh(dbfile, "mesh3D", coordnames, (float**) coordsD, dims, ndims,
        DB_DOUBLE, DB_NONCOLLINEAR, NULL);

Example File

The mesh variable in the file is mesh3D. There are other 2D meshes in the file also but you can ignore those.

quad_disk.silo.gz

Screen shot from (patched) VisIt

[markcmiller86]

Ok, here is another example based on VisIt's test data generator specmix.C. In this case, instead of z representing the "3rd" dimension, I made y be that dimension. And, when I do that, I begin to see that other operations in VisIt do not behave as nicely as they should. For example, VisIt's IndexSelect operator drops y-oriented edges.

Silo code

  // Create a 3D version of the same thing on the y=x+z plane
  {
      int dims3d[3] = {dims[0], 1, dims[1]};
      float f12[1000];
      float *coord3d[3] = {f1, f12, f2};
      char *coordnames3d[3] = {"f1","f1+f2","f2"};

      c=0;
      for (x=0;x<mesh.nx;x++) {
        for (y=0;y<mesh.ny;y++) {
          f12[c]=f1[c]+f2[c]; // y = x + z
          c++;
        }
      }
      
      DBPutQuadmesh(db, "Mesh_3d", coordnames3d, coord3d, dims3d, 3, 
                    DB_FLOAT, DB_NONCOLLINEAR, NULL);
  }

Data File

You want the Mesh_3d and p3d and u3d variables in the file...

specmix_quad.silo.gz

Screen shot from VisIt for two orientations

[patrickb314]

This is great, thanks. So, basically, this seems to be a VisIt issue, not a Silo issue. In particular, I was able to use Paraview to read the files I was already generating with Silo, which is something I hadn't tried before.

[patrickb314]

Silo can generate this type of data as shown in the examples Mark provides above, but VisIt requires a patch to be able to read them. Paraview can read these files without patching.

[markcmiller86]

FYI...this commit fixes the problem with VisIt.