1 /*
2 CalculiX - A 3-dimensional finite element program
3 Copyright (C) 1998-2021 Guido Dhondt
4
5 This program is free software; you can redistribute it and/or
6 modify it under the terms of the GNU General Public License as
7 published by the Free Software Foundation(version 2);
8
9
10 This program is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 GNU General Public License for more details.
14
15 You should have received a copy of the GNU General Public License
16 along with this program; if not, write to the Free Software
17 Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
18 */
19 /* */
20 /* author: Reinhold Fischer */
21
22 #include <stdio.h>
23 #include <math.h>
24 #include <stdlib.h>
25 #include <string.h>
26 #include "CalculiX.h"
27
dfdbj(double * bcont,double ** dbcontp,ITG * neq,ITG * nope,ITG * konl,ITG * nactdof,double * s,double * z,ITG * ikmpc,ITG * ilmpc,ITG * ipompc,ITG * nodempc,ITG * nmpc,double * coefmpc,double * fnl,ITG * nev,ITG ** ikactcontp,ITG ** ilactcontp,ITG * nactcont,ITG * nactcont_,ITG * mi,ITG * cyclicsymmetry,ITG * izdof,ITG * nzdof)28 void dfdbj(double *bcont,double **dbcontp,ITG *neq,ITG *nope,ITG *konl,
29 ITG* nactdof,double *s,double *z,ITG *ikmpc,ITG *ilmpc,
30 ITG *ipompc,ITG *nodempc,ITG *nmpc,double *coefmpc,
31 double *fnl,ITG *nev,ITG **ikactcontp,ITG **ilactcontp,
32 ITG *nactcont,ITG *nactcont_,ITG *mi, ITG *cyclicsymmetry,
33 ITG *izdof, ITG *nzdof){
34
35 ITG j,j1,jdof,kdof,k,k1,l,id,index,ist,id1,ist1,index1,id2,ist2,index2,
36 jdbcontcol,i1,i3,i4,mt=mi[1]+1,im,*ikactcont=*ikactcontp,
37 *ilactcont=*ilactcontp,kdofm1;
38
39 double d1,sl,*dbcont=*dbcontp;
40
41 for(j=0; j<*nope; j++){
42 i1=mt*(konl[j]-1)+1;
43 for(j1=0; j1<3; j1++){
44 jdof=nactdof[i1+j1];
45 if(jdof>0){
46 jdof--;
47 FORTRAN(nident,(ikactcont,&jdof,nactcont,&id));
48 do{
49 if(id>0){
50 if(ikactcont[id-1]==jdof){
51 jdbcontcol=ilactcont[id-1];
52 break;
53 }
54 }
55 (*nactcont)++;
56 if(*nactcont>*nactcont_){
57 *nactcont_=(ITG)(1.1**nactcont_);
58 RENEW(ikactcont,ITG,*nactcont_);
59 RENEW(ilactcont,ITG,*nactcont_);
60 RENEW(dbcont,double,*nev**nactcont_);
61 }
62 k=*nactcont-1;
63 l=k-1;
64 while(k>id){
65 ikactcont[k]=ikactcont[l];
66 ilactcont[k--]=ilactcont[l--];
67 }
68 jdbcontcol=*nactcont;
69 ikactcont[id]=jdof;
70 ilactcont[id]=*nactcont;
71 // memset(&dbcont[(*nactcont-1)**nev],0,sizeof(double)**nev);
72 DMEMSET(dbcont,(*nactcont-1)**nev,*nactcont**nev,0.);
73 break;
74 }while(1);
75 bcont[jdof]-=fnl[j*3+j1];
76 i4=(jdbcontcol-1)**nev;
77 i3=(3*j+j1);
78 for(k=0; k<*nope; k++){
79 for(k1=0; k1<3; k1++){
80 sl=s[(3*k+k1)*60+i3];
81 kdof=nactdof[mt*(konl[k]-1)+k1+1];
82 if(kdof>0){
83 if(!(*cyclicsymmetry)){
84 for(l=0; l<*nev; l++){
85 dbcont[i4+l]-=sl*z[(long long)l**neq+kdof-1];
86 }
87 }else{
88 kdofm1=kdof-1;
89 FORTRAN(nident,(izdof,&kdofm1,nzdof,&id));
90 if(id!=0){
91 if(izdof[id-1]==kdofm1){
92 for(l=0; l<*nev; l++){
93 dbcont[i4+l]-=sl*z[l**nzdof+id-1];
94 }
95 }else{printf("*ERROR in dfdbj\n");FORTRAN(stop,());}
96 }else{printf("*ERROR in dfdbj\n");FORTRAN(stop,());}
97 }
98 }
99 else{
100 kdof=8*(konl[k]-1)+k1+1;
101 FORTRAN(nident,(ikmpc,&kdof,nmpc,&id));
102 if(id>0){
103 id--;
104 if(ikmpc[id]==kdof){
105 id=ilmpc[id];
106 ist=ipompc[id-1];
107 ist--;
108 index=nodempc[ist*3+2];
109 if(index==0) continue;
110 index--;
111 do{
112 kdof=nactdof[mt*(nodempc[index*3]-1)+nodempc[index*3+1]];
113 d1=sl*coefmpc[index]/coefmpc[ist];
114 if(kdof>0){
115 if(!(*cyclicsymmetry)){
116 for(l=0; l<*nev; l++){
117 dbcont[i4+l]+=d1*z[(long long)l**neq+kdof-1];
118 }
119 }
120 }else{
121 kdofm1=kdof-1;
122 FORTRAN(nident,(izdof,&kdofm1,nzdof,&id));
123 if(id!=0){
124 if(izdof[id-1]==kdofm1){
125 for(l=0; l<*nev; l++){
126 dbcont[i4+l]+=d1*z[l**nzdof+id-1];
127 }
128 }else{printf("*ERROR in dfdbj\n");FORTRAN(stop,());}
129 }else{printf("*ERROR in dfdbj\n");FORTRAN(stop,());}
130 }
131 index=nodempc[index*3+2];
132 if(index==0) break;
133 index--;
134 }while(1);
135 }
136 }
137 }
138 }
139 }
140 }
141 else{
142 jdof=8*(konl[j]-1)+j1+1;
143 FORTRAN(nident,(ikmpc,&jdof,nmpc,&id1));
144 if(id1>0){
145 id1--;
146 if(ikmpc[id1]==jdof){
147 id1=ilmpc[id1];
148 ist1=ipompc[id1-1];
149 ist1--;
150 index1=nodempc[ist1*3+2];
151 if(index1==0) continue;
152 index1--;
153 do{
154 jdof=nactdof[mt*(nodempc[index1*3]-1)+nodempc[index1*3+1]];
155 if(jdof>0){
156 jdof--;
157 FORTRAN(nident,(ikactcont,&jdof,nactcont,&id));
158 do{
159 if(id>0){
160 if(ikactcont[id-1]==jdof){
161 jdbcontcol=ilactcont[id-1];
162 }
163 }
164 (*nactcont)++;
165 if(*nactcont>*nactcont_){
166 *nactcont_=(ITG)(1.1**nactcont_);
167 RENEW(ikactcont,ITG,*nactcont_);
168 RENEW(ilactcont,ITG,*nactcont_);
169 RENEW(dbcont,double,*nev**nactcont_);
170 }
171 k=*nactcont-1;
172 l=k-1;
173 do{
174 ikactcont[k]=ikactcont[l];
175 ilactcont[k--]=ilactcont[l--];
176 }while(k>id);
177 jdbcontcol=*nactcont;
178 ikactcont[id]=jdof;
179 ilactcont[id]=*nactcont;
180 // memset(&dbcont[(*nactcont-1)**nev],0,sizeof(double)**nev);
181 DMEMSET(dbcont,(*nactcont-1)**nev,*nactcont**nev,0.);
182 break;
183 }while(1);
184 bcont[jdof]+=coefmpc[index1]*fnl[j*3+j1]/coefmpc[ist1];
185 i4=(jdbcontcol-1)**nev;
186 i3=(3*j+j1);
187 for(k=0; k<*nope; k++){
188 for(k1=0; k1<3; k1++){
189 sl=s[(3*k+k1)*60+i3];
190 kdof=nactdof[mt*(konl[k]-1)+k1+1];
191 if(kdof>0){
192 d1=sl*coefmpc[index1]/coefmpc[ist1];
193 if(!(*cyclicsymmetry)){
194 for(l=0; l<*nev; l++){
195 dbcont[i4+l]+=d1*z[(long long)l**neq+kdof-1];
196 }
197 }else{
198 kdofm1=kdof-1;
199 FORTRAN(nident,(izdof,&kdofm1,nzdof,&id));
200 if(id!=0){
201 if(izdof[id-1]==kdofm1){
202 for(l=0; l<*nev; l++){
203 dbcont[i4+l]+=d1*z[l**nzdof+id-1];
204 }
205 }else{printf("*ERROR in dfdbj\n");FORTRAN(stop,());}
206 }else{printf("*ERROR in dfdbj\n");FORTRAN(stop,());}
207 }
208 }
209 else{
210 kdof=8*(konl[k]-1)+k1+1;
211 FORTRAN(nident,(ikmpc,&kdof,nmpc,&id2));
212 if(id2>0){
213 id2--;
214 if(ikmpc[id2]==kdof){
215 id2=ilmpc[id2];
216 ist2=ipompc[id2-1];
217 ist2--;
218 index2=nodempc[ist2*3+2];
219 if(index2==0) continue;
220 index2--;
221 do{
222 kdof=nactdof[mt*(nodempc[index2*3]-1)+nodempc[index2*3+1]];
223 if(kdof>0){
224 d1=sl*coefmpc[index1]*coefmpc[index2]/(coefmpc[ist1]*coefmpc[ist2]);
225 if(!(*cyclicsymmetry)){
226 for(l=0; l<*nev; l++){
227 dbcont[i4+l]-=d1*z[(long long)l**neq+kdof-1];
228 }
229 }else{
230 kdofm1=kdof-1;
231 FORTRAN(nident,(izdof,&kdofm1,nzdof,&id));
232 if(id!=0){
233 if(izdof[id-1]==kdofm1){
234 for(l=0; l<*nev; l++){
235 dbcont[i4+l]-=d1*z[l**nzdof+id-1];
236 }
237 }else{printf("*ERROR in dfdbj\n");FORTRAN(stop,());}
238 }else{printf("*ERROR in dfdbj\n");FORTRAN(stop,());}
239 }
240 }
241 index2=nodempc[index2*3+2];
242 if(index2==0) break;
243 index2--;
244 }while(1);
245 }
246 }
247 }
248 }
249 }
250 }
251 index1=nodempc[index1*3+2];
252 if(index1==0) break;
253 index1--;
254 }while(1);
255 }
256 }
257 }
258 }
259 }
260 *dbcontp=dbcont;
261 *ikactcontp=ikactcont;
262 *ilactcontp=ilactcont;
263 }
264
265 /*!
266 ! CalculiX - A 3-dimensional finite element program
267 ! Copyright (C) 1998-2021 Guido Dhondt
268 !
269 ! This program is free software; you can redistribute it and/or
270 ! modify it under the terms of the GNU General Public License as
271 ! published by the Free Software Foundation(version 2);
272 !
273 !
274 ! This program is distributed in the hope that it will be useful,
275 ! but WITHOUT ANY WARRANTY; without even the implied warranty of
276 ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
277 ! GNU General Public License for more details.
278 !
279 ! You should have received a copy of the GNU General Public License
280 ! along with this program; if not, write to the Free Software
281 ! Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
282 !
283 subroutine dfdbj(bcont,dbcont,neq,nope,konl,nactdof,s,z,
284 & ikmpc,ilmpc,ipompc,nodempc,nmpc,coefmpc,fnl,nev,iactcont,
285 & nactcont)
286 !
287 ! calculates the derivative of the contact forces with respect
288 ! to the modal variables
289 !
290 implicit none
291 !
292 integer j,j1,neq,nope,konl(*),nactdof(0:3,*),jdof,kdof,
293 & k,k1,l,id,ikmpc(*),ilmpc(*),ipompc(*),nodempc(3,*),nmpc,
294 & index,ist,id1,ist1,index1,id2,ist2,index2,nev ,iactcont(*),
295 & nactcont,jdofcont
296 !
297 real*8 bcont(*),dbcont(nev,*),s(60,60),z(neq,*),coefmpc(*),
298 & fnl(3,9)
299 !
300 do j=1,nope
301 do j1=1,3
302 jdof=nactdof(j1,konl(j))
303 if(jdof.ne.0) then
304 call nident(iactcont,jdof,nactcont,id)
305 jdofcont=0
306 if(id.gt.0)then
307 if(iactcont(id).eq.jdof) then
308 jdofcont=id
309 endif
310 endif
311 if(jdofcont.eq.0) then
312 nactcont=nactcont+1
313 do k=nactcont,id+2,-1
314 iactcont(k)=iactcont(k-1)
315 do l=1,nev
316 dbcont(l,k)=dbcont(l,k-1)
317 enddo
318 enddo
319 jdofcont=id+1
320 iactcont(jdofcont)=jdof
321 do l=1,nev
322 dbcont(l,jdofcont)=0.d0
323 enddo
324 endif
325 bcont(jdof)=bcont(jdof)-fnl(j1,j)
326 do k=1,nope
327 do k1=1,3
328 kdof=nactdof(k1,konl(k))
329 if(kdof.ne.0) then
330 do l=1,nev
331 dbcont(l,jdofcont)=dbcont(l,jdofcont)-
332 & s(3*(j-1)+j1,3*(k-1)+k1)*z(kdof,l)
333 enddo
334 else
335 kdof=8*(konl(k)-1)+k1
336 call nident(ikmpc,kdof,nmpc,id)
337 if(id.gt.0) then
338 if(ikmpc(id).eq.kdof) then
339 id=ilmpc(id)
340 ist=ipompc(id)
341 index=nodempc(3,ist)
342 if(index.eq.0) cycle
343 do
344 kdof=nactdof(nodempc(2,index),
345 & nodempc(1,index))
346 if(kdof.ne.0) then
347 do l=1,nev
348 dbcont(l,jdofcont)=
349 & dbcont(l,jdofcont)+
350 & s(3*(j-1)+j1,3*(k-1)+k1)*
351 & coefmpc(index)*z(kdof,l)/
352 & coefmpc(ist)
353 enddo
354 endif
355 index=nodempc(3,index)
356 if(index.eq.0) exit
357 enddo
358 endif
359 endif
360 endif
361 enddo
362 enddo
363 else
364 jdof=8*(konl(j)-1)+j1
365 call nident(ikmpc,jdof,nmpc,id1)
366 if(id1.gt.0) then
367 if(ikmpc(id1).eq.jdof) then
368 id1=ilmpc(id1)
369 ist1=ipompc(id1)
370 index1=nodempc(3,ist1)
371 if(index1.eq.0) cycle
372 do
373 jdof=nactdof(nodempc(2,index1),
374 & nodempc(1,index1))
375 if(jdof.ne.0) then
376 call nident(iactcont,jdof,nactcont,id)
377 jdofcont=0
378 if(id.gt.0)then
379 if(iactcont(id).eq.jdof) then
380 jdofcont=id
381 endif
382 endif
383 if(jdofcont.eq.0) then
384 nactcont=nactcont+1
385 do k=nactcont,id+2,-1
386 iactcont(k)=iactcont(k-1)
387 do l=1,nev
388 dbcont(l,k)=dbcont(l,k-1)
389 enddo
390 enddo
391 jdofcont=id+1
392 iactcont(jdofcont)=jdof
393 do l=1,nev
394 dbcont(l,jdofcont)=0.d0
395 enddo
396 endif
397 bcont(jdofcont)=bcont(jdofcont)+
398 & coefmpc(index1)*
399 & fnl(j1,j)/coefmpc(ist1)
400 do k=1,nope
401 do k1=1,3
402 kdof=nactdof(k1,konl(k))
403 if(kdof.ne.0) then
404 do l=1,nev
405 dbcont(l,jdofcont)=
406 & dbcont(l,jdofcont)
407 & +s(3*(j-1)+j1,3*(k-1)+k1)
408 & *coefmpc(index1)*z(kdof,l)/
409 & coefmpc(ist1)
410 enddo
411 else
412 kdof=8*(konl(k)-1)+k1
413 call nident(ikmpc,kdof,nmpc,id2)
414 if(id2.gt.0) then
415 if(ikmpc(id2).eq.kdof) then
416 id2=ilmpc(id2)
417 ist2=ipompc(id2)
418 index2=nodempc(3,ist2)
419 if(index2.eq.0) cycle
420 do
421 !
422 ! translated to the left to avoid exceedance
423 ! of 72 columns
424 !
425 kdof=nactdof(nodempc(2,index2),
426 & nodempc(1,index2))
427 if(kdof.ne.0) then
428 do l=1,nev
429 dbcont(l,jdofcont)=dbcont(l,jdofcont)
430 & -s(3*(j-1)+j1,3*(k-1)+k1)
431 & *coefmpc(index1)
432 & *coefmpc(index2)*z(kdof,l)/
433 & (coefmpc(ist1)*coefmpc(ist2))
434 enddo
435 endif
436 index2=nodempc(3,index2)
437 if(index2.eq.0) exit
438 !
439 ! end of translation
440 !
441 enddo
442 endif
443 endif
444 endif
445 enddo
446 enddo
447 endif
448 index1=nodempc(3,index1)
449 if(index1.eq.0) exit
450 enddo
451 endif
452 endif
453 endif
454 enddo
455 enddo
456 !
457 return
458 end
459 */
460