CalculiX  2.13
A Free Software Three-Dimensional Structural Finite Element Program
Functions/Subroutines
generatecycmpcs.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine generatecycmpcs (tolloc, co, nk, ipompc, nodempc, coefmpc, nmpc, ikmpc, ilmpc, mpcfree, rcs, zcs, ics, nr, nz, rcs0, zcs0, labmpc, mcs, triangulation, csab, xn, yn, zn, phi, noded, ncsnodes, rcscg, rcs0cg, zcscg, zcs0cg, nrcg, nzcg, jcs, lcs, kontri, straight, ne, ipkon, kon, lakon, ifacetet, inodface, ncounter, jobnamec, vold, cfd, mi, indepset)
 

Function/Subroutine Documentation

◆ generatecycmpcs()

subroutine generatecycmpcs ( real*8  tolloc,
real*8, dimension(3,*)  co,
integer  nk,
integer, dimension(*)  ipompc,
integer, dimension(3,*)  nodempc,
real*8, dimension(*)  coefmpc,
integer  nmpc,
integer, dimension(*)  ikmpc,
integer, dimension(*)  ilmpc,
integer  mpcfree,
real*8, dimension(*)  rcs,
real*8, dimension(*)  zcs,
integer, dimension(*)  ics,
integer, dimension(*)  nr,
integer, dimension(*)  nz,
real*8, dimension(*)  rcs0,
real*8, dimension(*)  zcs0,
character*20, dimension(*)  labmpc,
integer  mcs,
logical  triangulation,
real*8, dimension(7)  csab,
real*8  xn,
real*8  yn,
real*8  zn,
real*8  phi,
integer  noded,
integer  ncsnodes,
real*8, dimension(*)  rcscg,
real*8, dimension(*)  rcs0cg,
real*8, dimension(*)  zcscg,
real*8, dimension(*)  zcs0cg,
integer, dimension(*)  nrcg,
integer, dimension(*)  nzcg,
integer, dimension(*)  jcs,
integer, dimension(*)  lcs,
integer, dimension(3,*)  kontri,
real*8, dimension(9,*)  straight,
integer  ne,
integer, dimension(*)  ipkon,
integer, dimension(*)  kon,
character*8, dimension(*)  lakon,
integer, dimension(*)  ifacetet,
integer, dimension(*)  inodface,
integer  ncounter,
character*132, dimension(*)  jobnamec,
real*8, dimension(0:mi(2),*)  vold,
integer  cfd,
integer, dimension(*)  mi,
character*81  indepset 
)
26 !
27 ! generate cyclic mpc's
28 !
29  implicit none
30 !
31  logical triangulation,interpolation,multistage
32 !
33  character*1 c
34  character*3 m1,m2,m3
35  character*5 p0,p1,p2,p3,p7,p9999
36  character*8 lakon(*)
37  character*20 labmpc(*),label
38  character*81 indepset
39  character*132 jobnamec(*),fntria
40 !
41  integer ipompc(*),nodempc(3,*),nneigh,ne,ipkon(*),kon(*),
42  & j,k,nk,nmpc,mpcfree,ics(*),nterms,ncyclicsymmetrymodel,
43  & nr(*),nz(*),noded,nodei,ikmpc(*),ilmpc(*),kontri(3,*),
44  & number,idof,ndir,node,ncsnodes,id,mpcfreeold,
45  & mcs,nrcg(*),nzcg(*),jcs(*),lcs(*),nodef(8),
46  & netri,ifacetet(*),inodface(*),lathyp(3,6),inum,one,i,
47  & noden(10),ncounter,ier,ipos,cfd,mi(*),ilen
48 !
49  real*8 tolloc,co(3,*),coefmpc(*),rcs(*),zcs(*),rcs0(*),zcs0(*),
50  & csab(7),xn,yn,zn,xap,yap,zap,rp,zp,al(3,3),ar(3,3),phi,
51  & x2,y2,z2,x3,y3,z3,rcscg(*),rcs0cg(*),zcscg(*),zcs0cg(*),
52  & straight(9,*),ratio(8),vold(0:mi(2),*)
53 !
54  save netri,ncyclicsymmetrymodel
55 !
56  data ncyclicsymmetrymodel /0/
57 !
58 ! latin hypercube positions in a 3 x 3 matrix
59 !
60  data lathyp /1,2,3,1,3,2,2,1,3,2,3,1,3,1,2,3,2,1/
61 !
62  multistage=.false.
63  nneigh=10
64 !
65  xap=co(1,noded)-csab(1)
66  yap=co(2,noded)-csab(2)
67  zap=co(3,noded)-csab(3)
68 !
69  zp=xap*xn+yap*yn+zap*zn
70  rp=dsqrt((xap-zp*xn)**2+(yap-zp*yn)**2+(zap-zp*zn)**2)
71 !
72  call near2d(rcs0,zcs0,rcs,zcs,nr,nz,rp,zp,ncsnodes,noden,nneigh)
73  node=noden(1)
74  nodei=abs(ics(noden(1)))
75 !
76 ! check whether node is on axis
77 !
78  if(nodei.eq.noded) then
79  return
80  endif
81 !
82  interpolation=.false.
83 !
84  if(rp.gt.1.d-10) then
85  x2=(xap-zp*xn)/rp
86  y2=(yap-zp*yn)/rp
87  z2=(zap-zp*zn)/rp
88  x3=yn*z2-y2*zn
89  y3=x2*zn-xn*z2
90  z3=xn*y2-x2*yn
91  endif
92 !
93  if((tolloc.ge.0.d0).and.
94  & (tolloc.le.dsqrt((rp-rcs0(node))**2+(zp-zcs0(node))**2)))
95  & then
96 !
97 ! the nodal positions on the dependent and independent
98 ! sides of the mpc's do no agree: interpolation is
99 ! necessary.
100 !
101 c write(*,*) '*WARNING in generatecycmpcs: no cyclic'
102 c write(*,*) ' symmetric partner found for'
103 c write(*,*) ' dependent node ',noded,'.'
104 c write(*,*) ' allowed tolerance:',tolloc
105 c write(*,*) ' best partner node number:',nodei
106 c write(*,*) ' actual distance in a radial plane: ',
107 c & dsqrt((rp-rcs0(node))**2+(zp-zcs0(node))**2)
108 c write(*,*) ' Remedy: the node is connected to an'
109 c write(*,*) ' independent element side.'
110 c write(*,*)
111 !
112  interpolation=.true.
113 !
114  if(.not.triangulation) then
115  call triangulate(ics,rcs0,zcs0,ncsnodes,
116  & rcscg,rcs0cg,zcscg,zcs0cg,nrcg,nzcg,jcs,kontri,
117  & straight,ne,ipkon,kon,lakon,lcs,netri,ifacetet,
118  & inodface)
119  triangulation=.true.
120 !
121 c fntria(1:28)='TriMasterCyclicSymmetryModel'
122 c ncyclicsymmetrymodel=ncyclicsymmetrymodel+1
123 c if(ncyclicsymmetrymodel.lt.10) then
124 c write(fntria(29:29),'(i1)')ncyclicsymmetrymodel
125 c fntria(30:33)='.frd'
126 c ipos=34
127 c elseif(ncyclicsymmetrymodel.lt.100) then
128 c write(fntria(29:30),'(i2)')ncyclicsymmetrymodel
129 c fntria(31:34)='.frd'
130 c ipos=35
131 c else
132 c write(*,*) '*ERROR in generatecycmpcs: no more than'
133 c write(*,*) ' 99 cyclic symmetry model cards'
134 c write(*,*) ' allowed'
135 c call exit(201)
136 c endif
137 c do i=ipos,132
138 c fntria(i:i)=' '
139 c enddo
140 c
141  ilen=index(indepset,' ')
142  fntria(1:3)='Tri'
143  do j=4,ilen+2
144  fntria(j:j)=indepset(j-3:j-3)
145  enddo
146  fntria(ilen+3:ilen+6)='.frd'
147  do j=ilen+7,132
148  fntria(j:j)=' '
149  enddo
150 !
151 c open(70,file=fntria,status='unknown')
152 c c='C'
153 c m1=' -1'
154 c m2=' -2'
155 c m3=' -3'
156 c p0=' 0'
157 c p1=' 1'
158 c p2=' 2'
159 c p3=' 3'
160 c p7=' 7'
161 c p9999=' 9999'
162 c one=1
163 c write(70,'(a5,a1)') p1,c
164 c write(70,'(a5,a1,67x,i1)') p2,c,one
165 c do i=1,nk
166 c write(70,'(a3,i10,1p,3e12.5)') m1,i,(co(j,i),j=1,3)
167 c enddo
168 c write(70,'(a3)') m3
169 c write(70,'(a5,a1,67x,i1)') p3,c,one
170 c do i=1,netri
171 c write(70,'(a3,i10,2a5)')m1,i,p7,p0
172 c write(70,'(a3,3i10)') m2,(kontri(j,i),j=1,3)
173 c enddo
174 c write(70,'(a3)') m3
175 c write(70,'(a5)') p9999
176 c close(70)
177 !
178  endif
179 !
180  label='CYCLIC '
181  if(mcs.lt.10) then
182  write(label(7:7),'(i1)') mcs
183  elseif(mcs.lt.100) then
184  write(label(7:8),'(i2)') mcs
185  else
186  write(*,*)'*ERROR in generatecycmpcs: no more than 99'
187  write(*,*)' cyclic symmetry definitions allowed'
188  call exit(201)
189  endif
190 !
191  nodei=nk+1
192 !
193 ! copying the initial conditions for the new node
194 !
195  do i=0,mi(2)
196  vold(i,nodei)=vold(i,noded)
197  enddo
198 !
199  co(1,nodei)=csab(1)+zp*xn+rp*(x2*dcos(phi)+x3*dsin(phi))
200  co(2,nodei)=csab(2)+zp*yn+rp*(y2*dcos(phi)+y3*dsin(phi))
201  co(3,nodei)=csab(3)+zp*zn+rp*(z2*dcos(phi)+z3*dsin(phi))
202 !
203  ier=0
204 !
205  call linkdissimilar(co,csab,
206  & rcscg,rcs0cg,zcscg,zcs0cg,nrcg,nzcg,straight,
207  & nodef,ratio,nterms,rp,zp,netri,
208  & nodei,ifacetet,inodface,noded,xn,yn,
209  & zn,ier,multistage)
210 !
211  if(ier.ne.0) then
212  ncounter=ncounter+1
213  return
214  endif
215 !
216  else
217  if(ics(node).lt.0) return
218  endif
219 !
220 ! generating the mechanical MPC's; the generated MPC's are for
221 ! nodal diameter 0. For other nodal diameters the MPC's are
222 ! changed implicitly in mastructcs and mafillsmcs
223 !
224  call transformatrix(csab,co(1,noded),al)
225  call transformatrix(csab,co(1,nodei),ar)
226 !
227 ! checking for latin hypercube positions in matrix al none of
228 ! which are zero
229 !
230  do inum=1,6
231  if((dabs(al(lathyp(1,inum),1)).gt.1.d-3).and.
232  & (dabs(al(lathyp(2,inum),2)).gt.1.d-3).and.
233  & (dabs(al(lathyp(3,inum),3)).gt.1.d-3)) exit
234  enddo
235 !
236  do ndir=1,3
237 !
238 ! determining which direction to use for the
239 ! dependent side: should not occur on the dependent
240 ! side in another MPC and should have a nonzero
241 ! coefficient
242 !
243  number=lathyp(ndir,inum)
244  idof=8*(noded-1)+number
245  call nident(ikmpc,idof,nmpc,id)
246  if(id.gt.0) then
247  if(ikmpc(id).eq.idof) then
248  write(*,*) '*WARNING in generatecycmpcs: cyclic MPC in no
249  &de'
250  write(*,*) ' ',noded,' and direction ',ndir
251  write(*,*) ' cannot be created: the'
252  write(*,*) ' DOF in this node is already used'
253  cycle
254  endif
255  endif
256  number=number-1
257 !
258  nmpc=nmpc+1
259  labmpc(nmpc)='CYCLIC '
260  if(mcs.lt.10) then
261  write(labmpc(nmpc)(7:7),'(i1)') mcs
262  elseif(mcs.lt.100) then
263  write(labmpc(nmpc)(7:8),'(i2)') mcs
264  else
265  write(*,*)'*ERROR in generatecycmpcs: no more than 99'
266  write(*,*)' cyclic symmetry definitions allowed'
267  call exit(201)
268  endif
269  ipompc(nmpc)=mpcfree
270 !
271 ! updating ikmpc and ilmpc
272 !
273  do j=nmpc,id+2,-1
274  ikmpc(j)=ikmpc(j-1)
275  ilmpc(j)=ilmpc(j-1)
276  enddo
277  ikmpc(id+1)=idof
278  ilmpc(id+1)=nmpc
279 !
280  do j=1,3
281  number=number+1
282  if(number.gt.3) number=1
283  if(dabs(al(number,ndir)).lt.1.d-5) cycle
284  nodempc(1,mpcfree)=noded
285  nodempc(2,mpcfree)=number
286  coefmpc(mpcfree)=al(number,ndir)
287  mpcfree=nodempc(3,mpcfree)
288  if(mpcfree.eq.0) then
289  write(*,*)'*ERROR in generatecycmpcs: increase memmpc_'
290  call exit(201)
291  endif
292  enddo
293  do j=1,3
294  number=number+1
295  if(number.gt.3) number=1
296  if(dabs(ar(number,ndir)).lt.1.d-5) cycle
297  if(.not.interpolation) then
298  nodempc(1,mpcfree)=nodei
299  nodempc(2,mpcfree)=number
300  coefmpc(mpcfree)=-ar(number,ndir)
301  mpcfreeold=mpcfree
302  mpcfree=nodempc(3,mpcfree)
303  if(mpcfree.eq.0) then
304  write(*,*)
305  & '*ERROR in generatecycmpcs: increase memmpc_'
306  call exit(201)
307  endif
308  else
309  do k=1,nterms
310  nodempc(1,mpcfree)=nodef(k)
311  nodempc(2,mpcfree)=number
312  coefmpc(mpcfree)=-ar(number,ndir)*ratio(k)
313  mpcfreeold=mpcfree
314  mpcfree=nodempc(3,mpcfree)
315  if(mpcfree.eq.0) then
316  write(*,*) '*ERROR in generatecycmpcs: increase nmp
317  &c_'
318  call exit(201)
319  endif
320  enddo
321  endif
322  enddo
323  nodempc(3,mpcfreeold)=0
324  enddo
325 !
326 ! generating the thermal MPC's; the generated MPC's are for nodal
327 ! diameter 0. BETTER: based on ithermal(2), cf. gen3dfrom2d.f
328 !
329  nmpc=nmpc+1
330  labmpc(nmpc)='CYCLIC '
331  if(mcs.lt.10) then
332  write(labmpc(nmpc)(7:7),'(i1)') mcs
333  elseif(mcs.lt.100) then
334  write(labmpc(nmpc)(7:8),'(i2)') mcs
335  else
336  write(*,*)'*ERROR in generatecycmpcs: no more than 99'
337  write(*,*)' cyclic symmetry definitions allowed'
338  call exit(201)
339  endif
340  ipompc(nmpc)=mpcfree
341  idof=8*(noded-1)
342  call nident(ikmpc,idof,nmpc-1,id)
343  if(id.gt.0) then
344  if(ikmpc(id).eq.idof) then
345  write(*,*) '*ERROR in generatecycmpcs: temperature'
346  write(*,*) ' in node',noded,'is already used'
347  call exit(201)
348  endif
349  endif
350 !
351 ! updating ikmpc and ilmpc
352 !
353  do j=nmpc,id+2,-1
354  ikmpc(j)=ikmpc(j-1)
355  ilmpc(j)=ilmpc(j-1)
356  enddo
357  ikmpc(id+1)=idof
358  ilmpc(id+1)=nmpc
359 !
360  nodempc(1,mpcfree)=noded
361  nodempc(2,mpcfree)=0
362  coefmpc(mpcfree)=1.d0
363  mpcfree=nodempc(3,mpcfree)
364  if(mpcfree.eq.0) then
365  write(*,*)'*ERROR in generatecycmpcs: increase memmpc_'
366  call exit(201)
367  endif
368  if(.not.interpolation) then
369  nodempc(1,mpcfree)=nodei
370  nodempc(2,mpcfree)=0
371  coefmpc(mpcfree)=-1.d0
372  mpcfreeold=mpcfree
373  mpcfree=nodempc(3,mpcfree)
374  if(mpcfree.eq.0) then
375  write(*,*)'*ERROR in generatecycmpcs: increase memmpc_'
376  call exit(201)
377  endif
378  else
379  do k=1,nterms
380  nodempc(1,mpcfree)=nodef(k)
381  nodempc(2,mpcfree)=0
382  coefmpc(mpcfree)=-ratio(k)
383  mpcfreeold=mpcfree
384  mpcfree=nodempc(3,mpcfree)
385  if(mpcfree.eq.0) then
386  write(*,*)'*ERROR in generatecycmpcs: increase memmpc_'
387  call exit(201)
388  endif
389  enddo
390  endif
391  nodempc(3,mpcfreeold)=0
392 !
393 ! generating the static pressure MPC's for 3D fluid calculations;
394 ! the generated MPC's are for nodal diameter 0.
395 !
396  if(cfd.eq.1) then
397  nmpc=nmpc+1
398  labmpc(nmpc)='CYCLIC '
399  if(mcs.lt.10) then
400  write(labmpc(nmpc)(7:7),'(i1)') mcs
401  elseif(mcs.lt.100) then
402  write(labmpc(nmpc)(7:8),'(i2)') mcs
403  else
404  write(*,*)'*ERROR in generatecycmpcs: no more than 99'
405  write(*,*)' cyclic symmetry definitions allowed'
406  call exit(201)
407  endif
408  ipompc(nmpc)=mpcfree
409  idof=8*(noded-1)+4
410  call nident(ikmpc,idof,nmpc-1,id)
411  if(id.gt.0) then
412  if(ikmpc(id).eq.idof) then
413  write(*,*) '*ERROR in generatecycmpcs: pressure'
414  write(*,*) ' in node',noded,'is already used'
415  call exit(201)
416  endif
417  endif
418 !
419 ! updating ikmpc and ilmpc
420 !
421  do j=nmpc,id+2,-1
422  ikmpc(j)=ikmpc(j-1)
423  ilmpc(j)=ilmpc(j-1)
424  enddo
425  ikmpc(id+1)=idof
426  ilmpc(id+1)=nmpc
427 !
428  nodempc(1,mpcfree)=noded
429  nodempc(2,mpcfree)=4
430  coefmpc(mpcfree)=1.d0
431  mpcfree=nodempc(3,mpcfree)
432  if(mpcfree.eq.0) then
433  write(*,*)'*ERROR in generatecycmpcs: increase memmpc_'
434  call exit(201)
435  endif
436  if(.not.interpolation) then
437  nodempc(1,mpcfree)=nodei
438  nodempc(2,mpcfree)=4
439  coefmpc(mpcfree)=-1.d0
440  mpcfreeold=mpcfree
441  mpcfree=nodempc(3,mpcfree)
442  if(mpcfree.eq.0) then
443  write(*,*)'*ERROR in generatecycmpcs: increase memmpc_'
444  call exit(201)
445  endif
446  else
447  do k=1,nterms
448  nodempc(1,mpcfree)=nodef(k)
449  nodempc(2,mpcfree)=4
450  coefmpc(mpcfree)=-ratio(k)
451  mpcfreeold=mpcfree
452  mpcfree=nodempc(3,mpcfree)
453  if(mpcfree.eq.0) then
454  write(*,*)
455  & '*ERROR in generatecycmpcs: increase memmpc_'
456  call exit(201)
457  endif
458  enddo
459  endif
460  nodempc(3,mpcfreeold)=0
461  endif
462 !
463  return
triangulate
subroutine triangulate(ics, rcs0, zcs0, ncsnodes, rcscg, rcs0cg, zcscg, zcs0cg, nrcg, nzcg, jcs, kontri, straight, ne, ipkon, kon, lakon, lcs, netri, ifacetet, inodface)
Definition: triangulate.f:22
transformatrix
subroutine transformatrix(xab, p, a)
Definition: transformatrix.f:20
nident
subroutine nident(x, px, n, id)
Definition: nident.f:26
linkdissimilar
subroutine linkdissimilar(co, csab, rcscg, rcs0cg, zcscg, zcs0cg, nrcg, nzcg, straight, nodef, ratio, nterms, rp, zp, netri, nodei, ifacetet, inodface, noded, xn, yn, zn, ier, multistage)
Definition: linkdissimilar.f:23
near2d
subroutine near2d(xo, yo, x, y, nx, ny, xp, yp, n, neighbor, k)
Definition: near2d.f:20
Hosted by OpenAircraft.com, (Michigan UAV, LLC)