1 /* Siconos is a program dedicated to modeling, simulation and control
2 * of non smooth dynamical systems.
3 *
4 * Copyright 2021 INRIA.
5 *
6 * Licensed under the Apache License, Version 2.0 (the "License");
7 * you may not use this file except in compliance with the License.
8 * You may obtain a copy of the License at
9 *
10 * http://www.apache.org/licenses/LICENSE-2.0
11 *
12 * Unless required by applicable law or agreed to in writing, software
13 * distributed under the License is distributed on an "AS IS" BASIS,
14 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
15 * See the License for the specific language governing permissions and
16 * limitations under the License.
17 */
18 #include <float.h> // for DBL_EPSILON
19 #include <math.h> // for sqrt, fabs, isinf
20 #include <stdio.h> // for NULL, printf
21 #include <stdlib.h> // for calloc, realloc
22 #include "AlartCurnierGenerated.h" // for fc3d_AlartCurn...
23 #include "FrictionContactProblem.h" // for FrictionContac...
24 #include "Friction_cst.h" // for SICONOS_FRICTI...
25 #include "NonSmoothNewton.h" // for nonSmoothDirec...
26 #include "NumericsFwd.h" // for SolverOptions
27 #include "NumericsMatrix.h" // for NumericsMatrix
28 #include "SolverOptions.h" // for SolverOptions
29 #include "fc3d_2NCP_Glocker.h" // for NCPGlocker_ini...
30 #include "fc3d_AlartCurnier_functions.h" // for compute_rho_sp...
31 #include "fc3d_GlockerFischerBurmeister_functions.h" // for F_GlockerFisch...
32 #include "fc3d_Solvers.h" // for FreeSolverNSGSPtr
33 #include "fc3d_compute_error.h" // for fc3d_unitary_c...
34 #include "fc3d_local_problem_tools.h" // for fc3d_local_pro...
35 #include "fc3d_onecontact_nonsmooth_Newton_solvers.h" // for computeNonsmoo...
36 #include "fc3d_projection.h" // for fc3d_projectio...
37 #include "numerics_verbose.h" // for numerics_print...
38 #include "op3x3.h" // for cpy3, mvp3x3
39 #include "SiconosBlas.h" // for cblas_ddot
40 #include "NSSTools.h" // for max
41
42 /* #define DEBUG_CHECK */
43 /* #define DEBUG_NOCOLOR */
44 /* #define DEBUG_MESSAGES */
45 /* #define DEBUG_STDOUT */
46 #include "siconos_debug.h" // for DEBUG_PRINTF
47
48 #ifdef DEBUG_MESSAGES
49 #include "NumericsVector.h"
50 #endif
51 static computeNonsmoothFunction Function = NULL;
52 static NewtonFunctionPtr F = NULL;
53 static NewtonFunctionPtr jacobianF = NULL;
54 static UpdateSolverPtr updateSolver = NULL;
55 static PostSolverPtr postSolver = NULL;
56 static FreeSolverNSGSPtr freeSolver = NULL;
57
58 /* size of a block */
59 static int Fsize;
fc3d_AC_initialize(FrictionContactProblem * problem,FrictionContactProblem * localproblem,SolverOptions * options)60 static void fc3d_AC_initialize(FrictionContactProblem* problem,
61 FrictionContactProblem* localproblem,
62 SolverOptions * options)
63 {
64 /** In initialize, these operators are "connected" to their corresponding static variables,
65 * that will be used to build local problem for each considered contact.
66 * Local problem is built during call to update (which depends on the storage type for M).
67 */
68
69 DEBUG_PRINTF("fc3d_AC_initialize starts with options->iparam[SICONOS_FRICTION_3D_NSN_FORMULATION] = %i\n",
70 options->iparam[SICONOS_FRICTION_3D_NSN_FORMULATION]);
71
72 if(options->iparam[SICONOS_FRICTION_3D_NSN_FORMULATION] ==
73 SICONOS_FRICTION_3D_NSN_FORMULATION_ALARTCURNIER_STD)
74 {
75 Function = &(computeAlartCurnierSTD);
76 }
77 else if(options->iparam[SICONOS_FRICTION_3D_NSN_FORMULATION] ==
78 SICONOS_FRICTION_3D_NSN_FORMULATION_JEANMOREAU_STD)
79 {
80 Function = &(computeAlartCurnierJeanMoreau);
81 }
82 else if(options->iparam[SICONOS_FRICTION_3D_NSN_FORMULATION] ==
83 SICONOS_FRICTION_3D_NSN_FORMULATION_ALARTCURNIER_GENERATED)
84 {
85 Function = &(fc3d_AlartCurnierFunctionGenerated);
86 }
87 else if(options->iparam[SICONOS_FRICTION_3D_NSN_FORMULATION] ==
88 SICONOS_FRICTION_3D_NSN_FORMULATION_JEANMOREAU_GENERATED)
89 {
90 Function = &fc3d_AlartCurnierJeanMoreauFunctionGenerated;;
91 }
92 else if(options->iparam[SICONOS_FRICTION_3D_NSN_FORMULATION] ==
93 SICONOS_FRICTION_3D_NSN_FORMULATION_NULL)
94 {
95 Function = NULL;
96 }
97
98 /* Compute and store default value of rho value */
99 size_t nc = problem->numberOfContacts;
100
101 double avg_rho[3] = {0.0, 0.0, 0.0};
102
103 if(options->solverId == SICONOS_FRICTION_3D_ONECONTACT_NSN ||
104 options->solverId == SICONOS_FRICTION_3D_ONECONTACT_NSN_GP)
105 {
106 if(!options->dWork ||
107 options->dWorkSize < 3*nc)
108 {
109 options->dWork = (double *)realloc(options->dWork,
110 3*nc * sizeof(double));
111 options->dWorkSize = 3*nc ;
112 }
113 }
114 else if(options->solverId == SICONOS_FRICTION_3D_ONECONTACT_NSN_GP_HYBRID)
115 {
116 if(!options->dWork ||
117 options->dWorkSize < 4*nc)
118 {
119 options->dWork = (double *)realloc(options->dWork,
120 4*nc * sizeof(double));
121 options->dWorkSize = 4*nc ;
122 }
123 }
124
125
126
127 double* rho=0;
128 for(size_t contact =0; contact <nc ; contact++)
129 {
130 if(options->solverId == SICONOS_FRICTION_3D_ONECONTACT_NSN ||
131 options->solverId == SICONOS_FRICTION_3D_ONECONTACT_NSN_GP)
132 {
133 rho = &options->dWork[3*contact];
134 }
135 else if(options->solverId == SICONOS_FRICTION_3D_ONECONTACT_NSN_GP_HYBRID)
136 {
137 options->dWork[contact] = 1.0; // for PLI algorithm.
138 rho = &options->dWork[3*contact+nc];
139 }
140 numerics_printf("fc3d_AC_initialize "" compute rho for contact = %i",contact);
141
142 if(options->iparam[SICONOS_FRICTION_3D_NSN_RHO_STRATEGY] == SICONOS_FRICTION_3D_NSN_FORMULATION_RHO_STRATEGY_SPLIT_SPECTRAL_NORM_COND)
143 {
144 fc3d_local_problem_fill_M(problem, localproblem, contact);
145 compute_rho_split_spectral_norm_cond(localproblem, rho);
146 }
147 else if(options->iparam[SICONOS_FRICTION_3D_NSN_RHO_STRATEGY] == SICONOS_FRICTION_3D_NSN_FORMULATION_RHO_STRATEGY_SPLIT_SPECTRAL_NORM)
148 {
149 fc3d_local_problem_fill_M(problem, localproblem, contact);
150 compute_rho_split_spectral_norm(localproblem, rho);
151 }
152 else if(options->iparam[SICONOS_FRICTION_3D_NSN_RHO_STRATEGY] == SICONOS_FRICTION_3D_NSN_FORMULATION_RHO_STRATEGY_SPECTRAL_NORM)
153 {
154 fc3d_local_problem_fill_M(problem, localproblem, contact);
155 compute_rho_spectral_norm(localproblem, rho);
156 }
157 else if(options->iparam[SICONOS_FRICTION_3D_NSN_RHO_STRATEGY] == SICONOS_FRICTION_3D_NSN_FORMULATION_RHO_STRATEGY_CONSTANT)
158 {
159 rho[0]=options->dparam[SICONOS_FRICTION_3D_NSN_RHO];
160 rho[1]=options->dparam[SICONOS_FRICTION_3D_NSN_RHO];
161 rho[2]=options->dparam[SICONOS_FRICTION_3D_NSN_RHO];
162 }
163 else if(options->iparam[SICONOS_FRICTION_3D_NSN_RHO_STRATEGY] == SICONOS_FRICTION_3D_NSN_FORMULATION_RHO_STRATEGY_ADAPTIVE)
164 {
165 numerics_error("fc3d_AC_initialize", "Adaptive strategy for computing rho not yet implemented");
166 }
167 else
168 numerics_error("fc3d_AC_initialize", "unknown strategy for computing rho");
169
170
171 if(verbose >0)
172 {
173 avg_rho[0] += rho[0];
174 avg_rho[1] += rho[1];
175 avg_rho[2] += rho[2];
176 }
177 numerics_printf("fc3d_AC_initialize""contact = %i, rho[0] = %4.2e, rho[1] = %4.2e, rho[2] = %4.2e", contact, rho[0], rho[1], rho[2]);
178
179 fc3d_local_problem_fill_M(problem, localproblem, contact);
180 double m_row_norm = 0.0, sum;
181 for(int i =0; i<3; i++)
182 {
183 sum =0.0;
184 for(int j =0; j<3; j++)
185 {
186 sum += fabs(localproblem->M->matrix0[i+j*3]);
187 }
188 m_row_norm = max(sum, m_row_norm);
189 }
190 numerics_printf("fc3d_AC_initialize" " inverse of norm of M = %e", 1.0/hypot9(localproblem->M->matrix0));
191 numerics_printf("fc3d_AC_initialize" " inverse of row norm of M = %e", 1.0/m_row_norm);
192
193 DEBUG_EXPR(NM_display(localproblem->M););
194
195 }
196 numerics_printf("fc3d_AC_initialize" " Avg. rho value = %e\t%e\t%e\t",avg_rho[0]/nc,avg_rho[1]/nc,avg_rho[2]/nc);
197
198 }
199
fc3d_AC_free(FrictionContactProblem * problem,FrictionContactProblem * localproblem,SolverOptions * localsolver_options)200 static void fc3d_AC_free(FrictionContactProblem * problem, FrictionContactProblem * localproblem, SolverOptions* localsolver_options)
201 {
202 free(localsolver_options->dWork);
203 localsolver_options->dWork=NULL;
204 }
205
206
fc3d_AC_post(int contact,double * reaction)207 static void fc3d_AC_post(int contact, double* reaction)
208 {
209 /* This function is required in the interface but useless in Alart-Curnier case */
210 }
211
fc3d_onecontact_nonsmooth_Newton_solvers_initialize(FrictionContactProblem * problem,FrictionContactProblem * localproblem,SolverOptions * localsolver_options)212 void fc3d_onecontact_nonsmooth_Newton_solvers_initialize(FrictionContactProblem* problem,
213 FrictionContactProblem* localproblem,
214 SolverOptions * localsolver_options)
215 {
216
217 /* Initialize solver (Connect F and its jacobian, set local size ...) according to the chosen formulation. */
218
219 /* Alart-Curnier formulation */
220 if(localsolver_options->solverId == SICONOS_FRICTION_3D_ONECONTACT_NSN)
221 {
222 fc3d_AC_initialize(problem, localproblem,localsolver_options);
223 postSolver = &fc3d_AC_post;
224 freeSolver = &fc3d_AC_free;
225 }
226 else if(localsolver_options->solverId == SICONOS_FRICTION_3D_ONECONTACT_NSN_GP)
227 {
228 fc3d_AC_initialize(problem, localproblem,localsolver_options);
229 postSolver = &fc3d_AC_post;
230 freeSolver = &fc3d_AC_free;
231
232 }
233 else if(localsolver_options->solverId == SICONOS_FRICTION_3D_ONECONTACT_NSN_GP_HYBRID)
234 {
235
236 fc3d_AC_initialize(problem, localproblem,localsolver_options);
237 postSolver = &fc3d_AC_post;
238 freeSolver = &fc3d_AC_free;
239 }
240 /* Glocker formulation - Fischer-Burmeister function used in Newton */
241 else if(localsolver_options->solverId == SICONOS_FRICTION_3D_NCPGlockerFBNewton)
242 {
243 Fsize = 5;
244 NCPGlocker_initialize(problem, localproblem);
245 F = &F_GlockerFischerBurmeister;
246 jacobianF = &jacobianF_GlockerFischerBurmeister;
247 /* updateSolver = &NCPGlocker_update; */
248 postSolver = &NCPGlocker_post;
249 freeSolver = (FreeSolverNSGSPtr)&NCPGlocker_free;
250 }
251 else
252 {
253 numerics_error("fc3d_onecontact_nonsmooth_Newton_solvers_initialize", "Unknown formulation type.");
254 }
255 }
256
257
258
259
fc3d_onecontact_nonsmooth_Newton_solvers_solve(FrictionContactProblem * localproblem,double * local_reaction,SolverOptions * options)260 int fc3d_onecontact_nonsmooth_Newton_solvers_solve(FrictionContactProblem* localproblem, double* local_reaction, SolverOptions * options)
261 {
262 numerics_printf_verbose(2, "--------------- fc3d_onecontact_nonsmooth_Newton_solvers_solve starts");
263 numerics_printf_verbose(2, "-- contact %i", options->iparam[SICONOS_FRICTION_3D_CURRENT_CONTACT_NUMBER]);
264
265 /* (*updateSolver)(contact, local_reaction); */
266 int info =1;
267
268 /* check trivial solution */
269
270 double* q = localproblem->q;
271 if(q[0] > -DBL_EPSILON)
272 {
273 local_reaction[0] = 0.0;
274 local_reaction[1] = 0.0;
275 local_reaction[2] = 0.0;
276 numerics_printf_verbose(2,"take off, trivial solution reaction = 0, velocity = q.");
277 info=0;
278 options->iparam[SICONOS_IPARAM_ITER_DONE] = 0;
279 options->dparam[SICONOS_DPARAM_RESIDU] = 0.0;
280 numerics_printf_verbose(2, "--------------- fc3d_onecontact_nonsmooth_Newton_solvers_solve ends");
281 return info;
282 }
283
284 if(options->solverId == SICONOS_FRICTION_3D_ONECONTACT_NSN)
285 {
286 info = fc3d_onecontact_nonsmooth_Newton_solvers_solve_direct(localproblem, local_reaction, options);
287 }
288 else if(options->solverId == SICONOS_FRICTION_3D_ONECONTACT_NSN_GP)
289 {
290 info = fc3d_onecontact_nonsmooth_Newton_solvers_solve_damped(localproblem, local_reaction, options);
291 }
292 else if(options->solverId == SICONOS_FRICTION_3D_ONECONTACT_NSN_GP_HYBRID)
293 {
294 if(options->iparam[SICONOS_FRICTION_3D_NSN_HYBRID_STRATEGY] == SICONOS_FRICTION_3D_NSN_HYBRID_STRATEGY_PLI_NSN_LOOP ||
295 options->iparam[SICONOS_FRICTION_3D_NSN_HYBRID_STRATEGY] == SICONOS_FRICTION_3D_NSN_HYBRID_STRATEGY_NSN_AND_PLI_NSN_LOOP)
296 {
297 info = fc3d_onecontact_nonsmooth_Newton_solvers_solve_hybrid(localproblem, local_reaction, options);
298 }
299 else
300 {
301 numerics_error("fc3d_onecontact_nonsmooth_Newton_solvers_solve","Unknown local nsn hybrid solver");
302 }
303 }
304 else
305 {
306 info = nonSmoothDirectNewton(Fsize, local_reaction, &F, &jacobianF, options);
307 }
308 if(info > 0)
309 {
310 if(verbose > 0)
311 {
312 if(options->iparam[SICONOS_IPARAM_MAX_ITER] == options->iparam[SICONOS_IPARAM_ITER_DONE])
313 {
314 numerics_warning("fc3d_onecontact_nonsmooth_Newton_solvers_solve",
315 "reached max. number of iterations (%i) for contact %i. Residual = %12.8e",
316 options->iparam[SICONOS_IPARAM_MAX_ITER], options->iparam[SICONOS_FRICTION_3D_CURRENT_CONTACT_NUMBER], options->dparam[SICONOS_DPARAM_RESIDU]);
317 }
318 else
319 {
320 numerics_warning("fc3d_onecontact_nonsmooth_Newton_solvers_solve",
321 "no convergence for contact %i with error = %12.8e",
322 options->iparam[SICONOS_IPARAM_MAX_ITER], options->iparam[SICONOS_FRICTION_3D_CURRENT_CONTACT_NUMBER], options->dparam[SICONOS_DPARAM_RESIDU]);
323 }
324 /* note : exit on failure should be done in DefaultCheckSolverOutput */
325 }
326 }
327 numerics_printf_verbose(2, "--------------- fc3d_onecontact_nonsmooth_Newton_solvers_solve ends");
328 return info;
329 /* (*postSolver)(contact,reaction); */
330 }
331
fc3d_onecontact_nonsmooth_Newton_solvers_free(FrictionContactProblem * problem,FrictionContactProblem * localproblem,SolverOptions * localsolver_options)332 void fc3d_onecontact_nonsmooth_Newton_solvers_free(FrictionContactProblem * problem, FrictionContactProblem * localproblem, SolverOptions* localsolver_options)
333 {
334 F = NULL;
335 jacobianF = NULL;
336 updateSolver = NULL;
337 postSolver = NULL;
338 (*freeSolver)(problem, localproblem, localsolver_options);
339 }
340
341
fc3d_onecontact_nonsmooth_Newton_solvers_computeError(int n,double * velocity,double * reaction,double * error)342 void fc3d_onecontact_nonsmooth_Newton_solvers_computeError(int n, double* velocity, double*reaction, double * error)
343 {
344 /* int numberOfContacts = n/3; */
345 /* int sizeGlobal = numberOfContacts*FSize; */
346 /* // double * FGlobal = (double*)malloc(sizeGlobal*sizeof(*FGlobal)); */
347 /* (*computeFGlobal)(reaction,velocity); */
348 /* int i; */
349 /* double Fz; */
350 /* *error = 0; */
351 /* for(i=0;i<sizeGlobal;++i) */
352 /* { */
353 /* Fz = velocity[i]*reaction[i]; */
354 /* if(Fz>0) */
355 /* *error+=Fz; */
356 /* if(reaction[i]<0) */
357 /* *error+=reaction[i]; */
358 /* if(velocity[i]<0) */
359 /* *error+=velocity[i]; */
360 /* } */
361
362 /* // (*computeVelocity)(FGlobal); */
363
364 /* free(FGlobal); */
365
366 }
367 #ifdef DEBUG_CHECK
fc3d_onecontact_nonsmooth_Newton_AC_debug(double * R,double * velocity,double mu,double * rho,double * MLocal,double * F,double * A,double * B,double * AWplusB,int * iparam)368 static int fc3d_onecontact_nonsmooth_Newton_AC_debug(double *R, double * velocity, double mu, double * rho, double * MLocal,
369 double * F, double * A, double *B, double * AWplusB, int * iparam)
370 {
371 double AWpB[9];
372 if(iparam[SICONOS_FRICTION_3D_NSN_FORMULATION] !=
373 SICONOS_FRICTION_3D_NSN_FORMULATION_JEANMOREAU_GENERATED
374 && iparam[SICONOS_FRICTION_3D_NSN_FORMULATION] !=
375 SICONOS_FRICTION_3D_NSN_FORMULATION_NULL)
376 {
377 double Fg[3] = {0., 0., 0.};
378 double Ag[9] = {0., 0., 0., 0., 0., 0., 0., 0., 0.};
379 double Bg[9] = {0., 0., 0., 0., 0., 0., 0., 0., 0.};
380 assert(*rho > 0. && *(rho + 1) > 0. && *(rho + 2) > 0.);
381 if(iparam[SICONOS_FRICTION_3D_NSN_FORMULATION] ==
382 SICONOS_FRICTION_3D_NSN_FORMULATION_ALARTCURNIER_STD)
383 {
384 fc3d_AlartCurnierFunctionGenerated(R, velocity, mu, rho, Fg, Ag, Bg);
385 }
386
387 if(iparam[SICONOS_FRICTION_3D_NSN_FORMULATION] ==
388 SICONOS_FRICTION_3D_NSN_FORMULATION_JEANMOREAU_STD)
389 {
390 fc3d_AlartCurnierJeanMoreauFunctionGenerated(R, velocity, mu, rho, Fg, Ag, Bg);
391
392 }
393
394 sub3(F, Fg);
395 sub3x3(A, Ag);
396 sub3x3(B, Bg);
397
398 assert(hypot3(Fg) <= 1e-7);
399 assert(hypot9(Ag) <= 1e-7);
400 assert(hypot9(Bg) <= 1e-7);
401 cpy3x3(A, Ag);
402 cpy3x3(B, Bg);
403 mm3x3(A, MLocal, AWpB);
404 add3x3(B, AWpB);
405
406 scal3x3(-1., AWpB);
407 sub3x3(AWplusB, AWpB);
408 assert(hypot9(AWpB) <= 1e-7);
409 }
410
411 return 0;
412 }
413 #endif
414
fc3d_onecontact_nonsmooth_Newton_AC_update(int contact,FrictionContactProblem * problem,FrictionContactProblem * localproblem,double * reaction,SolverOptions * options)415 void fc3d_onecontact_nonsmooth_Newton_AC_update(int contact, FrictionContactProblem* problem, FrictionContactProblem* localproblem, double * reaction, SolverOptions* options)
416 {
417 /* Build a local problem for a specific contact
418 reaction corresponds to the global vector (size n) of the global problem.
419 */
420 /* Call the update function which depends on the storage for MGlobal/MBGlobal */
421 /* Build a local problem for a specific contact
422 reaction corresponds to the global vector (size n) of the global problem.
423 */
424
425 /* The part of MGlobal which corresponds to the current block is copied into MLocal */
426 fc3d_local_problem_fill_M(problem, localproblem, contact);
427
428 /**** Computation of qLocal = qBlock + sum over a row of blocks in MGlobal of the products MLocal.reactionBlock,
429 excluding the block corresponding to the current contact. ****/
430 fc3d_local_problem_compute_q(problem, localproblem, reaction, contact);
431
432 /* Friction coefficient for current block*/
433 localproblem->mu[0] = problem->mu[contact];
434
435
436 }
437
438
fc3d_onecontact_nonsmooth_Newton_solvers_solve_direct(FrictionContactProblem * localproblem,double * R,SolverOptions * options)439 int fc3d_onecontact_nonsmooth_Newton_solvers_solve_direct(FrictionContactProblem* localproblem,
440 double * R, SolverOptions * options)
441 {
442
443 int * iparam = options->iparam;
444 double * dparam = options->dparam;
445
446 numerics_printf_verbose(2, "--------------- fc3d_onecontact_nonsmooth_Newton_solvers_solve_direct starts");
447
448 double mu = localproblem->mu[0];
449 double * qLocal = localproblem->q;
450
451 double norm_qLocal = sqrt(qLocal[0]*qLocal[0]+qLocal[1]*qLocal[1]+qLocal[2]*qLocal[2]);
452 double norm_relative = 1.0;
453 if(norm_qLocal> DBL_EPSILON)
454 {
455 norm_relative /= norm_qLocal;
456 }
457
458 double * MLocal = localproblem->M->matrix0;
459
460
461 /* store the increment */
462 double dR[3] = {0., 0., 0.};
463 /* store the value fo the function */
464 double F[3] = {0., 0., 0.};
465 /* Store the (sub)-gradient of the function */
466 double A[9] = {0., 0., 0., 0., 0., 0., 0., 0., 0.};
467 double B[9] = {0., 0., 0., 0., 0., 0., 0., 0., 0.};
468 /* Value of AW+B */
469 double AWplusB[9] = {0., 0., 0., 0., 0., 0., 0., 0., 0.};
470
471 /* retrieve value of rho */
472 double * rho;
473 if(options->solverId == SICONOS_FRICTION_3D_ONECONTACT_NSN ||
474 options->solverId == SICONOS_FRICTION_3D_ONECONTACT_NSN_GP)
475 {
476 rho = &options->dWork[3*iparam[SICONOS_FRICTION_3D_CURRENT_CONTACT_NUMBER]];
477 }
478 else
479 {
480 int nc = options->dWorkSize/4;
481 rho = &options->dWork[3*iparam[SICONOS_FRICTION_3D_CURRENT_CONTACT_NUMBER]+nc];
482 }
483
484 /* compute the velocity */
485 double velocity[3] = {0., 0., 0.};
486 cpy3(qLocal,velocity);
487 mvp3x3(MLocal,R,velocity);
488
489
490 int itermax = iparam[SICONOS_IPARAM_MAX_ITER];
491
492 /* Newton iteration */
493 int inew;
494 int info_solv3x3;
495 for(inew = 0 ; inew < itermax ; ++inew)
496 {
497 /* Update function and gradient */
498 Function(R, velocity, mu, rho, F, A, B);
499
500 /* compute -(A MLocal +B) */
501 mm3x3(A, MLocal, AWplusB);
502 add3x3(B, AWplusB);
503 scal3x3(-1., AWplusB);
504
505 #ifdef DEBUG_CHECK
506 fc3d_onecontact_nonsmooth_Newton_AC_debug(R, velocity, mu, rho, MLocal,
507 F, A, B, AWplusB, iparam);
508 #endif
509
510 /* Solve the linear system */
511 cpy3(F,dR);
512 info_solv3x3 = solve_3x3_gepp(AWplusB, dR);
513
514 /* if determinant is zero, replace dR=NaN with zero (i.e. don't modify R) and return early */
515 if(info_solv3x3)
516 {
517 dR[0] = 0;
518 dR[1] = 0;
519 dR[2] = 0;
520 }
521
522 /* update iterates */
523 add3(dR, R);
524
525 /* compute new residue */
526 cpy3(qLocal,velocity);
527 mvp3x3(MLocal,R,velocity);
528 Function(R, velocity, mu, rho, F, NULL, NULL);
529 dparam[SICONOS_DPARAM_RESIDU] = 0.5 * (F[0] * F[0] + F[1] * F[1] + F[2] * F[2])* norm_relative ; // improve with relative tolerance
530
531 if(info_solv3x3)
532 {
533 if(verbose > 0)
534 numerics_warning("fc3d_onecontact_nonsmooth_Newton_solvers_solve_direct",
535 "contact %i 3x3 linear system is irregular # iteration = %"
536 " error = %.10e \n",
537 iparam[SICONOS_FRICTION_3D_CURRENT_CONTACT_NUMBER],
538 inew, dparam[SICONOS_DPARAM_RESIDU]);
539 break;
540 }
541
542 if(verbose > 1) printf("--------------- fc3d_onecontact_nonsmooth_Newton_solvers_solve_direct -- contact %i # iteration = %i error = %.10e \n", iparam[SICONOS_FRICTION_3D_CURRENT_CONTACT_NUMBER], inew, dparam[SICONOS_DPARAM_RESIDU]);
543
544 if(dparam[SICONOS_DPARAM_RESIDU] < dparam[SICONOS_DPARAM_TOL])
545 {
546 iparam[SICONOS_IPARAM_ITER_DONE]=inew;
547 return 0;
548 }
549 }// End of the Newton iteration
550
551 iparam[SICONOS_IPARAM_ITER_DONE]=inew;
552 return 1;
553 }
554
555 #undef DEBUG_MESSAGES
LineSearchGP(FrictionContactProblem * localproblem,computeNonsmoothFunction Function,double * t_opt,double R[3],double dR[3],double * rho,int LSitermax,double * F,double * A,double * B,double * velocity)556 static int LineSearchGP(FrictionContactProblem* localproblem,
557 computeNonsmoothFunction Function,
558 double * t_opt,
559 double R[3],
560 double dR[3],
561 double *rho,
562 int LSitermax,
563 double * F,
564 double * A,
565 double * B,
566 double * velocity)
567 {
568 numerics_printf_verbose(2, "-- LineSearchGP starts");
569
570 double alpha = *t_opt;
571
572 double inf = 1e20;
573
574 double alphamin = 0.0;
575 double alphamax = inf;
576
577 double m1 = 0.1, m2 = 0.9;
578
579
580 /* // store the value fo the function */
581 /* double F[3]={0.,0.,0.}; */
582
583 /* // Store the (sub)-gradient of the function */
584 /* double A[9]={0.,0.,0.,0.,0.,0.,0.,0.,0.}; */
585 /* double B[9]={0.,0.,0.,0.,0.,0.,0.,0.,0.}; */
586
587 /* double velocity[3]={0.,0.,0.}; */
588
589 double mu = localproblem->mu[0];
590 double * qLocal = localproblem->q;
591 double * MLocal = localproblem->M->matrix0;
592
593 /* for (int i=0; i<3; i++) velocity[i] = MLocal[i+0*3]*R[0] + qLocal[i] */
594 /* + MLocal[i+1*3]*R[1] + */
595 /* + MLocal[i+2*3]*R[2] ; */
596
597 /* Function(R,velocity,mu,rho,F,A,B); */
598
599
600 // Computation of q(t) and q'(t) for t =0
601
602 double q0 = 0.5 * cblas_ddot(3, F, 1, F, 1);
603
604 double tmp[3] = {0., 0., 0.};
605
606 // Value of AW+B
607 double AWplusB[9] = {0., 0., 0., 0., 0., 0., 0., 0., 0.};
608
609 // compute A MLocal +B
610 for(int i = 0; i < 3; i++)
611 {
612 for(int j = 0; j < 3; j++)
613 {
614 AWplusB[i + 3 * j] = 0.0;
615 for(int k = 0; k < 3; k++)
616 {
617 AWplusB[i + 3 * j] += A[i + 3 * k] * MLocal[k + j * 3];
618 }
619 AWplusB[i + 3 * j] += B[i + 3 * j];
620 }
621 }
622
623 #ifdef DEBUG_MESSAGES
624 for(int l = 0; l < 3; l++)
625 {
626 for(int k = 0; k < 3; k++)
627 {
628 printf("AWplusB[%i+3*%i] = %le\t", l, k, AWplusB[l + 3 * k]);
629 }
630 printf("\n");
631 }
632 #endif
633
634 for(int i = 0; i < 3; i++)
635 {
636 tmp[i] = 0.0;
637 for(int j = 0; j < 3; j++)
638 {
639 tmp[i] += AWplusB[i + 3 * j] * dR[j] ;
640 }
641 }
642
643
644
645
646 double dqdt0 = 0.0;
647 for(int i = 0; i < 3; i++)
648 {
649 dqdt0 += F[i] * tmp[i];
650 }
651 #ifdef DEBUG_MESSAGES
652 printf("q0 = %12.8e \n", q0);
653 printf("dqdt0 = %12.8e \n", dqdt0);
654 for(int i = 0; i < 3; i++)
655 {
656 printf("tmp[%i] = %12.8e \t", i, tmp[i]);
657 }
658 printf("\n");
659 for(int i = 0; i < 3; i++)
660 {
661 printf("dR[%i] = %12.8e \t", i, dR[i]);
662 }
663 printf("\n");
664 #endif
665
666 for(int iter = 0; iter < LSitermax; iter++)
667 {
668
669 for(int i = 0; i < 3; i++) tmp[i] = R[i] + alpha * dR[i];
670
671 for(int i = 0; i < 3; i++) velocity[i] = MLocal[i + 0 * 3] * tmp[0] + qLocal[i]
672 + MLocal[i + 1 * 3] * tmp[1] +
673 + MLocal[i + 2 * 3] * tmp[2] ;
674
675 Function(tmp, velocity, mu, rho, F, NULL, NULL);
676
677 double q = 0.5 * cblas_ddot(3, F, 1, F, 1);
678
679 double slope = (q - q0) / alpha;
680
681 #ifdef DEBUG_MESSAGES
682 printf("q = %12.8e \n", q);
683 printf("slope = %12.8e \n", slope);
684 #endif
685
686
687 int C1 = (slope >= m2 * dqdt0);
688 int C2 = (slope <= m1 * dqdt0);
689
690 if(C1 && C2)
691 {
692 DEBUG_PRINTF("Success in LS: alpha = %12.8e\n", alpha);
693 *t_opt = alpha;
694
695 numerics_printf_verbose(2, "-- LineSearchGP success number of iteration = %i alpha = %.10e", iter, alpha);
696 return 0;
697
698 }
699 else if(!C1)
700 {
701 #ifdef DEBUG_MESSAGES
702 printf("LS: alpha too small = %12.8e\t, slope =%12.8e\n", alpha, slope);
703 printf(" m1*dqdt0 =%12.8e\t, m2*dqdt0 =%12.8e\n ", m1 * dqdt0, m2 * dqdt0);
704 #endif
705 //std::cout << "t = " << t << " is too small : slope = " << slope << ", m2*qp0 = " << m2*qp0 << std::endl;
706 alphamin = alpha;
707 }
708 else // not(C2)
709 {
710 #ifdef DEBUG_MESSAGES
711 printf("LS: alpha too big = %12.8e\t, slope =%12.8e\n", alpha, slope);
712 printf(" m1*dqdt0 =%12.8e\t, m2*dqdt0 =%12.8e\n ", m1 * dqdt0, m2 * dqdt0);
713 #endif
714 //std::cout << "t = " << t << " is too big : slope = " << slope << ", m1*qp0 = " << m1*qp0 << std::endl;
715 alphamax = alpha;
716 }
717 if(alpha < inf)
718 {
719 alpha = 0.5 * (alphamin + alphamax);
720 }
721 else
722 {
723 alpha = 10 * alpha;
724 }
725
726
727 }
728 numerics_printf_verbose(2,"-- LineSearchGP not succeed max number of iteration reached = %i with alpha = %.10e", LSitermax, alpha);
729 *t_opt = alpha;
730 return -1;
731 }
732 #define DEBUG_MESSAGES
fc3d_onecontact_nonsmooth_Newton_solvers_solve_damped(FrictionContactProblem * localproblem,double * R,SolverOptions * options)733 int fc3d_onecontact_nonsmooth_Newton_solvers_solve_damped(FrictionContactProblem* localproblem,
734 double * R, SolverOptions * options)
735 {
736
737 int * iparam = options->iparam;
738 double * dparam = options->dparam;
739
740 numerics_printf_verbose(2, "--------------- fc3d_onecontact_nonsmooth_Newton_solvers_solve_damped starts");
741
742 double mu = localproblem->mu[0];
743 double * qLocal = localproblem->q;
744
745 double norm_qLocal = sqrt(qLocal[0]*qLocal[0]+qLocal[1]*qLocal[1]+qLocal[2]*qLocal[2]);
746 double norm_relative = 1.0;
747 if(norm_qLocal> DBL_EPSILON)
748 {
749 norm_relative /= norm_qLocal;
750 }
751
752 double * MLocal = localproblem->M->matrix0;
753
754 /* store the increment */
755 double dR[3] = {0., 0., 0.};
756 /* store the value fo the function */
757 double F[3] = {0., 0., 0.};
758 /* Store the (sub)-gradient of the function */
759 double A[9] = {0., 0., 0., 0., 0., 0., 0., 0., 0.};
760 double B[9] = {0., 0., 0., 0., 0., 0., 0., 0., 0.};
761 /* Value of AW+B */
762 double AWplusB[9] = {0., 0., 0., 0., 0., 0., 0., 0., 0.};
763
764 /* retrieve value of rho */
765 double * rho;
766 if(options->solverId == SICONOS_FRICTION_3D_ONECONTACT_NSN ||
767 options->solverId == SICONOS_FRICTION_3D_ONECONTACT_NSN_GP)
768 {
769 rho = &options->dWork[3*iparam[SICONOS_FRICTION_3D_CURRENT_CONTACT_NUMBER]];
770 }
771 else
772 {
773 int nc = options->dWorkSize/4;
774 rho = &options->dWork[3*iparam[SICONOS_FRICTION_3D_CURRENT_CONTACT_NUMBER]+nc];
775 }
776
777 /* compute the velocity */
778 double velocity[3] = {0., 0., 0.};
779 cpy3(qLocal,velocity);
780 mvp3x3(MLocal,R,velocity);
781
782
783 int itermax = iparam[SICONOS_IPARAM_MAX_ITER];
784
785 /* Newton iteration */
786 int inew;
787 int info_solv3x3;
788 double t = 1.;
789 double t_opt = 1.;
790 double t_init = 1.;
791
792 int LSitermax = iparam[SICONOS_FRICTION_3D_NSN_LINESEARCH_MAX_ITER];
793 for(inew = 0 ; inew < itermax ; ++inew)
794 {
795 /* Update function and gradient */
796 Function(R, velocity, mu, rho, F, A, B);
797
798 /* compute -(A MLocal +B) */
799 mm3x3(A, MLocal, AWplusB);
800 add3x3(B, AWplusB);
801 scal3x3(-1., AWplusB);
802
803 #ifdef DEBUG_CHECK
804 fc3d_onecontact_nonsmooth_Newton_AC_debug(R, velocity, mu, rho, MLocal,
805 F, A, B, AWplusB, iparam);
806 #endif
807
808 /* Solve the linear system */
809 cpy3(F,dR);
810 info_solv3x3 = solve_3x3_gepp(AWplusB, dR);
811
812 /* if determinant is zero, replace dR=NaN with zero (i.e. don't modify R) and return early */
813 if(info_solv3x3)
814 {
815 dR[0] = 0;
816 dR[1] = 0;
817 dR[2] = 0;
818 }
819 if(!info_solv3x3)
820 {
821 /* Perform Line Search */
822
823 t_opt = t_init;
824 LineSearchGP(localproblem, Function, &t_opt, R, dR,
825 rho, LSitermax, F, A, B, velocity);
826 t = t_opt;
827 }
828
829 /* update iterates */
830 R[0] = R[0] + t * dR[0];
831 R[1] = R[1] + t * dR[1];
832 R[2] = R[2] + t * dR[2];
833
834 /* compute new residue */
835 cpy3(qLocal,velocity);
836 mvp3x3(MLocal,R,velocity);
837 Function(R, velocity, mu, rho, F, NULL, NULL);
838 dparam[SICONOS_DPARAM_RESIDU] = 0.5 * (F[0] * F[0] + F[1] * F[1] + F[2] * F[2])* norm_relative ; // improve with relative tolerance
839
840 if(info_solv3x3)
841 {
842 if(verbose > 0)
843 numerics_warning("fc3d_onecontact_nonsmooth_Newton_solvers_solve_damped",
844 "contact %i 3x3 linear system is irregular # iteration = %"
845 " error = %.10e \n",
846 iparam[SICONOS_FRICTION_3D_CURRENT_CONTACT_NUMBER],
847 inew, dparam[SICONOS_DPARAM_RESIDU]);
848 break;
849 }
850 numerics_printf_verbose(2,"-- contact %i # iteration = %i error = %.10e", iparam[SICONOS_FRICTION_3D_CURRENT_CONTACT_NUMBER], inew, dparam[SICONOS_DPARAM_RESIDU]);
851
852 if(dparam[SICONOS_DPARAM_RESIDU] < dparam[SICONOS_DPARAM_TOL])
853 {
854 iparam[SICONOS_IPARAM_ITER_DONE]=inew;
855 return 0;
856 }
857 }// End of the Newton iteration
858
859 iparam[SICONOS_IPARAM_ITER_DONE]=inew;
860 return 1;
861 }
862
compute_local_error(FrictionContactProblem * localproblem,double * local_reaction)863 static double compute_local_error(FrictionContactProblem* localproblem, double * local_reaction)
864 {
865 double * local_q = localproblem->q;
866 double norm_q = cblas_dnrm2(3, localproblem->q, 1);
867 double * local_M = localproblem->M->matrix0;
868 double worktmp[3];
869 double local_velocity[3] = {0., 0., 0.};
870 cpy3(local_q,local_velocity);
871 mvp3x3(local_M,local_reaction,local_velocity);
872 DEBUG_EXPR(NV_display(local_q,3););
873 DEBUG_EXPR(NV_display(local_velocity,3););
874 DEBUG_EXPR(NV_display(local_reaction,3););
875
876
877
878 double current_error = 0.0;
879
880 fc3d_unitary_compute_and_add_error(local_reaction,local_velocity, localproblem->mu[0], ¤t_error, worktmp);
881 current_error = sqrt(current_error);
882 DEBUG_PRINTF("absolute local error = %e", current_error);
883 if(fabs(norm_q) > DBL_EPSILON)
884 current_error /= norm_q;
885 DEBUG_PRINTF("relative local error = %e", current_error);
886 return current_error;
887 }
888
keep_or_discard_solution(FrictionContactProblem * localproblem,double * local_reaction,double * local_reaction_backup,SolverOptions * options,double * current_error)889 static void keep_or_discard_solution(FrictionContactProblem* localproblem, double * local_reaction,
890 double * local_reaction_backup, SolverOptions * options,
891 double * current_error)
892 {
893
894 double error=0.0;
895 /* error = options->dparam[SICONOS_DPARAM_RESIDU]; */
896
897 error = compute_local_error(localproblem, local_reaction);
898 DEBUG_PRINTF("New local error = %e\n", error);
899 int nan = isnan(options->dparam[SICONOS_DPARAM_RESIDU]) || isinf(options->dparam[SICONOS_DPARAM_RESIDU]);
900 if(nan)
901 {
902 DEBUG_PRINT("Residu is equal to nan or inf\n");
903 DEBUG_EXPR(NM_display(localproblem->M));
904 DEBUG_EXPR(NM_vector_display(localproblem->q,3));
905 DEBUG_EXPR(NM_vector_display(local_reaction,3));
906 DEBUG_EXPR(NM_vector_display(local_reaction_backup,3));
907
908 /* DEBUG_PRINTF("No hope for contact %d, setting to zero.\n", */
909 /* options->iparam[SICONOS_FRICTION_3D_CURRENT_CONTACT_NUMBER]); */
910 /* local_reaction[0] = 0; */
911 /* local_reaction[1] = 0; */
912 /* local_reaction[2] = 0; */
913 DEBUG_PRINTF("Discard the new local solution with error = %e\n", error);
914 DEBUG_PRINTF("Get back to the local backup solution = %e\n", *current_error);
915 cpy3(local_reaction_backup, local_reaction);
916 //memcpy(local_reaction, local_reaction_backup, sizeof(double)*3);
917 }
918 else
919 {
920 if(error <= options->dparam[SICONOS_DPARAM_TOL] || error <= *current_error)
921 {
922 DEBUG_PRINTF("Keep the new local solution with error = %e\n", error);
923 *current_error = error;
924 cpy3(local_reaction, local_reaction_backup);
925 //memcpy(local_reaction_backup, local_reaction, sizeof(double)*3);
926 }
927 else
928 {
929 DEBUG_PRINTF("Discard the new local solution with error = %e\n", error);
930 DEBUG_PRINTF("Get back to the local backup solution = %e\n", *current_error);
931 cpy3(local_reaction_backup, local_reaction);
932 //memcpy(local_reaction, local_reaction_backup, sizeof(double)*3);
933 }
934 }
935 }
936
fc3d_onecontact_nonsmooth_Newton_solvers_solve_hybrid(FrictionContactProblem * localproblem,double * local_reaction,SolverOptions * options)937 int fc3d_onecontact_nonsmooth_Newton_solvers_solve_hybrid(FrictionContactProblem* localproblem, double * local_reaction, SolverOptions * options)
938 {
939 numerics_printf_verbose(2, "--------------- fc3d_onecontact_nonsmooth_Newton_solvers_solve_hybrid starts");
940
941 int info = -1;
942 double local_reaction_backup[3] = {local_reaction[0], local_reaction[1], local_reaction[2]};
943
944 int max_loop = options->iparam[SICONOS_FRICTION_3D_NSN_HYBRID_MAX_LOOP];
945
946 int newton_iteration_number = options->iparam[SICONOS_IPARAM_MAX_ITER];
947 int pli_iteration_number = options->iparam[SICONOS_FRICTION_3D_NSN_HYBRID_MAX_ITER];
948
949 /* compute current_error */
950 double current_error = compute_local_error(localproblem, local_reaction);
951
952 DEBUG_PRINTF("fc3d_onecontact_nonsmooth_Newton_solvers_solve_hybrid current_error= %e\n", current_error);
953
954 if(!(options->iparam[SICONOS_FRICTION_3D_NSN_HYBRID_STRATEGY] == SICONOS_FRICTION_3D_NSN_HYBRID_STRATEGY_PLI_NSN_LOOP ||
955 options->iparam[SICONOS_FRICTION_3D_NSN_HYBRID_STRATEGY] == SICONOS_FRICTION_3D_NSN_HYBRID_STRATEGY_NSN_AND_PLI_NSN_LOOP))
956 {
957 numerics_error("fc3d_onecontact_nonsmooth_Newton_solvers_solve_hybrid","Unknown local nsn hybrid solver");
958 }
959
960 /* 0 - Perform a first call to NSN solver to see if it succeeds quickly */
961
962 if(options->iparam[SICONOS_FRICTION_3D_NSN_HYBRID_STRATEGY ] == SICONOS_FRICTION_3D_NSN_HYBRID_STRATEGY_NSN_AND_PLI_NSN_LOOP)
963 {
964
965 options->iparam[SICONOS_IPARAM_MAX_ITER]= newton_iteration_number;
966 info = fc3d_onecontact_nonsmooth_Newton_solvers_solve_damped(localproblem, local_reaction, options);
967
968 DEBUG_PRINTF("NSN solver ended with residual = %e\n", options->dparam[SICONOS_DPARAM_RESIDU]);
969
970 keep_or_discard_solution(localproblem, local_reaction,
971 local_reaction_backup, options,
972 & current_error);
973
974
975 if(current_error <= options->dparam[SICONOS_DPARAM_TOL])
976 {
977 options->dparam[SICONOS_DPARAM_RESIDU] = current_error;
978 DEBUG_PRINTF("First call of NSN solver ends with current_error= %e\n", current_error);
979 DEBUG_END("fc3d_onecontact_nonsmooth_Newton_solvers_solve_hybrid\n");
980 return info;
981 }
982 }
983
984 if(options->iparam[SICONOS_FRICTION_3D_NSN_HYBRID_STRATEGY] == SICONOS_FRICTION_3D_NSN_HYBRID_STRATEGY_PLI_NSN_LOOP ||
985 options->iparam[SICONOS_FRICTION_3D_NSN_HYBRID_STRATEGY] == SICONOS_FRICTION_3D_NSN_HYBRID_STRATEGY_NSN_AND_PLI_NSN_LOOP)
986 {
987 int loop = 0 ;
988 while(loop < max_loop && current_error >= options->dparam[SICONOS_DPARAM_TOL])
989 {
990 loop++;
991 DEBUG_PRINTF("SICONOS_FRICTION_3D_ONECONTACT_NSN_GP_HYBRID: loop = %i\n", loop);
992
993 /* 1 - fixed point projection solver */
994
995 options->iparam[SICONOS_IPARAM_MAX_ITER]= pli_iteration_number;
996 info = fc3d_projectionOnConeWithLocalIteration_solve(localproblem, local_reaction, options);
997 DEBUG_PRINTF("PLI solver ended with residual = %e\n", options->dparam[SICONOS_DPARAM_RESIDU]);
998
999 keep_or_discard_solution(localproblem, local_reaction,
1000 local_reaction_backup, options,
1001 & current_error);
1002 if(current_error < options->dparam[SICONOS_DPARAM_TOL])
1003 {
1004 options->iparam[SICONOS_IPARAM_MAX_ITER]= newton_iteration_number;
1005 break;
1006 }
1007
1008 /* 2 - nonsmooth Newton solver */
1009
1010 options->iparam[SICONOS_IPARAM_MAX_ITER]= newton_iteration_number;
1011 info = fc3d_onecontact_nonsmooth_Newton_solvers_solve_damped(localproblem, local_reaction, options);
1012 DEBUG_PRINTF("NSN solver ended with residual = %e\n", options->dparam[SICONOS_DPARAM_RESIDU]);
1013
1014 keep_or_discard_solution(localproblem, local_reaction,
1015 local_reaction_backup, options,
1016 & current_error);
1017 }
1018
1019 if(loop == max_loop && max_loop != 1)
1020 {
1021 if(verbose >0)
1022 {
1023 printf("Maximum number of loop (%i) in SICONOS_FRICTION_3D_ONECONTACT_NSN_GP_HYBRID has been reached for contact %i with error = %e \n",max_loop,options->iparam[SICONOS_FRICTION_3D_CURRENT_CONTACT_NUMBER],current_error);
1024 }
1025 }
1026 }
1027
1028 options->dparam[SICONOS_DPARAM_RESIDU] = current_error;
1029 DEBUG_PRINTF("fc3d_onecontact_nonsmooth_Newton_solvers_solve_hybrid END current_error= %e\n", current_error);
1030 DEBUG_END("fc3d_onecontact_nonsmooth_Newton_solvers_solve_hybrid\n")
1031 return info;
1032 }
1033
1034
1035
fc3d_onecontact_nsn_set_default(SolverOptions * options)1036 void fc3d_onecontact_nsn_set_default(SolverOptions* options)
1037 {
1038 /* Value of rho parameter */
1039 options->iparam[SICONOS_FRICTION_3D_NSN_RHO_STRATEGY] = SICONOS_FRICTION_3D_NSN_FORMULATION_RHO_STRATEGY_SPLIT_SPECTRAL_NORM_COND;
1040 options->dparam[SICONOS_FRICTION_3D_NSN_RHO] =1.0;
1041
1042 /* Choice of formulation */
1043 options->iparam[SICONOS_FRICTION_3D_NSN_FORMULATION] = SICONOS_FRICTION_3D_NSN_FORMULATION_ALARTCURNIER_STD ;
1044 /* Choice of line -search method */
1045 options->iparam[SICONOS_FRICTION_3D_NSN_LINESEARCH] = SICONOS_FRICTION_3D_NSN_LINESEARCH_NO;
1046
1047 /* parameters for hybrid solvers */
1048 options->iparam[SICONOS_FRICTION_3D_NSN_HYBRID_STRATEGY] = SICONOS_FRICTION_3D_NSN_HYBRID_STRATEGY_PLI_NSN_LOOP;
1049
1050 options->iparam[SICONOS_FRICTION_3D_NSN_HYBRID_MAX_LOOP] = 1;
1051 options->iparam[SICONOS_FRICTION_3D_NSN_HYBRID_MAX_ITER] = 10;
1052 }
1053
fc3d_onecontact_nsn_gp_set_default(SolverOptions * options)1054 void fc3d_onecontact_nsn_gp_set_default(SolverOptions* options)
1055 {
1056 /* Value of rho parameter */
1057 options->iparam[SICONOS_FRICTION_3D_NSN_RHO_STRATEGY] = SICONOS_FRICTION_3D_NSN_FORMULATION_RHO_STRATEGY_SPLIT_SPECTRAL_NORM;
1058 options->dparam[SICONOS_FRICTION_3D_NSN_RHO] =1.0;
1059
1060 /* Choice of formulation */
1061 options->iparam[SICONOS_FRICTION_3D_NSN_FORMULATION] = SICONOS_FRICTION_3D_NSN_FORMULATION_ALARTCURNIER_STD ;
1062 /* Choice of line -search method */
1063 options->iparam[SICONOS_FRICTION_3D_NSN_LINESEARCH] = SICONOS_FRICTION_3D_NSN_LINESEARCH_GOLDSTEINPRICE;
1064 options->iparam[SICONOS_FRICTION_3D_NSN_LINESEARCH_MAX_ITER] = 10;
1065
1066 /* parameters for hybrid solvers */
1067 options->iparam[SICONOS_FRICTION_3D_NSN_HYBRID_STRATEGY] = SICONOS_FRICTION_3D_NSN_HYBRID_STRATEGY_NSN_AND_PLI_NSN_LOOP;
1068 options->iparam[SICONOS_FRICTION_3D_NSN_HYBRID_MAX_LOOP] = 1;
1069 options->iparam[SICONOS_FRICTION_3D_NSN_HYBRID_MAX_ITER] = 100;
1070 }
1071