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 projectmidnodes(nktet_,ipoed,iedgmid,iexternedg, 20 & iedgext,cotet,nktet,iedg,iexternfa,ifacext,itreated, 21 & ilist,isharp,ipofa,ifac,iedgextfa,ifacexted,jfix,co,idimsh, 22 & ipoeled,ieled,kontet,c1,jflag,iedtet,ibadnodes,nbadnodes, 23 & iwrite) 24! 25! 1. projects all midnodes on external edges of the unrefined mesh 26! on the parent external edges if they are sharp 27! 28! 2. projects all other midnodes belonging to external faces on 29! neighboring external faces of the unrefined mesh 30! 31 implicit none 32! 33 integer nktet_,ipoed(*),index,i,j,k,iedgmid(*),iexternedg(*), 34 & imasted,iedgext(3,*),nterms,node1,node2,iedg(3,*), 35 & nktet,imastfa,ifacext(6,*),iexternfa(*),itreated(*),ii,node, 36 & isharp(*),ilist(*),nlist,ifac(4,*),iedgextfa(2,*), 37 & ifacexted(3,*),id,indexe,isol,ipofa(*),jfix(*),idimsh(*), 38 & ipoeled(*),ieled(2,*),kontet(4,*),jflag,iedtet(6,*), 39 & ibadnodes(*),nbadnodes,iwrite 40! 41 real*8 pneigh(3,9),cotet(3,*),pnode(3),ratio(9),dist,xi,et, 42 & pnodeproj(3),co(3,*),c1 43! 44 nbadnodes=0 45! 46! set all nodes on "non-treated" 47! 48 do i=1,nktet 49 itreated(i)=0 50 enddo 51! 52! loop over all edges: projection on sharp edges from the unrefined mesh 53! 54 loop1: do i=1,nktet_ 55 index=ipoed(i) 56! 57 do 58 if(index.eq.0) cycle loop1 59! 60 if(iexternedg(index).gt.0) then 61! 62! recovering master edge 63! 64 imasted=iexternedg(index) 65! 66! if the master edge is not sharp: loop 67! 68 if(isharp(imasted).ne.1) then 69 index=iedg(3,index) 70 cycle 71 endif 72! 73! edge is a subset of an external sharp edge 74! of the unrefined mesh 75! 76! end nodes belonging to the edge 77! 78 node1=iedg(1,index) 79 node2=iedg(2,index) 80! 81 if(iedgext(2,imasted).ne.0) then 82 nterms=3 83 do j=1,3 84 do k=1,3 85 pneigh(k,j)=cotet(k,iedgext(j,imasted)) 86 enddo 87 enddo 88 else 89 nterms=2 90 endif 91! 92! edge is a subset of an external sharp edge 93! of the unrefined mesh 94! 95 if((jfix(node1).ne.1).or.(jfix(node2).ne.1)) then 96 node=iedgmid(index) 97! 98 if(nterms.eq.3) then 99 do k=1,3 100 pnode(k)=cotet(k,node) 101 enddo 102! 103! projection for quadratic master edges 104! 105 call attach_1d(pneigh,pnode,nterms,ratio,dist,xi) 106 call checkjac(cotet,node,pnode,kontet,c1,jflag, 107 & iedtet,iedgmid,ipoeled,ieled,index,ibadnodes, 108 & nbadnodes,iwrite) 109! 110 itreated(node)=1 111 endif 112 endif 113 endif 114! 115 index=iedg(3,index) 116 enddo 117 enddo loop1 118! 119! projection of the external nodes not treated yet onto the 120! faces of the unrefined mesh 121! 122 loop2: do i=1,nktet_ 123! 124! loop over all faces of the refined mesh 125! 126 index=ipofa(i) 127 do 128 if(index.eq.0) cycle loop2 129! 130! if no external face: loop 131! 132 if(iexternfa(index).le.0) then 133 index=ifac(4,index) 134 cycle 135 endif 136! 137! external face; treat the midnodes of the face 138! 139 do ii=1,3 140 if(ii.eq.1) then 141 node1=ifac(1,index) 142 node2=ifac(2,index) 143 elseif(ii.eq.2) then 144 node1=ifac(2,index) 145 node2=ifac(3,index) 146 else 147 node1=ifac(1,index) 148 node2=ifac(3,index) 149 endif 150! 151! determine the number of the edge 152! 153 if(node1.gt.node2) then 154 node=node2 155 node2=node1 156 node1=node 157 endif 158 indexe=ipoed(node1) 159 do 160 if(iedg(2,indexe).eq.node2) exit 161 indexe=iedg(3,indexe) 162 enddo 163! 164! check whether a middle node was treated 165! if so, nothing to do -> check the next edge 166! 167 if(itreated(iedgmid(indexe)).ne.0) cycle 168! 169! project the middle node 170! 171 if((jfix(node1).eq.1).and.(jfix(node2).eq.1)) then 172 cycle 173 else 174 node=iedgmid(indexe) 175 do k=1,3 176 pnode(k)=cotet(k,node) 177 enddo 178 endif 179! 180! parent face 181! 182 imastfa=iexternfa(index) 183 ilist(1)=imastfa 184 nlist=1 185! 186! start the loop looking for the correct face; 187! starting with the parent face 188! 189 do 190 if(ifacext(4,imastfa).ne.0) then 191 nterms=6 192 do j=1,6 193 do k=1,3 194 pneigh(k,j)=co(k,ifacext(j,imastfa)) 195 enddo 196 enddo 197 else 198 nterms=3 199 do j=1,3 200 do k=1,3 201 pneigh(k,j)=co(k,ifacext(j,imastfa)) 202 enddo 203 enddo 204 endif 205! 206 do k=1,3 207 pnodeproj(k)=pnode(k) 208 enddo 209 call attach_2d(pneigh,pnodeproj,nterms,ratio,dist, 210 & xi,et) 211! 212! check whether this face is the correct one; 213! if dabs(xi)=1 or dabs(et)=1 or xi+et=0 this 214! may not be the case 215! 216! the solution is found (isol=1) unless proved 217! otherwise 218! 219 isol=1 220! 221 if((dabs(et).lt.1.d-10).and. 222 & (dabs(xi).gt.1.d-10)) then 223! 224! take neighboring face across edge 1-2 unless sharp 225! 226 imasted=ifacexted(1,imastfa) 227 if(isharp(imasted).eq.0) then 228 if(iedgextfa(1,imasted).eq.imastfa) then 229 imastfa=iedgextfa(2,imasted) 230 else 231 imastfa=iedgextfa(1,imasted) 232 endif 233 isol=0 234 endif 235! 236 elseif((dabs(xi+et-1.d0).lt.1.d-10).and. 237 & (dabs(et).gt.1.d-10)) then 238! 239! take neighboring face across edge 2-3 unless sharp 240! 241 imasted=ifacexted(2,imastfa) 242 if(isharp(imasted).eq.0) then 243 if(iedgextfa(1,imasted).eq.imastfa) then 244 imastfa=iedgextfa(2,imasted) 245 else 246 imastfa=iedgextfa(1,imasted) 247 endif 248 isol=0 249 endif 250! 251 elseif((dabs(xi).lt.1.d-10).and. 252 & (dabs(xi+et-1.d0).gt.1.d-10)) then 253! 254! take neighboring face across edge 3-1 unless sharp 255! 256 imasted=ifacexted(3,imastfa) 257 if(isharp(imasted).eq.0) then 258 if(iedgextfa(1,imasted).eq.imastfa) then 259 imastfa=iedgextfa(2,imasted) 260 else 261 imastfa=iedgextfa(1,imasted) 262 endif 263 isol=0 264 endif 265 endif 266! 267! if solution is found: copy projected coordinates 268! else continue with a neighbor 269! 270 if(isol.eq.1) then 271 call checkjac(cotet,node,pnodeproj,kontet,c1,jflag, 272 & iedtet,iedgmid,ipoeled,ieled,indexe,ibadnodes, 273 & nbadnodes,iwrite) 274 itreated(node)=2 275 exit 276 else 277! 278! update list; exit if an element in the list is 279! revisited 280! 281 call nident(ilist,imastfa,nlist,id) 282 if(id.gt.0) then 283 if(ilist(id).eq.imastfa) then 284 if(nlist.eq.2) then 285 isol=1 286 call checkjac(cotet,node,pnodeproj,kontet,c1, 287 & jflag,iedtet,iedgmid,ipoeled,ieled,indexe, 288 & ibadnodes,nbadnodes,iwrite) 289 itreated(node)=2 290 exit 291 endif 292! 293 write(*,*) '*WARNING in projectmidnodes:' 294 write(*,*) ' face in list revisited' 295 write(*,*) ' no projection for node', 296 & node 297 write(*,*) 298 exit 299 endif 300 endif 301 nlist=nlist+1 302 do k=nlist,id+2,-1 303 ilist(k)=ilist(k-1) 304 enddo 305 ilist(id+1)=imastfa 306! 307 endif 308 enddo 309! 310 enddo 311! 312! treat the next face 313! 314 index=ifac(4,index) 315 enddo 316 enddo loop2 317! 318 return 319 end 320