Polyhedral mesh renders as split, negative volumes, and ordering questions

[mlohry]

The following code for a single tetrahedron mesh generates a blueprint hdf5 mesh file and render,

  std::vector<double> xvert = {0.0, 1.0, 0.0, 0.0};
  std::vector<double> yvert = {0.0, 0.0, 0.0, 1.0};
  std::vector<double> zvert = {1.0, 0.0, 0.0, 0.0};
  std::vector<int>    conn  = {0, 1, 2, 3};
  conduit::Node       mesh;
  mesh["coordsets/coords/type"] = "explicit";
  mesh["coordsets/coords/values/x"].set_external(xvert);
  mesh["coordsets/coords/values/y"].set_external(yvert);
  mesh["coordsets/coords/values/z"].set_external(zvert);
  mesh["topologies/topo/type"]           = "unstructured";
  mesh["topologies/topo/elements/shape"] = "tet";
  mesh["topologies/topo/coordset"]       = "coords";
  mesh["topologies/topo/elements/connectivity"].set_external(conn);

  ascent::Ascent ascent_viewer;
  conduit::Node  ascent_opts;
  ascent_opts["mpi_comm"] = MPI_Comm_c2f(MPI_COMM_WORLD);
  ascent_viewer.open(ascent_opts);
  ascent_viewer.publish(mesh);
  conduit::Node actions;
  conduit::Node& add_actions = actions.append();
  add_actions["action"]      = "add_scenes";
  conduit::Node& scenes      = add_actions["scenes"];
  scenes["s1/plots/p1/type"] = "mesh";
  scenes["s1/image_prefix"]  = "single_tet_mesh_render";
  conduit::Node& add_extracts    = actions.append();
  add_extracts["action"]         = "add_extracts";
  conduit::Node& extracts        = add_extracts["extracts"];
  extracts["e0/type"]            = "relay";
  extracts["e0/params/path"]     = "single_tet_mesh_extract";
  extracts["e0/params/protocol"] = "blueprint/mesh/hdf5";
  ascent_viewer.execute(actions);
  ascent_viewer.close();

The rendered image and the mesh viewed in visit are as expected,

Attempting to render the same thing, but using the general polyhedral format,

  std::vector<double> xvert = {0.0, 1.0, 0.0, 0.0};
  std::vector<double> yvert = {0.0, 0.0, 0.0, 1.0};
  std::vector<double> zvert = {1.0, 0.0, 0.0, 0.0};
  conduit::Node       mesh;
  mesh["coordsets/coords/type"] = "explicit";
  mesh["coordsets/coords/values/x"].set_external(xvert);
  mesh["coordsets/coords/values/y"].set_external(yvert);
  mesh["coordsets/coords/values/z"].set_external(zvert);
  std::vector<std::int32_t> poly_face_conn = {
      0,1,3,  // face 0
      1,2,3,  // face 1
      2,0,3,  // face 2
      0,2,1   // face 3
  };
  std::vector<std::int32_t> poly_face_sizes   = {3, 3, 3, 3};
  std::vector<std::int32_t> poly_face_offsets = {0, 3, 6, 9};
  std::vector<std::int32_t> poly_cell_conn    = {0, 1, 2, 3};
  std::vector<std::int32_t> poly_cell_sizes   = {4};
  std::vector<std::int32_t> poly_cell_offsets = {0};

  mesh["topologies/topo/type"]              = "unstructured";
  mesh["topologies/topo/elements/shape"]    = "polyhedral";
  mesh["topologies/topo/subelements/shape"] = "polygonal";
  mesh["topologies/topo/coordset"]          = "coords";
  mesh["topologies/topo/elements/connectivity"].set_external(poly_cell_conn);
  mesh["topologies/topo/elements/sizes"].set_external(poly_cell_sizes);
  mesh["topologies/topo/elements/offsets"].set_external(poly_cell_offsets);
  mesh["topologies/topo/subelements/connectivity"].set_external(poly_face_conn);
  mesh["topologies/topo/subelements/sizes"].set_external(poly_face_sizes);
  mesh["topologies/topo/subelements/offsets"].set_external(poly_face_offsets);
  std::vector<double> ones(xvert.size(), 1.0);
  mesh["fields/dummy/association"] = "vertex";
  mesh["fields/dummy/topology"] = "topo";
  mesh["fields/dummy/values"].set_external(ones); // needed to work around https://github.com/Alpine-DAV/ascent/issues/1260

  ascent::Ascent ascent_viewer;
  conduit::Node  ascent_opts;
  ascent_opts["mpi_comm"] = MPI_Comm_c2f(MPI_COMM_WORLD);
  ascent_viewer.open(ascent_opts);
  ascent_viewer.publish(mesh);
  conduit::Node actions;
  conduit::Node& add_actions = actions.append();
  add_actions["action"]      = "add_scenes";
  conduit::Node& scenes      = add_actions["scenes"];
  scenes["s1/plots/p1/type"] = "mesh";
  scenes["s1/image_prefix"]  = "single_tet_poly_mesh_render";
  conduit::Node& add_extracts    = actions.append();
  add_extracts["action"]         = "add_extracts";
  conduit::Node& extracts        = add_extracts["extracts"];
  extracts["e0/type"]            = "relay";
  extracts["e0/params/path"]     = "single_tet_poly_mesh_extract";
  extracts["e0/params/protocol"] = "blueprint/mesh/hdf5";

  std::cout << actions.to_yaml() << std::endl;
  ascent_viewer.execute(actions);
  ascent_viewer.close();

In the render it appears the tetrahedron has been subdivided,

and when loading the mesh into visit, what appears to be subdivided cells, which additionally have all negative volumes by visit's calculation.

The polygonal face ordering used here is (I believe) the same as that specified in the VTK source definition for a tetra 4.

// vtk tetra ordering
  std::vector<std::int32_t> poly_face_conn = {
      0,1,3,  // face 0
      1,2,3,  // face 1
      2,0,3,  // face 2
      0,2,1   // face 3
  };

Looking at the VTK-m source is a different winding equivalent to

//vtk-m ordering
std::vector<std::int32_t> poly_face_conn = {
    0,3,1,  // face 0
    1,2,3,  // face 1
    0,2,3,  // face 2
    0,2,1   // face 3
};

again the render is split but this time some sub-cells are listed with negative volumes of -0.01389 and positive +0.01389.

If I use the SIDS tetra_4 ordering for faces equivalent to

// sids/cgns ordering
std::vector<std::int32_t> poly_face_conn = {
    0,2,1,  // face 0
    0,1,3,  // face 1
    1,2,3,  // face 2
    2,0,3   // face 3
};

again the cell is subdivided, this time with all negative -0.01389 volumes showing in visit:

Questions:

1. Is the subdivision of the poly cell expected behavior or a result of bad ordering?
2. If the subdivision is expected, can I avoid this so the exported mesh "looks like" my source mesh?
3. What is the appropriate ordering for polyhedral cell types? Specifically interested in TETRA_4, PYRA_5, PENTA_6, and HEXA_8 cells for now.

[mlohry]

Another example, this time using the two pyramid connectivity example from the blueprint docs

  const std::string   image_filename = "two_pyr_poly_mesh_render";
  const std::string   mesh_filename  = "two_pyr_poly_mesh_extract";
  std::vector<double> xvert          = {0.5, 0.0, 1.0, 0.0, 1.0, 0.5};
  std::vector<double> yvert          = {1.0, 0.0, 0.0, 0.0, 0.0, -1.0};
  std::vector<double> zvert          = {0.5, 1.0, 1.0, 0.0, 0.0, 0.5};
  conduit::Node       mesh;
  mesh["coordsets/coords/type"] = "explicit";
  mesh["coordsets/coords/values/x"].set_external(xvert);
  mesh["coordsets/coords/values/y"].set_external(yvert);
  mesh["coordsets/coords/values/z"].set_external(zvert);

  // connectivity from example https://llnl-conduit.readthedocs.io/en/latest/blueprint_mesh.html
  std::vector<std::int32_t> poly_face_conn = {1, 2, 4, 3, 1, 2, 0, 2, 4, 0, 4, 3, 0, 3, 1, 0, 1, 2, 5, 2, 4, 5, 4, 3, 5, 3, 1, 5};
  std::vector<std::int32_t> poly_face_sizes   = {4, 3, 3, 3, 3, 3, 3, 3, 3};
  std::vector<std::int32_t> poly_face_offsets = {0, 4, 7, 10, 13, 16, 19, 22, 25};
  std::vector<std::int32_t> poly_cell_conn    = {0, 1, 2, 3, 4, 0, 5, 6, 7, 8};
  std::vector<std::int32_t> poly_cell_sizes   = {5, 5};
  std::vector<std::int32_t> poly_cell_offsets = {0, 5};

  mesh["topologies/topo/type"]              = "unstructured";
  mesh["topologies/topo/elements/shape"]    = "polyhedral";
  mesh["topologies/topo/subelements/shape"] = "polygonal";
  mesh["topologies/topo/coordset"]          = "coords";
  mesh["topologies/topo/elements/connectivity"].set_external(poly_cell_conn);
  mesh["topologies/topo/elements/sizes"].set_external(poly_cell_sizes);
  mesh["topologies/topo/elements/offsets"].set_external(poly_cell_offsets);
  mesh["topologies/topo/subelements/connectivity"].set_external(poly_face_conn);
  mesh["topologies/topo/subelements/sizes"].set_external(poly_face_sizes);
  mesh["topologies/topo/subelements/offsets"].set_external(poly_face_offsets);
  // dummy field needed to work around https://github.com/Alpine-DAV/ascent/issues/1260
  std::vector<double> ones(xvert.size(), 1.0);
  mesh["fields/dummy/association"] = "vertex";
  mesh["fields/dummy/topology"]    = "topo";
  mesh["fields/dummy/values"].set_external(ones);

  conduit::Node verify_info;
  const bool    verified = conduit::blueprint::mesh::verify(mesh, verify_info);
  ASSERT_TRUE(verified) << verify_info.to_yaml() << std::endl;

  ascent::Ascent ascent_viewer;
  conduit::Node  ascent_opts;
  ascent_opts["mpi_comm"] = MPI_Comm_c2f(MPI_COMM_WORLD);
  ascent_viewer.open(ascent_opts);
  ascent_viewer.publish(mesh);

  conduit::Node actions;

  // declare an action to create a .png render
  conduit::Node& add_actions = actions.append();
  add_actions["action"]      = "add_scenes";
  conduit::Node& scenes      = add_actions["scenes"];
  scenes["s1/plots/p1/type"] = "mesh";
  scenes["s1/image_prefix"]  = image_filename;

  // declare an action to create a mesh extraction to hdf5
  conduit::Node& add_extracts    = actions.append();
  add_extracts["action"]         = "add_extracts";
  conduit::Node& extracts        = add_extracts["extracts"];
  extracts["e0/type"]            = "relay";
  extracts["e0/params/path"]     = mesh_filename;
  extracts["e0/params/protocol"] = "blueprint/mesh/hdf5";

  std::cout << actions.to_yaml() << std::endl;
  ascent_viewer.execute(actions);
  //  conduit::Node info;
  //  ascent_viewer.info(info);
  //  info.print();
  ascent_viewer.close();

  std::cout << "Image file should have been written to disk: " << image_filename << std::endl;
  std::cout << "Mesh file should have been written to disk: " << mesh_filename << std::endl;

This also shows negative volumes, and oddly the mesh edge lines don't appear (or are white and opaque)

[cyrush]

@mlohry Thank you for all of these examples.

We have some face ordering constraints in another set of codes that I think explain why we have negative volumes in VisIt. We will look at addressing those by making VisIt's volume expressions more robust.

For the internal mesh lines, this happens b/c we do subdivide the mesh to render using simpler primitives.

VisIt has logic to remove the mesh lines of automatically subdivided mesh (when it is presented directly with a polytopal mesh). Ascent doesn't have the same logic, but needs it.

Ascent subdividing is impacting the extract (you are exporting the subdivided result instead of the original published mesh)

Both versions are useful, so we will look at adding options that allow both flavors to be exported.