# How to use the geometry within an element

(Difference between revisions)
 Revision as of 20:25, 9 July 2014 (view source)Rrossi (Talk | contribs) (Created page with "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 derivat...") Revision as of 20:32, 9 July 2014 (view source)Rrossi (Talk | contribs) Newer edit → Line 15: Line 15: rGeom.ShapeFunctionsIntegrationPointsGradients(DN_DX,DetJ,GeometryData::GI_GAUSS_2); 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 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. + 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, or to a similar precision for the case of triangles and tetrahedra. GI_GAUSS_1,2,3,4 are available, providing an increasing precision at the price of using more integration points 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" + + '''DN_DX[IntPointNumber](I,k)''' contains 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 + One may obtain the Shape functions at the gauss points by using + Matrix NContainer = rGeom.ShapeFunctionsValues(GeometryData::GI_GAUSS_2); + + Ncontainer is defined such that Ncontainer(IntPointNumber, I) contains the shape functions of node I evaluated at the integration point "IntPointNumber" shape functions of node I evaluated at the integration point "IntPointNumber" − NContainer = rGeom.ShapeFunctionsValues(GeometryData::GI_GAUSS_2); + − + Note that if the parameter "GeometryData::GI_GAUSS_2" can be omitted in all of the functions. In this case the default integration rule will be used. Such integration rule will provide exact results − obtain the list of integration points + for the computation of a laplacian or of a the stiffness matrix for a linear structural problem. + + the list of integration points can be queried by const GeometryType::IntegrationPointsArrayType& IntegrationPoints = rGeom.IntegrationPoints(GeometryData::GI_GAUSS_2); const GeometryType::IntegrationPointsArrayType& IntegrationPoints = rGeom.IntegrationPoints(GeometryData::GI_GAUSS_2); Line 32: Line 39: for (unsigned int g = 0; g < IntegrationPoints.size(); g++) for (unsigned int g = 0; g < IntegrationPoints.size(); g++) area += DetJ[g] * IntegrationPoints[g].Weight(); area += DetJ[g] * IntegrationPoints[g].Weight(); + + similarly the geometry allows the computation of the Jacobians on all of the gauss points as well as of the "Local" ShapeFunctionDerivatives, expressed in terms of parent coordinates.

## Revision as of 20:32, 9 July 2014

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;
```

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, or to a similar precision for the case of triangles and tetrahedra. 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 the derivative of the shape function for node I, with respect to the "k" component, evaluated at the gauss point "IntPointNumber"

One may obtain the Shape functions at the gauss points by using

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

Ncontainer is defined such that Ncontainer(IntPointNumber, I) contains the shape functions of node I evaluated at the integration point "IntPointNumber"

Note that if the parameter "GeometryData::GI_GAUSS_2" can be omitted in all of the functions. In this case the default integration rule will be used. Such integration rule will provide exact results for the computation of a laplacian or of a the stiffness matrix for a linear structural problem.

the list of integration points can be queried by

```  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();
```

similarly the geometry allows the computation of the Jacobians on all of the gauss points as well as of the "Local" ShapeFunctionDerivatives, expressed in terms of parent coordinates.