Overview

LehrFEM++ provides an interface for geometry data: lf::geometry::Geometry. Geometry data for entities is stored as a mapping \( \Phi_K(\hat{x}) \) from the reference element (unit triangle, square, unit segment) to any general element. \( \Phi_K^{*}(\hat{x}) \) represents the inverse mapping from the general element to the reference element.

For any one element, it is enough to store the mapping \( \Phi_K(\hat{x}) \) from the reference element and the type of the reference element.

(Affine) Geometry Mapping for Triangles

More details and mathematical background can be found in Lecture Document Section 2.8.

Geometry Interface

To get the geometry of an entity:

for (const lf::mesh::Entity* entity : mesh_p->Entities(codim)) {
  // Get the geometry of an entity
  const lf::geometry::Geometry* geometry = entity->Geometry();
}

Utility Functions

A number of convenience functions are provided by the Geometry class.

const lf::geometry::Geometry* geometry = entity->Geometry();
 
// Get the volume of an entity (for edges length, for triangles area, etc.)
double v = lf::geometry::Volume(*geometry);
 
// Get the corners of an entity as a dxn matrix
// where d is the dimension of physical space, and n the number of corners.
// (for edges 2 corners, for triangles 3 corners, etc.)
Eigen::MatrixXd corners = lf::geometry::Corners(*geometry);

In addition, we can test for some properties of the geometry:

// Test if geometry is a non-degenerate bilinear quadrilateral (corners
// matrix with 4 cols)
bool is_not_deg_q = lf::geometry::assertNonDegenerateQuad(corners);
 
// Test if geometry is a non-degenerate triangle (corners matrix with 3
// cols)
bool is_not_deg_t = lf::geometry::assertNonDegenerateTriangle(corners);

Geometry Methods

Objects implementing the lf::geometry::Geometry interface provide the following methods:

// Dimension of the local coordinate system
unsigned dim_local = geometry->DimLocal();
 
// Dimension of the physical coordinate system
unsigned dim_global = geometry->DimGlobal();
 
// Get the reference element for the geometry
lf::base::RefEl ref_el = geometry->RefEl();
 
// Check if mapping is affine
bool is_affine = geometry->isAffine();

Transformations and Mappings

This part is discussed in detail in Lecture Document Section 2.8. Usage of the following methods is discussed in the Assembly Quick Reference.

Note: The following methods accept a matrix of points in the local coordinate system as input. The points are stored in a matrix where each column represents a point. The number of rows corresponds to the dimension of the local coordinate system.

Global

The lf::geometry::Geometry class provides the lf::geometry::Geometry::Global method to map points from local to global coordinates. In LehrFEM++, local points generally refer to points on the reference element.

// Define two point on the reference element (local).
Eigen::MatrixXd local_points(2, 2);
local_points << 0.1, 0.5, 0.6, 0.2;
 
// Map the points from the local into the global space.
Eigen::MatrixXd global = geometry->Global(local_points);

Mapping of points from local to global coordinates

IntegrationElement

The method lf::geometry::Geometry::IntegrationElement computes the integration element

\[ g(\xi) := \sqrt{\mathrm{det}\left|D\Phi^T(\xi) D\Phi(\xi) \right|} \]

for each point in points as an Eigen::VectorXd.

// Compute the integration element for each point in points
Eigen::VectorXd integration_element =
    geometry->IntegrationElement(local_points);

Note: The following two methods return a matrix for each evaluated point. For efficiency reasons, the returned matrices are horizontally stacked. The number of rows corresponds to the dimension of the global coordinate system.

Jacobian

The method lf::geometry::Geometry::Jacobian computes the Jacobian matrix \( D\Phi(\hat{x}) \) for each point in points.

// Compute the Jacobian matrix for each point in points

Eigen::MatrixXd jacobian = geometry->Jacobian(local_points);

The Jacobian evaluated at every point is itself a matrix of size dim_global x num_points. The JacobianInverseGramian evaluated at every point is a matrix of size dim_local x (dim_global * num_points). To access the Jacobian at a specific point, use the following code:

unsigned dim_global = geometry->DimGlobal();
unsigned dim_local = geometry->DimLocal();
 
// Access the Jacobian matrix for the n-th point (starting at 0)
Eigen::MatrixXd jacobian_n =
    jacobian.block(0, n * dim_local, dim_global, dim_local);

If dim_local == dim_global == 2 and we pass three points for evaluation, Geometry::Jacobian returns a \( 2 \times 6 \) matrix:

Block matrix returned by Geometry::Jacobian

To access the Jacobian for the second point, we can use Eigen block access:

// Access 2x2 Jacobian matrix starting at row 0 and column 2

Eigen::MatrixXd jacobian_2 = jacobian.block(0, 2, 2, 2);

JacobianInverseGramian

The method lf::geometry::Geometry::JacobianInverseGramian computes the inverse of the JacobianInverseGramian matrix for each point in points. Details can be found in Lecture Document Paragraph 2.8.3.14 and the method docs. Similarly to Geometry::Jacobian, it also returns a horizontally stacked matrix. Individual matrices can be accessed using Eigen block.

// Compute the inverse of the Jacobian matrix for each point in points
Eigen::MatrixXd jacobian_inv =
    geometry->JacobianInverseGramian(local_points);
 
// Access the JacobianInverseGramian matrix for the n-th point (starting
// at 0)
Eigen::MatrixXd jacobian_inv_n =
    jacobian_inv.block(0, n * dim_local, dim_global, dim_local);

Previous	Next
Mesh	DOFHandlers

Table of Contents