xujun_zhao_test.cc

#include "libmesh/libmesh.h"
#include "libmesh/mesh.h"
#include "libmesh/mesh_generation.h"
#include "libmesh/elem.h"
#include "libmesh/fe.h"
#include "libmesh/quadrature_gauss.h"
#include "libmesh/getpot.h"
#include "libmesh/mesh_modification.h"
#include "libmesh/cell_hex27.h"

#include <iomanip>

using namespace libMesh;

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

  // Supposedly this point and the following element produces a zero
  // determinant?  I am highly skeptical of this.

  // pt[186] = (-11.545647,22.442960,12.336682)
  // Elem id = 449, and volume is 383.684010
  // Element Nodes are :

  Mesh mesh(init.comm());
  unsigned cnt = 0;
  mesh.add_point(Point(-18.750000,7.911975,7.911975), cnt++);
  mesh.add_point(Point(-18.750000,3.750062,7.806584), cnt++);
  mesh.add_point(Point(-18.750000,3.834322,12.613942), cnt++);
  mesh.add_point(Point(-18.750000,8.099410,12.760054), cnt++);
  mesh.add_point(Point(-30.163472,12.728141,12.728141), cnt++);
  mesh.add_point(Point(-31.084410,6.216985,12.942030), cnt++);
  mesh.add_point(Point(-29.400453,6.012310,19.778966), cnt++);
  mesh.add_point(Point(-28.660566,12.380462,19.504553), cnt++);
  mesh.add_point(Point(-18.750000,5.831019,7.859280), cnt++);
  mesh.add_point(Point(-18.750000,3.792192,10.210263), cnt++);
  mesh.add_point(Point(-18.750000,5.966866,12.686998), cnt++);
  mesh.add_point(Point(-18.750000,8.005692,10.336015), cnt++);
  mesh.add_point(Point(-24.456736,10.320058,10.320058), cnt++);
  mesh.add_point(Point(-24.917205,4.983523,10.374307), cnt++);
  mesh.add_point(Point(-24.075226,4.923316,16.196454), cnt++);
  mesh.add_point(Point(-23.705283,10.239936,16.132304), cnt++);
  mesh.add_point(Point(-30.623941,9.472563,12.835085), cnt++);
  mesh.add_point(Point(-30.242432,6.114647,16.360498), cnt++);
  mesh.add_point(Point(-29.030509,9.196386,19.641760), cnt++);
  mesh.add_point(Point(-29.412019,12.554301,16.116347), cnt++);
  mesh.add_point(Point(-18.750000,5.898942,10.273139), cnt++);
  mesh.add_point(Point(-24.686971,7.651791,10.347183), cnt++);
  mesh.add_point(Point(-24.496216,4.953420,13.285381), cnt++);
  mesh.add_point(Point(-23.890255,7.581626,16.164379), cnt++);
  mesh.add_point(Point(-24.081009,10.279997,13.226181), cnt++);
  mesh.add_point(Point(-29.827225,9.334474,16.238423), cnt++);
  mesh.add_point(Point(-24.288613,7.616708,13.255781), cnt++);

  Elem * hex27 = mesh.add_elem(new Hex27);
  for (unsigned i=0; i<27; ++i)
    hex27->set_node(i) = mesh.node_ptr(i);

  mesh.prepare_for_use();
  mesh.write("xujun_zhao_test.e");

  // This should produce a zero determinant, according to Xujun Zhao.
  bool contains_point = hex27->contains_point(Point(-11.545647,22.442960,12.336682));
  libMesh::out << "contains_point=" << contains_point << std::endl;

  return 0;
}