Evaluate solution at given point in Hex20 element.

andrew_yt.cc

#include "libmesh/libmesh.h"
#include "libmesh/mesh.h"
#include "libmesh/elem.h"
#include "libmesh/cell_hex20.h"
#include "libmesh/fe.h"
#include "libmesh/fe_interface.h"

using namespace libMesh;

// This is an example hopfully showing YT people how to evaluate
// solution values on a single Hex20 element.
int main (int argc, char** argv)
{
  LibMeshInit init(argc, argv);

  // Node positions from the original email:
  const Real vertices[20][3] =
    {
      {3.00608789e-03,   4.64941000e-02,  -3.95758979e-04},
      {3.03202730e-03,   4.64941000e-02,   0.00000000e+00},
      {3.03202730e-03,   4.70402000e-02,   3.34389809e-20},
      {3.00608789e-03,   4.70402000e-02,  -3.95758979e-04},
      {2.45511948e-03,   4.64941000e-02,  -3.23222611e-04},
      {2.47630461e-03,   4.64941000e-02,   1.20370622e-35},
      {2.47630461e-03,   4.70402000e-02,   3.34389809e-20},
      {2.45511948e-03,   4.70402000e-02,  -3.23222611e-04},
      {3.01905760e-03,   4.64941000e-02,  -1.97879489e-04},
      {3.03202730e-03,   4.67671500e-02,   3.34389809e-20},
      {3.01905760e-03,   4.70402000e-02,  -1.97879489e-04},
      {3.00608789e-03,   4.67671500e-02,  -3.95758979e-04},
      {2.73060368e-03,   4.64941000e-02,  -3.59490795e-04},
      {2.75416596e-03,   4.64941000e-02,  -1.86574463e-34},
      {2.75416596e-03,   4.70402000e-02,   6.68779617e-20},
      {2.73060368e-03,   4.70402000e-02,  -3.59490795e-04},
      {2.47100265e-03,   4.64941000e-02,  -1.61958070e-04},
      {2.47630461e-03,   4.67671500e-02,   1.67194904e-20},
      {2.47100265e-03,   4.70402000e-02,  -1.61958070e-04},
      {2.45511948e-03,   4.67671500e-02,  -3.23222611e-04}
    };

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

  // Create a Hex20 element and assign points to it.
  Elem * elem = mesh.add_elem(new Hex20);

  for (unsigned i = 0; i<20; ++i)
    {
      Node * node = mesh.add_point(Point(vertices[i][0],
                                         vertices[i][1],
                                         vertices[i][2]), i);
      elem->set_node(i) = node;
    }

  // Write to file, make sure it looks OK.
  // mesh.prepare_for_use();
  // mesh.write("andrew_yt.e");

  // Field values from the original email:
  const Real field_values[20] = {659.80151432,  650.95995348,  650.02809796,  658.81589888,
                                 659.77560908,  650.93582507,  649.99987015,  658.78508795,
                                 655.3807339 ,  650.49402572,  654.42199842,  659.3087066 ,
                                 659.7885617 ,  650.94788928,  650.01398406,  658.80049342,
                                 655.35571708,  650.46784761,  654.39247905,  659.28034852};

  // Point to evaluate the solution (physical space) from email:
  Point physical_x = Point(0.00274291,
                           0.04676715,
                           -0.00017978);

  // A second-order Lagrange FE type object.
  FEType fe_type = FEType(SECOND, LAGRANGE);

  // Map physical_x into reference coordinates
  Point reference_x = FEInterface::inverse_map(/*dim=*/3,
                                               fe_type,
                                               elem,
                                               physical_x);

  // Debugging: check that the reference point is in [-1,1]^3
  // libMesh::out << "reference_x = " << reference_x << std::endl;

  // Construct an FE object -- there does not seem to be an
  // FEInterface function that just computes the solution at a custom
  // point, so we will use the standard reinit() method.
  UniquePtr<FEBase> fe (FEBase::build(/*dim=*/3, fe_type));

  // Request stuff we need for evaluating the solution field so that
  // reinit() knows to compute it.
  const std::vector<std::vector<Real> >& phi = fe->get_phi();

  // Create custom vector of reference points where shape functions
  // should be evaluated.  (The reinit() interface requires this.)
  std::vector<Point> ref_points(1, reference_x);
  
  // Recomputes shape function values, etc. at the requested point.
  fe->reinit(elem, &ref_points);

  // Evaluate solution value at physical_x.  u = sum_i u_i * phi_i
  Real field_val = 0.;
  for (unsigned i=0; i<phi.size(); ++i)
    {
      // libMesh::out << phi[i][0] << std::endl;
      field_val += field_values[i]*phi[i][0];
    }
  
  libMesh::out << "field_val=" << field_val << std::endl;
  
  return 0;
}