1 /* CalculiX - A 3-dimensional finite element program */
2 /* Copyright (C) 1998-2021 Guido Dhondt */
3
4 /* This program is free software; you can redistribute it and/or */
5 /* modify it under the terms of the GNU General Public License as */
6 /* published by the Free Software Foundation(version 2); */
7 /* */
8
9 /* This program is distributed in the hope that it will be useful, */
10 /* but WITHOUT ANY WARRANTY; without even the implied warranty of */
11 /* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the */
12 /* GNU General Public License for more details. */
13
14 /* You should have received a copy of the GNU General Public License */
15 /* along with this program; if not, write to the Free Software */
16 /* Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */
17
18 #include <stdio.h>
19 #include <math.h>
20 #include <stdlib.h>
21 #include "CalculiX.h"
22 #include "spooles.h"
23
24 /**
25 * @brief *MASSLESS DYNAMIC CONTACT*: Main computation of step for dynamic massless contact
26 *
27 * calling also functions:
28 * + createcontactdofs():
29 * + expand_auw:
30 * + extract_matrices:
31 * + detectactivecont1
32 * + detectactivecont2
33 *
34 * @param au LOWER triangle of STIFFNESS matrix of size: neq[0]
35 * @param ad DIAGONAL of STIFFNESS matrix of size: neq[0]
36 * @param aub UPPER triangle of MASS matrix of size: neq[0]
37 * @param adb DIAGONAL of MASS matrix of size: neq[0]
38 * @param jq Location in field **irow** of the first subdiagonal nonzero in column i (only for symmetric matrices)
39 * @param irow Row of element i in field au (i.e. au(i))
40 * @param neq NTOT of equations:
41 * + neq[0]: mechanical TOTDOF,
42 * + neq[1]: TOTDOF_mech + TOTDOF_thermal
43 * + neq[2]: neq[1] + # of single point constraints (only for modal calculations)
44 * @param nzs projected nonzeros:
45 * + nzs[0]: sum of projected nonzero mechanical off-diagonal terms
46 * + nzs[1]: nzs[0]+sum of projected nonzero thermal off-diagonal terms
47 * + nzs[2]: nzs[1] + sum of nonzero coefficients of SPC DOFs (only for modal calculations)+
48 * @param auwp The auwp for W_b matrix??
49 * @param jqwp The jqwp for W_b matrix??
50 * @param irowwp The irowwp for W_b matrix??
51 * @param nzsw The nzsw for W_b matrix??
52 * @param islavnode Slave nodes are stored tie by tie
53 * @param nslavnode Position in islavnode before the first slave node of tie i.
54 * @param nslavs The total number of slave nodes
55 * @param imastnode Master nodes catalogued as in **islavnode** of size nmastnode(ntie+1)
56 * @param nmastnode The nmastnode(i) points towards the master node within tie i with the largest number
57 * @param ntie Total number of tie constraints
58 * @param nactdof Active degrees of freedom are stored in a two-dimensional field.
59 * + It has as many rows as there are nodes in the model and four columns since each node has
60 * one temperature degree of freedom and three translational degrees.
61 * + In C this field is mapped into a one-dimensional field starting
62 * with the degrees of freedom of node 1, then those of node 2, and so on.
63 * @param mi Maximums for:
64 * + mi(1): max # of integration points per element (max over all elements)
65 * + mi(2): max degree of freedom per node (max over all nodes)
66 * @param vold vold(0..4,1..nk) solution field in all nodes
67 * + 0: temperature
68 * + 1: displacement in global x-direction
69 * + 2: displacement in global y-direction
70 * + 3: displacement in global z-direction
71 * + 4: static pressure
72 * @param nk highest node number
73 * @param fext external mechanical forces in DOF i (due to point loads and
74 distributed loads, including centrifugal and gravity loads,
75 but excluding temperature loading and displacement loading)
76 * @param isolver Sparse linear solver selector
77 * + 'SPOOLES' 0
78 * + 'ITERATIVESCALING' 2
79 * + 'ITERATIVECHOLESKY' 3
80 * + 'SGI' 4
81 * + 'TAUCS' 5
82 * + 'PARDISO' 7
83 * + 'PASTIX' 8
84 * @param iperturb Perturbation?
85 * + iperturb(1)
86 * + = -1 : linear iteration in a nonlinear calculation
87 * + = 0 : linear
88 * + = 1 : second order theory for frequency/buckling/Green calculations following a static step (PERTURBATION selected)
89 * + $ \ge$ 2 : Newton-Raphson iterative procedure is active
90 * + = 3 : nonlinear material (linear or nonlinear geometric and/or heat transfer)
91 * + iperturb(2)
92 * +0 : linear geometric (NLGEOM not selected)
93 * +1 : nonlinear geometric (NLGEOM
94 *
95 */
massless(double * au,double * ad,double * aub,double * adb,ITG * jq,ITG * irow,ITG * neq,ITG * nzs,double ** auwp,ITG ** jqwp,ITG ** irowwp,ITG * nzsw,ITG * islavnode,ITG * nslavnode,ITG * nslavs,ITG * imastnode,ITG * nmastnode,ITG * ntie,ITG * nactdof,ITG * mi,double * vold,ITG * nk,double * fext,ITG * isolver,ITG * iperturb,double * co,double * springarea)96 void massless(double *au, double *ad, double *aub, double *adb, ITG *jq,
97 ITG *irow, ITG *neq, ITG *nzs, double **auwp, ITG **jqwp,
98 ITG **irowwp, ITG *nzsw, ITG *islavnode, ITG *nslavnode,
99 ITG *nslavs, ITG *imastnode, ITG *nmastnode, ITG *ntie,
100 ITG *nactdof, ITG *mi, double *vold, ITG *nk,
101 double *fext, ITG *isolver, ITG *iperturb, double *co, double *springarea) {
102
103 ITG nzswnew, *jqwnew = NULL, *irowwnew = NULL, *jqw = NULL, *iroww = NULL, nmasts,
104 neqbb, *jqbb = NULL, *irowbb = NULL, *kslav = NULL, *lslav = NULL, *ktot = NULL,
105 *ltot = NULL, neqslav, neqtot, *jqbi = NULL, *jqib = NULL, *irowbi = NULL,
106 *irowib = NULL, nzsbi, nzsib, symmetryflag = 0, inputformat = 0,nrhs=1,
107 nzsbb, *icolbb = NULL,*iacti = NULL, nacti=0;
108
109 double *auwnew = NULL, *auw = NULL, *aubb = NULL, *adbb = NULL, *aubi = NULL, *auib = NULL,
110 sigma = 0.0, *gapdof = NULL, *gapnorm = NULL, *qtmp = NULL, *cvec = NULL, sum,
111 *adbbb = NULL, *aubbb = NULL , *gcontvec=NULL, *gcontfull=NULL;
112
113 auw = *auwp;
114 jqw = *jqwp;
115 iroww = *irowwp;
116 nmasts = nmastnode[*ntie];
117
118 /* catalogue slave dofs in kslav,lslav and
119 slave and master dofs in ktot,ltot */
120
121 NNEW(kslav, ITG, 3 * *nslavs);
122 NNEW(lslav, ITG, 3 * *nslavs);
123 NNEW(ktot, ITG, 3 * *nslavs + 3 * nmasts);
124 NNEW(ltot, ITG, 3 * *nslavs + 3 * nmasts);
125
126 // Create set of slave and slave+master contact DOFS. Sorted.
127 FORTRAN(create_contactdofs, (kslav, lslav, ktot, ltot, nslavs, islavnode, &nmasts,
128 imastnode, nactdof, mi, &neqslav, &neqtot)); // TODO: nactdof, mi, to declare
129
130 /* expanding the matrix Wb according to the number of degrees
131 of freedom */
132
133 NNEW(jqwnew, ITG, 3 * *nslavs + 1);
134 NNEW(auwnew, double, 3 * *nzsw);
135 NNEW(irowwnew, ITG, 3 * *nzsw); //TODO: declaring irowwnew?? typo?
136
137 // Expand contact direction matrix W_b from nodal size to DOF size (3x)
138 FORTRAN(expand_auw, (auw, jqw, iroww, nslavs, nzsw, auwnew, jqwnew, irowwnew,
139 &nzswnew, &neqslav, lslav, ntie, nactdof, mi, ktot, &neqtot,
140 islavnode, nslavnode, imastnode));
141
142 *nzsw = nzswnew;
143
144 RENEW(jqw, ITG, 3 * *nslavs + 1);
145 RENEW(auw, double, *nzsw);
146 RENEW(iroww, ITG, *nzsw);
147
148 memcpy(jqw, jqwnew, sizeof(ITG) * (3 * *nslavs + 1));
149 memcpy(auw, auwnew, sizeof(double) * *nzsw); // TODO: aunew auwnew?? declare
150 memcpy(iroww, irowwnew, sizeof(ITG) * *nzsw);
151
152 SFREE(jqwnew);
153 SFREE(auwnew);
154 SFREE(irowwnew);
155
156 /* extracting Kbb,Kbi,Kib,Kii from the stiffness matrix */
157
158 NNEW(jqbb, ITG, neqtot + 1);
159 NNEW(aubb, double, nzs[0]); // QUESTION: here we are allocation for the size of the whole K matrix?
160 NNEW(adbb, double, neqtot);
161 NNEW(irowbb, ITG, nzs[0]);
162
163 NNEW(jqbi, ITG, neq[0] + 1);
164 NNEW(aubi, double, nzs[0]);
165 NNEW(irowbi, ITG, nzs[0]);
166
167 NNEW(jqib, ITG, neqtot + 1);// QUESTION: do we need also the IB partition? can we just use the transposed BI?
168 NNEW(auib, double, nzs[0]);
169 NNEW(irowib, ITG, nzs[0]);
170
171 NNEW(icolbb, ITG, neqtot);
172 // NNEW(icontactv, ITG, neqtot + 1); // is this size good?
173
174 // extracting the submatrices from the global stiffness matrix
175 FORTRAN(extract_matrices, (au, ad, jq, irow, neq, nzs, aubb, adbb, jqbb, irowbb, &neqtot,
176 &nzsbb, islavnode, nslavs, imastnode, &nmasts, aubi, jqbi,
177 irowbi, &nzsbi, auib, jqib, irowib, &nzsib, kslav, ktot, icolbb));
178
179 RENEW(aubb, double, nzsbb);
180 RENEW(irowbb, ITG, nzsbb);
181 // RENEW(icolbb, ITG, neqtot); /* length of each column icolbb(i)=jqbb(i+1)-jqbb(i) */
182 RENEW(aubi, double, nzsbi);
183 RENEW(irowbi, ITG, nzsbi);
184 RENEW(auib, double, nzsib);
185 RENEW(irowib, ITG, nzsib);
186
187 /* factorize kbb */
188 NNEW(gapdof, double, neqtot);
189 NNEW(gapnorm, double, *nslavs);
190
191 MNEW(qtmp, double, neq[0]);
192
193 FORTRAN(detectactivecont1, (vold, nk, mi, aubi, irowbi, jqbi, &neqtot, fext, aubb, adbb, ltot,
194 irowbb, jqbb, auw, iroww, jqw, &neqslav, gapdof,
195 auib, irowib, jqib, icolbb, nactdof, qtmp, neq, co));
196
197 /* call spooles_solve..... */
198
199 if (*isolver == 0) {
200 #ifdef SPOOLES
201 // spooles(adbb,aubb,adbbb,aubbb,&sigma,gapdof,icolbb,irowbb,&neqtot,&nzsbb,&symmetryflag,
202 // &inputformat,&nzsbb);
203 spooles_factor(adbb,aubb,adbbb,aubbb,&sigma,icolbb,irowbb,
204 &neqtot,&nzsbb,&symmetryflag, &inputformat,&nzsbb);
205
206 spooles_solve(gapdof,&neqtot);
207 #else
208 printf("*ERROR in linstatic: the SPOOLES library is not linked\n\n");
209 FORTRAN(stop, ());
210 #endif
211 } else if ((*isolver == 2) || (*isolver == 3)) {
212 preiter(ad, &au, gapnorm, &icolbb, &irow, neq, nzs, isolver, iperturb);
213 } else if (*isolver == 4) {
214 #ifdef SGI
215 token=1;
216 sgapdof_main(ad,au,adb,aub,&sigma,g,icolbb,irow,neq,nzs,token);
217 #else
218 printf("*ERROR in linstatic: the SGI library is not linked\n\n");
219 FORTRAN(stop, ());
220 #endif
221 } else if (*isolver == 5) {
222 #ifdef TAUCS
223 tau(ad,&au,adb,aub,&sigma,g,icolbb,&irow,neq,nzs);
224 #else
225 printf("*ERROR in linstatic: the TAUCS library is not linked\n\n");
226 FORTRAN(stop, ());
227 #endif
228 } else if (*isolver == 7) {
229 #ifdef PARDISO
230 pardiso_main(adbb,aubb,adbbb,aubbb,&sigma,gapdof,icolbb,irowbb,&neqtot,nzsbb,
231 &symmetryflag,&inputformat,jqbb,&nzs[2],&nrhs); // TODO check PARDISO!
232 #else
233 printf("*ERROR in linstatic: the PARDISO library is not linked\n\n");
234 FORTRAN(stop, ());
235 #endif
236 } else if (*isolver == 8) {
237 #ifdef PASTIX
238 pastix_main(ad,au,adb,aub,&sigma,gapdof,icolbb,irow,neq,nzs,
239 &symmetryflag,&inputformat,jq,&nzs[2],&nrhs);
240 #else
241 printf("*ERROR in linstatic: the PASTIX library is not linked\n\n");
242 FORTRAN(stop, ());
243 #endif
244 }
245
246 NNEW(iacti, ITG, *nslavs);
247 /* premultiply g by Wb^T and add g0 => determine active degrees => reduce g to c */
248 FORTRAN(detectactivecont2, (gapnorm, gapdof, auw, iroww, jqw, &neqtot, nslavs, springarea, iacti,&nacti));
249
250
251 // ***************** G ("DELASSUS") MATRIX + Cvector NORMAL CONTACT ***********************++
252 int ncdim = 1;
253
254
255 // initializing cvec with values of gnorm
256 NNEW(cvec,double,nacti);
257 for(int i=0;i<*nslavs;i++){
258 if(iacti[i]!=0){
259 cvec[iacti[i]-1]=gapnorm[i];
260 }
261 }
262
263 NNEW(gcontfull, double, nacti * nacti);
264
265 // calculate Wb^T.Kbb^(-1).Wb
266 int indexnorm ;
267 for (int i = 0; i < *nslavs; i++) {
268 if(iacti[i]!=0){
269 indexnorm= 3*i ; // normal contact index
270 //initialize in zeros vector
271 NNEW(gcontvec, double, neqtot);
272 for(int j=jqw[indexnorm]-1;j<jqw[indexnorm+1]-1;j++){
273 gcontvec[iroww[j]-1]=auw[j];
274 }
275 }
276
277 spooles_solve(gcontvec,&neqtot); // TODO add other solvers!!!!
278
279 for (int j = 0; j < *nslavs; ++j) {
280 if(iacti[j]!=0){
281 indexnorm= 3*j ;
282 sum=0.0;
283 for (int k = jqw[indexnorm]-1; k < jqw[indexnorm+1]-1 ; k++) {
284 sum+=auw[k]*gcontvec[iroww[k]-1];
285 }
286 gcontfull[i*nacti+j]=sum;
287 }
288 }
289 }
290
291 FORTRAN(relaxval_al, (gcontfull,&nacti, &ncdim )); // last parameter: nCdim: 1:NORMALonly, 3:NTT
292
293
294 // conttype: contact type
295 int conttype = 1 ;//'NORMAL'
296 double mufric = 0.0;
297 double atol = 1.0e-8;
298 double rtol = 1.0e-8;
299 double *pkvec = NULL;
300 int kitermax = 1000;
301 int timek = 0.0;
302
303 FORTRAN( auglag_inclusion ,(conttype, gcontfull,&nacti, &ncdim, mufric, atol,rtol, pkvec, kitermax, timek ));
304
305
306
307 if (*isolver == 0) {
308 #ifdef SPOOLES
309 spooles_cleanup();
310 #endif
311 }
312
313 SFREE(kslav);
314 SFREE(lslav);
315 SFREE(ktot);
316 SFREE(ltot);
317
318 *auwp = auw;
319 *jqwp = jqw;
320 *irowwp = iroww; // QUESTION: why ?
321
322 SFREE(gcontvec);
323 SFREE(gcontfull);
324
325 return;
326 }
327