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 <stdio.h> // for NULL, fprintf, stderr
20 #include <stdlib.h> // for exit, EXIT_FAILURE
21 #include "Friction_cst.h" // for SICONOS_GLOBAL_FRICTION_3D...
22 #include "GlobalRollingFrictionContactProblem.h" // for GlobalFrictionContactProblem
23 #include "NonSmoothDrivers.h" // for gfc3d_driver
24 #include "NumericsFwd.h" // for SolverOptions, GlobalFrict...
25 #include "SolverOptions.h" // for SolverOptions, solver_opti...
26 #include "NumericsMatrix.h" //
27
28 #include "global_rolling_fc_Solvers.h" // for gfc3d_ACLMFixedPoint, gfc3...
29 #include "numerics_verbose.h" // for numerics_printf_verbose
30 //#include "gfc3d_compute_error.h"
31 //#include "SiconosBlas.h" // for cblas_dcopy, cblas_dscal
32
33 #include "siconos_debug.h" // for DEBUG_EXPR
34 #ifdef DEBUG_MESSAGES
35 #include "NumericsVector.h"
36 #include "NumericsMatrix.h"
37 #endif
38
39 const char* const SICONOS_GLOBAL_ROLLING_FRICTION_3D_NSGS_WR_STR = "GFC3D_NSGS_WR";
40 #ifdef WITH_FCLIB
41 #include "string.h" // for strcpy, strcat
42 #include "fclib_interface.h" // for frictionContact_fclib_write, fri...
43 #endif
44
45 //#define FCLIB_OUTPUT
46
47 #ifdef FCLIB_OUTPUT
48 #ifdef WITH_FCLIB
49 #include "string.h" // for strcpy, strcat
50 #include "fclib_interface.h" // for frictionContact_fclib_write, fri...
51 #endif
52 static int fccounter = -1;
53 #endif
54
g_rolling_fc3d_driver(GlobalRollingFrictionContactProblem * problem,double * reaction,double * velocity,double * globalVelocity,SolverOptions * options)55 int g_rolling_fc3d_driver(GlobalRollingFrictionContactProblem* problem, double *reaction, double *velocity,
56 double* globalVelocity, SolverOptions* options)
57 {
58 #ifdef FCLIB_OUTPUT
59 #ifdef WITH_FCLIB
60 fccounter ++;
61 int freq_output=1;
62 int nc = problem->numberOfContacts;
63 if(nc >0)
64 {
65 if(fccounter % freq_output == 0)
66 {
67 char fname[256];
68 sprintf(fname, "GRFC3D-%.5d-%.5d.hdf5", (int)nc, fccounter);
69 printf("Dump GRFC3D-%.5d-%.5d.hdf5.\n", (int)nc, fccounter);
70 /* printf("ndof = %i.\n", ndof); */
71 /* FILE * foutput = fopen(fname, "w"); */
72 int n = 100;
73 char * title = (char *)malloc(n * sizeof(char));
74 strcpy(title, "GRFC3 dump in hdf5");
75 char * description = (char *)malloc(n * sizeof(char));
76 strcpy(description, "Rewriting in hdf5 through siconos of ");
77 strcat(description, fname);
78 strcat(description, " in FCLIB format");
79 char * mathInfo = (char *)malloc(n * sizeof(char));
80 strcpy(mathInfo, "unknown");
81 globalRollingFrictionContact_fclib_write(problem,
82 title,
83 description,
84 mathInfo,
85 fname);
86 }
87 /* fclose(foutput); */
88 }
89 #else
90 printf("Fclib is not available ...\n");
91 #endif
92 #endif
93
94
95
96
97 assert(options->isSet);
98 DEBUG_EXPR(NV_display(globalVelocity,problem_ori->M->size0););
99 if(verbose > 0)
100 solver_options_print(options);
101
102 /* Solver name */
103 /* const char* const name = options->solverName;*/
104
105 int info = -1 ;
106
107 if(problem->dimension != 5)
108 numerics_error("grfc3d_driver", "Dimension of the problem : problem-> dimension is not compatible or is not set");
109
110 /* if there is no contact, we compute directly the global velocity as M^{-1}q */
111 int m = problem->H->size1;
112 if(m ==0)
113 {
114 numerics_printf_verbose(1,"---- GFC3D - DRIVER . No contact case. Direct computation of global velocity");
115 globalRollingFrictionContact_computeGlobalVelocity(problem, reaction, globalVelocity);
116 return 0;
117 }
118
119 /* Non Smooth Gauss Seidel (NSGS) with reduction*/
120 switch(options->solverId)
121 {
122 case SICONOS_GLOBAL_ROLLING_FRICTION_3D_NSGS_WR:
123 {
124
125 numerics_printf_verbose(1," ========================== Call NSGS_WR solver with reformulation into Rolling Friction-Contact 3D problem ==========================\n");
126 grfc3d_nsgs_wr(problem, reaction, velocity, globalVelocity, &info, options);
127 break;
128 }
129 /* case SICONOS_GLOBAL_ROLLING_FRICTION_3D_IPM: */
130 /* { */
131 /* grfc3d_IPM(problem, reaction, velocity, */
132 /* globalVelocity, &info, options); */
133 /* break; */
134
135 /* } */
136 default:
137 {
138 fprintf(stderr, "Numerics, grfc3d_driver failed. Unknown solver %d.\n", options->solverId);
139 exit(EXIT_FAILURE);
140
141 }
142 }
143
144 return info;
145
146 }
147
grfc3d_checkTrivialCaseGlobal(int n,double * q,double * velocity,double * reaction,double * globalVelocity,SolverOptions * options)148 int grfc3d_checkTrivialCaseGlobal(int n, double* q, double* velocity, double* reaction, double * globalVelocity, SolverOptions* options)
149 {
150 /* norm of vector q */
151 /* double qs = cblas_dnrm2( n , q , 1 ); */
152 /* int i; */
153 int info = -1;
154 /* if( qs <= DBL_EPSILON ) */
155 /* { */
156 /* // q norm equal to zero (less than DBL_EPSILON) */
157 /* // -> trivial solution: reaction = 0 and velocity = q */
158 /* for( i = 0 ; i < n ; ++i ) */
159 /* { */
160 /* velocity[i] = q[i]; */
161 /* reaction[i] = 0.; */
162 /* } */
163 /* iparam[2] = 0; */
164 /* iparam[4]= 0; */
165 /* dparam[1] = 0.0; */
166 /* dparam[3] = 0.0; */
167 /* info = 0; */
168 /* if(iparam[1]>0) */
169 /* printf("fc3d driver, norm(q) = 0, trivial solution reaction = 0, velocity = q.\n"); */
170 /* } */
171 return info;
172 }
173