=Paper=
{{Paper
|id=Vol-2930/paper12
|storemode=property
|title=Variational Method for Solving Contact Problem of Elasticity
|pdfUrl=https://ceur-ws.org/Vol-2930/paper12.pdf
|volume=Vol-2930
|authors=Robert Namm,Georgiy Tsoy,Ellina Vikhtenko,Gyungsoo Woo
}}
==Variational Method for Solving Contact Problem of Elasticity==
https://ceur-ws.org/Vol-2930/paper12.pdf
Variational Method for Solving Contact Problem of Elasticity
Robert Namma, Georgiy Tsoya, Ellina Vikhtenkob and Gyungsoo Wooc
a
Computing Center of Far Eastern Branch Russian Academy of Sciences, 65 Kim Yu Chen Street, Khabarovsk,
680000, Russia
b
Pacific National University, 136 Tikhookeanskaya Street, Khabarovsk, 680035, Russia
c
Changwon National University, Changwon, 51140, Republic of Korea
Abstract
Variational inequalities corresponding to nonlinear contact problems in mechanics often arise
in engineering practice. To solve them, duality methods are widely used. As a rule, they are
based on the classical methods of constructing Lagrange functionals with linear dependence in
dual variables. This approach is typical for determining the saddle point β the displacement
vector and normal stress in the contact area. The linear dependence in the dual variables does
not allow to prove the theoretical convergence to the saddle point of the well-known iterative
methods. It is possible to justify the convergence only in the primal variable under the
condition that the shift in the dual variable is sufficiently small.
Keywords 1
Variational inequality, contact problem, modified Lagrange functional, saddle point, Uzawa
algorithm.
1. Introduction
The contact problem for an elastic body with a rigid support is represented in the form of a
variational inequality or an equivalent constrained minimization problem of a convex energy functional.
Using the classical duality scheme, this problem can be reduced to the problem of finding a saddle point
for the Lagrange functional [1-2]. Saddle point search methods for classical Lagrange functionals have
been deeply and in detail investigated in many works [3-5], but as a rule, they guarantee convergence
only in a primal variable. The question of convergence in a dual variable remains open. To overcome
this drawback, a modified duality scheme is considered in the paper, on the basis of which a saddle
point search method is constructed, which guarantees convergence in the dual variable as well.
2. Two-dimensional contact problem of elasticity
Let Ξ© β π
2 be a bounded domain with Lipschitz boundary Ξ.
Figure 1: Elastic body with the contact zone
VI International Conference Information Technologies and High-Performance Computing (ITHPC-2021),
September 14β16, 2021, Khabarovsk, Russia
EMAIL: rnamm@yandex.ru (A. 1); tsoy.dv@mail.ru (A. 2); vikht.el@gmail.com (A. 3); gswoo@changwon.ac.kr (A. 4)
ORCID: 0000-0001-7226-3166 (A. 1); 0000-0002-4209-1284 (A. 2); 0000-0002-7152-2311 (A. 3)
Β©οΈ 2020 Copyright for this paper by its authors.
Use permitted under Creative Commons License Attribution 4.0 International (CC BY 4.0).
CEUR Workshop Proceedings (CEUR-WS.org)
98
We assume that on part of the boundary ΞD the body is rigidly fixed, on part Ξπ the surface forces
are given. The contact zone of an elastic body with a rigid foundation will be denoted by ΞπΎ .
2
For the displacement vector π£ = (π£1 , π£2 ), define the strain tensor π = {πππ }π,π=1
1 ππ£π ππ£π
πππ (π£) = ( + ), π, π = 1,2,
2 ππ₯π ππ₯π
2
and the stress tensor π = {πππ }π,π=1
πππ (π£) = πππππ πππ (π£),
where πΆ = {πππππ} is a given elasticity tensor with the usual properties of positive definiteness and
symmetry πππππ = πππππ = πππππ , π, π, π, π=1,2; πππππ πΌππ πΌππ β₯ π0 πΌππ πΌππ , π0 > 0 β const.
Summation over repeated indices is assumed.
Let us specify vector-functions of the body and surface forces π = (π1 , π2 ) and π = (π1 , π2 ),
respectively. The boundary value problem is formulated as follows
βπππ£ π(π’) = π in Ξ©, (1)
π’ = 0 on Ξπ· , (2)
Ο(π’)π = π on Ξπ , (3)
π’π β€ 0, ππ (π’) β€ 0, ππ (π’)π’π = 0 on ΞπΎ , (4)
ππ = 0 ππ ΞπΎ , (5)
where π = (π1 , π2 ) is the unit outward normal vector to Ξπ , π = (π1 , π2 ) is the unit outward normal
vector to ΞπΎ , ππ (π’) = πππ (π’)ππ ππ , ππ (π’) = π(π’)π β ππ (π’) β π, where ππ (π’) and ππ (π’) are normal and
tangential components of the surface traction on ΞπΎ , respectively.
The boundary value problem (1)-(5) belongs to the class of problems with a free boundary since the
adhesion region (π’π = 0) on ΞπΎ is not known in advance and is found simultaneously with the desired
solution of the problem. Condition (4) means non-penetration of an elastic body into a rigid foundation.
The problem (1)-(5) has a variational formulation. Let π β [πΏ2 (Ξ©)]2 . Define the set of admissible
displacements
2
πΎ = {π£ β [π»Ξ1π· (Ξ©)] : π£π β€ 0 on ΞπΎ },
where π»Ξ1π· (Ξ©) = {π£ β π»1 (Ξ©): π£ = 0 on Ξπ· }.
The problem (1)-(5) corresponds to the variational inequality [6]
π(π’, π£ β π’) β₯ β« π β
(π£ β π’) πΞ© + β« π β
(π£ β π’) πΞ β π£ β πΎ, (6)
Ξ© Ξπ
here π(π’, π£) = β«Ξ© πππ (π£)πππ (π£) πΞ© = β«Ξ© πππππ πππ (π£)πππ (π£) πΞ©.
Variational inequality (6) is equivalent to the minimization problem
π½(π£) β min,
{ (7)
π£ β πΎ,
1
where π½(π£) = 2 π(π’, π£) β β«Ξ© π β
π£ πΞ© β β«Ξ π β
π£ πΞ.
π
It is known that the solution π’ β πΎ to the problem (7) exists and is unique, and it satisfies the
equilibrium equation (1) and boundary conditions (2)-(5) in the generalized sense [6].
Let us formulate the dual problem for the problem (7) using the classical duality scheme. For this,
we define the Lagrange functional
2
πΏ(π£, π) = π½(π£) + β« π π£π πΞ β (π£, π) β [π»Ξ1π· (Ξ©)] Γ πΏ2 (ΞπΎ ) (8)
ΞπΎ
and the corresponding dual functional
π(π) = inf πΏ(π£, π).
π£β[π»Ξ1 (Ξ©)]
2 (9)
π·
Problem
π(π) β sup,
{
π β πΏ+2 (Ξπ ),
+ (Ξ )
where πΏ2 π = {π β πΏ2 (ΞπΎ ): π β₯ 0 on ΞπΎ }, is called dual to the problem (7).
99
2
A pair (π£ β , π β ) β [π»Ξ1π· (Ξ©)] Γ πΏ+2 (ΞπΎ ) is called a saddle point of the Lagrange functional πΏ(π£, π) if
the following two-sided inequality takes place
2
πΏ(π£ β , π) β€ πΏ(π£ β , π β ) β€ πΏ(π£, π β ) β(π£, π) β [π»Ξ1π· (Ξ©)] Γ πΏ+2 (ΞπΎ )
In this case, π£ β is the desired solution π’ to the problem (7), and π β is the solution of the dual problem (9)
and coincides on ΞπΎ with the normal stress ππ (π’).
3. Modified Lagrange functional
As already noted, the solution of the contact problem of the theory of elasticity is closely related to
the search for the saddle point of the classical Lagrange functional. The well-known saddle point search
algorithms for classical Lagrange functionals do not guarantee convergence in dual variables. This
situation occurs, for example, in the well-known Uzawa method [3-4]. To overcome this serious
drawback, let us consider a modified duality scheme that allows one to construct algorithms for finding
saddle points that provide convergence in both primal and dual variables.
Consider the modified Lagrange functional [7-8]
1 2
π(π£, π) = π½(π£) + β« (((π + ππ£π )+ )2 β π 2 )πΞ β (π£, π) β [π»Ξ1π· (Ξ©)] Γ πΏ2 (ΞπΎ ). (10)
2π ΞπΎ
Here (π + ππ£π )+ = πππ₯{0; π + ππ£π }, π > 0 is arbitrary positive constant.
For the modified functional π(π£, π), we define a saddle point as follows.
2
Definition. A pair (π£ β , π β ) β [π»Ξ1π· (Ξ©)] Γ πΏ2 (ΞπΎ ) is called a saddle point of the modified Lagrange
functional π(π£, π) if the following two-sided inequality takes place
2
π(π£ β , π) β€ π(π£ β , π β ) β€ π(π£, π β ) β(π£, π) β [π»Ξ1π· (Ξ©)] Γ πΏ2 (ΞπΎ ). (11)
The definition of the saddle point for the modified Lagrange functional π(π£, π) differs from the
definition of the saddle point for the classical one in that in the two-sided inequality (11) the domain of
variation of the dual variable π coincides with the entire functional space πΏ2 (ΞπΎ ), in contrast to the
corresponding inequality for the classical analogue, whereas the domain variation of the dual variable π
is taken by πΏ+2 (ΞπΎ ). Despite this, the sets of saddle points for πΏ(π£, π) and π(π£, π) coincide. This
important property of the modified Lagrange functionals is provided by the nonlinear dependence of the
π(π£, π) on the dual variable in formula (10).
Introduce the dual functional
π(π) = inf 2 π(π£, π)
π£β[π»Ξ1 (Ξ©)]
π·
and the corresponding dual problem
π(π) β sup,
{ (12)
π β πΏ2 (ΞπΎ ).
The following statement holds [7, 8].
Theorem. The dual functional π(π) is Gateaux differentiable in πΏ2 (ΞπΎ ) and its derivative βπ(π)
satisfies the Lipschitz condition with the constant 1/π, that is
1
ββπ(π1 ) β βπ(π2 ) βπΏ (Ξ ) β€ βπ1 β π2 βπΏ2 (ΞπΎ ) β π1 , π2 β πΏ2 (Ξπ ).
2 πΎ π
It can be shown that βπ(π) = max{π’π , βπ/π} β π β πΏ2 (Ξπ ), where
1
π’ = arg min 2 {π½(π£) + β« (((π + ππ£π )+ )2 β π 2 )πΞ} .
π£β[π» 1 (Ξ©)]
2π ΞπΎ
Ξπ·
To solve the dual problem (12), taking into account the above theorem, we can consider the gradient
method
π π+1 = π π + π βπ(π π ), π = 1,2, β¦
with any initial π 0 β πΏ2 (ΞπΎ ).
100
Since βπ(π π ) = max{π’ππ+1 , βπ π /π}, the gradient method is transformed into the following Uzawa
algorithm
1 + 2 2
(π) π’π+1 =arg min 2 {π½(π£) + β« (((π π + ππ£π ) ) β (π π ) ) πΞ} ,
π£β[π» 1 (Ξ©)]
2π ΞπΎ (13)
Ξπ·
+
(ππ) π π+1 = (π π + π π’ππ+1 ) .
Under the condition of solvability of the dual problem (12), it is possible to prove the weak
convergence of the sequence {π π } generated by the Uzawa algorithm to the solution π β of the dual
problem. In this case, the sequence {π’π } converges to the desired solution π’β with respect to the
minimized functional, that is lim π½(π’π ) = π½(π’β ) [7]. It can be proved, if the sequence {π’π } belongs to
πββ
[π» 2 (Ξ©)]2 and is bounded, then the sequence {(π’π , π π )} converges to (π’β , π β ) at the norm of
2
[π»Ξ1π· (Ξ©)] Γ πΏ2 (ΞπΎ ) and, at the same time, π β = βππ (π’β ).
4. Numerical solution of the contact problem of elasticity
Research on the numerical analysis of variational inequalities is carried out, as a rule, on the basis of
the finite element method [2, 3].
Let Ξ© β π
2 be a bounded polygonal domain. Let us carry out a finite element approximation of
problem (7) using piecewise bilinear basis functions [9]. By standard transformations, problem (7) is
transformed into a finite-dimensional quadratic programming problem
1
π₯(π₯) = < π΄π₯, π₯ > β β (ππ1 π₯π1 + ππ2 π₯π2 ) β β (ππ1 π₯π1 + ππ2 π₯π2 ) β min ,
2 (14)
πβπ© πββ³
{ π(π₯ β
π)π β‘ π₯π1 π1 + π₯π2 π2 β€ 0 π β π«.
where π© is the set of indices of quadrilateral mesh nodes, |π©| is the cardinality of the set π©, β³ is the
set of indices of mesh nodes on Ξπ , π« is the set of mesh nodes indices on ΞπΎ , π΄ = (πππ ), π, π =
1, β¦ ,2|π©| is the stiffness matrix, (ππ1 , ππ2 ), (ππ1 , ππ2 ) are the coordinates of the expansion of the
vectors of volume and surface forces in the finite element basis of each node π, respectively.
Under the natural condition that ππππ Ξπ· > 0, the matrix π΄ is symmetric and positive definite. The
function to be minimized in problem (14) corresponds to a finite-dimensional approximation of the
functional π½(π£) in problem (7). To solve the problem (7), we use the Uzawa algorithm with a modified
Lagrange functional in a finite-dimensional version. Let us apply one of the quadrature formulas for the
finite-dimensional approximation of the expression
1
β« (((π + ππ£π )+ )2 β π 2 )πΞ.
2π ΞπΎ
As a result, we obtain a continuously differentiable piecewise-quadratic function in the variables π₯π1 ,
π₯π2 , π β π«, of the form
1 + 2 2
β (((ππ + π(π₯ β
π)π ) ) β (ππ ) ) βπ ,
2π
πβπ«
where ππ , βπ are known quantities.
Let us set an arbitrary π 0 β π
|π«| , where |π«| is the cardinality of the set π«. Uzawa algorithm (13) in
the finite-dimensional case has the form
1 + 2 2
(π)β² π₯ π+1 =arg min {π₯(π₯) + β (((πππ + π(π₯ β
π)π ) ) β (πππ ) ) βπ } ,
2|π©|
π₯βπ
2π (15)
πβπ«
+
(ππ)β²β² πππ+1 = (πππ + π(π₯ β
π)π ) , π β π«.
Let us consider the step (π)β². It is the problem of minimizing a continuously differentiable piecewise
quadratic function. A feature of the function to be minimized is that its Hessian has discontinuities on
101
some linear manifolds and, at the same time, the gradient of the function is continuous. To minimize
function (16) we apply a natural generalization of Newtonβs method [10], [12]
1 + 2 2
π₯(π₯) + β (((πππ + π(π₯ β
π)π ) ) β (πππ ) ) βπ (16)
2π
πβπ«
under a fixed π π .
Let us assume that at some step π the monotony of generalized Newton method is broken, i.e.
1 + 2 2
π₯(π₯ πβ1 ) + β (((πππ + π(π₯ πβ1 β
π)π ) ) β (πππ ) ) <
2π
πβπ«
1 + 2 2
< π₯(π₯ π ) + β (((πππ + π(π₯ π β
π)π ) ) β (πππ ) ).
2π
πβπ«
+
Let π«Μ β π« and be such that (πππ + π(π₯ πβ1 β
π)π ) = πππ + π(π₯ πβ1 β
π)π for all π β π«Μ . Consider the
convex quadratic function
1 2 2
π₯(π₯) + β ((πππ + π(π₯ β
π)π ) β (πππ ) ). (17)
2π
πβπ«Μ
Its Hessian and gradient will be denoted as π΄πβ1 and π΄πβ1 π₯ β π΅πβ1 , respectively, where π΅πβ1 is
some vector. Taking into account that according to Newtonβs method
π₯ π = π₯ πβ1 β π΄β1 πβ1 (π΄πβ1 π₯
πβ1
β π΅πβ1 ),
then
ββ©(π΄πβ1 π₯ πβ1 β π΅πβ1 ), π₯ π β π₯ πβ1 βͺ = β©π΄πβ1 (π₯ π β π₯ πβ1 ), π₯ π β π₯ πβ1 βͺ > 0
or
β©π΄πβ1 π₯ πβ1 β π΅πβ1 , π₯ π β π₯ πβ1 βͺ < 0.
Thus, in the vicinity of the point π₯ πβ1 , the minimized function (16) decreases locally in the direction
π₯ π β π₯ πβ1. The point π₯ π is the minimum point of the quadratic function (17), but is not the desired
minimum point of the convex function (16). Using the well-known rule of Armijo [11] we find a
number πΌπ > 0 such that the value of the minimized function (17) at the element π₯ πβ1 +
πΌπ (π₯ π β π₯ πβ1 ) will be less that at the element π₯ πβ1. Next, we take π₯ πβ1 + πΌπ (π₯ π β π₯ πβ1 ) as a
new π₯Μ π and return to the usual Newtonian method. We ensured monotonicity in the process of
minimizing the piecewise quadratic continuously differentiable function (16). This ensures the
convergence in the minimized function of the generalized Newtonβs methods, and from strong
convexity, the convergence in the arguments π₯ in a finite number of steps (due to the finiteness of the
set π«). The criterion that the minimum point of function (16) is found is the repeatability of the set π«Μ at
two successive steps of the generalized Newton method. The use of the generalized Newton method at
step (π)β² of the Uzawa algorithm with a modified Lagrange function significantly accelerates the search
for a solution to the problem (14). As a rule, the monotonicity of the minimization process at the (π)β²
step is achieved automatically after the first or second application of the Armijo rule.
Overall, the Uzawa algorithm (15) rapidly converges to a saddle point due to the fast stabilization of
the sequence {π π }.
Note that, within the framework of solving the linear programming problem, several authors
previously investigated similar algorithms for minimizing piecewise quadratic functions [12].
5. Numerical experiments
We shall consider the following model example: the body β¦ = (0, 3) Γ (0, 1) (in m) is made of an
elastic isotropic, homogeneous material characterized by Youngβs modulus E = 21.19x104 MPa and
Poissonβs ratio ΞΌ = 0.277. It is fixed along Ξπ· = {0} Γ (0, 1) and linearly distributed surface tractions
of density π = (π1 , π2 ) are applied of Ξπ = Ξπ1 βͺ Ξπ2 , where Ξπ1 = {3} Γ (0, 1) and Ξπ2 = (0, 3) Γ {1}.
We consider the following traction forces: π1 (x) = 0, π2 (x) = 1 MPa on Ξπ1 , π1 (x) = 0, π2 (x) =
β 2β3 (3 β x1 ) MPa on Ξπ2 , parameter π = 108 .
102
The body is discretized into ππ₯ Γ ππ¦ 4-node quadrilateral finite elements, where ππ₯ (= 3ππ¦ ) is
varied to generate problem instances with different sizes. The number of degrees of freedom of
displacements is ππ = 2(ππ₯ + 1)(ππ¦ + 1) and the number of contact candidate nodes is ππ = ππ₯ + 1.
We assume the small deformation and solve the finite-dimensional optimization problem using the
generalized Newton method (GNM). Finding an inverse matrix is a very computationally expensive
operation. Therefore, rather than computing the inverse of the generalized Hessian matrix, one may
save time and increase numerical stability by solving the system of linear equations. For this purpose,
we use the conjugate gradient method implemented in the SciPy [13] and CuPy [14] packages and
compare the computation time. In all numerical experiments considered below, we choose the following
condition:
βπ₯ π β π₯ πβ1 β2
max (βπΊ(π₯ π )β2 , ) < 10β10
βπ₯ π β2
as a stopping criterion for the generalized Newton method. Uzawa algorithm terminates if
βπ π β π πβ1 β2
< 10β8 .
βπ π β2
All experiments are implemented in Python, using the scikit-fem library [15] for performing finite
element assembly and CuPy library for GPU-accelerated computing. Computation was carried out on
IBM Power Systems S822LC 8335-GTB server, which is based on two 10-core IBM POWER8
processors with a maximum operating frequency of 4.023 GHz and two NVIDIA Tesla P100 GPU
accelerators.
Table 1 shows how the total number of Uzawa method iterations and the number of generalized
Newton method iterations depend on ππ ΠΈ ππ . The number of GNM iterations on the first step of the
Uzawa method is represented by the first integer in the respective column. The second integer
represents the number of iterations in subsequent steps of the Uzawa method. We can see that the
number of iterations slightly increases with the increasing number of primal and dual variables.
Calculations show that the use of CuPy library can speed up the execution time for problems of large
size up to 12.6 times.
Table 1
Computational results
ππ₯ Γ ππ¦ ππ ππ GNM it Uzawa it Time CPU, s Time GPU, s M(uβ , lβ )
60x20 2562 61 7/2 6 1.89 5.06 -6.700472e-05
120x40 9922 121 8/2 7 7.70 8.76 -6.711822e-05
240x80 39042 241 9/2 8 56.35 17.24 -6.714898e-05
480x160 154882 481 10/2 11 480.60 38.06 -6.715744e-05
Figure 2 shows normal and tangential displacements of the body on the contact zone for ππ = 481.
Figure 2: Normal and tangential contact displacements
103
The value of the dual variable (normal contact stress) is depicted in Figure 3. We see that the
penetration of the elastic body into a rigid foundation does not occur and at the points where the body is
in contact, the dual variable is positive. The resulting domain deformation in Lagrange coordinates
x+1000u(x) with an amplification factor 1000 and Von Mises stresses are presented in Figure 4.
Figure 3: Normal contact stress
Figure 4: Distribution of the von Mises stresses
6. Conclusion
In the paper, the numerical algorithm of solving the contact problem was proposed. The algorithm
based on the modified Lagrange functionals and Uzawa method. The algorithm was implemented by
using the finite element method. The numerical experiments illustrating the fast convergence of the
algorithm by primal and dual variables were presented. This circumstance can be explained by the good
differential properties of the modified dual functional, which makes it possible to implement the
gradient method for solving the dual problem. After a discretization, we used the generalized Newton
method with Armijo line search to solve the minimization problem of piecewise quadratic functional.
7. Acknowledgements
This study was supported by the Russian Foundation for Basic Research (Project 20-01-00450 A),
by the Ministry of Science and Higher Education of the Russian Federation, Supplementary Agreement
No. 075-02-2020-1529/1 dated April 21, 2020.
104
The studies were carried out using the resources of the Center for Shared Use of Scientific
Equipment "Center for Processing and Storage of Scientific Data of the Far Eastern Branch of the
Russian Academy of Sciences", funded by the Russian Federation represented by the Ministry of
Science and Higher Education of the Russian Federation under project No. 075-15-2021-663.
8. References
[1] I. HlavaΔhek, Y. Haslinger, I. Nechas, Y. LoviΕ‘hek, Solution of Variational Inequalities in
Mechanics, Springer-Verlag, New York, 1986.
[2] N. Kikuchi, T. Oden, Contact Problem in Elasticity: a Study of Variational Inequalities and Finite
Element Methods, SIAM, Philadelphia, 1988.
[3] R. Glowinski, Numerical Methods for Nonlinear Variational Problems, Springer, New York, 1984.
[4] I. Ekland, R. Temam, Convex Analysis and Variational Problems, North-Holland, Amsterdam,
1976.
[5] R. Glowinski R., J. L. Lions, R. Tremolier, Numerical Analysis of Variational Inequalities, North-
Holland, Amsterdam, 1981.
[6] A. M. Khludnev, Elasticity Problems in Nonsmooth Domains, Fizmatlit, Moscow, 2010.
[7] R. V. Namm, G. I. Tsoy, A modified dual scheme for solving an elastic crack problem, Numer.
Anal. Appl. 10(1) (2017) 37β46. doi.org/10.1134/S1995423917010050.
[8] R. V. Namm, E. M. Vikhtenko, N. N. Maksimova, Sensitivity functionals in variational
inequalities of mechanics and their applications to duality schemes, Numer. Anal. Appl. 7(1)
(2014) 36β44. doi.org/10.1134/S1995423914010042.
[9] P. Ciarlet, The Finite Element Method for Elliptic Problems, North Holland, Amsterdam, 1978.
[10] O. L. Mangasarian, A generalized Newton method for absolute value equations, Optim. Lett. 3(1)
(2009) 101-108. doi.org/10.1007/s11590-008-0094-5.
[11] L. Armijo, Minimization of functions having Lipschitz continuous first partial derivatives, Pacific
J. Math. 16(1) (1966) 1-3. doi.org/10.2140/pjm.1966.16.1.
[12] A. I. Golikov, Yu. G. Evtushenko, Generalized Newtonβs method for linear optimization problems
with inequality constraints, Proc. Steklov Inst. Math. 284 (2014) 96β107.
doi.org/10.1134/S0081543814020096.
[13] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, et al., SciPy
1.0: fundamental algorithms for scientific computing in Python, Nat. Methods 17(3) (2020) 261β
272. doi.org/10.1038/s41592-019-0686-2.
[14] R. Okuta, Y. Unno, D. Nishino, S. Hido, and C. Loomis, CuPy: A NumPy-Compatible Library for
NVIDIA GPU Calculations, in: Proceedings of Workshop on Machine Learning Systems
(LearningSys) in The Thirty-first Annual Conference on Neural Information Processing
Systems(NIPS), ), Long Beach, CA, USA, 2017.
[15] T. Gustafsson, G. D. McBain, scikit-fem: A Python package for finite element assembly, J. Open
Source Softw. 5(52) (2020) 2369. doi.org/10.21105/joss.02369.
105