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 attach_3d(pneigh,pnode,nterms,ratio,dist,xil,etl,zel,
20     &  loopa)
21!
22!     ataches node with coordinates in "pnode" to the face containing
23!     "nterms" nodes with coordinates in field "pneigh" (nterms < 20).
24!     cave: the coordinates are stored in pneigh(1..3,*)
25!
26      implicit none
27!
28      integer nterms,i,j,k,m,imin,jmin,kmin,im,jm,km,n,loopa
29!
30      real*8 ratio(20),pneigh(3,20),pnode(3),a,xi(-1:1,-1:1,-1:1),
31     &  et(-1:1,-1:1,-1:1),ze(-1:1,-1:1,-1:1),p(3),distmin,d1,dist,
32     &  xil,etl,zel,dx(3),al,dummy
33!
34      n=3
35!
36      d1=1.d0
37!
38      xi(0,0,0)=0.d0
39      et(0,0,0)=0.d0
40      ze(0,0,0)=0.d0
41      call distattach_3d(xi(0,0,0),et(0,0,0),ze(0,0,0),pneigh,pnode,a,p,
42     &     ratio,nterms)
43      distmin=a
44      imin=0
45      jmin=0
46      kmin=0
47!
48      do m=1,loopa
49!
50!     initialisation
51!
52         d1=d1/10.d0
53!
54         do i=-1,1
55            do j=-1,1
56               do k=-1,1
57                  if((i.eq.0).and.(j.eq.0).and.(k.eq.0)) cycle
58!
59                  xi(i,j,k)=xi(0,0,0)+i*d1
60                  et(i,j,k)=et(0,0,0)+j*d1
61                  ze(i,j,k)=ze(0,0,0)+k*d1
62!
63!              check whether inside the (-1,1)x(-1,1)x(-1,1) domain
64!
65                  if((xi(i,j,k).le.1.d0).and.
66     &               (xi(i,j,k).ge.-1.d0).and.
67     &               (et(i,j,k).le.1.d0).and.
68     &               (et(i,j,k).ge.-1.d0).and.
69     &               (ze(i,j,k).le.1.d0).and.
70     &               (ze(i,j,k).ge.-1.d0)) then
71                     call distattach_3d(xi(i,j,k),et(i,j,k),ze(i,j,k),
72     &                    pneigh,pnode,a,p,ratio,nterms)
73!
74!     checking for smallest initial distance
75!
76                     if(a.lt.distmin) then
77                        distmin=a
78                        imin=i
79                        jmin=j
80                        kmin=k
81                     endif
82                  endif
83!
84               enddo
85            enddo
86         enddo
87!
88!     minimizing the distance from the face to the node
89!
90         do
91!
92!     exit if minimum found
93!
94c            write(*,*) 'attach_3d',m,imin,jmin,kmin
95            if((imin.eq.0).and.(jmin.eq.0).and.(kmin.eq.0)) exit
96!
97!           new center of 3x3x3 matrix
98!
99            xi(0,0,0)=xi(imin,jmin,kmin)
100            et(0,0,0)=et(imin,jmin,kmin)
101            ze(0,0,0)=ze(imin,jmin,kmin)
102!
103            im=imin
104            jm=jmin
105            km=kmin
106!
107            imin=0
108            jmin=0
109            kmin=0
110!
111            do i=-1,1
112               do j=-1,1
113                  do k=-1,1
114                     if((i+im.lt.-1).or.(i+im.gt.1).or.
115     &                  (j+jm.lt.-1).or.(j+jm.gt.1).or.
116     &                  (k+km.lt.-1).or.(k+km.gt.1)) then
117!
118                        xi(i,j,k)=xi(0,0,0)+i*d1
119                        et(i,j,k)=et(0,0,0)+j*d1
120                        ze(i,j,k)=ze(0,0,0)+k*d1
121!
122!              check whether inside the (-1,1)x(-1,1)x(-1,1) domain
123!
124                        if((xi(i,j,k).le.1.d0).and.
125     &                     (xi(i,j,k).ge.-1.d0).and.
126     &                     (et(i,j,k).le.1.d0).and.
127     &                     (et(i,j,k).ge.-1.d0).and.
128     &                     (ze(i,j,k).le.1.d0).and.
129     &                     (ze(i,j,k).ge.-1.d0)) then
130                           call distattach_3d(xi(i,j,k),et(i,j,k),
131     &                          ze(i,j,k),pneigh,pnode,a,p,ratio,nterms)
132!
133!                       check for new minimum
134!
135                           if(a.lt.distmin) then
136                              distmin=a
137                              imin=i
138                              jmin=j
139                              kmin=k
140                           endif
141                        endif
142!
143                     endif
144                  enddo
145               enddo
146            enddo
147         enddo
148      enddo
149!
150      call distattach_3d(xi(0,0,0),et(0,0,0),ze(0,0,0),pneigh,pnode,a,p,
151     &     ratio,nterms)
152!
153      do i=1,3
154        pnode(i)=p(i)
155      enddo
156!
157      dist=dsqrt(a)
158!
159      if((nterms.eq.20).or.(nterms.eq.8)) then
160         xil=xi(0,0,0)
161         etl=et(0,0,0)
162         zel=ze(0,0,0)
163      elseif((nterms.eq.4).or.(nterms.eq.10)) then
164         xil=(xi(0,0,0)+1.d0)/2.d0
165         etl=(et(0,0,0)+1.d0)/2.d0
166         zel=(ze(0,0,0)+1.d0)/2.d0
167         dx(1)=xil
168         dx(2)=etl
169         dx(3)=zel
170         call insertsortd(dx,n)
171c         call dsort(dx,iy,n,kflag)
172         if(dx(3).gt.1.d-30) then
173            al=dx(3)/(xil+etl+zel)
174            xil=al*xil
175            etl=al*etl
176            zel=al*zel
177         endif
178      elseif((nterms.eq.6).or.(nterms.eq.15)) then
179         xil=(xi(0,0,0)+1.d0)/2.d0
180         etl=(et(0,0,0)+1.d0)/2.d0
181         if(xil+etl.gt.1.d0) then
182            dummy=xil
183            xil=1.d0-etl
184            etl=1.d0-dummy
185         endif
186         zel=ze(0,0,0)
187      endif
188!
189      return
190      end
191
192