!
!     CalculiX - A 3-dimensional finite element program
!              Copyright (C) 1998-2007 Guido Dhondt
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation(version 2);
!     
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of 
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with this program; if not, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
      subroutine extrapolate(yi,yn,ipkon,inum,kon,lakon,nfield,nk,
     &  ne,mint_,ndim,orab,ielorien,co,iorienglob,cflag,nelemload,
     &  nload,nodeboun,nboun)
!
!     extrapolates field values at the integration points to the 
!     nodes
!
!     the number of internal state variables is limited to 999
!     (cfr. array field)
!
      implicit none
!
      character*1 cflag
      character*8 lakon(*),lakonl
!
      integer ipkon(*),inum(*),kon(*),mint_,ne,indexe,nope,
     &  nonei20(3,12),nfield,nonei10(3,6),nk,i,j,k,l,ndim,
     &  nonei15(3,9),iorienglob,iorien,ielorien(*),konl,
     &  mint3d,m,iflag,nelemload(2,*),nload,node,nboun,
     &  nodeboun(*)
!
      real*8 yi(ndim,mint_,*),yn(nfield,*),field(999,20),a8(8,8),
     &  a4(4,4),a27(20,27),a9(6,9),a2(6,2),orab(7,*),co(3,*),
     &  coords(3,27),xi,et,ze,xl(3,20),xsj,shp(4,20),weight,
     &  yiloc(6,27),a(3,3),b(3,3),c(3,3)
!
      real*8 gauss2d1(2,1),gauss2d2(2,4),gauss2d3(2,9),gauss2d4(2,1),
     &  gauss2d5(2,3),gauss3d1(3,1),gauss3d2(3,8),gauss3d3(3,27),
     &  gauss3d4(3,1),gauss3d5(3,4),gauss3d6(3,15),gauss3d7(3,2),
     &  gauss3d8(3,9),gauss3d9(3,18),weight2d1(1),weight2d2(4),
     &  weight2d3(9),weight2d4(1),weight2d5(3),weight3d1(1),
     &  weight3d2(8),weight3d3(27),weight3d4(1),weight3d5(4),
     &  weight3d6(15),weight3d7(2),weight3d8(9),weight3d9(18)
!
      include "gauss.f"
!
      data nonei10 /5,1,2,6,2,3,7,3,1,8,1,4,9,2,4,10,3,4/
!
      data nonei15 /7,1,2,8,2,3,9,3,1,10,4,5,11,5,6,12,6,4,
     &  13,1,4,14,2,5,15,3,6/
!
      data nonei20 /9,1,2,10,2,3,11,3,4,12,4,1,
     &  13,5,6,14,6,7,15,7,8,16,8,5,
     &  17,1,5,18,2,6,19,3,7,20,4,8/
!
      data a2 /  1.1455,-0.1455,1.1455,-0.1455,1.1455,-0.1455,
     &           -0.1455,1.1455,-0.1455,1.1455,-0.1455,1.1455/
      data a4 /  1.92705, -0.30902, -0.30902, -0.30902,
     &          -0.30902,  1.92705, -0.30902, -0.30902,
     &          -0.30902, -0.30902,  1.92705, -0.30902,
     &          -0.30902, -0.30902, -0.30902,  1.92705/
      data a9 / 1.63138,-0.32628,-0.32628,-0.52027, 0.10405, 0.10405,
     &         -0.32628, 1.63138,-0.32628, 0.10405,-0.52027, 0.10405,
     &         -0.32628,-0.32628, 1.63138, 0.10405, 0.10405,-0.52027,
     &          0.55556,-0.11111,-0.11111, 0.55556,-0.11111,-0.11111,
     &         -0.11111, 0.55556,-0.11111,-0.11111, 0.55556,-0.11111,
     &         -0.11111,-0.11111, 0.55556,-0.11111,-0.11111, 0.55556,
     &         -0.52027, 0.10405, 0.10405, 1.63138,-0.32628,-0.32628,
     &          0.10405,-0.52027, 0.10405,-0.32628, 1.63138,-0.32628,
     &          0.10405, 0.10405,-0.52027,-0.32628,-0.32628, 1.63138/
      data a8 /2.549,-.683,.183,-.683,-.683,.183,
     &        -.04904,.183,-.683,2.549,-.683,.183,
     &        .183,-.683,.183,-.04904,-.683,.183,
     &        -.683,2.549,.183,-.04904,.183,-.683,
     &        .183,-.683,2.549,-.683,-.04904,.183,
     &        -.683,.183,-.683,.183,-.04904,.183,
     &        2.549,-.683,.183,-.683,.183,-.683,
     &        .183,-.04904,-.683,2.549,-.683,.183,
     &        .183,-.04904,.183,-.683,-.683,.183,
     &        -.683,2.549,-.04904,.183,-.683,.183,
     &        .183,-.683,2.549,-.683/      
      data a27 /
     &  2.37499,-0.12559,-0.16145,-0.12559,-0.12559,-0.16145, 0.11575,
     & -0.16145, 0.32628, 0.11111, 0.11111, 0.32628, 0.11111,-0.10405,
     & -0.10405, 0.11111, 0.32628, 0.11111,-0.10405, 0.11111,-0.31246,
     & -0.31246, 0.31481, 0.31481, 0.31481, 0.31481,-0.16902,-0.16902,
     &  1.28439,-0.27072,-0.19444,-0.27072,-0.19444, 0.15961,-0.00661,
     &  0.15961,-0.27072,-0.27072, 0.15961, 0.15961,-0.12559, 2.37499,
     & -0.12559,-0.16145,-0.16145,-0.12559,-0.16145, 0.11575, 0.32628,
     &  0.32628, 0.11111, 0.11111, 0.11111, 0.11111,-0.10405,-0.10405,
     &  0.11111, 0.32628, 0.11111,-0.10405,-0.31246, 0.31481, 0.31481,
     & -0.31246, 0.31481,-0.16902,-0.16902, 0.31481,-0.27072,-0.19444,
     & -0.27072, 1.28439, 0.15961,-0.00661, 0.15961,-0.19444,-0.27072,
     &  0.15961, 0.15961,-0.27072,-0.48824,-0.48824,-0.48824,-0.48824,
     &  0.22898, 0.22898, 0.22898, 0.22898, 0.05556, 0.05556, 0.05556,
     &  0.05556, 0.05556, 0.05556, 0.05556, 0.05556,-0.22222,-0.22222,
     & -0.22222,-0.22222, 0.31481,-0.31246,-0.31246, 0.31481,-0.16902,
     &  0.31481, 0.31481,-0.16902,-0.27072, 1.28439,-0.27072,-0.19444,
     &  0.15961,-0.19444, 0.15961,-0.00661, 0.15961,-0.27072,-0.27072,
     &  0.15961,-0.12559,-0.16145,-0.12559, 2.37499,-0.16145, 0.11575,
     & -0.16145,-0.12559, 0.11111, 0.11111, 0.32628, 0.32628,-0.10405,
     & -0.10405, 0.11111, 0.11111, 0.11111,-0.10405, 0.11111, 0.32628,
     &  0.31481, 0.31481,-0.31246,-0.31246,-0.16902,-0.16902, 0.31481,
     &  0.31481,-0.19444,-0.27072, 1.28439,-0.27072,-0.00661, 0.15961,
     & -0.19444, 0.15961, 0.15961, 0.15961,-0.27072,-0.27072,-0.16145,
     & -0.12559, 2.37499,-0.12559, 0.11575,-0.16145,-0.12559,-0.16145,
     &  0.11111, 0.32628, 0.32628, 0.11111,-0.10405, 0.11111, 0.11111,
     & -0.10405,-0.10405, 0.11111, 0.32628, 0.11111,-0.31246, 0.31481,
     & -0.16902, 0.31481,-0.31246, 0.31481,-0.16902, 0.31481,-0.27072,
     &  0.15961, 0.15961,-0.27072,-0.27072, 0.15961, 0.15961,-0.27072,
     &  1.28439,-0.19444,-0.00661,-0.19444,-0.48824,-0.48824, 0.22898,
     &  0.22898,-0.48824,-0.48824, 0.22898, 0.22898, 0.05556,-0.22222,
     &  0.05556,-0.22222, 0.05556,-0.22222, 0.05556,-0.22222, 0.05556,
     &  0.05556, 0.05556, 0.05556, 0.31481,-0.31246, 0.31481,-0.16902,
     &  0.31481,-0.31246, 0.31481,-0.16902,-0.27072,-0.27072, 0.15961,
     &  0.15961,-0.27072,-0.27072, 0.15961, 0.15961,-0.19444, 1.28439,
     & -0.19444,-0.00661,-0.48824, 0.22898, 0.22898,-0.48824,-0.48824,
     &  0.22898, 0.22898,-0.48824,-0.22222, 0.05556,-0.22222, 0.05556,
     & -0.22222, 0.05556,-0.22222, 0.05556, 0.05556, 0.05556, 0.05556,
     &  0.05556,-0.29630,-0.29630,-0.29630,-0.29630,-0.29630,-0.29630,
     & -0.29630,-0.29630,-0.11111,-0.11111,-0.11111,-0.11111,-0.11111,
     & -0.11111,-0.11111,-0.11111,-0.11111,-0.11111,-0.11111,-0.11111,
     &  0.22898,-0.48824,-0.48824, 0.22898, 0.22898,-0.48824,-0.48824,
     &  0.22898,-0.22222, 0.05556,-0.22222, 0.05556,-0.22222, 0.05556,
     & -0.22222, 0.05556, 0.05556, 0.05556, 0.05556, 0.05556, 0.31481,
     & -0.16902, 0.31481,-0.31246, 0.31481,-0.16902, 0.31481,-0.31246,
     &  0.15961, 0.15961,-0.27072,-0.27072, 0.15961, 0.15961,-0.27072,
     & -0.27072,-0.19444,-0.00661,-0.19444, 1.28439, 0.22898, 0.22898,
     & -0.48824,-0.48824, 0.22898, 0.22898,-0.48824,-0.48824, 0.05556,
     & -0.22222, 0.05556,-0.22222, 0.05556,-0.22222, 0.05556,-0.22222,
     &  0.05556, 0.05556, 0.05556, 0.05556,-0.16902, 0.31481,-0.31246,
     &  0.31481,-0.16902, 0.31481,-0.31246, 0.31481, 0.15961,-0.27072,
     & -0.27072, 0.15961, 0.15961,-0.27072,-0.27072, 0.15961,-0.00661,
     & -0.19444, 1.28439,-0.19444,-0.12559,-0.16145, 0.11575,-0.16145,
     &  2.37499,-0.12559,-0.16145,-0.12559, 0.11111,-0.10405,-0.10405,
     &  0.11111, 0.32628, 0.11111, 0.11111, 0.32628, 0.32628, 0.11111,
     & -0.10405, 0.11111, 0.31481, 0.31481,-0.16902,-0.16902,-0.31246,
     & -0.31246, 0.31481, 0.31481,-0.19444, 0.15961,-0.00661, 0.15961,
     &  1.28439,-0.27072,-0.19444,-0.27072,-0.27072,-0.27072, 0.15961,
     &  0.15961,-0.16145,-0.12559,-0.16145, 0.11575,-0.12559, 2.37499,
     & -0.12559,-0.16145, 0.11111, 0.11111,-0.10405,-0.10405, 0.32628,
     &  0.32628, 0.11111, 0.11111, 0.11111, 0.32628, 0.11111,-0.10405,
     &  0.31481,-0.16902,-0.16902, 0.31481,-0.31246, 0.31481, 0.31481,
     & -0.31246, 0.15961,-0.00661, 0.15961,-0.19444,-0.27072,-0.19444,
     & -0.27072, 1.28439,-0.27072, 0.15961, 0.15961,-0.27072, 0.22898,
     &  0.22898, 0.22898, 0.22898,-0.48824,-0.48824,-0.48824,-0.48824,
     &  0.05556, 0.05556, 0.05556, 0.05556, 0.05556, 0.05556, 0.05556,
     &  0.05556,-0.22222,-0.22222,-0.22222,-0.22222,-0.16902, 0.31481,
     &  0.31481,-0.16902, 0.31481,-0.31246,-0.31246, 0.31481, 0.15961,
     & -0.19444, 0.15961,-0.00661,-0.27072, 1.28439,-0.27072,-0.19444,
     &  0.15961,-0.27072,-0.27072, 0.15961,-0.16145, 0.11575,-0.16145,
     & -0.12559,-0.12559,-0.16145,-0.12559, 2.37499,-0.10405,-0.10405,
     &  0.11111, 0.11111, 0.11111, 0.11111, 0.32628, 0.32628, 0.11111,
     & -0.10405, 0.11111, 0.32628,-0.16902,-0.16902, 0.31481, 0.31481,
     &  0.31481, 0.31481,-0.31246,-0.31246,-0.00661, 0.15961,-0.19444,
     &  0.15961,-0.19444,-0.27072, 1.28439,-0.27072, 0.15961, 0.15961,
     & -0.27072,-0.27072, 0.11575,-0.16145,-0.12559,-0.16145,-0.16145,
     & -0.12559, 2.37499,-0.12559,-0.10405, 0.11111, 0.11111,-0.10405,
     &  0.11111, 0.32628, 0.32628, 0.11111,-0.10405, 0.11111, 0.32628,
     &  0.11111/
      data iflag /1/
!
      do i=1,nk
         inum(i)=0
      enddo
!
      do i=1,nk
         do j=1,nfield
            yn(j,i)=0.d0
         enddo
      enddo
!
      do i=1,ne
!
         if(ipkon(i).lt.0) cycle
         indexe=ipkon(i)
         lakonl=lakon(i)
!
         if(lakonl(4:4).eq.'2') then
            nope=20
         elseif(lakonl(4:4).eq.'8') then
            nope=8
         elseif(lakonl(4:5).eq.'10') then
            nope=10
         elseif(lakonl(4:4).eq.'4') then
            nope=4
         elseif(lakonl(4:5).eq.'15') then
            nope=15
         elseif(lakonl(4:4).eq.'6') then
            nope=6
         elseif(lakonl(1:1).eq.'D') then
            if(kon(indexe+1).ne.0) 
     &           inum(kon(indexe+1))=inum(kon(indexe+1))+1
            inum(kon(indexe+2))=inum(kon(indexe+2))+1
            if(kon(indexe+3).ne.0)
     &           inum(kon(indexe+3))=inum(kon(indexe+3))+1
            cycle
         else
            cycle
         endif
!
!     calculation of the integration point coordinates for
!     output in the local system
!
         if((iorienglob.ne.0).and.(ielorien(i).ne.0)) then
!
            iorien=ielorien(i)
!
            if(lakon(i)(4:5).eq.'8R') then
               mint3d=1
            elseif((lakon(i)(4:4).eq.'8').or.
     &              (lakon(i)(4:6).eq.'20R')) then
               mint3d=8
            elseif(lakon(i)(4:4).eq.'2') then
               mint3d=27
            elseif(lakon(i)(4:5).eq.'10') then
               mint3d=4
            elseif(lakon(i)(4:4).eq.'4') then
               mint3d=1
            elseif(lakon(i)(4:5).eq.'15') then
               mint3d=9
            elseif(lakon(i)(4:4).eq.'6') then
               mint3d=2
            endif
!
            do j=1,nope
               konl=kon(indexe+j)
               do k=1,3
                  xl(k,j)=co(k,konl)
               enddo
            enddo
!
            do j=1,mint3d
               if(lakon(i)(4:5).eq.'8R') then
                  xi=gauss3d1(1,j)
                  et=gauss3d1(2,j)
                  ze=gauss3d1(3,j)
                  weight=weight3d1(j)
               elseif((lakon(i)(4:4).eq.'8').or.
     &                 (lakon(i)(4:6).eq.'20R'))
     &                 then
                  xi=gauss3d2(1,j)
                  et=gauss3d2(2,j)
                  ze=gauss3d2(3,j)
                  weight=weight3d2(j)
               elseif(lakon(i)(4:4).eq.'2') then
                  xi=gauss3d3(1,j)
                  et=gauss3d3(2,j)
                  ze=gauss3d3(3,j)
                  weight=weight3d3(j)
               elseif(lakon(i)(4:5).eq.'10') then
                  xi=gauss3d5(1,j)
                  et=gauss3d5(2,j)
                  ze=gauss3d5(3,j)
                  weight=weight3d5(j)
               elseif(lakon(i)(4:4).eq.'4') then
                  xi=gauss3d4(1,j)
                  et=gauss3d4(2,j)
                  ze=gauss3d4(3,j)
                  weight=weight3d4(j)
               elseif(lakon(i)(4:5).eq.'15') then
                  xi=gauss3d8(1,j)
                  et=gauss3d8(2,j)
                  ze=gauss3d8(3,j)
                  weight=weight3d8(j)
               elseif(lakon(i)(4:4).eq.'6') then
                  xi=gauss3d7(1,j)
                  et=gauss3d7(2,j)
                  ze=gauss3d7(3,j)
                  weight=weight3d7(j)
               endif
!
               if(nope.eq.20) then
                  call shape20h(xi,et,ze,xl,xsj,shp,iflag)
               elseif(nope.eq.8) then
                  call shape8h(xi,et,ze,xl,xsj,shp,iflag)
               elseif(nope.eq.10) then
                  call shape10tet(xi,et,ze,xl,xsj,shp,iflag)
               elseif(nope.eq.4) then
                  call shape4tet(xi,et,ze,xl,xsj,shp,iflag)
               elseif(nope.eq.15) then
                  call shape15w(xi,et,ze,xl,xsj,shp,iflag)
               else
                  call shape6w(xi,et,ze,xl,xsj,shp,iflag)
               endif
!
               do k=1,3
                  coords(k,j)=0.d0
                  do l=1,nope
                     coords(k,j)=coords(k,j)+xl(k,l)*shp(4,l)
                  enddo
               enddo
!
               if(nfield.eq.3) then
                  call transformatrix(orab(1,iorien),coords(1,j),a)
                  yiloc(1,j)=yi(1,j,i)*a(1,1)+yi(2,j,i)*a(2,1)+
     &                     yi(3,j,i)*a(3,1)
                  yiloc(2,j)=yi(1,j,i)*a(1,2)+yi(2,j,i)*a(2,2)+
     &                     yi(3,j,i)*a(3,2)
                  yiloc(3,j)=yi(1,j,i)*a(1,3)+yi(2,j,i)*a(2,3)+
     &                     yi(3,j,i)*a(3,3)
               elseif(nfield.eq.6) then
                  call transformatrix(orab(1,iorien),coords(1,j),a)
                  b(1,1)=yi(1,j,i)
                  b(2,2)=yi(2,j,i)
                  b(3,3)=yi(3,j,i)
                  b(1,2)=yi(4,j,i)
                  b(1,3)=yi(5,j,i)
                  b(2,3)=yi(6,j,i)
                  b(2,1)=b(1,2)
                  b(3,1)=b(1,3)
                  b(3,2)=b(2,3)
                  do k=1,3
                     do l=1,3
                        c(k,l)=0.d0
                        do m=1,3
                           c(k,l)=c(k,l)+b(k,m)*a(m,l)
                        enddo
                     enddo
                  enddo
                  do k=1,3
                     do l=k,3
                        b(k,l)=0.d0
                        do m=1,3
                           b(k,l)=b(k,l)+a(m,k)*c(m,l)
                        enddo
                     enddo
                  enddo
                  yiloc(1,j)=b(1,1)
                  yiloc(2,j)=b(2,2)
                  yiloc(3,j)=b(3,3)
                  yiloc(4,j)=b(1,2)
                  yiloc(5,j)=b(1,3)
                  yiloc(6,j)=b(2,3)
               endif
            enddo
!
            if((lakonl(4:6).eq.'20R').or.(lakonl(4:5).eq.'8 ')) then
               do j=1,8
                  do k=1,nfield
                     field(k,j)=0.d0
                     do l=1,8
                        field(k,j)=field(k,j)+a8(j,l)*yiloc(k,l)
                     enddo
                  enddo
               enddo
            elseif(lakonl(4:4).eq.'8') then
               do j=1,8
                  do k=1,nfield
                     field(k,j)=yiloc(k,1)
                  enddo
               enddo
            elseif(lakonl(4:5).eq.'10') then
               do j=1,4
                  do k=1,nfield
                     field(k,j)=0.d0
                     do l=1,4
                        field(k,j)=field(k,j)+a4(j,l)*yiloc(k,l)
                     enddo
                  enddo
               enddo
            elseif(lakonl(4:4).eq.'2') then
               do j=1,20
                  do k=1,nfield
                     field(k,j)=0.d0
                     do l=1,27
                        field(k,j)=field(k,j)+a27(j,l)*yiloc(k,l)
                     enddo
                  enddo
               enddo
            elseif(lakonl(4:4).eq.'4') then
               do j=1,4
                  do k=1,nfield
                     field(k,j)=yiloc(k,1)
                  enddo
               enddo
            elseif(lakonl(4:4).eq.'1') then
               do j=1,6
                  do k=1,nfield
                     field(k,j)=0.d0
                     do l=1,9
                        field(k,j)=field(k,j)+a9(j,l)*yiloc(k,l)
                     enddo
                  enddo
               enddo
            else
               do j=1,6
                  do k=1,nfield
                     field(k,j)=0.d0
                     do l=1,2
                        field(k,j)=field(k,j)+a2(j,l)*yiloc(k,l)
                     enddo
                  enddo
               enddo
            endif
         else
!
!        determining the field values in the vertex nodes
!        for C3D20R and C3D8: trilinear extrapolation (= use of the
!                             C3D8 linear interpolation functions)
!        for C3D8R: constant field value in each element
!        for C3D10: use of the C3D4 linear interpolation functions
!        for C3D4: constant field value in each element
!        for C3D15: use of the C3D6 linear interpolation functions
!        for C3D6: use of a linear interpolation function
!
            if((lakonl(4:6).eq.'20R').or.(lakonl(4:5).eq.'8 ')) then
               do j=1,8
                  do k=1,nfield
                     field(k,j)=0.d0
                     do l=1,8
                        field(k,j)=field(k,j)+a8(j,l)*yi(k,l,i)
                     enddo
                  enddo
               enddo
            elseif(lakonl(4:4).eq.'8') then
               do j=1,8
                  do k=1,nfield
                     field(k,j)=yi(k,1,i)
                  enddo
               enddo
            elseif(lakonl(4:5).eq.'10') then
               do j=1,4
                  do k=1,nfield
                     field(k,j)=0.d0
                     do l=1,4
                        field(k,j)=field(k,j)+a4(j,l)*yi(k,l,i)
                     enddo
                  enddo
               enddo
            elseif(lakonl(4:4).eq.'2') then
               do j=1,20
                  do k=1,nfield
                     field(k,j)=0.d0
                     do l=1,27
                        field(k,j)=field(k,j)+a27(j,l)*yi(k,l,i)
                     enddo
                  enddo
               enddo
            elseif(lakonl(4:4).eq.'4') then
               do j=1,4
                  do k=1,nfield
                     field(k,j)=yi(k,1,i)
                  enddo
               enddo
            elseif(lakonl(4:4).eq.'1') then
               do j=1,6
                  do k=1,nfield
                     field(k,j)=0.d0
                     do l=1,9
                        field(k,j)=field(k,j)+a9(j,l)*yi(k,l,i)
                     enddo
                  enddo
               enddo
            else
               do j=1,6
                  do k=1,nfield
                     field(k,j)=0.d0
                     do l=1,2
                        field(k,j)=field(k,j)+a2(j,l)*yi(k,l,i)
                     enddo
                  enddo
               enddo
            endif
         endif
!
!        determining the field values in the midside nodes
!
         if(lakonl(4:6).eq.'20R') then
            do j=9,20
               do k=1,nfield
                  field(k,j)=(field(k,nonei20(2,j-8))+
     &                 field(k,nonei20(3,j-8)))/2.d0
               enddo
            enddo
         elseif(lakonl(4:5).eq.'10') then
            do j=5,10
               do k=1,nfield
                  field(k,j)=(field(k,nonei10(2,j-4))+
     &                 field(k,nonei10(3,j-4)))/2.d0
               enddo
            enddo
         elseif(lakonl(4:5).eq.'15') then
            do j=7,15
               do k=1,nfield
                  field(k,j)=(field(k,nonei15(2,j-6))+
     &                 field(k,nonei15(3,j-6)))/2.d0
               enddo
            enddo
         endif
!
!        transferring the field values into yn
!
         do j=1,nope
            do k=1,nfield
               yn(k,kon(indexe+j))=yn(k,kon(indexe+j))+
     &              field(k,j)
            enddo
            inum(kon(indexe+j))=inum(kon(indexe+j))+1
         enddo
!
      enddo
!
!     taking the mean of nodal contributions coming from different
!     elements having the node in common
!
      do i=1,nk
         if(inum(i).gt.0) then
            do j=1,nfield
               yn(j,i)=yn(j,i)/inum(i)
            enddo
         endif
      enddo
!
!     for 1d and 2d elements only:
!     finding the solution in the original nodes
!
      if((cflag.ne.' ').and.(cflag.ne.'E')) then
         call map3dto1d2d(yn,ipkon,inum,kon,lakon,nfield,nk,ne,cflag,co)
      endif
!
!     printing values for environmental film and radiation nodes
!
      do i=1,nload
         node=nelemload(2,i)
         if(node.gt.0) then
            if(inum(node).gt.0) cycle
            inum(node)=1
         endif
      enddo
!
!     printing values of prescribed boundary conditions
!
      do i=1,nboun
         node=nodeboun(i)
         if(inum(node).gt.0) cycle
         inum(node)=1
      enddo
!
      return
      end