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 <time.h>
22 #include "CalculiX.h"
23 #include "mortar.h"
24
25 /**
26 * function called before solver transforming the SPCs/MPCs, the matrix and the right hand side for quadratic elements
27 (see Phd-thesis Sitzmann,Chapter 4)
28
29 * * Author: Saskia Sitzmann
30 *
31 * [out] iflagact flag indicating, whether the coupling matrices have to
32 be recalculated
33 * [out] nzsc2 number of nonzero,nondiagonal entries of intermediate
34 system matrix
35 * [out] auc2p intermediate system matrix
36 * [out] adc2p intermediate system matrix, diagonal terms
37 * [out] irowc2p rows for intermediate system matrix
38 * [out] icolc2p columns for intermediate system matrix
39 * [out] jqc2p pointer to irowc
40 * [out] aubdp coupling matrix \f$ B_d[nactdof(i,p),nactdof(j,q)]\f$
41 for all active degrees od freedoms
42 * [out] irowbdp field containing row numbers of aubd
43 * [out] jqbdp pointer into field irowbd
44 * [out] aubdtilp matrix \f$ \tilde{D}^{-1}\tilde{B}_d[nactdof(i,p),
45 nactdof(j,q)]\f$ for all active degrees od freedoms
46 * [out] irowbdtilp field containing row numbers of aubd
47 * [out] jqbdtilp pointer into field irowbdtil
48 * [out] aubdtil2p coupling matrix \f$ \tilde{D}$ and
49 $\tilde{B}^2_d[nactdof(i,p),nactdof(j,q)]\f$ for all
50 active degrees od freedoms
51 * [out] irowbdtil2p field containing row numbers of aubdtil2
52 * [out] jqbdtil2p pointer into field irowbdtil2
53 * [out] auddp coupling matrix \f$ D_d[nactdof(i,p),nactdof(j,q)]\f$
54 for all active degrees od freedoms
55 * [out] irowddp field containing row numbers of audd
56 * [out] jqddp pointer into field irowdd
57 * [out] auddtilp coupling matrix \f$ \tilde{D}_d[nactdof(i,p),nactdof(j,
58 q)]\f$ for all active degrees od freedoms
59 * [out] irowddtilp field containing row numbers of audd
60 * [out] jqddtilp pointer into field irowdd
61 * [out] auddtil2p matrix \f$ Id_d[nactdof(i,p),nactdof(j,q)]\f$ for all
62 active degrees od freedoms
63 * [out] irowddtil2p field containing row numbers of audd
64 * [out] jqddtil2p pointer into field irowdd
65 * [out] auddinvp coupling matrix \f$ \tilde{D}^{-1}_d[nactdof(i,p),
66 nactdof(j,q)]\f$ for all active degrees od freedoms
67 * [out] irowddinvp field containing row numbers of auddinv
68 * [out] jqddinvp pointer into field irowddinv
69 * [out] jqtempp field storing the untransformed stiffness matrix
70 representation
71 * [out] irowtempp field storing the untransformed stiffness matrix
72 representation
73 * [out] icoltempp field storing the untransformed stiffness matrix
74 representation
75 * [out] nzstemp field storing the untransformed stiffness matrix size
76 * [out] slavnor slave normal
77 * [out] slavtan slave tangent
78 * [out] nboun2 number of transformed SPCs
79 * [out] ndirboun2p (i) direction of transformed SPC i
80 * [out] nodeboun2p (i) node of transformed SPC i
81 * [out] xboun2p (i) value of transformed SPC i
82 * [out] nmpc2 number of transformed mpcs
83 * [out] ipompc2p (i) pointer to nodempc and coeffmpc for transformed
84 MPC i
85 * [out] nodempc2p i and directions of transformed MPCs
86 * [out] coefmpc2p coefficients of transformed MPCs
87 * [out] labmpc2p transformed mpc labels
88 * [out] ikboun2p sorted dofs idof=8*(node-1)+dir for transformed SPCs
89 * [out] ilboun2p transformed SPC numbers for sorted dofs
90 * [out] ikmpc2p sorted dofs idof=8*(node-1)+dir for transformed MPCs
91 * [out] ilmpc2p transformed SPC numbers for sorted dofs
92 * [out] nslavspcp (2*i) pointer to islavspc...
93 * [out] islavspcp ... which stores SPCs for slave node i
94 * [out] nsspc number of SPC for slave i
95 * [out] nslavmpcp (2*i) pointer to islavmpc...
96 * [out] islavmpcp ... which stores MPCs for slave node i
97 * [out] nsmpc number of MPC for slave i
98 * [out] nslavspc2p (2*i) pointer to islavspc2...
99 * [out] islavspc2p ... which stores transformed SPCs for slave node i
100 * [out] nsspc2 number of transformed SPC for slave i
101 * [out] nslavmpc2p (2*i) pointer to islavmpc2...
102 * [out] islavmpc2p ... which stores transformed MPCs for slave node i
103 * [out] nsmpc2 number of transformed MPC for slave i
104 * [in] imastnode field storing the i of the master surfaces
105 * [in] nmastnode (i)pointer into field imastnode for contact tie i
106 * [out] nmastspcp (2*i) pointer to imastspc...
107 * [out] imastspcp ... which stores SPCs for master node i
108 * [out] nmspc number of SPC for master i
109 * [out] nmastmpcp (2*i) pointer to imastmpc...
110 * [out] imastmpcp ... which stores MPCs for master node i
111 * [out] nmmpc number of MPC for master i
112 * [in] islavelinv (i)==0 if there is no slave node in the element,
113 >0 otherwise
114 * [in] islavsurf islavsurf(1,i) slaveface i islavsurf(2,i) # integration
115 points generated before looking at face i
116 * [in] autloc transformation matrix \f$ T[p,q]\f$ for slave i \f$ p,q
117 \f$
118 * [in] irowtloc field containing row numbers of autloc
119 * [in] jqtloc pointer into field irowtloc
120 * [in] autlocinv transformation matrix \f$ T^{-1}[p,q]\f$ for slave i
121 \f$ p,q \f$
122 * [in] irowtlocinv field containing row numbers of autlocinv
123 * [in] jqtlocinv pointer into field irowtlocinv
124 * [in] nk2 number or generated points needed for transformed SPCs
125 * [in] iflagdualquad flag indicating what mortar contact is used
126 (=1 quad-lin, =2 quad-quad, =3 PG quad-lin, =4 PG
127 quad-quad)
128 * [out] f_cmp not used any more
129 * [out] f_csp contact force for active degrees of freedom
130 * [in] auxtil2 auxilary field
131 * [in] pslavsurf field storing position xil, etal and weight for
132 integration point on slave side
133 * [in] pmastsurf field storing position and etal for integration points
134 on master side
135 * [in] islavact (i) indicates, if slave node i is active
136 (=-3 no-slave-node, =-2 no-LM-node, =-1 no-gap-node,
137 =0 inactive node, =1 sticky node, =2 slipping/active
138 node)
139 * [in] cdn ?
140 * [in] cvtilini C*v at start of the increment
141 * [in] cvtil C*v
142 * [in] idamping flag indicating whether damping is used
143 * [in] ilin flag indicating wheter fist iteration is calculated
144 linear geometrically
145 * [in] iperturb_sav saved iperturb values
146 * [out] nodeforc2p transformed point force, node
147 * [out] ndirforc2p transformed point force, dir
148 * [out] xforc2p transformed point force, value
149 * [out] nforc2 number of transformed point forces
150 **/
151
premortar(ITG * iflagact,ITG * ismallsliding,ITG * nzs,ITG * nzsc2,double ** auc2p,double ** adc2p,ITG ** irowc2p,ITG ** icolc2p,ITG ** jqc2p,double ** aubdp,ITG ** irowbdp,ITG ** jqbdp,double ** aubdtilp,ITG ** irowbdtilp,ITG ** jqbdtilp,double ** aubdtil2p,ITG ** irowbdtil2p,ITG ** jqbdtil2p,double ** auddp,ITG ** irowddp,ITG ** jqddp,double ** auddtilp,ITG ** irowddtilp,ITG ** jqddtilp,double ** auddtil2p,ITG ** irowddtil2p,ITG ** jqddtil2p,double ** auddinvp,ITG ** irowddinvp,ITG ** jqddinvp,ITG ** jqtempp,ITG ** irowtempp,ITG ** icoltempp,ITG * nzstemp,ITG * iit,double * slavnor,double * slavtan,ITG * icol,ITG * irow,ITG * jq,ITG * ikboun,ITG * ilboun,ITG * ikmpc,ITG * ilmpc,ITG * nboun2,ITG ** ndirboun2p,ITG ** nodeboun2p,double ** xboun2p,ITG * nmpc2,ITG ** ipompc2p,ITG ** nodempc2p,double ** coefmpc2p,char ** labmpc2p,ITG ** ikboun2p,ITG ** ilboun2p,ITG ** ikmpc2p,ITG ** ilmpc2p,ITG ** nslavspcp,ITG ** islavspcp,ITG ** nslavmpcp,ITG ** islavmpcp,ITG ** nslavspc2p,ITG ** islavspc2p,ITG ** nslavmpc2p,ITG ** islavmpc2p,ITG ** nmastspcp,ITG ** imastspcp,ITG ** nmastmpcp,ITG ** imastmpcp,ITG ** nmastmpc2p,ITG ** imastmpc2p,ITG * nmmpc2,ITG * nsspc,ITG * nsspc2,ITG * nsmpc,ITG * nsmpc2,ITG * imastnode,ITG * nmastnode,ITG * nmspc,ITG * nmmpc,double * co,ITG * nk,ITG * kon,ITG * ipkon,char * lakon,ITG * ne,double * stn,double * elcon,ITG * nelcon,double * rhcon,ITG * nrhcon,double * alcon,ITG * nalcon,double * alzero,ITG * ielmat,ITG * ielorien,ITG * norien,double * orab,ITG * ntmat_,double * t0,double * t1,ITG * ithermal,double * prestr,ITG * iprestr,char * filab,double * eme,double * emn,double * een,ITG * iperturb,double * f,ITG * nactdof,ITG * iout,double * qa,double * vold,double * b,ITG * nodeboun,ITG * ndirboun,double * xbounact,double * xboun,ITG * nboun,ITG * ipompc,ITG * nodempc,double * coefmpc,char * labmpc,ITG * nmpc,ITG * nmethod,ITG * neq,double * veold,double * accold,double * dtime,double * time,double * ttime,double * plicon,ITG * nplicon,double * plkcon,ITG * nplkcon,double * xstateini,double * xstiff,double * xstate,ITG * npmat_,char * matname,ITG * mi,ITG * ielas,ITG * icmd,ITG * ncmat_,ITG * nstate_,double * stiini,double * vini,double * ener,double * enern,double * emeini,double * xstaten,double * eei,double * enerini,double * cocon,ITG * ncocon,char * set,ITG * nset,ITG * istartset,ITG * iendset,ITG * ialset,ITG * nprint,char * prlab,char * prset,double * qfx,double * qfn,double * trab,ITG * inotr,ITG * ntrans,ITG * nelemload,ITG * nload,ITG * istep,ITG * iinc,double * springarea,double * reltime,ITG * ne0,double * xforc,ITG * nforc,double * thicke,double * shcon,ITG * nshcon,char * sideload,double * xload,double * xloadold,ITG * icfd,ITG * inomat,ITG * islavelinv,ITG * islavsurf,ITG * iponoels,ITG * inoels,ITG * mortar,ITG * nslavnode,ITG * islavnode,ITG * nslavs,ITG * ntie,double * autloc,ITG * irowtloc,ITG * jqtloc,double * autlocinv,ITG * irowtlocinv,ITG * jqtlocinv,ITG * nk2,ITG * iflagdualquad,char * tieset,ITG * itiefac,ITG * rhsi,double * au,double * ad,double ** f_cmp,double ** f_csp,double * t1act,double * cam,double * bet,double * gam,double * epn,double * xloadact,ITG * nodeforc,ITG * ndirforc,double * xforcact,double * xbodyact,ITG * ipobody,ITG * nbody,double * cgr,ITG * nzl,double * sti,ITG * iexpl,ITG * mass,ITG * buckling,ITG * stiffness,ITG * intscheme,double * physcon,ITG * coriolis,ITG * ibody,ITG * integerglob,double * doubleglob,ITG * nasym,double * alpham,double * betam,double * auxtil2,double * pslavsurf,double * pmastsurf,double * clearini,ITG * ielprop,double * prop,ITG * islavact,double * cdn,ITG * memmpc_,double * cvinitil,double * cvtil,ITG * idamping,ITG * ilin,ITG * iperturb_sav,double * adb,double * aub,ITG ** nodeforc2p,ITG ** ndirforc2p,double ** xforc2p,ITG * nforc2,ITG * itietri,double * cg,double * straight,ITG * koncont,double * energyini,double * energy,ITG * kscale,ITG * iponoel,ITG * inoel,ITG * nener,char * orname,ITG * network,char * typeboun,ITG * num_cpus,double * t0g,double * t1g,double * smscale,ITG * mscalmethod)152 void premortar(ITG *iflagact,ITG *ismallsliding,ITG *nzs,ITG *nzsc2,
153 double **auc2p,double **adc2p,ITG **irowc2p,ITG **icolc2p,
154 ITG **jqc2p,
155 double **aubdp,ITG **irowbdp,ITG **jqbdp,
156 double **aubdtilp,ITG **irowbdtilp,ITG **jqbdtilp,
157 double **aubdtil2p,ITG **irowbdtil2p,ITG **jqbdtil2p,
158 double **auddp,ITG **irowddp,ITG **jqddp,
159 double **auddtilp,ITG **irowddtilp,ITG **jqddtilp,
160 double **auddtil2p,ITG **irowddtil2p,ITG **jqddtil2p,
161 double **auddinvp,ITG **irowddinvp,ITG **jqddinvp,
162 ITG **jqtempp,ITG **irowtempp,ITG **icoltempp,ITG *nzstemp,
163 ITG *iit,double *slavnor,double *slavtan,
164 ITG *icol,ITG *irow,ITG *jq,
165 ITG *ikboun,ITG *ilboun,ITG *ikmpc,ITG *ilmpc,
166 ITG *nboun2,ITG **ndirboun2p,ITG **nodeboun2p,
167 double **xboun2p,
168 ITG *nmpc2,ITG **ipompc2p,ITG **nodempc2p,double **coefmpc2p,
169 char **labmpc2p,
170 ITG **ikboun2p,ITG **ilboun2p,ITG **ikmpc2p,ITG **ilmpc2p,
171 ITG **nslavspcp,ITG **islavspcp,ITG **nslavmpcp,
172 ITG **islavmpcp,
173 ITG **nslavspc2p,ITG **islavspc2p,ITG **nslavmpc2p,
174 ITG **islavmpc2p,
175 ITG **nmastspcp,ITG **imastspcp,ITG **nmastmpcp,
176 ITG **imastmpcp,
177 ITG **nmastmpc2p,ITG **imastmpc2p,ITG *nmmpc2,
178 ITG *nsspc,ITG *nsspc2,ITG *nsmpc,ITG *nsmpc2,
179 ITG *imastnode,ITG *nmastnode,ITG *nmspc,ITG *nmmpc,
180 double *co,ITG *nk,ITG *kon,ITG *ipkon,char *lakon,
181 ITG *ne,double *stn,
182 double *elcon,ITG *nelcon,double *rhcon,ITG *nrhcon,
183 double *alcon,ITG *nalcon,double *alzero,ITG *ielmat,
184 ITG *ielorien,ITG *norien,double *orab,ITG *ntmat_,
185 double *t0,double *t1,ITG *ithermal,double *prestr,
186 ITG *iprestr,char *filab,double *eme,double *emn,
187 double *een,ITG *iperturb,double *f,ITG *nactdof,
188 ITG *iout,double *qa,
189 double *vold,double *b,ITG *nodeboun,ITG *ndirboun,
190 double *xbounact,double *xboun,ITG *nboun,ITG *ipompc,
191 ITG *nodempc,
192 double *coefmpc,char *labmpc,ITG *nmpc,ITG *nmethod,
193 ITG *neq,double *veold,double *accold,
194 double *dtime,double *time,
195 double *ttime,double *plicon,
196 ITG *nplicon,double *plkcon,ITG *nplkcon,
197 double *xstateini,double *xstiff,double *xstate,ITG *npmat_,
198 char *matname,ITG *mi,ITG *ielas,
199 ITG *icmd,ITG *ncmat_,ITG *nstate_,double *stiini,
200 double *vini,double *ener,
201 double *enern,double *emeini,double *xstaten,double *eei,
202 double *enerini,double *cocon,ITG *ncocon,char *set,
203 ITG *nset,ITG *istartset,
204 ITG *iendset,ITG *ialset,ITG *nprint,char *prlab,
205 char *prset,double *qfx,double *qfn,double *trab,
206 ITG *inotr,ITG *ntrans,ITG *nelemload,
207 ITG *nload,ITG *istep,ITG *iinc,
208 double *springarea,double *reltime,ITG *ne0,double *xforc,
209 ITG *nforc,double *thicke,
210 double *shcon,ITG *nshcon,char *sideload,double *xload,
211 double *xloadold,ITG *icfd,ITG *inomat,
212 ITG *islavelinv,ITG *islavsurf,
213 ITG *iponoels,ITG *inoels,
214 ITG *mortar,ITG *nslavnode,ITG *islavnode,ITG *nslavs,
215 ITG *ntie,
216 double *autloc,ITG *irowtloc,ITG *jqtloc,
217 double *autlocinv,ITG *irowtlocinv,ITG *jqtlocinv,
218 ITG *nk2,ITG *iflagdualquad,
219 char *tieset,ITG *itiefac ,ITG *rhsi,
220 double *au,double *ad,double **f_cmp,double **f_csp,
221 double *t1act,double *cam,double *bet,double *gam,
222 double *epn,
223 double *xloadact,ITG *nodeforc,ITG *ndirforc,double *xforcact,
224 double *xbodyact,ITG *ipobody,ITG *nbody,double *cgr,
225 ITG *nzl,double *sti,ITG *iexpl,ITG *mass,ITG *buckling,
226 ITG *stiffness,
227 ITG *intscheme,double *physcon,ITG *coriolis,ITG *ibody,
228 ITG *integerglob,double *doubleglob,ITG *nasym,
229 double *alpham,double *betam,double *auxtil2,
230 double *pslavsurf,double *pmastsurf,
231 double *clearini,ITG *ielprop,double *prop,
232 ITG *islavact,double *cdn,ITG *memmpc_,
233 double *cvinitil,double *cvtil,ITG *idamping,
234 ITG *ilin,ITG *iperturb_sav,double *adb,double *aub,
235 ITG **nodeforc2p,ITG **ndirforc2p,double **xforc2p,
236 ITG *nforc2,
237 ITG *itietri,double *cg,double *straight,ITG *koncont,
238 double *energyini,
239 double *energy,ITG *kscale,ITG *iponoel,ITG *inoel,ITG *nener,
240 char *orname,ITG *network,
241 char *typeboun,ITG *num_cpus,double *t0g,double *t1g,
242 double *smscale,ITG *mscalmethod){
243
244 ITG im,i,ii,j,jj,k,l,mt=mi[1]+1,node,jfaces,nelems,ifaces,nope,nopes,idummy,
245 *ndirboun2=NULL,*nodeboun2=NULL,*ipompc2=NULL,*nodempc2=NULL,nodes[8],
246 konl[20],jj2,ifac,debug,
247 *ikboun2=NULL,*ilboun2=NULL,*ikmpc2=NULL,*ilmpc2=NULL,
248 *nslavspc=NULL,*islavspc=NULL,*nslavmpc=NULL,*islavmpc=NULL,
249 *nslavspc2=NULL,*islavspc2=NULL,*nslavmpc2=NULL,*islavmpc2=NULL,
250 *nmastspc=NULL,*imastspc=NULL,*nmastmpc=NULL,*imastmpc2=NULL,
251 *nmastmpc2=NULL,*imastmpc=NULL,
252 *irowc2=NULL,*icolc2=NULL,*jqc2=NULL,*irowbd=NULL,*jqbd=NULL,
253 *irowbdtil=NULL,*jqbdtil=NULL,
254 *irowbdtil2=NULL,*jqbdtil2=NULL,*irowdd=NULL,*jqdd=NULL,*irowddtil=NULL,
255 *jqddtil=NULL,
256 *irowddinv=NULL,*jqddinv=NULL,*irowddtil2=NULL,*jqddtil2=NULL,
257 *irowtemp=NULL,*icoltemp=NULL,*jqtemp=NULL,*inum=NULL,*icoltil=NULL,
258 *irowtil=NULL,*jqtil=NULL,
259 *nodeforc2=NULL,*ndirforc2=NULL,mortartrafoflag=1;
260
261 double alpha,*xboun2=NULL,*coefmpc2=NULL,*auc2=NULL,*adc2=NULL,*aubd=NULL,
262 *aubdtil=NULL,*aubdtil2=NULL,*ftil=NULL,*fexttil=NULL,
263 *audd=NULL,*auddtil=NULL,*auddinv=NULL,*auddtil2=NULL,*v=NULL,
264 *stx=NULL,*fn=NULL,*autil=NULL,*finitil=NULL,*fextinitil=NULL,
265 *adtil=NULL,*fmpc2=NULL,*voldtil=NULL,*vinitil=NULL,*accoldtil=NULL,
266 *btil=NULL,*aubtil=NULL,*adbtil=NULL,
267 *adctil=NULL,*auctil=NULL,*veoldtil=NULL,*volddummy=NULL,
268 *vectornull=NULL,*f_cs=NULL,*f_cm=NULL,
269 *xforc2=NULL,*fnext=NULL,*veolddummy=NULL,*accolddummy=NULL;
270
271 char *labmpc2=NULL;
272
273 debug=0;
274
275 alpha=1-2*sqrt(*bet);
276
277 ndirboun2=*ndirboun2p;nodeboun2=*nodeboun2p;xboun2=*xboun2p;
278 ipompc2=*ipompc2p;nodempc2=*nodempc2p;coefmpc2=*coefmpc2p;
279 labmpc2=*labmpc2p;ikboun2=*ikboun2p;ilboun2=*ilboun2p;ikmpc2=*ikmpc2p;
280 ilmpc2=*ilmpc2p;
281 nslavspc=*nslavspcp;islavspc=*islavspcp;nslavmpc=*nslavmpcp;
282 islavmpc=*islavmpcp;
283 nslavspc2=*nslavspc2p;islavspc2=*islavspc2p;nslavmpc2=*nslavmpc2p;
284 islavmpc2=*islavmpc2p;
285 nmastspc=*nmastspcp;imastspc=*imastspcp;nmastmpc=*nmastmpcp;
286 imastmpc=*imastmpcp;nmastmpc2=*nmastmpc2p;imastmpc2=*imastmpc2p;
287 auc2=*auc2p;adc2=*adc2p;irowc2=*irowc2p;icolc2=*icolc2p;jqc2=*jqc2p;
288 aubd=*aubdp;irowbd=*irowbdp;jqbd=*jqbdp;
289 aubdtil=*aubdtilp;irowbdtil=*irowbdtilp;jqbdtil=*jqbdtilp;
290 aubdtil2=*aubdtil2p;irowbdtil2=*irowbdtil2p;jqbdtil2=*jqbdtil2p;
291 audd=*auddp;irowdd=*irowddp;jqdd=*jqddp;
292 auddtil=*auddtilp;irowddtil=*irowddtilp;jqddtil=*jqddtilp;
293 auddtil2=*auddtil2p;irowddtil2=*irowddtil2p;jqddtil2=*jqddtil2p;
294 auddinv=*auddinvp;irowddinv=*irowddinvp;jqddinv=*jqddinvp;
295 irowtemp=*irowtempp;icoltemp=*icoltempp;jqtemp=*jqtempp;
296 f_cs=*f_csp;f_cm=*f_cmp;
297 ndirforc2=*ndirforc2p;nodeforc2=*nodeforc2p;xforc2=*xforc2p;
298
299 fflush(stdout);
300
301 NNEW(f_cs,double,neq[1]);
302 NNEW(f_cm,double,neq[1]);
303
304 // check for coupled thermo-mechanical calculation
305
306 if(ithermal[0]>1){
307 printf("\tprecontactmortar: coupled thermo-mechanical calculations NOT");
308 printf("supported yet!\n \tPlease use surface-to-surface penalty contact");
309 printf("instead.\n\n STOP!\n");
310 fflush(stdout);
311 FORTRAN(stop,());
312 }
313
314 // fix for linear calculation in first iteration of first increment
315
316 if((*ilin==1)&&(*iit==1)&&(*iinc==1)){
317 *ielas=1;
318 iperturb[0]=-1;
319 iperturb[1]=0;
320 }
321
322 /* small sliding is autometically set active due to combined fix-point
323 Newton approach
324 do NOT change this unless the additional derivates neglected here
325 have been implemented */
326
327 *iflagact=0;
328 *ismallsliding=1;
329 *nzsc2=nzs[1];
330 NNEW(auc2,double,*nzsc2);
331 NNEW(adc2,double,neq[1]);
332 NNEW(irowc2,ITG,*nzsc2);
333 NNEW(icolc2,ITG,neq[1]);
334 NNEW(jqc2,ITG,neq[1]+1);
335 if(*iit==1 || *ismallsliding==0){
336 NNEW(aubd,double,6*nslavnode[*ntie]);
337 NNEW(irowbd,ITG,6*nslavnode[*ntie]);
338 NNEW(jqbd,ITG,neq[1]+1);
339 NNEW(aubdtil,double,6*nslavnode[*ntie]);
340 NNEW(irowbdtil,ITG,6*nslavnode[*ntie]);
341 NNEW(jqbdtil,ITG,neq[1]+1);
342 NNEW(aubdtil2,double,6*nslavnode[*ntie]);
343 NNEW(irowbdtil2,ITG,6*nslavnode[*ntie]);
344 NNEW(jqbdtil2,ITG,neq[1]+1);
345 NNEW(audd,double,3*nslavnode[*ntie]);
346 NNEW(irowdd,ITG,3*nslavnode[*ntie]);
347 NNEW(jqdd,ITG,neq[1]+1);
348 NNEW(auddtil,double,3*nslavnode[*ntie]);
349 NNEW(irowddtil,ITG,3*nslavnode[*ntie]);
350 NNEW(jqddtil,ITG,neq[1]+1);
351 NNEW(auddtil2,double,3*nslavnode[*ntie]);
352 NNEW(irowddtil2,ITG,3*nslavnode[*ntie]);
353 NNEW(jqddtil2,ITG,neq[1]+1);
354 NNEW(auddinv,double,3*nslavnode[*ntie]);
355 NNEW(irowddinv,ITG,3*nslavnode[*ntie]);
356 NNEW(jqddinv,ITG,neq[1]+1);
357 }
358
359 /* storing the original stiffness matrix */
360
361 NNEW(jqtemp,ITG,neq[1]+1);
362 NNEW(irowtemp,ITG,nzs[1]);
363 NNEW(icoltemp,ITG,neq[1]);
364 for(i=0;i<3;i++){
365 nzstemp[i]=nzs[i];}
366 for (i=0;i<neq[1];i++){jqtemp[i]=jq[i];icoltemp[i]=icol[i];}
367 jqtemp[neq[1]]=jq[neq[1]];
368 for (i=0;i<nzs[1];i++){irowtemp[i]=irow[i];}
369
370 if(*iit==1 || *ismallsliding==0){
371 DMEMSET(slavnor,0,3*nslavnode[*ntie],0.);
372 DMEMSET(slavtan,0,6*nslavnode[*ntie],0.);
373 }
374
375 /* setting iflagact=1 before calling contactmortar invokes
376 combined fix-point Newton approach */
377
378 if(*iit>1 && *ismallsliding==1){*iflagact=1;}
379
380 /* transform SPCs/MPCs in case of quadratic finite elements
381 Caution: there is still a problem with MPCs on slave mid nodes,
382 avoid this!!! */
383
384 //wird immer aufgerufen,da sich xbounact geaendert haben kann
385
386 transformspcsmpcs_quad(nboun,ndirboun,nodeboun,xbounact,
387 nmpc,ipompc,nodempc,coefmpc,labmpc,
388 ikboun,ilboun,ikmpc,ilmpc,
389 nboun2,&ndirboun2,&nodeboun2,&xboun2,
390 nmpc2,&ipompc2,&nodempc2,&coefmpc2,&labmpc2,
391 &ikboun2,&ilboun2,&ikmpc2,&ilmpc2,
392 irowtlocinv,jqtlocinv,autlocinv,
393 nk,nk2,iflagdualquad,
394 ntie,tieset,itiefac,islavsurf,
395 lakon,ipkon,kon,&mt,memmpc_,
396 nodeforc,ndirforc,xforcact,nforc,
397 &nodeforc2,&ndirforc2,&xforc2,nforc2);
398
399 RENEW(islavspc,ITG,2**nboun);
400 RENEW(islavmpc,ITG,2**nmpc);
401 RENEW(islavspc2,ITG,2**nboun2);
402 RENEW(islavmpc2,ITG,2**nmpc2);
403 RENEW(imastspc,ITG,2**nboun);
404 RENEW(imastmpc,ITG,2**nmpc);
405 RENEW(imastmpc2,ITG,2**nmpc);
406
407 /* cataloque SPCs/MPCs */
408
409 FORTRAN(catsmpcslavno,(ntie,islavnode,imastnode,nslavnode,
410 nmastnode,nboun,ndirboun,nodeboun,nmpc,
411 ipompc,nodempc,ikboun,ilboun,ikmpc,ilmpc,
412 nboun2,nmpc2,ipompc2,nodempc2,ikboun2,
413 ilboun2,ikmpc2,ilmpc2,nslavspc,islavspc,
414 nsspc,nslavmpc,islavmpc,nsmpc,nslavspc2,
415 islavspc2,nsspc2,nslavmpc2,islavmpc2,
416 nsmpc2,nmastspc,imastspc,nmspc,nmastmpc,
417 imastmpc,nmmpc,nmastmpc2,imastmpc2,
418 nmmpc2));
419
420 RENEW(islavspc,ITG,2**nsspc+1);
421 RENEW(islavmpc,ITG,2**nsmpc+1);
422 RENEW(islavspc2,ITG,2**nsspc2+1);
423 RENEW(islavmpc2,ITG,2**nsmpc2+1);
424 RENEW(imastspc,ITG,2**nmspc+1);
425 RENEW(imastmpc,ITG,2**nmmpc+1);
426 RENEW(imastmpc2,ITG,2**nmmpc2+1);
427
428 /* Check for additional MPCs on slave mid nodes */
429
430 if(*iit==1 && *iinc==1){
431
432 FORTRAN(nortanslav,(tieset,ntie,ipkon,kon,lakon,set,co,vold,nset,
433 islavsurf,itiefac,islavnode,nslavnode,slavnor,slavtan,
434 mi));
435
436 // call checkspsmpc
437
438 FORTRAN(checkspcmpc,(ntie,tieset,islavnode,imastnode,nslavnode,nmastnode,
439 slavnor,islavact,nboun,ndirboun,xboun,
440 nodempc,coefmpc,ikboun,ilboun,nmpc2,ipompc2,nodempc2,
441 nslavspc,islavspc,nsspc,nslavmpc,islavmpc,nsmpc,
442 nmspc,nmastmpc,imastmpc,nmmpc));
443
444 for (i=0;i<*ntie;i++){
445 if(tieset[i*(81*3)+80]=='C'){
446 for(j=nslavnode[i];j<nslavnode[i+1];j++){
447 node=islavnode[j];
448 int checkformidnode=0;
449 for(jj=nslavmpc[2*(j)];jj<nslavmpc[2*(j)+1];jj++){
450 if(islavmpc[2*jj+1]==-2){
451
452 // MPC cannot be identified
453
454 checkformidnode=1;
455 if(debug==1)printf(" check for mid node %"ITGFORMAT"\n",node);
456 }
457 }
458
459 if(checkformidnode==1){
460
461 //check for mid node
462
463 do{
464
465 for(l=itiefac[2*i];l<=itiefac[2*i+1];l++){
466 ifaces = islavsurf[2*(l-1)+0];
467 nelems = (ITG)(ifaces/10);
468 jfaces = ifaces - nelems*10;
469 FORTRAN(getnumberofnodes,(&nelems,&jfaces,lakon,&nope,&nopes,
470 &idummy));
471 for(jj=0;jj<nope;jj++){
472 konl[jj]=kon[ipkon[nelems-1]+jj];
473 }
474 for(jj=0;jj<nopes;jj++){
475 jj2=jj+1;
476 ifac=FORTRAN(getlocno,(&jj2,&jfaces,&nope));
477 nodes[jj]=konl[ifac-1];
478 }
479 ii=-1;
480 for(jj=0;jj<nopes;jj++){
481 if(nodes[jj]==node){ii=jj;}
482 }
483 if(ii>-1){
484 break;
485 }
486 }
487
488 break;
489 }while(1);
490
491 if((ii>2 && nopes==6)||(ii>3 && nopes==8)){
492
493 // mid node found with extra not supported mpc ->error
494
495 printf(" precontactmortar: Problem with slave mid node \n\n");
496 printf(" *ERROR: Slave mid node %"ITGFORMAT" has",node);
497 printf(" additional MPC which is not a directional blocking ");
498 printf("MPC or a 1-to-1 cyclic symmetry MPC. \n");
499 printf("\t\t This is not supported yet!!!!!!!!!!!!!\n");
500 fflush(stdout);
501 FORTRAN(stop,());
502 }
503 }
504 }
505 }
506 }
507
508 if(*iit==1 || *ismallsliding==0){
509 DMEMSET(slavnor,0,3*nslavnode[*ntie],0.);
510 DMEMSET(slavtan,0,6*nslavnode[*ntie],0.);
511 }
512 if(debug==1)printf(" precontactmortar: MPC-check OK\n\n");
513 }
514
515 /* fix for quadratic FE */
516
517 NNEW(v,double,mt**nk);
518 NNEW(stx,double,6*mi[0]**ne);
519 NNEW(fn,double,mt*(*nk+*nk2));
520 NNEW(fmpc2,double,*nmpc2);
521 memcpy(&v[0],&vold[0],sizeof(double)*mt**nk);
522 *iout=-1;
523 NNEW(ftil,double,neq[1]);
524 NNEW(fexttil,double,neq[1]);
525 NNEW(inum,ITG,*nk);
526
527 /* calculating the internal forces and tangent stiffness using
528 modified shape functions for quadratic elements */
529
530 if(debug==1)printf(" precontactmortar: call results\n");
531 results(co,nk,kon,ipkon,lakon,ne,v,stn,inum,stx,
532 elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,
533 ielorien,norien,orab,ntmat_,t0,t1act,ithermal,
534 prestr,iprestr,filab,eme,emn,een,iperturb,
535 ftil,fn,nactdof,iout,qa,vold,b,nodeboun,
536 ndirboun,xbounact,nboun,ipompc,
537 nodempc,coefmpc,labmpc,nmpc,nmethod,cam,&neq[1],veold,accold,
538 bet,gam,dtime,time,ttime,plicon,nplicon,plkcon,nplkcon,
539 xstateini,xstiff,xstate,npmat_,epn,matname,mi,ielas,icmd,
540 ncmat_,nstate_,stiini,vini,ikboun,ilboun,ener,enern,emeini,
541 xstaten,eei,enerini,cocon,ncocon,set,nset,istartset,iendset,
542 ialset,nprint,prlab,prset,qfx,qfn,trab,inotr,ntrans,fmpc2,
543 nelemload,nload,ikmpc,ilmpc,istep,iinc,springarea,
544 reltime,ne0,thicke,shcon,nshcon,
545 sideload,xloadact,xloadold,icfd,inomat,pslavsurf,pmastsurf,
546 mortar,islavact,cdn,islavnode,nslavnode,ntie,clearini,
547 islavsurf,ielprop,prop,energyini,energy,kscale,iponoel,
548 inoel,nener,orname,network,ipobody,xbodyact,ibody,typeboun,
549 itiefac,tieset,smscale,mscalmethod,nbody,t0g,t1g,
550 islavelinv,autloc,irowtloc,jqtloc,nboun2,
551 ndirboun2,nodeboun2,xboun2,nmpc2,ipompc2,nodempc2,coefmpc2,
552 labmpc2,ikboun2,ilboun2,ikmpc2,ilmpc2,&mortartrafoflag,
553 intscheme);
554
555 if(debug==1)printf(" precontactmortar: results_dstil finished\n");
556 fflush(stdout);
557
558 SFREE(v);SFREE(stx);SFREE(fn);SFREE(inum);SFREE(fmpc2);
559 *iout=0;
560
561 *rhsi=1;
562 NNEW(adtil,double,neq[1]);
563 NNEW(autil,double,nzs[1]);
564 NNEW(irowtil,ITG,nzs[1]);
565 NNEW(jqtil,ITG,neq[1]+1);
566 NNEW(icoltil,ITG,neq[1]);
567 for(i=0;i<neq[1]+1;i++){
568 jqtil[i]=jq[i];}
569 for(i=0;i<neq[1];i++){
570 icoltil[i]=icol[i];}
571 for(i=0;i<nzs[1];i++){
572 irowtil[i]=irow[i];}
573
574 if(debug==1)printf(" precontactmortar: call mafillsmmain\n");
575 fflush(stdout);
576
577 /* calculating the external forces fext and stiffness matrix au/ad using
578 modified shape functions for quadratic elements */
579
580 mafillsmmain(co,nk,kon,ipkon,lakon,ne,nodeboun2,ndirboun2,xboun2,nboun2,
581 ipompc2,nodempc2,coefmpc2,nmpc2,nodeforc2,ndirforc2,xforc2,
582 nforc2,nelemload,sideload,xloadact,nload,xbodyact,ipobody,
583 nbody,cgr,adtil,autil,fexttil,nactdof,icol,jqtil,irowtil,
584 neq,nzl,
585 nmethod,ikmpc2,ilmpc2,ikboun2,ilboun2,
586 elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,
587 ielmat,ielorien,norien,orab,ntmat_,
588 t0,t1act,ithermal,prestr,iprestr,vold,iperturb,sti,
589 nzs,stx,adb,aub,iexpl,plicon,nplicon,plkcon,nplkcon,
590 xstiff,npmat_,dtime,matname,mi,
591 ncmat_,mass,stiffness,buckling,rhsi,intscheme,
592 physcon,shcon,nshcon,cocon,ncocon,ttime,time,istep,iinc,
593 coriolis,ibody,xloadold,reltime,veold,springarea,nstate_,
594 xstateini,xstate,thicke,integerglob,doubleglob,
595 tieset,istartset,iendset,ialset,ntie,nasym,pslavsurf,
596 pmastsurf,mortar,clearini,ielprop,prop,ne0,fnext,kscale,
597 iponoel,inoel,network,ntrans,inotr,trab,smscale,
598 mscalmethod,set,nset,islavelinv,autloc,irowtloc,jqtloc,
599 &mortartrafoflag);
600
601 if(debug==1)printf(" precontactmortar: mafillsmmain_dstil finished\n");
602 fflush(stdout);
603
604 for(i=0;i<neq[1];i++){
605 ad[i]=adtil[i];
606 }
607 for(i=0;i<nzs[1];i++){
608 au[i]=autil[i];
609 }
610
611 /* transform vold,vini,accold for dynamic calculations */
612
613 NNEW(voldtil,double,mt**nk);
614 NNEW(veoldtil,double,mt**nk);
615 NNEW(vinitil,double,mt**nk);
616 NNEW(accoldtil,double,mt**nk);
617 NNEW(btil,double,neq[1]);
618
619 if(debug==1)printf(" precontactmortar: call calcresidual\n");
620
621 /* calculating the residual b */
622
623 calcresidual(nmethod,neq,btil,fexttil,ftil,iexpl,nactdof,auxtil2,voldtil,
624 vinitil,dtime,accoldtil,nk,adbtil,aubtil,jqtil,irowtil,nzl,
625 &alpha,fextinitil,finitil,islavnode,nslavnode,mortar,ntie,f_cm,
626 f_cs,mi,nzs,nasym,idamping,veoldtil,adctil,auctil,cvinitil,cvtil,
627 alpham,num_cpus);
628
629 if(debug==1)printf(" precontactmortar: calcresidual finished\n");
630 fflush(stdout);
631
632 for(k=0;k<neq[1];k++){
633 b[k]=btil[k];}
634
635 SFREE(btil);SFREE(ftil);SFREE(fexttil);
636
637 SFREE(voldtil);SFREE(vinitil);SFREE(accoldtil);SFREE(veoldtil);
638
639 SFREE(adtil);SFREE(autil);SFREE(irowtil);SFREE(jqtil);SFREE(icoltil);
640
641 /* update vold due to spcs to get gap right for rigid body movements */
642
643 if(*iinc==1 && *iit==1 &&*nmethod!=4){
644 NNEW(v,double,mt**nk);
645 NNEW(volddummy,double,mt**nk);
646 NNEW(veolddummy,double,mt**nk);
647 NNEW(accolddummy,double,mt**nk);
648 for(k=0;k<mt**nk;k++){
649 volddummy[k]=0.0;
650 veolddummy[k]=0.0;
651 accolddummy[k]=0.0;}
652 memcpy(&v[0],&vold[0],sizeof(double)*mt**nk);
653 NNEW(vectornull,double,neq[1]);
654 *iout=-1;
655
656 FORTRAN(resultsini_mortar,(nk,v,ithermal,iperturb,
657 nactdof,iout,volddummy,vectornull,nodeboun,ndirboun,
658 xbounact,nboun,ipompc,nodempc,coefmpc,labmpc,nmpc,nmethod,cam,
659 bet,gam,dtime,mi));
660
661 memcpy(&vold[0],&v[0],sizeof(double)*mt**nk);
662 SFREE(v);SFREE(vectornull);SFREE(volddummy);SFREE(veolddummy);
663 SFREE(accolddummy);
664 }
665
666 *ielas=0;
667 *iout=0;
668 iperturb[0]=iperturb_sav[0];
669 iperturb[1]=iperturb_sav[1];
670 mortartrafoflag=0;
671
672 *ndirboun2p=ndirboun2;*nodeboun2p=nodeboun2;*xboun2p=xboun2;
673 *ipompc2p=ipompc2;*nodempc2p=nodempc2;*coefmpc2p=coefmpc2;
674 *labmpc2p=labmpc2;*ikboun2p=ikboun2;*ilboun2p=ilboun2;*ikmpc2p=ikmpc2;
675 *ilmpc2p=ilmpc2;
676 *nslavspcp=nslavspc;*islavspcp=islavspc;*nslavmpcp=nslavmpc;
677 *islavmpcp=islavmpc;
678 *nslavspc2p=nslavspc2;*islavspc2p=islavspc2;*nslavmpc2p=nslavmpc2;
679 *islavmpc2p=islavmpc2;
680 *nmastspcp=nmastspc;*imastspcp=imastspc;*nmastmpcp=nmastmpc;
681 *nmastmpc2p=nmastmpc2;*imastmpcp=imastmpc;*imastmpc2p=imastmpc2;
682 *auc2p=auc2;*adc2p=adc2;*irowc2p=irowc2;*icolc2p=icolc2;*jqc2p=jqc2;
683 *aubdp=aubd;*irowbdp=irowbd;*jqbdp=jqbd;
684 *aubdtilp=aubdtil;*irowbdtilp=irowbdtil;*jqbdtilp=jqbdtil;
685 *aubdtil2p=aubdtil2;*irowbdtil2p=irowbdtil2;*jqbdtil2p=jqbdtil2;
686 *auddp=audd;*irowddp=irowdd;*jqddp=jqdd;
687 *auddtilp=auddtil;*irowddtilp=irowddtil;*jqddtilp=jqddtil;
688 *auddtil2p=auddtil2;*irowddtil2p=irowddtil2;*jqddtil2p=jqddtil2;
689 *auddinvp=auddinv;*irowddinvp=irowddinv;*jqddinvp=jqddinv;
690 *irowtempp=irowtemp;*icoltempp=icoltemp;*jqtempp=jqtemp;
691 *f_csp=f_cs;*f_cmp=f_cm;
692 *nodeforc2p=nodeforc2;*ndirforc2p=ndirforc2;
693 *xforc2p=xforc2;
694
695 return;
696 }
697