Rocstar  1.0
Rocstar multiphysics simulation application
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
RFLU_ModCOBALT.F90
Go to the documentation of this file.
1 ! *********************************************************************
2 ! * Rocstar Simulation Suite *
3 ! * Copyright@2015, Illinois Rocstar LLC. All rights reserved. *
4 ! * *
5 ! * Illinois Rocstar LLC *
6 ! * Champaign, IL *
7 ! * www.illinoisrocstar.com *
8 ! * sales@illinoisrocstar.com *
9 ! * *
10 ! * License: See LICENSE file in top level of distribution package or *
11 ! * http://opensource.org/licenses/NCSA *
12 ! *********************************************************************
13 ! *********************************************************************
14 ! * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, *
15 ! * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES *
16 ! * OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND *
17 ! * NONINFRINGEMENT. IN NO EVENT SHALL THE CONTRIBUTORS OR *
18 ! * COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
19 ! * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, *
20 ! * Arising FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
21 ! * USE OR OTHER DEALINGS WITH THE SOFTWARE. *
22 ! *********************************************************************
23 ! ******************************************************************************
24 !
25 ! Purpose: Collection of routines to read and convert COBALT grids.
26 !
27 ! Description: None.
28 !
29 ! Notes: None.
30 !
31 ! ******************************************************************************
32 !
33 ! $Id: RFLU_ModCOBALT.F90,v 1.5 2008/12/06 08:45:03 mtcampbe Exp $
34 !
35 ! Copyright: (c) 2004-2006 by the University of Illinois
36 !
37 ! ******************************************************************************
38 
39 MODULE rflu_modcobalt
40 
41  USE modparameters
42  USE moddatatypes
43  USE moderror
44  USE modglobal, ONLY: t_global
45  USE modgrid, ONLY: t_grid
46  USE modbndpatch, ONLY: t_patch
47  USE moddatastruct, ONLY: t_region
48  USE modmpi
49 
50  USE modbuildfilenames, ONLY: buildfilenameplain
51 
52  USE rflu_moddimensions, ONLY: rflu_setmaxdimension
53 
54  IMPLICIT NONE
55 
56 ! ******************************************************************************
57 ! Declarations and definitions
58 ! ******************************************************************************
59 
60  PRIVATE
61 
62 ! ==============================================================================
63 ! Data
64 ! ==============================================================================
65 
66  CHARACTER(CHRLEN) :: &
67  RCSIdentString = '$RCSfile: RFLU_ModCOBALT.F90,v $ $Revision: 1.5 $'
68 
69  TYPE t_gridcobalt
70  INTEGER :: nFaces,nFacesPerCellMax,nMappings,nPatches,nQuads,nTris, &
71  nVertPerFaceMax
72  INTEGER, DIMENSION(:), POINTER :: nvpf
73  INTEGER, DIMENSION(:,:), POINTER :: f2c,f2v,patch2bc
74  END TYPE t_gridcobalt
75 
76  TYPE(t_gridcobalt) :: gridCOBALT
77 
78 ! ==============================================================================
79 ! Public procedures
80 ! ==============================================================================
81 
82  PUBLIC :: rflu_convcobalt2rocflu, &
83  rflu_readgridcobalt
84 
85 ! ******************************************************************************
86 ! Routines
87 ! ******************************************************************************
88 
89  CONTAINS
90 
91 
92 
93 
94 
95 
96 
97 
98 ! ******************************************************************************
99 !
100 ! Purpose: Convert grid format from COBALT to ROCFLU.
101 !
102 ! Description: None.
103 !
104 ! Input:
105 ! pRegion Pointer to region
106 !
107 ! Output: None.
108 !
109 ! Notes:
110 ! 1. This routine is very long, but easily understood due to structure and
111 ! extensive comments. Breaking it into pieces might be more confusing.
112 !
113 ! ******************************************************************************
114 
115  SUBROUTINE rflu_convcobalt2rocflu(pRegion)
116 
117  USE modsortsearch
118 
119  USE rflu_modcellmapping
120  USE rflu_moddimensions, ONLY: rflu_setmaxdimensions
121  USE rflu_modgeometry, ONLY: rflu_computeapproxcentroids, &
122  rflu_createapproxcentroids, &
123  rflu_destroyapproxcentroids
124  USE rflu_modgrid
125 
126  USE modinterfaces, ONLY: facecentroidtria, &
127  facecentroidquad, &
128  facevectortria, &
129  facevectorquad
130 
131  IMPLICIT NONE
132 
133 ! ******************************************************************************
134 ! Declarations and definitions
135 ! ******************************************************************************
136 
137 ! ==============================================================================
138 ! Arguments
139 ! ==============================================================================
140 
141  TYPE(t_region), POINTER :: pregion
142 
143 ! ==============================================================================
144 ! Local variables
145 ! ==============================================================================
146 
147  CHARACTER(CHRLEN) :: ifilename
148  INTEGER :: cntr,c1,c2,errorflag,ibegmax,ibegmin,ibeg1,ibeg2,icg,icl, &
149  iendmax,iendmin,iend1,iend2,ifg,ifile,ifl,iloc,imap,imap2, &
150  ipatch,ivg,ivl,j,v1,v2,v3,v4
151  INTEGER :: flag(4),key(4),vert(4),verttemp(4)
152  INTEGER, DIMENSION(:), ALLOCATABLE :: cellmapp,celltype
153  INTEGER, DIMENSION(:,:), ALLOCATABLE :: celldegr
154  REAL(RFREAL) :: cofgappx,cofgappy,cofgappz,dotprod,fcenx,fceny,fcenz, &
155  fvecx,fvecy,fvecz
156  REAL(RFREAL) :: xyznodes(xcoord:zcoord,4)
157  TYPE(t_grid), POINTER :: pgrid
158  TYPE(t_patch), POINTER :: ppatch
159  TYPE(t_global), POINTER :: global
160 
161 ! ******************************************************************************
162 ! Start
163 ! ******************************************************************************
164 
165  global => pregion%global
166 
167  CALL registerfunction(global,'RFLU_ConvCOBALT2ROCFLU', &
168  'RFLU_ModCOBALT.F90')
169 
170  IF ( global%verbLevel > verbose_none ) THEN
171  WRITE(stdout,'(A,1X,A)') solver_name, &
172  'Converting from COBALT to ROCFLU format...'
173  END IF ! global%verbLevel
174 
175 ! ******************************************************************************
176 ! Set grid pointer and initialize variables
177 ! ******************************************************************************
178 
179  pgrid => pregion%grid
180 
181  pgrid%nEdges = 0
182  pgrid%nEdgesTot = 0
183 
184  pgrid%nFaces = 0
185  pgrid%nFacesTot = 0
186 
187 ! ******************************************************************************
188 ! Allocate temporary memory
189 ! ******************************************************************************
190 
191  ALLOCATE(celldegr(face_type_tri:face_type_quad,gridcobalt%nFaces), &
192  stat=errorflag)
193  global%error = errorflag
194  IF ( global%error /= err_none ) THEN
195  CALL errorstop(global,err_allocate,__line__,'cellDegr')
196  END IF ! global%error
197 
198  DO icg = 1,pgrid%nCellsTot
199  celldegr(face_type_tri,icg) = 0
200  celldegr(face_type_quad,icg) = 0
201  END DO ! icg
202 
203  ALLOCATE(celltype(pgrid%nCellsTot),stat=errorflag)
204  global%error = errorflag
205  IF ( global%error /= err_none ) THEN
206  CALL errorstop(global,err_allocate,__line__,'cellType')
207  END IF ! global%error
208 
209  ALLOCATE(cellmapp(pgrid%nCellsTot),stat=errorflag)
210  global%error = errorflag
211  IF ( global%error /= err_none ) THEN
212  CALL errorstop(global,err_allocate,__line__,'cellMapp')
213  END IF ! global%error
214 
215 ! ******************************************************************************
216 ! Determine number of cells for each cell type
217 ! ******************************************************************************
218 
219  IF ( global%verbLevel > verbose_none ) THEN
220  WRITE(stdout,'(A,3X,A)') solver_name,'Determining cell types...'
221  END IF ! global%verbLevel
222 
223 ! ==============================================================================
224 ! Determine number of triangular and quadrilateral faces per cell
225 ! ==============================================================================
226 
227  DO ifg = 1,gridcobalt%nFaces
228  c1 = gridcobalt%f2c(1,ifg)
229  c2 = gridcobalt%f2c(2,ifg)
230 
231  IF ( gridcobalt%nvpf(ifg) == 3 ) THEN
232  IF ( c1 > 0 ) THEN
233  celldegr(face_type_tri,c1) = celldegr(face_type_tri,c1) + 1
234  END IF ! c1
235 
236  IF ( c2 > 0 ) THEN
237  celldegr(face_type_tri,c2) = celldegr(face_type_tri,c2) + 1
238  END IF ! c2
239  ELSE IF ( gridcobalt%nvpf(ifg) == 4 ) THEN
240  IF ( c1 > 0 ) THEN
241  celldegr(face_type_quad,c1) = celldegr(face_type_quad,c1) + 1
242  END IF ! c1
243 
244  IF ( c2 > 0 ) THEN
245  celldegr(face_type_quad,c2) = celldegr(face_type_quad,c2) + 1
246  END IF ! c2
247  ELSE ! Should be impossible
248  CALL errorstop(global,err_reached_default,__line__)
249  END IF ! gridCOBALT%nvpf
250  END DO ! ifg
251 
252 ! ==============================================================================
253 ! Determine cell types using number of triangular and quadrilateral faces
254 ! ==============================================================================
255 
256  pgrid%nTetsTot = 0
257  pgrid%nHexsTot = 0
258  pgrid%nPrisTot = 0
259  pgrid%nPyrsTot = 0
260 
261  DO icg = 1,pgrid%nCellsTot
262  IF ( celldegr(face_type_tri ,icg) == 4 .AND. &
263  celldegr(face_type_quad,icg) == 0 ) THEN ! Tetrahedron
264  pgrid%nTetsTot = pgrid%nTetsTot + 1
265  celltype(icg) = cell_type_tet
266  cellmapp(icg) = pgrid%nTetsTot
267  ELSE IF ( celldegr(face_type_tri ,icg) == 0 .AND. &
268  celldegr(face_type_quad,icg) == 6 ) THEN ! Hexahedron
269  pgrid%nHexsTot = pgrid%nHexsTot + 1
270  celltype(icg) = cell_type_hex
271  cellmapp(icg) = pgrid%nHexsTot
272  ELSE IF ( celldegr(face_type_tri ,icg) == 2 .AND. &
273  celldegr(face_type_quad,icg) == 3 ) THEN ! Prism
274  pgrid%nPrisTot = pgrid%nPrisTot + 1
275  celltype(icg) = cell_type_pri
276  cellmapp(icg) = pgrid%nPrisTot
277  ELSE IF ( celldegr(face_type_tri ,icg) == 4 .AND. &
278  celldegr(face_type_quad,icg) == 1 ) THEN ! Pyramid
279  pgrid%nPyrsTot = pgrid%nPyrsTot + 1
280  celltype(icg) = cell_type_pyr
281  cellmapp(icg) = pgrid%nPyrsTot
282  ELSE
283  CALL errorstop(global,err_reached_default,__line__)
284  END IF ! cellDegr
285  END DO ! icg
286 
287  pgrid%nTets = pgrid%nTetsTot
288  pgrid%nHexs = pgrid%nHexsTot
289  pgrid%nPris = pgrid%nPrisTot
290  pgrid%nPyrs = pgrid%nPyrsTot
291 
292  pgrid%nTetsMax = rflu_setmaxdimension(global,pgrid%nTetsTot)
293  pgrid%nHexsMax = rflu_setmaxdimension(global,pgrid%nHexsTot)
294  pgrid%nPrisMax = rflu_setmaxdimension(global,pgrid%nPrisTot)
295  pgrid%nPyrsMax = rflu_setmaxdimension(global,pgrid%nPyrsTot)
296 
297 ! ===============================================================================
298 ! Write statistics
299 ! ===============================================================================
300 
301  IF ( global%verbLevel > verbose_low ) THEN
302  WRITE(stdout,'(A,5X,A)') solver_name,'Cell statistics:'
303  WRITE(stdout,'(A,7X,A,1X,I9)') solver_name,'Tetrahedra:',pgrid%nTetsTot
304  WRITE(stdout,'(A,7X,A,1X,I9)') solver_name,'Hexahedra: ',pgrid%nHexsTot
305  WRITE(stdout,'(A,7X,A,1X,I9)') solver_name,'Prisms: ',pgrid%nPrisTot
306  WRITE(stdout,'(A,7X,A,1X,I9)') solver_name,'Pyramids: ',pgrid%nPyrsTot
307  END IF ! global%verbLevel
308 
309  IF ( global%verbLevel > verbose_none ) THEN
310  WRITE(stdout,'(A,3X,A)') solver_name,'Determining cell types done.'
311  END IF ! global%verbLevel
312 
313 ! ******************************************************************************
314 ! Set maximum dimensions. NOTE this is necessary because dimensioning of
315 ! arrays is in terms of maximum dimensions.
316 ! ******************************************************************************
317 
318  CALL rflu_setmaxdimensions(pregion)
319 
320 ! ******************************************************************************
321 ! Renumber cells. NOTE this is necessary because Rocflu assumes that cells are
322 ! numbered in the order tetrahedra, hexahedra, prisms, and pyramids, but
323 ! Cobalt violates this assumption.
324 ! ******************************************************************************
325 
326  IF ( global%verbLevel > verbose_none ) THEN
327  WRITE(stdout,'(A,3X,A)') solver_name,'Renumbering cells...'
328  END IF ! global%verbLevel
329 
330 ! ==============================================================================
331 ! Build cell mapping
332 ! ==============================================================================
333 
334  CALL rflu_createcellmapping(pregion)
335  CALL rflu_buildloc2globcellmapping(pregion)
336  CALL rflu_buildglob2loccellmapping(pregion)
337 
338 ! ==============================================================================
339 ! Renumber cells, make face list consistent
340 ! ==============================================================================
341 
342  DO ifg = 1,gridcobalt%nFaces
343  DO j = 1,2
344  c1 = gridcobalt%f2c(j,ifg)
345 
346  IF ( c1 > 0 ) THEN
347  icl = cellmapp(c1)
348 
349  SELECT CASE ( celltype(c1) )
350  CASE ( cell_type_tet )
351  icg = pgrid%tet2CellGlob(icl)
352  CASE ( cell_type_hex )
353  icg = pgrid%hex2CellGlob(icl)
354  CASE ( cell_type_pri )
355  icg = pgrid%pri2CellGlob(icl)
356  CASE ( cell_type_pyr )
357  icg = pgrid%pyr2CellGlob(icl)
358  CASE default
359  CALL errorstop(global,err_reached_default,__line__)
360  END SELECT ! cellType
361 
362  gridcobalt%f2c(j,ifg) = icg
363  END IF ! c1
364  END DO ! j
365  END DO ! ifg
366 
367 ! ==============================================================================
368 ! Deallocate part of temporary memory
369 ! ==============================================================================
370 
371  DEALLOCATE(cellmapp,stat=errorflag)
372  global%error = errorflag
373  IF ( global%error /= err_none ) THEN
374  CALL errorstop(global,err_deallocate,__line__,'cellMapp')
375  END IF ! global%error
376 
377  DEALLOCATE(celltype,stat=errorflag)
378  global%error = errorflag
379  IF ( global%error /= err_none ) THEN
380  CALL errorstop(global,err_deallocate,__line__,'cellType')
381  END IF ! global%error
382 
383  IF ( global%verbLevel > verbose_none ) THEN
384  WRITE(stdout,'(A,3X,A)') solver_name,'Renumbering cells done.'
385  END IF ! global%verbLevel
386 
387 ! ******************************************************************************
388 ! Building cell connectivity. NOTE in the following, the cell connectivity is
389 ! built first without checking for correctness for convenience. The
390 ! connectivity will be checked later in two passes.
391 ! ******************************************************************************
392 
393  IF ( global%verbLevel > verbose_none ) THEN
394  WRITE(stdout,'(A,3X,A)') solver_name,'Building cell connectivity...'
395  END IF ! global%verbLevel
396 
397  IF ( pgrid%nTetsMax > 0 ) THEN
398  ALLOCATE(pgrid%tet2v(4,pgrid%nTetsMax),stat=errorflag)
399  global%error = errorflag
400  IF ( global%error /= err_none ) THEN
401  CALL errorstop(global,err_allocate,__line__,'pGrid%tet2v')
402  END IF ! global%error
403 
404  DO icl = 1,pgrid%nTetsMax
405  pgrid%tet2v(1,icl) = c2v_init
406  pgrid%tet2v(2,icl) = c2v_init
407  pgrid%tet2v(3,icl) = c2v_init
408  pgrid%tet2v(4,icl) = c2v_init
409  END DO ! icl
410  END IF ! pGrid%nTetsMax
411 
412  IF ( pgrid%nHexsMax > 0 ) THEN
413  ALLOCATE(pgrid%hex2v(8,pgrid%nHexsMax),stat=errorflag)
414  global%error = errorflag
415  IF ( global%error /= err_none ) THEN
416  CALL errorstop(global,err_allocate,__line__,'pGrid%hex2v')
417  END IF ! global%error
418 
419  DO icl = 1,pgrid%nHexsMax
420  pgrid%hex2v(1,icl) = c2v_init
421  pgrid%hex2v(2,icl) = c2v_init
422  pgrid%hex2v(3,icl) = c2v_init
423  pgrid%hex2v(4,icl) = c2v_init
424  pgrid%hex2v(5,icl) = c2v_init
425  pgrid%hex2v(6,icl) = c2v_init
426  pgrid%hex2v(7,icl) = c2v_init
427  pgrid%hex2v(8,icl) = c2v_init
428  END DO ! icl
429  END IF ! pGrid%nHexsMax
430 
431  IF ( pgrid%nPrisMax > 0 ) THEN
432  ALLOCATE(pgrid%pri2v(6,pgrid%nPrisMax),stat=errorflag)
433  global%error = errorflag
434  IF ( global%error /= err_none ) THEN
435  CALL errorstop(global,err_allocate,__line__,'pGrid%pri2v')
436  END IF ! global%error
437 
438  DO icl = 1,pgrid%nPrisMax
439  pgrid%pri2v(1,icl) = c2v_init
440  pgrid%pri2v(2,icl) = c2v_init
441  pgrid%pri2v(3,icl) = c2v_init
442  pgrid%pri2v(4,icl) = c2v_init
443  pgrid%pri2v(5,icl) = c2v_init
444  pgrid%pri2v(6,icl) = c2v_init
445  END DO ! icl
446  END IF ! pGrid%nPrisMax
447 
448  IF ( pgrid%nPyrsMax > 0 ) THEN
449  ALLOCATE(pgrid%pyr2v(5,pgrid%nPyrsMax),stat=errorflag)
450  global%error = errorflag
451  IF ( global%error /= err_none ) THEN
452  CALL errorstop(global,err_allocate,__line__,'pGrid%pyr2v')
453  END IF ! global%error
454 
455  DO icl = 1,pgrid%nPyrsMax
456  pgrid%pyr2v(1,icl) = c2v_init
457  pgrid%pyr2v(2,icl) = c2v_init
458  pgrid%pyr2v(3,icl) = c2v_init
459  pgrid%pyr2v(4,icl) = c2v_init
460  pgrid%pyr2v(5,icl) = c2v_init
461  END DO ! icl
462  END IF ! pGrid%nPyrsMax
463 
464 ! ==============================================================================
465 ! Loop over faces and visit interior cells for each face. NOTE do not use
466 ! quadrilateral faces to define prism because it is easier to use the
467 ! triangular faces (no overlap of vertices) and then use the quadrilateral
468 ! faces to get the correct relative orientation. NOTE for pyramids, if first
469 ! face to be added is a triangular face, store in positions 1 to 3 of pyr2v
470 ! array despite these being the wrong locations. When quadrilateral face is
471 ! added to existing triangular face, make sure that these 3 vertices are
472 ! taken into account.
473 ! ==============================================================================
474 
475  DO icg = 1,pgrid%nCellsTot ! Reinitialize
476  celldegr(face_type_tri,icg) = 0
477  celldegr(face_type_quad,icg) = 0
478  END DO ! icg
479 
480  DO ifg = 1,gridcobalt%nFaces
481  DO j = 1,2
482  c1 = gridcobalt%f2c(j,ifg)
483 
484  IF ( c1 > 0 ) THEN ! interior cell
485 
486 ! ------------------------------------------------------------------------------
487 ! Triangular face
488 ! ------------------------------------------------------------------------------
489 
490  IF ( gridcobalt%nvpf(ifg) == 3 ) THEN
491  SELECT CASE ( pgrid%cellGlob2Loc(1,c1) )
492 
493 ! ----------- Tetrahedron ------------------------------------------------------
494 
495  CASE ( cell_type_tet )
496  IF ( celldegr(face_type_tri,c1) < 2 ) THEN
497  celldegr(face_type_tri,c1) = celldegr(face_type_tri,c1) + 1
498 
499  icl = pgrid%cellGlob2Loc(2,c1)
500 
501  IF ( celldegr(face_type_tri,c1) == 1 ) THEN
502  pgrid%tet2v(1:3,icl) = gridcobalt%f2v(1:3,ifg)
503  ELSE
504  vert(1:3) = pgrid%tet2v(1:3,icl)
505  CALL quicksortinteger(vert(1:3),3)
506 
507  DO ivl = 1,3
508  CALL binarysearchinteger(vert,3,gridcobalt%f2v(ivl,ifg), &
509  iloc)
510 
511  IF ( iloc == element_not_found ) THEN
512  pgrid%tet2v(4,icl) = gridcobalt%f2v(ivl,ifg)
513  EXIT
514  END IF ! iloc
515  END DO ! ivl
516 
517  IF ( pgrid%tet2v(4,icl) == c2v_init ) THEN
518  CALL errorstop(global,err_c2vlist_invalid,__line__)
519  END IF ! pGrid%tet2v
520  END IF ! cellDegr
521  END IF ! cellDegr
522 
523 ! ----------- Prism ------------------------------------------------------------
524 
525  CASE ( cell_type_pri )
526  IF ( celldegr(face_type_tri,c1) < 2 ) THEN
527  celldegr(face_type_tri,c1) = celldegr(face_type_tri,c1) + 1
528 
529  icl = pgrid%cellGlob2Loc(2,c1)
530 
531  IF ( celldegr(face_type_tri,c1) == 1 ) THEN
532  pgrid%pri2v(1:3,icl) = gridcobalt%f2v(1:3,ifg)
533  ELSE
534  pgrid%pri2v(4:6,icl) = gridcobalt%f2v(1:3,ifg)
535  END IF ! cellDegr
536  END IF ! cellDegr
537 
538 ! ----------- Pyramid ----------------------------------------------------------
539 
540  CASE ( cell_type_pyr )
541  IF ( celldegr(face_type_tri,c1) < 1 ) THEN
542  celldegr(face_type_tri,c1) = celldegr(face_type_tri,c1) + 1
543 
544  icl = pgrid%cellGlob2Loc(2,c1)
545 
546  IF ( celldegr(face_type_quad,c1) == 0 ) THEN
547  pgrid%pyr2v(1,icl) = gridcobalt%f2v(1,ifg)
548  pgrid%pyr2v(2,icl) = gridcobalt%f2v(2,ifg)
549  pgrid%pyr2v(3,icl) = gridcobalt%f2v(3,ifg)
550  ELSE
551  vert(1:4) = pgrid%pyr2v(1:4,icl)
552  CALL quicksortinteger(vert(1:4),4)
553 
554  DO ivl = 1,3
555  CALL binarysearchinteger(vert(1:4),4, &
556  gridcobalt%f2v(ivl,ifg),iloc)
557 
558  IF ( iloc == element_not_found ) THEN
559  pgrid%pyr2v(5,icl) = gridcobalt%f2v(ivl,ifg)
560  END IF ! iloc
561  END DO ! ivl
562 
563  IF ( pgrid%pyr2v(5,icl) == c2v_init ) THEN
564  CALL errorstop(global,err_c2vlist_invalid,__line__)
565  END IF ! pGrid%pyr2v
566  END IF ! cellDegr
567  END IF ! cellDegr
568  END SELECT ! pGrid%cellGlob2Loc
569 
570 ! ------------------------------------------------------------------------------
571 ! Quadrilateral face
572 ! ------------------------------------------------------------------------------
573 
574  ELSE IF ( gridcobalt%nvpf(ifg) == 4 ) THEN
575  SELECT CASE ( pgrid%cellGlob2Loc(1,c1) )
576 
577 ! ----------- Hexahedron -------------------------------------------------------
578 
579  CASE ( cell_type_hex )
580  IF ( celldegr(face_type_quad,c1) < 2 ) THEN
581  icl = pgrid%cellGlob2Loc(2,c1)
582 
583  IF ( celldegr(face_type_quad,c1) == 0 ) THEN ! Face 1
584  celldegr(face_type_quad,c1) = celldegr(face_type_quad,c1) &
585  + 1
586  pgrid%hex2v(1:4,icl) = gridcobalt%f2v(1:4,ifg)
587  ELSE
588  vert(1:4) = pgrid%hex2v(1:4,icl)
589  CALL quicksortinteger(vert(1:4),4)
590 
591  cntr = 0
592 
593  DO ivl = 1,4
594  ivg = gridcobalt%f2v(ivl,ifg)
595 
596  CALL binarysearchinteger(vert(1:4),4,ivg,iloc)
597 
598  IF ( iloc == element_not_found ) THEN
599  cntr = cntr + 1
600  END IF ! iloc
601  END DO ! ivl
602 
603  IF ( cntr == 4 ) THEN ! Face 6
604  celldegr(face_type_quad,c1) = &
605  celldegr(face_type_quad,c1) + 1
606  pgrid%hex2v(5:8,icl) = gridcobalt%f2v(1:4,ifg)
607  END IF ! cntr
608  END IF ! cellDegr
609  END IF ! cellDegr
610 
611 ! ----------- Pyramid ----------------------------------------------------------
612 
613  CASE ( cell_type_pyr )
614  IF ( celldegr(face_type_quad,c1) < 1 ) THEN
615  celldegr(face_type_quad,c1) = celldegr(face_type_quad,c1) + 1
616 
617  icl = pgrid%cellGlob2Loc(2,c1)
618 
619  IF ( celldegr(face_type_tri,c1) == 0 ) THEN
620  pgrid%pyr2v(1:4,icl) = gridcobalt%f2v(1:4,ifg)
621  ELSE
622  verttemp(1:3) = pgrid%pyr2v(1:3,icl)
623 
624  pgrid%pyr2v(1:4,icl) = gridcobalt%f2v(1:4,ifg)
625  vert(1:4) = pgrid%pyr2v(1:4,icl)
626  CALL quicksortinteger(vert(1:4),4)
627 
628  DO ivl = 1,3
629  CALL binarysearchinteger(vert(1:4),4, &
630  verttemp(ivl),iloc)
631 
632  IF ( iloc == element_not_found ) THEN
633  pgrid%pyr2v(5,icl) = verttemp(ivl)
634  END IF ! iloc
635  END DO ! ivl
636 
637  IF ( pgrid%pyr2v(5,icl) == c2v_init ) THEN
638  CALL errorstop(global,err_c2vlist_invalid,__line__)
639  END IF ! pGrid%pyr2v
640  END IF ! cellDegr
641  END IF ! cellDegr
642  END SELECT ! pGrid%cellGlob2Loc
643 
644  ELSE ! Should be impossible by now, but leave for defensive coding
645  CALL errorstop(global,err_reached_default,__line__)
646  END IF ! gridCOBALT%nvpf
647  END IF ! c1
648  END DO ! j
649  END DO ! ifg
650 
651  IF ( global%verbLevel > verbose_none ) THEN
652  WRITE(stdout,'(A,3X,A)') solver_name,'Building cell connectivity done.'
653  END IF ! global%verbLevel
654 
655 ! ******************************************************************************
656 ! Checking connectivity, pass 1 (checking orientation). For cells other than
657 ! tetrahedra, this is only done for those faces which were added above, and
658 ! for pyramids, it is only done for the quadrilateral face. It is not
659 ! necessary for the triangular faces of the pyramids to be checked because
660 ! reversing the order would mean that the order of the quadrilateral face is
661 ! wrong.
662 ! ******************************************************************************
663 
664 ! ==============================================================================
665 ! Compute approximate centroids. NOTE cell mapping (built above) required for
666 ! computation of approximate centroids.
667 ! ==============================================================================
668 
669  CALL rflu_createapproxcentroids(pregion)
670  CALL rflu_computeapproxcentroids(pregion)
671 
672 ! ==============================================================================
673 ! Check connectivity by checking orientation of faces. This is done by
674 ! computing the dot product between the face normal and the relative position
675 ! vector between the cell centroid and the face centroid. This must be
676 ! positive if the face normal is outward pointing, which means that the face
677 ! vertices are oriented properly.
678 ! ==============================================================================
679 
680  IF ( global%verbLevel > verbose_none ) THEN
681  WRITE(stdout,'(A,3X,A)') solver_name,'Checking cell connectivity...'
682 
683  IF ( global%verbLevel > verbose_low ) THEN
684  WRITE(stdout,'(A,5X,A)') solver_name,'Pass 1...'
685  END IF ! global%verbLevel
686  END IF ! global%verbLevel
687 
688 ! ------------------------------------------------------------------------------
689 ! Tetrahedra
690 ! ------------------------------------------------------------------------------
691 
692  IF ( pgrid%nTetsTot > 0 ) THEN
693  IF ( global%verbLevel > verbose_low ) THEN
694  WRITE(stdout,'(A,7X,A)') solver_name,'Tetrahedra...'
695  END IF ! global%verbLevel
696 
697  DO icl = 1,pgrid%nTetsTot
698  icg = pgrid%tet2CellGlob(icl)
699 
700  cofgappx = pgrid%cofgApp(xcoord,icg)
701  cofgappy = pgrid%cofgApp(ycoord,icg)
702  cofgappz = pgrid%cofgApp(zcoord,icg)
703 
704  cntr = 0
705 
706  DO ifl = 1,4
707  v1 = pgrid%tet2v(f2vtet(1,ifl),icl)
708  v2 = pgrid%tet2v(f2vtet(2,ifl),icl)
709  v3 = pgrid%tet2v(f2vtet(3,ifl),icl)
710 
711  xyznodes(xcoord:zcoord,1) = pgrid%xyz(xcoord:zcoord,v1)
712  xyznodes(xcoord:zcoord,2) = pgrid%xyz(xcoord:zcoord,v2)
713  xyznodes(xcoord:zcoord,3) = pgrid%xyz(xcoord:zcoord,v3)
714 
715  CALL facevectortria(xyznodes(xcoord:zcoord,1:3),fvecx,fvecy,fvecz)
716  CALL facecentroidtria(xyznodes(xcoord:zcoord,1:3),fcenx,fceny,fcenz)
717 
718  dotprod = fvecx*(fcenx - cofgappx) &
719  + fvecy*(fceny - cofgappy) &
720  + fvecz*(fcenz - cofgappz)
721 
722  IF ( dotprod < 0.0_rfreal ) THEN ! Only one face can be wrong
723  cntr = cntr + 1
724 
725  IF ( cntr > 1 ) THEN
726  CALL errorstop(global,err_face_orient,__line__)
727  END IF ! cntr
728 
729  pgrid%tet2v(f2vtet(1,ifl),icl) = v1
730  pgrid%tet2v(f2vtet(2,ifl),icl) = v3
731  pgrid%tet2v(f2vtet(3,ifl),icl) = v2
732  END IF ! dotProd
733  END DO ! ifl
734  END DO ! icl
735  END IF ! pGrid%nTetsTot
736 
737 ! ------------------------------------------------------------------------------
738 ! Hexahedra
739 ! ------------------------------------------------------------------------------
740 
741  IF ( pgrid%nHexsTot > 0 ) THEN
742  IF ( global%verbLevel > verbose_low ) THEN
743  WRITE(stdout,'(A,7X,A)') solver_name,'Hexahedra...'
744  END IF ! global%verbLevel
745 
746  DO icl = 1,pgrid%nHexsTot
747  icg = pgrid%hex2CellGlob(icl)
748 
749  cofgappx = pgrid%cofgApp(xcoord,icg)
750  cofgappy = pgrid%cofgApp(ycoord,icg)
751  cofgappz = pgrid%cofgApp(zcoord,icg)
752 
753  DO ifl = 1,6,5 ! Note loop over opposing faces
754  v1 = pgrid%hex2v(f2vhex(1,ifl),icl)
755  v2 = pgrid%hex2v(f2vhex(2,ifl),icl)
756  v3 = pgrid%hex2v(f2vhex(3,ifl),icl)
757  v4 = pgrid%hex2v(f2vhex(4,ifl),icl)
758 
759  xyznodes(xcoord:zcoord,1) = pgrid%xyz(xcoord:zcoord,v1)
760  xyznodes(xcoord:zcoord,2) = pgrid%xyz(xcoord:zcoord,v2)
761  xyznodes(xcoord:zcoord,3) = pgrid%xyz(xcoord:zcoord,v3)
762  xyznodes(xcoord:zcoord,4) = pgrid%xyz(xcoord:zcoord,v4)
763 
764  CALL facevectorquad(xyznodes(xcoord:zcoord,1:4),fvecx,fvecy,fvecz)
765  CALL facecentroidquad(xyznodes(xcoord:zcoord,1:4),fcenx,fceny,fcenz)
766 
767  dotprod = fvecx*(fcenx - cofgappx) &
768  + fvecy*(fceny - cofgappy) &
769  + fvecz*(fcenz - cofgappz)
770 
771  IF ( dotprod < 0.0_rfreal ) THEN
772  pgrid%hex2v(f2vhex(1,ifl),icl) = v1
773  pgrid%hex2v(f2vhex(2,ifl),icl) = v4
774  pgrid%hex2v(f2vhex(3,ifl),icl) = v3
775  pgrid%hex2v(f2vhex(4,ifl),icl) = v2
776  END IF ! dotProd
777  END DO ! ifl
778  END DO ! icl
779  END IF ! pGrid%nHexsTot
780 
781 ! ------------------------------------------------------------------------------
782 ! Prisms
783 ! ------------------------------------------------------------------------------
784 
785  IF ( pgrid%nPrisTot > 0 ) THEN
786  IF ( global%verbLevel > verbose_low ) THEN
787  WRITE(stdout,'(A,7X,A)') solver_name,'Prisms...'
788  END IF ! global%verbLevel
789 
790  DO icl = 1,pgrid%nPrisTot
791  icg = pgrid%pri2CellGlob(icl)
792 
793  cofgappx = pgrid%cofgApp(xcoord,icg)
794  cofgappy = pgrid%cofgApp(ycoord,icg)
795  cofgappz = pgrid%cofgApp(zcoord,icg)
796 
797  DO ifl = 1,5,4 ! NOTE loop only over triangular faces
798  v1 = pgrid%pri2v(f2vpri(1,ifl),icl)
799  v2 = pgrid%pri2v(f2vpri(2,ifl),icl)
800  v3 = pgrid%pri2v(f2vpri(3,ifl),icl)
801 
802  xyznodes(xcoord:zcoord,1) = pgrid%xyz(xcoord:zcoord,v1)
803  xyznodes(xcoord:zcoord,2) = pgrid%xyz(xcoord:zcoord,v2)
804  xyznodes(xcoord:zcoord,3) = pgrid%xyz(xcoord:zcoord,v3)
805 
806  CALL facevectortria(xyznodes(xcoord:zcoord,1:3),fvecx,fvecy,fvecz)
807  CALL facecentroidtria(xyznodes(xcoord:zcoord,1:3),fcenx,fceny,fcenz)
808 
809  dotprod = fvecx*(fcenx - cofgappx) &
810  + fvecy*(fceny - cofgappy) &
811  + fvecz*(fcenz - cofgappz)
812 
813  IF ( dotprod < 0.0_rfreal ) THEN
814  pgrid%pri2v(f2vpri(1,ifl),icl) = v1
815  pgrid%pri2v(f2vpri(2,ifl),icl) = v3
816  pgrid%pri2v(f2vpri(3,ifl),icl) = v2
817  END IF ! dotProd
818  END DO ! ifl
819  END DO ! icl
820  END IF ! pGrid%nPrisTot
821 
822 ! ------------------------------------------------------------------------------
823 ! Pyramids
824 ! ------------------------------------------------------------------------------
825 
826  IF ( pgrid%nPyrsTot > 0 ) THEN
827  IF ( global%verbLevel > verbose_low ) THEN
828  WRITE(stdout,'(A,7X,A)') solver_name,'Pyramids...'
829  END IF ! global%verbLevel
830 
831  DO icl = 1,pgrid%nPyrsTot
832  icg = pgrid%pyr2CellGlob(icl)
833 
834  cofgappx = pgrid%cofgApp(xcoord,icg)
835  cofgappy = pgrid%cofgApp(ycoord,icg)
836  cofgappz = pgrid%cofgApp(zcoord,icg)
837 
838  DO ifl = 1,1 ! Note loop only over first face
839  v1 = pgrid%pyr2v(f2vpyr(1,ifl),icl)
840  v2 = pgrid%pyr2v(f2vpyr(2,ifl),icl)
841  v3 = pgrid%pyr2v(f2vpyr(3,ifl),icl)
842  v4 = pgrid%pyr2v(f2vpyr(4,ifl),icl)
843 
844  xyznodes(xcoord:zcoord,1) = pgrid%xyz(xcoord:zcoord,v1)
845  xyznodes(xcoord:zcoord,2) = pgrid%xyz(xcoord:zcoord,v2)
846  xyznodes(xcoord:zcoord,3) = pgrid%xyz(xcoord:zcoord,v3)
847  xyznodes(xcoord:zcoord,4) = pgrid%xyz(xcoord:zcoord,v4)
848 
849  CALL facevectorquad(xyznodes(xcoord:zcoord,1:4),fvecx,fvecy,fvecz)
850  CALL facecentroidquad(xyznodes(xcoord:zcoord,1:4),fcenx,fceny,fcenz)
851 
852  dotprod = fvecx*(fcenx - cofgappx) &
853  + fvecy*(fceny - cofgappy) &
854  + fvecz*(fcenz - cofgappz)
855 
856  IF ( dotprod < 0.0_rfreal ) THEN
857  pgrid%pyr2v(f2vpyr(1,ifl),icl) = v1
858  pgrid%pyr2v(f2vpyr(2,ifl),icl) = v4
859  pgrid%pyr2v(f2vpyr(3,ifl),icl) = v3
860  pgrid%pyr2v(f2vpyr(4,ifl),icl) = v2
861  END IF ! dotProd
862  END DO ! ifl
863  END DO ! icl
864  END IF ! pGrid%nPyrsTot
865 
866 ! ******************************************************************************
867 ! Checking orientation, pass 2. At this point, tetrahedra and pyramids are
868 ! oriented properly. However, hexahedra and prisms are not necessarily
869 ! oriented properly, because although the two individual faces are oriented
870 ! properly, they may not be aligned correctly relative to each other. To find
871 ! correct alignment, use third face.
872 ! ******************************************************************************
873 
874  IF ( global%verbLevel > verbose_low ) THEN
875  WRITE(stdout,'(A,5X,A)') solver_name,'Pass 2...'
876  END IF ! global%verbLevel
877 
878 ! ==============================================================================
879 ! Loop over faces
880 ! ==============================================================================
881 
882  DO ifg = 1,gridcobalt%nFaces
883  DO j = 1,2
884  c1 = gridcobalt%f2c(j,ifg)
885 
886  IF ( c1 > 0 ) THEN
887  IF ( gridcobalt%nvpf(ifg) == 4 ) THEN
888  SELECT CASE ( pgrid%cellGlob2Loc(1,c1) )
889 
890 ! ------------------------------------------------------------------------------
891 ! Hexahedra
892 ! ------------------------------------------------------------------------------
893 
894  CASE ( cell_type_hex )
895  IF ( celldegr(face_type_quad,c1) < 3 ) THEN
896  icl = pgrid%cellGlob2Loc(2,c1)
897 
898 ! --------------- Find faces which can be used to determine orientation --------
899 
900  verttemp(1:4) = pgrid%hex2v(1:4,icl)
901  CALL quicksortinteger(verttemp(1:4),4)
902 
903  cntr = 0
904 
905  DO ivl = 1,4
906  ivg = gridcobalt%f2v(ivl,ifg)
907 
908  CALL binarysearchinteger(verttemp(1:4),4,ivg,iloc)
909 
910  IF ( iloc /= element_not_found ) THEN
911  cntr = cntr + 1
912  END IF ! iloc
913  END DO ! ivl
914 
915 ! --------------- Need faces which share 2 vertices with face 1 ----------------
916 
917  IF ( cntr == 2 ) THEN ! Have proper candidate face
918  celldegr(face_type_quad,c1) = celldegr(face_type_quad,c1) &
919  + 1
920 
921  v1 = gridcobalt%f2v(1,ifg)
922  v2 = gridcobalt%f2v(2,ifg)
923  v3 = gridcobalt%f2v(3,ifg)
924  v4 = gridcobalt%f2v(4,ifg)
925 
926  xyznodes(xcoord:zcoord,1) = pgrid%xyz(xcoord:zcoord,v1)
927  xyznodes(xcoord:zcoord,2) = pgrid%xyz(xcoord:zcoord,v2)
928  xyznodes(xcoord:zcoord,3) = pgrid%xyz(xcoord:zcoord,v3)
929  xyznodes(xcoord:zcoord,4) = pgrid%xyz(xcoord:zcoord,v4)
930 
931  CALL facevectorquad(xyznodes(xcoord:zcoord,1:4), &
932  fvecx,fvecy,fvecz)
933  CALL facecentroidquad(xyznodes(xcoord:zcoord,1:4), &
934  fcenx,fceny,fcenz)
935 
936 ! ----------------- Check orientation, reverse if necessary --------------------
937 
938  dotprod = fvecx*(fcenx - pgrid%cofgApp(xcoord,c1)) &
939  + fvecy*(fceny - pgrid%cofgApp(ycoord,c1)) &
940  + fvecz*(fcenz - pgrid%cofgApp(zcoord,c1))
941 
942  IF ( dotprod < 0.0_rfreal ) THEN
943  vert(1) = v1
944  vert(2) = v4
945  vert(3) = v3
946  vert(4) = v2
947  ELSE
948  vert(1) = v1
949  vert(2) = v2
950  vert(3) = v3
951  vert(4) = v4
952  END IF ! dotProd
953 
954 ! ----------------- Determine local face number --------------------------------
955 
956  verttemp(1:4) = pgrid%hex2v(1:4,icl)
957  CALL quicksortinteger(verttemp(1:4),4)
958 
959  DO ivl = 1,4
960  ivg = pgrid%hex2v(ivl,icl)
961 
962  CALL binarysearchinteger(verttemp(1:4),4,ivg,iloc)
963 
964  IF ( iloc /= element_not_found ) THEN
965  key(iloc) = ivl
966  ELSE
967  CALL errorstop(global,err_reached_default,__line__)
968  END IF ! iloc
969  END DO ! ivl
970 
971  cntr = 0
972  flag(1:4) = 0
973 
974  DO ivl = 1,4
975  ivg = vert(ivl)
976 
977  CALL binarysearchinteger(verttemp(1:4),4,ivg,iloc)
978 
979  IF ( iloc /= element_not_found ) THEN
980  flag(key(iloc)) = 1
981  ELSE
982  cntr = cntr + 1
983  END IF ! iloc
984  END DO ! ivl
985 
986 ! ----------------- Cycle vertices of face 6 to match ordering on face 1 -------
987 
988  IF ( flag(1) == 1 .AND. flag(2) == 1 ) THEN ! Face 2
989  CALL cyclelist(vert(1:4),4,1,pgrid%hex2v(1,icl))
990  CALL cyclelist(pgrid%hex2v(5:8,icl),4,1,vert(4))
991  ELSE IF ( flag(2) == 1 .AND. flag(3) == 1 ) THEN ! Face 3
992  CALL cyclelist(vert(1:4),4,1,pgrid%hex2v(2,icl))
993  CALL cyclelist(pgrid%hex2v(5:8,icl),4,2,vert(4))
994  ELSE IF ( flag(3) == 1 .AND. flag(4) == 1 ) THEN ! Face 4
995  CALL cyclelist(vert(1:4),4,1,pgrid%hex2v(3,icl))
996  CALL cyclelist(pgrid%hex2v(5:8,icl),4,3,vert(4))
997  ELSE IF ( flag(4) == 1 .AND. flag(1) == 1 ) THEN ! Face 5
998  CALL cyclelist(vert(1:4),4,1,pgrid%hex2v(4,icl))
999  CALL cyclelist(pgrid%hex2v(5:8,icl),4,4,vert(4))
1000  ELSE
1001  CALL errorstop(global,err_reached_default,__line__)
1002  END IF ! flag
1003  ELSE IF ( cntr /= 4 .AND. cntr /= 0 ) THEN
1004  CALL errorstop(global,err_reached_default,__line__)
1005  END IF ! cntr
1006  END IF ! cellDegr
1007 
1008 ! ------------------------------------------------------------------------------
1009 ! Prism
1010 ! ------------------------------------------------------------------------------
1011 
1012  CASE ( cell_type_pri )
1013  IF ( celldegr(face_type_quad,c1) < 1 ) THEN
1014  icl = pgrid%cellGlob2Loc(2,c1)
1015 
1016 ! --------------- Find faces which can be used to determine orientation --------
1017 
1018  verttemp(1:4) = pgrid%pri2v(1:4,icl)
1019  CALL quicksortinteger(verttemp(1:4),4)
1020 
1021  cntr = 0
1022 
1023  DO ivl = 1,4
1024  ivg = gridcobalt%f2v(ivl,ifg)
1025 
1026  CALL binarysearchinteger(verttemp(1:4),4,ivg,iloc)
1027 
1028  IF ( iloc /= element_not_found ) THEN
1029  cntr = cntr + 1
1030  END IF ! iloc
1031  END DO ! ivl
1032 
1033 ! --------------- Need faces which share 2 vertices with face 1 ----------------
1034 
1035  IF ( cntr == 2 ) THEN ! Have proper candidate face
1036  celldegr(face_type_quad,c1) = celldegr(face_type_quad,c1) &
1037  + 1
1038 
1039  v1 = gridcobalt%f2v(1,ifg)
1040  v2 = gridcobalt%f2v(2,ifg)
1041  v3 = gridcobalt%f2v(3,ifg)
1042  v4 = gridcobalt%f2v(4,ifg)
1043 
1044  xyznodes(xcoord:zcoord,1) = pgrid%xyz(xcoord:zcoord,v1)
1045  xyznodes(xcoord:zcoord,2) = pgrid%xyz(xcoord:zcoord,v2)
1046  xyznodes(xcoord:zcoord,3) = pgrid%xyz(xcoord:zcoord,v3)
1047  xyznodes(xcoord:zcoord,4) = pgrid%xyz(xcoord:zcoord,v4)
1048 
1049  CALL facevectorquad(xyznodes(xcoord:zcoord,1:4), &
1050  fvecx,fvecy,fvecz)
1051  CALL facecentroidquad(xyznodes(xcoord:zcoord,1:4), &
1052  fcenx,fceny,fcenz)
1053 
1054 ! ----------------- Check orientation, reverse if necessary --------------------
1055 
1056  dotprod = fvecx*(fcenx - pgrid%cofgApp(xcoord,c1)) &
1057  + fvecy*(fceny - pgrid%cofgApp(ycoord,c1)) &
1058  + fvecz*(fcenz - pgrid%cofgApp(zcoord,c1))
1059 
1060  IF ( dotprod < 0.0_rfreal ) THEN
1061  vert(1) = v1
1062  vert(2) = v4
1063  vert(3) = v3
1064  vert(4) = v2
1065  ELSE
1066  vert(1) = v1
1067  vert(2) = v2
1068  vert(3) = v3
1069  vert(4) = v4
1070  END IF ! dotProd
1071 
1072 ! ----------------- Determine local face number --------------------------------
1073 
1074  verttemp(1:3) = pgrid%pri2v(1:3,icl)
1075  CALL quicksortinteger(verttemp(1:3),3)
1076 
1077  DO ivl = 1,3
1078  ivg = pgrid%pri2v(ivl,icl)
1079 
1080  CALL binarysearchinteger(verttemp(1:3),3,ivg,iloc)
1081 
1082  IF ( iloc /= element_not_found ) THEN
1083  key(iloc) = ivl
1084  ELSE
1085  CALL errorstop(global,err_reached_default,__line__)
1086  END IF ! iloc
1087  END DO ! ivl
1088 
1089  cntr = 0
1090  flag(1:3) = 0
1091 
1092  DO ivl = 1,4
1093  ivg = vert(ivl)
1094 
1095  CALL binarysearchinteger(verttemp(1:3),3,ivg,iloc)
1096 
1097  IF ( iloc /= element_not_found ) THEN
1098  flag(key(iloc)) = 1
1099  ELSE
1100  cntr = cntr + 1
1101  END IF ! iloc
1102  END DO ! ivl
1103 
1104 ! ----------------- Cycle vertices of face 5 to match ordering on face 1 -------
1105 
1106  IF ( flag(1) == 1 .AND. flag(2) == 1 ) THEN ! Face 2
1107  CALL cyclelist(vert(1:4),4,1,pgrid%pri2v(1,icl))
1108  CALL cyclelist(pgrid%pri2v(4:6,icl),3,1,vert(4))
1109  ELSE IF ( flag(2) == 1 .AND. flag(3) == 1 ) THEN ! Face 3
1110  CALL cyclelist(vert(1:4),4,1,pgrid%pri2v(2,icl))
1111  CALL cyclelist(pgrid%pri2v(4:6,icl),3,2,vert(4))
1112  ELSE IF ( flag(3) == 1 .AND. flag(1) == 1 ) THEN ! Face 4
1113  CALL cyclelist(vert(1:4),4,1,pgrid%pri2v(3,icl))
1114  CALL cyclelist(pgrid%pri2v(4:6,icl),3,3,vert(4))
1115  ELSE
1116  CALL errorstop(global,err_reached_default,__line__)
1117  END IF ! flag
1118  ELSE IF ( cntr /= 3 .AND. cntr /= 0 ) THEN
1119  CALL errorstop(global,err_reached_default,__line__)
1120  END IF ! cntr
1121  END IF ! cellDegr
1122 
1123  END SELECT ! pGrid%cellGlob2Loc
1124  END IF ! gridCOBALT%nvpf
1125  END IF ! c1
1126 
1127  END DO ! j
1128  END DO ! ifg
1129 
1130  IF ( global%verbLevel > verbose_none ) THEN
1131  WRITE(stdout,'(A,3X,A)') solver_name,'Checking cell connectivity done.'
1132  END IF ! global%verbLevel
1133 
1134 ! ==============================================================================
1135 ! Deallocate remaining temporary memory
1136 ! ==============================================================================
1137 
1138  DEALLOCATE(celldegr,stat=errorflag)
1139  global%error = errorflag
1140  IF ( global%error /= err_none ) THEN
1141  CALL errorstop(global,err_deallocate,__line__,'cellDegr')
1142  END IF ! global%error
1143 
1144 ! ******************************************************************************
1145 ! Building boundary data structure. Cobalt has boundary flags in the face
1146 ! lists as negative cell indices. The boundary flags do not have to be
1147 ! contiguous. Use mapping file to map the boundary flags to Rocflu patches.
1148 ! ******************************************************************************
1149 
1150  IF ( global%verbLevel > verbose_none ) THEN
1151  WRITE(stdout,'(A,3X,A)') solver_name,'Building boundary face lists...'
1152  END IF ! global%verbLevel
1153 
1154 ! ==============================================================================
1155 ! Read additional info on grid, required for mapping of patches
1156 ! ==============================================================================
1157 
1158  IF ( global%verbLevel > verbose_none ) THEN
1159  WRITE(stdout,'(A,5X,A)') solver_name,'Reading information file...'
1160  END IF ! global%verbLevel
1161 
1162 ! ------------------------------------------------------------------------------
1163 ! Open file
1164 ! ------------------------------------------------------------------------------
1165 
1166  ifile = if_grid
1167 
1168  CALL buildfilenameplain(global,filedest_indir,'.cgi',ifilename)
1169 
1170  OPEN(ifile,file=ifilename,form="FORMATTED",status="OLD",iostat=errorflag)
1171  global%error = errorflag
1172  IF ( global%error /= err_none ) THEN
1173  CALL errorstop(global,err_file_open,__line__,ifilename)
1174  END IF ! global%error
1175 
1176 ! ------------------------------------------------------------------------------
1177 ! Read file
1178 ! ------------------------------------------------------------------------------
1179 
1180  READ(ifile,*) pgrid%nPatches
1181  READ(ifile,*) gridcobalt%nMappings
1182 
1183  ALLOCATE(gridcobalt%patch2bc(3,gridcobalt%nMappings),stat=errorflag)
1184  global%error = errorflag
1185  IF ( global%error /= err_none ) THEN
1186  CALL errorstop(global,err_allocate,__line__,'gridCOBALT%patch2bc')
1187  END IF ! global%error
1188 
1189  DO imap = 1,gridcobalt%nMappings
1190  READ(ifile,*) (gridcobalt%patch2bc(j,imap),j=1,3)
1191  END DO ! iMap
1192 
1193 ! ------------------------------------------------------------------------------
1194 ! Close file
1195 ! ------------------------------------------------------------------------------
1196 
1197  CLOSE(ifile,iostat=errorflag)
1198  global%error = errorflag
1199  IF ( global%error /= err_none ) THEN
1200  CALL errorstop(global,err_file_close,__line__,ifilename)
1201  END IF ! global%error
1202 
1203  IF ( global%verbLevel > verbose_none ) THEN
1204  WRITE(stdout,'(A,5X,A)') solver_name,'Reading information file done.'
1205  END IF ! global%verbLevel
1206 
1207 ! ==============================================================================
1208 ! Check for consistent input - somewhat complicated...
1209 ! ==============================================================================
1210 
1211  IF ( global%checkLevel > check_none ) THEN
1212  IF ( global%verbLevel > verbose_none ) THEN
1213  WRITE(stdout,'(A,5X,A)') solver_name,'Checking patch mapping entries...'
1214  END IF ! global%verbLevel
1215 
1216  DO imap = 1,gridcobalt%nMappings
1217  IF ( gridcobalt%patch2bc(2,imap) < gridcobalt%patch2bc(1,imap) ) THEN
1218  IF ( global%verbLevel > verbose_none ) THEN
1219  WRITE(stdout,'(A,7X,A)') solver_name,'Check failed.'
1220  END IF ! global%verbLevel
1221  CALL errorstop(global,err_patch_numbering,__line__)
1222  END IF ! gridCOBALT
1223  END DO ! iMap
1224 
1225  IF ( minval(gridcobalt%patch2bc(3,:)) /= 1 .OR. &
1226  maxval(gridcobalt%patch2bc(3,:)) /= pgrid%nPatches ) THEN
1227  IF ( global%verbLevel > verbose_none ) THEN
1228  WRITE(stdout,'(A,7X,A)') solver_name,'Check failed.'
1229  END IF ! global%verbLevel
1230  CALL errorstop(global,err_patch_numbering,__line__)
1231  END IF ! gridCOBALT
1232 
1233  DO imap = 1,gridcobalt%nMappings
1234  DO imap2 = 1,gridcobalt%nMappings
1235 
1236  IF ( imap /= imap2 ) THEN
1237  ibeg1 = gridcobalt%patch2bc(1,imap)
1238  iend1 = gridcobalt%patch2bc(2,imap)
1239 
1240  ibeg2 = gridcobalt%patch2bc(1,imap2)
1241  iend2 = gridcobalt%patch2bc(2,imap2)
1242 
1243  IF ( ibeg1 < ibeg2 ) THEN
1244  ibegmin = ibeg1
1245  iendmin = iend1
1246  ibegmax = ibeg2
1247  iendmax = iend2
1248  ELSE IF ( ibeg1 > ibeg2 ) THEN
1249  ibegmin = ibeg2
1250  iendmin = iend2
1251  ibegmax = ibeg1
1252  iendmax = iend1
1253  ELSE ! iBeg1 and iBeg2 have the same value
1254  IF ( global%verbLevel > verbose_none ) THEN
1255  WRITE(stdout,'(A,7X,A)') solver_name,'Check failed.'
1256  END IF ! global%verbLevel
1257  CALL errorstop(global,err_patch_numbering,__line__)
1258  END IF ! iBeg1
1259 
1260  IF ( iendmin >= ibegmax ) THEN
1261  IF ( global%verbLevel > verbose_none ) THEN
1262  WRITE(stdout,'(A,7X,A)') solver_name,'Check failed.'
1263  END IF ! global%verbLevel
1264  CALL errorstop(global,err_patch_numbering,__line__)
1265  END IF ! iebmax
1266  END IF ! iMap
1267 
1268  END DO ! iMap2
1269  END DO ! iMap
1270 
1271  IF ( global%verbLevel > verbose_none ) THEN
1272  WRITE(stdout,'(A,5X,A)') solver_name, &
1273  'Checking patch mapping entries done.'
1274  END IF ! global%verbLevel
1275  END IF ! global%checkLevel
1276 
1277 ! ******************************************************************************
1278 ! Convert boundary information to ROCFLU format
1279 ! ******************************************************************************
1280 
1281 ! ==============================================================================
1282 ! Allocate patch memory and initialize patch structure
1283 ! ==============================================================================
1284 
1285  ALLOCATE(pregion%patches(pgrid%nPatches),stat=errorflag)
1286  global%error = errorflag
1287  IF ( global%error /= err_none ) THEN
1288  CALL errorstop(global,err_allocate,__line__,'pRegion%patches')
1289  END IF ! global%error
1290 
1291  DO ipatch = 1,pgrid%nPatches
1292  ppatch => pregion%patches(ipatch)
1293 
1294  ppatch%nBTris = 0
1295  ppatch%nBQuads = 0
1296  ppatch%nBVert = 0
1297 
1298  ppatch%iPatchGlobal = ipatch
1299  ppatch%iBorder = patch_iborder_default
1300  ppatch%renumFlag = .false.
1301  END DO ! iPatch
1302 
1303  global%nPatches = pgrid%nPatches
1304 
1305 ! ==============================================================================
1306 ! Determine number of triangular and quadrilateral faces on each patch
1307 ! ==============================================================================
1308 
1309  DO ifg = 1,gridcobalt%nFaces
1310  DO j = 1,2
1311  c1 = gridcobalt%f2c(j,ifg)
1312 
1313  IF ( c1 < 0 ) THEN ! boundary cell
1314  nullify(ppatch)
1315 
1316  DO imap = 1,gridcobalt%nMappings
1317  IF ( abs(c1) >= gridcobalt%patch2bc(1,imap) .AND. &
1318  abs(c1) <= gridcobalt%patch2bc(2,imap) ) THEN
1319  ipatch = gridcobalt%patch2bc(3,imap)
1320  ppatch => pregion%patches(ipatch)
1321  EXIT
1322  END IF ! ABS
1323  END DO ! iMap
1324 
1325  IF ( ASSOCIATED(ppatch) .EQV. .true. ) THEN
1326  IF ( gridcobalt%nvpf(ifg) == 3 ) THEN
1327  ppatch%nBTris = ppatch%nBTris + 1
1328  ELSE IF ( gridcobalt%nvpf(ifg) == 4 ) THEN
1329  ppatch%nBQuads = ppatch%nBQuads + 1
1330  ELSE ! Should be impossible, keep for defensive programming
1331  CALL errorstop(global,err_reached_default,__line__)
1332  END IF ! gridCOBALT%nvpf
1333  ELSE
1334  CALL errorstop(global,err_reached_default,__line__)
1335  END IF ! ASSOCIATED
1336  END IF ! c1
1337  END DO ! j
1338  END DO ! ifg
1339 
1340 ! ==============================================================================
1341 ! Set total boundary patch quantities and number of boundary faces
1342 ! ==============================================================================
1343 
1344  pgrid%nBFaces = 0
1345 
1346  DO ipatch = 1,pgrid%nPatches
1347  ppatch => pregion%patches(ipatch)
1348 
1349  ppatch%nBFaces = ppatch%nBTris + ppatch%nBQuads
1350  pgrid%nBFaces = pgrid%nBFaces + ppatch%nBFaces
1351 
1352  ppatch%nBFacesTot = ppatch%nBFaces
1353  ppatch%nBQuadsTot = ppatch%nBQuads
1354  ppatch%nBTrisTot = ppatch%nBTris
1355  ppatch%nBVertTot = ppatch%nBVert
1356 
1357  ppatch%nBTrisMax = rflu_setmaxdimension(global,ppatch%nBTrisTot)
1358  ppatch%nBQuadsMax = rflu_setmaxdimension(global,ppatch%nBQuadsTot)
1359  ppatch%nBFacesMax = rflu_setmaxdimension(global,ppatch%nBFacesTot)
1360  ppatch%nBVertMax = rflu_setmaxdimension(global,ppatch%nBVertTot)
1361 
1362  ppatch%nBCellsVirt = 0
1363  END DO ! iPatch
1364 
1365  pgrid%nBFacesTot = pgrid%nBFaces
1366 
1367 ! ==============================================================================
1368 ! Allocate memory for boundary face lists
1369 ! ==============================================================================
1370 
1371  DO ipatch = 1,pgrid%nPatches
1372  ppatch => pregion%patches(ipatch)
1373 
1374  IF ( ppatch%nBTrisMax > 0 ) THEN
1375  ALLOCATE(ppatch%bTri2v(3,ppatch%nBTrisMax),stat=errorflag)
1376  global%error = errorflag
1377  IF ( global%error /= err_none ) THEN
1378  CALL errorstop(global,err_allocate,__line__,'pPatch%bTri2v')
1379  END IF ! global%error
1380 
1381  ppatch%nBTris = 0 ! Reused as counter below
1382  END IF ! pPatch%nBTrisMax
1383 
1384  IF ( ppatch%nBQuadsMax > 0 ) THEN
1385  ALLOCATE(ppatch%bQuad2v(4,ppatch%nBQuadsMax),stat=errorflag)
1386  global%error = errorflag
1387  IF ( global%error /= err_none ) THEN
1388  CALL errorstop(global,err_allocate,__line__,'pPatch%bQuad2v')
1389  END IF ! global%error
1390 
1391  ppatch%nBQuads = 0 ! Reused as counter below
1392  END IF ! pPatch%nBQuadsMax
1393  END DO ! iPatch
1394 
1395 ! ==============================================================================
1396 ! Extract boundary faces and check orientation (best done here because have
1397 ! access to other cell index)
1398 ! ==============================================================================
1399 
1400  DO ifg = 1,gridcobalt%nFaces
1401  DO j = 1,2
1402  c1 = gridcobalt%f2c(j,ifg)
1403 
1404  IF ( c1 < 0 ) THEN ! boundary cell
1405  DO imap = 1,gridcobalt%nMappings
1406  IF ( abs(c1) >= gridcobalt%patch2bc(1,imap) .AND. &
1407  abs(c1) <= gridcobalt%patch2bc(2,imap) ) THEN
1408  ipatch = gridcobalt%patch2bc(3,imap)
1409  ppatch => pregion%patches(ipatch)
1410  EXIT
1411  END IF ! ABS
1412  END DO ! iMap
1413 
1414  IF ( ASSOCIATED(ppatch) .EQV. .true. ) THEN
1415  c2 = gridcobalt%f2c(3-j,ifg)
1416 
1417  IF ( gridcobalt%nvpf(ifg) == 3 ) THEN
1418  ppatch%nBTris = ppatch%nBTris + 1
1419 
1420  v1 = gridcobalt%f2v(1,ifg)
1421  v2 = gridcobalt%f2v(2,ifg)
1422  v3 = gridcobalt%f2v(3,ifg)
1423 
1424  xyznodes(xcoord:zcoord,1) = pgrid%xyz(xcoord:zcoord,v1)
1425  xyznodes(xcoord:zcoord,2) = pgrid%xyz(xcoord:zcoord,v2)
1426  xyznodes(xcoord:zcoord,3) = pgrid%xyz(xcoord:zcoord,v3)
1427 
1428  CALL facevectortria(xyznodes(xcoord:zcoord,1:3), &
1429  fvecx,fvecy,fvecz)
1430  CALL facecentroidtria(xyznodes(xcoord:zcoord,1:3), &
1431  fcenx,fceny,fcenz)
1432 
1433  dotprod = fvecx*(fcenx - pgrid%cofgApp(xcoord,c2)) &
1434  + fvecy*(fceny - pgrid%cofgApp(ycoord,c2)) &
1435  + fvecz*(fcenz - pgrid%cofgApp(zcoord,c2))
1436 
1437  IF ( dotprod < 0.0_rfreal ) THEN
1438  ppatch%bTri2v(1,ppatch%nBTris) = v1
1439  ppatch%bTri2v(2,ppatch%nBTris) = v3
1440  ppatch%bTri2v(3,ppatch%nBTris) = v2
1441  ELSE
1442  ppatch%bTri2v(1,ppatch%nBTris) = v1
1443  ppatch%bTri2v(2,ppatch%nBTris) = v2
1444  ppatch%bTri2v(3,ppatch%nBTris) = v3
1445  END IF ! dotProd
1446  ELSE IF ( gridcobalt%nvpf(ifg) == 4 ) THEN
1447  ppatch%nBQuads = ppatch%nBQuads + 1
1448 
1449  v1 = gridcobalt%f2v(1,ifg)
1450  v2 = gridcobalt%f2v(2,ifg)
1451  v3 = gridcobalt%f2v(3,ifg)
1452  v4 = gridcobalt%f2v(4,ifg)
1453 
1454  xyznodes(xcoord:zcoord,1) = pgrid%xyz(xcoord:zcoord,v1)
1455  xyznodes(xcoord:zcoord,2) = pgrid%xyz(xcoord:zcoord,v2)
1456  xyznodes(xcoord:zcoord,3) = pgrid%xyz(xcoord:zcoord,v3)
1457  xyznodes(xcoord:zcoord,4) = pgrid%xyz(xcoord:zcoord,v4)
1458 
1459  CALL facevectorquad(xyznodes(xcoord:zcoord,1:4), &
1460  fvecx,fvecy,fvecz)
1461  CALL facecentroidquad(xyznodes(xcoord:zcoord,1:4), &
1462  fcenx,fceny,fcenz)
1463 
1464  dotprod = fvecx*(fcenx - pgrid%cofgApp(xcoord,c2)) &
1465  + fvecy*(fceny - pgrid%cofgApp(ycoord,c2)) &
1466  + fvecz*(fcenz - pgrid%cofgApp(zcoord,c2))
1467 
1468  IF ( dotprod < 0.0_rfreal ) THEN
1469  ppatch%bQuad2v(1,ppatch%nBQuads) = v1
1470  ppatch%bQuad2v(2,ppatch%nBQuads) = v4
1471  ppatch%bQuad2v(3,ppatch%nBQuads) = v3
1472  ppatch%bQuad2v(4,ppatch%nBQuads) = v2
1473  ELSE
1474  ppatch%bQuad2v(1,ppatch%nBQuads) = v1
1475  ppatch%bQuad2v(2,ppatch%nBQuads) = v2
1476  ppatch%bQuad2v(3,ppatch%nBQuads) = v3
1477  ppatch%bQuad2v(4,ppatch%nBQuads) = v4
1478  END IF ! dotProd
1479  ELSE ! Should be impossible, keep for defensive programming
1480  CALL errorstop(global,err_reached_default,__line__)
1481  END IF ! gridCOBALT%nvpf
1482  ELSE
1483  CALL errorstop(global,err_reached_default,__line__)
1484  END IF ! ASSOCIATED
1485  END IF ! c1
1486 
1487  END DO ! j
1488  END DO ! ifg
1489 
1490  IF ( global%verbLevel > verbose_none ) THEN
1491  WRITE(stdout,'(A,3X,A)') solver_name,'Building boundary face lists done.'
1492  END IF ! global%verbLevel
1493 
1494 ! ==============================================================================
1495 ! Deallocate temporary memory
1496 ! ==============================================================================
1497 
1498  DEALLOCATE(gridcobalt%patch2bc,stat=errorflag)
1499  global%error = errorflag
1500  IF ( global%error /= err_none ) THEN
1501  CALL errorstop(global,err_deallocate,__line__,'gridCOBALT%patch2bc')
1502  END IF ! global%error
1503 
1504 ! ==============================================================================
1505 ! Destroy approximate centroids and cell mapping.
1506 ! ==============================================================================
1507 
1508  CALL rflu_destroyapproxcentroids(pregion)
1509  CALL rflu_destroycellmapping(pregion)
1510 
1511 ! *****************************************************************************
1512 ! Allocate memory for boundary face lists bf2c and bf2v
1513 ! *****************************************************************************
1514 
1515  DO ipatch = 1,pgrid%nPatches
1516  ppatch => pregion%patches(ipatch)
1517 
1518  ALLOCATE(ppatch%bf2c(ppatch%nBFacesMax),stat=errorflag)
1519  global%error = errorflag
1520  IF ( global%error /= err_none ) THEN
1521  CALL errorstop(global,err_allocate,__line__,'pPatch%bf2c')
1522  END IF ! global%error
1523 
1524  ALLOCATE(ppatch%bf2v(4,ppatch%nBFacesMax),stat=errorflag)
1525  global%error = errorflag
1526  IF ( global%error /= err_none ) THEN
1527  CALL errorstop(global,err_allocate,__line__,'pPatch%bf2v')
1528  END IF ! global%error
1529 
1530  DO ifl = 1,ppatch%nBFacesMax
1531  ppatch%bf2v(1,ifl) = vert_none
1532  ppatch%bf2v(2,ifl) = vert_none
1533  ppatch%bf2v(3,ifl) = vert_none
1534  ppatch%bf2v(4,ifl) = vert_none
1535  END DO ! ifl
1536  END DO ! iPatch
1537 
1538 ! ******************************************************************************
1539 ! End
1540 ! ******************************************************************************
1541 
1542  IF ( global%verbLevel > verbose_none ) THEN
1543  WRITE(stdout,'(A,1X,A)') solver_name, &
1544  'Converting from COBALT to ROCFLU format done.'
1545  END IF ! global%verbLevel
1546 
1547  CALL deregisterfunction(global)
1548 
1549  END SUBROUTINE rflu_convcobalt2rocflu
1550 
1551 
1552 
1553 
1554 
1555 
1556 
1557 
1558 
1559 
1560 ! ******************************************************************************
1561 !
1562 ! Purpose: Read grid file from COBALT in ASCII format.
1563 !
1564 ! Description: None.
1565 !
1566 ! Input:
1567 ! pRegion Pointer to region
1568 !
1569 ! Output: None.
1570 !
1571 ! Notes:
1572 ! 1. COBALT grid file is usually called <something>.inp, which conflicts
1573 ! with the Rocflu input file, so the COBALT grid file needs to be renamed
1574 ! before being run through rfluprep. The same applies to the COBALT
1575 ! boundary condition file, which is called <something>.bc.
1576 ! 2. Part of the file is read directly into the grid data structure, while
1577 ! the remainder needs to be processed in RFLU_ConvCOBALT2ROCFLU.F90.
1578 !
1579 ! ******************************************************************************
1580 
1581  SUBROUTINE rflu_readgridcobalt(pRegion)
1582 
1583  IMPLICIT NONE
1584 
1585 ! ******************************************************************************
1586 ! Declarations and definitions
1587 ! ******************************************************************************
1588 
1589 ! ==============================================================================
1590 ! Arguments
1591 ! ==============================================================================
1592 
1593  TYPE(t_region), POINTER :: pregion
1594 
1595 ! ==============================================================================
1596 ! Local variables
1597 ! ==============================================================================
1598 
1599  CHARACTER(CHRLEN) :: dummystring,ifilename
1600  INTEGER :: errorflag,ic,ifc,ifile,iv,j,ndim,nvertperface,nzones
1601  TYPE(t_grid), POINTER :: pgrid
1602  TYPE(t_global), POINTER :: global
1603 
1604 ! ******************************************************************************
1605 ! Start
1606 ! ******************************************************************************
1607 
1608  global => pregion%global
1609 
1610  CALL registerfunction(global,'RFLU_ReadGridCOBALT', &
1611  'RFLU_ModCOBALT.F90')
1612 
1613  IF ( global%verbLevel > verbose_none ) THEN
1614  WRITE(stdout,'(A,1X,A)') solver_name,'Reading COBALT grid file...'
1615  END IF ! global%verbLevel
1616 
1617 ! ==============================================================================
1618 ! Open grid file
1619 ! ==============================================================================
1620 
1621  ifile = if_grid
1622 
1623  CALL buildfilenameplain(global,filedest_indir,'.cgr',ifilename)
1624 
1625  OPEN(ifile,file=ifilename,form="FORMATTED",status="OLD",iostat=errorflag)
1626  global%error = errorflag
1627  IF ( global%error /= err_none ) THEN
1628  CALL errorstop(global,err_file_open,__line__,ifilename)
1629  END IF ! global%error
1630 
1631 ! ==============================================================================
1632 ! Header
1633 ! ==============================================================================
1634 
1635  IF ( global%verbLevel > verbose_none ) THEN
1636  WRITE(stdout,'(A,3X,A)') solver_name,'Header...'
1637  END IF ! global%verbLevel
1638 
1639  pgrid => pregion%grid
1640 
1641  READ(ifile,*) ndim,nzones,gridcobalt%nPatches
1642 
1643  IF ( ndim /= 3 ) THEN
1644  CALL errorstop(global,err_ndimens_invalid,__line__)
1645  END IF ! nDim
1646 
1647  IF ( nzones /= 1 ) THEN
1648  CALL errorstop(global,err_nzones_invalid,__line__)
1649  END IF ! nZones
1650 
1651  READ(ifile,*) pgrid%nVertTot,gridcobalt%nFaces,pgrid%nCellsTot, &
1652  gridcobalt%nVertPerFaceMax,gridcobalt%nFacesPerCellMax
1653 
1654  pgrid%nVert = pgrid%nVertTot
1655  pgrid%nCells = pgrid%nCellsTot
1656 
1657  pgrid%nVertMax = rflu_setmaxdimension(global,pgrid%nVertTot)
1658  pgrid%nCellsMax = rflu_setmaxdimension(global,pgrid%nCellsTot)
1659 
1660 ! ==============================================================================
1661 ! Coordinates
1662 ! ==============================================================================
1663 
1664  IF ( global%verbLevel > verbose_none ) THEN
1665  WRITE(stdout,'(A,3X,A)') solver_name,'Coordinates...'
1666  END IF ! global%verbLevel
1667 
1668  ALLOCATE(pgrid%xyz(3,pgrid%nVertMax),stat=errorflag)
1669  global%error = errorflag
1670  IF ( global%error /= err_none ) THEN
1671  CALL errorstop(global,err_allocate,__line__,'pGrid%xyz')
1672  END IF ! global%error
1673 
1674  DO iv = 1,pgrid%nVertTot
1675  READ(ifile,*) (pgrid%xyz(j,iv),j=1,3)
1676  END DO ! iv
1677 
1678 ! ==============================================================================
1679 ! Connectivity
1680 ! ==============================================================================
1681 
1682  IF ( global%verbLevel > verbose_none ) THEN
1683  WRITE(stdout,'(A,3X,A)') solver_name,'Connectivity...'
1684  END IF ! global%verbLevel
1685 
1686  ALLOCATE(gridcobalt%f2v(gridcobalt%nVertPerFaceMax,gridcobalt%nFaces), &
1687  stat=errorflag)
1688  global%error = errorflag
1689  IF ( global%error /= err_none ) THEN
1690  CALL errorstop(global,err_allocate,__line__,'gridCOBALT%f2v')
1691  END IF ! global%error
1692 
1693  DO ifc = 1,gridcobalt%nFaces
1694  DO j = 1,gridcobalt%nVertPerFaceMax
1695  gridcobalt%f2v(j,ifc) = 0
1696  END DO ! j
1697  END DO ! ifc
1698 
1699  ALLOCATE(gridcobalt%f2c(2,gridcobalt%nFaces),stat=errorflag)
1700  global%error = errorflag
1701  IF ( global%error /= err_none ) THEN
1702  CALL errorstop(global,err_allocate,__line__,'gridCOBALT%f2c')
1703  END IF ! errorFlag
1704 
1705  ALLOCATE(gridcobalt%nvpf(gridcobalt%nFaces),stat=errorflag)
1706  global%error = errorflag
1707  IF ( global%error /= err_none ) THEN
1708  CALL errorstop(global,err_allocate,__line__,'gridCOBALT%nvpf')
1709  END IF ! errorFlag
1710 
1711  gridcobalt%nTris = 0
1712  gridcobalt%nQuads = 0
1713 
1714  DO ifc = 1,gridcobalt%nFaces
1715  READ(ifile,*) nvertperface
1716  backspace(ifile)
1717 
1718  IF ( nvertperface == 3 ) THEN
1719  gridcobalt%nTris = gridcobalt%nTris + 1
1720  ELSE IF ( nvertperface == 4 ) THEN
1721  gridcobalt%nQuads = gridcobalt%nQuads + 1
1722  ELSE
1723  CALL errorstop(global,err_facetype_invalid,__line__)
1724  END IF ! nVertPerFace
1725 
1726  READ(ifile,*) gridcobalt%nvpf(ifc), &
1727  (gridcobalt%f2v(j,ifc),j=1,nvertperface), &
1728  (gridcobalt%f2c(j,ifc),j=1,2)
1729  END DO ! ifc
1730 
1731 ! ==============================================================================
1732 ! Close grid file
1733 ! ==============================================================================
1734 
1735  CLOSE(ifile,iostat=errorflag)
1736  global%error = errorflag
1737  IF ( global%error /= err_none ) THEN
1738  CALL errorstop(global,err_file_close,__line__,ifilename)
1739  END IF ! global%error
1740 
1741 ! ******************************************************************************
1742 ! Print grid statistics
1743 ! ******************************************************************************
1744 
1745  IF ( global%verbLevel > verbose_none ) THEN
1746  WRITE(stdout,'(A,3X,A)') solver_name,'Grid Statistics:'
1747  WRITE(stdout,'(A,5X,A,11X,I9)') solver_name,'Vertices: ', &
1748  pgrid%nVertTot
1749  WRITE(stdout,'(A,5X,A,11X,I9)') solver_name,'Cells: ', &
1750  pgrid%nCellsTot
1751  WRITE(stdout,'(A,5X,A,11X,I9)') solver_name,'Faces: ', &
1752  gridcobalt%nFaces
1753  WRITE(stdout,'(A,7X,A,9X,I9)') solver_name,'Triangular faces: ', &
1754  gridcobalt%nTris
1755  WRITE(stdout,'(A,7X,A,9X,I9)') solver_name,'Quadrilateral faces:', &
1756  gridcobalt%nQuads
1757  WRITE(stdout,'(A,5X,A,11X,I9)') solver_name,'Patches: ', &
1758  gridcobalt%nPatches
1759  END IF ! global%verbLevel
1760 
1761 ! ******************************************************************************
1762 ! End
1763 ! ******************************************************************************
1764 
1765  IF ( global%verbLevel > verbose_none ) THEN
1766  WRITE(stdout,'(A,1X,A)') solver_name,'Reading COBALT grid file done.'
1767  END IF ! global%verbLevel
1768 
1769  CALL deregisterfunction(global)
1770 
1771  END SUBROUTINE rflu_readgridcobalt
1772 
1773 
1774 
1775 
1776 
1777 
1778 
1779 
1780 
1781 ! ******************************************************************************
1782 ! End
1783 ! ******************************************************************************
1784 
1785 
1786 END MODULE rflu_modcobalt
1787 
1788 ! ******************************************************************************
1789 !
1790 ! RCS Revision history:
1791 !
1792 ! $Log: RFLU_ModCOBALT.F90,v $
1793 ! Revision 1.5 2008/12/06 08:45:03 mtcampbe
1794 ! Updated license.
1795 !
1796 ! Revision 1.4 2008/11/19 22:18:14 mtcampbe
1797 ! Added Illinois Open Source License/Copyright
1798 !
1799 ! Revision 1.3 2006/04/28 19:48:58 haselbac
1800 ! Bug fix: Cell connectivity arrays were only dimensioned to Tot, not Max
1801 !
1802 ! Revision 1.2 2006/03/25 22:04:29 haselbac
1803 ! Changes because of sype patches
1804 !
1805 ! Revision 1.1 2005/04/15 15:09:08 haselbac
1806 ! Initial revision
1807 !
1808 ! Revision 1.7 2005/01/20 14:54:56 haselbac
1809 ! Added setting of nBFaces and nBFacesTot
1810 !
1811 ! Revision 1.6 2004/12/10 15:22:14 haselbac
1812 ! Bug fix: Added missing RFLU_SetMaxDimensions call
1813 !
1814 ! Revision 1.5 2004/12/04 03:37:35 haselbac
1815 ! Adapted to changes in RFLU_ModCellMapping
1816 !
1817 ! Revision 1.4 2004/11/03 17:09:24 haselbac
1818 ! Removed setting of vertex and cell flags
1819 !
1820 ! Revision 1.3 2004/11/03 15:06:08 haselbac
1821 ! Bug fix: Corrected loop limits
1822 !
1823 ! Revision 1.2 2004/10/19 19:31:08 haselbac
1824 ! Removed renumbering of bface lists
1825 !
1826 ! Revision 1.1 2004/07/06 15:15:47 haselbac
1827 ! Initial revision
1828 !
1829 ! ******************************************************************************
1830 
1831 
1832 
1833 
1834 
1835 
1836 
1837 
rflu_modgrid
Definition: RFLU_ModGrid.F90:39
facevectortria
subroutine facevectortria(xyzNodes, fVecX, fVecY, fVecZ)
Definition: FaceVector.F90:49
rflu_modcobalt
Definition: RFLU_ModCOBALT.F90:39
rflu_modgeometry::rflu_destroyapproxcentroids
subroutine, public rflu_destroyapproxcentroids(pRegion)
Definition: RFLU_ModGeometry.F90:1960
rflu_modcobalt::t_gridcobalt
Definition: RFLU_ModCOBALT.F90:69
moderror::registerfunction
subroutine registerfunction(global, funName, fileName)
Definition: ModError.F90:449
Attribute::status
int status() const
Obtain the status of the attribute.
Definition: Attribute.h:240
ic
**********************************************************************Rocstar Simulation Suite Illinois Rocstar LLC All rights reserved ****Illinois Rocstar LLC IL **www illinoisrocstar com **sales illinoisrocstar com WITHOUT WARRANTY OF ANY **EXPRESS OR INCLUDING BUT NOT LIMITED TO THE WARRANTIES **OF FITNESS FOR A PARTICULAR PURPOSE AND **NONINFRINGEMENT IN NO EVENT SHALL THE CONTRIBUTORS OR **COPYRIGHT HOLDERS BE LIABLE FOR ANY DAMAGES OR OTHER WHETHER IN AN ACTION OF TORT OR **Arising OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE **USE OR OTHER DEALINGS WITH THE SOFTWARE **********************************************************************INTERFACE SUBROUTINE ic
Definition: RFLO_FindSourceCell.h:27
rflu_moddimensions::rflu_setmaxdimensions
subroutine, public rflu_setmaxdimensions(pRegion)
Definition: RFLU_ModDimensions.F90:613
modgrid::t_grid
Definition: ModGrid.F90:56
modsortsearch::cyclelist
subroutine cyclelist(a, na, i, v)
Definition: ModSortSearch.F90:188
modsortsearch::quicksortinteger
subroutine quicksortinteger(a, n)
Definition: ModSortSearch.F90:482
rflu_modcellmapping
Definition: RFLU_ModCellMapping.F90:39
modsortsearch::binarysearchinteger
subroutine binarysearchinteger(a, n, v, i, j)
Definition: ModSortSearch.F90:98
buildfilenameplain
subroutine buildfilenameplain(global, dest, ext, fileName)
Definition: BuildFileNamePlain.F90:48
rflu_modcobalt::rflu_convcobalt2rocflu
subroutine, public rflu_convcobalt2rocflu(pRegion)
Definition: RFLU_ModCOBALT.F90:115
rflu_modcobalt::rflu_readgridcobalt
subroutine, public rflu_readgridcobalt(pRegion)
Definition: RFLU_ModCOBALT.F90:1581
rflu_modcellmapping::rflu_buildglob2loccellmapping
subroutine, public rflu_buildglob2loccellmapping(pRegion)
Definition: RFLU_ModCellMapping.F90:307
form
**********************************************************************Rocstar Simulation Suite Illinois Rocstar LLC All rights reserved ****Illinois Rocstar LLC IL **www illinoisrocstar com **sales illinoisrocstar com WITHOUT WARRANTY OF ANY **EXPRESS OR INCLUDING BUT NOT LIMITED TO THE WARRANTIES **OF FITNESS FOR A PARTICULAR PURPOSE AND **NONINFRINGEMENT IN NO EVENT SHALL THE CONTRIBUTORS OR **COPYRIGHT HOLDERS BE LIABLE FOR ANY DAMAGES OR OTHER WHETHER IN AN ACTION OF TORT OR **Arising OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE **USE OR OTHER DEALINGS WITH THE SOFTWARE **********************************************************************INTERFACE SUBROUTINE form
Definition: RFLO_ReadDataFileInt.h:26
rflu_modcellmapping::rflu_destroycellmapping
subroutine, public rflu_destroycellmapping(pRegion)
Definition: RFLU_ModCellMapping.F90:648
rflu_moddimensions::rflu_setmaxdimension
INTEGER function, public rflu_setmaxdimension(global, nXyzTot)
Definition: RFLU_ModDimensions.F90:559
j
j indices j
Definition: Indexing.h:6
rflu_moddimensions
Definition: RFLU_ModDimensions.F90:39
modglobal::t_global
Definition: ModGlobal.F90:50
moderror::errorstop
subroutine errorstop(global, errorCode, errorLine, addMessage)
Definition: ModError.F90:483
rflu_modgeometry::rflu_computeapproxcentroids
subroutine, public rflu_computeapproxcentroids(pRegion)
Definition: RFLU_ModGeometry.F90:1208
rflu_modcellmapping::rflu_createcellmapping
subroutine, public rflu_createcellmapping(pRegion)
Definition: RFLU_ModCellMapping.F90:181
facevectorquad
subroutine facevectorquad(xyzNodes, fVecX, fVecY, fVecZ)
Definition: FaceVector.F90:80
key
static T_Key key
Definition: vinci_lass.c:76
moderror::deregisterfunction
subroutine deregisterfunction(global)
Definition: ModError.F90:469
rflu_modgeometry::rflu_createapproxcentroids
subroutine, public rflu_createapproxcentroids(pRegion)
Definition: RFLU_ModGeometry.F90:1558
facecentroidtria
subroutine facecentroidtria(xyz, fCenX, fCenY, fCenZ)
Definition: FaceCentroid.F90:47
facecentroidquad
subroutine facecentroidquad(xyz, fCenX, fCenY, fCenZ)
Definition: FaceCentroid.F90:73
rflu_modgeometry
Definition: RFLU_ModGeometry.F90:46
modbndpatch::t_patch
Definition: ModBndPatch.F90:91
rflu_modcellmapping::rflu_buildloc2globcellmapping
subroutine, public rflu_buildloc2globcellmapping(pRegion)
Definition: RFLU_ModCellMapping.F90:463