Next Article in Journal
Sustainability Viewed from Farmers’ Perspectives in a Resource-Constrained Environment
Next Article in Special Issue
Durability Characteristics of Concrete Mixture Based on Red Ceramic Waste Aggregate
Previous Article in Journal
Did Safe Cycling Infrastructure Still Matter During a COVID-19 Lockdown?
Previous Article in Special Issue
An Analysis of Traffic Conflicts as a Tool for Sustainable Road Transport
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Non-Monotone Projected Gradient Method in Linear Elasticity Contact Problems with Given Friction

1
Department of Mathematics, Faculty of Civil Engineering, VSB-TU Ostrava, Ludvíka Podéště 1875/17, 708 00 Ostrava, Czech Republic
2
Department of Applied Mathematics, Faculty of Electrical Engineering and Computer Science, VSB-TU Ostrava, 17 Listopadu 2172/15, 708 00 Ostrava, Czech Republic
3
Institute of Geonics of the Czech Academy of Sciences, Studentská 1768, 708 00 Ostrava, Czech Republic
*
Author to whom correspondence should be addressed.
Sustainability 2020, 12(20), 8674; https://doi.org/10.3390/su12208674
Submission received: 29 September 2020 / Revised: 13 October 2020 / Accepted: 15 October 2020 / Published: 19 October 2020
(This article belongs to the Special Issue Civil Engineering as a Tool for Developing a Sustainable Society)

Abstract

:
We are focusing on the algorithms for solving the large-scale convex optimization problem in linear elasticity contact problems discretized by Finite Element method (FEM). The unknowns of the problem are the displacements of the FEM nodes, the corresponding objective function is defined as a convex quadratic function with symmetric positive definite stiffness matrix and additional non-linear term representing the friction in contact. The feasible set constraints the displacement subject to non-penetration conditions. The dual formulation of this optimization problem is well-known as a Quadratic Programming (QP) problem and can be considered as a most basic non-linear optimization problem. Understanding these problems and the development of efficient algorithms for solving them play the crucial role in the large-scale problems in practical applications. We shortly review the theory and examine the behavior and the efficiency of Spectral Projected Gradient method modified for QP problems (SPG-QP) on the solution of a toy example in MATLAB environment.

1. Introduction

One of the key components of the structural design is the contact problem. This problem is crucial, for example, during the design process of steel components, namely welded and screwed contacts. From the mathematical point of view, we are dealing with variational equalities and variational inequalities [1]. Since the problem is analytically solvable only for simple geometries, the typical approach is to discretize the problem and using the Finite Element Method (FEM, [2,3]). Using the reformulation of the problem by the so-called Ritz–Galerkin method and under certain assumptions, one can mathematically prove that the solution of the discretized problem is the approximation of the contiguous one. Additionally, the error converges asymptotically to zero with the density increase in the used mesh [1,3]. The advantage of the method is the possibility of the formulation of the problem as the (constrained) energy functional minimization problem. In the case of the linear elasticity, this corresponding problem is the Quadratic programming (QP) problem [4].
QP problems belong among to the most important optimization problems in practical applications; for example, in the case of the least square methods for regression analysis with the Euclidean distance [5], the solution of the problems with symmetric positive definite system matrix [6,7], Support Vector Machine (SVM, [8,9]), the inner subproblem of more complex (non-linear) optimization problems [10,11], and the numerical solution of variational inequalities in the linear elasticity contact problems [1,4,12,13].
In our case, the Hessian matrix represents the stiffness matrix and the linear term of the objective function is the force density vector. The unknown is the vector of the displacements in the FEM nodes. In the case of contact problems, the non-penetration conditions of the bodies and obstacles are defining the feasible set of the corresponding optimization problem with the linear inequality constraints. Additionally, if we consider also a friction between the bodies and/or the obstacles, the original contiguous problem captures this property by the new term influencing the solution displacement. For example, in the case of the linear elasticity contact problem with the Tresca (given) friction, the discretized minimization problem of the energy functional with the additional non-linear term representing the response of displacement caused by friction and constraints representing the non-penetration, can be reformulated using the dualization into the QP problem with bound and separable quadratic inequality constraints [1,12,14].
Although the QP problems belong to the basic non-linear optimization problems necessary for understanding non-linear optimization theory in general, there does not exist the best algorithm for solving problems from different applications, especially when one introduces the Dirichlet boundary conditions by equality constraints and/or uses non-overlapping domain decomposition methods, such as the Finite Element Tearing and Interconnecting Method (FETI, [4,15,16,17,18,19,20,21]). The contact problems of linear elasticity are characterized by the bounded spectrum independent of the used mesh [22]. The utilization of this property has been proven and practically applied/demonstrated only in the case of active-set algorithms, namely the MPRGP algorithm (Modified Proportioning with Reduced Gradient Projections, [4,8,12,23,24]). Another approach is to use Interior-Point (IP, [25]) method—we eliminate the inequality constraints using barrier functions and solve the sequence of resulting saddle-point problems. In this paper, we examine the projected gradient descent method [26,27,28], especially SPG-QP (Spectral Projected Gradients [29] for QP [5]).

2. Problem Definition

As a simple benchmark, we consider the block of homogeneous material Ω with fixed zero displacement on boundary Γ D and imposed traction forces F on Γ F . The part Γ C denotes the part of boundary that may get into contact with a rigid obstacle. The block is attracted to the obstacle as consequence of the gravity force F g , see Figure 1.
We solve the discretized form of the problem using FEM (see [2,3,30]), which leads to the optimization problem
u * : = arg min u Ω C f ( u ) + j h ( u ) ,
where N N is a number of used FEM nodes, n = 3 N is a number of primal variables, m c N is a number of FEM nodes in Γ C , the terms of the objective function are defined as
f ( u ) : = 1 2 u T K u f T u , j h ( u ) : = i = 1 m c ψ i T i u 2 ,
and
  • u R n is a vector of unknown displacements in FEM nodes;
  • Ω C R n is a set of feasible (non-penetrated) displacements;
  • K R n , n is a symmetric positive definite (SPD) stiffness matrix (the Dirichlet boundary condition is implemented by the modification of the stiffness matrix);
  • f R n is a vector of forces density resulting from the stresses imposed on the structure during a displacement (the discretized form of F and F g );
  • j h : R n R is a numerical integration of functional describing the friction forces in the weak formulation of the problem;
  • T i R 2 , n , i = 1 , , m c are basis vectors formed by appropriately placed multiples of the unit tangential vectors in such a way that the jump of tangential displacement in i-th FEM node is given by T i u ;
  • ψ i R + , i = 1 , , m are slip bound coefficients associated with T i .
Since our problem has a simple geometry, see Figure 2, we set n : = [ 0 , 0 , 1 ] as the unit normal contact vector and t 1 : = [ 1 , 0 , 0 ] , t 2 : = [ 0 , 1 , 0 ] as the tangential vectors for every FEM node in Γ C .
For every contact node (i-th node of Γ C ) T i R 2 , n is given by a zero matrix with 1 in the first row on appropriate x-coordinate of i-th node and in the second row on appropriate y-coordinate of i-th node. Consequently, the matrix composed as
T : = T 1 T m c R 2 m c , n
is a full rank matrix.
We can modify the non-differentiable term j h in (1) into equivalent form (see [1])
j h ( u ) = i = 1 m c max τ i 2 ψ i τ i T T i u ,
where τ i R 2 are regulation variables. We denote the function and the vector
L ( u , τ ) : = f ( u ) + τ T T u , τ = τ 1 τ m c R 2 m c ,
and we write constraints τ i 2 ψ i in (2) in the form
[ τ i ] 1 2 + [ τ i ] 2 2 ψ i , i = 1 , , m c .
After the substitution into (1), we get
min u Ω C f ( u ) + j h ( u ) = min u Ω C max τ Λ τ L ( u , τ ) ,
where we denoted Λ τ R 2 m c as a feasible set defined by constraints (4). We consider L ( u , τ ) as a Lagrange function and τ as a vector of the Lagrange multipliers (in notation (3)) and we use the duality theorem (see Dostál [4]) to reformulate problem (5) into
min u Ω C max τ Λ τ L ( u , τ ) = max τ Λ τ min u Ω C L ( u , τ ) .
We include condition u Ω C by creating new Lagrange multipliers
max τ Λ τ min u Ω C L ( u , τ ) = max τ Λ τ , λ C 0 min u R n L ( u , τ ) + λ C T ( B u c ) ,
where matrix B R m c , n and vector c R m c are constructed in such way that
Ω C = { u R n : B u c } .
Due to geometry in our problem we can construct B simply; B is zero matrix with 1 in every i-th row (which is corresponding to i-th node in Γ C ) on appropriate z-coordinate of i-th node (see the choice of normal vectors for nodes in Γ C in Figure 2). Using these values, the set (8) consists of the feasible displacements which do not penetrate the rigid obstacle.
Problem (6) is equivalent to the saddle point problem
( u * , λ * ) : = arg max λ Λ min u R n f ( u ) + λ T ( B ^ u c ^ ) ,
where
λ : = τ λ C , B ^ : = T B , c ^ : = 0 c
and
Λ : = { [ τ , λ C ] R 3 m c : τ 2 i 1 2 + τ 2 i 2 ψ i , i = 1 , , m c , λ C 0 } .
We derive the dual problem corresponding to (9). From the first KKT condition (since stiffness matrix is SPD, we can use standard inversion K 1 ) we get
K u f + B ^ T λ = 0 u = K 1 f B ^ T λ
and substitute into Lagrange function (7). After the simplifications, we get
L ( u , λ ) = L K 1 ( f B ^ T λ ) , λ = 1 2 λ T B ^ K 1 T T λ + λ T B ^ K 1 f 1 2 f T K 1 f .
We obtain the function of only one variable λ . Since we want to find maximizer (see saddle-point problem (9)), we omit the constant term and change signs and λ * solves equivalent minimization problem
λ * = arg min λ Λ Θ ( λ ) , Θ ( λ ) : = 1 2 λ T A λ λ T b ,
where we denoted
A : = B ^ K 1 B ^ T , b : = B ^ K 1 f .
After solving the minimization problem (12), the corresponding solution u * of primal problem (1) can be evaluated using (11).
Obviously, A R 3 m c , 3 m c is SPD and problem (12) is the QP problem with separable quadratic constraints combinated with bound constraints. Such a problem always has a unique solution [4].

3. Numerical Solution

To solve the problem (12), we use the Spectral Projected Gradient method (SPG, [29]), which is a projection gradient method for solving the general minimization problems
x ¯ : = arg min x Θ f ( x ) .
The method is based on the Barzilai–Borwein (BB, [31,32]) step-size
x it + 1 2 = P Θ ( x it α it f ( x it ) ) , α it : = x it x it 1 , x it x it 1 x i t x it 1 , f ( x it ) f ( x it 1 ) ,
where f ( x ) R n is a gradient of objective function f : R n R and P Θ : R n Θ is the projection onto closed convex feasible set Θ R n . However, since this procedure does not necessarily provide the decreasing sequence of function values (see [33]), an additional step is introduced
x it + 1 = x it + β it d it , d it = x it + 1 2 x it ,
where d it R n is called spectral projected gradient. The appropriate step-size β it ( 0 , ) is computed using the iterative Grippo–Lampariello–Lucidi line-search method (GLL, [34]) to satisfy so-called generalized Armijo condition
f ( x it + β it d it ) < f max + τ β it f ( x it ) , d it .
Here, τ ( 0 , 1 ) represents a safeguarding parameter and
f max : = max { f ( x it j ) : 0 j min { it , m 1 } }
with predefined constant m N . Because of the enforcement of condition (14), the algorithm generates the sequence of approximations with a decrease in objective function in every m iterations and the algorithm is converging to the minimizer.
Recently, [5] introduced the simplification of the method in the case of the quadratic optimization problem. If the objective function is quadratic, one can derive the analytical formula for the computation of step-size, which satisfies (14)
β it σ 1 , min { σ 2 , β ^ }
with safeguarding parameters 0 < σ 1 < σ 2 < 1 and
β ^ = ( 1 τ ) β ¯ + ( 1 τ ) 2 β ¯ 2 + 2 ξ , ξ = f max f ( x it ) A d it , d it , β ¯ = f ( x it ) , d it A d it , d it ,
where A = 2 f ( x it ) R n , n is a Hessian matrix of quadratic objective function f : R n R . The algorithm can be written with only single matrix-vector multiplication during one iteration, see Algorithm 1.
Algorithm 1: SPG for QP problems (SPG-QP, [5])
Given a quadratic objective function f ( x ) = 1 2 x T A x b T x with A R n , n SPD and b R n , closed convex feasible set Θ R n , initial approximation x 0 Θ , projection onto feasible set P Θ ( x ) , parameters m N , τ ( 0 , 1 ) , safeguarding parameters σ 1 , σ 2 R : 0 < σ 1 < σ 2 < 1 , precision ε > 0 , and initial step-size α 0 > 0 .

   g 0 : = A x 0 b
   f 0 : = 1 / 2 g 0 b , x 0
  set index of iteration i t : = 0

  while stopping criterium is not safisfied
    d i t : = P Θ ( x i t α i t g i t ) x i t

   compute matrix-vector multiplication A d i t
   compute multiple dot-product d i t , { d i t , A d i t , g i t }

    f max : = max { f ( x i t j ) : 0 j min { i t , m 1 } }
    ξ : = ( f max f ( x i t ) ) / d i t , A d i t
    β ¯ : = g i t , d i t / d i t , A d i t
    β ^ : = τ β ¯ + τ 2 β ¯ 2 + 2 ξ
   choose β i t [ σ 1 , min { σ 2 , β ^ } ]

    x i t + 1 : = x i t + β i t d i t
    g i t + 1 : = g i t + β i t A d i t
    f i t + 1 : = f ( x i t ) + β i t d i t , g i t + 1 2 β i t 2 d i t , A d i t

    α i t + 1 : = d i t , d i t / d i t , A d i t

    i t : = i t + 1
  endwhile
Return the approximation of solution x i t and the number of performed iterations it .
We used the presented Spectral Projected Gradient method for QP for solving the dual problem (12) with variable λ R m c , objective function Θ : R 3 m c R , and feasible set (10) defined by separable quadratic constraints and bound constraints.

4. Results

We implement and solve the problem introduced in Section 2 using MATLAB. For the FEM discretization, we adopt the open-source library presented in [35]. Instead of evaluating and computing the inverse of the matrix in the dualization (13), we solve the system of linear equations with multiple right-hand side vectors. In MATLAB, we evaluate
Y = K \ B ^ ; A = B ^ Y ; b = Y T f .
In our numerical experiment, we consider steel brick with E = 2 × 10 5 [ MPa ] , μ = 0.33 , ρ = 7.85 × 10 2 [ kg · m 3 ] of size 2 × 1 × 0.25 [ m ] in contact with the rigid obstacle. The displacement and the friction are caused by the force F = [ 450 , 0 , 0 ] [ MPa ] . We consider given (Tresca) friction between the brick and the obstacle. To discretize the problem, we construct regular grid with the number of FEM nodes in axis equal to N x = 25 , N y = 13 , N z = 4 . In total, we obtain N = 1300 FEM nodes, the primal problem of dimension n = 3900 , the number of FEM nodes in contact m c = 325 , and dual problem of 975 unknowns. The feasible set of dual problem is composed from 325 bound constraints and 325 separable quadratic constraints.
We implemented SPG-QP in MATLAB and set the parameters of the Armijo criterion to recommended values m = 10 , τ = 0.9 and initial approximation equal to the (feasible) zero vector. The first BB step-size is equal to steepest descent (Cauchy) step-size. We solved the problem (see Figure 3 and Figure 4) with the coefficient of friction (slip bound) ψ i = 100 S i [ MPa ] , where S i is the contact surface corresponding to i-th node. Figure 5 demonstrates the decrease in the Euclidean norm of scaled projected gradient
g α ¯ P x : = 1 α ¯ P Θ ( x α ¯ f ( x ) ) x ,
where α ¯ 0 , 1 / λ max A (see [36]) and λ max A is an upper estimation of the largest eigenvalue of the Hessian matrix A , which we compute using Gershgorin circle theorem. The norm of the projected gradient is used as a stopping criterion of the algorithm—i.e., we stop the algorithm if
d it , d it < ε .
In the SPG-QP algorithm, we decided to choose the largest possible β it to perform the steps in as large a manner as possible. This approach has been suggested by Fletcher [33].
Figure 6 shows the dependency of the number of iterations on the value of friction coefficient (the size of quadratic constraints) and the size of the problem. In the first case, we solve the problem of dimension N = 1300 and absolute error ε = 10 6 for the norm of scaled projected gradient. In the second case, we are solving the problem with friction coefficient ψ ^ = 100 [ MPa ] while varying the problem size using the same stopping criterion.
To improve the performance of the solution process, we utilize the GPU capability of the MATLAB. In our case, we compute on GeForce RTX 2080 Ti (4352 cuda cores, 11 G B GDDR6 memory) and compare with Intel Core i9-7900X in MATLAB R2017a. Since our code is highly vectorized, the only prerequisite to perform the computation on GPU is to transfer the data using MATLAB command gpuArray. The overloaded MATLAB operations perform the matrix and vector operations on GPU. To analyze the performance of this approach, we fix the number of iterations of the SPG-QP algorithm to be equal to 10 4 and measure the computational time. From the results presented in Figure 7, it is clear that the main computational bottleneck is not the solution of dual QP problem, but the dualization—i.e., the solution of the system of linear Equations (15). The results suggest the suitability of used gradient projected method on GPU; however, the problem has to be sufficiently large to efficiently use this architecture. The performance of SPG-QP is limited by the matrix–vector multiplication with dense dual Hessian matrix, which is the most time-consuming operation. Consequently, the time complexity of the method scales quadratically with the dimension of dual problem, see Figure 7.

5. Conclusions

We examined SPG-QP algorithm for solving the strictly convex QPQC in the solution of linear elasticity contact problem with given friction. The results demonstrate the suitability of the algorithm for this type of problems. However, with the increasing size of the primal problem, we observe the increase in the computational complexity of the corresponding Hessian matrix inverse (which has polynomial complexity, in general). We suggested using domain decomposition methods, for instance the Finite Element Tearing and Interconnecting method (FETI). This extension will be the topic of our upcoming research.

Author Contributions

Conceptualization, L.P. and D.H.; methodology, L.P., M.Č., J.K. and D.H.; software, L.P. and M.Č.; validation, D.H. and J.K.; formal analysis, L.P. and D.H.; investigation, L.P. and J.K.; resources, D.H., L.P. and M.Č.; data curation, L.P. and M.Č.; writing—original draft preparation, L.P.; writing—review and editing, L.P. and M.Č.; visualization, L.P.; supervision, D.H. and M.Č.; project administration, M.Č. and D.H.; funding acquisition, M.C. and D.H. All authors have read and agreed to the published version of the manuscript.

Funding

This paper has been completed thanks to the financial support provided to VSB Technical University of Ostrava by the Czech Ministry of Education, Youth and Sports from the budget for conceptual development of science, research and innovations for the 2020 year.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Hlaváček, I.; Haslinger, J.; Nečas, J.; Lovíšek, J. Solution of Variational Inequalities in Mechanics; Springer: Berlin, Germany, 1988. [Google Scholar]
  2. Reddy, J.N. An Introduction to the Finite Element Method; McGraw-Hill Education: New York, NY, USA, 1993. [Google Scholar]
  3. Brenner, S.; Scott, R. The Mathematical Theory of Finite Element Methods, 3rd ed.; Springer: Berlin/Heidelberg, Germany, 2008. [Google Scholar]
  4. Dostál, Z. Optimal Quadratic Programming Algorithms, with Applications to Variational Inequalities; SOIA, Springer: New York, NY, USA, 2009; Volume 23. [Google Scholar]
  5. Pospíšil, L.; Gagliardini, P.; Sawyer, W.; Horenko, I. On a scalable nonparametric denoising of time series signals. Commun. Appl. Math. Comput. Sci. 2018, 13, 107–138. [Google Scholar] [CrossRef]
  6. Hestenes, M.R.; Stiefel, E. Methods of conjugate gradients for solving linear systems. J. Res. Natl. Bur. Stand. 1952, 49, 409–436. [Google Scholar] [CrossRef]
  7. Dostál, Z.; Pospíšil, L. Conjugate gradients for symmetric positive semidefinite least-squares problems. Int. J. Comput. Math. 2018, 95, 2229–2239. [Google Scholar] [CrossRef]
  8. Kružík, J.; Pecha, M.; Hapla, V.; Horák, D.; Čermák, M. Investigating Convergence of Linear SVM Implemented in PermonSVM Employing MPRGP Algorithm. In Proceedings of the International Conference on High Performance Computing in Science and Engineering, Karolinka, Czech Republic, 22–25 May 2017; pp. 115–129. [Google Scholar] [CrossRef]
  9. Pecha, M.; Čermák, M. Analyzing l1-loss and l2-loss Support Vector Machines Implemented in PERMON Toolbox. Int. Conf. Adv. Eng. Theory Appl. 2018, 13–23. [Google Scholar] [CrossRef]
  10. Nocedal, J.; Wright, S.J. Numerical Optimization, 2nd ed.; Springer: New York, NY, USA, 2006. [Google Scholar]
  11. Boyd, S.; Vandenberghe, L. Convex Optimization; Cambridge University Press: Cambridge, UK, 2004. [Google Scholar]
  12. Dostál, Z.; Pospíšil, L. Optimal iterative QP and QPQC algorithms. Ann. Oper. Res. 2016, 243, 5–18. [Google Scholar] [CrossRef]
  13. Hapla, V.; Horák, D.; Pospíšil, L.; Čermák, M.; Vašatová, A.; Sojka, R. Solving Contact Mechanics Problems with PERMON. In High Performance Computing in Science and Engineering; Kozubek, T., Blaheta, R., Šístek, J., Rozložník, M., Čermák, M., Eds.; Springer International Publishing: Zurich, Switzerland, 2016; pp. 101–115. [Google Scholar]
  14. Dostál, Z.; Horák, D.; Kučera, R.; Vondrák, V.; Haslinger, J.; Dobiáš, J.; Pták, S. FETI based algorithms for contact problems: Scalability, large displacements and 3D coulomb friction. Comput. Methods Appl. Mech. Eng. 2005, 194, 395–409. [Google Scholar] [CrossRef]
  15. Farhat, C.; Roux, F.X. A method of finite element tearing and interconnecting and its parallel solution algorithm. Int. J. Numer. Methods Eng. 1991, 32, 1205–1227. [Google Scholar] [CrossRef]
  16. Dostál, Z.; Horák, D.; Kučera, R. Total FETI—An easier implementable variant of the FETI method for numerical solution of elliptic PDE. Commun. Numer. Methods Eng. 2006, 22, 1155–1162. [Google Scholar] [CrossRef]
  17. Kruis, J. Domain Decomposition Methods on Parallel Computers. In Progress in Engineering Computational Technology; Topping, B.H.V., Mota Soares, C.A., Eds.; Saxe-Coburg Publications: Stirling, UK, 2004; pp. 299–322. ISBN 1-874672-22-9. [Google Scholar]
  18. Medek, O.; Kruis, J.; Bittnar, Z.; Tvrdík, P. Static Load Balancing Applied to Schur Complement Method. Comput. Struct. 2007, 85, 489–498. [Google Scholar] [CrossRef]
  19. Kruis, J. The FETI Method and its Applications: A Review. In Parallel, Distributed and Grid Computing for Engineering; Topping, B., Iványi, P., Eds.; Saxe-Coburg Publications: Stirling, UK, 2009; pp. 199–216. ISBN 978-1-874672-41-8/1759-3158. [Google Scholar]
  20. Kruis, J.; Matouš, K.; Dostál, Z. Solving laminated plates by domain decomposition. Adv. Eng. Softw. 2002, 33, 445–452. [Google Scholar] [CrossRef]
  21. Kruis, J. Domain Decomposition Methods for Distributed Computing; Saxe-Coburg Publications: Stirling, UK, 2006; ISBN-13: 978-1-874672-23-4. [Google Scholar]
  22. Dostál, Z.; Kučera, R. An Optimal Algorithm for Minimization of Quadratic Functions with Bounded Spectrum Subject to Separable Convex Inequality and Linear Equality Constraints. SIAM J. Optim. 2010, 20, 2913–2938. [Google Scholar] [CrossRef]
  23. Bouchala, J.; Dostál, Z.; Kozubek, T.; Pospíšil, L.; Vodstrčil, P. On the solution of convex QPQC problems with elliptic and other separable constraints. Appl. Math. Comput. 2014, 247, 848–864. [Google Scholar] [CrossRef]
  24. Dostál, Z.; Kozubek, T.; Markopoulos, A.; Brzobohatý, T.; Vondrák, V.; Horyl, P. A theoretically supported scalable TFETI algorithm for the solution of multibody 3D contact problems with friction. Comput. Methods Appl. Mech. Eng. 2012, 110–120. [Google Scholar] [CrossRef]
  25. Wright, S.J. Primal-Dual Interior-Point Methods; Society for Industrial and Applied Mathematics: Philadelphia, PA, USA, 1997. [Google Scholar]
  26. Pospíšil, L.; Dostál, Z. The Projected Barzilai-Borwein Method with Fall-Back for Strictly Convex QCQP Problems with Separable Constraints. Math. Comput. Simul. 2018, 145, 79–89. [Google Scholar] [CrossRef]
  27. Dai, Y.H.; Fletcher, R. Projected Barzilai-Borwein Methods for Large-Scale Box-Constrained Quadratic Programming. Numer. Math. 2005, 100, 21–47. [Google Scholar] [CrossRef]
  28. Nesterov, Y.E. A method for solving the convex programming problem with convergence rate O(1/k2). Dokl. Akad. Nauk SSSR 1983, 269, 543–547. [Google Scholar]
  29. Birgin, E.G.; Martínez, J.M.; Raydan, M.M. Nonmonotone spectral projected gradient methods on convex sets. SIAM J. Optim. 2000, 10, 1196–1211. [Google Scholar] [CrossRef]
  30. Alberty, J.; Carstensen, C.; Funken, S.A. Remarks around 50 line of Matlab: Short finite element implementation. Numer. Algorithms 1999, 20, 117–137. [Google Scholar] [CrossRef]
  31. Barzilai, J.; Borwein, J.M. Two point step size gradient methods. IMA J. Numer. Anal. 1988, 8, 141–148. [Google Scholar] [CrossRef]
  32. Raydan, M.M. Convergence Properties of the Barzilai and Borwein Gradient Method. Ph.D. Thesis, Rice University, Houston, TX, USA, 1991. [Google Scholar]
  33. Fletcher, R. On the Barzilai-Borwein Method. In Optimization and Control with Applications; Springer: Boston, MA, USA, 2005; pp. 235–256. [Google Scholar]
  34. Grippo, L.; Lampariello, F.; Lucidi, S. A Nonmonotone Line Search Technique for Newton’s Method. SIAM J. Numer. Anal. 1986, 23, 707–716. [Google Scholar] [CrossRef]
  35. Čermák, M.; Sysala, S.; Valdman, J. Efficient and flexible MATLAB implementation of 2D and 3D elastoplastic problems. Appl. Math. Comput. 2019, 355, 595–614. [Google Scholar] [CrossRef] [Green Version]
  36. Dostál, Z. On the decrease of a quadratic function along the projected-gradient path. Electron. Trans. Numer. Anal. 2008, 31, 25–29. [Google Scholar]
Figure 1. Contact problem with rigid obstacle and given friction.
Figure 1. Contact problem with rigid obstacle and given friction.
Sustainability 12 08674 g001
Figure 2. Normal and tangential vectors of FEM nodes in Γ C .
Figure 2. Normal and tangential vectors of FEM nodes in Γ C .
Sustainability 12 08674 g002
Figure 3. The solution of the problem: the displacement.
Figure 3. The solution of the problem: the displacement.
Sustainability 12 08674 g003
Figure 4. The solution of the problem: the Eulidean norm of Friction forces computed from the corresponding Lagrange multipliers of dual quadratic constraints. The arrows represent the relative size and the direction of the force in FEM nodes of the contact zone.
Figure 4. The solution of the problem: the Eulidean norm of Friction forces computed from the corresponding Lagrange multipliers of dual quadratic constraints. The arrows represent the relative size and the direction of the force in FEM nodes of the contact zone.
Sustainability 12 08674 g004
Figure 5. The decrease in objective function (left) and the norm of scaled projected gradient (right) during the solution process. The notation x * is used for the solution obtained in the last iteration.
Figure 5. The decrease in objective function (left) and the norm of scaled projected gradient (right) during the solution process. The notation x * is used for the solution obtained in the last iteration.
Sustainability 12 08674 g005
Figure 6. The number of iterations for the problems with different friction coefficients; (left) the number of iterations and the solution time depending on the problem size (right).
Figure 6. The number of iterations for the problems with different friction coefficients; (left) the number of iterations and the solution time depending on the problem size (right).
Sustainability 12 08674 g006
Figure 7. The computational time: the dualization by evaluation of (15) (left), and the performance of 10 4 iterations of SPG-QP on one CPU and one GPU (right).
Figure 7. The computational time: the dualization by evaluation of (15) (left), and the performance of 10 4 iterations of SPG-QP on one CPU and one GPU (right).
Sustainability 12 08674 g007
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Pospíšil, L.; Čermák, M.; Horák, D.; Kružík, J. Non-Monotone Projected Gradient Method in Linear Elasticity Contact Problems with Given Friction. Sustainability 2020, 12, 8674. https://doi.org/10.3390/su12208674

AMA Style

Pospíšil L, Čermák M, Horák D, Kružík J. Non-Monotone Projected Gradient Method in Linear Elasticity Contact Problems with Given Friction. Sustainability. 2020; 12(20):8674. https://doi.org/10.3390/su12208674

Chicago/Turabian Style

Pospíšil, Lukáš, Martin Čermák, David Horák, and Jakub Kružík. 2020. "Non-Monotone Projected Gradient Method in Linear Elasticity Contact Problems with Given Friction" Sustainability 12, no. 20: 8674. https://doi.org/10.3390/su12208674

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop