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 <string.h>
22 #include "CalculiX.h"
23
remastruct(ITG * ipompc,double ** coefmpcp,ITG ** nodempcp,ITG * nmpc,ITG * mpcfree,ITG * nodeboun,ITG * ndirboun,ITG * nboun,ITG * ikmpc,ITG * ilmpc,ITG * ikboun,ITG * ilboun,char * labmpc,ITG * nk,ITG * memmpc_,ITG * icascade,ITG * maxlenmpc,ITG * kon,ITG * ipkon,char * lakon,ITG * ne,ITG * nactdof,ITG * icol,ITG * jq,ITG ** irowp,ITG * isolver,ITG * neq,ITG * nzs,ITG * nmethod,double ** fp,double ** fextp,double ** bp,double ** aux2p,double ** finip,double ** fextinip,double ** adbp,double ** aubp,ITG * ithermal,ITG * iperturb,ITG * mass,ITG * mi,ITG * iexpl,ITG * mortar,char * typeboun,double ** cvp,double ** cvinip,ITG * iit,ITG * network,ITG * itiefac,ITG * ne0,ITG * nkon0,ITG * nintpoint,ITG * islavsurf,double * pmastsurf,char * tieset,ITG * ntie,ITG * num_cpus)24 void remastruct(ITG *ipompc, double **coefmpcp, ITG **nodempcp, ITG *nmpc,
25 ITG *mpcfree, ITG *nodeboun, ITG *ndirboun, ITG *nboun,
26 ITG *ikmpc, ITG *ilmpc, ITG *ikboun, ITG *ilboun,
27 char *labmpc, ITG *nk,
28 ITG *memmpc_, ITG *icascade, ITG *maxlenmpc,
29 ITG *kon, ITG *ipkon, char *lakon, ITG *ne,
30 ITG *nactdof, ITG *icol, ITG *jq, ITG **irowp, ITG *isolver,
31 ITG *neq, ITG *nzs,ITG *nmethod, double **fp,
32 double **fextp, double **bp, double **aux2p, double **finip,
33 double **fextinip,double **adbp, double **aubp, ITG *ithermal,
34 ITG *iperturb, ITG *mass, ITG *mi,ITG *iexpl,ITG *mortar,
35 char *typeboun,double **cvp,double **cvinip,ITG *iit,
36 ITG *network,ITG *itiefac,ITG *ne0,ITG *nkon0,ITG *nintpoint,
37 ITG *islavsurf,double *pmastsurf,char*tieset,ITG *ntie,
38 ITG *num_cpus){
39
40 /* reconstructs the nonzero locations in the stiffness and mass
41 matrix after a change in MPC's */
42
43 char *lakontot=NULL;
44
45 ITG *nodempc=NULL,*mast1=NULL,*ipointer=NULL,mpcend,isiz,
46 callfrommain,i,*irow=NULL,mt,im,*ipkontot=NULL,*kontot=NULL,
47 netot;
48
49 double *coefmpc=NULL,*f=NULL,*fext=NULL,*b=NULL,*aux2=NULL,
50 *fini=NULL,*fextini=NULL,*adb=NULL,*aub=NULL,*cv=NULL,*cvini=NULL;
51
52 nodempc=*nodempcp;coefmpc=*coefmpcp;irow=*irowp;
53 f=*fp;fext=*fextp;b=*bp;aux2=*aux2p;fini=*finip;
54 fextini=*fextinip;adb=*adbp;aub=*aubp;cv=*cvp;cvini=*cvinip;
55
56 mt=mi[1]+1;
57
58 /* decascading the MPC's */
59
60 if(*icascade>0){
61 if(*iexpl<=1) printf(" Decascading the MPC's\n\n");
62
63 callfrommain=0;
64 cascade(ipompc,&coefmpc,&nodempc,nmpc,
65 mpcfree,nodeboun,ndirboun,nboun,ikmpc,
66 ilmpc,ikboun,ilboun,&mpcend,
67 labmpc,nk,memmpc_,icascade,maxlenmpc,
68 &callfrommain,iperturb,ithermal);
69 }
70
71 /* determining the matrix structure */
72
73 if(*iexpl<=1) printf(" Determining the structure of the matrix:\n");
74
75 if(nzs[1]<10) nzs[1]=10;
76 NNEW(mast1,ITG,nzs[1]);
77 NNEW(ipointer,ITG,mt**nk);
78
79 if(*mortar==1){
80
81 /* for face-to-face penalty: take all possible contact
82 elements into account: this number is constant during an
83 increment. Therefore, the location of nonzero's does not change
84 and the renumbering during the equation solution (PaStiX) has only
85 to be done once*/
86
87 NNEW(ipkontot,ITG,*ne0+*nintpoint);
88 NNEW(kontot,ITG,*nkon0+22**nintpoint);
89 NNEW(lakontot,char,8*(*ne0+*nintpoint));
90
91 isiz=*ne0;cpyparitg(ipkontot,ipkon,&isiz,num_cpus);
92 isiz=*nkon0;cpyparitg(kontot,kon,&isiz,num_cpus);
93 memcpy(&lakontot[0],&lakon[0],sizeof(char)*8**ne0);
94
95 FORTRAN(totalcontact,(tieset,ntie,&netot,ipkontot,kontot,lakontot,
96 islavsurf,itiefac,pmastsurf,ne0,nkon0));
97
98 printf("maximal possible contact elements = \n%d\n\n",netot-*ne0);
99
100 mastruct(nk,kontot,ipkontot,lakontot,&netot,nodeboun,ndirboun,nboun,
101 ipompc,nodempc,nmpc,nactdof,icol,jq,&mast1,&irow,isolver,neq,
102 ikmpc,ilmpc,ipointer,nzs,nmethod,ithermal,
103 ikboun,ilboun,iperturb,mi,mortar,typeboun,labmpc,
104 iit,icascade,network,iexpl);
105
106 SFREE(ipkontot);SFREE(kontot);SFREE(lakontot);
107
108 }else{
109
110 /* take actual contact elements into account */
111
112 mastruct(nk,kon,ipkon,lakon,ne,nodeboun,ndirboun,nboun,ipompc,
113 nodempc,nmpc,nactdof,icol,jq,&mast1,&irow,isolver,neq,
114 ikmpc,ilmpc,ipointer,nzs,nmethod,ithermal,
115 ikboun,ilboun,iperturb,mi,mortar,typeboun,labmpc,
116 iit,icascade,network,iexpl);
117 }
118
119 SFREE(ipointer);SFREE(mast1);
120 RENEW(irow,ITG,nzs[2]);
121
122 *nodempcp=nodempc;*coefmpcp=coefmpc;*irowp=irow;
123
124 /* reallocating fields the size of which depends on neq[1] or *nzs */
125
126 RENEW(f,double,neq[1]);DMEMSET(f,0,neq[1],0.);
127 RENEW(fext,double,neq[1]);DMEMSET(fext,0,neq[1],0.);
128 RENEW(b,double,neq[1]);DMEMSET(b,0,neq[1],0.);
129 RENEW(fini,double,neq[1]);
130
131 /* for static calculations fini has to be set to f at the
132 start of the calculation; in dynamic calculations this is
133 not needed, since the initial accelerations has already
134 been calculated */
135
136 if((*nmethod!=4)&&(*iit==-1)) DMEMSET(fini,0,neq[1],0.);
137
138 if(*nmethod==4){
139 RENEW(aux2,double,neq[1]);DMEMSET(aux2,0,neq[1],0.);
140 RENEW(cv,double,neq[1]);
141 RENEW(cvini,double,neq[1]);
142 RENEW(fextini,double,neq[1]);
143
144 /* the mass matrix is diagonal in an explicit dynamic
145 calculation and is not changed by contact; this
146 assumes that the number of degrees of freedom does
147 not change */
148
149 if(*iexpl<=1){
150 RENEW(adb,double,neq[1]);for(i=0;i<neq[1];i++) adb[i]=0.;
151 RENEW(aub,double,nzs[1]);for(i=0;i<nzs[1];i++) aub[i]=0.;
152 mass[0]=1;
153 }
154 }
155
156 *fp=f;*fextp=fext;*bp=b;*aux2p=aux2;*finip=fini;
157 *fextinip=fextini;*adbp=adb;*aubp=aub;*cvp=cv;*cvinip=cvini;
158
159 return;
160 }
161
162
163