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 map3dtolayer(yn,ipkon,kon,lakon,nfield,
20     &  ne,co,ielmat,mi)
21!
22!     interpolates 3d field nodal values to nodal values in the
23!     layers of composite materials
24!
25      implicit none
26!
27      character*8 lakon(*)
28!
29      integer ipkon(*),kon(*),ne,nfield,i,j,k,l,mi(*),nelem,
30     &  ielmat(mi(3),*),nlayer,iflag,nbot20(8),nmid20(8),ntop20(8),
31     &  nbot15(6),nmid15(6),ntop15(6),ibot,itop,nope,nopes,nopeexp,
32     &  indexe,konl(20),imid,node,nodebot,nodetop
33!
34      real*8 yn(nfield,*),shp(4,20),xsj(3),co(3,*),xl(3,20),
35     &  xi,et,ze,dd,dt,xi20(8),et20(8),xi15(6),et15(6)
36!
37!
38!
39      data iflag /1/
40!
41      data nbot20 /1,2,3,4,9,10,11,12/
42      data nmid20 /17,18,19,20,0,0,0,0/
43      data ntop20 /5,6,7,8,13,14,15,16/
44!
45      data nbot15 /1,2,3,7,8,9/
46      data nmid15 /13,14,15,0,0,0/
47      data ntop15 /4,5,6,10,11,12/
48!
49      data xi20 /-1.d0,1.d0,1.d0,-1.d0,0.d0,1.d0,0.d0,-1.d0/
50      data et20 /-1.d0,-1.d0,1.d0,1.d0,-1.d0,0.d0,1.d0,0.d0/
51!
52      data xi15 /0.d0,1.d0,0.d0,0.5d0,0.5d0,0.d0/
53      data et15 /0.d0,0.d0,1.d0,0.d0,0.5d0,0.5d0/
54!
55      include "gauss.f"
56!
57      do nelem=1,ne
58         if(lakon(nelem)(7:8).eq.'LC') then
59!
60!        composite materials
61!
62!        determining the number of layers
63!
64            nlayer=0
65            do k=1,mi(3)
66               if(ielmat(k,nelem).ne.0) then
67                  nlayer=nlayer+1
68               endif
69            enddo
70!
71            indexe=ipkon(nelem)
72!
73            if(lakon(nelem)(4:5).eq.'20') then
74               nope=20
75               nopes=8
76               nopeexp=28
77            elseif(lakon(nelem)(4:5).eq.'15') then
78               nope=15
79               nopes=6
80               nopeexp=21
81            endif
82!
83            do i=1,nope
84               konl(i)=kon(indexe+i)
85               do j=1,3
86                  xl(j,i)=co(j,konl(i))
87               enddo
88            enddo
89!
90            do i=1,nopes
91               if(lakon(nelem)(4:5).eq.'20') then
92                  ibot=nbot20(i)
93                  imid=nmid20(i)
94                  itop=ntop20(i)
95                  xi=xi20(i)
96                  et=et20(i)
97               else
98                  ibot=nbot15(i)
99                  imid=nmid15(i)
100                  itop=ntop15(i)
101                  xi=xi15(i)
102                  et=et15(i)
103               endif
104!
105               nodebot=konl(ibot)
106               nodetop=konl(itop)
107               dd=sqrt((co(1,nodebot)-co(1,nodetop))**2+
108     &                 (co(2,nodebot)-co(2,nodetop))**2+
109     &                 (co(3,nodebot)-co(3,nodetop))**2)
110               do j=0,nlayer-1
111!
112!                 bottom node
113!
114                  node=kon(indexe+nopeexp+j*nope+ibot)
115                  dt=sqrt((co(1,nodebot)-co(1,node))**2+
116     &                 (co(2,nodebot)-co(2,node))**2+
117     &                 (co(3,nodebot)-co(3,node))**2)
118                  ze=2.d0*dt/dd-1.d0
119c                  write(*,*) 'map3dtolayer',node,xi,et,ze
120!
121!                 determining the value of the shape functions
122!
123                  if(lakon(nelem)(4:5).eq.'20') then
124                     call shape20h(xi,et,ze,xl,xsj,shp,iflag)
125                  elseif(lakon(nelem)(4:5).eq.'15') then
126                     call shape15w(xi,et,ze,xl,xsj,shp,iflag)
127                  endif
128!
129                  do k=1,nfield
130                     yn(k,node)=0.d0
131                     do l=1,nope
132                        yn(k,node)=yn(k,node)+
133     &                        shp(4,l)*yn(k,konl(l))
134c                        write(*,*) 'xi',xi,et,ze,node
135c                        write(*,*) l,shp(4,l),yn(k,konl(l))
136                     enddo
137                  enddo
138!
139!                 top node
140!
141                  node=kon(indexe+nopeexp+j*nope+itop)
142                  dt=sqrt((co(1,nodebot)-co(1,node))**2+
143     &                 (co(2,nodebot)-co(2,node))**2+
144     &                 (co(3,nodebot)-co(3,node))**2)
145                  ze=2.d0*dt/dd-1.d0
146c                  write(*,*) 'map3dtolayer',node,xi,et,ze
147!
148!                 determining the value of the shape functions
149!
150                  if(lakon(nelem)(4:5).eq.'20') then
151                     call shape20h(xi,et,ze,xl,xsj,shp,iflag)
152                  elseif(lakon(nelem)(4:5).eq.'15') then
153                     call shape15w(xi,et,ze,xl,xsj,shp,iflag)
154                  endif
155!
156                  do k=1,nfield
157                     yn(k,node)=0.d0
158                     do l=1,nope
159                        yn(k,node)=yn(k,node)+
160     &                        shp(4,l)*yn(k,konl(l))
161                     enddo
162                  enddo
163!
164!                 middle node, if any
165!
166                  if(i.gt.nopes/2) cycle
167!
168                  node=kon(indexe+nopeexp+j*nope+imid)
169                  dt=sqrt((co(1,nodebot)-co(1,node))**2+
170     &                 (co(2,nodebot)-co(2,node))**2+
171     &                 (co(3,nodebot)-co(3,node))**2)
172                  ze=2.d0*dt/dd-1.d0
173c                  write(*,*) 'map3dtolayer',node,xi,et,ze
174!
175!                 determining the value of the shape functions
176!
177                  if(lakon(nelem)(4:5).eq.'20') then
178                     call shape20h(xi,et,ze,xl,xsj,shp,iflag)
179                  elseif(lakon(nelem)(4:5).eq.'15') then
180                     call shape15w(xi,et,ze,xl,xsj,shp,iflag)
181                  endif
182!
183                  do k=1,nfield
184                     yn(k,node)=0.d0
185                     do l=1,nope
186                        yn(k,node)=yn(k,node)+
187     &                        shp(4,l)*yn(k,konl(l))
188                     enddo
189                  enddo
190               enddo
191            enddo
192         endif
193      enddo
194!
195      return
196      end
197