Newer
Older
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
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! ||===\\ |
! || \\ |
! || || //==\\ || || //==|| ||/==\\ |
! || || || || || || || || || || |
! || // || || || || || || || |
! ||===// \\==// \\==\\ \\==\\ || |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! CALCULATE_OCTREE Nov. 2006 |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
subroutine embed_surface_in_octree (osolve,params,surface,is,istep,iter,threadinfo)
!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine ))))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! this routine embed a surface in a new octree, computes a lsf, and then
! refines osolve accordingly and interpolates the lsf on osolve
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine arguments ))))))))))))))))))))
!------------------------------------------------------------------------------|
use threads
use definitions
implicit none
type(octreesolve) osolve
type(parameters) params
type (octreelsf) olsf
type (sheet) surface
integer is
integer istep
integer iter
type (thread) threadinfo
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|
integer,dimension(:),allocatable::levs
type (octreesolve) oint
integer err,ierr,iproc,nproc
integer ie,ilsf,k,i
double precision xxx,yyy,zzz,dist
integer ioctree_number_of_elements
external ioctree_number_of_elements
character*72 :: shift
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
INCLUDE 'mpif.h'
call mpi_comm_size (mpi_comm_world,nproc,ierr)
call mpi_comm_rank (mpi_comm_world,iproc,ierr)
shift=' '
if(iproc.eq.0) write(8,*) 'embedding surface',is,'/',params%ns
! we want to create the octree olsf in which the surface under consideration will be imbedded.
! we first create a uniform octree at level leveluniform_oct
olsf%noctree=params%noctreemax
allocate (olsf%octree(olsf%noctree),stat=threadinfo%err) ; call heap (threadinfo,'olsf%octree', 'embed_surf...',size(olsf%octree),'int',+1)
call octree_init (olsf%octree,olsf%noctree)
call octree_create_uniform (olsf%octree,olsf%noctree,params%leveluniform_oct)
! we then improve the octree to best fit/represent the surface
call build_surface_octree (surface,olsf,params,threadinfo)
if (iproc.eq.0) write (8,*) olsf%nleaves,' leaves in lsf octree'
! we then calculate the connectivity and node coordinates of the surface octree
olsf%nnode=olsf%nleaves*3
allocate (olsf%icon(8,olsf%nleaves),stat=threadinfo%err) ; call heap (threadinfo,'olsf%icon','embed_surf...',size(olsf%icon),'int',+1)
allocate (olsf%x(olsf%nnode),stat=threadinfo%err) ; call heap (threadinfo,'olsf%x','embed_surf...',size(olsf%x),'dp',+1)
allocate (olsf%y(olsf%nnode),stat=threadinfo%err) ; call heap (threadinfo,'olsf%y','embed_surf...',size(olsf%y),'dp',+1)
allocate (olsf%z(olsf%nnode),stat=threadinfo%err) ; call heap (threadinfo,'olsf%z','embed_surf...',size(olsf%z),'dp',+1)
call octree_find_node_connectivity (olsf%octree,olsf%noctree,olsf%icon,olsf%nleaves,olsf%x,olsf%y,olsf%z,olsf%nnode)
write(threadinfo%Logunit,*) '[same] embed_surface_in_octree: olsf%nnode ',olsf%nnode
call flush(threadinfo%Logunit)
allocate (olsf%lsf(olsf%nnode),stat=threadinfo%err) ; call heap (threadinfo,'olsf%lsf', 'embed_surf...',size(olsf%lsf),'dp',+1)
! we then calculate the lsf on the olsf octree
call calculate_lsf (olsf%lsf,olsf%octree,olsf%noctree, &
olsf%icon,olsf%nleaves,olsf%nnode, &
surface%x,surface%y,surface%z, &
surface%xn,surface%yn,surface%zn, &
surface%icon,surface%nt,surface%nsurface, &
params%levelmax_oct,params%normaladvect,is)
if (params%debug>=2) call output_octree_lsf (olsf,surface,is,istep,iter)
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! we create and allocate the oint octree that will receive the current osolve octree.
! this will allow us to resize osolve so that the new surface can be embedded in it.
oint%noctree=osolve%noctree
oint%nleaves=osolve%nleaves
oint%nnode=osolve%nnode
oint%nlsf=osolve%nlsf
allocate (oint%octree(oint%noctree),stat=threadinfo%err) ; call heap (threadinfo,'oint%octree', 'embed_surf...',size(oint%octree),'int',+1)
allocate (oint%icon(8,oint%nleaves),stat=threadinfo%err) ; call heap (threadinfo,'oint%icon', 'embed_surf...',size(oint%icon),'int',+1)
allocate (oint%lsf(oint%nnode,oint%nlsf),stat=threadinfo%err) ; call heap (threadinfo,'oint%lsf', 'embed_surf...',size(oint%lsf),'dp',+1)
allocate (oint%x(oint%nnode),stat=threadinfo%err) ; call heap (threadinfo,'oint%x', 'embed_surf...',size(oint%x),'dp',+1)
allocate (oint%y(oint%nnode),stat=threadinfo%err) ; call heap (threadinfo,'oint%y', 'embed_surf...',size(oint%y),'dp',+1)
allocate (oint%z(oint%nnode),stat=threadinfo%err) ; call heap (threadinfo,'oint%z', 'embed_surf...',size(oint%z),'dp',+1)
oint%octree=osolve%octree
oint%icon=osolve%icon
oint%lsf=osolve%lsf
oint%x=osolve%x
oint%y=osolve%y
oint%z=osolve%z
! having replicated osolve, we can deallocate it
call heap (threadinfo,'osolve%icon','embed_surf...',size(osolve%icon),'int',-1) ; deallocate (osolve%icon)
call heap (threadinfo,'osolve%x','embed_surf...',size(osolve%x),'dp',-1) ; deallocate (osolve%x)
call heap (threadinfo,'osolve%y','embed_surf...',size(osolve%y),'dp',-1) ; deallocate (osolve%y)
call heap (threadinfo,'osolve%z','embed_surf...',size(osolve%z),'dp',-1) ; deallocate (osolve%z)
call heap (threadinfo,'osolve%lsf','embed_surf...',size(osolve%lsf),'dp',-1) ; deallocate (osolve%lsf)
! we reinitialise it
call octree_init (osolve%octree,osolve%noctree)
osolve%octree=oint%octree
! we modify the osolve%octree by embedding the olsf in it
allocate (levs(olsf%nleaves),stat=threadinfo%err) ; call heap (threadinfo,'levs', 'embed_surf...',size(levs),'dp',+1)
call octree_find_element_level (olsf%octree,olsf%noctree,levs,olsf%nleaves)
do ie=1,olsf%nleaves
xxx=sum(olsf%x(olsf%icon(:,ie)))/8.d0
yyy=sum(olsf%y(olsf%icon(:,ie)))/8.d0
zzz=sum(olsf%z(olsf%icon(:,ie)))/8.d0
call octree_create_from_particles (osolve%octree,osolve%noctree,xxx,yyy,zzz,1,levs(ie))
enddo
call heap (threadinfo,'levs', 'embed_surf...',size(levs),'dp',-1) ; deallocate (levs)
! we smoothen osolve
call octree_smoothen (osolve%octree,osolve%noctree)
! the number of leaves in osolve is computed by means of an octree function, and the
! number of nodes is over estimated
osolve%nleaves = ioctree_number_of_elements(osolve%octree,osolve%noctree)
osolve%nnode = osolve%nleaves*3
osolve%nlsf = params%ns
allocate (osolve%icon(8,osolve%nleaves),stat=threadinfo%err) ; call heap (threadinfo,'osolve%icon','embed_surf...',size(osolve%icon),'int',+1)
allocate (osolve%x(osolve%nnode),stat=threadinfo%err) ; call heap (threadinfo,'osolve%x', 'embed_surf...',size(osolve%x),'dp',+1)
allocate (osolve%y(osolve%nnode),stat=threadinfo%err) ; call heap (threadinfo,'osolve%y', 'embed_surf...',size(osolve%y),'dp',+1)
allocate (osolve%z(osolve%nnode),stat=threadinfo%err) ; call heap (threadinfo,'osolve%z', 'embed_surf...',size(osolve%z),'dp',+1)
! we build the icon array of osolve. the routine also outputs the real number of nodes osolve%nnode
call octree_find_node_connectivity (osolve%octree,osolve%noctree, &
osolve%icon,osolve%nleaves, &
osolve%x,osolve%y,osolve%z,osolve%nnode)
! now that osolve%nnode is known we can allocate osolve%lsf
allocate (osolve%lsf(osolve%nnode,osolve%nlsf),stat=threadinfo%err) ; call heap (threadinfo,'osolve%lsf', 'embed_surf...',size(osolve%lsf),'dp',+1)
! if we are dealing with the last surface, nodes are renumbered
if (params%renumber_nodes .and. is.eq.params%ns) then
call octree_renumber_nodes (osolve%icon,osolve%nleaves, &
osolve%x,osolve%y,osolve%z,osolve%nnode)
end if
! the lsf values of the previous is-1 surfaces are known on the oint octree which is no more
! the same as osolve, so we have to interpolates these lsfs on final octree osolve
do ilsf=1,is-1
do k=1,osolve%nnode
call octree_interpolate (oint%octree,oint%noctree,oint%icon,oint%nleaves,oint%lsf(1,ilsf),oint%nnode, &
osolve%x(k),osolve%y(k),osolve%z(k),osolve%lsf(k,ilsf))
! call octree_interpolate3 (oint%octree,oint%noctree,oint%icon,oint%nleaves,oint%lsf(1,ilsf), &
! oint%x(1:oint%nnode),oint%y(1:oint%nnode),oint%z(1:oint%nnode),oint%nnode, &
! osolve%x(k),osolve%y(k),osolve%z(k),osolve%lsf(k,ilsf))
enddo
enddo
! finally the lsf corresponding to the current surface also has to be interpolated onto osolve
! since it was computed on a dedicated olsf
do k=1,osolve%nnode
call octree_interpolate (olsf%octree,olsf%noctree,olsf%icon,olsf%nleaves,olsf%lsf,olsf%nnode, &
osolve%x(k),osolve%y(k),osolve%z(k),osolve%lsf(k,is))
! call octree_interpolate3 (olsf%octree,olsf%noctree,olsf%icon,olsf%nleaves,olsf%lsf,olsf%x(1:olsf%nnode),olsf%y(1:olsf%nnode),olsf%z(1:olsf%nnode),olsf%nnode, &
! osolve%x(k),osolve%y(k),osolve%z(k),osolve%lsf(k,is))
enddo
call heap (threadinfo,'olsf%x','embed_surf...',size(olsf%x),'dp',-1) ; deallocate (olsf%x)
call heap (threadinfo,'olsf%y','embed_surf...',size(olsf%y),'dp',-1) ; deallocate (olsf%y)
call heap (threadinfo,'olsf%z','embed_surf...',size(olsf%z),'dp',-1) ; deallocate (olsf%z)
call heap (threadinfo,'olsf%lsf','embed_surf...',size(olsf%lsf),'dp',-1) ; deallocate (olsf%lsf)
call heap (threadinfo,'olsf%icon','embed_surf...',size(olsf%icon),'int',-1) ; deallocate (olsf%icon)
call heap (threadinfo,'olsf%octree','embed_surf...',size(olsf%octree),'int',-1) ; deallocate (olsf%octree)
call heap (threadinfo,'oint%x','embed_surf...',size(oint%x),'dp',-1) ; deallocate (oint%x)
call heap (threadinfo,'oint%y','embed_surf...',size(oint%y),'dp',-1) ; deallocate (oint%y)
call heap (threadinfo,'oint%z','embed_surf...',size(oint%z),'dp',-1) ; deallocate (oint%z)
call heap (threadinfo,'oint%lsf','embed_surf...',size(oint%lsf),'dp',-1) ; deallocate (oint%lsf)
call heap (threadinfo,'oint%icon','embed_surf...',size(oint%icon),'int',-1) ; deallocate (oint%icon)
call heap (threadinfo,'oint%octree','embed_surf...',size(oint%octree),'int',-1) ; deallocate (oint%octree)
if (params%debug>=1 .and. iproc==0) then
write(*,'(a,i8,a)') shift//'osolve%nleaves=',osolve%nleaves
end if
write(threadinfo%Logunit,*) '[same] build_surface_octree: osolve%nleaves ', osolve%nleaves
call flush(threadinfo%Logunit)
return
end
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|