-
Notifications
You must be signed in to change notification settings - Fork 38
/
connectivity.f90
3970 lines (3970 loc) · 135 KB
/
connectivity.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
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
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
module connectivity_mod
!
!.. Global Use Statements ..
use module_kind_types
!
implicit none
!
private
!
!
! ##################################
! ### Module Public Procedures ###
! ##################################
!
public :: get_grid_connectivity
public :: find_cells_surrounding_each_node
public :: find_cells_surrounding_each_cell
public :: find_nodes_on_each_edge
public :: find_cells_surrounding_each_edge
public :: create_face_derived_type
!
public :: get_unit_cell_subconnectivity
!
public :: GetFaceNodes_for_QuadFace
public :: GetFaceNodes_for_TriaFace
public :: RotateFluxPoints_Quad
public :: RotateFluxPoints_Tria
!
!
! #########################################
! ##### MODULE GENERIC INTERFACES #####
! #########################################
!
interface RotateFluxPoints_Edge
module procedure RotateFluxPoints_Edge, &
FluxPointRotation_Edge
end interface RotateFluxPoints_Edge
!
interface RotateFluxPoints_Quad
module procedure RotateFluxPoints_Quad, &
FluxPointRotation_Quad
end interface RotateFluxPoints_Quad
!
interface RotateFluxPoints_Tria
module procedure RotateFluxPoints_Tria, &
FluxPointRotation_Tria
end interface RotateFluxPoints_Tria
!
contains
!
!###############################################################################
!
subroutine get_grid_connectivity(bface,cell_geom,cell_order,xyz_nodes, &
nodes_of_cell,nodes_of_cell_ptr)
!
!.. Use Statements ..
use ovar, only : continuous_output
use geovar, only : nface,face,fp
!
!.. Formal Arguments ..
integer, dimension(:,:), intent(in) :: bface
integer, dimension(:), intent(in) :: cell_geom
integer, dimension(:), intent(in) :: cell_order
integer, dimension(:), intent(in) :: nodes_of_cell
integer, dimension(:), intent(in) :: nodes_of_cell_ptr
real(wp), dimension(:,:), intent(in) :: xyz_nodes
!
!.. Local Scalars ..
integer :: ierr
!
!.. Local Allocatable Arrays ..
integer, allocatable, dimension(:) :: cells_with_node
integer, allocatable, dimension(:) :: cells_with_node_ptr
integer, allocatable, dimension(:) :: cells_surr_cell
integer, allocatable, dimension(:) :: cells_surr_cell_ptr
integer, allocatable, dimension(:) :: cells_with_edge
integer, allocatable, dimension(:) :: cells_with_edge_ptr
integer, allocatable, dimension(:) :: edges_with_node
integer, allocatable, dimension(:) :: edges_with_node_ptr
integer, allocatable, dimension(:,:) :: nodes_on_edge
!
!.. Local Parameters ..
character(len=*), parameter :: pname = "get_grid_connectivity"
!
continue
!
call debug_timer(entering_procedure,pname)
!
! Get the cells surrounding each node
!
call memory_pause("entering grid connectivity master subroutine", &
"finding cells surrounding each node")
call find_cells_surrounding_each_node(bface, &
nodes_of_cell,nodes_of_cell_ptr, &
cells_with_node,cells_with_node_ptr)
!
! Get the cells surrounding each cell
!
call memory_pause("finding cells surrounding each node", &
"finding the cells surrounding each cell")
call find_cells_surrounding_each_cell(bface,cell_geom, &
nodes_of_cell,nodes_of_cell_ptr, &
cells_with_node,cells_with_node_ptr, &
cells_surr_cell,cells_surr_cell_ptr)
!
! If using the continuous output option, find the
! node and cell connectivity for all the grid edges
!
if (continuous_output) then
!
! Get the nodes defining each edge as well as the edges
! surrounding each node
!
call find_nodes_on_each_edge(cell_geom, &
nodes_of_cell,nodes_of_cell_ptr, &
cells_with_node,cells_with_node_ptr, &
nodes_on_edge)
!nodes_on_edge, &
!edges_with_node,edges_with_node_ptr)
!
! Get the cells that surround each edge
!
call find_cells_surrounding_each_edge(nodes_on_edge, &
nodes_of_cell,nodes_of_cell_ptr, &
cells_with_node,cells_with_node_ptr, &
cells_with_edge,cells_with_edge_ptr)
!
end if
!
! Get the cells on each face
!
call memory_pause("finding the cells surrounding each cell", &
"finding the cells on each face")
call create_face_derived_type(bface,cell_geom,cell_order,xyz_nodes, &
nodes_of_cell,nodes_of_cell_ptr, &
cells_surr_cell,cells_surr_cell_ptr, &
face,fp)
!
! Save the total number of faces to nface
!
nface = size(face)
!
! Get the indices for the flux points for each cell
!
call memory_pause("finding the cells on each face", &
"finding the flux points of each cell")
call create_cell_derived_type(nodes_of_cell,nodes_of_cell_ptr, &
cells_with_node_ptr,nodes_on_edge, &
cells_with_edge,cells_with_edge_ptr)
!
! If using the continuous output option, create the edge derived type array
!
if (continuous_output) then
call create_edge_derived_type(cells_with_edge,cells_with_edge_ptr, &
nodes_on_edge)
end if
!
! Get the unit cell subconnectivities for use when writing solution file
!
call memory_pause("finding the flux points of each cell", &
"getting the unit cell subconnectivities")
call get_unit_cell_subconnectivity
call memory_pause("getting the unit cell subconnectivities", &
"leaving the grid connectivity master subroutine")
!
! Deallocate the local arrays before leaving
!
if (allocated(cells_with_node)) then
deallocate ( cells_with_node , stat=ierr , errmsg=error_message )
call alloc_error(pname,"cells_with_node",2,__LINE__,__FILE__,ierr, &
error_message)
end if
!
if (allocated(cells_with_node_ptr)) then
deallocate ( cells_with_node_ptr , stat=ierr , errmsg=error_message )
call alloc_error(pname,"cells_with_node_ptr",2,__LINE__,__FILE__,ierr, &
error_message)
end if
!
if (allocated(cells_surr_cell)) then
deallocate ( cells_surr_cell , stat=ierr , errmsg=error_message )
call alloc_error(pname,"cells_surr_cell",2,__LINE__,__FILE__,ierr, &
error_message)
end if
!
if (allocated(cells_surr_cell_ptr)) then
deallocate ( cells_surr_cell_ptr , stat=ierr , errmsg=error_message )
call alloc_error(pname,"cells_surr_cell_ptr",2,__LINE__,__FILE__,ierr, &
error_message)
end if
!
if (allocated(edges_with_node)) then
deallocate ( edges_with_node , stat=ierr , errmsg=error_message )
call alloc_error(pname,"edges_with_node",2,__LINE__,__FILE__,ierr, &
error_message)
end if
!
if (allocated(edges_with_node_ptr)) then
deallocate ( edges_with_node_ptr , stat=ierr , errmsg=error_message )
call alloc_error(pname,"edges_with_node_ptr",2,__LINE__,__FILE__,ierr, &
error_message)
end if
!
if (allocated(cells_with_edge)) then
deallocate ( cells_with_edge , stat=ierr , errmsg=error_message )
call alloc_error(pname,"cells_with_edge",2,__LINE__,__FILE__,ierr, &
error_message)
end if
!
if (allocated(cells_with_edge_ptr)) then
deallocate ( cells_with_edge_ptr , stat=ierr , errmsg=error_message )
call alloc_error(pname,"cells_with_edge_ptr",2,__LINE__,__FILE__,ierr, &
error_message)
end if
!
if (allocated(nodes_on_edge)) then
deallocate ( nodes_on_edge , stat=ierr , errmsg=error_message )
call alloc_error(pname,"nodes_on_edge",2,__LINE__,__FILE__,ierr, &
error_message)
end if
!
call debug_timer(leaving_procedure,pname)
!
end subroutine get_grid_connectivity
!
!###############################################################################
!
subroutine find_cells_surrounding_each_node(bface, &
nodes_of_cell,nodes_of_cell_ptr, &
cells_with_node,cells_with_node_ptr, &
include_bnd_faces, &
make_bnd_faces_negative)
!
! Find the cells that surround each node
!
! For node i_node,
! cells_with_node(ibeg:iend) gives all of the grid cells that contain i_node
! where : ibeg = cells_with_node_ptr(i_node)+1
! iend = cells_with_node_ptr(i_node+1)
!
!.. Formal Arguments ..
integer, dimension(:,:), intent(in) :: bface
integer, dimension(:), intent(in) :: nodes_of_cell
integer, dimension(:), intent(in) :: nodes_of_cell_ptr
integer, allocatable, dimension(:), intent(inout) :: cells_with_node
integer, allocatable, dimension(:), intent(inout) :: cells_with_node_ptr
!
!.. Optional Arguments ..
logical(lk), optional, intent(in) :: include_bnd_faces
logical(lk), optional, intent(in) :: make_bnd_faces_negative
!
!.. Local Scalars ..
integer :: i,ip,is,ibeg,iend,n,nc,ierr
integer :: n1,n2,np,lnode,lfbnd,lcell
integer :: bnd_face_coef
logical(lk) :: include_bnd_face_ghost_cells
!
!.. Local Allocatable Arrays ..
integer, allocatable, dimension(:) :: cwn
!
!.. Local Parameters ..
character(len=*), parameter :: pname = "find_cells_surrounding_each_node"
!
continue
!
call debug_timer(entering_procedure,pname)
!
! Get the grid size from the input arrays
!
lnode = maxval(nodes_of_cell)
lfbnd = size(bface,dim=2)
lcell = size(nodes_of_cell_ptr)-1 - lfbnd
!
! Determine whether to include the boundary face ghost cells
! when finding the cells surrounding each node
!
include_bnd_face_ghost_cells = true
if (present(include_bnd_faces)) then
include_bnd_face_ghost_cells = include_bnd_faces
end if
!
! Determine the coefficient for the boundary face ghost cell indices
! based on the optional argument make_bnd_faces_negative
!
bnd_face_coef = 1
if (present(make_bnd_faces_negative)) then
bnd_face_coef = merge(-1,1,make_bnd_faces_negative)
end if
!
! Allocate cells_with_node_ptr
!
allocate ( cells_with_node_ptr(1:lnode+1) , source=0 , &
stat=ierr , errmsg=error_message )
call alloc_error(pname,"cells_with_node_ptr",1,__LINE__,__FILE__,ierr, &
error_message)
!
! Loop over the cells, storing 'ahead'
!
do n = 1,lcell
!
ibeg = nodes_of_cell_ptr(n)+1
iend = nodes_of_cell_ptr(n+1)
!
do i = ibeg,iend
ip = nodes_of_cell(i) + 1
cells_with_node_ptr(ip) = cells_with_node_ptr(ip) + 1
end do
!
end do
!
! Add in the nodes of the ghost faces
!
if (include_bnd_face_ghost_cells) then
!
do n = 1,lfbnd
!
nc = bface(10,n)
!
ibeg = nodes_of_cell_ptr(nc)+1
iend = nodes_of_cell_ptr(nc+1)
!
do i = ibeg,iend
ip = nodes_of_cell(i)
cells_with_node_ptr(ip+1) = cells_with_node_ptr(ip+1) + 1
end do
!
end do
!
end if
!
! Reshuffle cells_with_node_ptr
!
do n = 2,size(cells_with_node_ptr)
cells_with_node_ptr(n) = cells_with_node_ptr(n) + cells_with_node_ptr(n-1)
end do
!
allocate ( cells_with_node(1:last(cells_with_node_ptr)) , source=0 , &
stat=ierr , errmsg=error_message )
call alloc_error(pname,"cells_with_node",1,__LINE__,__FILE__,ierr, &
error_message)
!
! Now store the cells surrounding each node in cells_with_node
!
do n = 1,lcell
!
ibeg = nodes_of_cell_ptr(n)+1
iend = nodes_of_cell_ptr(n+1)
!
do i = ibeg,iend
ip = nodes_of_cell(i)
is = cells_with_node_ptr(ip) + 1
cells_with_node_ptr(ip) = is
cells_with_node(is) = n
end do
!
end do
!
! Add in the nodes of the ghost faces
!
if (include_bnd_face_ghost_cells) then
!
do n = 1,lfbnd
!
nc = bface(10,n)
!
ibeg = nodes_of_cell_ptr(nc)+1
iend = nodes_of_cell_ptr(nc+1)
!
do i = ibeg,iend
ip = nodes_of_cell(i)
is = cells_with_node_ptr(ip) + 1
cells_with_node_ptr(ip) = is
!
! Make the index for this boundary face ghost cell negative
! if the optional argument make_bnd_faces_negative is present
! with the value .TRUE.
!
cells_with_node(is) = bnd_face_coef * nc
!
end do
!
end do
!
end if
!
! Finally, reorder cells_with_node_ptr
!
do n = size(cells_with_node_ptr),2,-1
cells_with_node_ptr(n) = cells_with_node_ptr(n-1)
end do
cells_with_node_ptr(1) = 0
!
! Perform one last error check to make sure that no cells have been added
! more than once to a given node. This was known to happen when the ghost
! cells were treated as actual cells but shouldnt happen anymore now that
! ghost cells are actually ghost faces. I have removed the old code that
! prevented this from happening in the above loops since it required a mask
! array that had the potential to consume a lot of memory. This error check
! has been added JUST-IN-CASE I have overlooked a scenario in which this
! could possible happen.
!
nc = 0
do n = 1,lnode
nc = max( nc , cells_with_node_ptr(n+1) - cells_with_node_ptr(n) )
end do
!
allocate ( cwn(1:nc) , source=0 , stat=ierr , errmsg=error_message )
call alloc_error(pname,"cwn",1,__LINE__,__FILE__,ierr,error_message, &
skip_alloc_pause)
!
do n = 1,lnode
!
n1 = cells_with_node_ptr(n)+1
n2 = cells_with_node_ptr(n+1)
np = n2-n1+1
!
cwn(1:np) = cells_with_node(n1:n2)
call make_unique(cwn(1:np),ip)
if (ip /= np) then
!
write (iout,9023) n,(cells_with_node(i),i=n1,n2)
9023 format ("Found a duplicate cell around node ",i0,/, &
" Cell list: ",100(2x,i0))
!" Cell list: ",*(2x,i0))
!
! Compute the difference from the original number of cells found for
! this node and the new unique number of cells that surround this node
nc = np - ip
! Adjust the pointer index for the remaining nodes
cells_with_node_ptr(n+1:) = cells_with_node_ptr(n+1:) - nc
! Get the new value for n2
n2 = n1 + ip - 1
! Copy the unique cell list back into cells_with_node
cells_with_node(n1:n2) = cwn(1:ip)
! Perform end-off-shifts to move the remaining parts of cells_with_node
! to the end of the list of cells for the current node
cells_with_node(n2+1:) = eoshift(cells_with_node(n2+1:),shift=nc)
!
end if
!
end do
!
if (allocated(cwn)) then
deallocate ( cwn , stat=ierr , errmsg=error_message )
call alloc_error(pname,"cwn",2,__LINE__,__FILE__,ierr,error_message, &
skip_alloc_pause)
end if
!
! Reallocate cells_with_node if it was changed
!
if (last(cells_with_node_ptr) /= size(cells_with_node)) then
!
n1 = last(cells_with_node_ptr) + 1
n2 = ubound(cells_with_node,dim=1)
!
if (any(cells_with_node(n1:n2) /= 0)) then
write (error_message,10)
call stop_gfr(abort,pname,__LINE__,__FILE__,error_message)
else
call reallocate(cells_with_node,last(cells_with_node_ptr))
end if
!
end if
!
call debug_timer(leaving_procedure,pname)
!
! Format Statements
!
10 format ("Error: The end of cells_with_node should be zero!")
!
contains
!
pure subroutine make_unique(cell_list,nc)
!.. Formal Arguments ..
integer, intent(out) :: nc
integer, dimension(:), intent(inout) :: cell_list
!.. Local Scalars ..
integer :: i
!.. Local Arrays ..
integer, dimension(1:size(cell_list)) :: array_of_zeros
continue
!
array_of_zeros = 0
!
nc = size(cell_list)
do i = 1,nc-1
where (cell_list(i+1:nc) == cell_list(i)) cell_list(i+1:nc) = 0
end do
!
nc = count( cell_list /= 0 )
cell_list = pack( cell_list , cell_list /= 0 , array_of_zeros )
!
end subroutine make_unique
!
end subroutine find_cells_surrounding_each_node
!
!###############################################################################
!
subroutine find_cells_surrounding_each_cell(bface,cell_geom, &
nodes_of_cell,nodes_of_cell_ptr, &
cells_with_node,cells_with_node_ptr, &
cells_surr_cell,cells_surr_cell_ptr)
!
! Find the cells that surround each cell
!
! For cell i_cell, ...
! cells_surr_cell(ibeg:iend) gives all of the grid cells
! that are adjacent to i_cell
! where : ibeg = cells_surr_cell_ptr(i_cell)+1
! iend = cells_surr_cell_ptr(i_cell+1)
!
!.. Formal Arguments ..
integer, dimension(:,:), intent(in) :: bface
integer, dimension(:), intent(in) :: cell_geom
integer, dimension(:), intent(in) :: nodes_of_cell
integer, dimension(:), intent(in) :: nodes_of_cell_ptr
integer, dimension(:), intent(in) :: cells_with_node
integer, dimension(:), intent(in) :: cells_with_node_ptr
integer, allocatable, dimension(:), intent(inout) :: cells_surr_cell
integer, allocatable, dimension(:), intent(inout) :: cells_surr_cell_ptr
!
!.. Local Scalars ..
integer :: i,i1,i2,ieadj,csc_ptr
integer :: j,j1,j2,n,n1,n2,ierr
integer :: this_geom,this_side,nfnods
integer :: lnode,lcell
!
!.. Local Arrays ..
integer, dimension(1:4) :: ip
integer, dimension(1:4) :: face_nodes
!
!.. Local Allocatable Arrays ..
integer, allocatable, dimension(:) :: lpoin
!
!.. Local Parameters ..
character(len=*), parameter :: pname = "find_cells_surrounding_each_cell"
!
continue
!
call debug_timer(entering_procedure,pname)
!
lnode = maxval(nodes_of_cell)
lcell = size(nodes_of_cell_ptr)-1 - size(bface,dim=2)
!
allocate ( lpoin(1:lnode) , source=0 , stat=ierr , errmsg=error_message )
call alloc_error(pname,"lpoin",1,__LINE__,__FILE__,ierr,error_message)
!
allocate ( cells_surr_cell_ptr(1:lcell+1) , source=0 , &
stat=ierr , errmsg=error_message )
call alloc_error(pname,"cells_surr_cell_ptr",1,__LINE__,__FILE__,ierr, &
error_message)
!
! Loop to set up cells_surr_cell_ptr
!
do n = 1,lcell
cells_surr_cell_ptr(n+1) = cells_surr_cell_ptr(n) + geom_faces(cell_geom(n))
end do
!
allocate ( cells_surr_cell(1:last(cells_surr_cell_ptr)) , source=0 , &
stat=ierr , errmsg=error_message )
call alloc_error(pname,"cells_surr_cell",1,__LINE__,__FILE__,ierr, &
error_message)
!
! Loop to find cells_surr_cell
!
do n = 1,lcell
!
! Get the surrounding cells based on the geometry of the current cell
!
n1 = nodes_of_cell_ptr(n)+1
n2 = nodes_of_cell_ptr(n+1)
!
! Get the pointer index of cells_surr_cell for this cell
!
csc_ptr = cells_surr_cell_ptr(n)
!
! Geometry of the current cell
!
this_geom = cell_geom(n)
!
! Loop over the number of cell faces for the current geometry
!
do this_side = 1,geom_faces(this_geom)
!
! Find the number of nodes that need to be matched for this face
!
nfnods = num_face_nodes(this_side,this_geom)
!
! Node offsets for this face of the current cell
!
ip(1:nfnods) = side_nodes(1:nfnods,this_side,this_geom)
!
! Store nodes of the current face for the current cell
!
face_nodes(1:nfnods) = nodes_of_cell( n1 + ip(1:nfnods) )
!
! Mark each of these nodes on the current face
!
lpoin( face_nodes(1:nfnods) ) = 1
!
! Get the beginning and ending indices for the
! cells surrounding the first node on the face
!
i1 = cells_with_node_ptr( face_nodes(1) )+1
i2 = cells_with_node_ptr( face_nodes(1) +1)
!
! Initialize the adjacent face marker to zero
!
ieadj = 0
!
! Loop over the cells that contain the first node on the face
!
find_cell: do i = i1,i2
!
! Current cell containing this node
! NOTE: We need the absolute value of this because the cell index could
! be negative if it is for a boundary face ghost cell and the
! optional argument make_bnd_faces_negative was present with the
! value .TRUE. when subroutine find_cells_surrounding_each_node
! was called
!
j = abs(cells_with_node(i))
!
! Skip to the next cell if this is the current cell we are working on
!
if (j == n) cycle find_cell
!
! Get the beginning and ending indices of the nodes defining the
! possible adjacent cell
!
j1 = nodes_of_cell_ptr(j)+1
j2 = nodes_of_cell_ptr(j+1)
!
! If summing the lpoin array using the current nodes of the possible
! adjacent cell is equal to the number of nodes we need to match, then
! this is the adjacent cell to the current face.
! NOTE: We need to use the actual value of cells_with_node(i) instead
! of j because we want to keep the cell index negative if this
! is a boundary face ghost cell
!
if (sum(lpoin(nodes_of_cell(j1:j2))) == nfnods) then
ieadj = cells_with_node(i)
exit find_cell
end if
!
end do find_cell
!
! Store the adjacent cell index
!
cells_surr_cell(csc_ptr+this_side) = ieadj
!
! Re-mark the nodes in lpoin to zero
!
lpoin( face_nodes(1:nfnods) ) = 0
!
end do
!
end do
!
deallocate ( lpoin , stat=ierr , errmsg=error_message )
call alloc_error(pname,"lpoin",2,__LINE__,__FILE__,ierr,error_message)
!
call debug_timer(leaving_procedure,pname)
!
end subroutine find_cells_surrounding_each_cell
!
!###############################################################################
!
subroutine initialize_flxpts(fp)
!
!.. Use Statements ..
use order_mod, only : geom_solpts
use order_mod, only : n_min_order,n_max_order
use geovar, only : fp_t
!
!.. Formal Arguments ..
type(fp_t), allocatable, intent(inout) :: fp
!
!.. Local Scalars ..
integer :: gmin,gmax,omin,omax,ierr
integer :: n,nep,nfp,this_order
!
character(len=100) :: array_name
!
!.. Local Allocatable Arrays ..
integer, allocatable :: fp_idx(:,:)
!
!.. Local Parameters ..
character(len=*), parameter :: pname = "initialize_flxpts"
!
continue
!
call debug_timer(entering_procedure,pname)
!
gmin = min(Geom_Edge,Geom_Tria,Geom_Quad)
gmax = max(Geom_Edge,Geom_Tria,Geom_Quad)
!
omin = n_min_order
omax = n_max_order
!
allocate ( fp , stat=ierr , errmsg=error_message )
call alloc_error(pname,"fp",1,__LINE__,__FILE__,ierr,error_message)
!
allocate ( fp%geord(gmin:gmax,omin:omax) , stat=ierr , errmsg=error_message )
call alloc_error(pname,"fp%geord",1,__LINE__,__FILE__,ierr,error_message)
!
do this_order = omin,omax
!
! Edge face
!
nep = geom_solpts(Geom_Edge,this_order)
!
allocate ( fp%geord(Geom_Edge,this_order)%rot(0:2) , &
stat=ierr , errmsg=error_message )
write (array_name,1) "Geom_Edge",this_order
call alloc_error(pname,array_name,1,__LINE__,__FILE__,ierr,error_message)
!
allocate ( fp%geord(Geom_Edge,this_order)%rot(0)%pts(1:nep) , source=0 , &
stat=ierr , errmsg=error_message )
write (array_name,1) "Geom_Edge",this_order,1
call alloc_error(pname,array_name,1,__LINE__,__FILE__,ierr,error_message)
!
allocate ( fp%geord(Geom_Edge,this_order)%rot(1)%pts(1:nep) , source=0 , &
stat=ierr , errmsg=error_message )
write (array_name,1) "Geom_Edge",this_order,1
call alloc_error(pname,array_name,1,__LINE__,__FILE__,ierr,error_message)
!
allocate ( fp%geord(Geom_Edge,this_order)%rot(2)%pts(1:nep) , source=0 , &
stat=ierr , errmsg=error_message )
write (array_name,1) "Geom_Edge",this_order,2
call alloc_error(pname,array_name,1,__LINE__,__FILE__,ierr,error_message)
!
fp%geord(Geom_Edge,this_order)%rot(0)%pts(1:nep) = 0
fp%geord(Geom_Edge,this_order)%rot(1)%pts(1:nep) = intseq(0,nep-1)
fp%geord(Geom_Edge,this_order)%rot(2)%pts(1:nep) = intseq(nep-1,0,-1)
!
! Triangular Face
!
nfp = geom_solpts(Geom_Tria,this_order)
!
allocate ( fp%geord(Geom_Tria,this_order)%rot(0:6) , &
stat=ierr , errmsg=error_message )
write (array_name,1) "Geom_Tria",this_order
call alloc_error(pname,array_name,1,__LINE__,__FILE__,ierr,error_message)
!
do n = lbound(fp%geord(Geom_Tria,this_order)%rot,dim=1), &
ubound(fp%geord(Geom_Tria,this_order)%rot,dim=1)
!
allocate ( fp%geord(Geom_Tria,this_order)%rot(n)%pts(1:nfp) , source=0 , &
stat=ierr , errmsg=error_message )
write (array_name,1) "Geom_Tria",this_order,n
call alloc_error(pname,array_name,1,__LINE__,__FILE__,ierr,error_message)
!
end do
!
! Quadrilateral Face
!
nfp = geom_solpts(Geom_Quad,this_order)
!
allocate ( fp%geord(Geom_Quad,this_order)%rot(0:8) , &
stat=ierr , errmsg=error_message )
write (array_name,1) "Geom_Quad",this_order
call alloc_error(pname,array_name,1,__LINE__,__FILE__,ierr,error_message)
!
do n = lbound(fp%geord(Geom_Quad,this_order)%rot,dim=1), &
ubound(fp%geord(Geom_Quad,this_order)%rot,dim=1)
!
allocate ( fp%geord(Geom_Quad,this_order)%rot(n)%pts(1:nfp) , source=0 , &
stat=ierr , errmsg=error_message )
write (array_name,1) "Geom_Quad",this_order,n
call alloc_error(pname,array_name,1,__LINE__,__FILE__,ierr,error_message)
!
end do
!
! See the function RotateFluxPoints_Quad for a complete description of
! how these rotations are defined
!
if (allocated(fp_idx)) then
deallocate ( fp_idx , stat=ierr , errmsg=error_message )
call alloc_error(pname,"fp_idx",2,__LINE__,__FILE__,ierr,error_message)
end if
!
! Create fp_idx as a square 2D array containing
! the consecutive indices from 0 to nfp-1
!
allocate ( fp_idx(1:nep,1:nep) , &
source=reshape( intseq(0,nfp-1) , [nep,nep] ) , &
stat=ierr , errmsg=error_message )
call alloc_error(pname,"fp_idx",1,__LINE__,__FILE__,ierr,error_message)
!
! Set all the flux point indices for the error rotation to 0
!
fp%geord(Geom_Quad,this_order)%rot(0)%pts(1:nfp) = 0
!
! Create the flux point indices for the first four rotations where the
! "a" and "b" vectors are aligned in the same coordinate directions
!
fp%geord(Geom_Quad,this_order)%rot(1)%pts(1:nfp) = reshape( fp_idx , [nfp] )
!
fp%geord(Geom_Quad,this_order)%rot(2)%pts(1:nfp) = &
reshape( fp_idx(:,nep:1:-1) , [nfp] )
!
fp%geord(Geom_Quad,this_order)%rot(3)%pts(1:nfp) = &
reshape( fp_idx(nep:1:-1,:) , [nfp] )
!
fp%geord(Geom_Quad,this_order)%rot(4)%pts(1:nfp) = &
reshape( fp_idx(nep:1:-1,nep:1:-1) , [nfp] )
!
! Create the flux point indices for the last four rotations where the
! "a" and "b" vectors are aligned in opposite coordinate directions,
! which requires that we first transpose fp_idx
!
fp_idx = transpose(fp_idx)
!
fp%geord(Geom_Quad,this_order)%rot(5)%pts(1:nfp) = reshape( fp_idx , [nfp] )
!
fp%geord(Geom_Quad,this_order)%rot(6)%pts(1:nfp) = &
reshape( fp_idx(:,nep:1:-1) , [nfp] )
!
fp%geord(Geom_Quad,this_order)%rot(7)%pts(1:nfp) = &
reshape( fp_idx(nep:1:-1,:) , [nfp] )
!
fp%geord(Geom_Quad,this_order)%rot(8)%pts(1:nfp) = &
reshape( fp_idx(nep:1:-1,nep:1:-1) , [nfp] )
!
end do
!
if (allocated(fp_idx)) then
deallocate ( fp_idx , stat=ierr , errmsg=error_message )
call alloc_error(pname,"fp_idx",2,__LINE__,__FILE__,ierr,error_message)
end if
!
call debug_timer(leaving_procedure,pname)
!
! Format Statements
!
1 format ("fp%geord(",a,",",i0,")%rot",:,"%(",i0,")%pts")
!
end subroutine initialize_flxpts
!
!###############################################################################
!
subroutine create_face_derived_type(bface,cell_geom,cell_order,xyz_nodes, &
nodes_of_cell,nodes_of_cell_ptr, &
cells_surr_cell,cells_surr_cell_ptr, &
face,fp,force_flxpt_offset_to_zero)
!
! Procedure that creates the face array.
! Look in the module geovar for the definition of the derived type.
!
!.. Use Statements ..
use order_mod, only : maxFP
use order_mod, only : geom_solpts
use order_mod, only : find_geom_max_interp_order
use geovar, only : face_t,fp_t
!
!.. Formal Arguments ..
integer, intent(in) :: bface(:,:)
integer, intent(in) :: cell_geom(:)
integer, intent(in) :: cell_order(:)
integer, intent(in) :: nodes_of_cell(:)
integer, intent(in) :: nodes_of_cell_ptr(:)
real(wp), intent(in) :: xyz_nodes(:,:)
integer, intent(in) :: cells_surr_cell(:)
integer, intent(in) :: cells_surr_cell_ptr(:)
type(face_t), allocatable, intent(inout) :: face(:)
type(fp_t), allocatable, intent(inout) :: fp
!
!.. Optional Arguments ..
logical(lk), optional, intent(in) :: force_flxpt_offset_to_zero
!
!.. Local Scalars ..
integer :: n,n1,n2,nf,nc,ierr
integer :: i,ic,ifa,ip1,ip2,ip3,jc
integer :: ibeg,iend,interval
integer :: gmin,gmax,omin,omax
integer :: face_counter,nfp,l1,r1,l2,r2
integer :: left_cell,left_geom
integer :: host_cell,ghost_cell,host_geom
integer :: host_side,nfnods,ncnods
integer :: right_cell,right_geom,ncf
integer :: left_side,right_side
integer :: face_geom,face_order,face_pts
integer :: c1,c2,this_cell,this_geom,adj_cell
integer :: left_order,right_order
integer :: rotation
character(len=200) :: array_name
logical(lk) :: error_finding_rotation
logical(lk) :: flxpt_offset_is_zero
!
integer :: nr,lnode,lfbnd,lface,lcell
!
real(wp) :: dir
!
!.. Local Arrays ..
integer, dimension(1:4) :: ip
integer, dimension(1:4) :: left_nodes
integer, dimension(1:4) :: right_nodes
!integer, dimension(1:4) :: face_nodes
real(wp), dimension(1:3) :: rnrm
real(wp), dimension(1:3) :: face_center
real(wp), dimension(1:3) :: cell_center
real(wp), dimension(1:3,1:4) :: xyz_face
real(wp), dimension(1:3,1:8) :: xyz_cell
!integer, dimension(1:maxFP) :: rot_fpts
!
!.. Local Allocatable Arrays ..
integer, allocatable, dimension(:) :: lpoin
integer, allocatable, dimension(:,:) :: faces_of_cell
!
!.. Local Parameters ..
character(len=*), parameter :: pname = "create_face_derived_type"
!
continue
!
call debug_timer(entering_procedure,pname)
!
! Get the grid dimensions from the input arrays
!
nr = size(xyz_nodes,dim=1)
lnode = size(xyz_nodes,dim=2)
lfbnd = size(bface,dim=2)
lcell = size(nodes_of_cell_ptr)-1 - lfbnd
!
! Initialize the flxpts array
!
call initialize_flxpts(fp)
!
! ###########
! ###########
! ### 1 ###
! ###########
! ###########
call debug_timer(start_timer,"section 1")
!
! Find the total number of faces
! and then allocate the face array
!
lface = lfbnd
do n = 1,lcell
!
n1 = cells_surr_cell_ptr(n)+1
n2 = cells_surr_cell_ptr(n+1)
!
lface = lface + count( cells_surr_cell(n1:n2)>n .and. &
cells_surr_cell(n1:n2)<=lcell )
!
end do
!
allocate ( face(1:lface) , stat=ierr , errmsg=error_message )
call alloc_error(pname,"face",1,__LINE__,__FILE__,ierr,error_message)
!
allocate ( faces_of_cell(1:6,1:lcell) , source=-99 , &
stat=ierr , errmsg=error_message )
call alloc_error(pname,"faces_of_cell",1,__LINE__,__FILE__,ierr,error_message)
!
call debug_timer(stop_timer,"section 1")
! ###########
! ###########
! ### 2 ###
! ###########
! ###########
call debug_timer(start_timer,"section 2")
!
! Get the cells on each boundary face
!
do nf = 1,lfbnd
!
host_cell = bface(2,nf) ! Host cell on boundary face
ghost_cell = bface(10,nf) ! Ghost cell on boundary face
!
! Store the cells on this boundary face
!
face(nf)%left%cell = host_cell
face(nf)%right%cell = ghost_cell
!
host_geom = cell_geom(host_cell) ! geometry of host cell
host_side = bface(4,nf) ! face of host cell on boundary
nfnods = num_face_nodes(host_side,host_geom) ! # of nodes defining bnd face
!
! beginning index of nodes defining host cell
n1 = nodes_of_cell_ptr(host_cell)+1
! ending index of nodes defining host cell
n2 = nodes_of_cell_ptr(host_cell+1)
!
! Node offsets for this face
ip(1:nfnods) = side_nodes(1:nfnods,host_side,host_geom)
!
! Store the nodes on this boundary face
!
face(nf)%nodes(1:nfnods) = nodes_of_cell( n1 + ip(1:nfnods) )
!
! Store the geometry of this face
!
face(nf)%geom = geom_of_face(nfnods)
!
! Store the side of the left cell that is on this face since we have this
! information available right now. Store the side of the right/ghost cell
! as well since this will always be 1.
!