1!
2!     CalculiX - A 3-dimensional finite element program
3!              Copyright (C) 1998 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 checkjac(cotet,node,pnew,kontet,c1,jflag,
20     &     iedtet,iedgmid,ipoeled,ieled,iedge,ibadnodes,nbadnodes,
21     &     iwrite)
22!
23!     check the Jacobian in all integration points of all elements
24!     belonging to midnode "node" if the coordinates of "node" are
25!     changed from cotet(1..3,node) to pnew(1..3); if the Jacobian
26!     is strictly positive, the coordinates are changed, if not, they
27!     are left unchanged
28!     "node" belongs to edge "iedge"
29!
30      implicit none
31!
32      integer node,kontet(4,*),i,ielem,indexe,iflag,jflag,iwrite,
33     &     j,n,iedtet(6,*),iedgmid(*),ipoeled(*),ieled(2,*),
34     &     iedge,k,m,ibadnodes(*),nbadnodes,id,listed,neighedge
35!
36      real*8 cotet(3,*),pnew(3),pold(3),xsj,c1,c2,xi,et,ze,xl(3,10),
37     &     shp(4,10)
38!
39      include "gauss.f"
40!
41      iflag=2
42!
43!     backup of the old coordinates
44!
45      do i=1,3
46        pold(i)=cotet(i,node)
47      enddo
48!
49!     storing the new coordinates in cotet
50!
51      c2=c1
52      do i=1,3
53        cotet(i,node)=c2*pnew(i)+(1.d0-c2)*pold(i)
54      enddo
55!
56!     trial loops; if bad elements the projection is reduced and the
57!     elements are rechecked
58!
59      n=3
60      loop1: do j=1,n
61      indexe=ipoeled(iedge)
62!
63!       loop over all elements belonging to edge "iedge"
64!
65        loop2: do
66          if(indexe.eq.0) exit loop1
67          ielem=ieled(1,indexe)
68!
69!         vertex nodes
70!
71          do k=1,4
72            do m=1,3
73              xl(m,k)=cotet(m,kontet(k,ielem))
74            enddo
75          enddo
76!
77!         middle nodes
78!
79          do k=1,6
80            do m=1,3
81              xl(m,k+4)=cotet(m,iedgmid(iedtet(k,ielem)))
82            enddo
83          enddo
84!
85          do k=1,4
86            xi=gauss3d5(1,k)
87            et=gauss3d5(2,k)
88            ze=gauss3d5(3,k)
89            call shape10tet(xi,et,ze,xl,xsj,shp,iflag)
90            if(xsj.le.0.d0) then
91c              write(*,*) 'checkjac bad element ',ielem,' node ',node
92              exit loop2
93            endif
94          enddo
95          indexe=ieled(2,indexe)
96        enddo loop2
97!
98!       last loop: revert to non-projected coordinates
99!
100        if(j.eq.n) then
101          do i=1,3
102            cotet(i,node)=pold(i)
103          enddo
104          exit
105        endif
106!
107!       not last loop: reduce the projection and retry
108!
109        c2=c2/2.d0
110        do i=1,3
111          cotet(i,node)=c2*pnew(i)+(1.d0-c2)*pold(i)
112        enddo
113      enddo loop1
114!
115!     if node was not completely projected: list all edges of
116!     neighboring elements (corresponds 1-to-1 to midnodes)
117!     -> application of fminsirefine for all neighboring edges/midnodes
118!        which are not on the free surface
119!
120      if(j.ne.1) then
121        indexe=ipoeled(iedge)
122        do
123          if(indexe.eq.0) exit
124!
125!         adjacent element
126!
127          ielem=ieled(1,indexe)
128          do k=1,6
129            neighedge=iedtet(k,ielem)
130!
131!           edge of adjacent element
132!
133            listed=0
134            call nident(ibadnodes,neighedge,nbadnodes,id)
135            if(id.gt.0) then
136              if(ibadnodes(id).eq.neighedge) then
137                listed=1
138              endif
139            endif
140            if(listed.eq.0) then
141              nbadnodes=nbadnodes+1
142              do m=nbadnodes,id+2,-1
143                ibadnodes(m)=ibadnodes(m-1)
144              enddo
145              ibadnodes(id+1)=neighedge
146            endif
147          enddo
148          indexe=ieled(2,indexe)
149        enddo
150      endif
151!
152!     if jflag=1 and j>1: last chance for a complete projection
153!     did not work
154!
155      if((jflag.eq.1).and.(j.ne.1)) then
156        write(*,*) '*WARNING in checkjac: projection of midnode ',node
157        write(*,*) '         had to be reduced to keep the adjacent'
158        write(*,*) '         elements regular'
159        write(*,*)
160        write(40,*) node
161        iwrite=1
162      endif
163!
164      return
165      end
166