david_tet10.cc

// This test code is to investigate David Knezevic's bad Tet10
#include "libmesh/libmesh.h"
#include "libmesh/mesh.h"
#include "libmesh/elem.h"
#include "libmesh/cell_tet10.h"

// Bring in the libmesh namespace
using namespace libMesh;

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

    {
      Mesh mesh(init.comm(), 3);

      Node* nodes[10];
      unsigned ctr=0;
      nodes[ctr] = mesh.add_point(Point(-0.0677617, 0.0669608, 0.404735), ctr); ctr++;
      nodes[ctr] = mesh.add_point(Point(-1.74574e-15, -6.7675e-32,      0.5), ctr); ctr++;
      nodes[ctr] = mesh.add_point(Point(-0.0555556, 1.93971e-16, 0.444444), ctr); ctr++;
      nodes[ctr] = mesh.add_point(Point(-2.42464e-15, 0.0555556, 0.444444), ctr); ctr++;
      nodes[ctr] = mesh.add_point(Point(-0.0338809, 0.0334804, 0.452368), ctr); ctr++;
      nodes[ctr] = mesh.add_point(Point(-0.0277778, 9.69856e-17, 0.472222), ctr); ctr++;
      nodes[ctr] = mesh.add_point(Point(-0.0639647, 0.0347326, 0.427214), ctr); ctr++;
      nodes[ctr] = mesh.add_point(Point(-0.0351893, 0.0636239, 0.427293), ctr); ctr++;
      nodes[ctr] = mesh.add_point(Point(-2.52163e-15, 0.0277778, 0.472222), ctr); ctr++;
      nodes[ctr] = mesh.add_point(Point(-0.0335307, 0.0335307,  0.45258), ctr); ctr++;

      Elem* elem = mesh.add_elem(new Tet10);
      for (unsigned i=0; i<10; ++i)
        elem->set_node(i) = nodes[i];

      mesh.write("david_tet10.e");
    }

    // The reference Tet10
    {
      Mesh ref_mesh(init.comm(), 3);
      std::string libmesh_dir = std::getenv("LIBMESH_ROOT");
      ref_mesh.read(libmesh_dir+std::string("/reference_elements/3D/one_tet10.xda"));
      ref_mesh.write("one_tet10.e");
    }
    
  }

  return 0;
}