Skip to content

Instantly share code, notes, and snippets.

@jwpeterson
Created January 11, 2021 16:36
Show Gist options
  • Save jwpeterson/f1a4518b748c45f04d85e1372335ebec to your computer and use it in GitHub Desktop.
Save jwpeterson/f1a4518b748c45f04d85e1372335ebec to your computer and use it in GitHub Desktop.
// 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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment