tet_test_2831.cc

// libmesh includes
#include "libmesh/libmesh.h"
#include "libmesh/mesh.h"
#include "libmesh/elem.h"
#include "libmesh/fe.h"
#include "libmesh/quadrature_gauss.h"
#include "libmesh/dense_vector.h"
#include "libmesh/dense_matrix.h"

using namespace libMesh;

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

  // Even if libmesh is configured with performance logging enabled,
  // let's turn it off for this simple example.
  perflog.disable_logging();

  // Create a mesh of the dimension specified on the command line
  Mesh mesh(init.comm());

  // Location of the reference_elements meshes:
  std::string libmesh_dir = std::getenv("LIBMESH_ROOT");

  // Test that we can correctly compute the volume of the reference element.
  mesh.read(libmesh_dir + std::string("/reference_elements/3D/one_tet10.xda"));

  const unsigned int dim = 3;

  // Create quadrature rule, set flag to only allow rules with non-negative weights.
  QGauss qrule(dim, /*order=*/THIRD);
  qrule.allow_rules_with_negative_weights = false;

  // Set up FE
  std::unique_ptr<FEBase> fe =
    FEBase::build(dim, FEType(SECOND, LAGRANGE));
  fe->attach_quadrature_rule(&qrule);

  // Pre-request FE data
  const auto & JxW = fe->get_JxW();
  const auto & phi = fe->get_phi();

  // A pointer to the only element.
  Elem * elem = mesh.elem_ptr(0);

  // Reinit FE on the elem
  fe->reinit(elem);

  // There should be 8 qps for 2x2x2 conical product rule which is
  // substituted when negative weights are not allowed.
  const unsigned int n_qps = qrule.n_points();
  libMesh::out << "n_qps = " << n_qps << std::endl;

  // Compute mass matrix
  const unsigned int n_shapes = phi.size();

  DenseMatrix<Real> mass_matrix(n_shapes, n_shapes);

  for (unsigned int qp=0; qp<n_qps; qp++)
    for (unsigned int i=0; i<n_shapes; i++)
      for (unsigned int j=0; j<n_shapes; j++)
        mass_matrix(i,j) += JxW[qp] * phi[i][qp] * phi[j][qp];

  libMesh::out << "mass_matrix = " << mass_matrix << std::endl;

  // Compute eigenvalues of mass matrix. Eigenvalues of mass matrix
  // should be strictly real and positive...
  DenseMatrix<Real> A = mass_matrix;
  DenseVector<Real> lambda_real, lambda_imag;
  A.evd(lambda_real, lambda_imag);

  for (auto i : index_range(lambda_real))
    libMesh::out << "Eigenvalue " << i << ", (" << lambda_real(i) << ", " << lambda_imag(i) << ")" << std::endl;

  return 0;
}

// Output
// mass_matrix = 0.00158519 0.000185185 6.66667e-05 -0.00017037 -0.00133333 -0.0022963 -0.00115556 -0.000681481 -0.0022963 -0.00223704 
// 0.000185185 0.00111111 0.000185185 0.000185185 -0.000740741 -0.000740741 -0.00259259 -0.00259259 -0.000740741 -0.00259259 
// 6.66667e-05 0.000185185 0.00134815 6.66667e-05 -0.0022963 -0.00133333 -0.000681481 -0.00271111 -0.0022963 -0.000681481 
// -0.00017037 0.000185185 6.66667e-05 0.00158519 -0.0022963 -0.0022963 -0.00223704 -0.000681481 -0.00133333 -0.00115556 
// -0.00133333 -0.000740741 -0.0022963 -0.0022963 0.0118519 0.00592593 0.00651852 0.00651852 0.00592593 0.00325926 
// -0.0022963 -0.000740741 -0.00133333 -0.0022963 0.00592593 0.0118519 0.00651852 0.00325926 0.00592593 0.00651852 
// -0.00115556 -0.00259259 -0.000681481 -0.00223704 0.00651852 0.00651852 0.0113778 0.00663704 0.00325926 0.00568889 
// -0.000681481 -0.00259259 -0.00271111 -0.000681481 0.00651852 0.00325926 0.00663704 0.0104296 0.00651852 0.00663704 
// -0.0022963 -0.000740741 -0.0022963 -0.00133333 0.00592593 0.00592593 0.00325926 0.00651852 0.0118519 0.00651852 
// -0.00223704 -0.00259259 -0.000681481 -0.00115556 0.00325926 0.00651852 0.00568889 0.00663704 0.00651852 0.0113778 
// 
// Eigenvalue 0, (0.0416667, 0)
// Eigenvalue 1, (0.00934363, 0)
// Eigenvalue 2, (0.00881523, 0)
// Eigenvalue 3, (0.00833333, 0)
// Eigenvalue 4, (0.00146994, 0)
// Eigenvalue 5, (0.00218477, 0)
// Eigenvalue 6, (0.0025568, 0)
// Eigenvalue 7, (5.90387e-19, 0)
// Eigenvalue 8, (-1.96224e-18, 0)
// Eigenvalue 9, (2.36944e-19, 0)