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 <stdlib.h>
19 #include <math.h>
20 #include <stdio.h>
21 #include <time.h>
22 #include <string.h>
23 #include "CalculiX.h"
24 #include "mortar.h"
25
26 #define max(a,b) (((a) > (b)) ? (a) : (b))
27 #define min(a,b) (((a) < (b)) ? (a) : (b))
28 /**
29 * embedding the contact conditions into the matrix system
30 * and condening the Lagrange multiplier
31 * see phd-thesis Sitzmann equation (4.15) for quad-quad/quad-lin method or (4.24) for PG quad-lin method
32
33 * Authors: Samoela Rakotonanahary, Saskia Sitzmann
34 *
35 * [out] aucp nondiagonal entries of intermediate matrix \f$ K_c \f$
36 * [out] adc diagonal entries of intermediate matrix \f$ K_c \f$
37 * [out] irowcp field containing row number for entries in auc
38 * [out] jqc (i) first element in auc belonging to column i
39 * [out] nzsc size of auc
40 * [in] aubd coupling matrix \f$ B_d[nactdof(i,p),nactdof(j,q)]\f$ for all active degrees od freedoms
41 * [in] irowbd field containing row numbers of aubd
42 * [in] jqbd pointer into field irowbd
43 * [in] aubdtil matrix \f$ \tilde{D}^{-1}\tilde{B}_d[nactdof(i,p),nactdof(j,q)]\f$ for all active degrees od freedoms
44 * [in] irowbdtil field containing row numbers of aubd
45 * [in] jqbdtil pointer into field irowbdtil
46 * [in] aubdtil2 coupling matrix \f$ \tilde{D}$ and $\tilde{B}^2_d[nactdof(i,p),nactdof(j,q)]\f$ for all active degrees od freedoms
47 * [in] irowbdtil2 field containing row numbers of aubdtil2
48 * [in] jqbdtil2 pointer into field irowbdtil2
49 * [in] irowdd field containing row numbers of audd
50 * [in] jqdd pointer into field irowdd
51 * [in] audd coupling matrix \f$ D_d[nactdof(i,p),nactdof(j,q)]\f$ for all active degrees od freedoms
52 * [in] irowddtil2 field containing row numbers of audd
53 * [in] jqddtil2 pointer into field irowdd
54 * [in] auddtil2 matrix \f$ Id_d[nactdof(i,p),nactdof(j,q)]\f$ for all active degrees od freedoms
55 * [in] irowddinv field containing row numbers of auddinv
56 * [in] jqddinv pointer into field irowddinv
57 * [in] auddinv coupling matrix \f$ \tilde{D}^{-1}_d[nactdof(i,p),nactdof(j,q)]\f$ for all active degrees od freedoms
58 * [in] Bd coupling matrix \f$ B_d[p,q]=\int \psi_p \phi_q dS \f$, \f$ p \in S, q \in M \f$
59 * [in] irowb field containing row numbers of Bd
60 * [in] jqb pointer into field irowb
61 * [in] Dd coupling matrix \f$ D_d[p,q]=\int \psi_p \phi_q dS \f$, \f$ p,q \in S \f$
62 * [in] irowd field containing row numbers of Dd
63 * [in] jqd pointer into field irowd
64 * [in] Ddtil coupling matrix \f$ \tilde{D}_d[p,q]=\int \psi_p \tilde{\phi}_q dS \f$, \f$ p,q \in S \f$
65 * [in] irowdtil field containing row numbers of Ddtil
66 * [in] jqdtil pointer into field irowdtil
67 * [in] neq (0) # of mechanical equations (1) sum of mechanical and thermal equations (2) neq(1+ # of single point contraints)
68 * [in,out] b right hand side
69 * [out] bhat intermediate right hand side
70 * [in] islavnode field storing the nodes of the slave surface
71 * [in] imastnode field storing the nodes of the master surfaces
72 * [in] nslavnode (i)pointer into field isalvnode for contact tie i
73 * [in] nmastnode (i)pointer into field imastnode for contact tie i
74 * [in] islavact (i) indicates, if slave node i is active (=-3 no-slave-node, =-2 no-LM-node, =-1 no-gap-node, =0 inactive node, =1 sticky node, =2 slipping/active node)
75 * [in] islavactdof (i)=10*slavenodenumber+direction for active dof i
76 * [in] gap (i) \f$ g_i= <g, \Psi_i> \f$ for node i on slave surface
77 * [in] slavnor slave normal
78 * [in] slavtan slave tangent
79 * [in] cstress current Lagrange multiplier
80 * [in] cstressini Lagrange multiplier at start of the increment
81 * [in] bp_old old friction bounds
82 * [in] nslavspc (2*i) pointer to islavspc...
83 * [in] islavspc ... which stores SPCs for slave node i
84 * [in] nsspc number of SPC for slave nodes
85 * [in] nslavmpc (2*i) pointer to islavmpc...
86 * [in] islavmpc ... which stores MPCs for slave node i
87 * [in] nsmpc number of MPC for slave nodes
88 * [in] nmastspc (2*i) pointer to imastspc...
89 * [in] imastspc ... which stores SPCs for master node i
90 * [in] nmspc number of SPC for master nodes
91 * [in] nmastmpc (2*i) pointer to imastmpc...
92 * [in] imastmpc ... which stores MPCs for master node i
93 * [in] nmmpc number of MPC for master nodes
94 * [in] tieset (1,i) name of tie constraint (2,i) dependent surface (3,i) independent surface
95 * [in] islavactdoftie (i)=tie number for active dof i
96 * [in] irowtloc field containing row numbers of autloc
97 * [in] jqtloc pointer into field irowtloc
98 * [in] autloc transformation matrix \f$ T[p,q]\f$ for slave nodes \f$ p,q \f$
99 * [in] irowtlocinv field containing row numbers of autlocinv
100 * [in] jqtlocinv pointer into field irowtlocinv
101 * [in] autlocinv transformation matrix \f$ T^{-1}[p,q]\f$ for slave nodes \f$ p,q \f$
102 * [in] islavnodeinv (i) slave node index for node i
103 * [in] lambdaiwan Lagrange multiplier splitted to Iwan elements
104 * [in] lambdaiwanini Lagrange multiplier splitted to Iwan elements at start of increment
105 */
106
multimortar(double ** aup,double * ad,ITG ** irowp,ITG * jq,ITG * nzs,double ** aucp,double * adc,ITG ** irowcp,ITG * jqc,ITG * nzsc,double * aubd,ITG * irowbd,ITG * jqbd,double * aubdtil,ITG * irowbdtil,ITG * jqbdtil,double * aubdtil2,ITG * irowbdtil2,ITG * jqbdtil2,ITG * irowdd,ITG * jqdd,double * audd,ITG * irowddtil2,ITG * jqddtil2,double * auddtil2,ITG * irowddinv,ITG * jqddinv,double * auddinv,double * Bd,ITG * irowb,ITG * jqb,double * Dd,ITG * irowd,ITG * jqd,double * Ddtil,ITG * irowdtil,ITG * jqdtil,ITG * neq,double * b,double * bhat,ITG * islavnode,ITG * imastnode,ITG * nslavnode,ITG * nmastnode,ITG * islavact,ITG * islavactdof,double * gap,double * slavnor,double * slavtan,double * vold,double * vini,double * cstress,double * cstressini,double * bp_old,ITG * nactdof,ITG * ntie,ITG * mi,ITG * nk,ITG * nboun,ITG * ndirboun,ITG * nodeboun,double * xboun,ITG * nmpc,ITG * ipompc,ITG * nodempc,double * coefmpc,ITG * ikboun,ITG * ilboun,ITG * ikmpc,ITG * ilmpc,ITG * nslavspc,ITG * islavspc,ITG * nsspc,ITG * nslavmpc,ITG * islavmpc,ITG * nsmpc,ITG * nmastspc,ITG * imastspc,ITG * nmspc,ITG * nmastmpc,ITG * imastmpc,ITG * nmmpc,char * tieset,ITG * islavactdoftie,ITG * nelcon,double * elcon,double * tietol,ITG * ncmat_,ITG * ntmat_,double * plicon,ITG * nplicon,ITG * npmat_,double * dtime,ITG * irowtloc,ITG * jqtloc,double * autloc,ITG * irowtlocinv,ITG * jqtlocinv,double * autlocinv,ITG * islavnodeinv,double * lambdaiwan,double * lambdaiwanini,ITG * iit,ITG * nmethod,double * bet,ITG * ithermal,double * plkcon,ITG * nplkcon)107 void multimortar(double **aup,double *ad,ITG **irowp,ITG *jq,ITG *nzs,
108 double **aucp,double *adc,ITG **irowcp,ITG *jqc,ITG *nzsc,
109 double *aubd,ITG *irowbd,ITG *jqbd,double *aubdtil,
110 ITG *irowbdtil,ITG *jqbdtil,double *aubdtil2,ITG *irowbdtil2,
111 ITG *jqbdtil2,ITG *irowdd,ITG *jqdd,double *audd,
112 ITG *irowddtil2,ITG *jqddtil2,double *auddtil2,
113 ITG *irowddinv,ITG *jqddinv,double *auddinv,double *Bd,
114 ITG *irowb,ITG *jqb,double *Dd,ITG *irowd,ITG *jqd,
115 double *Ddtil,ITG *irowdtil,ITG *jqdtil,ITG *neq,double *b,
116 double *bhat,ITG *islavnode,ITG *imastnode,
117 ITG *nslavnode,ITG *nmastnode,ITG *islavact,ITG *islavactdof,
118 double *gap,double *slavnor,double *slavtan,double *vold,
119 double *vini,double *cstress,double *cstressini,
120 double *bp_old,ITG *nactdof,ITG *ntie,ITG *mi,ITG *nk,
121 ITG *nboun,ITG *ndirboun,ITG *nodeboun,double *xboun,
122 ITG *nmpc,ITG *ipompc,ITG *nodempc,double *coefmpc,
123 ITG *ikboun,ITG *ilboun,ITG *ikmpc,ITG *ilmpc,ITG *nslavspc,
124 ITG *islavspc,ITG *nsspc,ITG *nslavmpc,ITG *islavmpc,
125 ITG *nsmpc,ITG *nmastspc,ITG *imastspc,ITG *nmspc,
126 ITG *nmastmpc,ITG *imastmpc,ITG *nmmpc,char *tieset,
127 ITG *islavactdoftie,ITG *nelcon,double *elcon,
128 double *tietol,ITG *ncmat_,ITG *ntmat_,double *plicon,
129 ITG *nplicon,ITG *npmat_,double *dtime,ITG *irowtloc,
130 ITG *jqtloc,double *autloc, ITG *irowtlocinv,ITG *jqtlocinv,
131 double *autlocinv,ITG *islavnodeinv,double *lambdaiwan,
132 double *lambdaiwanini,ITG *iit,ITG *nmethod,double *bet,
133 ITG *ithermal,double *plkcon,ITG *nplkcon){
134
135 ITG i,j,k,l,m,mt=mi[1]+1,nodesf,nodem,irow_ln,irow_lm,irow_li,irow_la,debug,
136 impclack,dim,nzs_nn,nzs_nm,nzs_ni,nzs_na,nzs_mn,nzs_mm,nzs_mi,nzs_ma,
137 nzs_in,nzs_im,nzs_ii,nzs_ia,nzs_an,nzs_am,nzs_ai,nzs_aa,nzs_bd1,nzs_bdtil2,
138 nzs_ddtil2i,nzs_ddtil2a,nzsbdtil2,nzs_mmf,nzs_iif,nzs_aaf,nzs_intmn,
139 nzs_intmm,nzs_intmi,nzs_intma,nzs_mntil,nzs_mmtil,nzs_mitil,nzs_matil,
140 ifree_nn,ifree_nm,ifree_ni,ifree_na,ifree_mn,ifree_mm,ifree_mi,ifree_ma,
141 ifree_in,ifree_im,ifree_ii,ifree_ia,ifree_an,ifree_am,ifree_ai,ifree_aa,
142 ifree,*irow=NULL,*irow_iid=NULL,*irow_nn=NULL,*irow_nm=NULL,*irow_ni=NULL,
143 *irow_na=NULL,*irow_mn=NULL,*irow_mm=NULL,*irow_mi=NULL,*irow_ma=NULL,
144 *irow_in=NULL,*irow_im=NULL,*irow_ii=NULL,*irow_ia=NULL,*irow_an=NULL,
145 *irow_am=NULL,*irow_ai=NULL,*irow_aa=NULL,*jq_iid=NULL,*jq_nn=NULL,
146 *jq_nm=NULL,*jq_ni=NULL,*jq_na=NULL,*jq_mn=NULL,*jq_mm=NULL,*jq_mi=NULL,
147 *jq_ma=NULL,*jq_in=NULL,*jq_im=NULL,*jq_ii=NULL,*jq_ia=NULL,*jq_an=NULL,
148 *jq_am=NULL,*jq_ai=NULL,*jq_aa=NULL,*irowc=NULL,*irow_bdtil2=NULL,
149 *jq_bdtil2=NULL,*irow_ddtil2i=NULL,*jq_ddtil2i=NULL,*irow_ddtil2a=NULL,
150 *jq_ddtil2a=NULL,*irow_bd1=NULL,*jq_bd1=NULL,*irow_aad=NULL,*irow_mmd=NULL,
151 *jq_aad=NULL,*jq_mmd=NULL,*irow_mmf=NULL,*jq_mmf=NULL,*irow_mmtil=NULL,
152 *jq_mmtil=NULL,*irow_mntil=NULL,*jq_mntil=NULL,*irow_mitil=NULL,
153 *jq_mitil=NULL,*irow_matil=NULL,*jq_matil=NULL,*irow_intmn=NULL,
154 *jq_intmn=NULL,*irow_intmm=NULL,*jq_intmm=NULL,*irow_intmi=NULL,
155 *jq_intmi=NULL,*irow_intma=NULL,*jq_intma=NULL,*irow_iif=NULL,*jq_iif=NULL,
156 *irow_aaf=NULL,*jq_aaf=NULL,*irow_antil=NULL,
157 *jq_antil=NULL,*irow_amtil=NULL,*jq_amtil=NULL,*irow_aitil=NULL,
158 *jq_aitil=NULL,*irow_aatil=NULL,*jq_aatil=NULL,*irow_t=NULL,*jq_t=NULL,
159 nzs_t=*nzs,*mast1=NULL,*l_flag=NULL,*n_flag=NULL,*m_flag=NULL,*a_flag=NULL,
160 *i_flag=NULL,number=1,iact,*n_flagr=NULL,*m_flagr=NULL,*i_flagr=NULL,
161 *a_flagr=NULL,nzsddtil2a,nzsddtil2i,nzs_sym;
162
163 double
164 *auc=NULL,*au=NULL,*ad_nn=NULL,*ad_mm=NULL,*ad_ii=NULL,*ad_aa,*au_bd1=NULL,
165 *au_bdtil2=NULL,*au_ddtil2i=NULL,*au_ddtil2a=NULL,*au_nn=NULL,*au_nm=NULL,
166 *au_ni=NULL,*au_na=NULL,*au_mn=NULL,*au_mm=NULL,*au_mi=NULL,*au_ma=NULL,
167 *au_in=NULL,*au_im=NULL,*au_ii=NULL,*au_ia=NULL,*au_an=NULL,*au_am=NULL,
168 *au_ai=NULL,*au_aa=NULL,*au_mmf=NULL,*au_mmtil=NULL,*au_mntil=NULL,
169 *au_mitil=NULL,*au_matil=NULL,*au_aaf=NULL,*au_iif=NULL,*au_intmn=NULL,
170 *au_intmm=NULL,*au_intmi=NULL,*au_intma=NULL,*au_antil=NULL,
171 *au_amtil=NULL,*au_aitil=NULL,*au_aatil=NULL,*au_t=NULL,*f_a=NULL,
172 *f_m=NULL,*f_i=NULL,*v_r=NULL,*f_atil=NULL;
173
174 irowc=*irowcp; auc=*aucp;
175 irow=*irowp; au=*aup;
176 debug=0;
177
178 /* save full K matrix without contact conditions
179 needed for calculation of the Lagrange multiplier in stressmortar */
180
181 /* the symmetric storage of au is
182 converted into a asymmetric storage of the lower (au) and upper (au^T)
183 triangle in auc */
184
185 NNEW(au_t,double,nzs_t);
186 NNEW(irow_t,ITG,nzs_t);
187 NNEW(jq_t,ITG,neq[1]+1);
188 dim=neq[1];
189
190 transpose(au,jq,irow,&dim,
191 au_t,jq_t,irow_t);
192
193 add_rect(au,irow,jq,neq[1],neq[1],
194 au_t,irow_t,jq_t,neq[1],neq[1],
195 &auc,&irowc,jqc,nzsc);
196
197 /* guido : adc not further used??? */
198
199 for(i=0;i<neq[1];i++){
200 adc[i]=ad[i];}
201 *nzsc=jqc[neq[1]]-1;
202 RENEW(auc,double,*nzsc);
203 RENEW(irowc,ITG,*nzsc);
204
205 /* Flag to produce the bijection between local and global dof **/
206 /*
207 * l_flag[i]=1 for M Master dof
208 * l_flag[i]=2 for I Slave dof
209 * l_flag[i]=3 for A Slave dof
210 * l_flag[i]=0 for N rest of the dof
211 *
212 * n_flag contains local N_row number
213 * m_flag contains local M_row number
214 * i_flag contains local I_row number
215 * a_flag contains local A_row number
216 **/
217
218 NNEW(l_flag,ITG,neq[1]);
219 NNEW(n_flag,ITG,neq[1]);
220 NNEW(m_flag,ITG,neq[1]);
221 NNEW(i_flag,ITG,neq[1]);
222 NNEW(a_flag,ITG,neq[1]);
223
224 /* Fill l_flag*/
225 //Slave
226
227 for(i=0;i<*ntie;i++){
228 if(tieset[i*(81*3)+80]=='C'){
229 for(j=nslavnode[i];j<nslavnode[i+1];j++){
230 nodesf=islavnode[j];
231
232 /* right now temperature degrees of freedom are set to N,
233 since l_flag[k-1]==0 */
234
235 for(l=0;l<3;l++){
236 //test of dof
237 k=nactdof[mt*(nodesf)-3+l];
238 if(k>0&&islavact[j]==0){
239 l_flag[k-1]=2;}
240 if(k>0&&islavact[j]>0) {
241 l_flag[k-1]=3;
242 iact++;}
243
244 // no-gap and no-LM nodes must be treated as master nodes
245
246 if(k>0&&(islavact[j]==-1 ||islavact[j]==-2)){
247 l_flag[k-1]=1;}
248
249 // nodes with no integration points can be treated N-nodes
250
251 if(k>0&&islavact[j]==-3){
252 l_flag[k-1]=0;}
253 }
254 if (ithermal[0]==3){
255 //coupled thermo-mechanical calculation
256 k=nactdof[mt*(nodesf)-4];
257 if(k>0&&islavact[j]==0){
258 l_flag[k-1]=2;}
259 if(k>0&&islavact[j]>0) {
260 l_flag[k-1]=3;
261 iact++;}
262
263 // no-gap and no-LM nodes must be treated as master nodes
264
265 if(k>0&&(islavact[j]==-1 ||islavact[j]==-2)){
266 l_flag[k-1]=1;}
267
268 // nodes with no integration points can be treated N-nodes
269
270 if(k>0&&islavact[j]==-3){
271 l_flag[k-1]=0;}
272 }
273 }
274 }
275 }
276
277 //Master
278
279 for(i=0;i<*ntie;i++){
280 if(tieset[i*(81*3)+80]=='C'){
281 for(j=nmastnode[i];j<nmastnode[i+1];j++){
282 nodem=imastnode[j];
283 for(l=0;l<3;l++){
284 //test of dof
285 k=nactdof[mt*(nodem)-3+l];
286 if(k>0)
287 {if(jqbdtil[k]-jqbdtil[k-1]>0){
288 l_flag[k-1]=1;
289 }else{
290 l_flag[k-1]=0;
291 }
292 }
293 }
294 if (ithermal[0]==3){
295 //coupled thermo-mechanical calculation
296 k=nactdof[mt*(nodem)-4];
297 if(k>0){
298 if(jqbdtil[k]-jqbdtil[k-1]>0){
299 l_flag[k-1]=1;
300 }else{
301 l_flag[k-1]=0;
302 }
303 }
304 }
305 }
306 }
307 }
308
309 /*** Fill the local row ***/
310
311 irow_ln=0;
312 irow_lm=0;
313 irow_la=0;
314 irow_li=0;
315
316 /* Storage of the diagonal of auc and the rhs */
317
318 /* ad_nn just contains the diagonal in a single vector
319 ad_mm, jq_mm, irow_mm: storage of ad_mm as general matrix
320 ad_aa, jq_aa, irow_aa: storage of ad_aa as general matrix
321 ad_ii, jq_ii, irow_ii: storage of ad_mm as general matrix
322 storage of a diagonal matrix ad in the form of a general matrix
323 does not change ad, it only requires the additional calculation
324 of vectors jq and irow */
325
326 /* storage as a general matrix allows for the usage of add_rect and
327 multi_rect for the purpose of adding and multiplying the matrices
328 in sparse matrix format */
329
330 NNEW(ad_nn,double,neq[1]);
331 NNEW(ad_mm,double,neq[1]);
332 NNEW(irow_mmd,ITG,neq[1]);
333 NNEW(jq_mmd,ITG,neq[1]);
334 NNEW(ad_aa,double,neq[1]);
335 NNEW(irow_aad,ITG,neq[1]);
336 NNEW(jq_aad,ITG,neq[1]);
337 NNEW(ad_ii,double,neq[1]);
338 NNEW(irow_iid,ITG,neq[1]);
339 NNEW(jq_iid,ITG,neq[1]);
340
341 NNEW(f_a,double,neq[1]);
342 NNEW(f_i,double,neq[1]);
343 NNEW(f_m,double,neq[1]);
344
345 for(i=0;i<neq[1];i++){
346 bhat[i]=b[i];
347 switch(l_flag[i]){
348 case 0 : {
349 //N
350 ad_nn[irow_ln]=ad[i];
351 n_flag[i]=(++irow_ln);
352 break;
353 }
354 case 1 : {
355 //M
356 f_m[irow_lm]=b[i];
357 m_flag[i]=(++irow_lm);
358 ad_mm[irow_lm-1]=ad[i];
359 jq_mmd[irow_lm-1]=irow_lm;
360 irow_mmd[irow_lm-1]=irow_lm;
361 break;
362 }
363 case 2 : {
364 //I
365 ad_ii[irow_li]=ad[i];
366 f_i[irow_li]=b[i];
367 i_flag[i]=(++irow_li);
368 jq_iid[irow_li-1]=irow_li;
369 irow_iid[irow_li-1]=irow_li;
370 break;
371 }
372 case 3 : {
373 //A
374 f_a[irow_la]=b[i];
375 a_flag[i]=(++irow_la);
376 ad_aa[irow_la-1]=ad[i];
377 jq_aad[irow_la-1]=irow_la;
378 irow_aad[irow_la-1]=irow_la;
379 break;
380 }
381 default : {
382 printf("*ERROR in multimortar: error l_flag\n");
383 break;
384 }
385 }
386 }
387 jq_iid[irow_li]=irow_li+1;
388 jq_aad[irow_la]=irow_la+1;
389 jq_mmd[irow_lm]=irow_lm+1;
390 RENEW(f_a,double,irow_la);
391 RENEW(f_i,double,irow_li);
392 RENEW(f_m,double,irow_lm);
393 RENEW(ad_nn,double,irow_ln);
394 RENEW(ad_mm,double,irow_lm);
395 RENEW(irow_mmd,ITG,irow_lm);
396 RENEW(jq_mmd,ITG,irow_lm+1);
397 RENEW(ad_aa,double,irow_la);
398 RENEW(irow_aad,ITG,irow_la);
399 RENEW(jq_aad,ITG,irow_la+1);
400 RENEW(ad_ii,double,irow_li);
401 RENEW(irow_iid,ITG,irow_li);
402 RENEW(jq_iid,ITG,irow_li+1);
403
404 /* end of the storage of the diagonals and the rhs */
405
406 /* extracting the submatrices **/
407
408 NNEW(au_nn,double,*nzs/8);
409 NNEW(au_nm,double,*nzs/8);
410 NNEW(au_ni,double,*nzs/8);
411 NNEW(au_na,double,*nzs/8);
412 NNEW(au_mn,double,*nzs/8);
413 NNEW(au_mm,double,*nzs/8);
414 NNEW(au_mi,double,*nzs/8);
415 NNEW(au_ma,double,*nzs/8);
416 NNEW(au_in,double,*nzs/8);
417 NNEW(au_im,double,*nzs/8);
418 NNEW(au_ii,double,*nzs/8);
419 NNEW(au_ia,double,*nzs/8);
420 NNEW(au_an,double,*nzs/8);
421 NNEW(au_am,double,*nzs/8);
422 NNEW(au_ai,double,*nzs/8);
423 NNEW(au_aa,double,*nzs/8);
424
425 NNEW(irow_nm,ITG,*nzs/8);
426 NNEW(irow_nn,ITG,*nzs/8);
427 NNEW(irow_ni,ITG,*nzs/8);
428 NNEW(irow_na,ITG,*nzs/8);
429 NNEW(irow_mm,ITG,*nzs/8);
430 NNEW(irow_mn,ITG,*nzs/8);
431 NNEW(irow_mi,ITG,*nzs/8);
432 NNEW(irow_ma,ITG,*nzs/8);
433 NNEW(irow_im,ITG,*nzs/8);
434 NNEW(irow_in,ITG,*nzs/8);
435 NNEW(irow_ii,ITG,*nzs/8);
436 NNEW(irow_ia,ITG,*nzs/8);
437 NNEW(irow_an,ITG,*nzs/8);
438 NNEW(irow_am,ITG,*nzs/8);
439 NNEW(irow_ai,ITG,*nzs/8);
440 NNEW(irow_aa,ITG,*nzs/8);
441
442 NNEW(jq_nn,ITG,irow_ln+1);
443 NNEW(jq_nm,ITG,irow_lm+1);
444 NNEW(jq_ni,ITG,irow_li+1);
445 NNEW(jq_na,ITG,irow_la+1);
446 NNEW(jq_mn,ITG,irow_ln+1);
447 NNEW(jq_mm,ITG,irow_lm+1);
448 NNEW(jq_mi,ITG,irow_li+1);
449 NNEW(jq_ma,ITG,irow_la+1);
450 NNEW(jq_in,ITG,irow_ln+1);
451 NNEW(jq_im,ITG,irow_lm+1);
452 NNEW(jq_ii,ITG,irow_li+1);
453 NNEW(jq_ia,ITG,irow_la+1);
454 NNEW(jq_an,ITG,irow_ln+1);
455 NNEW(jq_am,ITG,irow_lm+1);
456 NNEW(jq_ai,ITG,irow_li+1);
457 NNEW(jq_aa,ITG,irow_la+1);
458
459 jq_nn[0]=1;
460 jq_nm[0]=1;
461 jq_ni[0]=1;
462 jq_na[0]=1;
463 jq_mn[0]=1;
464 jq_mm[0]=1;
465 jq_mi[0]=1;
466 jq_ma[0]=1;
467 jq_in[0]=1;
468 jq_im[0]=1;
469 jq_ii[0]=1;
470 jq_ia[0]=1;
471 jq_an[0]=1;
472 jq_am[0]=1;
473 jq_ai[0]=1;
474 jq_aa[0]=1;
475
476 nzs_nn=*nzs/8.0;
477 nzs_nm=*nzs/8.0;
478 nzs_ni=*nzs/8.0;
479 nzs_na=*nzs/8.0;
480 nzs_mn=*nzs/8.0;
481 nzs_mm=*nzs/8.0;
482 nzs_mi=*nzs/8.0;
483 nzs_ma=*nzs/8.0;
484 nzs_in=*nzs/8.0;
485 nzs_im=*nzs/8.0;
486 nzs_ii=*nzs/8.0;
487 nzs_ia=*nzs/8.0;
488 nzs_an=*nzs/8.0;
489 nzs_am=*nzs/8.0;
490 nzs_ai=*nzs/8.0;
491 nzs_aa=*nzs/8.0;
492
493 ifree_nn=1;
494 ifree_nm=1;
495 ifree_ni=1;
496 ifree_na=1;
497 ifree_mn=1;
498 ifree_mm=1;
499 ifree_mi=1;
500 ifree_ma=1;
501 ifree_in=1;
502 ifree_im=1;
503 ifree_ii=1;
504 ifree_ia=1;
505 ifree_an=1;
506 ifree_am=1;
507 ifree_ai=1;
508 ifree_aa=1;
509
510 for(j=0;j<neq[1];j++){
511 m=j+1;
512 for(i=jqc[j]-1;i<jqc[j+1]-1;i++){
513 switch(l_flag[j]){
514 case 0 : {
515 // matrices XN
516 m=n_flag[j]+1;
517 switch(l_flag[irowc[i]-1]){
518 case 0 : {
519 //NN
520 insertas_ws(&irow_nn,&(n_flag[irowc[i]-1]),&m,&ifree_nn,
521 &nzs_nn,&auc[i],&au_nn);
522 break;
523 }
524 case 1 :{
525 //MN
526 insertas_ws(&irow_mn,&(m_flag[irowc[i]-1]),&m,&ifree_mn,
527 &nzs_mn,&auc[i],&au_mn);
528 break; }
529 case 2 :{
530 //IN
531 insertas_ws(&irow_in,&(i_flag[irowc[i]-1]),&m,&ifree_in,
532 &nzs_in,&auc[i],&au_in);
533 break;}
534
535 case 3:{
536 //AN
537 insertas_ws(&irow_an,&(a_flag[irowc[i]-1]),&m,&ifree_an,
538 &nzs_an,&auc[i],&au_an);
539 break;}
540 default :{
541 break; }
542 }
543 break;
544 }
545 case 1 : {
546 // matrices XM
547 m=m_flag[j]+1;
548 switch(l_flag[irowc[i]-1]){
549 case 0 : {
550 //NM
551 insertas_ws(&irow_nm,&(n_flag[irowc[i]-1]),&m,&ifree_nm,
552 &nzs_nm,&auc[i],&au_nm);
553 break;
554 }
555 case 1 : {
556 //MM
557 insertas_ws(&irow_mm,&(m_flag[irowc[i]-1]),&m,&ifree_mm,
558 &nzs_mm,&auc[i],&au_mm);
559 break;
560 }
561 case 2 :{
562 //IM
563 insertas_ws(&irow_im,&(i_flag[irowc[i]-1]),&m,&ifree_im,
564 &nzs_im,&auc[i],&au_im);
565 break;}
566 case 3 : {
567 //AM
568 insertas_ws(&irow_am,&(a_flag[irowc[i]-1]),&m,&ifree_am,
569 &nzs_am,&auc[i],&au_am);
570 break;
571 }
572 default : {
573 break;}
574 }
575 break;
576 }
577 case 2 : {
578 // matrices XI
579 m=i_flag[j]+1;
580 switch(l_flag[irowc[i]-1]){
581 case 0 : {
582 //NI
583 insertas_ws(&irow_ni,&(n_flag[irowc[i]-1]),&m,&ifree_ni,
584 &nzs_ni,&auc[i],&au_ni);
585 break;
586 }
587 case 1 : {
588 //MI
589 insertas_ws(&irow_mi,&(m_flag[irowc[i]-1]),&m,&ifree_mi,
590 &nzs_mi,&auc[i],&au_mi);
591 break;
592 }
593 case 2 :{
594 //II
595 insertas_ws(&irow_ii,&(i_flag[irowc[i]-1]),&m,&ifree_ii,
596 &nzs_ii,&auc[i],&au_ii);
597 break;}
598 case 3 : {
599 //AI
600 insertas_ws(&irow_ai,&(a_flag[irowc[i]-1]),&m,&ifree_ai,
601 &nzs_ai,&auc[i],&au_ai);
602 break;
603 }
604 default :
605 break;
606 }
607 break;
608 }
609 case 3 : {
610 // matrices XA
611 m=a_flag[j]+1;
612 switch(l_flag[irowc[i]-1]){
613 case 0 : {
614 //NA
615 insertas_ws(&irow_na,&(n_flag[irowc[i]-1]),&m,&ifree_na,
616 &nzs_na,&auc[i],&au_na);
617 break;
618 }
619 case 1 : {
620 //MA
621 insertas_ws(&irow_ma,&(m_flag[irowc[i]-1]),&m,&ifree_ma,
622 &nzs_ma,&auc[i],&au_ma);
623 break;
624 }
625 case 2 :{
626 //IA
627 insertas_ws(&irow_ia,&(i_flag[irowc[i]-1]),&m,&ifree_ia,
628 &nzs_ia,&auc[i],&au_ia);
629 break;}
630 case 3 : {
631 //AA
632 insertas_ws(&irow_aa,&(a_flag[irowc[i]-1]),&m,&ifree_aa,
633 &nzs_aa,&auc[i],&au_aa);
634 break;
635 }
636 default : {
637 break; }
638 }
639 break;
640 }
641 default :
642 break;
643 }
644 // end switch column lflag
645 }
646 // end loop nzs column
647
648 /* actualisation of the jq **/
649
650 if (n_flag[j]!=0) jq_nn[n_flag[j]]=ifree_nn;
651 if (n_flag[j]!=0)jq_mn[n_flag[j]]=ifree_mn;
652 if (n_flag[j]!=0)jq_in[n_flag[j]]=ifree_in;
653 if (n_flag[j]!=0)jq_an[n_flag[j]]=ifree_an;
654 if (m_flag[j]!=0) jq_nm[m_flag[j]]=ifree_nm;
655 if (m_flag[j]!=0)jq_mm[m_flag[j]]=ifree_mm;
656 if (m_flag[j]!=0)jq_im[m_flag[j]]=ifree_im;
657 if (m_flag[j]!=0)jq_am[m_flag[j]]=ifree_am;
658 if (i_flag[j]!=0) jq_ni[i_flag[j]]=ifree_ni;
659 if (i_flag[j]!=0)jq_mi[i_flag[j]]=ifree_mi;
660 if (i_flag[j]!=0)jq_ii[i_flag[j]]=ifree_ii;
661 if (i_flag[j]!=0)jq_ai[i_flag[j]]=ifree_ai;
662 if (a_flag[j]!=0) jq_na[a_flag[j]]=ifree_na;
663 if (a_flag[j]!=0)jq_ma[a_flag[j]]=ifree_ma;
664 if (a_flag[j]!=0)jq_ia[a_flag[j]]=ifree_ia;
665 if (a_flag[j]!=0)jq_aa[a_flag[j]]=ifree_aa;
666 }
667 // end loop over the global column
668
669 if(debug==1)printf("\tmm2: N %" ITGFORMAT " M %" ITGFORMAT " I %" ITGFORMAT
670 " A %" ITGFORMAT " nzs %" ITGFORMAT " \n",irow_ln,irow_lm,irow_li,
671 irow_la,*nzs);
672
673 /* NN */
674
675 if ((ifree_nn-1)!=0){
676 nzs_nn=ifree_nn-1;
677 RENEW(au_nn,double,nzs_nn);
678 RENEW(irow_nn,ITG,nzs_nn);
679 }else{
680 nzs_nn=ifree_nn-1;
681 RENEW(au_nn,double,1);
682 RENEW(irow_nn,ITG,1);
683 }
684
685 /* NM */
686
687 if ((ifree_nm-1)!=0){
688 nzs_nm=ifree_nm-1;
689 RENEW(au_nm,double,nzs_nm);
690 RENEW(irow_nm,ITG,nzs_nm);
691 }else{
692 nzs_nm=ifree_nm-1;
693 RENEW(au_nm,double,1);
694 RENEW(irow_nm,ITG,1);
695 }
696
697 /* NI */
698
699 if ((ifree_ni-1)!=0){
700 nzs_ni=ifree_ni-1;
701 RENEW(au_ni,double,nzs_ni);
702 RENEW(irow_ni,ITG,nzs_ni);
703 }else{
704 nzs_ni=ifree_ni-1;
705 RENEW(au_ni,double,1);
706 RENEW(irow_ni,ITG,1);
707 }
708
709 /* NA */
710
711 if ((ifree_na-1)!=0){
712 nzs_na=ifree_na-1;
713 RENEW(au_na,double,nzs_na);
714 RENEW(irow_na,ITG,nzs_na);
715 }else{
716 nzs_na=ifree_na-1;
717 RENEW(au_na,double,1);
718 RENEW(irow_na,ITG,1);
719 }
720
721 /* MN */
722
723 if ((ifree_mn-1)!=0){
724 nzs_mn=ifree_mn-1;
725 RENEW(au_mn,double,nzs_mn);
726 RENEW(irow_mn,ITG,nzs_mn);
727 }else{
728 nzs_mn=ifree_mn-1;
729 RENEW(au_mn,double,1);
730 RENEW(irow_mn,ITG,1);
731 }
732
733 /* MM */
734
735 if ((ifree_mm-1)!=0){
736 nzs_mm=ifree_mm-1;
737 RENEW(au_mm,double,nzs_mm);
738 RENEW(irow_mm,ITG,nzs_mm);
739 }else{
740 nzs_mm=ifree_mm-1;
741 RENEW(au_mm,double,1);
742 RENEW(irow_mm,ITG,1);
743 }
744
745 /* MI */
746
747 if ((ifree_mi-1)!=0){
748 nzs_mi=ifree_mi-1;
749 RENEW(au_mi,double,nzs_mi);
750 RENEW(irow_mi,ITG,nzs_mi);
751 }else{
752 nzs_mi=ifree_mi-1;
753 RENEW(au_mi,double,1);
754 RENEW(irow_mi,ITG,1);
755 }
756
757 /* MA */
758
759 if ((ifree_ma-1)!=0){
760 nzs_ma=ifree_ma-1;
761 RENEW(au_ma,double,nzs_ma);
762 RENEW(irow_ma,ITG,nzs_ma);
763 }else{
764 nzs_ma=ifree_ma-1;
765 RENEW(au_ma,double,1);
766 RENEW(irow_ma,ITG,1);
767 }
768
769 /* IN */
770
771 if ((ifree_in-1)!=0){
772 nzs_in=ifree_in-1;
773 RENEW(au_in,double,nzs_in);
774 RENEW(irow_in,ITG,nzs_in);
775 }else{
776 nzs_in=ifree_in-1;
777 RENEW(au_in,double,1);
778 RENEW(irow_in,ITG,1);
779 }
780
781 /* IM */
782
783 if ((ifree_im-1)!=0){
784 nzs_im=ifree_im-1;
785 RENEW(au_im,double,nzs_im);
786 RENEW(irow_im,ITG,nzs_im);
787 }else{
788 nzs_im=ifree_im-1;
789 RENEW(au_im,double,1);
790 RENEW(irow_im,ITG,1);
791 }
792
793 /* II */
794
795 if ((ifree_ii-1)!=0){
796 nzs_ii=ifree_ii-1;
797 RENEW(au_ii,double,nzs_ii);
798 RENEW(irow_ii,ITG,nzs_ii);
799 }else{
800 nzs_ii=ifree_ii-1;
801 RENEW(au_ii,double,1);
802 RENEW(irow_ii,ITG,1);
803 }
804
805 /* IA */
806
807 if ((ifree_ia-1)!=0){
808 nzs_ia=ifree_ia-1;
809 RENEW(au_ia,double,nzs_ia);
810 RENEW(irow_ia,ITG,nzs_ia);
811 }else{
812 nzs_ia=ifree_ia-1;
813 RENEW(au_ia,double,1);
814 RENEW(irow_ia,ITG,1);
815 }
816
817 /* AN */
818
819 if ((ifree_an-1)!=0){
820 nzs_an=ifree_an-1;
821 RENEW(au_an,double,nzs_an);
822 RENEW(irow_an,ITG,nzs_an);
823 }else{
824 nzs_an=ifree_an-1;
825 RENEW(au_an,double,1);
826 RENEW(irow_an,ITG,1);
827 }
828
829 /* AM */
830
831 if ((ifree_am-1)!=0){
832 nzs_am=ifree_am-1;
833 RENEW(au_am,double,nzs_am);
834 RENEW(irow_am,ITG,nzs_am);
835 }else{
836 nzs_am=ifree_am-1;
837 RENEW(au_am,double,1);
838 RENEW(irow_am,ITG,1);
839 }
840
841 /* AI */
842
843 if ((ifree_ai-1)!=0){
844 nzs_ai=ifree_ai-1;
845 RENEW(au_ai,double,nzs_ai);
846 RENEW(irow_ai,ITG,nzs_ai);
847 }else{
848 nzs_ai=ifree_ai-1;
849 RENEW(au_ai,double,1);
850 RENEW(irow_ai,ITG,1);
851 }
852
853 /* AA */
854
855 if ((ifree_aa-1)!=0){
856 nzs_aa=ifree_aa-1;
857 RENEW(au_aa,double,nzs_aa);
858 RENEW(irow_aa,ITG,nzs_aa);
859 }else{
860 nzs_aa=ifree_aa-1;
861 RENEW(au_aa,double,1);
862 RENEW(irow_aa,ITG,1);
863 }
864
865 if(debug==1){
866 printf("\tmm2: MN %" ITGFORMAT " %" ITGFORMAT " IN %" ITGFORMAT " %"
867 ITGFORMAT " AN %" ITGFORMAT " %" ITGFORMAT " \n",nzs_mn,nzs_nm,
868 nzs_in,nzs_ni,nzs_an,nzs_na);
869 printf("\tmm2: IM %" ITGFORMAT " %" ITGFORMAT " AM %" ITGFORMAT " %"
870 ITGFORMAT " AI %" ITGFORMAT " %" ITGFORMAT " \n",nzs_im,nzs_mi,
871 nzs_am,nzs_ma,nzs_ai,nzs_ia);
872 }
873
874 /* generate D_d,D_d^-1,B_a ,B_d^til in local dofs **/
875
876 nzsbdtil2=jqbdtil[neq[1]]-1;
877 NNEW(au_bd1,double,nzsbdtil2);
878 NNEW(irow_bd1,ITG,nzsbdtil2);
879 NNEW(jq_bd1,ITG,irow_lm+1);
880 nzs_bd1=0;
881 jq_bd1[0]=1;
882
883 nzsbdtil2=jqbdtil2[neq[1]]-1;
884 NNEW(au_bdtil2,double,nzsbdtil2);
885 NNEW(irow_bdtil2,ITG,nzsbdtil2);
886 NNEW(jq_bdtil2,ITG,irow_lm+1);
887 nzs_bdtil2=0;
888 jq_bdtil2[0]=1;
889 impclack=0;
890
891 nzsddtil2a=jqbdtil2[neq[1]]-1;
892 NNEW(au_ddtil2a,double,nzsddtil2a);
893 NNEW(irow_ddtil2a,ITG,nzsddtil2a);
894 NNEW(jq_ddtil2a,ITG,irow_la+1);
895 nzs_ddtil2a=0;
896 jq_ddtil2a[0]=1;
897
898 nzsddtil2i=jqbdtil2[neq[1]]-1;
899 NNEW(au_ddtil2i,double,nzsddtil2i);
900 NNEW(irow_ddtil2i,ITG,nzsddtil2i);
901 NNEW(jq_ddtil2i,ITG,irow_li+1);
902 nzs_ddtil2i=0;
903 jq_ddtil2i[0]=1;
904
905 for(j=0;j<neq[1];j++){
906 for(i=jqbdtil2[j]-1;i<jqbdtil2[j+1]-1;i++){
907 switch(l_flag[j]){
908 case 1 :{
909 //Matrix XM (column is master, row is whatever)
910 switch(l_flag[irowbdtil2[i]-1]){
911 case 3 :{
912 //AM here Bdtild
913 irow_bdtil2[nzs_bdtil2]=a_flag[irowbdtil2[i]-1];
914 au_bdtil2[nzs_bdtil2]=aubdtil2[i];
915 // Needed for slave part
916 nzs_bdtil2++;
917 break;}
918 default : {
919 break;}
920 }
921 break;}
922 case 0 :{
923 impclack++;
924 break; }
925 case 3 :{
926 break;}
927 default :{
928 break;}
929 }
930 }
931 if (m_flag[j]!=0) jq_bdtil2[m_flag[j]]=nzs_bdtil2+1;
932 }
933 impclack=0;
934 for(j=0;j<neq[1];j++){
935 for(i=jqbdtil2[j]-1;i<jqbdtil2[j+1]-1;i++){
936 switch(l_flag[j]){
937 case 3 :{
938 //Matrix XA (column is active slave, row is whatever)
939 switch(l_flag[irowbdtil2[i]-1]){
940 case 3 :{
941 //AA here Ddtild
942 irow_ddtil2a[nzs_ddtil2a]=a_flag[irowbdtil2[i]-1];
943 au_ddtil2a[nzs_ddtil2a]=aubdtil2[i];
944 // Needed for slave part
945 nzs_ddtil2a++;
946 break;}
947 default : {
948 break;}
949 }
950 break;}
951 case 0 :{
952 impclack++;
953 break; }
954 default :{
955 break; }
956 }
957 }
958 if (a_flag[j]!=0) jq_ddtil2a[a_flag[j]]=nzs_ddtil2a+1;
959 }
960 impclack=0;
961 for(j=0;j<neq[1];j++){
962 for(i=jqbdtil2[j]-1;i<jqbdtil2[j+1]-1;i++){
963 switch(l_flag[j]){
964 case 2 :{
965 //Matrix XI (column is inactive slave, row is whatever)
966 switch(l_flag[irowbdtil2[i]-1]){
967 case 3 :{
968 //AI here Ddtild
969 irow_ddtil2i[nzs_ddtil2i]=a_flag[irowbdtil2[i]-1];
970 au_ddtil2i[nzs_ddtil2i]=aubdtil2[i];
971 // Needed for slave part
972 nzs_ddtil2i++;
973 break;}
974 default : {
975 break;}
976 }
977 break;}
978 case 0 :{
979 impclack++;
980 break;}
981 default :{
982 break; }
983 }
984 }
985 if (i_flag[j]!=0) jq_ddtil2i[i_flag[j]]=nzs_ddtil2i+1;
986 }
987
988 /* \tilde{D}^{-1}\tilde{B}_d */
989
990 for(j=0;j<neq[1];j++){
991 for(i=jqbdtil[j]-1;i<jqbdtil[j+1]-1;i++){
992 switch(l_flag[j]){
993 case 1 :{
994 //Matrix XM (column is master, row is whatever)
995 switch(l_flag[irowbdtil[i]-1]){
996 case 3 :{
997 //AM here Bdtild
998 irow_bd1[nzs_bd1]=a_flag[irowbdtil[i]-1];
999 au_bd1[nzs_bd1]=-aubdtil[i];
1000 // Needed for master part
1001 nzs_bd1++;
1002 break;}
1003 default : {
1004 break;}
1005 }
1006 break;}
1007 default : {
1008 break;}
1009 }
1010 }
1011 if (m_flag[j]!=0) jq_bd1[m_flag[j]]=nzs_bd1+1;
1012 }
1013 RENEW(irow_bd1,ITG,nzs_bd1);
1014 RENEW(au_bd1,double,nzs_bd1);
1015 RENEW(irow_bdtil2,ITG,nzs_bdtil2);
1016 RENEW(au_bdtil2,double,nzs_bdtil2);
1017 RENEW(irow_ddtil2i,ITG,nzs_ddtil2i);
1018 RENEW(au_ddtil2i,double,nzs_ddtil2i);
1019 RENEW(irow_ddtil2a,ITG,nzs_ddtil2a);
1020 RENEW(au_ddtil2a,double,nzs_ddtil2a);
1021
1022 if(debug==1){
1023 printf("\tmm2: nzs_ddtil2a %" ITGFORMAT " mpclack %" ITGFORMAT " \n",
1024 nzs_ddtil2a,impclack);
1025 printf("\tmm2: nzs_ddtil2i %" ITGFORMAT " mpclack %" ITGFORMAT " \n",
1026 nzs_ddtil2i,impclack);
1027 printf("\tmm2: nzs_bdtil %" ITGFORMAT " mpclack %" ITGFORMAT " \n",
1028 nzs_bdtil2,impclack);
1029 printf("\tmm2: nzs_bd1 %" ITGFORMAT "\n",nzs_bd1);
1030
1031 /* add diagonals in K for matrices mm, aa and ii
1032 matrix nn is not transformed, so the diagonal of nn
1033 can be kept in vector ad_nn **/
1034
1035 printf("\tmm2: add diagonals\n");
1036 }
1037
1038 /* ad_mm: diagonal of mm stored as general matrix
1039 au_mm: off-diagonal terms of mm stored as general matrix
1040 au_mmf: all terms of mm stored as general matrix */
1041
1042 nzs_mmf=nzs_mm+irow_lm;
1043 NNEW(jq_mmf,ITG,irow_lm+1);
1044 NNEW(au_mmf,double,nzs_mmf);
1045 NNEW(irow_mmf,ITG,nzs_mmf);
1046 add_rect(ad_mm,irow_mmd,jq_mmd,irow_lm,irow_lm,
1047 au_mm,irow_mm,jq_mm,irow_lm,irow_lm,
1048 &au_mmf,&irow_mmf,jq_mmf,&nzs_mmf);
1049
1050 /* ad_ii: diagonal of ii stored as general matrix
1051 au_ii: off-diagonal terms of ii stored as general matrix
1052 au_iif: all terms of ii stored as general matrix */
1053
1054 nzs_iif=nzs_ii+irow_li;
1055 NNEW(jq_iif,ITG,irow_li+1);
1056 NNEW(au_iif,double,nzs_iif);
1057 NNEW(irow_iif,ITG,nzs_iif);
1058 add_rect(ad_ii,irow_iid,jq_iid,irow_li,irow_li,
1059 au_ii,irow_ii,jq_ii,irow_li,irow_li,
1060 &au_iif,&irow_iif,jq_iif,&nzs_iif);
1061
1062 /* ad_aa: diagonal of aa stored as general matrix
1063 au_aa: off-diagonal terms of aa stored as general matrix
1064 au_aaf: all terms of aa stored as general matrix */
1065
1066 nzs_aaf=nzs_aa+irow_la;
1067 NNEW(jq_aaf,ITG,irow_la+1);
1068 NNEW(au_aaf,double,nzs_aaf);
1069 NNEW(irow_aaf,ITG,nzs_aaf);
1070 add_rect(ad_aa,irow_aad,jq_aad,irow_la,irow_la,
1071 au_aa,irow_aa,jq_aa,irow_la,irow_la,
1072 &au_aaf,&irow_aaf,jq_aaf,&nzs_aaf);
1073 if(debug==1)printf("\tmm2: NN %" ITGFORMAT " MM %" ITGFORMAT " II %" ITGFORMAT " AA %"
1074 ITGFORMAT "\n",nzs_nn,nzs_mmf,nzs_iif,nzs_aaf);
1075 SFREE(au_ii);SFREE(irow_ii);SFREE(jq_ii);
1076 SFREE(au_aa);SFREE(irow_aa);SFREE(jq_aa);
1077 SFREE(au_mm);SFREE(irow_mm);SFREE(jq_mm);
1078
1079 /* Local => Global topology **/
1080
1081 NNEW(n_flagr,ITG,irow_ln);
1082 NNEW(m_flagr,ITG,irow_lm);
1083 NNEW(i_flagr,ITG,irow_li);
1084 NNEW(a_flagr,ITG,irow_la);
1085 for(j=0;j<neq[1];j++){
1086 switch(l_flag[j]){
1087 case 0 : {
1088 n_flagr[n_flag[j]-1]=j+1;
1089 break;}
1090 //N
1091 case 1 : {
1092 m_flagr[m_flag[j]-1]=j+1;
1093 break;}
1094 //M
1095 case 2 : {
1096 i_flagr[i_flag[j]-1]=j+1;
1097 break;}
1098 //I
1099 case 3 : {
1100 a_flagr[a_flag[j]-1]=j+1;
1101 break;}
1102 //A
1103 default : break;
1104 }
1105 }
1106
1107 /* K_MX -> K_MX^til: K_MX^til=K_MX-Btil^T*K_AX
1108 entries in the second row of the lhs matrix: **/
1109
1110 if(debug==1)printf("\tmm2: alter K_MX\n");
1111
1112 /* K_MN^til **/
1113
1114 nzs_intmn=nzs_mn+1;
1115 NNEW(jq_intmn,ITG,irow_ln+1);
1116 NNEW(au_intmn,double,nzs_intmn);
1117 NNEW(irow_intmn,ITG,nzs_intmn);
1118 multi_rect(au_bd1,irow_bd1,jq_bd1,irow_la,irow_lm,
1119 au_an,irow_an,jq_an,irow_la,irow_ln,
1120 &au_intmn,&irow_intmn,jq_intmn,&nzs_intmn);
1121
1122 nzs_mntil=nzs_mn+1;
1123 NNEW(jq_mntil,ITG,irow_ln+1);
1124 NNEW(au_mntil,double,nzs_mntil);
1125 NNEW(irow_mntil,ITG,nzs_mntil);
1126 add_rect(au_intmn,irow_intmn,jq_intmn,irow_lm,irow_ln,
1127 au_mn,irow_mn,jq_mn,irow_lm,irow_ln,
1128 &au_mntil,&irow_mntil,jq_mntil,&nzs_mntil);
1129
1130 nzs_intmm=nzs_mmf;
1131 NNEW(jq_intmm,ITG,irow_lm+1);
1132 NNEW(au_intmm,double,nzs_intmm);
1133 NNEW(irow_intmm,ITG,nzs_intmm);
1134 multi_rect(au_bd1,irow_bd1,jq_bd1,irow_la,irow_lm,
1135 au_am,irow_am,jq_am,irow_la,irow_lm,
1136 &au_intmm,&irow_intmm,jq_intmm,&nzs_intmm);
1137
1138 nzs_mmtil=nzs_mmf;
1139 NNEW(jq_mmtil,ITG,irow_lm+1);
1140 NNEW(au_mmtil,double,nzs_mmtil);
1141 NNEW(irow_mmtil,ITG,nzs_mmtil);
1142 add_rect(au_intmm,irow_intmm,jq_intmm,irow_lm,irow_lm,
1143 au_mmf,irow_mmf,jq_mmf,irow_lm,irow_lm,
1144 &au_mmtil,&irow_mmtil,jq_mmtil,&nzs_mmtil);
1145
1146 /* K_MI^til **/
1147
1148 nzs_intmi=max(nzs_mi,nzs_bdtil2);
1149 NNEW(jq_intmi,ITG,irow_li+1);
1150 NNEW(au_intmi,double,nzs_intmi);
1151 NNEW(irow_intmi,ITG,nzs_intmi);
1152 multi_rect(au_bd1,irow_bd1,jq_bd1,irow_la,irow_lm,
1153 au_ai,irow_ai,jq_ai,irow_la,irow_li,
1154 &au_intmi,&irow_intmi,jq_intmi,&nzs_intmi);
1155
1156 nzs_mitil=max(nzs_mi,nzs_intmi);
1157 NNEW(jq_mitil,ITG,irow_li+1);
1158 NNEW(au_mitil,double,nzs_mitil);
1159 NNEW(irow_mitil,ITG,nzs_mitil);
1160 add_rect(au_intmi,irow_intmi,jq_intmi,irow_lm,irow_li,
1161 au_mi,irow_mi,jq_mi,irow_lm,irow_li,
1162 &au_mitil,&irow_mitil,jq_mitil,&nzs_mitil);
1163
1164 /* K_MA^til **/
1165
1166 nzs_intma=max(nzs_ma,nzs_bdtil2);
1167 NNEW(jq_intma,ITG,irow_la+1);
1168 NNEW(au_intma,double,nzs_intma);
1169 NNEW(irow_intma,ITG,nzs_intma);
1170 multi_rect(au_bd1,irow_bd1,jq_bd1,irow_la,irow_lm,
1171 au_aaf,irow_aaf,jq_aaf,irow_la,irow_la,
1172 &au_intma,&irow_intma,jq_intma,&nzs_intma);
1173
1174 nzs_matil=max(nzs_ma,nzs_intma);
1175 NNEW(jq_matil,ITG,irow_la+1);
1176 NNEW(au_matil,double,nzs_matil);
1177 NNEW(irow_matil,ITG,nzs_matil);
1178 add_rect(au_intma,irow_intma,jq_intma,irow_lm,irow_la,
1179 au_ma,irow_ma,jq_ma,irow_lm,irow_la,
1180 &au_matil,&irow_matil,jq_matil,&nzs_matil);
1181
1182 if(debug==1)printf("\tmm2: MNtil %" ITGFORMAT " %" ITGFORMAT " MMtil %" ITGFORMAT " %"
1183 ITGFORMAT " MItil %" ITGFORMAT " %" ITGFORMAT " MAtil %" ITGFORMAT
1184 " %" ITGFORMAT " \n",nzs_intmn,nzs_mntil,nzs_intmm,nzs_mmtil,
1185 nzs_intmi,nzs_mitil,nzs_intma,nzs_matil);
1186
1187 /* deleting the original 2nd row in the matrix (master row entries) */
1188
1189 SFREE(au_mmf);SFREE(irow_mmf);SFREE(jq_mmf);
1190 SFREE(au_mn);SFREE(irow_mn);SFREE(jq_mn);
1191 SFREE(au_mi);SFREE(irow_mi);SFREE(jq_mi);
1192 SFREE(au_ma);SFREE(irow_ma);SFREE(jq_ma);
1193
1194 /* call trafoNT: transformation of the last two rows of the matrix with the normal and
1195 tangential vectors at the contact faces **/
1196
1197 NNEW(jq_antil,ITG,irow_ln+1);
1198 NNEW(au_antil,double,jq_an[irow_ln]-1);
1199 NNEW(irow_antil,ITG,jq_an[irow_ln]-1);
1200 NNEW(jq_amtil,ITG,irow_lm+1);
1201 NNEW(au_amtil,double,jq_am[irow_lm]-1);
1202 NNEW(irow_amtil,ITG,jq_am[irow_lm]-1);
1203 NNEW(jq_aitil,ITG,irow_li+1);
1204 NNEW(au_aitil,double,jq_ai[irow_li]-1);
1205 NNEW(irow_aitil,ITG,jq_ai[irow_li]-1);
1206 NNEW(jq_aatil,ITG,irow_la+1);
1207 NNEW(au_aatil,double,jq_aaf[irow_la]-1);
1208 NNEW(irow_aatil,ITG,jq_aaf[irow_la]-1);
1209 NNEW(f_atil,double,irow_la);
1210 trafontmortar2(neq,nzs,islavactdof,islavact,nslavnode,nmastnode,f_a,f_atil,
1211 au_an,irow_an,jq_an,au_am,irow_am,jq_am,au_ai,irow_ai,jq_ai,
1212 au_aaf,irow_aaf,jq_aaf,&au_antil,&irow_antil,jq_antil,
1213 &au_amtil,&irow_amtil,jq_amtil,&au_aitil,&irow_aitil,
1214 jq_aitil,&au_aatil,&irow_aatil,jq_aatil,gap,Bd,irowb,jqb,Dd,
1215 irowd,jqd,Ddtil,irowdtil,jqdtil,au_bdtil2,irow_bdtil2,
1216 jq_bdtil2,au_ddtil2i,irow_ddtil2i,jq_ddtil2i,au_ddtil2a,
1217 irow_ddtil2a,jq_ddtil2a,m_flagr,i_flagr,a_flagr,a_flag,i_flag,
1218 m_flag,&irow_ln,&irow_lm,&irow_li,&irow_la,slavnor,slavtan,vold,
1219 vini,cstress,cstressini,bp_old,nactdof,islavnode,imastnode,
1220 ntie,mi,nk,nboun,ndirboun,nodeboun,xboun,nmpc,ipompc,nodempc,
1221 coefmpc,ikboun,ilboun,ikmpc,ilmpc,nslavspc,islavspc,nsspc,
1222 nslavmpc,islavmpc,nsmpc,nmastspc,imastspc,nmspc,nmastmpc,
1223 imastmpc,nmmpc,tieset,islavactdoftie,nelcon,elcon,tietol,
1224 ncmat_,ntmat_,plicon,nplicon,npmat_,dtime,irowtloc,jqtloc,
1225 autloc,irowtlocinv,jqtlocinv,autlocinv,islavnodeinv,
1226 lambdaiwan,lambdaiwanini,iit,nmethod,bet,ithermal,plkcon,
1227 nplkcon);
1228
1229 SFREE(au_an);SFREE(irow_an);SFREE(jq_an);
1230 SFREE(au_am);SFREE(irow_am);SFREE(jq_am);
1231 SFREE(au_ai);SFREE(irow_ai);SFREE(jq_ai);
1232 SFREE(au_aaf);SFREE(irow_aaf);SFREE(jq_aaf);
1233
1234 /* Calculation of f_mhat=f_m + Bd_T*f_a **/
1235
1236 multi_rectv(au_bd1,irow_bd1,jq_bd1,irow_la,irow_lm,f_a,&v_r);
1237 //Bd_T*f_a
1238
1239 // Just Master part is concerned
1240 for(i=0;i<irow_lm;i++) {
1241 b[m_flagr[i]-1]=f_m[i]+v_r[i];}
1242
1243 SFREE(v_r);
1244 for(i=0;i<irow_li;i++) {
1245 b[i_flagr[i]-1]=f_i[i];}
1246 for(i=0;i<irow_la;i++) {
1247 b[a_flagr[i]-1]=f_atil[i];}
1248
1249 /* up to this point nzs is still the number of nonzero off-diagonals in the
1250 lower triangle of the symmetric matrix (= upper triangle transpose) */
1251
1252 nzs_sym=*nzs;
1253
1254 /* ASSEMBLE **/
1255
1256 *nzs=*nzsc;
1257 ifree=1;
1258
1259 RENEW(au,double,*nzs);
1260
1261 /* CHANGE on 20190329: */
1262
1263 for(j=nzs_sym;j<*nzs;j++){
1264 au[j]=0.;}
1265
1266 /* END CHANGE */
1267
1268 RENEW(irow,ITG,*nzs+neq[1]);
1269
1270 /* CHANGE on 20190329: */
1271
1272 for(j=nzs_sym;j<*nzs;j++){
1273 irow[j]=0;}
1274
1275 /* END CHANGE */
1276
1277 NNEW(mast1,ITG,*nzs);
1278 ITG loc_n,loc_i,loc_a,loc_m;
1279 for(j=0;j<neq[1];j++){
1280 ad[j]=0.0;}
1281 for(j=0;j<neq[1];j++){
1282 m=j+1;
1283 switch(l_flag[j]){
1284 case 0 : {
1285 // insert Matrices NN,MN,IN,AN
1286 loc_n=n_flag[j];
1287 for(l=jq_nn[loc_n-1]-1;l<jq_nn[loc_n]-1;l++){
1288 //NN
1289 insertas(&irow,&mast1,&n_flagr[irow_nn[l]-1],&m,&ifree,nzs,&au_nn[l],
1290 &au);
1291 }
1292 for(l=jq_mntil[loc_n-1]-1;l<jq_mntil[loc_n]-1;l++){
1293 //MN
1294 insertas(&irow,&mast1,&m_flagr[irow_mntil[l]-1],&m,&ifree,nzs,
1295 &au_mntil[l],&au);
1296 }
1297 for(l=jq_in[loc_n-1]-1;l<jq_in[loc_n]-1;l++){
1298 //IN
1299 insertas(&irow,&mast1,&i_flagr[irow_in[l]-1],&m,&ifree,nzs,&au_in[l],
1300 &au);
1301 }
1302 for(l=jq_antil[loc_n-1]-1;l<jq_antil[loc_n]-1;l++){
1303 //AN
1304 insertas(&irow,&mast1,&a_flagr[irow_antil[l]-1],&m,&ifree,nzs,
1305 &au_antil[l],&au);
1306 }
1307 ad[j]=ad_nn[loc_n-1];
1308 break;}
1309 case 1 :{
1310 // insert Matrices NM,MM,IM,AM
1311 loc_m=m_flag[j];
1312 for(l=jq_nm[loc_m-1]-1;l<jq_nm[loc_m]-1;l++){
1313 //NM
1314 insertas(&irow,&mast1,&n_flagr[irow_nm[l]-1],&m,&ifree,nzs,&au_nm[l],
1315 &au);
1316 }
1317 for(l=jq_mmtil[loc_m-1]-1;l<jq_mmtil[loc_m]-1;l++){
1318 //MM
1319 if ((m==m_flagr[irow_mmtil[l]-1])){
1320 //diagonal of mm
1321 ad[j]=au_mmtil[l];
1322 }else{
1323 insertas(&irow,&mast1,&m_flagr[irow_mmtil[l]-1],&m,&ifree,nzs,
1324 &au_mmtil[l],&au);
1325 }
1326 }
1327 for(l=jq_im[loc_m-1]-1;l<jq_im[loc_m]-1;l++){
1328 //IM
1329 insertas(&irow,&mast1,&i_flagr[irow_im[l]-1],&m,&ifree,nzs,&au_im[l],
1330 &au);
1331 }
1332 for(l=jq_amtil[loc_m-1]-1;l<jq_amtil[loc_m]-1;l++){
1333 //AM
1334 insertas(&irow,&mast1,&a_flagr[irow_amtil[l]-1],&m,&ifree,nzs,
1335 &au_amtil[l],&au);
1336 }
1337 break;}
1338 case 2 :{
1339 // insert Matrix NI,MI,II,AI
1340 loc_i=i_flag[j];
1341 for(l=jq_ni[loc_i-1]-1;l<jq_ni[loc_i]-1;l++){
1342 //NI
1343 insertas(&irow,&mast1,&n_flagr[irow_ni[l]-1],&m,&ifree,nzs,&au_ni[l],
1344 &au);
1345 }
1346 for(l=jq_mitil[loc_i-1]-1;l<jq_mitil[loc_i]-1;l++){
1347 //MI
1348 insertas(&irow,&mast1,&m_flagr[irow_mitil[l]-1],&m,&ifree,nzs,
1349 &au_mitil[l],&au);
1350 }
1351 for(l=jq_iif[loc_i-1]-1;l<jq_iif[loc_i]-1;l++){
1352 //II
1353 if (m==i_flagr[irow_iif[l]-1]){
1354 ad[j]=au_iif[l];
1355 }else{
1356 insertas(&irow,&mast1,&i_flagr[irow_iif[l]-1],&m,&ifree,nzs,
1357 &au_iif[l],&au);
1358 }
1359 }
1360 for(l=jq_aitil[loc_i-1]-1;l<jq_aitil[loc_i]-1;l++){
1361 //AI
1362 insertas(&irow,&mast1,&a_flagr[irow_aitil[l]-1],&m,&ifree,nzs,
1363 &au_aitil[l],&au);
1364 }
1365 break;}
1366 case 3 :{
1367 // insert Matrix NA,MA,IA,AA
1368 loc_a=a_flag[j];
1369 for(l=jq_na[loc_a-1]-1;l<jq_na[loc_a]-1;l++){
1370 //NA
1371 insertas(&irow,&mast1,&n_flagr[irow_na[l]-1],&m,&ifree,nzs,&au_na[l],
1372 &au);
1373 }
1374 for(l=jq_matil[loc_a-1]-1;l<jq_matil[loc_a]-1;l++){
1375 //MA
1376 insertas(&irow,&mast1,&m_flagr[irow_matil[l]-1],&m,&ifree,nzs,
1377 &au_matil[l],&au);
1378 }
1379 for(l=jq_ia[loc_a-1]-1;l<jq_ia[loc_a]-1;l++){
1380 //IA
1381 insertas(&irow,&mast1,&i_flagr[irow_ia[l]-1],&m,&ifree,nzs,&au_ia[l],
1382 &au);
1383 }
1384 for(l=jq_aatil[loc_a-1]-1;l<jq_aatil[loc_a]-1;l++){
1385 //AA
1386 if (m==a_flagr[irow_aatil[l]-1]){
1387 ad[j]=au_aatil[l];
1388 }else{
1389 insertas(&irow,&mast1,&a_flagr[irow_aatil[l]-1],&m,&ifree,nzs,
1390 &au_aatil[l],&au);
1391 }
1392 }
1393 break; }
1394 default : {
1395 break;}
1396 }
1397 }
1398
1399 // sort pro column
1400
1401 *nzs=ifree-1;
1402 RENEW(au,double,*nzs);
1403 RENEW(irow,ITG,*nzs);
1404 dim=neq[1];
1405 matrixsort(au,mast1,irow,jq,nzs,&dim);
1406 SFREE(mast1);
1407
1408 /*********** Free the intermediate matrices ********/
1409
1410 SFREE(au_nn);SFREE(irow_nn);SFREE(jq_nn);
1411 SFREE(au_antil);SFREE(irow_antil);SFREE(jq_antil);
1412 SFREE(au_nm);SFREE(irow_nm);SFREE(jq_nm);
1413 SFREE(au_amtil);SFREE(irow_amtil);SFREE(jq_amtil);
1414 SFREE(au_ni);SFREE(irow_ni);SFREE(jq_ni);
1415 SFREE(au_aitil);SFREE(irow_aitil);SFREE(jq_aitil);
1416 SFREE(au_na);SFREE(irow_na);SFREE(jq_na);
1417 SFREE(au_aatil);SFREE(irow_aatil);SFREE(jq_aatil);
1418 SFREE(au_bd1);SFREE(irow_bd1);SFREE(jq_bd1);
1419 SFREE(au_bdtil2);SFREE(irow_bdtil2);SFREE(jq_bdtil2);
1420 SFREE(au_ddtil2i);SFREE(irow_ddtil2i);SFREE(jq_ddtil2i);
1421 SFREE(au_ddtil2a);SFREE(irow_ddtil2a);SFREE(jq_ddtil2a);
1422 SFREE(ad_ii);SFREE(irow_iid);SFREE(jq_iid);
1423 SFREE(ad_aa);SFREE(irow_aad);SFREE(jq_aad);
1424 SFREE(ad_mm);SFREE(irow_mmd);SFREE(jq_mmd);
1425 SFREE(au_intmn);SFREE(irow_intmn);SFREE(jq_intmn);
1426 SFREE(au_intmm);SFREE(irow_intmm);SFREE(jq_intmm);
1427 SFREE(au_intmi);SFREE(irow_intmi);SFREE(jq_intmi);
1428 SFREE(au_intma);SFREE(irow_intma);SFREE(jq_intma);
1429 SFREE(au_mntil);SFREE(irow_mntil);SFREE(jq_mntil);
1430 SFREE(au_mmtil);SFREE(irow_mmtil);SFREE(jq_mmtil);
1431 SFREE(au_mitil);SFREE(irow_mitil);SFREE(jq_mitil);
1432 SFREE(au_matil);SFREE(irow_matil);SFREE(jq_matil);
1433 SFREE(f_a); SFREE(f_i);
1434 SFREE(f_m); SFREE(f_atil);
1435 SFREE(ad_nn);
1436 SFREE(au_in);SFREE(irow_in);SFREE(jq_in);
1437 SFREE(au_im);SFREE(irow_im);SFREE(jq_im);
1438 SFREE(au_ia);SFREE(irow_ia);SFREE(jq_ia);
1439 SFREE(au_iif);SFREE(irow_iif);SFREE(jq_iif);
1440 SFREE(au_t);SFREE(jq_t);SFREE(irow_t);
1441 SFREE(l_flag);SFREE(n_flag);SFREE(i_flag);SFREE(a_flag);SFREE(m_flag);
1442 SFREE(n_flagr);SFREE(i_flagr);SFREE(a_flagr);SFREE(m_flagr);
1443
1444 /*************************/
1445 /*END transmit the new stiffness matrix*/
1446
1447 if(debug==1)printf("\tnzsc %" ITGFORMAT " nzs %" ITGFORMAT " \n",*nzsc,*nzs);
1448 if(debug==1){
1449 number=3;
1450 FORTRAN(writematrix,(auc,adc,irowc,jqc,&neq[1],&number));
1451
1452 number=4;
1453 FORTRAN(writematrix,(au,ad,irow,jq,&neq[1],&number));
1454 printf("\n");
1455 number=5;
1456 FORTRAN(writevector,(b,&neq[1],&number));
1457 number=6;
1458 FORTRAN(writevector,(bhat,&neq[1],&number));
1459 }
1460
1461 *irowcp=irowc; *aucp=auc;
1462 *irowp=irow; *aup=au;
1463 return;
1464 }
1465