FE-Project
Loading...
Searching...
No Matches
scale_meshutil_2d.F90
Go to the documentation of this file.
1!-------------------------------------------------------------------------------
9!-------------------------------------------------------------------------------
10#include "scaleFElib.h"
11module scale_meshutil_2d
12 !-----------------------------------------------------------------------------
13 !
14 !++ used modules
15 !
16 use scale_precision
17
18 !-----------------------------------------------------------------------------
19 implicit none
20 private
21
22 !-----------------------------------------------------------------------------
23 !
24 !++ Public procedure
25 !
26 public :: meshutil2d_genrectdomain
27 public :: meshutil2d_genconnectivity
28 public :: meshutil2d_buildinteriormap
29 public :: meshutil2d_genpatchboundarymap
30 public :: meshutil2d_buildglobalmap
31
32contains
33!OCL SERIAL
34 subroutine meshutil2d_genrectdomain( pos_v, EToV, &
35 Ke_x, xmin, xmax, Ke_y, ymin, ymax )
36
37 implicit none
38
39 integer, intent(in) :: ke_x
40 real(rp), intent(in) :: xmin, xmax
41 integer, intent(in) :: ke_y
42 real(rp), intent(in) :: ymin, ymax
43
44 real(rp), intent(out) :: pos_v((ke_x+1)*(ke_y+1),2)
45 integer, intent(out) :: etov(ke_x*ke_y,4)
46
47 integer :: i, j
48 integer :: n
49 integer :: nvx, nvy
50 !-----------------------------------------------------------------------------
51
52 nvx = ke_x + 1
53 nvy = ke_y + 1
54
55 do j=1, nvy
56 do i=1, nvx
57 n = i + (j-1)*nvx
58 pos_v(n,1) = (xmax - xmin)*dble(i - 1)/dble(nvx - 1) + xmin
59 pos_v(n,2) = (ymax - ymin)*dble(j - 1)/dble(nvy - 1) + ymin
60 end do
61 end do
62
63 do j=1, ke_y
64 do i=1, ke_x
65 n = i + (j-1)*ke_x
66 etov(n,1:4) = (j-1)*nvx + &
67 & (/ i, i+1, &
68 & i + nvx, i+1 + nvx /)
69 end do
70 end do
71
72 !---
73 !!$
74 !!$ write(*,*) "-- vx, vy --"
75 !!$ do j=1, NvY
76 !!$ do i=1, NvX
77 !!$ n = i + (j-1)*NvY
78 !!$ write(*,*) i, j, " :", mesh%vx(n), mesh%vy(n)
79 !!$ end do
80 !!$ end do
81 !!$
82 !!$ write(*,*) "-- EToV --"
83 !!$ do j=1, Ke_y
84 !!$ do i=1, Ke_x
85 !!$ n = i + (j-1)*Ke_y
86 !!$ write(*,*) i, j, " :", mesh%EToV(n,:)
87 !!$ end do
88 !!$ end do
89
90 return
91 end subroutine meshutil2d_genrectdomain
92
93!OCL SERIAL
94 subroutine meshutil2d_genconnectivity( EToE, EToF, &
95 EToV, Ne, Nfaces )
96
97 use scale_quicksort, only: quicksort_exec_with_idx
98 implicit none
99
100 integer, intent(in) :: ne
101 integer, intent(in) :: nfaces
102 integer, intent(out) :: etoe(ne, nfaces)
103 integer, intent(out) :: etof(ne, nfaces)
104 integer, intent(in) :: etov(ne,4)
105
106 integer :: nodes(ne*nfaces,2)
107 integer(kind=8) :: face_ids(ne*nfaces)
108 integer :: ke
109 integer :: f
110 integer :: n
111 integer :: n1, n2
112 integer :: nnodes
113 integer :: tmp
114 integer :: nnodes_row
115 integer :: matchl(2,3), matchr(2,3)
116
117 real(rp) :: etoe_1d(ne*nfaces)
118 real(rp) :: etof_1d(ne*nfaces)
119
120 integer :: spnodetonode(3,ne*nfaces)
121 integer :: spnodetonoderowtmp(3,ne*nfaces)
122 integer :: sort_indx(ne*nfaces)
123 integer(kind=8) :: sorted_faceid(ne*nfaces)
124 !-----------------------------------------------------------------------------
125
126 nnodes = maxval( etov )
127 nnodes_row = size(nodes,1)
128
129 !---------
130 do ke=1, ne
131 nodes(ke ,:) = etov(ke,(/ 1, 2 /))
132 nodes(ke+ne ,:) = etov(ke,(/ 2, 4 /))
133 nodes(ke+2*ne,:) = etov(ke,(/ 4, 3 /))
134 nodes(ke+3*ne,:) = etov(ke,(/ 3, 1 /))
135 end do
136
137 ! Sort
138 do n=1, nnodes_row
139 if (nodes(n,1) > nodes(n,2) ) then
140 tmp = nodes(n,1);
141 nodes(n,1) = nodes(n,2)
142 nodes(n,2) = tmp
143 end if
144 ! write(*,*) n, ":", nodes(n,:)
145 end do
146 nodes = nodes - 1
147 !---------
148
149 do ke=1, ne
150 etoe(ke,:) = ke
151 etof(ke,:) = (/ 1, 2, 3, 4 /)
152 end do
153
154 face_ids(:) = nodes(:,1)*nnodes + nodes(:,2) + 1
155
156 do f=1, nfaces
157 do ke=1, ne
158 n = ke + (f-1)*ne
159 spnodetonode(:,n) = (/ n, etoe(ke,f), etof(ke,f) /)
160 sorted_faceid(n) = face_ids(n)
161 sort_indx(n) = n
162 ! write(*,*) "face_id, n, EToE, EToF:", spNodeToNode(:,n)
163 end do
164 end do
165
166 !- sort row
167 call quicksort_exec_with_idx( ne*nfaces, sorted_faceid, sort_indx )
168 spnodetonoderowtmp(:,:) = spnodetonode(:,:)
169 do n=1, nnodes_row
170 spnodetonode(:,n) = spnodetonoderowtmp(:,sort_indx(n))
171 ! write(*,'(a,4i10)') "(sorted) face_id, n, EToE, EToF:", spNodeToNode(:,n)
172 end do
173
174 etoe_1d(:) = -1
175 etof_1d(:) = -1
176 do n=1, nnodes_row-1
177 if ( sorted_faceid(n) - sorted_faceid(n+1) == 0 ) then
178 matchl(:,:) = transpose( spnodetonode(:,(/ n, n+1 /)) )
179 matchr(:,:) = transpose( spnodetonode(:,(/ n+1, n /)) )
180
181 etoe_1d(matchl(:,1)) = matchr(:,2)
182 etof_1d(matchl(:,1)) = matchr(:,3)
183 end if
184 end do
185
186 do f=1, nfaces
187 do ke=1, ne
188 n = ke + (f-1)*ne
189 if ( etoe_1d(n) /= -1 ) then
190 etoe(ke,f) = etoe_1d(n)
191 etof(ke,f) = etof_1d(n)
192 end if
193 end do
194 end do
195
196 !-------------------
197
198! write(*,*) "EToE------"
199! do k=1, mesh%Ne
200! write(*,*) "k=", k, ":", mesh%EToE(k,:)
201! end do
202! write(*,*) "EToF------"
203! do k=1, mesh%Ne
204! write(*,*) "k=", k, ":", mesh%EToF(k,:)
205! end do
206
207 return
208 end subroutine meshutil2d_genconnectivity
209
210!OCL SERIAL
211 subroutine meshutil2d_buildinteriormap( VMapM, VMapP, MapM, MapP, &
212 pos_en, pos_ev, EtoE, EtoF, EtoV, Fmask, Ne, Np, Nfp, Nfaces, Nv)
213
214 implicit none
215
216 integer, intent(in) :: ne
217 integer, intent(in) :: np
218 integer, intent(in) :: nfp
219 integer, intent(in) :: nfaces
220 integer, intent(in) :: nv
221
222 integer, intent(out) :: vmapm(nfp,nfaces,ne)
223 integer, intent(out) :: vmapp(nfp,nfaces,ne)
224 integer, intent(out) :: mapm(nfp,nfaces,ne)
225 integer, intent(out) :: mapp(nfp,nfaces,ne)
226
227 real(rp), intent(in) :: pos_en(np,ne,2)
228 real(rp), intent(in) :: pos_ev(nv,2)
229 integer, intent(in) :: etoe(ne, nfaces)
230 integer, intent(in) :: etof(ne, nfaces)
231 integer, intent(in) :: etov(ne, 4)
232 integer, intent(in) :: fmask(nfp,4)
233
234
235 integer :: ke, ke1, ke2
236 integer :: f, f1, f2
237 integer :: p
238 integer :: n
239 integer :: idp, idm
240 integer :: v1, v2
241 integer :: nodeids(np,ne)
242 real(rp) :: x1(nfp,nfp), x2(nfp,nfp)
243 real(rp) :: y1(nfp,nfp), y2(nfp,nfp)
244 real(rp) :: dist(nfp,nfp)
245 real(rp) :: x(np*ne), y(np*ne)
246
247 integer :: mindist_indx(1)
248 !-----------------------------------------------------------------------------
249
250 !$omp parallel private(ke,f,p,n)
251 !$omp do
252 do ke=1, ne
253 do p=1, np
254 n = p + (ke-1)*np
255 nodeids(p,ke) = n
256 x(n) = pos_en(p,ke,1)
257 y(n) = pos_en(p,ke,2)
258 end do
259 end do
260 !$omp end do
261 !$omp do collapse(2)
262 do ke=1, ne
263 do f=1, nfaces
264 do p=1, nfp
265 n = p + (f-1)*nfp + (ke-1)*nfp*nfaces
266 mapm(p,f,ke) = n
267 mapp(p,f,ke) = n
268 vmapm(p,f,ke) = nodeids(fmask(p,f),ke)
269 end do
270 end do
271 end do
272
273 !$omp workshare
274 vmapp(:,:,:) = -1
275 !$omp end workshare
276 !$omp end parallel
277
278 !$omp parallel private( &
279 !$omp ke1, f1, ke2, f2, v1, v2, &
280 !$omp x1, x2, y1, y2, dist, mindist_indx, idP, idM )
281
282 !$omp do
283 do ke1=1, ne
284 do f1=1, nfaces
285 ke2 = etoe(ke1,f1); f2 = etof(ke1,f1)
286
287 v1 = etov(ke1,f1); v2 = etov(ke1,1+mod(f1,nfaces))
288
289 x1(:,:) = spread( x(vmapm(:,f1,ke1)), 2, nfp )
290 x2(:,:) = spread( x(vmapm(:,f2,ke2)), 1, nfp )
291 y1(:,:) = spread( y(vmapm(:,f1,ke1)), 2, nfp )
292 y2(:,:) = spread( y(vmapm(:,f2,ke2)), 1, nfp )
293
294 dist(:,:) = (x1(:,:) - x2(:,:))**2 &
295 + (y1(:,:) - y2(:,:))**2
296 do idm=1, nfp
297 mindist_indx(:) = minloc(dist(idm,:))
298 idp = mindist_indx(1)
299 vmapp(idm,f1,ke1) = vmapm(idp,f2,ke2)
300 mapp(idm,f1,ke1) = idp + (f2-1)*nfp + (ke2-1)*nfp*nfaces
301 end do
302 end do
303 end do
304 !$omp end do
305 !$omp end parallel
306
307 !-----
308 ! mapB_counter = 0
309 ! do k=1,mesh%Ne
310 ! do f=1,elem%Nfaces
311 ! do p=1,elem%Nfp
312 ! n = p + (f-1)*elem%Nfp + (k-1)*elem%Nfp*elem%Nfaces
313 ! if (mesh%VMapM(p,f,k) == mesh%VMapP(p,f,k)) then
314 ! mapB_counter = mapB_counter + 1
315 ! mapB_tmp(mapB_counter) = n
316 !
317 ! vmapB_tmp(mapB_counter) = mesh%VMapM(p,f,k)
318 ! mesh%VMapP(p,f,k) = elem%Np*mesh%NeE + mapB_counter
319 ! end if
320 ! end do
321 ! end do
322 ! end do
323 !
324 ! allocate( mesh%mapB(mapB_counter) )
325 ! allocate( mesh%vmapB(mapB_counter) )
326 ! mesh%mapB(:) = mapB_tmp(1:mapB_counter)
327 ! mesh%vmapB(:) = vmapB_tmp(1:mapB_counter)
328
329 !-------
330 ! write(*,*) "Build MapInfo: "
331 ! do k=mesh%NeS,mesh%NeE
332 ! write(*,*) "k=", k, "---"
333 ! write(*,*) " VMapM:", mesh%VMapM(:,:,k)
334 ! write(*,*) " VMapP:", mesh%VMapP(:,:,k)
335 ! end do
336
337 ! write(*,*) "mapB:", mesh%mapB(:)
338 ! write(*,*) "vmapB:", mesh%vmapB(:)
339
340 end subroutine meshutil2d_buildinteriormap
341
342!OCL SERIAL
343 subroutine meshutil2d_genpatchboundarymap( VMapB, MapB, VMapP, &
344 pos_en, xmin, xmax, ymin, ymax, Fmask, Ne, Np, Nfp, Nfaces, Nv)
345
346 implicit none
347
348 integer, intent(in) :: ne
349 integer, intent(in) :: np
350 integer, intent(in) :: nfp
351 integer, intent(in) :: nfaces
352 integer, intent(in) :: nv
353
354 integer, intent(inout), allocatable :: vmapb(:)
355 integer, intent(inout), allocatable :: mapb(:)
356 integer, intent(inout) :: vmapp(nfp,nfaces,ne)
357
358 real(rp), intent(in) :: pos_en(np,ne,2)
359 real(rp), intent(in) :: xmin, xmax
360 real(rp), intent(in) :: ymin, ymax
361 integer, intent(in) :: fmask(nfp,4)
362
363
364 integer :: ke
365 integer :: b
366 integer :: f
367 integer :: i, j
368 real(rp) :: x, y
369
370 real(rp), parameter :: node_tol = 1.0e-12_rp
371
372 real(rp) :: ordinfo(nfp*ne,4)
373 integer :: elemids(nfp*ne,4)
374 integer :: faceids(nfp*ne,4)
375 integer :: counterb(4)
376 integer :: mapb_counter
377 real(rp) :: rdomx, rdomy
378 !-----------------------------------------------------------------------------
379
380
381 counterb(:) = 0
382 rdomx = 1.0_rp/(xmax - xmin)
383 rdomy = 1.0_rp/(ymax - ymin)
384
385 do ke=1, ne
386 do f=1, nfaces
387 x = sum(pos_en(fmask(:,f),ke,1)/dble(nfp))
388 y = sum(pos_en(fmask(:,f),ke,2)/dble(nfp))
389
390 call eval_domain_boundary(1, y, ymin, x, ke, f, rdomy)
391 call eval_domain_boundary(2, x, xmax, y, ke, f, rdomx)
392 call eval_domain_boundary(3, y, ymax, x, ke, f, rdomy)
393 call eval_domain_boundary(4, x, xmin, y, ke, f, rdomx)
394 end do
395 end do
396
397 allocate( mapb(sum(counterb*nfp)) )
398 allocate( vmapb(size(mapb)) )
399
400 mapb_counter = 1
401 do b = 1, 4
402 ! write(*,*) "LocalMesh boundary ID:", b
403 ! write(*,*) counterB(b)
404 ! write(*,*) ordInfo(1:counterB(b),b)
405 ! write(*,*) elemIds(1:counterB(b),b)
406 ! write(*,*) faceIds(1:counterB(b),b)
407
408 do i=1, counterb(b)
409 ke = elemids(i,b); f = faceids(i,b)
410 do j=1, nfp
411 vmapp(j,f,ke) = np*ne + mapb_counter
412 vmapb(mapb_counter) = fmask(j,f) + (ke-1)*np
413 mapb_counter = mapb_counter + 1
414 end do
415 end do
416 end do
417
418 ! write(*,*) "VMapP:-----"
419 ! do b=1, 4
420 ! do i=1, counterB(b)
421 ! k = elemIds(i,b); f = faceIDs(i,b)
422 ! write(*,*) "bid=", b, ":", mesh%VmapP(:,f,k)
423 ! end do
424 ! end do
425 ! write(*,*) "-----"
426 ! write(*,*) "VMapB:", mesh%VmapB(:)
427 ! write(*,*) "NeA=", mesh%NeA
428
429 return
430 contains
431 subroutine eval_domain_boundary(domb_id, r, rbc, ord_info, ke_, f_, normalized_fac)
432 integer, intent(in) :: domb_id
433 real(rp), intent(in) :: r
434 real(rp), intent(in) :: rbc
435 real(rp), intent(in) :: ord_info
436 integer, intent(in) :: ke_, f_
437 real(rp), intent(in) :: normalized_fac
438 !-------------------------------------------------------------
439
440 if ( abs(r - rbc)*normalized_fac < node_tol ) then
441 counterb(domb_id) = counterb(domb_id) + 1
442 ordinfo(counterb(domb_id),domb_id) = ord_info
443 elemids(counterb(domb_id),domb_id) = ke_
444 faceids(counterb(domb_id),domb_id) = f_
445 end if
446
447 return
448 end subroutine eval_domain_boundary
449
450 end subroutine meshutil2d_genpatchboundarymap
451
452!OCL SERIAL
453 subroutine meshutil2d_buildglobalmap( &
454 panelID_table, pi_table, pj_table, &
455 tileID_map, tileFaceID_map, tilePanelID_map, &
456 Ntile, isPeriodicX, isPeriodicY, &
457 Ne_x, Ne_y )
458
459 implicit none
460
461 integer, intent(in) :: ntile
462 integer, intent(out) :: panelid_table(ntile)
463 integer, intent(out) :: pi_table(ntile)
464 integer, intent(out) :: pj_table(ntile)
465 integer, intent(out) :: tileid_map(4,ntile)
466 integer, intent(out) :: tilefaceid_map(4,ntile)
467 integer, intent(out) :: tilepanelid_map(4,ntile)
468 logical, intent(in) :: isperiodicx
469 logical, intent(in) :: isperiodicy
470 integer, intent(in) :: ne_x
471 integer, intent(in) :: ne_y
472
473 integer :: ntileperpanel
474 integer :: nvx, nvy
475 integer, allocatable :: nodesid_2d(:,:)
476 integer, allocatable :: etov(:,:)
477 integer, allocatable :: etoe(:,:)
478 integer, allocatable :: etof(:,:)
479 integer :: i, j, f
480 integer :: panelid
481 integer :: tileid, tileid_r
482 integer :: counter
483
484 !-----------------------------------------------------------------------------
485
486 ntileperpanel = ntile/1
487 nvx = ne_x + 1
488 nvy = ne_y + 1
489 allocate( nodesid_2d(nvx, nvy) )
490 allocate( etov(ntile,4), etoe(ntile,4), etof(ntile,4) )
491
492 counter = 0
493 do j = 1, nvy
494 do i = 1, nvx
495 counter = counter + 1
496 nodesid_2d(i,j) = counter
497 end do
498 end do
499
500 !----
501
502 tileid = 0
503 do j = 1, ne_y
504 do i = 1, ne_x
505 tileid = tileid + 1
506 panelid_table(tileid) = 1
507 pi_table(tileid) = i; pj_table(tileid) = j
508 etov(tileid,:) = (/ nodesid_2d(i,j ), nodesid_2d(i+1,j ), &
509 & nodesid_2d(i,j+1), nodesid_2d(i+1,j+1)/)
510 end do
511 end do
512
513 call meshutil2d_genconnectivity( etoe, etof, &
514 & etov, ntile, 4 )
515 tileid_map(:,:) = transpose(etoe)
516 tilefaceid_map(:,:) = transpose(etof)
517
518 do tileid=1, ntile
519 do f=1, 4
520 tileid_r = tileid_map(f,tileid)
521 tilepanelid_map(f,tileid) = panelid_table(tileid_r)
522 end do
523 end do
524
525 if (isperiodicx) then
526 do tileid=1, ntile
527 if (pi_table(tileid) == 1 .and. tilefaceid_map(4,tileid) == 4) then
528 tileid_map(4,tileid) = ne_x + (pj_table(tileid) - 1)*ne_x
529 tilefaceid_map(4,tileid) = 2
530 end if
531 if (pi_table(tileid) == ne_x .and. tilefaceid_map(2,tileid) == 2) then
532 tileid_map(2,tileid) = 1 + (pj_table(tileid) - 1)*ne_x
533 tilefaceid_map(2,tileid) = 4
534 end if
535 end do
536 end if
537
538 if (isperiodicy) then
539 do tileid=1, ntile
540 if (pj_table(tileid) == 1 .and. tilefaceid_map(1,tileid) == 1) then
541 tileid_map(1,tileid) = pi_table(tileid) + (ne_y - 1)*ne_x
542 tilefaceid_map(1,tileid) = 3
543 end if
544 if (pj_table(tileid) == ne_y .and. tilefaceid_map(3,tileid) == 3) then
545 tileid_map(3,tileid) = pi_table(tileid)
546 tilefaceid_map(3,tileid) = 1
547 end if
548 end do
549 end if
550
551 return
552
553 !--
554
555 ! if (PRC_isMaster) then
556 ! write(*,*) "TotTile", Ntile
557 ! do tileID = 1, Ntile
558 ! write(*,*) "tileID:", tileID, ", EtoV:", EtoV(tileID,:)
559 ! end do
560 ! write(*,*) "-----------"
561 ! do tileID = 1, Ntile
562 ! write(*,*) "tileID:", tileID, ", EtoE:", EtoE(tileID,:)
563 ! end do
564 ! write(*,*) "-----------"
565 ! do tileID = 1, Ntile
566 ! write(*,*) "tileID:", tileID, ", EtoF:", EtoF(tileID,:)
567 ! end do
568 ! end if
569 end subroutine meshutil2d_buildglobalmap
570
571end module scale_meshutil_2d
scale_quicksort::quicksort_exec_with_idx
Definition scale_quicksort.F90:31
scale_meshutil_2d
module FElib / Mesh / utility for 2D mesh
Definition scale_meshutil_2d.F90:11
scale_meshutil_2d::meshutil2d_buildinteriormap
subroutine, public meshutil2d_buildinteriormap(vmapm, vmapp, mapm, mapp, pos_en, pos_ev, etoe, etof, etov, fmask, ne, np, nfp, nfaces, nv)
Definition scale_meshutil_2d.F90:213
scale_meshutil_2d::meshutil2d_genconnectivity
subroutine, public meshutil2d_genconnectivity(etoe, etof, etov, ne, nfaces)
Definition scale_meshutil_2d.F90:96
scale_meshutil_2d::meshutil2d_genrectdomain
subroutine, public meshutil2d_genrectdomain(pos_v, etov, ke_x, xmin, xmax, ke_y, ymin, ymax)
Definition scale_meshutil_2d.F90:36
scale_meshutil_2d::meshutil2d_buildglobalmap
subroutine, public meshutil2d_buildglobalmap(panelid_table, pi_table, pj_table, tileid_map, tilefaceid_map, tilepanelid_map, ntile, isperiodicx, isperiodicy, ne_x, ne_y)
Definition scale_meshutil_2d.F90:458
scale_meshutil_2d::meshutil2d_genpatchboundarymap
subroutine, public meshutil2d_genpatchboundarymap(vmapb, mapb, vmapp, pos_en, xmin, xmax, ymin, ymax, fmask, ne, np, nfp, nfaces, nv)
Definition scale_meshutil_2d.F90:345
scale_quicksort
module common / sort algorithm
Definition scale_quicksort.F90:15
Generated by doxygen 1.10.0