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! 20! Cuts a convex part of the master surface with a slave surface 21! inserts new active edges into iactiveline for current triangle 22! calling sutherland-hodgeman algorithm 23! 24! [in] nopes number of slave nodes for current slave surface 25! [in] slavstraight plane equations for mean slave plane 26! [in] xn mean slave normal 27! [in] xns slave normals in nodes of slave surface 28! [in] xl2s current positions of lsave nodes 29! [in] xl2sp projected positions of slave nodes 30! [in] ipe (i) pointer to ime for node i 31! [in] ime ... cataloging the edges with node i 32! [in,out] iactiveline field storing the active master triange lines for the active line search 33! [in,out] nactiveline number of active lines 34! [in] nelemm current master element 35! [in,out] nintpoint number of generated integration points 36! [in,out] pslavsurf field storing position xil, etal and weight for integration point on slave side 37! [in,out] imastsurf pointer into pmastsurf 38! [in,out] pmastsurf field storing position and weights for integration points on master side 39! [in,out] xl2m coordinates of master nodes for current master face 40! [in,out] nnodelem number of nodes in curretn master element 41! [in,out] xl2m2 resorted coordinates of master nodes for current master face 42! [in] nmp number of master nodes for current master face 43! [in] nodem node numbers for current master face 44! [in,out] gapmints gap evaluated at the integration points 45! [in] issurf current slave surface index 46! [in,out] areaslav current covering of the slave surface (in the reference element) 47! [in] debug debug output parameter 48! 49 subroutine treatmasterface_mortar( 50 & nopes,slavstraight,xn,xns,xl2s,xl2sp, 51 & ipe,ime,iactiveline,nactiveline, 52 & nelemm,nintpoint,pslavsurf, 53 & imastsurf,pmastsurf,xl2m,nnodelem,xl2m2,nmp, 54 & nodem,gapmints,issurf,areaslav,debug,clearance, 55 & cl2s,cl2m,shrink,reltime) 56! 57! Author: Saskia Sitzmann 58! 59 implicit none 60! 61 logical debug,shrink 62! 63 integer nvertex,nopes,ipe(*),ime(4,*),iactiveline(3,*), 64 & nactiveline,ifreeintersec,nmp,i,j,k,nintpoint,imastsurf(*), 65 & nnodelem,ijk,nodem(*),modf,nelemm,k_max,issurf,ii,jj,nipold 66! 67 real*8 pvertex(3,13),slavstraight(36),xn(3), 68 & xilm,etlm,xnl(3),clearance,p(3),dist, 69 & xl2s(3,*),p1(3),p2(3),pslavsurf(3,*), 70 & xil,etl,area,areax,areay,areaz,pmastsurf(2,*), 71 & xl2m(3,8),xl2m2(3,8),al,gapmints(*),err,xns(3,8), 72 & xl2sp(3,*),xl2mp(3,8),cgp(3),spm,spmc, 73 & pm(3),ph(3),reltime,ps(3),xit(3),etat(3),phc(3), 74 & areaslav,cl2s(3,8),cl2m(3,8),psc(3),pmc(3) 75! 76! 77! 78 data ijk /0/ 79 save ijk 80! 81 include "gauss.f" 82! 83 ifreeintersec=0 84 err=1.d-6 85 nvertex=0 86 nipold=nintpoint 87! 88 if(debug) write(20,*) 'TT:slavsurf',issurf,' melem',nelemm 89 if(debug) write(20,*) 'TT:xn',(xn(1:3)) 90! 91! Project master nodes to meanplane, needed for Sutherland-Hodgman 92! 93 do j=1,nmp 94 al=-xn(1)*xl2m2(1,j)-xn(2)* 95 & xl2m2(2,j)-xn(3)* 96 & xl2m2(3,j)-slavstraight(nopes*4+4) 97 do k=1,3 98 xl2mp(k,j)=xl2m2(k,j)+al*xn(k) 99 enddo 100 enddo 101! 102 if(debug) then 103 write(20,*) 'TT: xl2sp' 104 do j=1,nopes 105 write(20,*)(xl2sp(k,j),k=1,3) 106 enddo 107 write(20,*)'TT:xl2mp' 108 do j=1,nmp 109 write(20,*)(xl2mp(k,j),k=1,3) 110 enddo 111 endif 112! 113 111 format(3(e27.20)) 114! 115! call Sutherland-Hodgman Algo 116! 117 call sutherland_hodgman(nopes,xn,xl2sp,xl2mp,nodem, 118 & ipe,ime,iactiveline,nactiveline, 119 & ifreeintersec,nelemm,nmp, 120 & nvertex,pvertex) 121! 122! do we have a degenerated triangle? 123! 124 if(debug)then 125 write(20,*) 'nelemm',nelemm ,'p_new' 126 do k=1,nvertex 127 write(20,111)(pvertex(i,k),i=1,3) 128 enddo 129 endif 130! 131 do k=1,3 132 cgp(k)=0.0 133 enddo 134 if(nvertex.lt.3) return 135! 136 if(nvertex.eq.3)then 137 do k=1,3 138 cgp(k)=pvertex(k,nvertex) 139 enddo 140 nvertex=nvertex-1 141 k_max=1 142 else 143 do i=1,nvertex 144 do k=1,3 145 cgp(k)=cgp(k)+pvertex(k,i)/nvertex 146 enddo 147 enddo 148 k_max=nvertex 149 endif 150! 151 if(debug)then 152 write(20,*) 'TT: nactiveline',nactiveline,'iactiveline' 153 write(20,*) (iactiveline(1,k),k=1,nactiveline) 154 write(20,*) 'cg' ,(cgp(k),k=1,3) 155 write(20,*)'********************************************' 156 endif 157! 158! Project center point back on slave face 159! 160 call attachline(xl2s,cgp,nopes,xit(3),etat(3),xn,p,dist) 161! 162! generating integration points on the slave surface S 163! 164 do k=1,k_max 165! 166! storing the triangulation of the slave surfaces 167! 168 if(debug) then 169 ijk=ijk+1 170 write(40,100) ijk,(cgp(i),i=1,3) 171 ijk=ijk+1 172 write(40,100) ijk,(pvertex(i,modf(nvertex,k)),i=1,3) 173 ijk=ijk+1 174 write(40,100) ijk,(pvertex(i,modf(nvertex,k+1)),i=1,3) 175 write(40,101) ijk-2,ijk-2,ijk-1 176 write(40,101) ijk-1,ijk-1,ijk 177 write(40,101) ijk,ijk,ijk-2 178 endif 179 100 format('PNT ',i10,'P',3(1x,e21.14)) 180 101 format('LINE ',i10,'L',i10,'P ',i10,'P') 181! 182! Project back on slave surface 183! 184 call attachline(xl2s,pvertex(1:3,modf(nvertex,k)), 185 & nopes,xit(1),etat(1),xn,p,dist) 186 call attachline(xl2s,pvertex(1:3,modf(nvertex,k+1)), 187 & nopes,xit(2),etat(2),xn,p,dist) 188 p1(1)=xit(1)-xit(3) 189 p1(2)=etat(1)-etat(3) 190 p1(3)=0 191 p2(1)=xit(2)-xit(3) 192 p2(2)=etat(2)-etat(3) 193 p2(3)=0 194 areax=((p1(2)*p2(3))-(p2(2)*p1(3)))**2 195 areay=(-(p1(1)*p2(3))+(p2(1)*p1(3)))**2 196 areaz=((p1(1)*p2(2))-(p2(1)*p1(2)))**2 197 area=dsqrt(areax+areay+areaz)/2. 198 if(area.lt.1.e-4) cycle 199 if((nopes.eq.4 .or. nopes.eq.8) 200 & .and.areaslav+area-4.0.gt.1.e-3 201 & .and.nactiveline.gt.0)then 202 write(*,*)'TT: face',issurf,'loop in slavintmortar' 203 write(*,*)'area',areaslav,'+',area,'.gt.4!',nactiveline 204 nactiveline=0 205 return 206 endif 207 if((nopes.eq.3.or.nopes.eq.6) 208 & .and.areaslav+area-0.5.gt.1.e-4 209 & .and.nactiveline.gt.0)then 210 write(*,*)'TT: face',issurf,'loop in slavintmortar' 211 write(*,*)'area',areaslav,'+',area,'.gt.0.5!' 212 write(20,*)'TT: face',issurf,'loop in slavintmortar' 213 write(20,*)'area',areaslav,'+',area,'.gt.0.5!' 214 nactiveline=0 215 return 216 endif 217 areaslav=areaslav+area 218 if(debug)then 219 write(20,106) k,area,areaslav 220 106 format('tri',i10,' area',e15.8,' atot',e15.8) 221 write(20,*) '# fuer itri',k,'werden 7intp gen.' 222 endif 223! 224! 7 points scheme 225! 226 do i=1,7 227 xil=xit(3)*gauss2d6(1,i)+ 228 & xit(1)*gauss2d6(2,i)+ 229 & xit(2)*(1-gauss2d6(1,i)-gauss2d6(2,i)) 230 231 etl=etat(3)*gauss2d6(1,i)+ 232 & etat(1)*gauss2d6(2,i)+ 233 & etat(2)*(1-gauss2d6(1,i)-gauss2d6(2,i)) 234! 235 call evalshapefunc(xil,etl,xns,nopes,xnl) 236 call evalshapefunc(xil,etl,xl2s,nopes,ps) 237! 238 nintpoint=nintpoint+1 239! 240! projection of the master integration point onto the 241! master surface in order to get the local coordinates 242! own xn for every integration point? 243! 244 call attachline(xl2m,ps,nnodelem,xilm,etlm,xn,p,dist) 245 call evalshapefunc(xilm,etlm,xl2m,nnodelem,pm) 246! 247 if(debug)then 248 ijk=ijk+1 249 write(40,100) ijk,(ps(jj),jj=1,3) 250 ijk=ijk+1 251 write(40,100) ijk,(pm(jj),jj=1,3) 252 write(40,101) ijk-1,ijk-1,ijk 253 endif 254! 255! Calculation of the gap function at the integration point 256! 257 do ii=1,3 258 ph(ii)=pm(ii)-ps(ii) 259 enddo 260 261 al=sqrt(ph(1)**2+ph(2)**2+ph(3)**2) 262 spm=ph(1)*xn(1) 263 & +ph(2)*xn(2) 264 & +ph(3)*xn(3) 265! 266 if(.not.((clearance.gt.1.2357111316d0).and. 267 & (clearance.lt.1.2357111318d0))) then 268! 269! assuming zero displacements at the beginning of the calculation 270! only valid for small tangential movement in the contact zone 271! 272 call evalshapefunc(xil,etl,cl2s,nopes,psc) 273 call evalshapefunc(xilm,etlm,cl2m,nnodelem,pmc) 274 do ii=1,3 275 phc(ii)=pmc(ii)-psc(ii) 276 enddo 277 spmc=phc(1)*xn(1) 278 & +phc(2)*xn(2) 279 & +phc(3)*xn(3) 280 spm=spm-spmc+clearance 281 endif 282 if(shrink)then 283 spm=reltime*spm 284 endif 285 gapmints(nintpoint)=spm 286 pslavsurf(1,nintpoint)=xil 287 pslavsurf(2,nintpoint)=etl 288! 289! weight add to 0.5 \hat{w}_p=A_triaref*w_p=0.5*w_p 290! division by 0.5 to get to original weights... 291! 292 pslavsurf(3,nintpoint)=area*weight2d6(i)/0.5 293 pmastsurf(1,nintpoint)=xilm 294 pmastsurf(2,nintpoint)=etlm 295 imastsurf(nintpoint)=nelemm 296 if(debug)then 297 write(30,201) xil,etl,pslavsurf(3,nintpoint),spm,nelemm 298 endif 299 201 format('xil ',1x,e15.8,2x,'etl',2x,e15.8,2x,'w ',e15.8,2x, 300 & e15.8,2x,'M',i10) 301 enddo 302! 303 enddo 304c write(20,*)'********************' 305 return 306 end 307