Rocstar  1.0
Rocstar multiphysics simulation application
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
RFLO_ModGridRegionShape.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: Suite for region shape related routines.
26 !
27 ! Description: None.
28 !
29 ! Notes: None.
30 !
31 ! ******************************************************************************
32 !
33 ! $Id: RFLO_ModGridRegionShape.F90,v 1.3 2008/12/06 08:44:16 mtcampbe Exp $
34 !
35 ! Copyright: (c) 2004 by the University of Illinois
36 !
37 ! ******************************************************************************
38 
39 MODULE rflo_modgridregionshape
40 
41  USE modglobal, ONLY : t_global
42  USE moddatastruct, ONLY: t_region
43  USE modgrid, ONLY : t_grid
44  USE modbndpatch, ONLY : t_patch
45  USE modparameters
46  USE moddatatypes
47  USE moderror
48  USE modmpi
49 
50  IMPLICIT NONE
51 
52  PRIVATE
53  PUBLIC :: rflo_gridflatpatch, &
54  rflo_findfunkyblocks
55 
56 ! private : RFLO_BlockSkewedCell
57 
58 ! ******************************************************************************
59 ! Declarations and definitions
60 ! ******************************************************************************
61 
62  CHARACTER(CHRLEN) :: RCSIdentString = &
63  '$RCSfile: RFLO_ModGridRegionShape.F90,v $ $Revision: 1.3 $'
64 
65 ! ******************************************************************************
66 ! Routines
67 ! ******************************************************************************
68 
69  CONTAINS
70 
71 !******************************************************************************
72 !
73 ! Purpose: identify flat patches
74 !
75 ! Description: none.
76 !
77 ! Input: region = data of current region, grid movements
78 ! patch = current patch.
79 !
80 ! Output: patch%bndFlat get value true or false
81 !
82 ! Notes: none.
83 !
84 !******************************************************************************
85 
86 SUBROUTINE rflo_gridflatpatch( region,patch )
87 
88  USE modinterfaces, ONLY : rflo_getpatchindicesnodes, rflo_getnodeoffset
89  USE rflo_modvectortensor, ONLY : rflo_normcrossprod
90 
91  IMPLICIT NONE
92 
93 #include "Indexing.h"
94 
95 ! ... parameters
96  TYPE(t_region) :: region
97  TYPE(t_patch) :: patch
98 
99 ! ... loop variables
100  INTEGER :: i, j, k
101 
102 ! ... local variables
103  INTEGER :: ilev, lbound, ihlf, jhlf, khlf
104  INTEGER :: ibeg, iend, jbeg, jend, kbeg, kend, inoff, ijnoff
105 
106  REAL(RFREAL) :: vmag, eps
107  REAL(RFREAL) :: v11(xcoord:zcoord), v12(xcoord:zcoord)
108  REAL(RFREAL) :: v21(xcoord:zcoord), v22(xcoord:zcoord)
109  REAL(RFREAL) :: v31(xcoord:zcoord), v32(xcoord:zcoord)
110  REAL(RFREAL) :: v41(xcoord:zcoord), v42(xcoord:zcoord)
111  REAL(RFREAL) :: v1(xcoord:zcoord), v2(xcoord:zcoord)
112  REAL(RFREAL) :: v3(xcoord:zcoord), v4(xcoord:zcoord), vec(xcoord:zcoord)
113  REAL(RFREAL), POINTER :: xyz(:,:)
114 
115 !******************************************************************************
116 
117  CALL registerfunction( region%global,'RFLO_GridFlatPatch',&
118  'RFLO_ModGridRegionShape.F90' )
119 
120 ! get dimensions and pointers -------------------------------------------------
121 
122  ilev = 1
123  lbound = patch%lbound
124 
125  CALL rflo_getpatchindicesnodes( region,patch,ilev,ibeg,iend, &
126  jbeg,jend,kbeg,kend )
127  CALL rflo_getnodeoffset( region,ilev,inoff,ijnoff )
128 
129  xyz => region%levels(ilev)%grid%xyz
130 
131  ihlf = ibeg+(iend-ibeg)/2
132  jhlf = jbeg+(jend-jbeg)/2
133  khlf = kbeg+(kend-kbeg)/2
134 
135 ! default values --------------------------------------------------------------
136 
137  patch%bndFlat = .false.
138  patch%dirFlat = -1
139 
140 ! four pair vectors with consistent orientation -------------------------------
141 
142  IF (lbound==1 .OR. lbound==2) THEN
143  v11(:) = xyz(:,indijk(ibeg,jhlf,kbeg,inoff,ijnoff)) - &
144  xyz(:,indijk(ibeg,jbeg,kbeg,inoff,ijnoff))
145  v12(:) = xyz(:,indijk(ibeg,jbeg,khlf,inoff,ijnoff)) - &
146  xyz(:,indijk(ibeg,jbeg,kbeg,inoff,ijnoff))
147  v21(:) = xyz(:,indijk(ibeg,jbeg,khlf,inoff,ijnoff)) - &
148  xyz(:,indijk(ibeg,jbeg,kend,inoff,ijnoff))
149  v22(:) = xyz(:,indijk(ibeg,jhlf,kend,inoff,ijnoff)) - &
150  xyz(:,indijk(ibeg,jbeg,kend,inoff,ijnoff))
151  v31(:) = xyz(:,indijk(ibeg,jhlf,kend,inoff,ijnoff)) - &
152  xyz(:,indijk(ibeg,jend,kend,inoff,ijnoff))
153  v32(:) = xyz(:,indijk(ibeg,jend,khlf,inoff,ijnoff)) - &
154  xyz(:,indijk(ibeg,jend,kend,inoff,ijnoff))
155  v41(:) = xyz(:,indijk(ibeg,jend,khlf,inoff,ijnoff)) - &
156  xyz(:,indijk(ibeg,jend,kbeg,inoff,ijnoff))
157  v42(:) = xyz(:,indijk(ibeg,jhlf,kbeg,inoff,ijnoff)) - &
158  xyz(:,indijk(ibeg,jend,kbeg,inoff,ijnoff))
159  ELSEIF (lbound==3 .OR. lbound==4) THEN
160  v11(:) = xyz(:,indijk(ihlf,jbeg,kbeg,inoff,ijnoff)) - &
161  xyz(:,indijk(ibeg,jbeg,kbeg,inoff,ijnoff))
162  v12(:) = xyz(:,indijk(ibeg,jbeg,khlf,inoff,ijnoff)) - &
163  xyz(:,indijk(ibeg,jbeg,kbeg,inoff,ijnoff))
164  v21(:) = xyz(:,indijk(ibeg,jbeg,khlf,inoff,ijnoff)) - &
165  xyz(:,indijk(ibeg,jbeg,kend,inoff,ijnoff))
166  v22(:) = xyz(:,indijk(ihlf,jbeg,kend,inoff,ijnoff)) - &
167  xyz(:,indijk(ibeg,jbeg,kend,inoff,ijnoff))
168  v31(:) = xyz(:,indijk(ihlf,jbeg,kend,inoff,ijnoff)) - &
169  xyz(:,indijk(iend,jbeg,kend,inoff,ijnoff))
170  v32(:) = xyz(:,indijk(iend,jbeg,khlf,inoff,ijnoff)) - &
171  xyz(:,indijk(iend,jbeg,kend,inoff,ijnoff))
172  v41(:) = xyz(:,indijk(iend,jbeg,khlf,inoff,ijnoff)) - &
173  xyz(:,indijk(iend,jbeg,kbeg,inoff,ijnoff))
174  v42(:) = xyz(:,indijk(ihlf,jbeg,kbeg,inoff,ijnoff)) - &
175  xyz(:,indijk(iend,jbeg,kbeg,inoff,ijnoff))
176  ELSEIF (lbound==5 .OR. lbound==6) THEN
177  v11(:) = xyz(:,indijk(ihlf,jbeg,kbeg,inoff,ijnoff)) - &
178  xyz(:,indijk(ibeg,jbeg,kbeg,inoff,ijnoff))
179  v12(:) = xyz(:,indijk(ibeg,jhlf,kbeg,inoff,ijnoff)) - &
180  xyz(:,indijk(ibeg,jbeg,kbeg,inoff,ijnoff))
181  v21(:) = xyz(:,indijk(ibeg,jhlf,kbeg,inoff,ijnoff)) - &
182  xyz(:,indijk(ibeg,jend,kbeg,inoff,ijnoff))
183  v22(:) = xyz(:,indijk(ihlf,jend,kbeg,inoff,ijnoff)) - &
184  xyz(:,indijk(ibeg,jend,kbeg,inoff,ijnoff))
185  v31(:) = xyz(:,indijk(ihlf,jend,kbeg,inoff,ijnoff)) - &
186  xyz(:,indijk(iend,jend,kbeg,inoff,ijnoff))
187  v32(:) = xyz(:,indijk(iend,jhlf,kbeg,inoff,ijnoff)) - &
188  xyz(:,indijk(iend,jend,kbeg,inoff,ijnoff))
189  v41(:) = xyz(:,indijk(iend,jhlf,kbeg,inoff,ijnoff)) - &
190  xyz(:,indijk(iend,jbeg,kbeg,inoff,ijnoff))
191  v42(:) = xyz(:,indijk(ihlf,jbeg,kbeg,inoff,ijnoff)) - &
192  xyz(:,indijk(iend,jbeg,kbeg,inoff,ijnoff))
193  ENDIF
194 
195 ! unit vectors orthogonal to the patch at four patch corners ------------------
196 
197  CALL rflo_normcrossprod( v11,v12,v1 )
198  CALL rflo_normcrossprod( v21,v22,v2 )
199  CALL rflo_normcrossprod( v31,v32,v3 )
200  CALL rflo_normcrossprod( v41,v42,v4 )
201 
202 ! their total length should equals 4. for flat patch --------------------------
203 
204  eps = 1.e-3_rfreal
205  vec = v1+v2+v3+v4
206  vmag = sqrt( vec(xcoord)**2 + vec(ycoord)**2 + vec(zcoord)**2 )
207  IF (abs( vmag-4._rfreal) < eps) THEN
208  patch%bndFlat = .true.
209  patch%dirFlat = off
210  ELSE
211  region%levels(ilev)%grid%boundFlat(lbound) = .false.
212  ENDIF
213 
214 ! search for bend edges -------------------------------------------------------
215 
216  IF (lbound==1) THEN
217  v1 = v12/sqrt( v12(xcoord)**2 + v12(ycoord)**2 + v12(zcoord)**2 )
218  v2 = v21/sqrt( v21(xcoord)**2 + v21(ycoord)**2 + v21(zcoord)**2 )
219  vec = v1 - v2
220  vmag = sqrt( vec(xcoord)**2 + vec(ycoord)**2 + vec(zcoord)**2 )
221  IF (abs( vmag-2._rfreal) > eps) THEN
222  region%levels(ilev)%grid%edgeStraight(1) = .false.
223  ENDIF
224 
225  v1 = v22/sqrt( v22(xcoord)**2 + v22(ycoord)**2 + v22(zcoord)**2 )
226  v2 = v31/sqrt( v31(xcoord)**2 + v31(ycoord)**2 + v31(zcoord)**2 )
227  vec = v1 - v2
228  vmag = sqrt( vec(xcoord)**2 + vec(ycoord)**2 + vec(zcoord)**2 )
229  IF (abs( vmag-2._rfreal) > eps) THEN
230  region%levels(ilev)%grid%edgeStraight(2) = .false.
231  ENDIF
232 
233  v1 = v32/sqrt( v32(xcoord)**2 + v32(ycoord)**2 + v32(zcoord)**2 )
234  v2 = v41/sqrt( v41(xcoord)**2 + v41(ycoord)**2 + v41(zcoord)**2 )
235  vec = v1 - v2
236  vmag = sqrt( vec(xcoord)**2 + vec(ycoord)**2 + vec(zcoord)**2 )
237  IF (abs( vmag-2._rfreal) > eps) THEN
238  region%levels(ilev)%grid%edgeStraight(3) = .false.
239  ENDIF
240 
241  v1 = v42/sqrt( v42(xcoord)**2 + v42(ycoord)**2 + v42(zcoord)**2 )
242  v2 = v11/sqrt( v11(xcoord)**2 + v11(ycoord)**2 + v11(zcoord)**2 )
243  vec = v1 - v2
244  vmag = sqrt( vec(xcoord)**2 + vec(ycoord)**2 + vec(zcoord)**2 )
245  IF (abs( vmag-2._rfreal) > eps) THEN
246  region%levels(ilev)%grid%edgeStraight(4) = .false.
247  ENDIF
248 
249  ELSEIF (lbound==2) THEN
250  v1 = v12/sqrt( v12(xcoord)**2 + v12(ycoord)**2 + v12(zcoord)**2 )
251  v2 = v21/sqrt( v21(xcoord)**2 + v21(ycoord)**2 + v21(zcoord)**2 )
252  vec = v1 - v2
253  vmag = sqrt( vec(xcoord)**2 + vec(ycoord)**2 + vec(zcoord)**2 )
254  IF (abs( vmag-2._rfreal) > eps) THEN
255  region%levels(ilev)%grid%edgeStraight(5) = .false.
256  ENDIF
257 
258  v1 = v22/sqrt( v22(xcoord)**2 + v22(ycoord)**2 + v22(zcoord)**2 )
259  v2 = v31/sqrt( v31(xcoord)**2 + v31(ycoord)**2 + v31(zcoord)**2 )
260  vec = v1 - v2
261  vmag = sqrt( vec(xcoord)**2 + vec(ycoord)**2 + vec(zcoord)**2 )
262  IF (abs( vmag-2._rfreal) > eps) THEN
263  region%levels(ilev)%grid%edgeStraight(6) = .false.
264  ENDIF
265 
266  v1 = v32/sqrt( v32(xcoord)**2 + v32(ycoord)**2 + v32(zcoord)**2 )
267  v2 = v41/sqrt( v41(xcoord)**2 + v41(ycoord)**2 + v41(zcoord)**2 )
268  vec = v1 - v2
269  vmag = sqrt( vec(xcoord)**2 + vec(ycoord)**2 + vec(zcoord)**2 )
270  IF (abs( vmag-2._rfreal) > eps) THEN
271  region%levels(ilev)%grid%edgeStraight(7) = .false.
272  ENDIF
273 
274  v1 = v42/sqrt( v42(xcoord)**2 + v42(ycoord)**2 + v42(zcoord)**2 )
275  v2 = v11/sqrt( v11(xcoord)**2 + v11(ycoord)**2 + v11(zcoord)**2 )
276  vec = v1 - v2
277  vmag = sqrt( vec(xcoord)**2 + vec(ycoord)**2 + vec(zcoord)**2 )
278  IF (abs( vmag-2._rfreal) > eps) THEN
279  region%levels(ilev)%grid%edgeStraight(8) = .false.
280  ENDIF
281 
282  ELSEIF (lbound==3) THEN
283  v1 = v11/sqrt( v11(xcoord)**2 + v11(ycoord)**2 + v11(zcoord)**2 )
284  v2 = v42/sqrt( v42(xcoord)**2 + v42(ycoord)**2 + v42(zcoord)**2 )
285  vec = v1 - v2
286  vmag = sqrt( vec(xcoord)**2 + vec(ycoord)**2 + vec(zcoord)**2 )
287  IF (abs( vmag-2._rfreal) > eps) THEN
288  region%levels(ilev)%grid%edgeStraight(9) = .false.
289  ENDIF
290 
291  v1 = v22/sqrt( v22(xcoord)**2 + v22(ycoord)**2 + v22(zcoord)**2 )
292  v2 = v31/sqrt( v31(xcoord)**2 + v31(ycoord)**2 + v31(zcoord)**2 )
293  vec = v1 - v2
294  vmag = sqrt( vec(xcoord)**2 + vec(ycoord)**2 + vec(zcoord)**2 )
295  IF (abs( vmag-2._rfreal) > eps) THEN
296  region%levels(ilev)%grid%edgeStraight(10) = .false.
297  ENDIF
298 
299  ELSEIF (lbound==4) THEN
300  v1 = v11/sqrt( v11(xcoord)**2 + v11(ycoord)**2 + v11(zcoord)**2 )
301  v2 = v42/sqrt( v42(xcoord)**2 + v42(ycoord)**2 + v42(zcoord)**2 )
302  vec = v1 - v2
303  vmag = sqrt( vec(xcoord)**2 + vec(ycoord)**2 + vec(zcoord)**2 )
304  IF (abs( vmag-2._rfreal) > eps) THEN
305  region%levels(ilev)%grid%edgeStraight(11) = .false. ! officially this
306  ENDIF ! should be edge 12
307 
308  v1 = v22/sqrt( v22(xcoord)**2 + v22(ycoord)**2 + v22(zcoord)**2 )
309  v2 = v31/sqrt( v31(xcoord)**2 + v31(ycoord)**2 + v31(zcoord)**2 )
310  vec = v1 - v2
311  vmag = sqrt( vec(xcoord)**2 + vec(ycoord)**2 + vec(zcoord)**2 )
312  IF (abs( vmag-2._rfreal) > eps) THEN
313  region%levels(ilev)%grid%edgeStraight(12) = .false. ! officially this
314  ENDIF ! should be edge 11
315 
316  ENDIF
317 
318 ! search for flat direction on curved patches ---------------------------------
319 
320  IF (.NOT. patch%bndFlat) THEN
321  IF (lbound==1) THEN
322  IF ((region%levels(ilev)%grid%edgeStraight(1) .EQV. .true.) .AND. &
323  (region%levels(ilev)%grid%edgeStraight(3) .EQV. .true.)) THEN
324  patch%dirFlat = kcoord
325  ELSEIF ((region%levels(ilev)%grid%edgeStraight(2) .EQV. .true.) .AND. &
326  (region%levels(ilev)%grid%edgeStraight(4) .EQV. .true.)) THEN
327  patch%dirFlat = jcoord
328  ENDIF
329  ELSEIF (lbound==2) THEN
330  IF ((region%levels(ilev)%grid%edgeStraight(5) .EQV. .true.) .AND. &
331  (region%levels(ilev)%grid%edgeStraight(7) .EQV. .true.)) THEN
332  patch%dirFlat = kcoord
333  ELSEIF ((region%levels(ilev)%grid%edgeStraight(6) .EQV. .true.) .AND. &
334  (region%levels(ilev)%grid%edgeStraight(8) .EQV. .true.)) THEN
335  patch%dirFlat = jcoord
336  ENDIF
337  ELSEIF (lbound==3) THEN
338  IF ((region%levels(ilev)%grid%edgeStraight(1) .EQV. .true.) .AND. &
339  (region%levels(ilev)%grid%edgeStraight(5) .EQV. .true.)) THEN
340  patch%dirFlat = kcoord
341  ELSEIF ((region%levels(ilev)%grid%edgeStraight(9) .EQV. .true.) .AND. &
342  (region%levels(ilev)%grid%edgeStraight(10) .EQV. .true.)) THEN
343  patch%dirFlat = icoord
344  ENDIF
345  ELSEIF (lbound==4) THEN
346  IF ((region%levels(ilev)%grid%edgeStraight(3) .EQV. .true.) .AND. &
347  (region%levels(ilev)%grid%edgeStraight(7) .EQV. .true.)) THEN
348  patch%dirFlat = kcoord
349  ELSEIF ((region%levels(ilev)%grid%edgeStraight(11) .EQV. .true.) .AND. &
350  (region%levels(ilev)%grid%edgeStraight(12) .EQV. .true.)) THEN
351  patch%dirFlat = icoord
352  ENDIF
353  ELSEIF (lbound==5) THEN
354  IF ((region%levels(ilev)%grid%edgeStraight(4) .EQV. .true.) .AND. &
355  (region%levels(ilev)%grid%edgeStraight(8) .EQV. .true.)) THEN
356  patch%dirFlat = jcoord
357  ELSEIF ((region%levels(ilev)%grid%edgeStraight(9) .EQV. .true.) .AND. &
358  (region%levels(ilev)%grid%edgeStraight(11) .EQV. .true.)) THEN
359  patch%dirFlat = icoord
360  ENDIF
361  ELSEIF (lbound==6) THEN
362  IF ((region%levels(ilev)%grid%edgeStraight(2) .EQV. .true.) .AND. &
363  (region%levels(ilev)%grid%edgeStraight(6) .EQV. .true.)) THEN
364  patch%dirFlat = jcoord
365  ELSEIF ((region%levels(ilev)%grid%edgeStraight(10) .EQV. .true.) .AND. &
366  (region%levels(ilev)%grid%edgeStraight(12) .EQV. .true.)) THEN
367  patch%dirFlat = icoord
368  ENDIF ! edgeStraight
369  ENDIF ! lbound
370  ENDIF ! patch not flat
371 
372 ! finalize --------------------------------------------------------------------
373 
374  CALL deregisterfunction( region%global )
375 
376 END SUBROUTINE rflo_gridflatpatch
377 
378 
379 !******************************************************************************
380 !
381 ! Purpose: identify blocks with unconventional shape, if any.
382 !
383 ! Description: none.
384 !
385 ! Input: regions%grid = dimensions, grid coordinates.
386 !
387 ! Output: regions%blockShape = normal or funky
388 !
389 ! Notes: none.
390 !
391 !******************************************************************************
392 
393 SUBROUTINE rflo_findfunkyblocks( regions )
394 
395  USE modinterfaces, ONLY : rflo_calcfacevectors
396 
397  IMPLICIT NONE
398 
399 ! ... parameters
400  TYPE (t_region), POINTER :: regions(:)
401 
402 ! ... loop variables
403  INTEGER :: ireg
404 
405 !******************************************************************************
406 
407  CALL registerfunction( regions(1)%global,'RFLO_FindFunkyBlocks',&
408  'RFLO_ModGridRegionShape.F90' )
409 
410 ! IF ( regions(1)%global%myProcid == MASTERPROC .AND. &
411 ! regions(1)%global%verbLevel > VERBOSE_NONE ) THEN
412 ! WRITE(STDOUT,'(A,1X,A)') SOLVER_NAME,'Entering RFLO_FindFunkyBlocks...'
413 ! END IF ! global%verbLevel
414 
415 ! loop over all regions
416 
417  DO ireg=1,regions(1)%global%nRegions
418  IF (regions(ireg)%procid==regions(ireg)%global%myProcid & ! region active and
419  .AND. regions(ireg)%active==active) THEN ! on my processor
420 
421  CALL rflo_calcfacevectors( regions(ireg) )
422  CALL rflo_blockskewedcell( regions(ireg) )
423 
424  ENDIF ! region on this processor and active
425  ENDDO ! iReg
426 
427 ! finalize --------------------------------------------------------------------
428 
429 ! IF ( regions(1)%global%myProcid == MASTERPROC .AND. &
430 ! regions(1)%global%verbLevel > VERBOSE_NONE ) THEN
431 ! WRITE(STDOUT,'(A,1X,A)') SOLVER_NAME,'Leaving RFLO_FindFunkyBlocks.'
432 ! END IF ! global%verbLevel
433 
434  CALL deregisterfunction( regions(1)%global )
435 
436 END SUBROUTINE rflo_findfunkyblocks
437 
438 
439 !******************************************************************************
440 !
441 ! Purpose: identify if this region contains skewed cells.
442 !
443 ! Description: if yes, mark it as funky block, otherwise normal block.
444 !
445 ! Input: region%levels%grid = dimensions, coordinates, face vectors
446 ! (current region)
447 !
448 ! Output: for finest grid only.
449 !
450 !******************************************************************************
451 
452 SUBROUTINE rflo_blockskewedcell( region )
453 
454  USE modinterfaces, ONLY : rflo_getdimensphys, rflo_getnodeoffset
455 
456  IMPLICIT NONE
457 
458 #include "Indexing.h"
459 
460 ! ... parameters
461  TYPE(t_region) :: region
462 
463 ! ... loop variables
464  INTEGER :: i, j, k, l
465 
466 ! ... local variables
467 
468  INTEGER :: ipcbeg, ipcend, jpcbeg, jpcend, kpcbeg, kpcend
469  INTEGER :: ilev, inoff, ijnoff, ijknode(8), face(6)
470  INTEGER :: im, jm, km
471 
472  REAL(RFREAL) :: edge(3,4), dprod(8)
473  REAL(RFREAL) :: dprodm, dpmax, emag, smag1, smag2
474  REAL(RFREAL), POINTER :: xyz(:,:), si(:,:), sj(:,:), sk(:,:)
475 
476 !******************************************************************************
477 
478  CALL registerfunction( region%global,'RFLO_BlockSkewedCell',&
479  'RFLO_ModGridRegionShape.F90' )
480 
481 ! compute cell skewness -------------------------------------------------------
482 
483  ilev=1
484 
485  CALL rflo_getdimensphys( region,ilev,ipcbeg,ipcend, &
486  jpcbeg,jpcend,kpcbeg,kpcend )
487  CALL rflo_getnodeoffset( region,ilev,inoff,ijnoff )
488 
489  xyz => region%levels(ilev)%grid%xyz
490  si => region%levels(ilev)%grid%si
491  sj => region%levels(ilev)%grid%sj
492  sk => region%levels(ilev)%grid%sk
493 
494 ! check face vectors: cell skewness
495 
496  dprodm = 1._rfreal
497  dpmax = -1._rfreal
498 
499  DO k=kpcbeg,kpcend
500  DO j=jpcbeg,jpcend
501  DO i=ipcbeg,ipcend
502  ijknode(1) = indijk(i ,j ,k ,inoff,ijnoff)
503  ijknode(2) = indijk(i+1,j ,k ,inoff,ijnoff)
504  ijknode(3) = indijk(i ,j+1,k ,inoff,ijnoff)
505  ijknode(4) = indijk(i ,j ,k+1,inoff,ijnoff)
506  ijknode(5) = indijk(i+1,j+1,k+1,inoff,ijnoff)
507  ijknode(6) = indijk(i+1,j+1,k ,inoff,ijnoff)
508  ijknode(7) = indijk(i+1,j ,k+1,inoff,ijnoff)
509  ijknode(8) = indijk(i ,j+1,k+1,inoff,ijnoff)
510 
511  face(1) = ijknode(1)
512  face(2) = ijknode(2)
513  face(4) = ijknode(3)
514  face(6) = ijknode(4)
515 
516 ! ----- check for skewed i-faces
517 
518  edge(1,1) = xyz(xcoord,ijknode(1))-xyz(xcoord,ijknode(2))
519  edge(2,1) = xyz(ycoord,ijknode(1))-xyz(ycoord,ijknode(2))
520  edge(3,1) = xyz(zcoord,ijknode(1))-xyz(zcoord,ijknode(2))
521 
522  edge(1,2) = xyz(xcoord,ijknode(3))-xyz(xcoord,ijknode(6))
523  edge(2,2) = xyz(ycoord,ijknode(3))-xyz(ycoord,ijknode(6))
524  edge(3,2) = xyz(zcoord,ijknode(3))-xyz(zcoord,ijknode(6))
525 
526  edge(1,3) = xyz(xcoord,ijknode(8))-xyz(xcoord,ijknode(5))
527  edge(2,3) = xyz(ycoord,ijknode(8))-xyz(ycoord,ijknode(5))
528  edge(3,3) = xyz(zcoord,ijknode(8))-xyz(zcoord,ijknode(5))
529 
530  edge(1,4) = xyz(xcoord,ijknode(4))-xyz(xcoord,ijknode(7))
531  edge(2,4) = xyz(ycoord,ijknode(4))-xyz(ycoord,ijknode(7))
532  edge(3,4) = xyz(zcoord,ijknode(4))-xyz(zcoord,ijknode(7))
533 
534  DO l=1,4
535  emag = sqrt( edge(1,l)*edge(1,l) + &
536  edge(2,l)*edge(2,l) + &
537  edge(3,l)*edge(3,l) )
538  smag1 = sqrt( si(1,face(1))*si(1,face(1)) + &
539  si(2,face(1))*si(2,face(1)) + &
540  si(3,face(1))*si(3,face(1)) )
541  smag2 = sqrt( si(1,face(2))*si(1,face(2)) + &
542  si(2,face(2))*si(2,face(2)) + &
543  si(3,face(2))*si(3,face(2)) )
544 
545  dprod(l) = (si(1,face(1))*edge(1,l) + &
546  si(2,face(1))*edge(2,l) + &
547  si(3,face(1))*edge(3,l))/(smag1*emag)
548 
549  dprod(l+4) = (si(1,face(2))*edge(1,l) + &
550  si(2,face(2))*edge(2,l) + &
551  si(3,face(2))*edge(3,l))/(smag2*emag)
552  ENDDO
553 
554  IF (minval(dprod) < dprodm) THEN
555  dprodm = minval( dprod )
556  im = i
557  jm = j
558  km = k
559  ENDIF
560 
561  IF (maxval(dprod) > dpmax) THEN
562  dpmax = maxval( dprod )
563  ENDIF
564 
565 ! ----- check for skewed j-faces
566 
567  edge(1,1) = xyz(xcoord,ijknode(1))-xyz(xcoord,ijknode(3))
568  edge(2,1) = xyz(ycoord,ijknode(1))-xyz(ycoord,ijknode(3))
569  edge(3,1) = xyz(zcoord,ijknode(1))-xyz(zcoord,ijknode(3))
570 
571  edge(1,2) = xyz(xcoord,ijknode(4))-xyz(xcoord,ijknode(8))
572  edge(2,2) = xyz(ycoord,ijknode(4))-xyz(ycoord,ijknode(8))
573  edge(3,2) = xyz(zcoord,ijknode(4))-xyz(zcoord,ijknode(8))
574 
575  edge(1,3) = xyz(xcoord,ijknode(7))-xyz(xcoord,ijknode(5))
576  edge(2,3) = xyz(ycoord,ijknode(7))-xyz(ycoord,ijknode(5))
577  edge(3,3) = xyz(zcoord,ijknode(7))-xyz(zcoord,ijknode(5))
578 
579  edge(1,4) = xyz(xcoord,ijknode(2))-xyz(xcoord,ijknode(6))
580  edge(2,4) = xyz(ycoord,ijknode(2))-xyz(ycoord,ijknode(6))
581  edge(3,4) = xyz(zcoord,ijknode(2))-xyz(zcoord,ijknode(6))
582 
583  DO l=1,4
584  emag = sqrt( edge(1,l)*edge(1,l) + &
585  edge(2,l)*edge(2,l) + &
586  edge(3,l)*edge(3,l) )
587  smag1 = sqrt( sj(1,face(1))*sj(1,face(1)) + &
588  sj(2,face(1))*sj(2,face(1)) + &
589  sj(3,face(1))*sj(3,face(1)) )
590  smag2 = sqrt( sj(1,face(4))*sj(1,face(4)) + &
591  sj(2,face(4))*sj(2,face(4)) + &
592  sj(3,face(4))*sj(3,face(4)) )
593 
594  dprod(l) = (sj(1,face(1))*edge(1,l) + &
595  sj(2,face(1))*edge(2,l) + &
596  sj(3,face(1))*edge(3,l))/(smag1*emag)
597 
598  dprod(l+4) = (sj(1,face(4))*edge(1,l) + &
599  sj(2,face(4))*edge(2,l) + &
600  sj(3,face(4))*edge(3,l))/(smag2*emag)
601  ENDDO
602 
603  IF (minval(dprod) < dprodm) THEN
604  dprodm = minval( dprod )
605  im = i
606  jm = j
607  km = k
608  ENDIF
609 
610  IF (maxval(dprod) > dpmax) THEN
611  dpmax = maxval( dprod )
612  ENDIF
613 
614 ! ----- check for skewed k-faces
615 
616  edge(1,1) = xyz(xcoord,ijknode(1))-xyz(xcoord,ijknode(4))
617  edge(2,1) = xyz(ycoord,ijknode(1))-xyz(ycoord,ijknode(4))
618  edge(3,1) = xyz(zcoord,ijknode(1))-xyz(zcoord,ijknode(4))
619 
620  edge(1,2) = xyz(xcoord,ijknode(2))-xyz(xcoord,ijknode(7))
621  edge(2,2) = xyz(ycoord,ijknode(2))-xyz(ycoord,ijknode(7))
622  edge(3,2) = xyz(zcoord,ijknode(2))-xyz(zcoord,ijknode(7))
623 
624  edge(1,3) = xyz(xcoord,ijknode(6))-xyz(xcoord,ijknode(5))
625  edge(2,3) = xyz(ycoord,ijknode(6))-xyz(ycoord,ijknode(5))
626  edge(3,3) = xyz(zcoord,ijknode(6))-xyz(zcoord,ijknode(5))
627 
628  edge(1,4) = xyz(xcoord,ijknode(3))-xyz(xcoord,ijknode(8))
629  edge(2,4) = xyz(ycoord,ijknode(3))-xyz(ycoord,ijknode(8))
630  edge(3,4) = xyz(zcoord,ijknode(3))-xyz(zcoord,ijknode(8))
631 
632  DO l=1,4
633  emag = sqrt( edge(1,l)*edge(1,l) + &
634  edge(2,l)*edge(2,l) + &
635  edge(3,l)*edge(3,l) )
636  smag1 = sqrt( sk(1,face(1))*sk(1,face(1)) + &
637  sk(2,face(1))*sk(2,face(1)) + &
638  sk(3,face(1))*sk(3,face(1)) )
639  smag2 = sqrt( sk(1,face(6))*sk(1,face(6)) + &
640  sk(2,face(6))*sk(2,face(6)) + &
641  sk(3,face(6))*sk(3,face(6)) )
642 
643  dprod(l) = (sk(1,face(1))*edge(1,l) + &
644  sk(2,face(1))*edge(2,l) + &
645  sk(3,face(1))*edge(3,l))/(smag1*emag)
646 
647  dprod(l+4) = (sk(1,face(6))*edge(1,l) + &
648  sk(2,face(6))*edge(2,l) + &
649  sk(3,face(6))*edge(3,l))/(smag2*emag)
650  ENDDO
651 
652  IF (minval(dprod) < dprodm) THEN
653  dprodm = minval( dprod )
654  im = i
655  jm = j
656  km = k
657  ENDIF
658 
659  IF (maxval(dprod) > dpmax) THEN
660  dpmax = maxval( dprod )
661  ENDIF
662 
663  ENDDO
664  ENDDO
665  ENDDO
666 
667  IF (dprodm < 1.e-1_rfreal) THEN
668  region%blockShape = region_shape_funky
669 
670  IF ( region%global%verbLevel >= verbose_none ) THEN
671  WRITE(stdout,'(A,1X,A,I6)') solver_name, &
672  'Warning: found funky block, region',region%iRegionGlobal
673  END IF ! global%verbLevel
674  ENDIF
675 
676 ! finalize --------------------------------------------------------------------
677 
678  CALL deregisterfunction( region%global )
679 
680 END SUBROUTINE rflo_blockskewedcell
681 
682 
683 ! ******************************************************************************
684 ! End
685 ! ******************************************************************************
686 
687 END MODULE rflo_modgridregionshape
688 
689 
690 ! ******************************************************************************
691 !
692 ! RCS Revision history:
693 !
694 ! $Log: RFLO_ModGridRegionShape.F90,v $
695 ! Revision 1.3 2008/12/06 08:44:16 mtcampbe
696 ! Updated license.
697 !
698 ! Revision 1.2 2008/11/19 22:17:27 mtcampbe
699 ! Added Illinois Open Source License/Copyright
700 !
701 ! Revision 1.1 2006/03/18 11:09:30 wasistho
702 ! initial import
703 !
704 !
705 !
706 ! ******************************************************************************
707 
708 
709 
710 
711 
712 
713 
714 
715 
ibeg
**********************************************************************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 ibeg
Definition: RFLO_GetPatchIndices.h:27
rflo_modgridregionshape
Definition: RFLO_ModGridRegionShape.F90:39
rflo_modvectortensor
Definition: RFLO_ModVectorTensor.F90:39
face
Definition: face.h:90
k
j indices k indices k
Definition: Indexing.h:6
kpcbeg
**********************************************************************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 kpcbeg
Definition: RFLO_GetDimensPhys.h:27
moderror::registerfunction
subroutine registerfunction(global, funName, fileName)
Definition: ModError.F90:449
rflo_modvectortensor::rflo_normcrossprod
subroutine, public rflo_normcrossprod(s1, s2, s3)
Definition: RFLO_ModVectorTensor.F90:84
rflo_modgridregionshape::rflo_findfunkyblocks
subroutine, public rflo_findfunkyblocks(regions)
Definition: RFLO_ModGridRegionShape.F90:393
modgrid::t_grid
Definition: ModGrid.F90:56
sqrt
double sqrt(double d)
Definition: double.h:73
jpcbeg
**********************************************************************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 jpcbeg
Definition: RFLO_GetDimensPhys.h:27
ipcend
**********************************************************************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 ipcend
Definition: RFLO_GetDimensPhys.h:27
rflo_getnodeoffset
subroutine rflo_getnodeoffset(region, iLev, iNodeOffset, ijNodeOffset)
Definition: RFLO_GetDimensions.F90:267
patch
Definition: patch.h:74
iend
**********************************************************************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 knode iend
Definition: RFLO_CalcGradDummySymm.h:32
rflo_calcfacevectors
subroutine rflo_calcfacevectors(region)
Definition: RFLO_CalcFaceVectors.F90:45
ipcbeg
**********************************************************************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 ipcbeg
Definition: RFLO_GetDimensPhys.h:27
rflo_getpatchindicesnodes
subroutine rflo_getpatchindicesnodes(region, patch, iLev, ibeg, iend, jbeg, jend, kbeg, kend)
Definition: RFLO_GetPatchDimens.F90:138
i
blockLoc i
Definition: read.cpp:79
j
j indices j
Definition: Indexing.h:6
jpcend
**********************************************************************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 jpcend
Definition: RFLO_GetDimensPhys.h:27
jend
**********************************************************************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 knode jend
Definition: RFLO_CalcGradDummySymm.h:32
modglobal::t_global
Definition: ModGlobal.F90:50
jbeg
**********************************************************************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 knode jbeg
Definition: RFLO_CalcGradDummySymm.h:32
rflo_modgridregionshape::rflo_blockskewedcell
subroutine rflo_blockskewedcell(region)
Definition: RFLO_ModGridRegionShape.F90:452
rflo_modgridregionshape::rflo_gridflatpatch
subroutine, public rflo_gridflatpatch(region, patch)
Definition: RFLO_ModGridRegionShape.F90:86
kbeg
**********************************************************************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 knode kbeg
Definition: RFLO_CalcGradDummySymm.h:32
moderror::deregisterfunction
subroutine deregisterfunction(global)
Definition: ModError.F90:469
rflo_getdimensphys
subroutine rflo_getdimensphys(region, iLev, ipcbeg, ipcend, jpcbeg, jpcend, kpcbeg, kpcend)
Definition: RFLO_GetDimensions.F90:109
modbndpatch::t_patch
Definition: ModBndPatch.F90:91