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 subroutine calculateh(nk,v,veold,stn,een,emn,epn,enern,qfn, 20 & errn,h,filab,mi,d,nh,dmin,ipoed,iedg,cotet,jfix) 21! 22! calculating the desired size of h in all nodes of the 23! original mesh 24! 25 implicit none 26! 27 character*4 label 28 character*87 filab(*) 29! 30 integer mi(*),i,nk,istat,nh(*),ipoed(*),iedg(3,*),index,n1,n2, 31 & jfix(*) 32! 33 real*8 v(0:mi(2),*),veold(0:mi(2),*),stn(6,*),een(6,*),emn(6,*), 34 & epn(*),enern(*),qfn(3,*),errn(6,*),h(*),targetsize,size,d(*), 35 & dmin,cotet(3,*) 36! 37! 38! 39 label(1:4)=filab(48)(3:6) 40! 41 read(filab(48)(7:26),'(f20.0)',iostat=istat) targetsize 42 if(istat.gt.0) then 43 write(*,*) '*ERROR in calculateh:' 44 write(*,*) ' targetsize not readable' 45 call exit(201) 46 endif 47! 48! determine the size of all edges 49! 50 dmin=1.d30 51! 52 loop: do i=1,nk 53 index=ipoed(i) 54 do 55 if(index.eq.0) cycle loop 56! 57 n1=iedg(1,index) 58 n2=iedg(2,index) 59! 60 d(index)=dsqrt((cotet(1,n1)-cotet(1,n2))**2+ 61 & (cotet(2,n1)-cotet(2,n2))**2+ 62 & (cotet(3,n1)-cotet(3,n2))**2) 63! 64 h(n1)=h(n1)+d(index) 65 h(n2)=h(n2)+d(index) 66 nh(n1)=nh(n1)+1 67 nh(n2)=nh(n2)+1 68! 69 if(d(index).lt.dmin) dmin=d(index) 70! 71 index=iedg(3,index) 72 enddo 73 enddo loop 74! 75 do i=1,nk 76! 77! take the mean at each node 78! 79 if(nh(i).le.0) cycle 80 h(i)=h(i)/nh(i) 81! 82! if the node is on the boundary with part of the mesh which 83! is not being refined: h is set to the mean edge length 84! 85 if(jfix(i).eq.1) cycle 86! 87 if(label.eq.'U ') then 88 size=dsqrt(v(1,i)**2+v(2,i)**2+v(3,i)**2) 89 elseif(label.eq.'V ') then 90 size=dsqrt(veold(1,i)**2+veold(2,i)**2+veold(3,i)**2) 91 elseif(label.eq.'S ') then 92 size=dsqrt(stn(1,i)**2+stn(2,i)**2+stn(3,i)**2+ 93 & 2.d0*(stn(4,i)**2+stn(5,i)**2+stn(6,i)**2)) 94 elseif(label.eq.'E ') then 95 size=dsqrt(een(1,i)**2+een(2,i)**2+een(3,i)**2+ 96 & 2.d0*(een(4,i)**2+een(5,i)**2+een(6,i)**2)) 97 elseif(label.eq.'ME ') then 98 size=dsqrt(emn(1,i)**2+emn(2,i)**2+emn(3,i)**2+ 99 & 2.d0*(emn(4,i)**2+emn(5,i)**2+emn(6,i)**2)) 100 elseif(label.eq.'PEEQ') then 101 size=dabs(epn(i)) 102 elseif(label.eq.'ENER') then 103 size=dabs(enern(i)) 104 elseif(label.eq.'HFL ') then 105 size=dsqrt(qfn(1,i)**2+qfn(2,i)**2+qfn(3,i)**2) 106 elseif(label.eq.'ERR ') then 107 size=dabs(errn(1,i)) 108 elseif(label.eq.'USER') then 109c call ucalculateh(v,veold,stn,een,emn,epn,enern,qfn, 110c & errn,size,mi) 111 endif 112! 113 if(size/targetsize.gt.1.d0) then 114 h(i)=h(i)*targetsize/size 115 endif 116 enddo 117! 118 return 119 end 120