10 integer,
parameter,
private :: kreal = kind(0.0d0)
11 integer,
parameter,
private :: l_max_surface_node = 20
12 integer,
parameter,
private :: l_max_elem_node = 100
30 real(kind=kreal) :: distance
31 real(kind=kreal) :: wkdist
32 real(kind=kreal) :: lpos(2)
33 real(kind=kreal) :: gpos(3)
34 real(kind=kreal) :: direction(3)
35 real(kind=kreal) :: multiplier(3)
37 real(kind=kreal) :: tangentforce(3)
38 real(kind=kreal) :: tangentforce1(3)
39 real(kind=kreal) :: tangentforce_trial(3)
40 real(kind=kreal) :: tangentforce_final(3)
41 real(kind=kreal) :: reldisp(3)
62 integer,
intent(in) :: fnum
64 write(fnum, *)
"--Contact state=",cstate%state
65 write(fnum, *) cstate%surface, cstate%distance
66 write(fnum, *) cstate%lpos
67 write(fnum, *) cstate%direction
68 write(fnum, *) cstate%multiplier
74 integer,
intent(in) :: etype
75 integer,
intent(in) :: nnode
76 real(kind=kreal),
intent(out) :: mpcval(nnode*3 + 4)
79 real(kind=kreal) :: shapefunc(nnode)
82 mpcval(1:3) = cstate%direction(1:3)
85 mpcval( i*3+j ) = -cstate%direction(j)*shapefunc(i)
88 mpcval( 3*nnode+4 )=cstate%distance
93 fcoeff, symm, stiff, force )
94 integer,
intent(in) :: flag
96 integer,
intent(in) :: etype
97 integer,
intent(in) :: nnode
98 real(kind=kreal),
intent(in) :: ele(3,nnode)
99 real(kind=kreal),
intent(in) :: mu
100 real(kind=kreal),
intent(in) :: mut
101 real(kind=kreal),
intent(in) :: fcoeff
102 logical,
intent(in) :: symm
103 real(kind=kreal),
intent(out) :: stiff(:,:)
104 real(kind=kreal),
intent(out) :: force(:)
107 real(kind=kreal) :: shapefunc(nnode)
108 real(kind=kreal) :: n(nnode*3+3), dispmat(2,nnode*3+3)
109 real(kind=kreal) :: metric(2,2)
110 real(kind=kreal) :: det, inverse(2,2), ff(2), cff(2)
111 real(kind=kreal) :: dum11(nnode*3+3,nnode*3+3), dum12(nnode*3+3,nnode*3+3)
112 real(kind=kreal) :: dum21(nnode*3+3,nnode*3+3), dum22(nnode*3+3,nnode*3+3)
113 real(kind=kreal) :: tangent(3,2)
116 n(1:3) = cstate%direction(1:3)
118 n( i*3+1:i*3+3 ) = -shapefunc(i)*cstate%direction(1:3)
120 forall( i=1:nnode*3+3, j=1:nnode*3+3 )
121 stiff(i,j) = mu* n(i)*n(j)
123 force(1:nnode*3+3) = n(:)
126 call dispincrematrix( cstate%lpos, etype, nnode, ele, tangent, metric, dispmat )
129 if( fcoeff/=0.d0 )
then
130 forall(i=1:nnode*3+3, j=1:nnode*3+3)
131 dum11(i,j) = mut*dispmat(1,i)*dispmat(1,j)
132 dum12(i,j) = mut*dispmat(1,i)*dispmat(2,j)
133 dum21(i,j) = mut*dispmat(2,i)*dispmat(1,j)
134 dum22(i,j) = mut*dispmat(2,i)*dispmat(2,j)
136 stiff(1:nnode*3+3,1:nnode*3+3) &
137 = stiff(1:nnode*3+3,1:nnode*3+3) &
138 + metric(1,1)*dum11 + metric(1,2)*dum12 &
139 + metric(2,1)*dum21 + metric(2,2)*dum22
142 det = metric(1,1)*metric(2,2)-metric(1,2)*metric(2,1)
143 if( det==0.d0 ) stop
"Math error in contact stiff calculation"
144 inverse(1,1) = metric(2,2)/det
145 inverse(2,1) = -metric(2,1)/det
146 inverse(1,2) = -metric(1,2)/det
147 inverse(2,2) = metric(1,1)/det
148 ff(:) = cstate%multiplier(2:3)
149 cff(:) = matmul( inverse, ff )
150 ff(:) = ff(:)/dsqrt( ff(1)*ff(1)+ff(2)*ff(2) )
151 cff(:) = cff(:)/dsqrt( cff(1)*cff(1)+cff(2)*cff(2) )
152 stiff(1:nnode*3+3,1:nnode*3+3) = stiff(1:nnode*3+3,1:nnode*3+3) - &
153 ( cff(1)*ff(1)*metric(1,1)+ cff(2)*ff(1)*metric(1,2) )*dum11 - &
154 ( cff(2)*ff(2)*metric(1,2)+ cff(1)*ff(2)*metric(1,1) )*dum21 - &
155 ( cff(1)*ff(1)*metric(1,2)+ cff(2)*ff(1)*metric(2,2) )*dum12 - &
156 ( cff(2)*ff(2)*metric(2,2)+ cff(1)*ff(2)*metric(1,2) )*dum22
164 real(kind=kreal),
intent(in) :: pos(2)
165 integer,
intent(in) :: etype
166 real(kind=kreal),
intent(in) :: ele(:,:)
167 real(kind=kreal),
intent(out) :: tensor(2,2)
170 real(kind=kreal) :: tangent(3,2)
173 tensor(1,1)= dot_product( tangent(:,1), tangent(:,1) )
174 tensor(1,2)= dot_product( tangent(:,1), tangent(:,2) )
175 tensor(2,1)= dot_product( tangent(:,2), tangent(:,1) )
176 tensor(2,2)= dot_product( tangent(:,2), tangent(:,2) )
182 real(kind=kreal),
intent(in) :: pos(2)
183 integer,
intent(in) :: etype
184 integer,
intent(in) :: nnode
185 real(kind=kreal),
intent(in) :: ele(3,nnode)
186 real(kind=kreal),
intent(out) :: tangent(3,2)
187 real(kind=kreal),
intent(out) :: tensor(2,2)
188 real(kind=kreal),
intent(out) :: matrix(2,nnode*3+3)
191 real(kind=kreal) :: det
192 real(kind=kreal) :: shapefunc(nnode), t1(nnode*3+3), t2(nnode*3+3)
193 call tangentbase( etype, nnode, pos, ele, tangent )
194 tensor(1,1)= dot_product( tangent(:,1), tangent(:,1) )
195 tensor(1,2)= dot_product( tangent(:,1), tangent(:,2) )
196 tensor(2,1)= dot_product( tangent(:,2), tangent(:,1) )
197 tensor(2,2)= dot_product( tangent(:,2), tangent(:,2) )
198 det = tensor(1,1)*tensor(2,2)-tensor(1,2)*tensor(2,1)
199 if( det==0.d0 ) stop
"Error in calculate DispIncreMatrix"
207 t1( j ) = tangent(j,1)
208 t2( j ) = tangent(j,2)
210 forall( i=1:nnode, j=1:3 )
211 t1( i*3+j ) = -tangent(j,1)*shapefunc(i)
212 t2( i*3+j ) = -tangent(j,2)*shapefunc(i)
215 matrix(1,:) = (tensor(2,2)*t1(:)-tensor(1,2)*t2(:))/det
216 matrix(2,:) = (tensor(1,1)*t2(:)-tensor(2,1)*t1(:))/det
217 tangent(:,1) = tangent(:,1)/dsqrt(dot_product(tangent(:,1),tangent(:,1)))
218 tangent(:,2) = tangent(:,2)/dsqrt(dot_product(tangent(:,2),tangent(:,2)))
223 real(kind=kreal),
intent(in) :: xyz(3)
224 integer,
intent(in) :: etype
225 integer,
intent(in) :: nn
226 real(kind=kreal),
intent(in) :: elemt(3,nn)
227 real(kind=kreal),
intent(in) :: reflen
229 logical,
intent(out) :: isin
230 real(kind=kreal),
intent(in) :: distclr
231 real(kind=kreal),
optional :: ctpos(2)
232 real(kind=kreal),
optional :: localclr
234 integer :: count,order, initstate
235 real(kind=kreal) :: determ, inverse(2,2)
236 real(kind=kreal) :: sfunc(nn), curv(3,2,2)
237 real(kind=kreal) :: r(2), dr(2), r_tmp(2)
238 real(kind=kreal) :: xyz_out(3)
239 real(kind=kreal) :: dist_last,dist_now, dxyz(3)
240 real(kind=kreal) :: tangent(3,2)
241 real(kind=kreal) :: df(2),d2f(2,2),normal(3)
242 real(kind=kreal),
parameter :: eps = 1.0d-8
243 real(kind=kreal) :: clr, tol, factor
245 initstate = cstate%state
247 if(
present( localclr ) ) clr=localclr
248 if(
present( ctpos ) )
then
257 xyz_out = matmul( elemt(:,:), sfunc )
258 dxyz(1:3) = xyz_out(1:3) - xyz(1:3)
259 dist_last = dot_product( dxyz, dxyz(:) )
262 call curvature( etype, nn, r, elemt, curv )
265 df(1:2) = -matmul( dxyz(:), tangent(:,:) )
267 d2f(1,1)= dot_product( tangent(:,1), tangent(:,1) ) - dot_product( dxyz, curv(:,1,1) )
268 d2f(1,2)= dot_product( tangent(:,1), tangent(:,2) ) - dot_product( dxyz, curv(:,1,2) )
269 d2f(2,1)= dot_product( tangent(:,2), tangent(:,1) ) - dot_product( dxyz, curv(:,2,1) )
270 d2f(2,2)= dot_product( tangent(:,2), tangent(:,2) ) - dot_product( dxyz, curv(:,2,2) )
273 determ = d2f(1,1)*d2f(2,2) - d2f(1,2)*d2f(2,1)
274 if( determ==0.d0 ) stop
"Math error in contact searching"
275 inverse(1,1) = d2f(2,2) / determ
276 inverse(2,2) = d2f(1,1) / determ
277 inverse(1,2) = -d2f(1,2) / determ
278 inverse(2,1) = -d2f(2,1) / determ
279 dr=matmul(inverse,df)
281 tol=dot_product(dr,dr)
282 if( dsqrt(tol)> 3.d0 )
then
288 r_tmp(1:2) = r(1:2) + factor*dr(1:2)
290 xyz_out(1:3) = matmul( elemt(:,:), sfunc(:) )
291 dxyz(1:3) = xyz(1:3)-xyz_out(:)
292 dist_now = dot_product( dxyz, dxyz )
293 if(dist_now <= dist_last)
exit
294 factor = factor*0.7d0
304 dxyz(:)=xyz_out(:)-xyz(:)
306 normal(:) = normal(:)/dsqrt( dot_product(normal, normal) )
308 if( dabs(normal(count))<1.d-10 ) normal(count) =0.d0
309 if( dabs(1.d0-dabs(normal(count)))<1.d-10 ) normal(count) =sign(1.d0, normal(count))
311 cstate%distance = dot_product( dxyz, normal )
313 if( cstate%distance < distclr*reflen .and. cstate%distance > -5.0d-01*reflen ) isin = .true.
319 cstate%state = initstate
321 cstate%gpos(:)=xyz_out(:)
323 cstate%direction(:) = normal(:)
324 cstate%wkdist = cstate%distance
331 integer,
intent(in) :: etype
332 integer,
intent(in) :: nn
333 real(kind=kreal),
intent(in) :: elemt0(3,nn)
334 real(kind=kreal),
intent(in) :: elemt(3,nn)
338 real(kind=kreal) :: tangent0(3,2), tangent(3,2)
339 real(kind=kreal) :: coeff(2), norm, norm_tmp
340 real(kind=kreal) :: tangentforce_tmp(3)
342 call tangentbase( etype, nn, cstate%lpos, elemt0, tangent0 )
343 call tangentbase( etype, nn, cstate%lpos, elemt, tangent )
347 coeff(i) = dot_product(cstate%tangentForce(1:3),tangent0(1:3,i))
348 coeff(i) = coeff(i)/dot_product(tangent0(1:3,i),tangent0(1:3,i))
350 tangentforce_tmp(1:3) = coeff(1)*tangent0(1:3,1) + coeff(2)*tangent0(1:3,2)
351 norm_tmp = dsqrt(dot_product(tangentforce_tmp,tangentforce_tmp))
353 if( norm_tmp > 1.d-6 )
then
354 norm = dsqrt(dot_product(cstate%tangentForce,cstate%tangentForce))
355 coeff(1:2) = (norm/norm_tmp)*coeff(1:2)
359 cstate%tangentForce1(1:3) = coeff(1)*tangent(1:3,1) + coeff(2)*tangent(1:3,2)
This module encapsulate the basic functions of all elements provide by this software.
subroutine getshapefunc(fetype, localcoord, func)
Calculate the shape function in natural coodinate system.
integer(kind=kind(2)) function getnumberofnodes(etype)
Obtain number of nodes of the element.
subroutine getelementcenter(fetype, localcoord)
Return natural coordinate of the center of surface element.
subroutine tangentbase(fetype, nn, localcoord, elecoord, tangent)
Calculate base vector of tangent space of 3d surface.
integer function isinsideelement(fetype, localcoord, clearance)
if a point is inside a surface element -1: No; 0: Yes; >0: Node's (vertex) number
real(kind=kreal) function, dimension(3) surfacenormal(fetype, nn, localcoord, elecoord)
Calculate normal of 3d-surface.
subroutine curvature(fetype, nn, localcoord, elecoord, l2ndderiv, normal, curv)
Calculate curvature tensor at a point along 3d surface.