Intrepid is a library of interoperable tools for compatible discretizations of Partial Differential Equations (PDEs). Included with the Trilinos 10.0 release is the "expert version" of Intrepid. This version is intended primarily for application developers who want to reuse large parts of their existing code frameworks such as I/O, data structures, assembly routines, etc. while gaining access to advanced discretization capabilities provided by Intrepid. In such cases the bulk of the data is owned and managed by the user rather than by Intrepid. To avoid unnecessary and possibly expensive copying of data to and from Intrepid, the expert version of the package comprises of mostly stateless classes operating on user-owned data. Virtually all numerical data required by PDE codes can be represented as a multi-dimensional array of scalar values. For this reason, and to enhance interoprability, Intrepid classes are templated on generic multi-dimensional arrays. The Shards package provides an implementation of a multi-dimensional array that can be used for that purpose, or users can write their own multi-dimensional arrays as long as a minimal interface is supported.
Current release of Intrepid includes the following features:
Familiarity with with the following concepts, objects, and tools is required:
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.
shards::CellTopology cellType = shards::getCellTopologyData< shards::Tetrahedron<> >(); // cell type: tetrahedron int spaceDim = cellType->getDimension(); // retrieve spatial dimension int numNodes = cellType->getNodeCount(); // retrieve number of 0-cells (nodes)
We additionally set the number of computational cells
DefaultCubatureFactory<double> cubFactory; // create cubature factory int cubDegree = 2; // set cubature degree, e.g. 2 Teuchos::RCP<Cubature<double> > myCub = cubFactory.create(cellType, cubDegree); // create default cubature int numCubPoints = myCub->getNumPoints(); // retrieve number of cubature points
Basis_HGRAD_TET_C1_FEM<double, FieldContainer<double> > tetBasis; // create tet basis int numFields = tetBasis.getCardinality(); // get basis cardinality
FieldContainer<double> cub_points(numCubPoints, spaceDim); FieldContainer<double> cub_weights(numCubPoints); FieldContainer<double> cell_nodes(numCells, numNodes, spaceDim); FieldContainer<double> jacobian(numCells, numCubPoints, spaceDim, spaceDim); FieldContainer<double> jacobian_inv(numCells, numCubPoints, spaceDim, spaceDim); FieldContainer<double> jacobian_det(numCells, numCubPoints); FieldContainer<double> weighted_measure(numCells, numCubPoints); FieldContainer<double> grad_at_cub_points(numFields, numCubPoints, spaceDim); FieldContainer<double> transformed_grad_at_cub_points(numCells, numFields, numCubPoints, spaceDim); FieldContainer<double> weighted_transformed_grad_at_cub_points(numCells, numFields, numCubPoints, spaceDim); FieldContainer<double> stiffness_matrices(numCells, numFields, numFields);
We assume that the array
cell_nodes is filled with nodes defining a set of computational (physical) cells.
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
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
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.
The next release of Intrepid is expected to support Finite Difference and Finite Volume discretizations on standard and non-standard (polygon and polyhedron) cell topologies. A "user-friendly" version for rapid development of PDE codes is also under development.