How to use the geometry within an element

The integration of FE functions, typically requires the evaluation of functions at given integration points as well as the computation of the shape functions and of their derivates on such points

the following code snippet is useful to compute some elementary FE operations, abstracting the computations from the actual choice of the geometry. The idea is that in general it is possible to write an element in a similar way for, say, a triangle and a quadrilatera.

If we are within an element, we can obtain a reference to the geometry (read as topology if you like) of the element by doing

  const GeometryType& rGeom = this->GetGeometry();

The geometry provides facilities to compute Jacobians, determinants, shape functions, etc. For example to compute the Determinant of the Jacobian AND the shape function derivatives at once for all of the gauss points of interest, one can do

  Vector DetJ;
  ShapeFunctionDerivativesArrayType DN_DX;
  rGeom.ShapeFunctionsIntegrationPointsGradients(DN_DX,DetJ,GeometryData::GI_GAUSS_2);

the vector "DetJ" will contain the determinant of the Jacobian for all of the integration points considered in the integration rule considered. The choice of the integration rule is controlled here by the parameter"GeometryData::GI_GAUSS_2" which here corresponds to using the 2nd order integration rule, which typically implies 2 gauss points per direction if kroneker product can be used. GI_GAUSS_1,2,3,4 are available, providing an increasing precision at the price of using more integration points DN_DX[IntPointNumber](I,k) contains here the derivative of the shape function for node I, with respect to the "k" component, evaluated at the gauss point "IntPointNumber"

here obtain a matrix Ncontainer such that Ncontainer(IntPointNumber, I) contains the shape functions of node I evaluated at the integration point "IntPointNumber"

  NContainer = rGeom.ShapeFunctionsValues(GeometryData::GI_GAUSS_2);
   

obtain the list of integration points

  const GeometryType::IntegrationPointsArrayType& IntegrationPoints = rGeom.IntegrationPoints(GeometryData::GI_GAUSS_2);

this can be used for example to compute the total area of the element as

   double area = 0.0;
   for (unsigned int g = 0; g < IntegrationPoints.size(); g++)
       area += DetJ[g] * IntegrationPoints[g].Weight();