SPEED
WRITE_FILE_DGFS.f90
Go to the documentation of this file.
1! Copyright (C) 2012 The SPEED FOUNDATION
2! Author: Ilario Mazzieri
3!
4! This file is part of SPEED.
5!
6! SPEED is free software; you can redistribute it and/or modify it
7! under the terms of the GNU Affero General Public License as
8! published by the Free Software Foundation, either version 3 of the
9! License, or (at your option) any later version.
10!
11! SPEED is distributed in the hope that it will be useful, but
12! WITHOUT ANY WARRANTY; without even the implied warranty of
13! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14! Affero General Public License for more details.
15!
16! You should have received a copy of the GNU Affero General Public License
17! along with SPEED. If not, see <http://www.gnu.org/licenses/>.
18
73
74
75 subroutine write_file_dgfs(nm, sd, tag_mat, cs_nnz_loc, cs_loc, &
76 nn_loc, local_n_num, ne_loc, local_el_num, &
77 xs,ys,zs,&
78 nel_dg_loc, nel_dg_glo, &
79 mpi_id, mpi_comm, mpi_np, &
80 alfa11, alfa12, alfa13, &
81 alfa21, alfa22, alfa23, &
82 alfa31, alfa32, alfa33, &
83 beta11, beta12, beta13, &
84 beta21, beta22, beta23, &
85 beta31, beta32, beta33, &
86 gamma1, gamma2, gamma3, &
87 delta1, delta2, delta3, &
88 faces, area_nodes, dg_els, scratch_dg_els, &
89 filename,mpi_file)
90
91 use max_var
92 use str_mesh
93 use str_mesh_scratch
94
95 implicit none
96
97
98 include 'SPEED.MPI'
99
100
101 type(element), dimension(nel_dg_loc), intent(inout) :: dg_els
102 type(scratch_element), dimension(nel_dg_loc), intent(inout) :: scratch_dg_els
103
104 character*70 :: filename, filempi, filenorm, mpi_file, filempi_new
105
106 integer*4 :: nm, cs_nnz_loc, nn_loc, ne_loc, nel_dg_loc, nel_dg_glo
107 integer*4 :: mpi_comm, mpi_id, mpi_np, dummy
108 integer*4 :: im, nn, ie, ned, yon
109 integer*4 :: ne1, ne2, ne3, ne4, ic1, ic2, ic3, ic4
110 integer*4 :: ne5, ne6, ne7, ne8, ic5, ic6, ic7, ic8
111 integer*4 :: el_conf, face_conf, face_found, unitmpi, unitname, unitnorm
112 integer*4 :: ip, k, j, i, it, ih, ik, tt, indice, node_not_ass, ic, dim1, dim2, jstart
113 integer*4 :: error, ncol
114 integer*4 :: mpierror, nofel
115
116 integer*4, dimension(2) :: dims, dimsfi
117 integer*4, dimension(nm) :: tag_mat, sd
118 integer*4, dimension(0:cs_nnz_loc) :: cs_loc
119 integer*4, dimension(nn_loc) :: local_n_num
120 integer*4, dimension(ne_loc) :: local_el_num
121 integer*4, dimension(mpi_np) :: count_dg
122
123 integer*4, dimension(:,:), allocatable :: DG_CON_ALL
124 integer*4, dimension(3,nel_dg_glo) :: faces
125 integer*4, dimension(nel_dg_glo,3) :: mat_el_fac
126 integer*4, dimension(nel_dg_glo) :: link_faces
127
128 real*8 :: normal_x, normal_y, normal_z, dummyreal
129 real*8 :: c_alfa11, c_alfa12, c_alfa13, c_alfa21, c_alfa22, c_alfa23, c_alfa31, c_alfa32, c_alfa33
130 real*8 :: c_beta11, c_beta12, c_beta13, c_beta21, c_beta22, c_beta23, c_beta31, c_beta32, c_beta33
131 real*8 :: c_gamma1, c_gamma2, c_gamma3, c_delta1, c_delta2, c_delta3
132 real*8 :: xnod, ynod, znod, csi, eta, zeta, xnod1, ynod1, znod1
133 real*8 :: val1, val2, val3, val4, val5, val6, valmin
134 real*8 :: coef_a, coef_b, coef_t, det_trasf
135
136 real*8, dimension(nn_loc) :: xs,ys,zs
137
138 real*8, dimension(:), allocatable :: ctgl,wwgl
139 real*8, dimension(:), allocatable :: ct, ww
140 real*8, dimension(ne_loc) :: alfa11,alfa12,alfa13
141 real*8, dimension(ne_loc) :: alfa21,alfa22,alfa23
142 real*8, dimension(ne_loc) :: alfa31,alfa32,alfa33
143 real*8, dimension(ne_loc) :: beta11,beta12,beta13
144 real*8, dimension(ne_loc) :: beta21,beta22,beta23
145 real*8, dimension(ne_loc) :: beta31,beta32,beta33
146 real*8, dimension(ne_loc) :: gamma1,gamma2,gamma3
147 real*8, dimension(ne_loc) :: delta1,delta2,delta3
148
149 real*8, dimension(:,:), allocatable :: dd
150 real*8, dimension(:,:), allocatable :: dg_pw_all
151 real*8, dimension(25,nel_dg_glo) :: area_nodes
152 real*8, dimension(3,3) :: rot_mat
153 real*8, dimension(nel_dg_glo,3) :: normalxyz
154
155 filempi = 'DGFS000000.mpi'
156 unitmpi = 50
157
158 if (mpi_id .lt. 10) then
159 write(filempi(10:10),'(i1)') mpi_id
160 else if (mpi_id .lt. 100) then
161 write(filempi(9:10),'(i2)') mpi_id
162 else if (mpi_id .lt. 1000) then
163 write(filempi(8:10),'(i3)') mpi_id
164 else if (mpi_id .lt. 10000) then
165 write(filempi(7:10),'(i4)') mpi_id
166 else if (mpi_id .lt. 100000) then
167 write(filempi(6:10),'(i5)') mpi_id
168 else if (mpi_id .lt. 1000000) then
169 write(filempi(5:10),'(i6)') mpi_id
170 endif
171
172 if(len_trim(mpi_file) .ne. 70) then
173 filempi_new = mpi_file(1:len_trim(mpi_file)) // '/' // filempi
174 else
175 filempi_new = filempi
176 endif
177
178 open(unitmpi,file=filempi_new)
179
180
181 filenorm = 'NORMALL.input'; unitnorm = 40
182
183 open(unitnorm,file=filenorm)
184 read(unitnorm,*) nofel
185 do i = 1, nofel
186 read(unitnorm,*) mat_el_fac(i,1),mat_el_fac(i,2),mat_el_fac(i,3) ,&
187 normalxyz(i,1), normalxyz(i,2), normalxyz(i,3), dummy, dummyreal, dummyreal, link_faces(i)
188 enddo
189
190 close(unitnorm)
191
192
193
194
195 ic = 0
196
197 ned = cs_loc(0) - 1
198 node_not_ass = 0
199
200 do it = 1, nel_dg_loc
201
202 if(dg_els(it)%proj_yn .eq. 1) then
203
204
205
206 do i = 1, (dg_els(it)%quad_rule)**2
207
208 tt = 0
209 ih = 1
210
211 do while (tt.eq.0 .and. ih.le. nel_dg_glo)
212
213 if( (faces(1,ih) .ne. dg_els(it)%mat) .and. (link_faces(ih) .eq. dg_els(it)%link) ) then
214
215 !ik = faces(2,ih)
216 !il = faces(3,ih)
217
218 !CHECK THE NORMAL !!!
219 call check_normal(dg_els(it)%nx, dg_els(it)%ny, dg_els(it)%nz, &
220 faces(1,ih), faces(2,ih), faces(3,ih), &
221 nel_dg_glo, normalxyz, mat_el_fac, yon)
222
223
224 if(yon .eq. 1) then
225
226 xnod = scratch_dg_els(it)%x_nq(i)
227 ynod = scratch_dg_els(it)%y_nq(i)
228 znod = scratch_dg_els(it)%z_nq(i)
229
230 call make_bilinear_map(area_nodes(2:25,ih), &
231 c_alfa11, c_alfa12, c_alfa13, &
232 c_alfa21, c_alfa22, c_alfa23, &
233 c_alfa31, c_alfa32, c_alfa33, &
234 c_beta11, c_beta12, c_beta13, &
235 c_beta21, c_beta22, c_beta23, &
236 c_beta31, c_beta32, c_beta33, &
237 c_gamma1, c_gamma2, c_gamma3, &
238 c_delta1, c_delta2, c_delta3)
239
240 call newton_rapson(xnod, ynod, znod, &
241 c_alfa11, c_alfa12, c_alfa13, &
242 c_alfa21, c_alfa22, c_alfa23, &
243 c_alfa31, c_alfa32, c_alfa33, &
244 c_beta11, c_beta12, c_beta13, &
245 c_beta21, c_beta22, c_beta23, &
246 c_beta31, c_beta32, c_beta33, &
247 c_gamma1, c_gamma2, c_gamma3, &
248 c_delta1, c_delta2, c_delta3, &
249 tt, csi, eta, zeta, nofinr, mpi_id,&
250 dg_els(it)%ind_el,faces(2,ih), 1.d-6, 1.01d0, 1)
251
252
253 if(tt == 1) then
254
255 ic = ic + 1
256
257 write(unitmpi,"(1I2,1X,1I12,1X,1I2,1X,1I2,1X,1I12,1X,1I2,3(1X,ES25.16),3(1X,ES25.16))") &
258 dg_els(it)%mat, dg_els(it)%ind_el, dg_els(it)%face_el, &
259 faces(1,ih), faces(2,ih), faces(3,ih), &
260 xnod, ynod, znod, &
261 dg_els(it)%wx_pl(i),dg_els(it)%wy_pl(i),dg_els(it)%wz_pl(i)
262
263 endif ! if (tt == 1)
264
265
266
267 else
268
269 tt = 0
270
271
272 endif ! (yon .eq. 1)
273
274 endif ! (faces(ih,1) .ne. dg_els(it)%mat)
275
276 ih = ih + 1
277
278 enddo ! while tt == 0
279
280 ! CHECK ON NODE NOT ASSIGNED
281 if(tt .eq. 0 .and. ih .gt. nel_dg_glo) then
282 node_not_ass = node_not_ass + 1
283 write(*,*) mpi_id ,'NODE', i,' NOT ASSIGNED', ' el', dg_els(it)%ind_el
284 endif
285
286
287 enddo ! i
288
289 endif !if proj_yn == 1
290
291 enddo ! it
292
293 close(unitmpi)
294
295
296 dims(1) = 6
297 dims(2) = ic
298
299 if(ic .ne. 0) then
300 write(*,'(A,I5,A,I6)') 'Proc. :', mpi_id, ' not assigned nodes : ', node_not_ass
301 write(*,'(A,I5,A,I6)') 'Proc. :', mpi_id, ' assigned nodes : ', ic
302 else
303 write(*,'(A,I5)') 'No projection for Proc. :', mpi_id
304 endif
305
306! write(unitmpi,*) 6, ic
307! do i = 1, ic
308! write(unitmpi,"(1I2,1X,1I12,1X,1I2,1X,1I2,1X,1I12,1X,1I2,3(1X,ES16.9),3(1X,ES16.9))") &
309! DG_CON(1,i), DG_CON(2,i), DG_CON(3,i), DG_CON(4,i), DG_CON(5,i), DG_CON(6,i), &
310! DG_PW(1,i), DG_PW(2,i), DG_PW(3,i), DG_PW(4,i), DG_PW(5,i), DG_PW(6,i)
311! write(unit_mpi,"(1I2,1X,1I12,1X,1I2,1X,1I2,1X,1I12,1X,1I2,3(1X,ES25.16),3(1X,ES25.16))") &
312! write(*,*)
313! enddo
314
315
316
317 call mpi_barrier(mpi_comm, mpierror)
318 call mpi_allgather(dims(2), 1, speed_integer, count_dg, 1, speed_integer, mpi_comm, mpierror)
319
320 dimsfi(1) = 6
321 dimsfi(2) = sum(count_dg)
322
323 ! write(*,*) 'qui in fondo'
324
325
326 if(mpi_id .eq. 0) then
327
328 allocate(dg_con_all(dimsfi(1),dimsfi(2)), dg_pw_all(dimsfi(1),dimsfi(2)))
329
330 do i = 1, mpi_np
331
332
333 if (count_dg(i) .ne. 0) then
334
335 filempi = 'DGFS000000.mpi'
336 unitmpi = 40
337 if (i-1 .lt. 10) then
338 write(filempi(10:10),'(i1)') i-1
339 else if (i-1 .lt. 100) then
340 write(filempi(9:10),'(i2)') i-1
341 else if (i-1 .lt. 1000) then
342 write(filempi(8:10),'(i3)') i-1
343 else if (i-1 .lt. 10000) then
344 write(filempi(7:10),'(i4)') i-1
345 else if (i-1 .lt. 100000) then
346 write(filempi(6:10),'(i5)') i-1
347 else if (i-1 .lt. 1000000) then
348 write(filempi(5:10),'(i6)') i-1
349 endif
350
351 if(len_trim(mpi_file) .ne. 70) then
352 filempi_new = mpi_file(1:len_trim(mpi_file)) // '/' // filempi
353 else
354 filempi_new = filempi
355 endif
356
357 open(unitmpi,file=filempi_new)
358 dim1 = 6
359 dim2 = count_dg(i)
360
361 if(i.eq. 1) then
362 jstart = 0
363 else
364 jstart = sum(count_dg(1:i-1))
365 endif
366
367 do j = 1, dim2
368 read(unitmpi,*) dg_con_all(1,j+jstart), dg_con_all(2,j+jstart), dg_con_all(3,j+jstart), &
369 dg_con_all(4,j+jstart), dg_con_all(5,j+jstart), dg_con_all(6,j+jstart), &
370 dg_pw_all(1,j+jstart), dg_pw_all(2,j+jstart), dg_pw_all(3,j+jstart), &
371 dg_pw_all(4,j+jstart), dg_pw_all(5,j+jstart), dg_pw_all(6,j+jstart)
372
373 ! write(*,*) DG_CON_ALL(1,j+jstart), DG_CON_ALL(2,j+jstart), DG_CON_ALL(3,j+jstart), &
374 ! DG_CON_ALL(4,j+jstart), DG_CON_ALL(5,j+jstart), DG_CON_ALl(6,j+jstart), &
375 ! DG_PW_ALL(1,j+jstart), DG_PW_ALL(2,j+jstart), DG_PW_ALL(3,j+jstart), &
376 ! DG_PW_ALL(4,j+jstart), DG_PW_ALL(5,j+jstart), DG_PW_ALL(6,j+jstart)
377
378
379
380 enddo
381
382 close(unitmpi)
383
384 endif
385
386 enddo
387
388 unitname = 500
389 open(unitname,file=filename)
390 write(unitname,"(1I15)") dimsfi(2)
391
392 do j = 1, dimsfi(2)
393
394 write(unitname,"(1I2,1X,1I12,1X,1I2,1X,1I2,1X,1I12,1X,1I2,3(1X,ES25.16),3(1X,ES25.16))") &
395 dg_con_all(1,j), dg_con_all(2,j), dg_con_all(3,j), &
396 dg_con_all(4,j), dg_con_all(5,j), dg_con_all(6,j), &
397 dg_pw_all(1,j), dg_pw_all(2,j), dg_pw_all(3,j), &
398 dg_pw_all(4,j), dg_pw_all(5,j), dg_pw_all(6,j)
399
400
401
402 enddo
403
404 close(unitname)
405
406
407 deallocate(dg_con_all, dg_pw_all)
408
409 endif
410
411 call mpi_barrier(mpi_comm, error)
412
413
414 return
415
416
417 end subroutine write_file_dgfs
check_normal
subroutine check_normal(nor_x, nor_y, nor_z, material, element, face, nel_glo, normalxyz, mat_el_fac, yon)
Verifies if 2 different faces are opposite or not.
Definition CHECK_NORMAL.f90:38
make_bilinear_map
subroutine make_bilinear_map(nodes, c_alfa11, c_alfa12, c_alfa13, c_alfa21, c_alfa22, c_alfa23, c_alfa31, c_alfa32, c_alfa33, c_beta11, c_beta12, c_beta13, c_beta21, c_beta22, c_beta23, c_beta31, c_beta32, c_beta33, c_gamma1, c_gamma2, c_gamma3, c_delta1, c_delta2, c_delta3)
Takes values from nodes array.
Definition MAKE_BILINEAR_MAP.f90:57
newton_rapson
subroutine newton_rapson(xnod, ynod, znod, alfa11, alfa12, alfa13, alfa21, alfa22, alfa23, alfa31, alfa32, alfa33, beta11, beta12, beta13, beta21, beta22, beta23, beta31, beta32, beta33, gamma1, gamma2, gamma3, delta1, delta2, delta3, tt, csi_s, eta_s, zeta_s, nofi, id_mpi, ipiu, imeno, toll1, toll2, flag)
Performs the Newton Rapson algorithm.
Definition NEWTON_RAPSON.f90:73
write_file_dgfs
subroutine write_file_dgfs(nm, sd, tag_mat, cs_nnz_loc, cs_loc, nn_loc, local_n_num, ne_loc, local_el_num, xs, ys, zs, nel_dg_loc, nel_dg_glo, mpi_id, mpi_comm, mpi_np, alfa11, alfa12, alfa13, alfa21, alfa22, alfa23, alfa31, alfa32, alfa33, beta11, beta12, beta13, beta21, beta22, beta23, beta31, beta32, beta33, gamma1, gamma2, gamma3, delta1, delta2, delta3, faces, area_nodes, dg_els, scratch_dg_els, filename, mpi_file)
Writes file DGFS.input, containing infos for computing integrals on DG interfaces.
Definition WRITE_FILE_DGFS.f90:90
max_var
Set maximal bounds.
Definition MODULES.f90:54
max_var::nofinr
integer, parameter nofinr
max number of newton rapson iterations
Definition MODULES.f90:58
str_mesh_scratch
Contains mesh structure (scratch)
Definition MODULES.f90:67
str_mesh
Contains mesh structure for DG interface elements.
Definition MODULES.f90:84
str_mesh::element
Definition MODULES.f90:88
str_mesh_scratch::scratch_element
Definition MODULES.f90:71