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
24 #define min(a,b) ((a) <= (b) ? (a) : (b))
25 #define max(a,b) ((a) >= (b) ? (a) : (b))
26
mastructse(ITG * kon,ITG * ipkon,char * lakon,ITG * ne,ITG * ipompc,ITG * nodempc,ITG * nmpc,ITG * nactdof,ITG * jqs,ITG ** mast1p,ITG ** irowsp,ITG * ipointer,ITG * nzss,ITG * mi,ITG * mortar,ITG * nodedesi,ITG * ndesi,ITG * icoordinate,ITG * ielorien,ITG * istartdesi,ITG * ialdesi)27 void mastructse(ITG *kon, ITG *ipkon, char *lakon, ITG *ne,
28 ITG *ipompc, ITG *nodempc, ITG *nmpc,
29 ITG *nactdof, ITG *jqs, ITG **mast1p, ITG **irowsp,
30 ITG *ipointer, ITG *nzss, ITG *mi, ITG *mortar,
31 ITG *nodedesi, ITG *ndesi,ITG *icoordinate,ITG *ielorien,
32 ITG *istartdesi,ITG *ialdesi){
33
34 /* determines the structure of the sensitivity matrix;
35 (i.e. the location of the nonzeros */
36
37 char lakonl[2]=" \0";
38
39 ITG i,j,k,id,index,jdof1,idof1,idof2,nmast,ifree,kdof1,
40 node,ist,kflag,indexe,nope,isize,*mast1=NULL,ii,
41 *irows=NULL,mt=mi[1]+1,*next=NULL,jstart,idesvar;
42
43 mast1=*mast1p;
44 irows=*irowsp;
45 ifree=0;
46 kflag=2;
47
48 NNEW(next,ITG,*nzss);
49
50 for(idesvar=0;idesvar<*ndesi;idesvar++){
51 idof2=idesvar+1;
52
53 for(ii=istartdesi[idesvar]-1;ii<istartdesi[idesvar+1]-1;ii++){
54 i=ialdesi[ii]-1;
55
56 if(ipkon[i]<0) continue;
57 if(strcmp1(&lakon[8*i],"F")==0)continue;
58 indexe=ipkon[i];
59 /* Bernhardi start */
60 if (strcmp1(&lakon[8*i+3],"8I")==0)nope=11;
61 else if(strcmp1(&lakon[8*i+3],"20")==0)nope=20;
62 /* Bernhardi end */
63 else if (strcmp1(&lakon[8*i+3],"8")==0)nope=8;
64 else if (strcmp1(&lakon[8*i+3],"10")==0)nope=10;
65 else if ((strcmp1(&lakon[8*i+3],"4")==0)||
66 (strcmp1(&lakon[8*i+2],"4")==0)) nope=4;
67 else if (strcmp1(&lakon[8*i+3],"15")==0)nope=15;
68 else if (strcmp1(&lakon[8*i+3],"6")==0)nope=6;
69 else if (strcmp1(&lakon[8*i],"E")==0){
70 if((strcmp1(&lakon[8*i+6],"C")==0)&&(*mortar==1)){
71
72 /* face-to-face contact (all nodes already belong
73 to other elements */
74
75 continue;
76 }else if(strcmp1(&lakon[8*i+6],"F")!=0){
77
78 /* node-to-face contact */
79
80 lakonl[0]=lakon[8*i+7];
81 nope=atoi(lakonl)+1;
82 }else{
83
84 /* advection elements */
85
86 continue;
87 }
88 }else continue;
89
90 /* displacement degrees of freedom */
91
92 for(j=0;j<nope;++j){
93 node=kon[indexe+j]-1;
94 for(k=1;k<4;k++){
95 idof1=nactdof[mt*node+k];
96 if(idof1>0){
97 insertfreq(ipointer,&mast1,&next,
98 &idof1,&idof2,&ifree,nzss);
99 }else if(*nmpc!=0){
100 if(idof1!=2*(idof1/2)){
101 id=(-idof1+1)/2-1;
102 ist=ipompc[id]-1;
103 index=nodempc[3*ist+2];
104 if(index==0) continue;
105 do{
106 jdof1=nactdof[mt*(nodempc[3*index-3]-1)
107 +nodempc[3*index-2]];
108 if(jdof1>0){
109 insertfreq(ipointer,&mast1,&next,
110 &jdof1,&idof2,&ifree,nzss);
111 }
112 index=nodempc[3*index-1];
113 if(index==0) break;
114 }while(1);
115 }
116 }
117 }
118 }
119 }
120 }
121
122 /* determine irows and jqs */
123
124 RENEW(irows,ITG,ifree);
125 nmast=0;
126 jqs[0]=1;
127 for(i=0;i<*ndesi;i++){
128 index=ipointer[i];
129 do{
130 if(index==0) break;
131 irows[nmast++]=mast1[index-1];
132 index=next[index-1];
133 }while(1);
134 jqs[i+1]=nmast+1;
135 }
136
137 /* sorting the row numbers within each column */
138
139 for(i=0;i<*ndesi;++i){
140 if(jqs[i+1]-jqs[i]>0){
141 isize=jqs[i+1]-jqs[i];
142 FORTRAN(isortii,(&irows[jqs[i]-1],&mast1[jqs[i]-1],&isize,&kflag));
143 }
144 }
145
146 /* removing duplicate entries */
147
148 nmast=0;
149 for(i=0;i<*ndesi;i++){
150 jstart=nmast+1;
151 if(jqs[i+1]-jqs[i]>0){
152 irows[nmast++]=irows[jqs[i]-1];
153 for(j=jqs[i];j<jqs[i+1]-1;j++){
154 if(irows[j]==irows[nmast-1])continue;
155 irows[nmast++]=irows[j];
156 }
157 }
158 jqs[i]=jstart;
159 }
160 jqs[*ndesi]=nmast+1;
161
162 *nzss=jqs[*ndesi]-1;
163
164 SFREE(next);
165
166 *mast1p=mast1;
167 *irowsp=irows;
168
169 return;
170
171 }
172