Issue loading quad4 nodal data from ExodusII

[drwells]

summary

Two nodes in a quadrilateral loaded from ExodusII have the wrong id()s, which causes the nodal data vectors to be in the wrong order.

description of problem

I checked this with both libMesh 1.7.8 and the development version (as of this morning).

I'm loading nodal ExodusII data and, at some point in the process, two nodes of a quadrilateral get each-other's node numbers. More precisely, I start with (via ncdump)

 coordx = 0, 1, 0, 1 ;

 coordy = 0, 0, 1, 1 ;

 eb_names =
  "" ;

 coor_names =
  "",
  "" ;

 connect1 =
  1, 2, 4, 3 ;

but later on libMesh appears to renumber things. The solitary element is read as

  Node id()=0, processor_id()=0, Point=(x,y,z)=(       0,        0,        0)
    DoFs=(0/0/0) (1/0/0) (2/0/0) 
  Node id()=1, processor_id()=0, Point=(x,y,z)=(       1,        0,        0)
    DoFs=(0/0/1) (1/0/1) (2/0/1) 
  Node id()=2, processor_id()=0, Point=(x,y,z)=(       1,        1,        0)
    DoFs=(0/0/2) (1/0/2) (2/0/2) 
  Node id()=3, processor_id()=0, Point=(x,y,z)=(       0,        1,        0)
    DoFs=(0/0/3) (1/0/3) (2/0/3) 

i.e., the nodes 2 and 3 have the wrong id()s, which causes ExodusII_IO::copy_nodal_solution() to read values in the wrong order.

I have absolutely no idea what causes this. Is this some sort of obscure issue related to not supplying the node number map, or something else in the ExodusII file?

exodus file

I narrowed this problem down to a one-element mesh (I started in 3d with hexahedra, the problem is the same there). Here's the test file (via ncdump)

netcdf test {
dimensions:
        len_name = 256 ;
        time_step = UNLIMITED ; // (1 currently)
        num_dim = 2 ;
        num_nodes = 4 ;
        num_elem = 1 ;
        num_el_blk = 1 ;
        num_el_in_blk1 = 1 ;
        num_nod_per_el1 = 4 ;
        num_nod_var = 4 ;
variables:
        double time_whole(time_step) ;
        int eb_status(num_el_blk) ;
        int eb_prop1(num_el_blk) ;
                eb_prop1:name = "ID" ;
        double coordx(num_nodes) ;
        double coordy(num_nodes) ;
        char eb_names(num_el_blk, len_name) ;
                eb_names:_FillValue = "" ;
        char coor_names(num_dim, len_name) ;
                coor_names:_FillValue = "" ;
        int connect1(num_el_in_blk1, num_nod_per_el1) ;
                connect1:elem_type = "QUAD4" ;
        char name_nod_var(num_nod_var, len_name) ;
                name_nod_var:_FillValue = "" ;
        double vals_nod_var1(time_step, num_nodes) ;
        double vals_nod_var2(time_step, num_nodes) ;
        double vals_nod_var3(time_step, num_nodes) ;
        double vals_nod_var4(time_step, num_nodes) ;

// global attributes:
                :api_version = 9.04f ;
                :version = 9.04f ;
                :floating_point_word_size = 8 ;
                :file_size = 1 ;
                :maximum_name_length = 32 ;
                :int64_status = 0 ;
                :title = "test" ;
data:

 time_whole = 0 ;

 eb_status = 1 ;

 eb_prop1 = 1 ;

 coordx = 0, 1, 0, 1 ;

 coordy = 0, 0, 1, 1 ;

 eb_names =
  "" ;

 coor_names =
  "",
  "" ;

 connect1 =
  1, 2, 4, 3 ;

 name_nod_var =
  "X1",
  "X2",
  "exp",
  "monomial" ;

 vals_nod_var1 =
  0, 1, 0, 1 ;

 vals_nod_var2 =
  0, 0, 1, 1 ;

 vals_nod_var3 =
  1, 2.71828182845905, 2.71828182845905, 7.38905609893065 ;

 vals_nod_var4 =
  0, 0, 0, 1 ;
}

source code

There's probably a better way to do this, but since ExodusII_IO::copy_to_nodal_solution() updates system.solution I created one System per vector:

#include <libmesh/equation_systems.h>
#include <libmesh/exodusII_io.h>
#include <libmesh/elem.h>
#include <libmesh/node.h>
#include <libmesh/numeric_vector.h>
#include <libmesh/replicated_mesh.h>

int main(int argc, char **argv) {
  using namespace libMesh;
  LibMeshInit init(argc, argv);

  ReplicatedMesh mesh(init.comm());
  ExodusII_IO exodus_io(mesh);
  exodus_io.read("test.e");
  mesh.prepare_for_use();

  EquationSystems equation_systems(mesh);
  System &X1_system = equation_systems.add_system<System>("X1");
  X1_system.add_variable("X1", FEType(1));
  System &X2_system = equation_systems.add_system<System>("X2");
  X2_system.add_variable("X2", FEType(1));
  equation_systems.init();

  exodus_io.copy_nodal_solution(X1_system, "X1", "X1");
  exodus_io.copy_nodal_solution(X2_system, "X2", "X2");

  const auto &X1 = *X1_system.solution;
  const auto &X2 = *X2_system.solution;

  const Elem *elem = mesh.elem_ptr(0);
  std::cout << "volume = " << elem->volume() << std::endl;
  for (unsigned int node_n = 0; node_n < elem->n_nodes(); ++node_n)
  {
      elem->node_ptr(node_n)->print_info(std::cout);
  }

  for (const auto *node : mesh.node_ptr_range()) {
      std::cout << "id = " << node->id() << std::endl;
    std::cout << "node  = " << (*node)(0) << ", " << (*node)(1) << std::endl;
    std::cout << "field = " << X1(node->dof_number(X1_system.number(), 0, 0))
              << ", " << X2(node->dof_number(X2_system.number(), 0, 0))
              << std::endl;
  }
}

executable output

volume = 1
  Node id()=0, processor_id()=0, Point=(x,y,z)=(       0,        0,        0)
    DoFs=(0/0/0) (1/0/0) 
  Node id()=1, processor_id()=0, Point=(x,y,z)=(       1,        0,        0)
    DoFs=(0/0/1) (1/0/1) 
  Node id()=2, processor_id()=0, Point=(x,y,z)=(       1,        1,        0)
    DoFs=(0/0/2) (1/0/2) 
  Node id()=3, processor_id()=0, Point=(x,y,z)=(       0,        1,        0)
    DoFs=(0/0/3) (1/0/3) 
id = 0
node  = 0, 0
field = 0, 0
id = 1
node  = 1, 0
field = 1, 0
id = 2
node  = 1, 1
field = 0, 1
id = 3
node  = 0, 1
field = 1, 1

[jwpeterson]

It looks to me like the mesh has been automatically renumbered, which libmesh does during prepare_for_use() unless you explicitly tell it not to. Did you try calling

mesh.allow_renumbering(false);

before reading in the exo file? I don't know exactly what version of libmesh we first settled on that API, but it should definitely work in the latest master.

[drwells]

Yup, that's exactly it! Thank you!

[drwells]

That was easier than I thought it would be - I didn't think libMesh would renumber nodes.

[jwpeterson]

The "allowed to renumber by default" behavior has unfortunately been around forever, and in hindsight it is obviously surprising and non-intuitive for most users when it happens. Our attempt to work around it in a backwards-compatible way is that MeshBase API, but we should probably also add a comment to the docs for ExodusII_IO::copy_nodal_solution() (and other APIs) that says something like "this method requires that the mesh has not been renumbered" and suggests the use of MeshBase::allow_renumbering(false);.

[drwells]

Those are good ideas - I'll add some more docs to #4341.

[roystgnr]

I disagree with "unfortunately" - renumbering by default was one of our last bastions of defense of "will handle transient adaptivity efficiently" applied mathematicians' libMesh, after we hit it big with engineers who unwittingly pushed it in the direction of "is this a nodal variable or an elemental variable? those are the two possibilities, right?" libMesh.

But I really love #4341. When we do conflict with the engineer types' code (both outside the library and inside) that definitely needs to result in friendly error messages, not corrupted data.