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! regularization function for normal contact mortar 21! see phd-thesis Sitzmann Chapter 3.2.1., semi-smooth Newton for normal contact 22! 23! [in]lambdap contact pressure in normal direction 24! [in]divmode indicates whether function or derivate 25! should be called 26! =0 function called 27! =1 derivative called 28! [in]regmode selects regularization funtion 29! =1 perturbed Lagrange 30! =2 piece wise linear with given data points 31! =3 expontential contact law 32! [out]gnc result regularization function 33! [in] aninvloc stiffness constant for perturbed Lagrange 34! [in] p0 parameter for exponential regularization 35! [in] beta parameter for exponential regularization 36! 37 subroutine regularization_gn_c(lambdap,divmode,regmode, 38 & gnc,aninvloc,p0,beta,elcon,nelcon,itie,ntmat_, 39 & plicon,nplicon,npmat_,ncmat_,tietol,scal) 40! 41! regularization function for normal contact 42! Author: Saskia Sitzmann 43! 44 implicit none 45! 46 integer divmode,regmode,i,ndata,kode,npmat_,ncmat_, 47 & itie,ntmat_,nplicon(0:ntmat_,*),nelcon(2,*), 48 & imat 49! 50 real*8 lambdap,gnc,pn_d(40),gn_d(40),aninv1(40),t(40), 51 & beta,p0,aninvloc,elconloc(21),plconloc(802),t1l, 52 & elcon(0:ncmat_,ntmat_,*),plicon(0:2*npmat_,ntmat_,*), 53 & tietol(3,*),scal 54! 55! 56! 57 kode=-51 58 t1l=0.0 59 imat=int(tietol(2,itie+1)) 60! 61 gnc=0.0 62! 63! perturbed Lagrange 64! 65 if(regmode.eq.1)then 66 if(divmode.eq.0)then 67 gnc=aninvloc*lambdap 68 elseif(divmode.eq.1)then 69 gnc=aninvloc 70 else 71 write(*,*)'error in regularzation_gn_c.f!' 72 call exit(201) 73 endif 74! 75! multiple perturbed Lagrange 76! 77 else if(regmode.eq.2)then 78 call materialdata_sp(elcon,nelcon,imat,ntmat_,i,t1l, 79 & elconloc,kode,plicon,nplicon,npmat_,plconloc,ncmat_) 80 ndata=int(plconloc(801)) 81! is this right? 82c ndata=int(plconloc(81)) 83 do i=1,ndata 84 gn_d(i)=plconloc(2*i-1)*scal 85 pn_d(i)=plconloc(2*i)*scal 86 enddo 87 do i=1,ndata-1 88 aninv1(i)=(gn_d(i+1)-gn_d(i))/(pn_d(i+1)-pn_d(i)) 89 t(i)=gn_d(i)-aninv1(i)*pn_d(i) 90 enddo 91 if(divmode.eq.0)then 92 if(lambdap.lt.pn_d(1))then 93 gnc=aninv1(1)*lambdap+t(1) 94 endif 95 do i=1,ndata-2 96 if(pn_d(i+1).gt.lambdap.and.pn_d(i).le.lambdap)then 97 gnc=aninv1(i)*lambdap+t(i) 98 endif 99 enddo 100 if(pn_d(ndata-1).le.lambdap)then 101 gnc=aninv1(ndata-1)*lambdap+t(ndata-1) 102 endif 103 elseif(divmode.eq.1)then 104 if(lambdap.lt.pn_d(1))then 105 gnc=aninv1(1) 106 endif 107 do i=1,ndata-2 108 if(pn_d(i+1).gt.lambdap.and.pn_d(i).le.lambdap)then 109 gnc=aninv1(i) 110 endif 111 enddo 112 if(pn_d(ndata-1).le.lambdap)then 113 gnc=aninv1(ndata-1) 114 endif 115 else 116 write(*,*)'error in regularzation_gn_c.f!' 117 call exit(201) 118 endif 119! 120! exponetial perturbed Lagrange 121! 122 else if (regmode.eq.3)then 123 if(divmode.eq.0)then 124 if(lambdap.gt.0.0)then 125 gnc=scal*beta*log(((lambdap/scal)+p0)/p0) 126 else 127 gnc=beta/(p0)*lambdap 128 endif 129 elseif(divmode.eq.1)then 130 if(lambdap.gt.0.0)then 131 gnc=beta/(lambdap+p0) 132 else 133 gnc=beta/(p0) 134 endif 135 else 136 write(*,*)'error in regularzation_gn_c.f!' 137 call exit(201) 138 endif 139 endif 140c call exit(201) 141 return 142 end 143