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!
20!     Subroutine x_interpolate.f
21!
22!     Triangulates the face and interpolates xstate variables according to
23!     the plane equations of these planes
24!
25!     by: Jaro Hokkanen
26!
27!
28      subroutine interpolateinface(kk,xstate,xstateini,numpts,nstate_,
29     &     mi,islavsurf,pslavsurf,
30     &     ne0,islavsurfold,pslavsurfold)
31!
32      implicit none
33!
34      integer numpts,i_int,n_int,i,j,k,kk,l,ll,nn,
35     &     koncont(3,2*numpts+1),itri,kflag,neigh(1),kneigh,
36     &     imastop(3,2*numpts+1),indexcj,nopespringj,list(numpts),
37     &     igauss,mi(*),nstate_,itriangle(100),itriold,
38     &     ifaceq(8,6),ip(numpts),ne0,itrinew,ntriangle,
39     &     ifacet(6,4),ifacew1(4,5),ifacew2(8,5),n,islavsurf(2,*),
40     &     ibin(numpts),ivert1,ntriangle_,nterms,m,islavsurfold(2,*),
41     &     nx(2*numpts+1),ny(2*numpts+1),isol,id
42!
43      real*8 xstate(nstate_,mi(1),*),p(3),pslavsurfold(3,*),
44     &     xstateini(nstate_,mi(1),*),coi(2,numpts+3),pneigh(3,3),
45     &     cg(2,2*numpts+1),pslavsurf(3,*),xil,etl,ratio(3),
46     &     z(9,3),x(2*numpts+1),xo(2*numpts+1),dis(3),
47     &     y(2*numpts+1),yo(2*numpts+1),straight(9,2*numpts+1),dist
48!
49!     nodes per face for hex elements
50!
51      data ifaceq /4,3,2,1,11,10,9,12,
52     &     5,6,7,8,13,14,15,16,
53     &     1,2,6,5,9,18,13,17,
54     &     2,3,7,6,10,19,14,18,
55     &     3,4,8,7,11,20,15,19,
56     &     4,1,5,8,12,17,16,20/
57!
58!     nodes per face for tet elements
59!
60      data ifacet /1,3,2,7,6,5,
61     &     1,2,4,5,9,8,
62     &     2,3,4,6,10,9,
63     &     1,4,3,8,10,7/
64!
65!     nodes per face for linear wedge elements
66!
67      data ifacew1 /1,3,2,0,
68     &     4,5,6,0,
69     &     1,2,5,4,
70     &     2,3,6,5,
71     &     3,1,4,6/
72!
73!     nodes per face for quadratic wedge elements
74!
75      data ifacew2 /1,3,2,9,8,7,0,0,
76     &     4,5,6,10,11,12,0,0,
77     &     1,2,5,4,7,14,10,13,
78     &     2,3,6,5,8,15,11,14,
79     &     3,1,4,6,9,13,12,15/
80!
81      kneigh=1
82!
83      do i=1,numpts
84        list(i)=i
85      enddo
86!
87!     Loop over the old integration points within the face
88!
89      ll=0
90      do l=islavsurfold(2,kk)+1,islavsurfold(2,kk+1)
91        ll=ll+1
92        coi(1,ll)=pslavsurfold(1,l)
93        coi(2,ll)=pslavsurfold(2,l)
94        ip(ll)=ne0+l
95      enddo
96!
97!     Calling deltri for triangulation
98!
99      kflag=2
100!
101      call deltri(numpts,numpts,coi(1,1:numpts+3),coi(2,1:numpts+3),
102     &     list,ibin,koncont,imastop,n)
103!
104!     rearranging imastop field: imastop(1,i) is the triangle opposite
105!     of local node 1 in triangle i etc..
106!
107      nn=0
108      do i=1,n
109        ivert1=imastop(1,i)
110        imastop(1,i)=imastop(2,i)
111        imastop(2,i)=imastop(3,i)
112        imastop(3,i)=ivert1
113      enddo
114!
115!     determining the center of gravity and the bounding planes
116!     of the triangles
117!
118      call updatecont2d(koncont,n,coi,cg,straight)
119!
120!     sorting the centers of gravity
121!
122      do l=1,n
123        xo(l)=cg(1,l)
124        x(l)=xo(l)
125        nx(l)=l
126        yo(l)=cg(2,l)
127        y(l)=yo(l)
128        ny(l)=l
129      enddo
130      kflag=2
131      call dsort(x,nx,n,kflag)
132      call dsort(y,ny,n,kflag)
133!
134!     Loop over the new integration points
135!
136      do igauss=islavsurf(2,kk)+1,islavsurf(2,kk+1)
137!
138!     coordinates of the new integration point
139!
140        p(1)=pslavsurf(1,igauss)
141        p(2)=pslavsurf(2,igauss)
142!
143!     closest triangle center of gravity
144!
145        call near2d(xo,yo,x,y,nx,ny,p(1),p(2),
146     &       n,neigh,kneigh)
147!
148        isol=0
149!
150        itriold=0
151        itri=neigh(1)
152        ntriangle=0
153        ntriangle_=100
154!
155        loop1: do
156        do l=1,3
157          ll=3*l-2
158          dist=straight(ll,itri)*p(1)+
159     &         straight(ll+1,itri)*p(2)+
160     &         straight(ll+2,itri)
161!
162          if(dist.gt.0.d0) then
163            itrinew=imastop(l,itri)
164            if(itrinew.eq.0) then
165c     write(*,*) '**border reached'
166              exit loop1
167            elseif(itrinew.eq.itriold) then
168c     write(*,*) '**solution in between triangles'
169              isol=itri
170              exit loop1
171            else
172              call nident(itriangle,itrinew,ntriangle,id)
173              if(id.gt.0) then
174                if(itriangle(id).eq.itrinew) then
175c     write(*,*) '**circular path;no solution'
176                  exit loop1
177                endif
178              endif
179              ntriangle=ntriangle+1
180              if(ntriangle.gt.ntriangle_) then
181c     write(*,*) '**too many iterations'
182                exit loop1
183              endif
184              do k=ntriangle,id+2,-1
185                itriangle(k)=itriangle(k-1)
186              enddo
187              itriangle(id+1)=itrinew
188              itriold=itri
189              itri=itrinew
190              cycle loop1
191            endif
192          elseif(l.eq.3) then
193c     write(*,*) '**regular solution'
194            isol=itri
195            exit loop1
196          endif
197        enddo
198      enddo loop1
199!
200      if(isol.eq.0) then
201!
202!     mapping the new integration point onto the nearest
203!     triangle
204!
205!     calculating the distance from each edge of the triangle
206!     dis(i) is the distance of the new integration point fromt
207!     the edge opposite to local node i of the triangle itri
208!
209        do l=1,3
210          ll=3*l-2
211          dis(l)=straight(ll,itri)*p(1)+
212     &         straight(ll+1,itri)*p(2)+
213     &         straight(ll+2,itri)
214        enddo
215!
216        if(dis(1).gt.0.d0) then
217          if(dis(2).gt.0.d0) then
218!
219!     projection is node 3 of the triangle
220!
221            p(1)=coi(1,koncont(3,itri))
222            p(2)=coi(2,koncont(3,itri))
223          elseif(dis(3).gt.0.d0) then
224!
225!     projection is node 2 of the triangle
226!
227            p(1)=coi(1,koncont(2,itri))
228            p(2)=coi(2,koncont(2,itri))
229          else
230!
231!     projection onto local edge 1 (opposite to local node 1)
232!
233            p(1)=p(1)-dis(1)*straight(1,itri)
234            p(2)=p(2)-dis(1)*straight(2,itri)
235          endif
236        elseif(dis(2).gt.0.d0) then
237          if(dis(3).gt.0.d0) then
238!
239!     projection is node 1 of the triangle
240!
241            p(1)=coi(1,koncont(1,itri))
242            p(2)=coi(2,koncont(1,itri))
243          else
244!
245!     projection onto local edge 2 (opposite to local node 2)
246!
247            p(1)=p(1)-dis(2)*straight(4,itri)
248            p(2)=p(2)-dis(2)*straight(5,itri)
249          endif
250        else
251!
252!     projection onto local edge 3 (opposite to local node 3)
253!
254          p(1)=p(1)-dis(3)*straight(7,itri)
255          p(2)=p(2)-dis(3)*straight(8,itri)
256        endif
257c     do k=1,3
258c     do m=1,2
259c     pneigh(m,k)=coi(m,koncont(k,itri))
260c     enddo
261c     pneigh(3,k)=0.d0
262c     enddo
263c     p(3)=0.d0
264c     nterms=3
265c     !
266c     call attach_2d(pneigh,p,nterms,ratio,dist,xil,etl)
267c     !
268c     do m=1,2
269c     p(m)=0.d0
270c     do k=1,3
271c     p(m)=p(m)+ratio(k)*pneigh(m,k)
272c     enddo
273c     enddo
274c
275      endif
276!
277!     Assigning xstate values from the vertices of the triangle
278!     to the field z (integration point values from xstateini)
279!
280      do k=1,3
281        do i=1,9
282          z(i,k)=xstateini(i,1,ip(koncont(k,itri)))
283        enddo
284      enddo
285!
286!     Calling plane_eq to interpolate values using plane equation
287!
288      do i=1,9
289        call plane_eq(
290     &       coi(1,koncont(1,itri)),coi(2,koncont(1,itri)),z(i,1),
291     &       coi(1,koncont(2,itri)),coi(2,koncont(2,itri)),z(i,2),
292     &       coi(1,koncont(3,itri)),coi(2,koncont(3,itri)),z(i,3),
293     &       p(1),p(2),xstate(i,1,ne0+igauss))
294      enddo
295!
296      enddo
297!
298      return
299      end
300
301