DOI:10.14311/APP.2021.30.0018
Acta Polytechnica CTU Proceedings 30:18–23, 2021 © Czech Technical University in Prague, 2021

available online at http://ojs.cvut.cz/ojs/index.php/app

FINITE ELEMENT IMPLEMENTATION OF GEOMETRICALLY
NONLINEAR CONTACT

Ondřej Faltus∗, Martin Horák

Czech Technical University in Prague, Faculty of Civil Engineering, Department of Mechanics, Thákurova 7,
166 29 Prague 6, Czech Republic

∗ corresponding author: ondrej.faltus@fsv.cvut.cz

Abstract. The OOFEM finite element software has been recently updated to include contact al-
gorithms for small strain applications. In this work, we attempt to extend the contact algorithms to
large strain problems. Reviewing the current code and comparing it with approaches encountered in
literature, we arrive at a specific algorithmic solution and integrate it into the current code base. The
current code is explained, the necessary extensions are derived and documented, and the algorith-
mic changes are described. Tests confirm the functionality and quadratic rate of convergence of the
proposed implementation.

Keywords: Contact, large strains, node-to-segment, OOFEM, penalty method.

1. Introduction

Having originated with Heinrich Hertz in 1881 [1],
contact mechanics has been a fringe discipline of
structural mechanics for most of the 20th century, ow-
ing mainly to its limited options for analytical analy-
sis. Only with the advent of numerical methods, such
as the finite element method (FEM), has development
in this field sprung forward [2]. This surge in new ad-
vances is only exacerbated by the ever more rapid evo-
lution of computational hardware capabilities experi-
enced today; ever more complicated contact problems
necessitating ever more complicated methods previ-
ously unimaginable due to hardware constraints are
coming into focus [2, 3].

The tool used for the purposes of this article is the
finite element computational software OOFEM [4].
OOFEM is a free, open-source C++ code developed
mainly at the Czech Technical University in Prague.
It is based on object-oriented design with emphasis on
variability and extendability. Currently, OOFEM’s
finite element code supports a number of modules for
solution of various types of numerical problems. Of
interest for our purposes is the sm module intended
for structural mechanics.

Basic introduction of contact algorithms into the
OOFEM software has been achieved in [5]. This
implementation, however, has certain limitations,
among them chiefly the lack of support for large
strain, i.e., geometrically non-linear problems. In this
work, we strive to remove this particular limitation
and demonstrate a simple method of extending exist-
ing basic penalty-based formulation to handle large
strain problems.

2. Algorithm
2.1. Basic equation
Without loss of generality, let us consider two de-
formable bodies Bγ , γ = 1, 2. Moreover, let Ωγ0 and
Ωγ represent the regions occupied by the body in the
reference and deformed configuration, respectively.
Position of a point in the reference configuration is
denoted by X γ , and a position xγ of a point in the
deformed configuration is obtained be means of the
one-to one deformation mapping ϕγ as

xγ = ϕγ (X γ ) (1)

The basic equations for each body under consider-
ation are the following:
• Kinematic equations

F γ (X γ ) =
∂ϕγ (X γ )

∂X γ
∀X γ ∈ Ωγ (2)

• Equilibrium equations

DivP γ (X γ ) + B̄γ (X γ ) = 0 ∀X γ ∈ Ωγ (3)

• Dirichlet boundary conditions

uγ (X γ ) = uγ (X γ ) ∀X γ ∈ Γγu (4)

• Neumann boundary conditions

t
γ (X γ ) = P γ (X γ ) · N̄ γ (X γ ) ∀X γ ∈ Γγσ (5)

• Constitutive equation

P γ (F ) = λγ Tr(Eγ )F γ + 2µγ F γ · Eγ (6)

Here, we use the so-called Saint-Venant Kirchhoff
material model. Note that the contact implementa-
tion is not affected by the choice of constitutive law
and extension to inelastic law is straightforward.

18

http://dx.doi.org/10.14311/APP.2021.30.0018
http://ojs.cvut.cz/ojs/index.php/app


vol. 30/2021 FEM Implementation of Geometrically Nonlinear Contact

In the equations above, F is the deformation gra-
dient, E is the Green-Lagrangian strain tensor, Div
is divergence operator in the reference configuration,
P is the first Piola-Kirchhoff stress, B̄ are prescribed
body forces in the reference configuration, u are pre-
scribed displacements, N̄ is a unit normal vector in
the reference configuration, Γu and Γσ are the Dirich-
let and Neumann boundaries which are mutually dis-
joint, i.e., Γγu ∩ Γγσ = ∅, and Γγu ∪ Γγσ = Γγ , λγ and µγ
are Lame’s parameters.

Moreover, contact conditions introduce essentially
a new set of boundary conditions into the task. Note
that only normal contact is considered in what fol-
lows. The core issue of contact problems can be de-
scribed by means of the contact constraints in the
form of Karush-Kuhn-Tucker optimality conditions
[6]:

tc ≤ 0 gc ≥ 0 gctc = 0 (7)

The first condition represents the fact that contact
forces tc are only compressive, the second condition
is characterized by a gap function g(u) and express a
non-penetration condition. The last equation is sim-
ply complementarity between two distinct and mutu-
ally exclusive states for the system; either contact oc-
curs, and therefore the gap function gc is zero, while a
non-zero compressive contact force is present, or the
contact force is zero and the gap function is larger
than zero.

2.2. Penalty Method
Due to the constraint conditions in form of ineqauli-
ties, see (7), weak form of the above introduced for-
mulation is represented by a variational inequality,
see [2]. Several methods how to handle the contact
constraints are possible, e.g., the Lagrangian multi-
plier method, the augmented Lagrangian method, or
the Nitsche method [3, 6]. Currently in OOFEM,
the Lagrangian multiplier method and the penalty
method are implemented [5]. It is the penalty method
which we have selected for the purpose of large strain
extension, as its formulation is more straightforward
and the necessary code complexity is notably lower.

For the penalty method, the weak form can be writ-
ten as

!

γ

"

Ωγ
δF γ : P γ − δuγ · B̄γ dV −

−
"

Γγσ
δuγ · t̄γ dΓ + W c = 0 (8)

where δ represents variation and W c is the penalty-
based contact contribution along contact interface Γc
defined as

W c(u) =
"

Γc

1
2

H(−g(u))pg2c (u) dΓ (9)

where H is the Heaviside function and p is the penalty
parameter. The main disadvantages of the penalty

method is the that the achieved solution is imprecise
(it would only be precise for p → ∞), but moreover,
any improvement in precision acts directly against
the stability of the computation, as larger penalty
parameters introduce larger members into computa-
tional stiffness matrices, which disturbs most numer-
ical solvers [2].

Those disadvantages are balanced by the relatively
simple formulation of the method and, in compari-
son with the Lagrangian multiplier method, by sim-
ple resulting equation systems, as no new variables
are introduced into the computation.

2.3. Node-to-Segment Discretization
An implementation of contact problems into the
framework of the finite element method necessi-
tates selecting a distinctive approach at contact dis-
cretization. As the physical world is described by
FEM in the form of finite elements and connect-
ing nodes, contact can only be considered between
those. Typical contact discretizations include node-
to-node, node-to-segment, segment-to-segment and
other approaches [3]. At current time, node-to-node
and node-to-segment discretizations are implemented
in OOFEM [5]. We are only extending the node-to-
segment discretization for large strains, as the node-
to-node discretization is by definition inadequate for
this purpose.

In node-to-segment discretization, a single node
comes into contact with a segment. Various objects
may be understood under the term "segment", most
prominently a set of element boundaries (element
edges in 2D domains or element surfaces in 3D do-
mains) or an arbitrarily defined rigid analytical sur-
face.

Contact detection in node-to-segment discretiza-
tion necessitates determining three things [3, 6]:
(1.) The gap function gc, essentially a projection of

the node onto the contact segment, which deter-
mines whether there is contact or not; analogical
to the penetration function in the analytical ap-
proach to contact problems.

(2.) The normal vector in the deformed configuration
n, which determines the orientation of the contact
surface at the contact point. Obtained for most
segments in OOFEM as vector perpendicular to
the contact segment’s surface at the contact point
[5].

(3.) The extended N-matrix N ∗, which ensures al-
location of the external forces and stiffness to the
proper degrees of freedom in the FEM formulation.
When the segment in question is an element bound-
ary, the matrix takes the form [5]

N ∗ =
#
I(d×d) −Ne(xc)

$
(10)

where I(d×d) is a square unit matrix of the appro-
priate dimension, which corresponds to degrees of
freedom of the contacting node, and Ne(xc) is the

19



Ondřej Faltus, Martin Horák Acta Polytechnica CTU Proceedings

matrix of element boundary interpolation functions
evaluated at the contact point.

2.4. FEM Formulation
To integrate contact constraints into the FEM frame-
work, considerations have to be made for the contact
forces and energy discussed above within the princi-
ples of FEM discretization.

For the penalty method, contact energy is ex-
pressed as

W c =
1
2

pg2c (11)

Note that the energy contact contribution above is
calculated only in the contact point, which can be
interpreted as Lobatto integration of (9), see, e.g., [6]
for more details.

Variation of this part of energy yields

δW c = pgcδgc (12)

It is worth noting that the term pg(u) represents the
contact force. Variation of the gap function can be
further expressed as

δgc =
%

N ∗T
n

||n||

&
δdc (13)

where δdc are the variational nodal displacements of
relevant FEM nodes and n is the normal vector in
deformed configuration.

The external forces of contact are therefore

fc = −pgc
%

N ∗T
n

||n||

&
(14)

To obtain the associated tangential stiffness ma-
trix. the internal forces have to be subjected to dif-
ferentiation with respect to nodal displacements dc
[5]. For small strains, it is generally assumed that
the direction of contact does not change throughout
the analysis, and therefore n remains independent of
dc [5, 6].

Kc = −
∂fc
∂dc

= p
%

N ∗T
n

||n||

&T %
N ∗T

n

||n||

&
(15)

2.5. Extension to Large Strain
Problems

Now we stand before a question: which of the previ-
ously described considerations change if we consider
large strains? The mechanics of introducing contact
constraints to FEM remain the same, as well as the
definition of contact energy.

Two critical differences appear. Firstly, it can no
longer be assumed that contact follows the same di-
rection throughout the whole computation. The nor-
mal vector has to be reassessed after each computa-
tional step and also the stiffness matrix has to reflect
this. Secondly, large deformations result in change of
contact point on the segment, which is reflected in
the evaluation of matrix N ∗.

The normal vector is deformation dependent and
can be expressed geometrically [2]

n = t1 × t2 (16)

where × is the vector cross product and ti are tangen-
tial vectors of the element boundary in the deformed
configuration (in the case of an element edge, i.e.,
contact formulation in 2D, one is the tangent and the
other is the e3, that is, a unit vector perpendicular
to the plane).

For the new expression of stiffness, we elect to uti-
lize a formula presented in [2]. Specifically for the
case of node-to-2D-segment penalty method, it reads

K = p
#
Nv N

T
v −

gc
l

#
Bv T

T
v + Tv BTv +

gc
l

Bv B
T
v

$$
(17)

with

Nv = N ∗T
n

||n||
Bv = B∗T

n

||n||
Tv = N ∗T t (18)

where B∗ is the extended B-matrix. Note that the
first member of the stiffness matrix corresponds with
the previously derived term for the small strain case.

For a linear edge in 2D, the extended matrices are
defined as

N ∗ =
'
−1 0 ξ − 1 0 −ξ 0
0 −1 0 ξ − 1 0 −ξ

(
(19)

B∗ = −
∂N ∗

∂ξ
=

'
0 0 −1 0 1 0
0 0 0 −1 0 1

(
(20)

where ξ is the parametric coordinate on the element
edge.

3. Implementation
3.1. Class Structure
The OOFEM code has a strictly object-oriented
structure [4]. The integration of contact into this
structure has been mainly accomplished in [5] and
only slight additions were necessary for the large
strain extension.

Figure 1 displays the relevant parts of OOFEM
class structure together with the classes used for con-
tact.

The general idea of class structure centers around
two main classes: ElementEdgeContactSegment and
Node2SegmentPenaltyContact.

ElementEdgeContactSegment is derived from the
general ContactSegment class and represents element
edges in 2D. When initialized from an OOFEM input
file, the class receives a set of element boundaries.
Internal processes ensure that contact is considered
always with the edge which is closest to the contacting
node [5].

Node2SegmentPenaltyContact is a subclass of
ActiveBoundaryCondition, an OOFEM feature

20



vol. 30/2021 FEM Implementation of Geometrically Nonlinear Contact

Figure 1. A diagram of OOFEM class structure, focused on contact-relevant parts, created with the help of [7].

which allows boundary conditions to contribute by
themselves to the internal forces, external forces, and
stiffness matrix. The utilization of this feature avoids
creating any form of a "contact element", otherwise
prevalent among FEM contact implementations [3].
The contact boundary condition is likewise initialized
from the OOFEM input file with information on the
size of the penalty parameter, a set of nodes and a set
of one or more contact segments. Any node is always
checked for contact with any segment [5].

3.2. Operations
The assembly of external forces and stiffness matrix
takes place in the Node2SegmentPenaltyContact
class, specifically in its overridden methods
assembleVector() and assemble(). For the
determination of gap gc, the contact segment
method computePenetration() is invoked, which
computes it from a purely geometrical consideration,
considering the node to have penetrated the segment
when its deformed coordinates are on a different side
of the segment than the undeformed ones.

The contact segment is also responsible for the
computation of the tangent and normal vectors. The
tangent vector is determined as

t =
∂x

∂ξ
= −B∗seg xN = (21)

=
'
−1 0 1 0
0 −1 0 1

(
)
**+

**,

x1
y1
x2
y2

-
**.

**/
=

0
x2 − x1
y2 − y1

1

where B∗seg is the segment part of the extended B-
matrix and xi, yi coordinates of the segment nodes.
Meanwhile, the normal vector is determined as

n = e3 × t (22)

for large strains and from the
Element::edgeEvalNormal() function for small
strains (which returns the undeformed normal
vector).

Neither of those vectors is normalized. The
normalization of the normal vector happens in
the contact condition, when the Nv , Bv and Tv
matrices are determined according to (18). In
the assembly of stiffness, the contact conditions
asks the contact segment whether the large strain
part shall be assembled. This is accomplished
by the ContactSegment::hasNonLinearGeometry()
function and based on the setting of the element in
contact.

4. Demonstration
To demonstrate the functionality of our approach, we
select a simple computational example in a 2D plane
strain idealization. The mesh used can be seen in
Figure 2a.

A double-clamped beam is positioned above a can-
tilever. The length of the beam is 5 m and length of
the cantilever is 3 m. Both are 0.25 m high with a
gap of 0.25 m between them. The top of the beam is
loaded in every step by an uniform distributed load
of 8000 kN/m, causing vertical deformation and thus

21



Ondřej Faltus, Martin Horák Acta Polytechnica CTU Proceedings

(a) . Undeformed mesh.

(b) . Deformed state and vertical component of stress in step 4 (initiation
of contact between beam and cantilever).

(c) . Deformed state and vertical component of stress in step 200 (end of
analysis).

Figure 2. A computational example demonstrating the code and its convergence.

Iteration Force err. x Disp. err. x Force err. y Disp. err. y

1 5.074e-05 0.000e+00 5.450e-03 0.000e+00
2 7.311e-03 8.044e-02 1.583e-01 4.515e-02
3 1.879e-04 2.920e-03 4.120e-03 2.014e-03
4 1.465e-04 2.844e-03 3.202e-03 1.062e-03
5 2.041e-07 5.071e-04 1.488e-06 1.809e-04
6 1.418e-12 1.863e-06 3.298e-11 7.549e-07

Table 1. Evolution of global force and displacement errors in iterations of loading step 14.

subsequently contact between the two bodies. Load-
ing is performed in 200 uniform loading steps.

Both beam and cantilever are modelled by rect-
angular plane strain elements (OOFEM element
Quad1PlaneStrain). In total, the task mesh con-
sists of 805 elements. Refinement of the mesh has
been accomplished with the help of a custom script
in the MATLAB software [8]. All elements are made
of the same Saint-Venant Kirchhoff material with the
parameters λ = 100 MPa and µ = 150 MPa.

Contact conditions are introduced between the

lower edge of the beam and the upper edge of the can-
tilever below. On the latter interface, two boundary
conditions are defined, one pairing the beam’s nodes
with the cantilever’s element edges and the other vice
versa. This ensures that no clipping caused by crude
or otherwise imperfect meshing can take place.

Figure 2b shows the deformed state of the mesh and
the vertical component of the stress tensor in step 4
of 200, when the beam and the cantilever first come
into contact. Likewise, Figure 2c shows the same in
step 200 at the end of the analysis. To render those

22



vol. 30/2021 FEM Implementation of Geometrically Nonlinear Contact

Figure 3. Evolution of error during the equilibrium
iteration process.

figures, the ParaView open-source software was used
[9]. Values of stress in the visualizations are given in
Pa.

Note that between steps 4 and 200, the contact
surface moves significantly to the left both on the
beam and the cantilever. This was accomplished by
alternating activation of the two opposite boundary
conditions on the interface, which managed to subse-
quently introduce contact between neighboring node-
segment pairs, while the contact gap between the
original pairs opened up again.

The main objective of this test was to prove
quadratic convergence of the new approach to large
strain contact stiffness matrix.

The vast majority of the 200 steps computed man-
aged to converge in 6 iterations or less. The only
exception is step 4, in which the contact occurs for
the first time. Table 1 shows error evolution in
a typical computational step (step 14 was chosen).
Moreover, the dependence of error (norm of residual
forces) on the number of iterations during the equi-
librium Newton-Raphson process for steps 4, 14, and
164 is visualized in Figure 3. Thanks to the semi-
logarithmic scale used, it is apparent that approxi-
mate quadratic convergence has been achieved.

5. Conclusions
The paper describes an implementation of a two
dimensional large strain node-to-segment penalty-
based approach to contact mechanics into an open-
source finite element code OOFEM. The emphasis
is given on proper design of the hierarchy of classes
with a focus on modularity and extensibility towards
different methods for handling contact constraints,
contact with friction, or three dimensional formula-
tion. The functionality of the presented approach
is successfully verified on a numerical example and a
quadratic rate of convergence of the Newton-Raphson
iteration procedure is demonstrated.

Acknowledgements
The authors would like to acknowledge support re-
ceived for work on this project from a CTU grant no.
SGS20/038/OHK1/1T/11.

References
[1] H. Hertz. On the contact of elastic solids. Z Reine

Angew Mathematik 92:156–171, 1881.
[2] P. Wriggers. Computational contact mechanics.

Springer, New York, 2nd edn., c2006.
[3] V. A. Yastrebov. Numerical methods in contact

mechanics. Wiley, Hoboken, NJ, 2013.
[4] B. Patzák. OOFEM home page, 2000.

Http://www.oofem.org.
[5] O. Faltus. Object-Oriented Design and

Implementation of Contact Mechanics into Finite
Element Code OOFEM. Master’s thesis, Czech
Technical University in Prague, 2020.

[6] A. Konyukhov, R. Izi. Introduction to computational
contact mechanics. Wiley, Chichester, West Sussex,
2015.

[7] StarUML. Version 3.2.2. MKLabs Co.,Ltd., Seoul,
Republic of Korea, 2020.

[8] MATLAB. Version 9.3.0 (R2017b). The MathWorks
Inc., Natick, Massachusetts, 2017.

[9] U. Ayachit. The ParaView Guide: A Parallel
Visualization Application. Kitware, 1st edn., 2015.

23