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 <assert.h> // for assert
19 #include <float.h> // for DBL_EPSILON
20 #include <stdio.h> // for NULL
21 #include <string.h> // for NULL
22 #include "FrictionContactProblem.h" // for FrictionContac...
23 #include "Friction_cst.h" // for SICONOS_FRICTI...
24 #include "NonSmoothDrivers.h" // for fc3d_driver
25 #include "NumericsFwd.h" // for SolverOptions
26 #include "SolverOptions.h" // for SolverOptions
27 #include "fc3d_Solvers.h" // for fc3d_ACLMFixed...
28 #include "fc3d_nonsmooth_Newton_AlartCurnier.h" // for fc3d_nonsmooth...
29 #include "fc3d_nonsmooth_Newton_FischerBurmeister.h" // for fc3d_nonsmooth...
30 #include "fc3d_nonsmooth_Newton_natural_map.h" // for fc3d_nonsmooth...
31 #include "fc3d_onecontact_nonsmooth_Newton_solvers.h" // for fc3d_onecontac...
32 #include "fc3d_projection.h" // for fc3d_projectio...
33 #include "fc3d_unitary_enumerative.h" // for fc3d_unitary_e...
34 #include "numerics_verbose.h" // for numerics_printf
35
36 const char* const SICONOS_FRICTION_3D_NSGS_STR = "FC3D_NSGS";
37 const char* const SICONOS_FRICTION_3D_NSGSV_STR = "FC3D_NSGSV";
38 const char* const SICONOS_FRICTION_3D_TFP_STR = "FC3D_TFP";
39 const char* const SICONOS_FRICTION_3D_PFP_STR = "FC3D_PFP";
40 const char* const SICONOS_FRICTION_3D_NSN_AC_STR = "FC3D_NSN_AC";
41 const char* const SICONOS_FRICTION_3D_NSN_AC_TEST_STR = "FC3D_NSN_AC_TEST";
42 const char* const SICONOS_FRICTION_3D_NSN_FB_STR = "FC3D_NSN_FB";
43 const char* const SICONOS_FRICTION_3D_NSN_NM_STR = "FC3D_NSN_NM";
44 const char* const SICONOS_FRICTION_3D_DSFP_STR = "FC3D_DeSaxceFixedPoint";
45 const char* const SICONOS_FRICTION_3D_NCPGlockerFBFixedPoint_STR = "FC3D_NCPGlockerFBFixedPoint";
46 const char* const SICONOS_FRICTION_3D_ONECONTACT_NSN_STR = "FC3D_ONECONTACT_NSN";
47 const char* const SICONOS_FRICTION_3D_ONECONTACT_NSN_GP_STR = "FC3D_ONECONTACT_NSN_GP";
48 const char* const SICONOS_FRICTION_3D_ONECONTACT_NSN_GP_HYBRID_STR = "FC3D_ONECONTACT_NSN_GP_HYBRID";
49 const char* const SICONOS_FRICTION_3D_NCPGlockerFBNewton_STR = "FC3D_NCPGlockerFBNewton";
50 const char* const SICONOS_FRICTION_3D_ONECONTACT_ProjectionOnConeWithDiagonalization_STR = "FC3D_ProjectionOnConeWithDiagonalization";
51 const char* const SICONOS_FRICTION_3D_ONECONTACT_ProjectionOnCone_STR = "FC3D_ProjectionOnCone";
52 const char* const SICONOS_FRICTION_3D_ONECONTACT_ProjectionOnConeWithLocalIteration_STR = "FC3D_ProjectionOnConeWithLocalIteration";
53 const char* const SICONOS_FRICTION_3D_ONECONTACT_ProjectionOnConeWithRegularization_STR = "FC3D_projectionOnConeWithRegularization";
54 const char* const SICONOS_FRICTION_3D_NCPGlockerFBPATH_STR = "FC3D_NCPGlockerFBPATH";
55 const char* const SICONOS_FRICTION_3D_ONECONTACT_ProjectionOnCylinder_STR = "FC3D_projectionOnCylinder";
56 const char* const SICONOS_FRICTION_3D_ONECONTACT_ProjectionOnCylinderWithLocalIteration_STR = "FC3D_projectionOnCylinderWithLocalIteration";
57 const char* const SICONOS_FRICTION_3D_ONECONTACT_ProjectionOnCone_velocity_STR = "FC3D_ProjectionOnCone_velocity";
58 const char* const SICONOS_FRICTION_3D_CONVEXQP_PG_CYLINDER_STR = "FC3D ConvexQP PG solver";
59 const char* const SICONOS_FRICTION_3D_VI_FPP_Cylinder_STR = "FC3D_VI_FPP_Cylinder";
60
61
62 const char* const SICONOS_FRICTION_3D_EG_STR = "FC3D_ExtraGradient";
63 const char* const SICONOS_FRICTION_3D_FPP_STR = "FC3D_FixedPointProjection";
64 const char* const SICONOS_FRICTION_3D_VI_EG_STR = "FC3D_VI_ExtraGradient";
65 const char* const SICONOS_FRICTION_3D_VI_FPP_STR = "FC3D_VI_FixedPointProjection";
66 const char* const SICONOS_FRICTION_3D_HP_STR = "FC3D_HyperplaneProjection";
67 const char* const SICONOS_FRICTION_3D_PROX_STR = "FC3D_PROX";
68 const char* const SICONOS_FRICTION_3D_GAMS_PATH_STR = "FC3D_GAMS_PATH";
69 const char* const SICONOS_FRICTION_3D_GAMS_PATHVI_STR = "FC3D_GAMS_PATHVI";
70 const char* const SICONOS_FRICTION_3D_GAMS_LCP_PATH_STR = "FC3D_GAMS_LCP_PATH";
71 const char* const SICONOS_FRICTION_3D_GAMS_LCP_PATHVI_STR = "FC3D_GAMS_LCP_PATHVI";
72 const char* const SICONOS_FRICTION_3D_ONECONTACT_QUARTIC_STR = "FC3D_QUARTIC";
73 const char* const SICONOS_FRICTION_3D_ONECONTACT_QUARTIC_NU_STR = "FC3D_QUARTIC_NU";
74 const char* const SICONOS_FRICTION_3D_ACLMFP_STR = "FC3D_ACLMFP";
75 const char* const SICONOS_FRICTION_3D_SOCLCP_STR = "FC3D_SOCLCP";
76
77
78
fc3d_driver(FrictionContactProblem * problem,double * reaction,double * velocity,SolverOptions * options)79 int fc3d_driver(FrictionContactProblem* problem,
80 double *reaction, double *velocity,
81 SolverOptions* options)
82 {
83 if(options == NULL)
84 numerics_error("fc3d_driver", "null input for solver options");
85
86 assert(options->isSet); /* true(1) if the SolverOptions structure has been filled in else false(0) */
87
88 if(verbose > 1)
89 solver_options_print(options);
90
91 int info = -1 ;
92
93 if(problem->dimension != 3)
94 numerics_error("fc3d_driver", "Dimension of the problem : problem-> dimension is not compatible or is not set");
95
96 /* Check for trivial case */
97 info = fc3d_checkTrivialCase(problem, velocity, reaction, options);
98 if(info == 0)
99 {
100 /* If a trivial solution is found, we set the number of iterations to 0
101 and the reached acuracy to 0.0 .
102 */
103 options->iparam[SICONOS_IPARAM_ITER_DONE] = 0;
104 options->dparam[SICONOS_DPARAM_RESIDU] = 0.0;
105 goto exit;
106 }
107
108
109 switch(options->solverId)
110 {
111 /* Non Smooth Gauss Seidel (NSGS) */
112 case SICONOS_FRICTION_3D_NSGS:
113 {
114 numerics_printf(" ========================== Call NSGS solver for Friction-Contact 3D problem ==========================\n");
115 fc3d_nsgs(problem, reaction, velocity, &info, options);
116 break;
117 }
118 case SICONOS_FRICTION_3D_NSGSV:
119 {
120 numerics_printf(" ========================== Call NSGSV solver for Friction-Contact 3D problem ==========================\n");
121 fc3d_nsgs_velocity(problem, reaction, velocity, &info, options);
122 break;
123 }
124 /* ADMM*/
125 case SICONOS_FRICTION_3D_ADMM:
126 {
127 numerics_printf(" ========================== Call ADMM solver for Friction-Contact 3D problem ==========================\n");
128 fc3d_admm(problem, reaction, velocity, &info, options);
129 break;
130 }
131 /* Proximal point algorithm */
132 case SICONOS_FRICTION_3D_PROX:
133 {
134 numerics_printf(" ========================== Call PROX (Proximal Point) solver for Friction-Contact 3D problem ==========================\n");
135 fc3d_proximal(problem, reaction, velocity, &info, options);
136 break;
137 }
138 /* Tresca Fixed point algorithm */
139 case SICONOS_FRICTION_3D_TFP:
140 {
141 numerics_printf(" ========================== Call TFP (Tresca Fixed Point) solver for Friction-Contact 3D problem ==========================\n");
142 fc3d_TrescaFixedPoint(problem, reaction, velocity, &info, options);
143 break;
144 }
145 /* Panagiotopoulos Fixed point algorithm */
146 case SICONOS_FRICTION_3D_PFP:
147 {
148 numerics_printf(" ========================== Call PFP (Panagiotopoulos Fixed Point) solver for Friction-Contact 3D problem ==========================\n");
149 fc3d_Panagiotopoulos_FixedPoint(problem, reaction, velocity, &info, options);
150 break;
151 }
152 /* ACLM Fixed point algorithm */
153 case SICONOS_FRICTION_3D_ACLMFP:
154 {
155 numerics_printf(" ========================== Call ACLM (Acary Cadoux Lemarechal Malick Fixed Point) solver for Friction-Contact 3D problem ==========================\n");
156 fc3d_ACLMFixedPoint(problem, reaction, velocity, &info, options);
157 break;
158 }
159 /* SOCLCP Fixed point algorithm */
160 case SICONOS_FRICTION_3D_SOCLCP:
161 {
162 numerics_printf(" ========================== Call SOCLCP solver for Friction-Contact 3D problem (Associated one) ==========================\n");
163 fc3d_SOCLCP(problem, reaction, velocity, &info, options);
164 break;
165 }
166 /* De Saxce Fixed point algorithm */
167 case SICONOS_FRICTION_3D_DSFP:
168 {
169 numerics_printf(" ========================== Call DeSaxce Fixed Point (DSFP) solver for Friction-Contact 3D problem ==========================\n");
170 fc3d_DeSaxceFixedPoint(problem, reaction, velocity, &info, options);
171 break;
172 }
173 /* Fixed point projection algorithm */
174 case SICONOS_FRICTION_3D_FPP:
175 {
176 numerics_printf(" ========================== Call Fixed Point Projection (FPP) solver for Friction-Contact 3D problem ==========================\n");
177 fc3d_fixedPointProjection(problem, reaction, velocity, &info, options);
178 break;
179 }
180
181 /* Extra Gradient algorithm */
182 case SICONOS_FRICTION_3D_EG:
183 {
184 numerics_printf(" ========================== Call ExtraGradient (EG) solver for Friction-Contact 3D problem ==========================\n");
185 fc3d_ExtraGradient(problem, reaction, velocity, &info, options);
186 break;
187 }
188 /* VI Fixed Point Projection algorithm */
189 case SICONOS_FRICTION_3D_VI_FPP:
190 {
191 numerics_printf(" ========================== Call VI_FixedPointProjection (VI_FPP) solver for Friction-Contact 3D problem ==========================\n");
192 fc3d_VI_FixedPointProjection(problem, reaction, velocity, &info, options);
193 break;
194 }
195 /* VI Extra Gradient algorithm */
196 case SICONOS_FRICTION_3D_VI_EG:
197 {
198 numerics_printf(" ========================== Call VI_ExtraGradient (VI_EG) solver for Friction-Contact 3D problem ==========================\n");
199 fc3d_VI_ExtraGradient(problem, reaction, velocity, &info, options);
200 break;
201 }
202 /* Hyperplane Projection algorithm */
203 case SICONOS_FRICTION_3D_HP:
204 {
205 numerics_printf(" ========================== Call Hyperplane Projection (HP) solver for Friction-Contact 3D problem ==========================\n");
206 fc3d_HyperplaneProjection(problem, reaction, velocity, &info, options);
207 break;
208 }
209 /* Alart Curnier in local coordinates */
210 case SICONOS_FRICTION_3D_NSN_AC:
211 {
212 numerics_printf(" ========================== Call Alart Curnier solver for Friction-Contact 3D problem ==========================\n");
213 fc3d_nonsmooth_Newton_AlartCurnier(problem, reaction, velocity, &info, options);
214 break;
215 }
216 /* XXX to delete */
217 case SICONOS_FRICTION_3D_NSN_AC_TEST:
218 {
219 numerics_printf(" ========================== Call Alart Curnier solver for Friction-Contact 3D problem ==========================\n");
220 fc3d_nonsmooth_Newton_AlartCurnier2(problem, reaction, velocity, &info, options);
221 break;
222 }
223 /* Fischer Burmeister in local coordinates */
224 case SICONOS_FRICTION_3D_NSN_FB:
225 {
226 numerics_printf(" ========================== Call Fischer Burmeister solver for Friction-Contact 3D problem ==========================\n");
227 fc3d_nonsmooth_Newton_FischerBurmeister(problem, reaction, velocity, &info, options);
228 break;
229 }
230 case SICONOS_FRICTION_3D_NSN_NM:
231 {
232 numerics_printf(" ========================== Call natural map solver for Friction-Contact 3D problem ==========================\n");
233 fc3d_nonsmooth_Newton_NaturalMap(problem, reaction, velocity, &info, options);
234 break;
235 }
236 case SICONOS_FRICTION_3D_ONECONTACT_QUARTIC_NU:
237 case SICONOS_FRICTION_3D_ONECONTACT_QUARTIC:
238 {
239 numerics_printf(" ========================== Call Quartic solver for Friction-Contact 3D problem ==========================\n");
240 fc3d_unitary_enumerative(problem, reaction, velocity, &info, options);
241 break;
242 }
243 case SICONOS_FRICTION_3D_ONECONTACT_NSN:
244 case SICONOS_FRICTION_3D_ONECONTACT_NSN_GP:
245 case SICONOS_FRICTION_3D_ONECONTACT_NSN_GP_HYBRID:
246 {
247 numerics_printf(" ========================== Call Newton-based solver for one contact Friction-Contact 3D problem ==========================\n");
248 fc3d_onecontact_nonsmooth_Newton_solvers_initialize(problem, problem, options);
249 info = fc3d_onecontact_nonsmooth_Newton_solvers_solve(problem, reaction, options);
250 break;
251 }
252 case SICONOS_FRICTION_3D_ONECONTACT_ProjectionOnConeWithLocalIteration:
253 {
254 numerics_printf(" ========================== Call Projection on cone solver for one contact Friction-Contact 3D problem ==========================\n");
255 fc3d_projectionOnConeWithLocalIteration_initialize(problem, problem, options);
256 info = fc3d_projectionOnConeWithLocalIteration_solve(problem, reaction, options);
257 fc3d_projectionOnConeWithLocalIteration_free(problem, problem, options);
258
259 break;
260 }
261 case SICONOS_FRICTION_3D_GAMS_PATH:
262 {
263 numerics_printf(" ========================== Call PATH solver via GAMS for an AVI Friction-Contact 3D problem ==========================\n");
264 fc3d_AVI_gams_path(problem, reaction, velocity, &info, options);
265 break;
266 }
267 case SICONOS_FRICTION_3D_GAMS_PATHVI:
268 {
269 numerics_printf(" ========================== Call PATHVI solver via GAMS for an AVI Friction-Contact 3D problem ==========================\n");
270 fc3d_AVI_gams_pathvi(problem, reaction, velocity, &info, options);
271 break;
272 }
273 case SICONOS_FRICTION_3D_GAMS_LCP_PATH:
274 {
275 numerics_printf(" ========================== Call PATH solver via GAMS for an LCP-based reformulation of the AVI Friction-Contact 3D problem ==========================\n");
276 fc3d_lcp_gams_path(problem, reaction, velocity, &info, options);
277 break;
278 }
279 case SICONOS_FRICTION_3D_GAMS_LCP_PATHVI:
280 {
281 numerics_printf(" ========================== Call PATHVI solver via GAMS for an LCP-based reformulation of the AVI Friction-Contact 3D problem ==========================\n");
282 fc3d_lcp_gams_pathvi(problem, reaction, velocity, &info, options);
283 break;
284 }
285 default:
286 {
287 char msg[200];
288 strcpy(msg, "Unknown solver : ");
289 strcat(msg, solver_options_id_to_name(options->solverId));
290 strcat(msg, "\n");
291 numerics_warning("fc3d_driver", msg);
292 numerics_error("fc3d_driver", msg);
293 info = 1;
294 }
295 }
296
297 exit:
298 return info;
299
300 }
301
fc3d_checkTrivialCase(FrictionContactProblem * problem,double * velocity,double * reaction,SolverOptions * options)302 int fc3d_checkTrivialCase(FrictionContactProblem* problem, double* velocity,
303 double* reaction, SolverOptions* options)
304 {
305 /* Number of contacts */
306 int nc = problem->numberOfContacts;
307 double* q = problem->q;
308 /* Dimension of the problem */
309 int n = 3 * nc;
310 int i = 0;
311 /*take off? R=0 ?*/
312 for(i = 0; i < nc; i++)
313 {
314 if(q[3 * i] < -DBL_EPSILON)
315 return -1;
316 }
317 for(i = 0 ; i < n ; ++i)
318 {
319 velocity[i] = q[i];
320 reaction[i] = 0.;
321 }
322
323 numerics_printf("fc3d fc3d_checkTrivialCase, take off, trivial solution reaction = 0, velocity = q.\n");
324 return 0;
325 }
326