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