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 <string.h>
22 #include "CalculiX.h"
23 
radcyc(ITG * nk,ITG * kon,ITG * ipkon,char * lakon,ITG * ne,double * cs,ITG * mcs,ITG * nkon,ITG * ialset,ITG * istartset,ITG * iendset,ITG ** kontrip,ITG * ntri,double ** cop,double ** voldp,ITG * ntrit,ITG * inocs,ITG * mi)24 void radcyc(ITG *nk,ITG *kon,ITG *ipkon,char *lakon,ITG *ne,
25 	    double *cs, ITG *mcs, ITG *nkon,ITG *ialset, ITG *istartset,
26             ITG *iendset,ITG **kontrip,ITG *ntri,
27             double **cop, double **voldp,ITG *ntrit, ITG *inocs,
28             ITG *mi){
29 
30   /* duplicates triangular faces for cyclic radiation conditions */
31 
32   char *filab=NULL;
33 
34   ITG i,is,nsegments,idtie,nkt,icntrl,imag=0,*kontri=NULL,mt=mi[1]+1,
35      node,i1,i2,nope,iel,indexe,j,k,ielset,node1,node2,node3,l,jj;
36 
37   double *vt=NULL,*fnt=NULL,*stnt=NULL,*eent=NULL,*qfnt=NULL,t[3],theta,
38      pi,*v=NULL,*fn=NULL,*stn=NULL,*een=NULL,*qfn=NULL,*co=NULL,
39      *vold=NULL,*emnt=NULL,*emn=NULL;
40 
41   pi=4.*atan(1.);
42 
43   kontri=*kontrip;co=*cop;vold=*voldp;
44 
45   /* determining the maximum number of sectors */
46 
47   nsegments=1;
48   for(j=0;j<*mcs;j++){
49       if(cs[17*j]>nsegments) nsegments=(ITG)(cs[17*j]);
50   }
51 
52   /* assigning nodes and elements to sectors */
53 
54   ielset=cs[12];
55   if((*mcs!=1)||(ielset!=0)){
56     for(i=0;i<*nk;i++) inocs[i]=-1;
57   }
58 
59   for(i=0;i<*mcs;i++){
60     is=cs[17*i+4];
61     if(is==1) continue;
62     ielset=cs[17*i+12];
63     if(ielset==0) continue;
64     for(i1=istartset[ielset-1]-1;i1<iendset[ielset-1];i1++){
65       if(ialset[i1]>0){
66         iel=ialset[i1]-1;
67         if(ipkon[iel]<0) continue;
68         indexe=ipkon[iel];
69         if(strcmp1(&lakon[8*iel+3],"2")==0)nope=20;
70         else if (strcmp1(&lakon[8*iel+3],"8")==0)nope=8;
71         else if (strcmp1(&lakon[8*iel+3],"10")==0)nope=10;
72         else if (strcmp1(&lakon[8*iel+3],"4")==0)nope=4;
73         else if (strcmp1(&lakon[8*iel+3],"15")==0)nope=15;
74         else {nope=6;}
75         for(i2=0;i2<nope;++i2){
76           node=kon[indexe+i2]-1;
77           inocs[node]=i;
78         }
79       }
80       else{
81         iel=ialset[i1-2]-1;
82         do{
83           iel=iel-ialset[i1];
84           if(iel>=ialset[i1-1]-1) break;
85           if(ipkon[iel]<0) continue;
86           indexe=ipkon[iel];
87           if(strcmp1(&lakon[8*iel+3],"2")==0)nope=20;
88           else if (strcmp1(&lakon[8*iel+3],"8")==0)nope=8;
89           else if (strcmp1(&lakon[8*iel+3],"10")==0)nope=10;
90           else if (strcmp1(&lakon[8*iel+3],"4")==0)nope=4;
91           else if (strcmp1(&lakon[8*iel+3],"15")==0)nope=15;
92           else {nope=6;}
93           for(i2=0;i2<nope;++i2){
94             node=kon[indexe+i2]-1;
95             inocs[node]=i;
96           }
97         }while(1);
98       }
99     }
100   }
101 
102   /* duplicating triangular faces
103      only those faces are duplicated the nodes of which belong to
104      the same cyclic symmetry. non-integer cyclic symmety numbers are
105      reduced to the next lower integer. */
106 
107   *ntrit=nsegments**ntri;
108   RENEW(kontri,ITG,4**ntrit);
109   for(i=4**ntri;i<4**ntrit;i++) kontri[i]=0;
110 
111   for(i=0;i<*ntri;i++){
112     node1=kontri[4*i];
113     if(inocs[node1-1]<0) continue;
114     idtie=inocs[node1-1];
115     node2=kontri[4*i+1];
116     if((inocs[node2-1]<0)||(inocs[node2-1]!=idtie)) continue;
117     node3=kontri[4*i+2];
118     if((inocs[node3-1]<0)||(inocs[node3-1]!=idtie)) continue;
119     idtie=cs[17*idtie];
120     for(k=1;k<idtie;k++){
121       j=i+k**ntri;
122       kontri[4*j]=node1+k**nk;
123       kontri[4*j+1]=node2+k**nk;
124       kontri[4*j+2]=node3+k**nk;
125       kontri[4*j+3]=kontri[4*i+3];
126     }
127   }
128 
129   RENEW(co,double,3**nk*nsegments);
130   RENEW(vold,double,mt**nk*nsegments);
131   nkt=*nk*nsegments;
132 
133   /* generating the coordinates for the other sectors */
134 
135   icntrl=1;
136 
137   FORTRAN(rectcyl,(co,v,fn,stn,qfn,een,cs,nk,&icntrl,t,filab,&imag,mi,emn));
138 
139   for(jj=0;jj<*mcs;jj++){
140     is=(ITG)(cs[17*jj]);
141     for(i=1;i<is;i++){
142 
143       theta=i*2.*pi/cs[17*jj];
144 
145       for(l=0;l<*nk;l++){
146         if(inocs[l]==jj){
147 	  co[3*l+i*3**nk]=co[3*l];
148 	  co[1+3*l+i*3**nk]=co[1+3*l]-theta;
149 	  co[2+3*l+i*3**nk]=co[2+3*l];
150         }
151       }
152     }
153   }
154 
155   icntrl=-1;
156 
157   FORTRAN(rectcyl,(co,vt,fnt,stnt,qfnt,eent,cs,&nkt,&icntrl,t,filab,
158 		   &imag,mi,emnt));
159 
160   *kontrip=kontri;*cop=co;*voldp=vold;
161 
162   return;
163 }
164 
165