# Building Blocks#

Have a look at the building blocks of FElupe.

create a meshed cube with hexahedron elements

setup your own numeric region with a mesh, an element and a quadrature

add a displacement field to a field container

define your own Neo-Hookean material formulation

apply your own boundary conditions

solve the problem (create your own Newton-Rhapson iteration loop)

export the displaced mesh along with cauchy stress projected to mesh-points

Start setting up a problem in FElupe by the creation of a numeric **Region** with a geometry (**Mesh**), a finite **Element** and a **Quadrature** rule, e.g. for hexahedrons or tetrahedrons. By using a template region like RegionHexahedron (see section **Region**), only the mesh has to be created.

```
import felupe as fem
mesh = fem.Cube(n=9)
element = fem.Hexahedron()
quadrature = fem.GaussLegendre(order=1, dim=3)
mesh.plot().show()
```

An overview of the mesh-metadata is available in the console.

```
mesh
```

```
<felupe Mesh object>
Number of points: 729
Number of cells:
hexahedron: 512
```

## Region#

A region essentially pre-calculates element shape functions and derivatives evaluated at every quadrature point of every cell w.r.t. the undeformed coordinates (as attribute dhdX). An array containing products of quadrature weights multiplied by the determinants of the (geometric) jacobians is stored as the differential volume. The sum of all differential volumes gives the total (undeformed) volume of the region. The attributes of a region are used in a Field.

```
region = fem.Region(mesh, element, quadrature)
# region = fem.RegionHexahedron(mesh)
region
```

```
<felupe Region object>
Element formulation: Hexahedron
Quadrature rule: GaussLegendre
Gradient evaluated: True
```

```
dV = region.dV
V = dV.sum()
```

The scheme of the region, the finite element formulation and the quadrature rule, may be visualized by its plot-method.

```
region.plot().show()
```

## Field#

In a second step fields are added to the Region which may be either scalar or vector fields. The values at mesh-points are obtained with the attribute `values`

. Interpolated field values at quadrature points are calculated with the `interpolate()`

method. Additionally, the displacement gradient w.r.t. the undeformed coordinates is calculated for every quadrature point of every cell in the region with the field method `grad()`

. A generalized extraction method `extract(grad=True, add_identity=True, sym=False)`

allows several arguments to be passed. This involves or whether the gradient or the values are extracted. If the gradient is extracted, the identity matrix may be added to the gradient (useful for the calculation of the deformation gradient). Optionally, the symmetric part is returned (small strain tensor).

```
displacement = fem.Field(region, dim=3)
u = displacement.values
ui = displacement.interpolate()
dudX = displacement.grad()
```

Next, the field is added to a field container, which handles one or several (vector) fields. Like a field, the field container also provides the `extract(grad=True, add_identity=True, sym=False)`

method, returning a list of interpolated field values or gradients.

```
field = fem.FieldContainer([displacement])
field
```

```
<felupe FieldContainer object>
Number of fields: 1
Dimension of fields:
Field: 3
```

The deformation gradient is obtained by a sum of the identity and the displacement gradient.

```
F = field.extract(grad=True, sym=False, add_identity=True)
```

## Constitution#

The material behavior has to be provided by the first Piola-Kirchhoff stress tensor as a function of the deformation gradient. FElupe provides a very basic hard-coded constitutive library (Neo-Hooke, linear elasticity and a generalized Hu-Washizu (u,p,J) three field variation). Alternatively, an isotropic material formulation is defined by a strain energy density function - both variation (stress) and linearization (elasticity) are carried out by automatic differentiation using tensortrax. The latter one is demonstrated here with a nearly-incompressible version of the Neo-Hookean material model.

```
import tensortrax.math as tm
def W(C, mu, bulk):
"Neo-Hooke"
J = tm.sqrt(tm.linalg.det(C))
return mu / 2 * (J ** (-2 / 3) * tm.trace(C) - 3) + bulk * (J - 1) ** 2 / 2
umat = fem.Hyperelastic(W, mu=1.0, bulk=2.0)
P = umat.gradient
A = umat.hessian
```

## Boundary Conditions#

Next we enforce boundary conditions on the displacement field. Boundaries are stored as a dictionary of multiple boundary instances. First, the left end of the cube is fixed. Displacements on the right end are fixed in directions y and z whereas displacements in direction x are prescribed with a user-defined value. A boundary instance hold useful attributes like `points`

or `dof`

.

```
import numpy as np
f0 = lambda x: np.isclose(x, 0)
f1 = lambda x: np.isclose(x, 1)
boundaries = {}
boundaries["left"] = fem.Boundary(displacement, fx=f0)
boundaries["right"] = fem.Boundary(displacement, fx=f1, skip=(1,0,0))
boundaries["move"] = fem.Boundary(displacement, fx=f1, skip=(0,1,1), value=0.5)
```

## Partition of deegrees of freedom#

The separation of active and inactive degrees of freedom is performed by a so-called **partition**. External values of prescribed displacement degrees of freedom are obtained by the application of the boundary values to the displacement field.

```
dof0, dof1 = fem.dof.partition(field, boundaries)
ext0 = fem.dof.apply(field, boundaries, dof0)
```

## Integral forms of equilibrium equations#

The integral (or weak) forms of equilibrium equations are defined by the `fem.IntegralForm`

class. The pre-evaluated function of interest has to be passed as the fun argument whereas the virtual field as the `v`

argument. By setting `grad_v=[True]`

(default), FElupe passes the gradient of the virtual field to the integral form. FElupe assumes a linear form if `u=None`

(default) or creates a bilinear form if a field is passed to the field argument `u`

.

```
linearform = fem.IntegralForm(P(F)[:-1], field, dV, grad_v=[True])
bilinearform = fem.IntegralForm(A(F), field, dV, u=field, grad_v=[True], grad_u=[True])
```

The assembly of both forms lead to the (point-based) internal force vector and the (sparse) stiffness matrix.

```
r = linearform.assemble()
K = bilinearform.assemble()
```

## Prepare (partition) and solve the linearized equation system#

In order to solve the linearized equation system a partition into active and inactive degrees of freedom has to be performed. This system may then be passed to the (sparse direct) solver. Given a set of nonlinear equilibrium equations \(\boldsymbol{g}\) the unknowns \(\boldsymbol{u}\) are found by linearization at a valid initial state of equilibrium and an iterative Newton-Rhapson solution prodecure. The incremental values of inactive degrees of freedom are given as the difference of external prescribed and current values of unknowns. The (linear) solution is equal to the first result of a Newton-Rhapson iterative solution procedure. The resulting point values `du`

are finally added to the displacement field.

The default solver of FElupe is SuperLU provided by the sparse package of SciPy. A significantly faster alternative is pypardiso which may be installed from PyPI with `pip install pypardiso`

(not included with FElupe). The optional argument `solver`

of `fem.solve.solve()`

accepts a user-defined solver.

```
from scipy.sparse.linalg import spsolve # default
# from pypardiso import spsolve
system = fem.solve.partition(field, K, dof1, dof0, r)
dfield = fem.solve.solve(*system, ext0, solver=spsolve)#.reshape(*u.shape)
# field += dfield
```

A very simple newton-rhapson code looks like this:

```
for iteration in range(8):
F = field.extract()
linearform = fem.IntegralForm(P(F)[:-1], field, dV)
bilinearform = fem.IntegralForm(A(F), field, dV, field)
r = linearform.assemble()
K = bilinearform.assemble()
system = fem.solve.partition(field, K, dof1, dof0, r)
dfield = fem.solve.solve(*system, ext0, solver=spsolve)
norm = np.linalg.norm(dfield)
print(iteration, norm)
field += dfield
if norm < 1e-12:
break
```

```
0 8.174180680860706
1 0.2940958778404007
2 0.02083230945148839
3 0.0001028992534421267
4 6.017153213511068e-09
5 5.675484825228616e-16
```

Alternatively, one may also use the Newton-Rhapson function of FElupe.

```
res = fem.newtonrhapson(
field, dof1=dof1, dof0=dof0, ext0=ext0, kwargs={"umat": umat}
)
field = res.x
```

All 3x3 components of the deformation gradient of integration point 1 of cell 1 (Python is 0-indexed) are obtained with

```
F = F[0]
F[:,:,0,0]
```

```
array([[ 1.49186831e+00, -1.17603278e-02, -1.17603278e-02],
[ 3.09611695e-01, 9.73138551e-01, 8.43648336e-04],
[ 3.09611695e-01, 8.43648336e-04, 9.73138551e-01]])
```

## Export of results#

Results are exported as VTK or XDMF files using meshio.

```
fem.save(region, field, filename="result.vtk")
```

Any tensor at quadrature points shifted or projected to, both averaged at mesh-points is evaluated for `quad`

and `hexahedron`

cell types by `fem.topoints`

or `fem.project`

, respectively. For example, the calculation of the cauchy stress involves the conversion from the first Piola-Kirchhoff stress to the Cauchy stress followed by the shift or the projection. The stress results at mesh points are passed as a dictionary to the `point_data`

argument.

```
from felupe.math import dot, det, transpose, tovoigt
s = dot(P([F])[0], transpose(F)) / det(F)
cauchy_shifted = fem.topoints(s, region)
cauchy_projected = fem.project(tovoigt(s), region)
fem.save(
region,
field,
filename="result_with_cauchy.vtk",
point_data={
"CauchyStressShifted": cauchy_shifted,
"CauchyStressProjected": cauchy_projected,
}
)
```