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 genmidnodes(nktet_,ipoed,iedgmid,iexternedg,
20     &     iedgext,cotet,nktet,iedg,jfix,ipoeled,ieled,kontet,
21     &     iedtet,iwrite)
22!
23!     generates midnodes in the middle between the neighboring vertex
24!     nodes
25!
26      implicit none
27!
28      integer nktet_,ipoed(*),index,i,k,iedgmid(*),iexternedg(*),
29     &     iedgext(3,*),node1,node2,iedg(3,*),nktet,jfix(*),ichange,
30     &     ielem,iflag,indexe,iwrite,m,ipoeled(*),ieled(2,*),
31     &     kontet(4,*),iedtet(6,*)
32!
33      real*8 cotet(3,*),xi,et,ze,shp(4,10),xsj,xl(3,10)
34!
35      include "gauss.f"
36!
37      iflag=2
38!
39!     generate middle nodes
40!
41      do i=1,nktet_
42        index=ipoed(i)
43!
44        do
45          if(index.eq.0) exit
46!
47          node1=iedg(1,index)
48          node2=iedg(2,index)
49!
50          if((jfix(node1).eq.1).and.(jfix(node2).eq.1).and.
51     &         (iexternedg(index).gt.0)) then
52            iedgmid(index)=iedgext(2,iexternedg(index))
53          else
54            nktet=nktet+1
55            iedgmid(index)=nktet
56            do k=1,3
57              cotet(k,nktet)=
58     &             (cotet(k,node1)+cotet(k,node2))/2.d0
59            enddo
60          endif
61!
62          index=iedg(3,index)
63        enddo
64      enddo
65!
66!     check quality of elements (may be bad due to the midnodes on
67!     the fixed edges; these are not necessarily in the middle between
68!     the vertex nodes)
69!
70      do
71        ichange=0
72        loop: do i=1,nktet_
73          index=ipoed(i)
74!
75          do
76            if(index.eq.0) exit
77!
78            node1=iedg(1,index)
79            node2=iedg(2,index)
80!
81            if((jfix(node1).eq.1).and.(jfix(node2).eq.1).and.
82     &           (iexternedg(index).gt.0)) then
83!
84!     check whether all elements in the shell of the edge
85!     have positive Jacobians
86!
87              indexe=ipoeled(index)
88              do
89                if(indexe.eq.0) exit
90                ielem=ieled(1,indexe)
91!
92!     vertex nodes
93!
94                do k=1,4
95                  do m=1,3
96                    xl(m,k)=cotet(m,kontet(k,ielem))
97                  enddo
98                enddo
99!
100!     middle nodes
101!
102                do k=1,6
103                  do m=1,3
104                    xl(m,k+4)=cotet(m,iedgmid(iedtet(k,ielem)))
105                  enddo
106                enddo
107!
108                do k=1,4
109                  xi=gauss3d5(1,k)
110                  et=gauss3d5(2,k)
111                  ze=gauss3d5(3,k)
112                  call shape10tet(xi,et,ze,xl,xsj,shp,iflag)
113                  if(xsj.le.0.d0) then
114                    do m=1,3
115                      cotet(m,iedgmid(index))=
116     &                     (cotet(m,node1)+cotet(m,node2))/2.d0
117                    enddo
118                    write(*,*) '*WARNING in genmidnodes: '
119                    write(*,*) '         fixed midnode ',iedgmid(index)
120                    write(*,*)
121     &                   '         had to be moved into the middle'
122                    write(*,*)
123     &                   '         of its neighboring vertex nodes'
124                    write(*,*) '         to keep the adjacent'
125                    write(*,*) '         elements regular'
126                    write(*,*)
127                    write(40,*) iedgmid(index)
128                    iwrite=1
129                    ichange=1
130                    cycle loop
131                  endif
132                enddo
133                indexe=ieled(2,indexe)
134              enddo
135            endif
136!
137            index=iedg(3,index)
138          enddo
139        enddo loop
140        if(ichange.eq.0) exit
141      enddo
142!
143      return
144      end
145