prism15_embedding.cc

// This utility code is used for generating the embedding matrix for the Prism15 element
#include <iomanip> // std::setw

#include "libmesh/libmesh.h"
#include "libmesh/mesh.h"
#include "libmesh/elem.h"
#include "libmesh/mesh_generation.h"
#include "libmesh/mesh_modification.h"
#include "libmesh/getpot.h"
#include "libmesh/fe_interface.h"
#include "libmesh/fe_type.h"

// Bring in the libmesh namespace
using namespace libMesh;


// Construct the nodes of a prism child with a base defined by the given vertices
void build_child(Point p0, Point p1, Point p2, std::vector<Point> & child_points)
{
  // Clear out any old data
  child_points.clear();
  child_points.resize(15);

  // Add the base vertex points
  child_points[0] = p0;
  child_points[1] = p1;
  child_points[2] = p2;

  // The reference element children are height 1.
  Real child_height = 1.;

  // Add the "roof" vertex points.
  child_points[3] = p0 + Point(0., 0., child_height);
  child_points[4] = p1 + Point(0., 0., child_height);
  child_points[5] = p2 + Point(0., 0., child_height);

  // Add base mid-edge points
  child_points[6] = 0.5*(p0+p1);
  child_points[7] = 0.5*(p1+p2);
  child_points[8] = 0.5*(p2+p0);

  // Add the vertical side mid-edge points
  child_points[ 9] = p0 + Point(0., 0., 0.5*child_height);
  child_points[10] = p1 + Point(0., 0., 0.5*child_height);
  child_points[11] = p2 + Point(0., 0., 0.5*child_height);

  // Add the "roof" mid-edge points
  child_points[12] = child_points[6] + Point(0., 0., child_height);
  child_points[13] = child_points[7] + Point(0., 0., child_height);
  child_points[14] = child_points[8] + Point(0., 0., child_height);
}



int main (int argc, char** argv)
{
  LibMeshInit init(argc, argv);
  {
    // Start by reading a reference element from file
    std::string libmesh_dir = std::getenv("LIBMESH_DIR");
    Mesh mesh(init.comm(), 3);
    // mesh.read(libmesh_dir+std::string("/share/reference_elements/3D/one_prism15.xda"));

    // We can use the nodes of the Prism18 to form the child "bases"
    mesh.read(libmesh_dir+std::string("/share/reference_elements/3D/one_prism18.xda"));

    Elem * prism18 = mesh.elem(0);

    // To be filled up with points
    std::vector<Point> child_points;

    // See the ASCII art in cell_prism18.h for details of this numbering
    const unsigned child_base_points[8][3] =
      {
        // bottom
        {0,6,8},
        {6,1,7},
        {8,7,2},
        {6,7,8},
        // top
        {9,15,17},
        {15,10,16},
        {17,16,11},
        {15,16,17},
      };

    // Print header so we can copy/paste this directly into source code.
    libMesh::out << "const float Prism15::_embedding_matrix[8][15][15] =" << std::endl;
    libMesh::out << '{' << std::endl;

    for (unsigned child=0; child<8; ++child)
      {
        libMesh::out << "\n// Embedding matrix for child " << child << std::endl;
        libMesh::out << "{" << std::endl;

        build_child(prism18->point(child_base_points[child][0]),
                    prism18->point(child_base_points[child][1]),
                    prism18->point(child_base_points[child][2]),
                    child_points);

        // Verification: For each child, construct a single-Prism15
        // element Mesh and verify the elements are valid.
        const bool do_verification = false;

        if (do_verification)
          {
            // Print points determined for the child
            for (unsigned i=0; i<child_points.size(); ++i)
              libMesh::out << child_points[i] << std::endl;

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

            Elem * prism15 = Elem::build(PRISM15).release();

            // Add all the points to the verification Mesh and set Node pointers in the element
            for (unsigned node=0; node<child_points.size(); ++node)
              prism15->set_node(node) = child_mesh.add_point(child_points[node]);

            // Finally add the Element to the Mesh.  We'll let its
            // destructor ensure that everything is cleaned up.
            prism15 = mesh.add_elem(prism15);

            // Compute the element volume -- note: since there is no
            // volume() specialization for Prism15, this actually reinits
            // a finite element and computes JxW to compute the volume.
            // If this produces a non-negative value, it implies the
            // element is not inverted.
            libMesh::out << "Child " << child << " volume: " << prism15->volume() << std::endl;
          }

        // Compute the embedding matrix entries: for each child,
        // evaluate each of the shape functions at each of the child
        // nodes.
        unsigned n_sf = FEInterface::n_shape_functions(/*dim=*/3,
                                                       FEType(SECOND, LAGRANGE),
                                                       PRISM15);

        // Write out a documentation row numbering the columns.
        libMesh::out << "  //";
        for (unsigned node=0; node<child_points.size(); ++node)
          {
            // The first column needs just slightly less space...
            unsigned width = node==0 ? 8 : 9;
            libMesh::out << std::setw(width) << node;
          }
        libMesh::out << std::endl;

        for (unsigned node=0; node<child_points.size(); ++node)
          {
            // Indent each row two spaces
            libMesh::out << "  {";

            for (unsigned sf=0; sf<n_sf; ++sf)
              {
                Real val = FEInterface::shape(/*dim=*/3,
                                              FEType(SECOND, LAGRANGE),
                                              PRISM15,
                                              sf,
                                              child_points[node]);

                // Avoid printing '-0'.  We know there are no really small shape function values.
                val = (std::abs(val) < 1.e-6) ? 0. : val;

                // Note: the default for std::fixed is to pad with zeros to the requested width?
                char col_comma = (sf+1 == n_sf) ? ' ' : ',';
                libMesh::out << std::setw(8) << val << col_comma;
              }

            char row_comma = (node+1 == child_points.size()) ? ' ' : ',';
            libMesh::out << '}' << row_comma << " // " << std::setw(2) << node << std::endl;
          }

        char child_comma = (child+1 == 8) ? ' ' : ',';
        libMesh::out << '}' << child_comma << std::endl;
      }

    // Print footer
    libMesh::out << "};" << std::endl;
  }

  return 0;
}