INTREPID2 Documentation (Development Version)

Introduction

Intrepid2 is an extension of Intrepid, a library of interoperable tools for compatible discretizations of Partial Differential Equations (PDEs). Intrepid2 utilizes Kokkos dynamic rank views as the default multidimensional array type, which enables the use of Intrepid2 on heterogeneous architectures.

Overview

Current release of Intrepid2 includes the following features:

• Default finite element basis functions for H(grad), H(curl), H(div) and H(vol) spaces of orders up to 2 on standard cell topologies in 1D, 2D and 3D
• High-order (up to 10) basis functions for H(grad), H(curl), H(div) and H(vol) spaces on select cell topologies
• Pullbacks (transformations) from reference coordinate frame of H(grad), H(curl), H(div) and H(vol) fields
• Pullbacks of gradient, curl and divergence of H(grad), H(curl), H(div) fields
• Cubature rules of orders up to 20 on most standard 1D, 2D and 3D cell topologies

Quick Start

Familiarity with with the following concepts, objects, and tools is required:

• Shards cell topologies,
• numerical integration / Intrepid2::Cubature,
• discrete (e.g. finite element) bases / Intrepid2::Basis / Intrepid2 basis class,
• cell mappings and transformations / Intrepid2::CellTools / Cell tools, and
• function mappings (pullbacks) / Intrepid2::FunctionSpaceTools / Function space tools.

The following example demonstrates, in 7 steps, the computation of finite element stiffness matrices on a set of tetrahedral cells using a piecewise linear basis and an appropriate integration rule.

Step 1: Select a cell topology

shards::CellTopology cellType = shards::getCellTopologyData< shards::Tetrahedron<> >(); // cell type: tetrahedron
ordinal_type spaceDim = cellType->getDimension();                                       // retrieve spatial dimension
ordinal_type numNodes = cellType->getNodeCount();                                       // retrieve number of 0-cells (nodes)

We additionally set the number of computational cells numCells.

Step 2: Select integration (cubature) rule

DefaultCubatureFactory<double> cubFactory;                                              // create cubature factory
ordinal_type cubDegree = 2;                                                             // set cubature degree, e.g. 2
Teuchos::RCP<Cubature<double> > myCub = cubFactory.create(cellType, cubDegree);         // create default cubature
ordinal_type numCubPoints = myCub->getNumPoints();                                      // retrieve number of cubature points

Step 3: Select discrete basis

Basis_HGRAD_TET_C1_FEM<double, DynRankView<double> > tetBasis;                          // create tet basis
ordinal_type numFields = tetBasis.getCardinality();                                     // get basis cardinality

Step 4: Format multi-dimensional arrays

DynRankView<double> cub_points(numCubPoints, spaceDim);
DynRankView<double> cub_weights(numCubPoints);

DynRankView<double> cell_nodes(numCells, numNodes, spaceDim);

DynRankView<double> jacobian(numCells, numCubPoints, spaceDim, spaceDim);
DynRankView<double> jacobian_inv(numCells, numCubPoints, spaceDim, spaceDim);
DynRankView<double> jacobian_det(numCells, numCubPoints);
DynRankView<double> weighted_measure(numCells, numCubPoints);

DynRankView<double> grad_at_cub_points(numFields, numCubPoints, spaceDim);
DynRankView<double> transformed_grad_at_cub_points(numCells, numFields, numCubPoints, spaceDim);
DynRankView<double> weighted_transformed_grad_at_cub_points(numCells, numFields, numCubPoints, spaceDim);
DynRankView<double> stiffness_matrices(numCells, numFields, numFields);

We assume that the array cell_nodes is filled with nodes defining a set of computational (physical) cells.

Step 5: Evaluate differential operator applied to basis at cubature points

myCub->getCubature(cub_points, cub_weights);                                          // retrieve cubature points and weights
tetBasis.getValues(grad_at_cub_points, cub_points, OPERATOR_GRAD);                    // evaluate grad operator at cubature points

Step 6: Apply cell tools

CellTools<double>::setJacobian(jacobian, cub_points, cell_nodes, cellType);           // compute cell Jacobians
CellTools<double>::setJacobianInv(jacobian_inv, jacobian);                            // compute inverses of cell Jacobians
CellTools<double>::setJacobianDet(jacobian_det, jacobian);                            // compute determinants of cell Jacobians

Step 7: Apply function space tools

FunctionSpaceTools::computeCellMeasure<double>(weighted_measure,                      // compute weighted cell measure
                                               jacobian_det,
                                               cub_weights);
FunctionSpaceTools::HGRADtransformGRAD<double>(transformed_grad_at_cub_points,        // transform reference gradients into physical space
                                               jacobian_inv,
                                               grad_at_cub_points);
FunctionSpaceTools::multiplyMeasure<double>(weighted_transformed_grad_at_cub_points,  // multiply with weighted measure
                                            weighted_measure,
                                            transformed_grad_at_cub_points);
FunctionSpaceTools::integrate<double>(stiffness_matrices,                             // compute stiffness matrices
                                      transformed_grad_at_cub_points,
                                      weighted_transformed_grad_at_cub_points,
                                      COMP_CPP);

The computed (local) stiffness matrices can now be used in the assembly of a (global) discrete differential operator, e.g. a discrete Laplacian.

Done!