-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathintegrate.f90
127 lines (93 loc) · 3.72 KB
/
integrate.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
module integrate
implicit none
private
public lockface
contains
! ######################################################################
subroutine lockface(itabsrt,itabsrtface,icoord,coord,facelst)
! ----------------------------------------------------------------------
! Given sorted notes itabsrt (element) and itabsrtface
! find icoord, coord indicating which face of itabsrt is the same
! as itabsrtface.
! Note that the elface below are in order so that parameters eg (xi,eta)
! on face increase from 1st to 3rd point. Thus compatible
! with face integration
! ----------------------------------------------------------------------
use iofile, only : prterr
integer itabsrt(8,3),itabsrtface(4,3),facelst(4),icoord
double precision coord
integer elface(6,4),icoords(6),itriplet(3),jtriplet(3)
integer i,j,ielface,starti,startj,tot,ind
double precision coords(6)
elface(1,:)=[2,3,7,6]; icoords(1)=1; coords(1)=+1.0
elface(2,:)=[1,4,8,5]; icoords(2)=1; coords(2)=-1.0
elface(3,:)=[4,3,7,8]; icoords(3)=2; coords(3)=+1.0
elface(4,:)=[1,2,6,5]; icoords(4)=2; coords(4)=-1.0
elface(5,:)=[5,6,7,8]; icoords(5)=3; coords(5)=+1.0
elface(6,:)=[1,2,3,4]; icoords(6)=3; coords(6)=-1.0
do ielface=1,6
! ----------------------------------------------------------------------
! for ielface, find first equivalent point
! ----------------------------------------------------------------------
starti=-1
outer: do i=1,4
itriplet=itabsrtface(i,:)
do j=1,4
jtriplet=itabsrt(elface(ielface,j),:)
if (all(itriplet.eq.jtriplet)) then
starti=i
startj=j
exit outer
endif
enddo
enddo outer
! ----------------------------------------------------------------------
! loop over a face and check if all triplets are the same
! between this face of the element and the face
! ----------------------------------------------------------------------
if (starti.ne.-1) then
tot=0
do ind=0,3
i=starti+ind
j=startj+ind
if (i.gt.4) i=i-4
if (i.lt.1) i=i+4
if (j.gt.4) j=j-4
if (j.lt.1) j=j+4
itriplet=itabsrtface(i,:)
jtriplet=itabsrt(elface(ielface,j),:)
if (all(itriplet.eq.jtriplet)) tot=tot+1
enddo
if (tot.eq.4) then
icoord=icoords(ielface)
coord = coords(ielface)
facelst=elface(ielface,:)
return
endif
! ----------------------------------------------------------------------
! Try it again in the opposite direction
! ----------------------------------------------------------------------
tot=0
do ind=0,3
i=starti+ind
j=startj-ind
if (i.gt.4) i=i-4
if (i.lt.1) i=i+4
if (j.gt.4) j=j-4
if (j.lt.1) j=j+4
itriplet=itabsrtface(i,:)
jtriplet=itabsrt(elface(ielface,j),:)
if (all(itriplet.eq.jtriplet)) tot=tot+1
enddo
if (tot.eq.4) then
icoord=icoords(ielface)
coord = coords(ielface)
facelst=elface(ielface,:)
return
endif
endif
enddo
call prterr('Did not find equivalent face: surface integrals must "point" into the simulation volume')
end subroutine lockface
! ######################################################################
end module integrate