Rocstar  1.0
Rocstar multiphysics simulation application
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
RFLO_ModLaplaceSmoothing.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 Laplacian-smoothing routines.
26 !
27 ! Description: None.
28 !
29 ! ******************************************************************************
30 !
31 ! $Id: RFLO_ModLaplaceSmoothing.F90,v 1.13 2009/08/27 14:04:50 mtcampbe Exp $
32 !
33 ! Copyright: (c) 2004 by the University of Illinois
34 !
35 ! ******************************************************************************
36 
37 MODULE rflo_modlaplacesmoothing
38 
39  USE modglobal, ONLY : t_global
40  USE moddatastruct, ONLY: t_region
41  USE modgrid, ONLY : t_grid
42  USE modbndpatch, ONLY : t_patch
43  USE modparameters
44  USE moddatatypes
45  USE moderror
46  USE modmpi
47 
48  IMPLICIT NONE
49 
50  PRIVATE
51  PUBLIC :: rflo_laplacegridsmoo, &
52  rflo_laplacegridsolve, &
53  rflo_laplacegridpatch, &
54  rflo_laplacegridjump
55 
56 ! private :
57 ! RFLO_LaplaceGridOrtho
58 ! RFLO_ProjectQuadCorner
59 
60 ! ******************************************************************************
61 ! Declarations and definitions
62 ! ******************************************************************************
63 
64  CHARACTER(CHRLEN) :: RCSIdentString = &
65  '$RCSfile: RFLO_ModLaplaceSmoothing.F90,v $ $Revision: 1.13 $'
66 
67 ! ******************************************************************************
68 ! Routines
69 ! ******************************************************************************
70 
71  CONTAINS
72 
73 !******************************************************************************
74 !
75 ! Purpose: smooth the distribution of grid points by solving simplified
76 ! Laplace equation in physical space.
77 !
78 ! Description: none.
79 !
80 ! Input: regions = data of all grid regions.
81 !
82 ! Output: regions%levels%grid%xyz = new grid coordinates
83 ! resid = convergence of the Jacobi iteration.
84 !
85 ! Notes: none.
86 !
87 !******************************************************************************
88 
89 SUBROUTINE rflo_laplacegridsmoo( regions,resid )
90 
91  USE modinterfaces, ONLY : rflo_generatecoarsegrids, &
92  rflo_copygeometrydummy, rflo_extrapolategeometry, &
93  rflo_exchangegeometry, rflo_calcfacevectors, rflo_exchangednodecopy, &
94  rflo_exchangednodesend, rflo_exchangednoderecv, rflo_clearsendrequests,&
95  rflo_getdimensphysnodes, rflo_getnodeoffset, &
96  rflo_calccellcentroids, rflo_calcfacecentroids, rflo_changeinteriorgrid
97  IMPLICIT NONE
98 
99 #include "Indexing.h"
100 
101 ! ... parameters
102  REAL(RFREAL) :: resid
103 
104  TYPE(t_region), POINTER :: regions(:)
105 
106 ! ... loop variables
107  INTEGER :: ireg, ipatch, ijk, i, j, k, iter
108 
109 ! ... local variables
110  INTEGER :: ipnbeg, ipnend, jpnbeg, jpnend, kpnbeg, kpnend, inoff, ijnoff
111  INTEGER :: bctype, iregsrc, ipatchsrc
112 
113  REAL(RFREAL) :: dx, dy, dz
114  REAL(RFREAL), POINTER :: xyz(:,:), xyzold(:,:), xyztemp(:,:)
115 
116  TYPE(t_global), POINTER :: global
117  TYPE(t_grid), POINTER :: grid, gridold, gridsrc
118  TYPE(t_patch), POINTER :: patch, patchsrc
119 
120 !******************************************************************************
121 
122  global => regions(1)%global
123 
124  CALL registerfunction( global,'RFLO_LaplaceGridSmoo',&
125  'RFLO_ModLaplaceSmoothing.F90' )
126 
127 ! smooth grid region-wise -----------------------------------------------------
128 
129  resid = 0._rfreal
130 
131  DO ireg=1,global%nRegions
132  IF (regions(ireg)%procid==global%myProcid .AND. & ! region active and
133  regions(ireg)%active==active .AND. & ! on my processor
134  regions(ireg)%mixtInput%moveGrid) THEN ! and moving
135 
136  IF (global%moveGridScheme==movegrid_foms) THEN
137 
138 ! ----- enforce grid orthogonality
139  CALL rflo_laplacegridortho( regions(ireg) )
140 
141  ELSE
142 ! ----- compute movements in the interior and along the boundaries
143  CALL rflo_laplacegridsolve( regions(ireg) )
144 
145  ENDIF
146 
147 ! --- zero out movements along certain boundaries
148 
149  DO ipatch=1,regions(ireg)%nPatches
150  patch => regions(ireg)%levels(1)%patches(ipatch)
151  bctype = patch%bcType
152  IF ((bctype>=bc_inflow .AND. bctype<=bc_inflow +bc_range) .OR. &
153  (bctype>=bc_outflow .AND. bctype<=bc_outflow +bc_range) .OR. &
154  (bctype>=bc_slipwall .AND. bctype<=bc_slipwall +bc_range) .OR. &
155  (bctype>=bc_noslipwall .AND. bctype<=bc_noslipwall+bc_range) .OR. &
156  (bctype>=bc_farfield .AND. bctype<=bc_farfield +bc_range) .OR. &
157  (bctype>=bc_injection .AND. bctype<=bc_injection +bc_range) .OR. &
158  (bctype>=bc_tra_peri .AND. bctype<=bc_tra_peri +bc_range) .OR. &
159  (bctype>=bc_rot_peri .AND. bctype<=bc_rot_peri +bc_range) .OR. &
160  (bctype>=bc_symmetry_free .AND. bctype<=bc_symmetry +bc_range)) THEN
161  CALL rflo_laplacegridpatch( regions(ireg),patch )
162  ENDIF ! bcType
163  ENDDO ! iPatch
164 
165  CALL rflo_getdimensphysnodes( regions(ireg),1,ipnbeg,ipnend, &
166  jpnbeg,jpnend,kpnbeg,kpnend )
167  CALL rflo_getnodeoffset( regions(ireg),1,inoff,ijnoff )
168 
169  xyz => regions(ireg)%levels(1)%grid%xyz
170  xyzold => regions(ireg)%levels(1)%grid%xyzOld
171 
172  DO k=kpnbeg,kpnend
173  DO j=jpnbeg,jpnend
174  DO i=ipnbeg,ipnend
175  ijk = indijk(i,j,k,inoff,ijnoff)
176  dx = xyz(xcoord,ijk) - xyzold(xcoord,ijk)
177  dy = xyz(ycoord,ijk) - xyzold(ycoord,ijk)
178  dz = xyz(zcoord,ijk) - xyzold(zcoord,ijk)
179  resid = resid + dx*dx + dy*dy +dz*dz
180  ENDDO
181  ENDDO
182  ENDDO
183 
184  ENDIF ! region on this processor and active, grid moving
185  ENDDO ! iReg
186 
187 ! fix interfaces between regions ----------------------------------------------
188 ! copy / send deformations
189 
190  DO ireg=1,global%nRegions
191  IF (regions(ireg)%procid==global%myProcid .AND. & ! region active and
192  regions(ireg)%active==active .AND. & ! on my processor
193  regions(ireg)%mixtInput%moveGrid) THEN ! and moving
194 
195  grid => regions(ireg)%levels(1)%grid
196  gridold => regions(ireg)%levels(1)%gridOld
197 
198  DO ipatch=1,regions(ireg)%nPatches
199  patch => regions(ireg)%levels(1)%patches(ipatch)
200  bctype = patch%bcType
201  IF ((bctype>=bc_regionconf .AND. bctype<=bc_regionconf+bc_range) .OR. &
202  (bctype>=bc_tra_peri .AND. bctype<=bc_tra_peri +bc_range) .OR. &
203  (bctype>=bc_rot_peri .AND. bctype<=bc_rot_peri +bc_range)) THEN
204  iregsrc = patch%srcRegion
205  ipatchsrc = patch%srcPatch
206  patchsrc => regions(iregsrc)%levels(1)%patches(ipatchsrc)
207  gridsrc => regions(iregsrc)%levels(1)%grid
208 
209  IF (regions(iregsrc)%procid == global%myProcid) THEN
210  CALL rflo_exchangednodecopy( regions(ireg),regions(iregsrc), &
211  patch,patchsrc,.true., &
212  grid%xyz,gridsrc%xyz )
213  ELSE
214  CALL rflo_exchangednodesend( regions(ireg),regions(iregsrc), &
215  patch,grid%xyz )
216  ENDIF
217  ENDIF ! bcType
218  ENDDO ! iPatch
219 
220  ENDIF ! region on this processor and active, grid moving
221  ENDDO ! iReg
222 
223 ! receive deformations
224 
225  DO ireg=1,global%nRegions
226  IF (regions(ireg)%procid==global%myProcid .AND. & ! region active and
227  regions(ireg)%active==active .AND. & ! on my processor
228  regions(ireg)%mixtInput%moveGrid) THEN ! and moving
229 
230  grid => regions(ireg)%levels(1)%grid
231  gridold => regions(ireg)%levels(1)%gridOld
232 
233  DO ipatch=1,regions(ireg)%nPatches
234  patch => regions(ireg)%levels(1)%patches(ipatch)
235  bctype = patch%bcType
236  IF ((bctype>=bc_regionconf .AND. bctype<=bc_regionconf+bc_range) .OR. &
237  (bctype>=bc_tra_peri .AND. bctype<=bc_tra_peri +bc_range) .OR. &
238  (bctype>=bc_rot_peri .AND. bctype<=bc_rot_peri +bc_range)) THEN
239  iregsrc = patch%srcRegion
240  ipatchsrc = patch%srcPatch
241  patchsrc => regions(iregsrc)%levels(1)%patches(ipatchsrc)
242  gridsrc => regions(iregsrc)%levels(1)%grid
243 
244  IF (regions(iregsrc)%procid /= global%myProcid) THEN
245  CALL rflo_exchangednoderecv( regions(ireg),regions(iregsrc), &
246  patch,patchsrc,.true.,grid%xyz )
247  ENDIF
248  ENDIF ! bcType
249  ENDDO ! iPatch
250 
251  ENDIF ! region on this processor and active, grid moving
252  ENDDO ! iReg
253 
254 ! clear send requests
255 
256  DO ireg=1,global%nRegions
257  IF (regions(ireg)%procid==global%myProcid .AND. & ! region active and
258  regions(ireg)%active==active .AND. & ! on my processor
259  regions(ireg)%mixtInput%moveGrid) THEN ! and moving
260  CALL rflo_clearsendrequests( regions,ireg,.true. )
261  ENDIF
262  ENDDO
263 
264 ! update grid, dummy, corner and edge cells -----------------------------------
265 
266  DO ireg=1,global%nRegions
267  IF (regions(ireg)%procid==global%myProcid .AND. & ! region active and
268  regions(ireg)%active==active .AND. & ! on my processor
269  regions(ireg)%mixtInput%moveGrid) THEN ! and moving
270 
271 ! --- change xyz from deformations to coordinates
272 
273  xyz => regions(ireg)%levels(1)%grid%xyz
274  xyzold => regions(ireg)%levels(1)%gridOld%xyz
275 
276  DO ijk=lbound(xyz,2),ubound(xyz,2)
277  xyz(xcoord,ijk) = xyz(xcoord,ijk) + xyzold(xcoord,ijk)
278  xyz(ycoord,ijk) = xyz(ycoord,ijk) + xyzold(ycoord,ijk)
279  xyz(zcoord,ijk) = xyz(zcoord,ijk) + xyzold(zcoord,ijk)
280  ENDDO
281 
282 ! --- update coarse grids and dummy cells
283 
284  CALL rflo_generatecoarsegrids( regions(ireg) ) ! coarsen finest grid
285  CALL rflo_copygeometrydummy( regions(ireg) ) ! copy to dummy nodes
286  CALL rflo_extrapolategeometry( regions(ireg) ) ! extrapolate
287  ENDIF ! region on this processor and active, grid moving
288  ENDDO ! iReg
289 
290  CALL rflo_exchangegeometry( regions ) ! exchange geometry
291 
292 ! recompute cell and face centroids for FOMS grid scheme ----------------------
293 
294 ! IF (global%moveGridScheme/=MOVEGRID_FOMS) GOTO 777
295  goto 777
296 
297  DO ireg=1,global%nRegions
298  IF (regions(ireg)%procid==global%myProcid .AND. & ! region active and
299  regions(ireg)%active==active .AND. & ! on my processor
300  regions(ireg)%mixtInput%moveGrid) THEN ! and moving
301  CALL rflo_calccellcentroids( regions(ireg) ) ! cell centroids
302  CALL rflo_calcfacecentroids( regions(ireg) ) ! face centroids
303  CALL rflo_calcfacevectors( regions(ireg) )
304  ENDIF ! region on this processor and active, grid moving
305  ENDDO ! iReg
306 
307  DO iter=1,1
308 
309  DO ireg=1,global%nRegions
310  IF (regions(ireg)%procid==global%myProcid .AND. & ! region active and
311  regions(ireg)%active==active .AND. & ! on my processor
312  regions(ireg)%mixtInput%moveGrid) THEN ! and moving
313 
314 ! --- enforce grid orthogonality
315 
316  CALL rflo_laplacegridortho( regions(ireg) )
317 
318 ! --- zero out movements and recored jumps along certain boundaries
319 
320 ! regions(iReg)%levels(1)%grid%xyzTemp = 0._RFREAL
321 
322  DO ipatch=1,regions(ireg)%nPatches
323  patch => regions(ireg)%levels(1)%patches(ipatch)
324  bctype = patch%bcType
325  IF ((bctype>=bc_inflow .AND. bctype<=bc_inflow +bc_range).OR. &
326  (bctype>=bc_outflow .AND. bctype<=bc_outflow +bc_range).OR. &
327  (bctype>=bc_slipwall .AND. bctype<=bc_slipwall +bc_range).OR. &
328  (bctype>=bc_noslipwall .AND. bctype<=bc_noslipwall+bc_range).OR. &
329  (bctype>=bc_farfield .AND. bctype<=bc_farfield +bc_range).OR. &
330  (bctype>=bc_injection .AND. bctype<=bc_injection +bc_range).OR. &
331  (bctype>=bc_tra_peri .AND. bctype<=bc_tra_peri +bc_range).OR. &
332  (bctype>=bc_rot_peri .AND. bctype<=bc_rot_peri +bc_range).OR. &
333  (bctype>=bc_symmetry_free .AND. bctype<=bc_symmetry +bc_range)) THEN
334  CALL rflo_laplacegridpatch( regions(ireg),patch )
335  ENDIF
336 ! CALL RFLO_LaplaceGridJump( regions(iReg),patch )
337  ENDDO ! iPatch
338 
339 ! --- change the interior movements
340 
341 ! grid => regions(iReg)%levels(1)%grid
342 ! gridOld => regions(iReg)%levels(1)%gridOld
343 
344 ! CALL RFLO_ChangeInteriorGrid( regions(iReg),grid%boundMoved, &
345 ! grid%edgeMoved,grid%arcLen12, &
346 ! grid%arcLen34,grid%arcLen56, &
347 ! gridOld%xyzOld,grid%xyzTemp )
348 
349 ! CALL RFLO_GetDimensPhysNodes( regions(iReg),1,ipnbeg,ipnend, &
350 ! jpnbeg,jpnend,kpnbeg,kpnend )
351 ! CALL RFLO_GetNodeOffset( regions(iReg),1,iNOff,ijNOff )
352 !IF (iReg==13) write(*,*)'ORTHO1',grid%xyz(:,IndIJK(10,10,10,iNOff,ijNOff))+ gridOld%xyz(:,IndIJK(10,10,10,iNOff,ijNOff)), grid%xyzOrth(:,IndIJK(10,10,10,iNOff,ijNOff))+ gridOld%xyz(:,IndIJK(10,10,10,iNOff,ijNOff))
353 
354 ! DO k=kpnbeg,kpnend
355 ! DO j=jpnbeg,jpnend
356 ! DO i=ipnbeg,ipnend
357 ! ijk = IndIJK(i,j,k,iNOff,ijNOff)
358 ! grid%xyz(:,ijk) = grid%xyzTemp(:,ijk) - &
359 ! gridOld%xyz(:,ijk) + grid%xyzOrth(:,ijk)
360 ! ENDDO
361 ! ENDDO
362 ! ENDDO
363 !IF (iReg==13) write(*,*)'ORTHO2',grid%xyz(:,IndIJK(10,10,10,iNOff,ijNOff)),grid%xyzOrth(:,IndIJK(10,10,10,iNOff,ijNOff))
364 
365 ! --- compute residuals
366 
367  CALL rflo_getdimensphysnodes( regions(ireg),1,ipnbeg,ipnend, &
368  jpnbeg,jpnend,kpnbeg,kpnend )
369  CALL rflo_getnodeoffset( regions(ireg),1,inoff,ijnoff )
370 
371  xyz => regions(ireg)%levels(1)%grid%xyz
372  xyzold => regions(ireg)%levels(1)%grid%xyzOld
373 
374  DO k=kpnbeg,kpnend
375  DO j=jpnbeg,jpnend
376  DO i=ipnbeg,ipnend
377  ijk = indijk(i,j,k,inoff,ijnoff)
378  dx = xyz(xcoord,ijk) - xyzold(xcoord,ijk)
379  dy = xyz(ycoord,ijk) - xyzold(ycoord,ijk)
380  dz = xyz(zcoord,ijk) - xyzold(zcoord,ijk)
381  resid = resid + dx*dx + dy*dy +dz*dz
382  ENDDO
383  ENDDO
384  ENDDO
385 
386  ENDIF ! region on this processor and active, grid moving
387  ENDDO ! iReg
388 
389 ! fix interfaces between regions ----------------------------------------------
390 ! copy / send deformations
391 
392  DO ireg=1,global%nRegions
393  IF (regions(ireg)%procid==global%myProcid .AND. & ! region active and
394  regions(ireg)%active==active .AND. & ! on my processor
395  regions(ireg)%mixtInput%moveGrid) THEN ! and moving
396 
397  grid => regions(ireg)%levels(1)%grid
398  gridold => regions(ireg)%levels(1)%gridOld
399 
400  DO ipatch=1,regions(ireg)%nPatches
401  patch => regions(ireg)%levels(1)%patches(ipatch)
402  bctype = patch%bcType
403  IF ((bctype>=bc_regionconf .AND. bctype<=bc_regionconf+bc_range) .OR. &
404  (bctype>=bc_tra_peri .AND. bctype<=bc_tra_peri +bc_range) .OR. &
405  (bctype>=bc_rot_peri .AND. bctype<=bc_rot_peri +bc_range)) THEN
406  iregsrc = patch%srcRegion
407  ipatchsrc = patch%srcPatch
408  patchsrc => regions(iregsrc)%levels(1)%patches(ipatchsrc)
409  gridsrc => regions(iregsrc)%levels(1)%grid
410 
411  IF (regions(iregsrc)%procid == global%myProcid) THEN
412  CALL rflo_exchangednodecopy( regions(ireg),regions(iregsrc), &
413  patch,patchsrc,.true., &
414  grid%xyz,gridsrc%xyz )
415  ELSE
416  CALL rflo_exchangednodesend( regions(ireg),regions(iregsrc), &
417  patch,grid%xyz )
418  ENDIF
419  ENDIF ! bcType
420  ENDDO ! iPatch
421 
422  ENDIF ! region on this processor and active, grid moving
423  ENDDO ! iReg
424 
425 ! receive deformations
426 
427  DO ireg=1,global%nRegions
428  IF (regions(ireg)%procid==global%myProcid .AND. & ! region active and
429  regions(ireg)%active==active .AND. & ! on my processor
430  regions(ireg)%mixtInput%moveGrid) THEN ! and moving
431 
432  grid => regions(ireg)%levels(1)%grid
433  gridold => regions(ireg)%levels(1)%gridOld
434 
435  DO ipatch=1,regions(ireg)%nPatches
436  patch => regions(ireg)%levels(1)%patches(ipatch)
437  bctype = patch%bcType
438  IF ((bctype>=bc_regionconf .AND. bctype<=bc_regionconf+bc_range) .OR. &
439  (bctype>=bc_tra_peri .AND. bctype<=bc_tra_peri +bc_range) .OR. &
440  (bctype>=bc_rot_peri .AND. bctype<=bc_rot_peri +bc_range)) THEN
441  iregsrc = patch%srcRegion
442  ipatchsrc = patch%srcPatch
443  patchsrc => regions(iregsrc)%levels(1)%patches(ipatchsrc)
444  gridsrc => regions(iregsrc)%levels(1)%grid
445 
446  IF (regions(iregsrc)%procid /= global%myProcid) THEN
447  CALL rflo_exchangednoderecv( regions(ireg),regions(iregsrc), &
448  patch,patchsrc,.true.,grid%xyz )
449  ENDIF
450  ENDIF ! bcType
451  ENDDO ! iPatch
452 
453  ENDIF ! region on this processor and active, grid moving
454  ENDDO ! iReg
455 
456 ! clear send requests
457 
458  DO ireg=1,global%nRegions
459  IF (regions(ireg)%procid==global%myProcid .AND. & ! region active and
460  regions(ireg)%active==active .AND. & ! on my processor
461  regions(ireg)%mixtInput%moveGrid) THEN ! and moving
462  CALL rflo_clearsendrequests( regions,ireg,.true. )
463  ENDIF
464  ENDDO
465 
466 ! update grid, dummy, corner and edge cells -----------------------------------
467 
468  DO ireg=1,global%nRegions
469  IF (regions(ireg)%procid==global%myProcid .AND. & ! region active and
470  regions(ireg)%active==active .AND. & ! on my processor
471  regions(ireg)%mixtInput%moveGrid) THEN ! and moving
472 
473 ! --- change xyz from deformations to coordinates
474 
475  xyz => regions(ireg)%levels(1)%grid%xyz
476  xyzold => regions(ireg)%levels(1)%gridOld%xyz
477 
478  DO ijk=lbound(xyz,2),ubound(xyz,2)
479  xyz(xcoord,ijk) = xyz(xcoord,ijk) + xyzold(xcoord,ijk)
480  xyz(ycoord,ijk) = xyz(ycoord,ijk) + xyzold(ycoord,ijk)
481  xyz(zcoord,ijk) = xyz(zcoord,ijk) + xyzold(zcoord,ijk)
482  ENDDO
483 
484 ! --- update coarse grids and dummy cells
485 
486  CALL rflo_generatecoarsegrids( regions(ireg) ) ! coarsen finest grid
487  CALL rflo_copygeometrydummy( regions(ireg) ) ! copy to dummy nodes
488  CALL rflo_extrapolategeometry( regions(ireg) ) ! extrapolate
489  ENDIF ! region on this processor and active, grid moving
490  ENDDO ! iReg
491 
492  CALL rflo_exchangegeometry( regions ) ! exchange geometry
493 
494  ENDDO ! iter
495 
496 ! finalize --------------------------------------------------------------------
497 
498 777 CONTINUE
499 
500  CALL deregisterfunction( global )
501 
502 END SUBROUTINE rflo_laplacegridsmoo
503 
504 
505 !******************************************************************************
506 !
507 ! Purpose: zero out movements obtained by LaplaceGridSolve along a boundary.
508 !
509 ! Description: none.
510 !
511 ! Input: region = data of current region, grid movements
512 ! patch = current patch.
513 !
514 ! Output: region%levels%grid%xyz = new grid movements.
515 !
516 ! Notes: none.
517 !
518 !******************************************************************************
519 
520 SUBROUTINE rflo_laplacegridpatch( region,patch )
521 
522  USE modinterfaces, ONLY : rflo_getpatchindicesnodes, rflo_getnodeoffset
523 
524  IMPLICIT NONE
525 
526 #include "Indexing.h"
527 
528 ! ... parameters
529  TYPE(t_region) :: region
530  TYPE(t_patch) :: patch
531 
532 ! ... loop variables
533  INTEGER :: i, j, k
534 
535 ! ... local variables
536  INTEGER :: ilev
537  INTEGER :: ibeg, iend, jbeg, jend, kbeg, kend, inoff, ijnoff, ijknb
538 
539  REAL(RFREAL), POINTER :: xyz(:,:), xyzold(:,:)
540 
541 !******************************************************************************
542 
543  CALL registerfunction( region%global,'RFLO_LaplaceGridPatch',&
544  'RFLO_ModLaplaceSmoothing.F90' )
545 
546 ! get dimensions and pointers
547 
548  ilev = 1
549 
550  CALL rflo_getpatchindicesnodes( region,patch,ilev,ibeg,iend, &
551  jbeg,jend,kbeg,kend )
552  CALL rflo_getnodeoffset( region,ilev,inoff,ijnoff )
553 
554  xyz => region%levels(ilev)%grid%xyz
555  xyzold => region%levels(ilev)%grid%xyzOld
556 
557 ! new = old
558 
559  DO k=kbeg,kend
560  DO j=jbeg,jend
561  DO i=ibeg,iend
562  ijknb = indijk(i,j,k,inoff,ijnoff)
563  xyz(xcoord,ijknb) = xyzold(xcoord,ijknb)
564  xyz(ycoord,ijknb) = xyzold(ycoord,ijknb)
565  xyz(zcoord,ijknb) = xyzold(zcoord,ijknb)
566  ENDDO
567  ENDDO
568  ENDDO
569 
570 ! finalize
571 
572  CALL deregisterfunction( region%global )
573 
574 END SUBROUTINE rflo_laplacegridpatch
575 
576 !******************************************************************************
577 !
578 ! Purpose: obtain jump of surface grid motion due to enforcing back boundary
579 ! surface motion.
580 !
581 ! Description: none.
582 !
583 ! Input: region = data of current region, grid movements
584 ! patch = current patch.
585 !
586 ! Output: region%levels%grid%xyz = new grid movements.
587 !
588 ! Notes: none.
589 !
590 !******************************************************************************
591 
592 SUBROUTINE rflo_laplacegridjump( region,patch )
593 
594  USE modinterfaces, ONLY : rflo_getpatchindicesnodes, rflo_getnodeoffset
595 
596  IMPLICIT NONE
597 
598 #include "Indexing.h"
599 
600 ! ... parameters
601  TYPE(t_region) :: region
602  TYPE(t_patch) :: patch
603 
604 ! ... loop variables
605  INTEGER :: i, j, k
606 
607 ! ... local variables
608  INTEGER :: ilev
609  INTEGER :: ibeg, iend, jbeg, jend, kbeg, kend, inoff, ijnoff, ijknb
610 
611  REAL(RFREAL), POINTER :: xyz(:,:), xyzorth(:,:), xyztemp(:,:)
612 
613 !******************************************************************************
614 
615  CALL registerfunction( region%global,'RFLO_LaplaceGridJump',&
616  'RFLO_ModLaplaceSmoothing.F90' )
617 
618 ! get dimensions and pointers
619 
620  ilev = 1
621 
622  CALL rflo_getpatchindicesnodes( region,patch,ilev,ibeg,iend, &
623  jbeg,jend,kbeg,kend )
624  CALL rflo_getnodeoffset( region,ilev,inoff,ijnoff )
625 
626  xyz => region%levels(ilev)%grid%xyz
627  xyzorth => region%levels(ilev)%grid%xyzOrth
628  xyztemp => region%levels(ilev)%grid%xyzTemp
629 
630 ! jump = new-(orth-orig)
631 
632  DO k=kbeg,kend
633  DO j=jbeg,jend
634  DO i=ibeg,iend
635  ijknb = indijk(i,j,k,inoff,ijnoff)
636  xyztemp(xcoord,ijknb) = xyz(xcoord,ijknb) - xyzorth(xcoord,ijknb)
637  xyztemp(ycoord,ijknb) = xyz(ycoord,ijknb) - xyzorth(ycoord,ijknb)
638  xyztemp(zcoord,ijknb) = xyz(zcoord,ijknb) - xyzorth(zcoord,ijknb)
639  ENDDO
640  ENDDO
641  ENDDO
642 
643 ! finalize
644 
645  CALL deregisterfunction( region%global )
646 
647 END SUBROUTINE rflo_laplacegridjump
648 
649 
650 !******************************************************************************
651 !
652 ! Purpose: conduct one Jacobi iteration to obtain new grid movements
653 ! in the flow domain (boundaries included).
654 !
655 ! Description: none.
656 !
657 ! Input: region = data of current region, old grid coordinates.
658 !
659 ! Output: region%levels%grid%xyz = grid movements.
660 !
661 ! Notes: on entry, xyz holds node coordinates from a previous smoothing
662 ! step. On exit however, xyz contains only the grid motion.
663 !
664 !******************************************************************************
665 
666 SUBROUTINE rflo_laplacegridsolve( region )
667 
668  USE modinterfaces, ONLY : rflo_getdimensdummynodes, rflo_getdimensphysnodes, &
669  rflo_getnodeoffset
670  IMPLICIT NONE
671 
672 #include "Indexing.h"
673 
674 ! ... parameters
675  TYPE(t_region) :: region
676 
677 ! ... loop variables
678  INTEGER :: i, j, k
679 
680 ! ... local variables
681  INTEGER :: idnbeg, idnend, jdnbeg, jdnend, kdnbeg, kdnend, ilev, inoff, ijnoff
682  INTEGER :: ipnbeg, ipnend, jpnbeg, jpnend, kpnbeg, kpnend, ibn, ien
683  INTEGER :: ijk, ip1, ip2, im1, im2, jp1, jp2, jm1, jm2, kp1, kp2, km1, km2
684  INTEGER :: moveblock, method, in, ico, ic(8)
685  INTEGER, POINTER :: idgen(:)
686 
687  LOGICAL, POINTER :: bndmoved(:)
688 
689  REAL(RFREAL) :: rxi, ryi, rzi, rxj, ryj, rzj, rxk, ryk, rzk, wi, wj, wk, &
690  si, sj, sk, sim, sjm, skm, sip, sjp, skp, d, p, eps
691  REAL(RFREAL) :: denom(8), dist(8,8)
692  REAL(RFREAL), POINTER :: xyz(:,:), xyzold(:,:), xyzorig(:,:)
693 
694 !******************************************************************************
695 
696  CALL registerfunction( region%global,'RFLO_LaplaceGridSolve',&
697  'RFLO_ModLaplaceSmoothing.F90' )
698 
699 ! get dimensions and pointers -------------------------------------------------
700 
701  ilev = 1
702 
703  CALL rflo_getdimensdummynodes( region,ilev,idnbeg,idnend, &
704  jdnbeg,jdnend,kdnbeg,kdnend )
705  CALL rflo_getdimensphysnodes( region,ilev,ipnbeg,ipnend,jpnbeg,jpnend, &
706  kpnbeg,kpnend )
707  CALL rflo_getnodeoffset( region,ilev,inoff,ijnoff )
708  ibn = indijk(idnbeg,jdnbeg,kdnbeg,inoff,ijnoff)
709  ien = indijk(idnend,jdnend,kdnend,inoff,ijnoff)
710 
711  xyz => region%levels(ilev)%grid%xyz
712  xyzold => region%levels(ilev)%grid%xyzOld
713  xyzorig => region%levels(ilev)%gridOld%xyz
714  bndmoved => region%levels(ilev)%grid%boundMoved
715  idgen => region%levels(ilev)%grid%ijkDgen
716 
717  p = region%global%moveGridPower
718  eps = epsilon( 1._rfreal )
719 
720 ! reset motion vectors --------------------------------------------------------
721 
722  DO ijk=ibn,ien
723  xyzold(xcoord,ijk) = xyz(xcoord,ijk) - xyzorig(xcoord,ijk)
724  xyzold(ycoord,ijk) = xyz(ycoord,ijk) - xyzorig(ycoord,ijk)
725  xyzold(zcoord,ijk) = xyz(zcoord,ijk) - xyzorig(zcoord,ijk)
726  ENDDO
727 
728 ! move block corners locally
729 
730  moveblock = 0
731  method = 1
732 
733  IF (moveblock==1) THEN
734  ic(1) = indijk(ipnbeg ,jpnbeg ,kpnbeg ,inoff,ijnoff)
735  ic(2) = indijk(ipnbeg ,jpnbeg ,kpnend ,inoff,ijnoff)
736  ic(3) = indijk(ipnbeg ,jpnend ,kpnend ,inoff,ijnoff)
737  ic(4) = indijk(ipnbeg ,jpnend ,kpnbeg ,inoff,ijnoff)
738  ic(5) = indijk(ipnend ,jpnbeg ,kpnbeg ,inoff,ijnoff)
739  ic(6) = indijk(ipnend ,jpnbeg ,kpnend ,inoff,ijnoff)
740  ic(7) = indijk(ipnend ,jpnend ,kpnend ,inoff,ijnoff)
741  ic(8) = indijk(ipnend ,jpnend ,kpnbeg ,inoff,ijnoff)
742 
743  DO ico = 1,8
744  denom(ico) = 0._rfreal
745  DO in = 1,8
746  IF (ico/=in) THEN
747  dist(in,ico) = sqrt( (xyz(xcoord,ic(in))-xyz(xcoord,ic(ico)))**2 + &
748  (xyz(ycoord,ic(in))-xyz(ycoord,ic(ico)))**2 + &
749  (xyz(zcoord,ic(in))-xyz(zcoord,ic(ico)))**2 )
750  dist(in,ico) = 1._rfreal/dist(in,ico)
751  denom(ico) = denom(ico) + dist(in,ico)
752  ENDIF
753  ENDDO
754  ENDDO
755 
756  IF (.NOT. bndmoved(1)) THEN
757  IF ((.NOT. bndmoved(3)).AND.(.NOT. bndmoved(5))) THEN
758  xyz(xcoord,ic(1)) = 0._rfreal
759  xyz(ycoord,ic(1)) = 0._rfreal
760  xyz(zcoord,ic(1)) = 0._rfreal
761  DO in=1,8
762  IF (in/=1) THEN
763  xyz(xcoord,ic(1))= xyz(xcoord,ic(1))+dist(in,1)*xyzold(xcoord,ic(in))
764  xyz(ycoord,ic(1))= xyz(ycoord,ic(1))+dist(in,1)*xyzold(ycoord,ic(in))
765  xyz(zcoord,ic(1))= xyz(zcoord,ic(1))+dist(in,1)*xyzold(zcoord,ic(in))
766  ENDIF
767  ENDDO
768  xyzold(xcoord,ic(1)) = xyz(xcoord,ic(1))/denom(1)
769  xyzold(ycoord,ic(1)) = xyz(ycoord,ic(1))/denom(1)
770  xyzold(zcoord,ic(1)) = xyz(zcoord,ic(1))/denom(1)
771  ENDIF
772  IF ((.NOT. bndmoved(3)).AND.(.NOT. bndmoved(6))) THEN
773  xyz(xcoord,ic(2)) = 0._rfreal
774  xyz(ycoord,ic(2)) = 0._rfreal
775  xyz(zcoord,ic(2)) = 0._rfreal
776  DO in=1,8
777  IF (in/=2) THEN
778  xyz(xcoord,ic(2))= xyz(xcoord,ic(2))+dist(in,2)*xyzold(xcoord,ic(in))
779  xyz(ycoord,ic(2))= xyz(ycoord,ic(2))+dist(in,2)*xyzold(ycoord,ic(in))
780  xyz(zcoord,ic(2))= xyz(zcoord,ic(2))+dist(in,2)*xyzold(zcoord,ic(in))
781  ENDIF
782  ENDDO
783  xyzold(xcoord,ic(2)) = xyz(xcoord,ic(2))/denom(2)
784  xyzold(ycoord,ic(2)) = xyz(ycoord,ic(2))/denom(2)
785  xyzold(zcoord,ic(2)) = xyz(zcoord,ic(2))/denom(2)
786  ENDIF
787  IF ((.NOT. bndmoved(4)).AND.(.NOT. bndmoved(6))) THEN
788  xyz(xcoord,ic(3)) = 0._rfreal
789  xyz(ycoord,ic(3)) = 0._rfreal
790  xyz(zcoord,ic(3)) = 0._rfreal
791  DO in=1,8
792  IF (in/=3) THEN
793  xyz(xcoord,ic(3))= xyz(xcoord,ic(3))+dist(in,3)*xyzold(xcoord,ic(in))
794  xyz(ycoord,ic(3))= xyz(ycoord,ic(3))+dist(in,3)*xyzold(ycoord,ic(in))
795  xyz(zcoord,ic(3))= xyz(zcoord,ic(3))+dist(in,3)*xyzold(zcoord,ic(in))
796  ENDIF
797  ENDDO
798  xyzold(xcoord,ic(3)) = xyz(xcoord,ic(3))/denom(3)
799  xyzold(ycoord,ic(3)) = xyz(ycoord,ic(3))/denom(3)
800  xyzold(zcoord,ic(3)) = xyz(zcoord,ic(3))/denom(3)
801  ENDIF
802  IF ((.NOT. bndmoved(4)).AND.(.NOT. bndmoved(5))) THEN
803  xyz(xcoord,ic(4)) = 0._rfreal
804  xyz(ycoord,ic(4)) = 0._rfreal
805  xyz(zcoord,ic(4)) = 0._rfreal
806  DO in=1,8
807  IF (in/=4) THEN
808  xyz(xcoord,ic(4))= xyz(xcoord,ic(4))+dist(in,4)*xyzold(xcoord,ic(in))
809  xyz(ycoord,ic(4))= xyz(ycoord,ic(4))+dist(in,4)*xyzold(ycoord,ic(in))
810  xyz(zcoord,ic(4))= xyz(zcoord,ic(4))+dist(in,4)*xyzold(zcoord,ic(in))
811  ENDIF
812  ENDDO
813  xyzold(xcoord,ic(4)) = xyz(xcoord,ic(4))/denom(4)
814  xyzold(ycoord,ic(4)) = xyz(ycoord,ic(4))/denom(4)
815  xyzold(zcoord,ic(4)) = xyz(zcoord,ic(4))/denom(4)
816  ENDIF
817  ENDIF
818 
819  IF (.NOT. bndmoved(2)) THEN
820  IF ((.NOT. bndmoved(3)).AND.(.NOT. bndmoved(5))) THEN
821  xyz(xcoord,ic(5)) = 0._rfreal
822  xyz(ycoord,ic(5)) = 0._rfreal
823  xyz(zcoord,ic(5)) = 0._rfreal
824  DO in=1,8
825  IF (in/=5) THEN
826  xyz(xcoord,ic(5))= xyz(xcoord,ic(5))+dist(in,5)*xyzold(xcoord,ic(in))
827  xyz(ycoord,ic(5))= xyz(ycoord,ic(5))+dist(in,5)*xyzold(ycoord,ic(in))
828  xyz(zcoord,ic(5))= xyz(zcoord,ic(5))+dist(in,5)*xyzold(zcoord,ic(in))
829  ENDIF
830  ENDDO
831  xyzold(xcoord,ic(5)) = xyz(xcoord,ic(5))/denom(5)
832  xyzold(ycoord,ic(5)) = xyz(ycoord,ic(5))/denom(5)
833  xyzold(zcoord,ic(5)) = xyz(zcoord,ic(5))/denom(5)
834  ENDIF
835  IF ((.NOT. bndmoved(3)).AND.(.NOT. bndmoved(6))) THEN
836  xyz(xcoord,ic(6)) = 0._rfreal
837  xyz(ycoord,ic(6)) = 0._rfreal
838  xyz(zcoord,ic(6)) = 0._rfreal
839  DO in=1,8
840  IF (in/=6) THEN
841  xyz(xcoord,ic(6))= xyz(xcoord,ic(6))+dist(in,6)*xyzold(xcoord,ic(in))
842  xyz(ycoord,ic(6))= xyz(ycoord,ic(6))+dist(in,6)*xyzold(ycoord,ic(in))
843  xyz(zcoord,ic(6))= xyz(zcoord,ic(6))+dist(in,6)*xyzold(zcoord,ic(in))
844  ENDIF
845  ENDDO
846  xyzold(xcoord,ic(6)) = xyz(xcoord,ic(6))/denom(6)
847  xyzold(ycoord,ic(6)) = xyz(ycoord,ic(6))/denom(6)
848  xyzold(zcoord,ic(6)) = xyz(zcoord,ic(6))/denom(6)
849  ENDIF
850  IF ((.NOT. bndmoved(4)).AND.(.NOT. bndmoved(6))) THEN
851  xyz(xcoord,ic(7)) = 0._rfreal
852  xyz(ycoord,ic(7)) = 0._rfreal
853  xyz(zcoord,ic(7)) = 0._rfreal
854  DO in=1,8
855  IF (in/=7) THEN
856  xyz(xcoord,ic(7))= xyz(xcoord,ic(7))+dist(in,7)*xyzold(xcoord,ic(in))
857  xyz(ycoord,ic(7))= xyz(ycoord,ic(7))+dist(in,7)*xyzold(ycoord,ic(in))
858  xyz(zcoord,ic(7))= xyz(zcoord,ic(7))+dist(in,7)*xyzold(zcoord,ic(in))
859  ENDIF
860  ENDDO
861  xyzold(xcoord,ic(7)) = xyz(xcoord,ic(7))/denom(7)
862  xyzold(ycoord,ic(7)) = xyz(ycoord,ic(7))/denom(7)
863  xyzold(zcoord,ic(7)) = xyz(zcoord,ic(7))/denom(7)
864  ENDIF
865  IF ((.NOT. bndmoved(4)).AND.(.NOT. bndmoved(5))) THEN
866  xyz(xcoord,ic(8)) = 0._rfreal
867  xyz(ycoord,ic(8)) = 0._rfreal
868  xyz(zcoord,ic(8)) = 0._rfreal
869  DO in=1,8
870  IF (in/=8) THEN
871  xyz(xcoord,ic(8))= xyz(xcoord,ic(8))+dist(in,8)*xyzold(xcoord,ic(in))
872  xyz(ycoord,ic(8))= xyz(ycoord,ic(8))+dist(in,8)*xyzold(ycoord,ic(in))
873  xyz(zcoord,ic(8))= xyz(zcoord,ic(8))+dist(in,8)*xyzold(zcoord,ic(in))
874  ENDIF
875  ENDDO
876  xyzold(xcoord,ic(8)) = xyz(xcoord,ic(8))/denom(8)
877  xyzold(ycoord,ic(8)) = xyz(ycoord,ic(8))/denom(8)
878  xyzold(zcoord,ic(8)) = xyz(zcoord,ic(8))/denom(8)
879  ENDIF
880  ENDIF ! bndMoved
881  ENDIF ! moveBlock corners
882 
883 ! compute new coordinates -----------------------------------------------------
884 
885  IF (method==1) THEN
886 
887 ! - New weighted Laplacian Smoothing of Bono Wasistho:
888 
889  DO k=kpnbeg,kpnend
890  DO j=jpnbeg,jpnend
891  DO i=ipnbeg,ipnend
892  ijk = indijk(i ,j ,k ,inoff,ijnoff)
893  ip1 = indijk(i+1,j ,k ,inoff,ijnoff)
894  ip2 = indijk(i+2,j ,k ,inoff,ijnoff)
895  im1 = indijk(i-1,j ,k ,inoff,ijnoff)
896  im2 = indijk(i-2,j ,k ,inoff,ijnoff)
897  jp1 = indijk(i ,j+1,k ,inoff,ijnoff)
898  jp2 = indijk(i ,j+2,k ,inoff,ijnoff)
899  jm1 = indijk(i ,j-1,k ,inoff,ijnoff)
900  jm2 = indijk(i ,j-2,k ,inoff,ijnoff)
901  kp1 = indijk(i ,j ,k+1,inoff,ijnoff)
902  kp2 = indijk(i ,j ,k+2,inoff,ijnoff)
903  km1 = indijk(i ,j ,k-1,inoff,ijnoff)
904  km2 = indijk(i ,j ,k-2,inoff,ijnoff)
905 
906  sim = sqrt( (xyzorig(xcoord,im1)+xyzold(xcoord,im1)-xyzorig(xcoord,ijk)- &
907  xyzold(xcoord,ijk))**2 + &
908  (xyzorig(ycoord,im1)+xyzold(ycoord,im1)-xyzorig(ycoord,ijk)- &
909  xyzold(ycoord,ijk))**2 + &
910  (xyzorig(zcoord,im1)+xyzold(zcoord,im1)-xyzorig(zcoord,ijk)- &
911  xyzold(zcoord,ijk))**2 )**p
912  sip = sqrt( (xyzorig(xcoord,ip1)+xyzold(xcoord,ip1)-xyzorig(xcoord,ijk)- &
913  xyzold(xcoord,ijk))**2 + &
914  (xyzorig(ycoord,ip1)+xyzold(ycoord,ip1)-xyzorig(ycoord,ijk)- &
915  xyzold(ycoord,ijk))**2 + &
916  (xyzorig(zcoord,ip1)+xyzold(zcoord,ip1)-xyzorig(zcoord,ijk)- &
917  xyzold(zcoord,ijk))**2 )**p
918 
919  sjm = sqrt( (xyzorig(xcoord,jm1)+xyzold(xcoord,jm1)-xyzorig(xcoord,ijk)- &
920  xyzold(xcoord,ijk))**2 + &
921  (xyzorig(ycoord,jm1)+xyzold(ycoord,jm1)-xyzorig(ycoord,ijk)- &
922  xyzold(ycoord,ijk))**2 + &
923  (xyzorig(zcoord,jm1)+xyzold(zcoord,jm1)-xyzorig(zcoord,ijk)- &
924  xyzold(zcoord,ijk))**2 )**p
925  sjp = sqrt( (xyzorig(xcoord,jp1)+xyzold(xcoord,jp1)-xyzorig(xcoord,ijk)- &
926  xyzold(xcoord,ijk))**2 + &
927  (xyzorig(ycoord,jp1)+xyzold(ycoord,jp1)-xyzorig(ycoord,ijk)- &
928  xyzold(ycoord,ijk))**2 + &
929  (xyzorig(zcoord,jp1)+xyzold(zcoord,jp1)-xyzorig(zcoord,ijk)- &
930  xyzold(zcoord,ijk))**2 )**p
931 
932  skm = sqrt( (xyzorig(xcoord,km1)+xyzold(xcoord,km1)-xyzorig(xcoord,ijk)- &
933  xyzold(xcoord,ijk))**2 + &
934  (xyzorig(ycoord,km1)+xyzold(ycoord,km1)-xyzorig(ycoord,ijk)- &
935  xyzold(ycoord,ijk))**2 + &
936  (xyzorig(zcoord,km1)+xyzold(zcoord,km1)-xyzorig(zcoord,ijk)- &
937  xyzold(zcoord,ijk))**2 )**p
938  skp = sqrt( (xyzorig(xcoord,kp1)+xyzold(xcoord,kp1)-xyzorig(xcoord,ijk)- &
939  xyzold(xcoord,ijk))**2 + &
940  (xyzorig(ycoord,kp1)+xyzold(ycoord,kp1)-xyzorig(ycoord,ijk)- &
941  xyzold(ycoord,ijk))**2 + &
942  (xyzorig(zcoord,kp1)+xyzold(zcoord,kp1)-xyzorig(zcoord,ijk)- &
943  xyzold(zcoord,ijk))**2 )**p
944 
945  sim = 1._rfreal/(sim+eps)
946  sip = 1._rfreal/(sip+eps)
947  sjm = 1._rfreal/(sjm+eps)
948  sjp = 1._rfreal/(sjp+eps)
949  skm = 1._rfreal/(skm+eps)
950  skp = 1._rfreal/(skp+eps)
951 
952  rxi = sim*xyzold(xcoord,im1) + sip*xyzold(xcoord,ip1)
953  ryi = sim*xyzold(ycoord,im1) + sip*xyzold(ycoord,ip1)
954  rzi = sim*xyzold(zcoord,im1) + sip*xyzold(zcoord,ip1)
955 
956  rxj = sjm*xyzold(xcoord,jm1) + sjp*xyzold(xcoord,jp1)
957  ryj = sjm*xyzold(ycoord,jm1) + sjp*xyzold(ycoord,jp1)
958  rzj = sjm*xyzold(zcoord,jm1) + sjp*xyzold(zcoord,jp1)
959 
960  rxk = skm*xyzold(xcoord,km1) + skp*xyzold(xcoord,kp1)
961  ryk = skm*xyzold(ycoord,km1) + skp*xyzold(ycoord,kp1)
962  rzk = skm*xyzold(zcoord,km1) + skp*xyzold(zcoord,kp1)
963 
964  d = 1._rfreal/(sim+sip+sjm+sjp+skm+skp)
965 
966 ! xyz(XCOORD,ijk) = (rxi+rxj+rxk)*d
967 ! xyz(YCOORD,ijk) = (ryi+ryj+ryk)*d
968 ! xyz(ZCOORD,ijk) = (rzi+rzj+rzk)*d
969 
970  xyz(xcoord,ijk) = (1._rfreal-REAL( idgen(ijk) ))*(rxi+rxj+rxk)*d + &
971  REAL( idgen(ijk) )*(xyz(xcoord,ijk)-xyzorig(xcoord,ijk))
972  xyz(ycoord,ijk) = (1._rfreal-REAL( idgen(ijk) ))*(ryi+ryj+ryk)*d + &
973  REAL( idgen(ijk) )*(xyz(ycoord,ijk)-xyzorig(ycoord,ijk))
974  xyz(zcoord,ijk) = (1._rfreal-REAL( idgen(ijk) ))*(rzi+rzj+rzk)*d + &
975  REAL( idgen(ijk) )*(xyz(zcoord,ijk)-xyzorig(zcoord,ijk))
976 
977  ENDDO ! i
978  ENDDO ! j
979  ENDDO ! k
980 
981  ELSEIF (method==2) THEN
982 
983 ! - Old method by Jiri Blazek
984 
985  DO k=kpnbeg,kpnend
986  DO j=jpnbeg,jpnend
987  DO i=ipnbeg,ipnend
988  ijk = indijk(i ,j ,k ,inoff,ijnoff)
989  ip1 = indijk(i+1,j ,k ,inoff,ijnoff)
990  ip2 = indijk(i+2,j ,k ,inoff,ijnoff)
991  im1 = indijk(i-1,j ,k ,inoff,ijnoff)
992  im2 = indijk(i-2,j ,k ,inoff,ijnoff)
993  jp1 = indijk(i ,j+1,k ,inoff,ijnoff)
994  jp2 = indijk(i ,j+2,k ,inoff,ijnoff)
995  jm1 = indijk(i ,j-1,k ,inoff,ijnoff)
996  jm2 = indijk(i ,j-2,k ,inoff,ijnoff)
997  kp1 = indijk(i ,j ,k+1,inoff,ijnoff)
998  kp2 = indijk(i ,j ,k+2,inoff,ijnoff)
999  km1 = indijk(i ,j ,k-1,inoff,ijnoff)
1000  km2 = indijk(i ,j ,k-2,inoff,ijnoff)
1001 
1002  rxi = xyzold(xcoord,im1) + xyzold(xcoord,ip1)
1003  ryi = xyzold(ycoord,im1) + xyzold(ycoord,ip1)
1004  rzi = xyzold(zcoord,im1) + xyzold(zcoord,ip1)
1005 
1006  rxj = xyzold(xcoord,jm1) + xyzold(xcoord,jp1)
1007  ryj = xyzold(ycoord,jm1) + xyzold(ycoord,jp1)
1008  rzj = xyzold(zcoord,jm1) + xyzold(zcoord,jp1)
1009 
1010  rxk = xyzold(xcoord,km1) + xyzold(xcoord,kp1)
1011  ryk = xyzold(ycoord,km1) + xyzold(ycoord,kp1)
1012  rzk = xyzold(zcoord,km1) + xyzold(zcoord,kp1)
1013 
1014  si = (xyzold(xcoord,ip1)+xyzorig(xcoord,ip1)- &
1015  xyzold(xcoord,ijk)-xyzorig(xcoord,ijk))**2 + &
1016  (xyzold(ycoord,ip1)+xyzorig(ycoord,ip1)- &
1017  xyzold(ycoord,ijk)-xyzorig(ycoord,ijk))**2 + &
1018  (xyzold(zcoord,ip1)+xyzorig(zcoord,ip1)- &
1019  xyzold(zcoord,ijk)-xyzorig(zcoord,ijk))**2 + &
1020  (xyzold(xcoord,im1)+xyzorig(xcoord,im1)- &
1021  xyzold(xcoord,ijk)-xyzorig(xcoord,ijk))**2 + &
1022  (xyzold(ycoord,im1)+xyzorig(ycoord,im1)- &
1023  xyzold(ycoord,ijk)-xyzorig(ycoord,ijk))**2 + &
1024  (xyzold(zcoord,im1)+xyzorig(zcoord,im1)- &
1025  xyzold(zcoord,ijk)-xyzorig(zcoord,ijk))**2 + &
1026  (xyzold(xcoord,ip2)+xyzorig(xcoord,ip2)- &
1027  xyzold(xcoord,ip1)-xyzorig(xcoord,ip1))**2 + &
1028  (xyzold(ycoord,ip2)+xyzorig(ycoord,ip2)- &
1029  xyzold(ycoord,ip1)-xyzorig(ycoord,ip1))**2 + &
1030  (xyzold(zcoord,ip2)+xyzorig(zcoord,ip2)- &
1031  xyzold(zcoord,ip1)-xyzorig(zcoord,ip1))**2 + &
1032  (xyzold(xcoord,im2)+xyzorig(xcoord,im2)- &
1033  xyzold(xcoord,im1)-xyzorig(xcoord,im1))**2 + &
1034  (xyzold(ycoord,im2)+xyzorig(ycoord,im2)- &
1035  xyzold(ycoord,im1)-xyzorig(ycoord,im1))**2 + &
1036  (xyzold(zcoord,im2)+xyzorig(zcoord,im2)- &
1037  xyzold(zcoord,im1)-xyzorig(zcoord,im1))**2
1038 
1039  sj = (xyzold(xcoord,jp1)+xyzorig(xcoord,jp1)- &
1040  xyzold(xcoord,ijk)-xyzorig(xcoord,ijk))**2 + &
1041  (xyzold(ycoord,jp1)+xyzorig(ycoord,jp1)- &
1042  xyzold(ycoord,ijk)-xyzorig(ycoord,ijk))**2 + &
1043  (xyzold(zcoord,jp1)+xyzorig(zcoord,jp1)- &
1044  xyzold(zcoord,ijk)-xyzorig(zcoord,ijk))**2 + &
1045  (xyzold(xcoord,jm1)+xyzorig(xcoord,jm1)- &
1046  xyzold(xcoord,ijk)-xyzorig(xcoord,ijk))**2 + &
1047  (xyzold(ycoord,jm1)+xyzorig(ycoord,jm1)- &
1048  xyzold(ycoord,ijk)-xyzorig(ycoord,ijk))**2 + &
1049  (xyzold(zcoord,jm1)+xyzorig(zcoord,jm1)- &
1050  xyzold(zcoord,ijk)-xyzorig(zcoord,ijk))**2 + &
1051  (xyzold(xcoord,jp2)+xyzorig(xcoord,jp2)- &
1052  xyzold(xcoord,jp1)-xyzorig(xcoord,jp1))**2 + &
1053  (xyzold(ycoord,jp2)+xyzorig(ycoord,jp2)- &
1054  xyzold(ycoord,jp1)-xyzorig(ycoord,jp1))**2 + &
1055  (xyzold(zcoord,jp2)+xyzorig(zcoord,jp2)- &
1056  xyzold(zcoord,jp1)-xyzorig(zcoord,jp1))**2 + &
1057  (xyzold(xcoord,jm2)+xyzorig(xcoord,jm2)- &
1058  xyzold(xcoord,jm1)-xyzorig(xcoord,jm1))**2 + &
1059  (xyzold(ycoord,jm2)+xyzorig(ycoord,jm2)- &
1060  xyzold(ycoord,jm1)-xyzorig(ycoord,jm1))**2 + &
1061  (xyzold(zcoord,jm2)+xyzorig(zcoord,jm2)- &
1062  xyzold(zcoord,jm1)-xyzorig(zcoord,jm1))**2
1063 
1064  sk = (xyzold(xcoord,kp1)+xyzorig(xcoord,kp1)- &
1065  xyzold(xcoord,ijk)-xyzorig(xcoord,ijk))**2 + &
1066  (xyzold(ycoord,kp1)+xyzorig(ycoord,kp1)- &
1067  xyzold(ycoord,ijk)-xyzorig(ycoord,ijk))**2 + &
1068  (xyzold(zcoord,kp1)+xyzorig(zcoord,kp1)- &
1069  xyzold(zcoord,ijk)-xyzorig(zcoord,ijk))**2 + &
1070  (xyzold(xcoord,km1)+xyzorig(xcoord,km1)- &
1071  xyzold(xcoord,ijk)-xyzorig(xcoord,ijk))**2 + &
1072  (xyzold(ycoord,km1)+xyzorig(ycoord,km1)- &
1073  xyzold(ycoord,ijk)-xyzorig(ycoord,ijk))**2 + &
1074  (xyzold(zcoord,km1)+xyzorig(zcoord,km1)- &
1075  xyzold(zcoord,ijk)-xyzorig(zcoord,ijk))**2 + &
1076  (xyzold(xcoord,kp2)+xyzorig(xcoord,kp2)- &
1077  xyzold(xcoord,kp1)-xyzorig(xcoord,kp1))**2 + &
1078  (xyzold(ycoord,kp2)+xyzorig(ycoord,kp2)- &
1079  xyzold(ycoord,kp1)-xyzorig(ycoord,kp1))**2 + &
1080  (xyzold(zcoord,kp2)+xyzorig(zcoord,kp2)- &
1081  xyzold(zcoord,kp1)-xyzorig(zcoord,kp1))**2 + &
1082  (xyzold(xcoord,km2)+xyzorig(xcoord,km2)- &
1083  xyzold(xcoord,km1)-xyzorig(xcoord,km1))**2 + &
1084  (xyzold(ycoord,km2)+xyzorig(ycoord,km2)- &
1085  xyzold(ycoord,km1)-xyzorig(ycoord,km1))**2 + &
1086  (xyzold(zcoord,km2)+xyzorig(zcoord,km2)- &
1087  xyzold(zcoord,km1)-xyzorig(zcoord,km1))**2
1088 
1089  wi = 1._rfreal/(si)**p
1090  wj = 1._rfreal/(sj)**p
1091  wk = 1._rfreal/(sk)**p
1092  d = 2._rfreal*(wi+wj+wk)
1093 
1094 ! xyz(XCOORD,ijk) = (wi*rxi+wj*rxj+wk*rxk)/d
1095 ! xyz(YCOORD,ijk) = (wi*ryi+wj*ryj+wk*ryk)/d
1096 ! xyz(ZCOORD,ijk) = (wi*rzi+wj*rzj+wk*rzk)/d
1097 
1098  xyz(xcoord,ijk) = (1._rfreal-REAL( idgen(ijk) ))* &
1099  (wi*rxi+wj*rxj+wk*rxk)/d + &
1100  REAL( idgen(ijk) )*(xyz(xcoord,ijk)-xyzorig(xcoord,ijk))
1101  xyz(ycoord,ijk) = (1._rfreal-REAL( idgen(ijk) ))* &
1102  (wi*ryi+wj*ryj+wk*ryk)/d + &
1103  REAL( idgen(ijk) )*(xyz(ycoord,ijk)-xyzorig(ycoord,ijk))
1104  xyz(zcoord,ijk) = (1._rfreal-REAL( idgen(ijk) ))* &
1105  (wi*rzi+wj*rzj+wk*rzk)/d + &
1106  REAL( idgen(ijk) )*(xyz(zcoord,ijk)-xyzorig(zcoord,ijk))
1107  ENDDO ! i
1108  ENDDO ! j
1109  ENDDO ! k
1110  ENDIF ! method
1111 !if (region%iRegionGlobal==13) write(*,*)'xyzOrth0',ipnbeg,ipnend,jpnbeg,jpnend,kpnbeg,kpnend,xyz(XCOORD,IndIJK(1 ,1 ,23 ,iNOff,ijNOff))
1112 
1113 ! finalize --------------------------------------------------------------------
1114 
1115  CALL deregisterfunction( region%global )
1116 
1117 END SUBROUTINE rflo_laplacegridsolve
1118 
1119 
1120 !******************************************************************************
1121 !
1122 ! Purpose: Enforce orthogonal Laplacian smoothing
1123 !
1124 ! Description: none.
1125 !
1126 ! Input: region = data of current region, old grid coordinates, and
1127 ! previous grid movements
1128 !
1129 ! Output: region%levels%grid%xyz = new grid movements.
1130 !
1131 ! Notes: on entry xyz holds the actual grid, on exit xyz holds the grid motion
1132 !
1133 !******************************************************************************
1134 
1135 SUBROUTINE rflo_laplacegridortho( region )
1136 
1137  USE modinterfaces, ONLY : rflo_getdimensdummy, rflo_getdimensdummynodes, &
1138  rflo_getdimensphysnodes, &
1139  rflo_getcelloffset, rflo_getnodeoffset, &
1140  rflo_calcfacevectors
1141  IMPLICIT NONE
1142 
1143 #include "Indexing.h"
1144 
1145 ! ... parameters
1146  TYPE(t_region) :: region
1147 
1148 ! ... loop variables
1149  INTEGER :: i, j, k
1150 
1151 ! ... local variables
1152  INTEGER :: ilev, ibn, ien, ibc, iec, icoff, ijcoff, inoff, ijnoff
1153  INTEGER :: idcbeg, idcend, jdcbeg, jdcend, kdcbeg, kdcend
1154  INTEGER :: idnbeg, idnend, jdnbeg, jdnend, kdnbeg, kdnend
1155  INTEGER :: ipnbeg, ipnend, jpnbeg, jpnend, kpnbeg, kpnend
1156  INTEGER :: ijk, ijkc, ijkn, ijkm, ijkp, ijka(3), errorflag
1157  INTEGER :: ijkn1, ijkn2, ijkn3, ijkn4, ijkn5, ijkn6, ijkn7, ijkn8
1158  INTEGER :: ijkp1, ijkm1, ijkp2, ijkm2, ijkp3, ijkm3, ijkp4, ijkm4
1159  INTEGER, POINTER :: idgen(:)
1160 
1161  REAL(RFREAL) :: rnfac, ocell, eps, oo6, oo8, smod, fac1, fac2
1162  REAL(RFREAL) :: usf(xcoord:zcoord), fc(xcoord:zcoord)
1163  REAL(RFREAL), POINTER :: xyz(:,:), xyzorig(:,:), xyzold(:,:), xyzorth(:,:)
1164  REAL(RFREAL), POINTER :: xyztemp(:,:), si(:,:), sj(:,:), sk(:,:)
1165 ! REAL(RFREAL), POINTER :: cofg(:,:) ! *
1166  REAL(RFREAL), ALLOCATABLE :: cofg(:,:)
1167  REAL(RFREAL), ALLOCATABLE :: xyzwork(:,:), xyzbuf(:,:)
1168  REAL(RFREAL), ALLOCATABLE :: fci(:,:), fcj(:,:), fck(:,:)
1169 
1170  TYPE(t_global), POINTER :: global
1171 
1172 !******************************************************************************
1173 
1174  global => region%global
1175  CALL registerfunction( global,'RFLO_LaplaceGridOrtho',&
1176  'RFLO_ModLaplaceSmoothing.F90' )
1177 
1178 ! get local parameters --------------------------------------------------------
1179 
1180  rnfac = 1._rfreal/12._rfreal
1181  oo6 = 1._rfreal/6._rfreal
1182  oo8 = 1._rfreal/8._rfreal
1183  ocell = global%moveGridOrthCell
1184  eps = 100._rfreal*epsilon(1._rfreal)
1185  fac1 = region%global%moveGridWeight
1186  fac2 = 1._rfreal-fac1
1187 
1188 ! get dimensions and pointers -------------------------------------------------
1189 
1190  ilev = 1
1191 
1192  CALL rflo_getdimensdummy( region,ilev,idcbeg,idcend, &
1193  jdcbeg,jdcend,kdcbeg,kdcend )
1194  CALL rflo_getdimensdummynodes( region,ilev,idnbeg,idnend, &
1195  jdnbeg,jdnend,kdnbeg,kdnend )
1196  CALL rflo_getdimensphysnodes( region,ilev,ipnbeg,ipnend, &
1197  jpnbeg,jpnend,kpnbeg,kpnend )
1198  CALL rflo_getcelloffset( region,ilev,icoff,ijcoff )
1199  CALL rflo_getnodeoffset( region,ilev,inoff,ijnoff )
1200  ibn = indijk(idnbeg,jdnbeg,kdnbeg,inoff,ijnoff)
1201  ien = indijk(idnend,jdnend,kdnend,inoff,ijnoff)
1202  ibc = indijk(idcbeg,jdcbeg,kdcbeg,icoff,ijcoff)
1203  iec = indijk(idcend,jdcend,kdcend,icoff,ijcoff)
1204 
1205  xyz => region%levels(ilev)%grid%xyz
1206  xyzorig => region%levels(ilev)%gridOld%xyz
1207  xyzold => region%levels(ilev)%grid%xyzOld
1208  xyzorth => region%levels(ilev)%grid%xyzOrth
1209  xyztemp => region%levels(ilev)%grid%xyzTemp
1210  si => region%levels(ilev)%grid%si
1211  sj => region%levels(ilev)%grid%sj
1212  sk => region%levels(ilev)%grid%sk
1213 ! cofg => region%levels(iLev)%grid%cofg ! *
1214  idgen => region%levels(ilev)%grid%ijkDgen
1215 
1216  ALLOCATE( xyzwork(3,ibn:ien),stat=errorflag )
1217  global%error = errorflag
1218  IF (global%error /= 0) CALL errorstop( global,err_allocate,&
1219  __line__ )
1220 
1221  ALLOCATE( xyzbuf(3,ibn:ien),stat=errorflag )
1222  global%error = errorflag
1223  IF (global%error /= 0) CALL errorstop( global,err_allocate,&
1224  __line__ )
1225 
1226  ALLOCATE( fci(3,ibn:ien),stat=errorflag )
1227  global%error = errorflag
1228  IF (global%error /= 0) CALL errorstop( global,err_allocate,&
1229  __line__ )
1230 
1231  ALLOCATE( fcj(3,ibn:ien),stat=errorflag )
1232  global%error = errorflag
1233  IF (global%error /= 0) CALL errorstop( global,err_allocate,&
1234  __line__ )
1235 
1236  ALLOCATE( fck(3,ibn:ien),stat=errorflag )
1237  global%error = errorflag
1238  IF (global%error /= 0) CALL errorstop( global,err_allocate,&
1239  __line__ )
1240 
1241  ALLOCATE( cofg(3,ibc:iec),stat=errorflag )
1242  global%error = errorflag
1243  IF (global%error /= 0) CALL errorstop( global,err_allocate,&
1244  __line__ )
1245 
1246  xyzwork = 0._rfreal
1247 
1248 ! reset motion vectors --------------------------------------------------------
1249 
1250  DO ijk=ibn,ien
1251  xyzold(xcoord,ijk) = xyz(xcoord,ijk) - xyzorig(xcoord,ijk)
1252  xyzold(ycoord,ijk) = xyz(ycoord,ijk) - xyzorig(ycoord,ijk)
1253  xyzold(zcoord,ijk) = xyz(zcoord,ijk) - xyzorig(zcoord,ijk)
1254  ENDDO
1255 
1256 ! xyzTemp = xyz
1257  xyztemp = xyzorig + xyzold
1258  xyzbuf = xyz
1259 
1260 ! get cell center -------------------------------------------------------------
1261 
1262  DO k=kdcbeg,kdcend
1263  DO j=jdcbeg,jdcend
1264  DO i=idcbeg,idcend
1265  ijkc = indijk(i ,j ,k ,icoff,ijcoff)
1266  ijkn1 = indijk(i ,j ,k ,inoff,ijnoff)
1267  ijkn2 = indijk(i ,j ,k+1,inoff,ijnoff)
1268  ijkn3 = indijk(i ,j+1,k+1,inoff,ijnoff)
1269  ijkn4 = indijk(i ,j+1,k ,inoff,ijnoff)
1270  ijkn5 = indijk(i+1,j ,k ,inoff,ijnoff)
1271  ijkn6 = indijk(i+1,j ,k+1,inoff,ijnoff)
1272  ijkn7 = indijk(i+1,j+1,k+1,inoff,ijnoff)
1273  ijkn8 = indijk(i+1,j+1,k ,inoff,ijnoff)
1274  cofg(:,ijkc) = xyztemp(:,ijkn1)+xyztemp(:,ijkn2)+xyztemp(:,ijkn3)+ &
1275  xyztemp(:,ijkn4)+xyztemp(:,ijkn5)+xyztemp(:,ijkn6)+ &
1276  xyztemp(:,ijkn7)+xyztemp(:,ijkn8)
1277  cofg(:,ijkc) = oo8*cofg(:,ijkc)
1278  ENDDO
1279  ENDDO
1280  ENDDO
1281 
1282 ! perturb grid ----------------------------------------------------------------
1283 
1284  DO k=kdnbeg+1,kdnend-1
1285  DO j=jdnbeg+1,jdnend-1
1286  DO i=idnbeg+1,idnend-1
1287  ijkn = indijk(i ,j ,k ,inoff,ijnoff)
1288  ijkn1 = indijk(i+1,j ,k ,inoff,ijnoff)
1289  ijkn2 = indijk(i ,j+1,k ,inoff,ijnoff)
1290  ijkn3 = indijk(i ,j ,k+1,inoff,ijnoff)
1291  ijkn4 = indijk(i-1,j ,k ,inoff,ijnoff)
1292  ijkn5 = indijk(i ,j-1,k ,inoff,ijnoff)
1293  ijkn6 = indijk(i ,j ,k-1,inoff,ijnoff)
1294  xyzorth(:,ijkn) = xyztemp(:,ijkn1)+xyztemp(:,ijkn2)+xyztemp(:,ijkn3)+ &
1295  xyztemp(:,ijkn4)+xyztemp(:,ijkn5)+xyztemp(:,ijkn6)
1296  xyz(:,ijkn) = oo6*xyzorth(:,ijkn)
1297  ENDDO
1298  ENDDO
1299  ENDDO
1300 
1301 ! get face vectors ------------------------------------------------------------
1302 
1303  CALL rflo_calcfacevectors( region )
1304 
1305 ! weighted face-centers -------------------------------------------------------
1306 
1307  DO k=kdnbeg+1,kdnend-1
1308  DO j=jdnbeg+1,jdnend-1
1309  DO i=idnbeg+1,idnbeg+1+(idnend-1-idnbeg-1)/2
1310  ijkn = indijk(i ,j ,k ,inoff,ijnoff)
1311  ijkp = indijk(i ,j ,k ,icoff,ijcoff)
1312  ijkm = indijk(i-1,j ,k ,icoff,ijcoff)
1313  fci(:,ijkn) = fac1*cofg(:,ijkm) + fac2*cofg(:,ijkp)
1314  ENDDO
1315  DO i=idnbeg+1+(idnend-1-idnbeg-1)/2+1,idnend-1
1316  ijkn = indijk(i ,j ,k ,inoff,ijnoff)
1317  ijkp = indijk(i ,j ,k ,icoff,ijcoff)
1318  ijkm = indijk(i-1,j ,k ,icoff,ijcoff)
1319  fci(:,ijkn) = fac2*cofg(:,ijkm) + fac1*cofg(:,ijkp)
1320  ENDDO
1321  ENDDO
1322  ENDDO
1323  DO k=kdnbeg+1,kdnend-1
1324  DO j=jdnbeg+1,jdnbeg+1+(jdnend-1-jdnbeg-1)/2
1325  DO i=idnbeg+1,idnend-1
1326  ijkn = indijk(i ,j ,k ,inoff,ijnoff)
1327  ijkp = indijk(i ,j ,k ,icoff,ijcoff)
1328  ijkm = indijk(i ,j-1,k ,icoff,ijcoff)
1329  fcj(:,ijkn) = fac1*cofg(:,ijkm) + fac2*cofg(:,ijkp)
1330  ENDDO
1331  ENDDO
1332  DO j=jdnbeg+1+(jdnend-1-jdnbeg-1)/2+1,jdnend-1
1333  DO i=idnbeg+1,idnend-1
1334  ijkn = indijk(i ,j ,k ,inoff,ijnoff)
1335  ijkp = indijk(i ,j ,k ,icoff,ijcoff)
1336  ijkm = indijk(i ,j-1,k ,icoff,ijcoff)
1337  fcj(:,ijkn) = fac2*cofg(:,ijkm) + fac1*cofg(:,ijkp)
1338  ENDDO
1339  ENDDO
1340  ENDDO
1341  DO k=kdnbeg+1,kdnbeg+1+(kdnend-1-kdnbeg-1)/2
1342  DO j=jdnbeg+1,jdnend-1
1343  DO i=idnbeg+1,idnend-1
1344  ijkn = indijk(i ,j ,k ,inoff,ijnoff)
1345  ijkp = indijk(i ,j ,k ,icoff,ijcoff)
1346  ijkm = indijk(i ,j ,k-1,icoff,ijcoff)
1347  fck(:,ijkn) = fac1*cofg(:,ijkm) + fac2*cofg(:,ijkp)
1348  ENDDO
1349  ENDDO
1350  ENDDO
1351  DO k=kdnbeg+1+(kdnend-1-kdnbeg-1)/2+1,kdnend-1
1352  DO j=jdnbeg+1,jdnend-1
1353  DO i=idnbeg+1,idnend-1
1354  ijkn = indijk(i ,j ,k ,inoff,ijnoff)
1355  ijkp = indijk(i ,j ,k ,icoff,ijcoff)
1356  ijkm = indijk(i ,j ,k-1,icoff,ijcoff)
1357  fck(:,ijkn) = fac2*cofg(:,ijkm) + fac1*cofg(:,ijkp)
1358  ENDDO
1359  ENDDO
1360  ENDDO
1361 
1362 ! orthogonal Laplacian --------------------------------------------------------
1363 
1364  DO k=kpnbeg,kpnend
1365  DO j=jpnbeg,jpnend
1366  DO i=ipnbeg,ipnend
1367 
1368 ! ----- I-Face
1369 
1370  ijkn1 = indijk(i ,j ,k ,inoff,ijnoff)
1371  ijkn2 = indijk(i ,j-1,k ,inoff,ijnoff)
1372  ijkn3 = indijk(i ,j-1,k-1,inoff,ijnoff)
1373  ijkn4 = indijk(i ,j ,k-1,inoff,ijnoff)
1374  ijkp1 = indijk(i ,j ,k ,icoff,ijcoff)
1375  ijkm1 = indijk(i-1,j ,k ,icoff,ijcoff)
1376  ijkp2 = indijk(i ,j-1,k ,icoff,ijcoff)
1377  ijkm2 = indijk(i-1,j-1,k ,icoff,ijcoff)
1378  ijkp3 = indijk(i ,j-1,k-1,icoff,ijcoff)
1379  ijkm3 = indijk(i-1,j-1,k-1,icoff,ijcoff)
1380  ijkp4 = indijk(i ,j ,k-1,icoff,ijcoff)
1381  ijkm4 = indijk(i-1,j ,k-1,icoff,ijcoff)
1382 !IF (i==2.AND.j==2.AND.k==2) write(*,*)'orthoI'
1383  CALL accumulatexyzwork( si, &
1384  fci(:,ijkn1),fci(:,ijkn2),fci(:,ijkn3),fci(:,ijkn4) )
1385 
1386 ! ----- J-Face
1387 
1388  ijkn1 = indijk(i ,j ,k ,inoff,ijnoff)
1389  ijkn2 = indijk(i ,j ,k-1,inoff,ijnoff)
1390  ijkn3 = indijk(i-1,j ,k-1,inoff,ijnoff)
1391  ijkn4 = indijk(i-1,j ,k ,inoff,ijnoff)
1392  ijkp1 = indijk(i ,j ,k ,icoff,ijcoff)
1393  ijkm1 = indijk(i ,j-1,k ,icoff,ijcoff)
1394  ijkp2 = indijk(i ,j ,k-1,icoff,ijcoff)
1395  ijkm2 = indijk(i ,j-1,k-1,icoff,ijcoff)
1396  ijkp3 = indijk(i-1,j ,k-1,icoff,ijcoff)
1397  ijkm3 = indijk(i-1,j-1,k-1,icoff,ijcoff)
1398  ijkp4 = indijk(i-1,j ,k ,icoff,ijcoff)
1399  ijkm4 = indijk(i-1,j-1,k ,icoff,ijcoff)
1400 !IF (i==2.AND.j==2.AND.k==2) write(*,*)'orthoJ'
1401  CALL accumulatexyzwork( sj, &
1402  fcj(:,ijkn1),fcj(:,ijkn2),fcj(:,ijkn3),fcj(:,ijkn4) )
1403 
1404 ! ----- K-Face
1405 
1406  ijkn1 = indijk(i ,j ,k ,inoff,ijnoff)
1407  ijkn2 = indijk(i-1,j ,k ,inoff,ijnoff)
1408  ijkn3 = indijk(i-1,j-1,k ,inoff,ijnoff)
1409  ijkn4 = indijk(i ,j-1,k ,inoff,ijnoff)
1410  ijkp1 = indijk(i ,j ,k ,icoff,ijcoff)
1411  ijkm1 = indijk(i ,j ,k-1,icoff,ijcoff)
1412  ijkp2 = indijk(i-1,j ,k ,icoff,ijcoff)
1413  ijkm2 = indijk(i-1,j ,k-1,icoff,ijcoff)
1414  ijkp3 = indijk(i-1,j-1,k ,icoff,ijcoff)
1415  ijkm3 = indijk(i-1,j-1,k-1,icoff,ijcoff)
1416  ijkp4 = indijk(i ,j-1,k ,icoff,ijcoff)
1417  ijkm4 = indijk(i ,j-1,k-1,icoff,ijcoff)
1418 !IF (i==2.AND.j==2.AND.k==2) write(*,*)'orthoK'
1419  CALL accumulatexyzwork( sk, &
1420  fck(:,ijkn1),fck(:,ijkn2),fck(:,ijkn3),fck(:,ijkn4) )
1421 
1422  ENDDO ! i
1423  ENDDO ! j
1424  ENDDO ! k
1425 
1426 ! average accumulated grid coordinates
1427 
1428  xyz = xyzbuf-xyzorig
1429  ijka(:) = 1
1430 
1431  DO k=kpnbeg+ijka(3),kpnend-ijka(3)
1432  DO j=jpnbeg+ijka(2),jpnend-ijka(2)
1433  DO i=ipnbeg+ijka(1),ipnend-ijka(1)
1434  ijkn = indijk(i ,j ,k ,inoff,ijnoff)
1435  xyzorth(:,ijkn) = ocell*rnfac*xyzwork(:,ijkn)
1436 ! +(1._RFREAL-ocell)*xyz(:,ijkN)
1437  xyz(:,ijkn) = xyzorth(:,ijkn)
1438 ! xyz(:,ijkN) = (1._RFREAL-REAL( idgen(ijkN) ))*xyzOrth(:,ijkN) + &
1439 ! REAL( idgen(ijkN) )*xyz(:,ijkN)
1440  ENDDO ! i
1441  ENDDO ! j
1442  ENDDO ! k
1443 
1444 ! deallocate temporary memory -------------------------------------------------
1445 
1446  DEALLOCATE( xyzwork,stat=errorflag )
1447  global%error = errorflag
1448  IF (global%error /= 0) CALL errorstop( global,err_deallocate,&
1449  __line__ )
1450 
1451  DEALLOCATE( xyzbuf,stat=errorflag )
1452  global%error = errorflag
1453  IF (global%error /= 0) CALL errorstop( global,err_deallocate,&
1454  __line__ )
1455 
1456  DEALLOCATE( fci,stat=errorflag )
1457  global%error = errorflag
1458  IF (global%error /= 0) CALL errorstop( global,err_deallocate,&
1459  __line__ )
1460 
1461  DEALLOCATE( fcj,stat=errorflag )
1462  global%error = errorflag
1463  IF (global%error /= 0) CALL errorstop( global,err_deallocate,&
1464  __line__ )
1465 
1466  DEALLOCATE( fck,stat=errorflag )
1467  global%error = errorflag
1468  IF (global%error /= 0) CALL errorstop( global,err_deallocate,&
1469  __line__ )
1470 
1471  DEALLOCATE( cofg,stat=errorflag )
1472  global%error = errorflag
1473  IF (global%error /= 0) CALL errorstop( global,err_deallocate,&
1474  __line__ )
1475 
1476 ! finalize --------------------------------------------------------------------
1477 
1478  CALL deregisterfunction( global )
1479 
1480 CONTAINS
1481 
1482 ! ==============================================================================
1483 ! subroutine for accumulation of projected cell corners coordinates
1484 ! ==============================================================================
1485 
1486  SUBROUTINE accumulatexyzwork( sf,fc1,fc2,fc3,fc4 )
1487 
1488  IMPLICIT NONE
1489  REAL(RFREAL), POINTER :: sf(:,:)
1490  REAL(RFREAL) :: fc1(xcoord:zcoord), fc2(xcoord:zcoord)
1491  REAL(RFREAL) :: fc3(xcoord:zcoord), fc4(xcoord:zcoord)
1492 
1493  smod = sqrt( sf(xcoord,ijkn1)*sf(xcoord,ijkn1) + &
1494  sf(ycoord,ijkn1)*sf(ycoord,ijkn1) + &
1495  sf(zcoord,ijkn1)*sf(zcoord,ijkn1) )
1496  usf(:) = sf(:,ijkn1)/(smod+eps)
1497 ! fc1(:) = 0.5_RFREAL*(cofg(:,ijkp1)+cofg(:,ijkm1))
1498  xyzwork(:,ijkn1) = xyzwork(:,ijkn1) + &
1499  dot_product( fc1(:)-xyztemp(:,ijkn1), usf(:) )*usf(:)
1500 !IF (i==2.AND.j==2.AND.k==2) write(*,*)'ortho1',xyzWork(:,ijkN1),usf(:)
1501 
1502  smod = sqrt( sf(xcoord,ijkn2)*sf(xcoord,ijkn2) + &
1503  sf(ycoord,ijkn2)*sf(ycoord,ijkn2) + &
1504  sf(zcoord,ijkn2)*sf(zcoord,ijkn2) )
1505  usf(:) = sf(:,ijkn2)/(smod+eps)
1506 ! fc2(:) = 0.5_RFREAL*(cofg(:,ijkp2)+cofg(:,ijkm2))
1507  xyzwork(:,ijkn1) = xyzwork(:,ijkn1) + &
1508  dot_product( fc2(:)-xyztemp(:,ijkn1), usf(:) )*usf(:)
1509 !IF (i==2.AND.j==2.AND.k==2) write(*,*)'ortho2',xyzWork(:,ijkN1),usf(:)
1510 
1511  smod = sqrt( sf(xcoord,ijkn3)*sf(xcoord,ijkn3) + &
1512  sf(ycoord,ijkn3)*sf(ycoord,ijkn3) + &
1513  sf(zcoord,ijkn3)*sf(zcoord,ijkn3) )
1514  usf(:) = sf(:,ijkn3)/(smod+eps)
1515 ! fc3(:) = 0.5_RFREAL*(cofg(:,ijkp3)+cofg(:,ijkm3))
1516  xyzwork(:,ijkn1) = xyzwork(:,ijkn1) + &
1517  dot_product( fc3(:)-xyztemp(:,ijkn1), usf(:) )*usf(:)
1518 !IF (i==2.AND.j==2.AND.k==2) write(*,*)'ortho3',xyzWork(:,ijkN1),usf(:)
1519 
1520  smod = sqrt( sf(xcoord,ijkn4)*sf(xcoord,ijkn4) + &
1521  sf(ycoord,ijkn4)*sf(ycoord,ijkn4) + &
1522  sf(zcoord,ijkn4)*sf(zcoord,ijkn4) )
1523  usf(:) = sf(:,ijkn4)/(smod+eps)
1524 ! fc4(:) = 0.5_RFREAL*(cofg(:,ijkp4)+cofg(:,ijkm4))
1525  xyzwork(:,ijkn1) = xyzwork(:,ijkn1) + &
1526  dot_product( fc4(:)-xyztemp(:,ijkn1), usf(:) )*usf(:)
1527 !IF (i==2.AND.j==2.AND.k==2) write(*,*)'ortho4',xyzWork(:,ijkN1),usf(:)
1528 
1529  END SUBROUTINE accumulatexyzwork
1530 
1531 END SUBROUTINE rflo_laplacegridortho
1532 
1533 !******************************************************************************
1534 !
1535 ! Purpose: Enforce orthogonality to Laplacian smoothing
1536 !
1537 ! Description: none.
1538 !
1539 ! Input: region = data of current region, old grid coordinates, and
1540 ! previous grid movements
1541 !
1542 ! Output: region%levels%grid%xyz = new grid movements.
1543 !
1544 ! Notes: on entry xyz holds the actual grid, on exit xyz holds the grid motion
1545 !
1546 !******************************************************************************
1547 
1548 SUBROUTINE rflo_laplacegridortho1( region )
1549 
1550  USE modinterfaces, ONLY : rflo_getdimensdummynodes, &
1551  rflo_getdimensphysnodes, &
1552  rflo_getcelloffset, rflo_getnodeoffset
1553  IMPLICIT NONE
1554 
1555 #include "Indexing.h"
1556 
1557 ! ... parameters
1558  TYPE(t_region) :: region
1559 
1560 ! ... loop variables
1561  INTEGER :: i, j, k
1562 
1563 ! ... local variables
1564  INTEGER :: ilev, ibn, ien, icoff, ijcoff, inoff, ijnoff
1565  INTEGER :: idnbeg, idnend, jdnbeg, jdnend, kdnbeg, kdnend
1566  INTEGER :: ipnbeg, ipnend, jpnbeg, jpnend, kpnbeg, kpnend
1567  INTEGER :: ijk, ijkc1, ijkc2, ijkc3, ijkc4, ijkn1, ijkn2, ijkn3, ijkn4
1568  INTEGER :: ijkc1m, ijkc2m, ijkc3m, ijkc4m, ia, ja, ka, errorflag
1569 
1570  REAL(RFREAL) :: rnfac, ocell, eps
1571  REAL(RFREAL) :: cofgp(xcoord:zcoord), cofgm(xcoord:zcoord), cfc(xcoord:zcoord)
1572  REAL(RFREAL) :: xyzc(xcoord:zcoord), xyzp(xcoord:zcoord)
1573  REAL(RFREAL), POINTER :: xyz(:,:), xyzorig(:,:), xyzold(:,:), xyzorth(:,:)
1574  REAL(RFREAL), POINTER :: cofg(:,:), cfci(:,:), cfcj(:,:), cfck(:,:)
1575  REAL(RFREAL), ALLOCATABLE :: xyzwork(:,:)
1576 
1577  TYPE(t_global), POINTER :: global
1578 
1579 !******************************************************************************
1580 
1581  global => region%global
1582  CALL registerfunction( global,'RFLO_LaplaceGridOrtho1',&
1583  'RFLO_ModLaplaceSmoothing.F90' )
1584 
1585 ! get local parameters --------------------------------------------------------
1586 
1587  rnfac = 1._rfreal/12._rfreal
1588  ocell = global%moveGridOrthCell
1589  eps = 100._rfreal*epsilon(1._rfreal)
1590 
1591 ! get dimensions and pointers -------------------------------------------------
1592 
1593  ilev = 1
1594  ia = 1
1595  ja = 1
1596  ka = 1
1597 
1598  CALL rflo_getdimensdummynodes( region,ilev,idnbeg,idnend, &
1599  jdnbeg,jdnend,kdnbeg,kdnend )
1600  CALL rflo_getdimensphysnodes( region,ilev,ipnbeg,ipnend, &
1601  jpnbeg,jpnend,kpnbeg,kpnend )
1602  CALL rflo_getcelloffset( region,ilev,icoff,ijcoff )
1603  CALL rflo_getnodeoffset( region,ilev,inoff,ijnoff )
1604  ibn = indijk(idnbeg,jdnbeg,kdnbeg,inoff,ijnoff)
1605  ien = indijk(idnend,jdnend,kdnend,inoff,ijnoff)
1606 
1607  xyz => region%levels(ilev)%grid%xyz
1608  xyzorig => region%levels(ilev)%gridOld%xyz
1609  xyzold => region%levels(ilev)%grid%xyzOld
1610  xyzorth => region%levels(ilev)%grid%xyzOrth
1611  cofg => region%levels(ilev)%grid%cofg
1612  cfci => region%levels(ilev)%grid%cfcI
1613  cfcj => region%levels(ilev)%grid%cfcJ
1614  cfck => region%levels(ilev)%grid%cfcK
1615 
1616  ALLOCATE( xyzwork(3,ibn:ien),stat=errorflag )
1617  global%error = errorflag
1618  IF (global%error /= 0) CALL errorstop( global,err_allocate,&
1619  __line__ )
1620 
1621  xyzwork = 0._rfreal
1622 
1623 ! reset motion vectors --------------------------------------------------------
1624 
1625 ! xyzOrth = xyz
1626  xyzorth = xyz + xyzorig
1627 
1628 ! DO ijk=ibn,ien
1629 ! xyzOld(XCOORD,ijk) = xyz(XCOORD,ijk) - xyzOrig(XCOORD,ijk)
1630 ! xyzOld(YCOORD,ijk) = xyz(YCOORD,ijk) - xyzOrig(YCOORD,ijk)
1631 ! xyzOld(ZCOORD,ijk) = xyz(ZCOORD,ijk) - xyzOrig(ZCOORD,ijk)
1632 ! xyz(XCOORD,ijk) = xyz(XCOORD,ijk) - xyzOrig(XCOORD,ijk)
1633 ! xyz(YCOORD,ijk) = xyz(YCOORD,ijk) - xyzOrig(YCOORD,ijk)
1634 ! xyz(ZCOORD,ijk) = xyz(ZCOORD,ijk) - xyzOrig(ZCOORD,ijk)
1635 ! ENDDO
1636 
1637 !if (region%iRegionGlobal==1) write(*,*)'xyzOrth0',xyzOrth(:,IndIJK(2,4,2 ,iNOff,ijNOff))-xyzOrig(:,IndIJK(2,4,2 ,iNOff,ijNOff))
1638 
1639 ! reset motion vectors --------------------------------------------------------
1640 
1641  DO k=kpnbeg+ka,kpnend-ka
1642  DO j=jpnbeg+ja,jpnend-ja
1643  DO i=ipnbeg+ia,ipnend-ia
1644 
1645 ! ----- I-Face
1646 
1647  ijkc1 = indijk(i ,j ,k ,icoff,ijcoff)
1648  ijkc1m = indijk(i-1,j ,k ,icoff,ijcoff)
1649  ijkc2 = indijk(i ,j ,k-1,icoff,ijcoff)
1650  ijkc2m = indijk(i-1,j ,k-1,icoff,ijcoff)
1651  ijkc3 = indijk(i ,j-1,k-1,icoff,ijcoff)
1652  ijkc3m = indijk(i-1,j-1,k-1,icoff,ijcoff)
1653  ijkc4 = indijk(i ,j-1,k ,icoff,ijcoff)
1654  ijkc4m = indijk(i-1,j-1,k ,icoff,ijcoff)
1655  ijkn1 = indijk(i ,j ,k ,inoff,ijnoff)
1656  ijkn2 = indijk(i ,j ,k-1,inoff,ijnoff)
1657  ijkn3 = indijk(i ,j-1,k-1,inoff,ijnoff)
1658  ijkn4 = indijk(i ,j-1,k ,inoff,ijnoff)
1659 
1660  CALL accumulatexyzwork( cfci )
1661 !if (i==2 .AND. j==4 .AND. k==2) write(*,*)'xyzOrthI',xyzWork(:,IndIJK(2,4,2 ,iNOff,ijNOff))/4-xyzOrig(:,IndIJK(2,4,2 ,iNOff,ijNOff))
1662 
1663 ! ----- J-Face
1664 
1665  ijkc1 = indijk(i ,j ,k ,icoff,ijcoff)
1666  ijkc1m = indijk(i ,j-1,k ,icoff,ijcoff)
1667  ijkc2 = indijk(i-1,j ,k ,icoff,ijcoff)
1668  ijkc2m = indijk(i-1,j-1,k ,icoff,ijcoff)
1669  ijkc3 = indijk(i-1,j ,k-1,icoff,ijcoff)
1670  ijkc3m = indijk(i-1,j-1,k-1,icoff,ijcoff)
1671  ijkc4 = indijk(i ,j ,k-1,icoff,ijcoff)
1672  ijkc4m = indijk(i ,j-1,k-1,icoff,ijcoff)
1673  ijkn1 = indijk(i ,j ,k ,inoff,ijnoff)
1674  ijkn2 = indijk(i-1,j ,k ,inoff,ijnoff)
1675  ijkn3 = indijk(i-1,j ,k-1,inoff,ijnoff)
1676  ijkn4 = indijk(i ,j ,k-1,inoff,ijnoff)
1677 
1678  CALL accumulatexyzwork( cfcj )
1679 !if (i==2 .AND. j==4 .AND. k==2) write(*,*)'xyzOrthJ',xyzWork(:,IndIJK(2,4,2 ,iNOff,ijNOff))/8-xyzOrig(:,IndIJK(2,4,2 ,iNOff,ijNOff))
1680 
1681 ! ----- K-Face
1682 
1683  ijkc1 = indijk(i ,j ,k ,icoff,ijcoff)
1684  ijkc1m = indijk(i ,j ,k-1,icoff,ijcoff)
1685  ijkc2 = indijk(i ,j-1,k ,icoff,ijcoff)
1686  ijkc2m = indijk(i ,j-1,k-1,icoff,ijcoff)
1687  ijkc3 = indijk(i-1,j-1,k ,icoff,ijcoff)
1688  ijkc3m = indijk(i-1,j-1,k-1,icoff,ijcoff)
1689  ijkc4 = indijk(i-1,j ,k ,icoff,ijcoff)
1690  ijkc4m = indijk(i-1,j ,k-1,icoff,ijcoff)
1691  ijkn1 = indijk(i ,j ,k ,inoff,ijnoff)
1692  ijkn2 = indijk(i ,j-1,k ,inoff,ijnoff)
1693  ijkn3 = indijk(i-1,j-1,k ,inoff,ijnoff)
1694  ijkn4 = indijk(i-1,j ,k ,inoff,ijnoff)
1695 
1696  CALL accumulatexyzwork( cfck )
1697 !if (i==2 .AND. j==4 .AND. k==2) write(*,*)'xyzOrthK',xyzWork(:,IndIJK(2,4,2 ,iNOff,ijNOff))/12-xyzOrig(:,IndIJK(2,4,2 ,iNOff,ijNOff))
1698 
1699 ! ----- average accumulated grid coordinates
1700 
1701 !if (region%iRegionGlobal==13 .AND. i==10 .AND. j==10 .AND. k==10) THEN
1702 ! xyzc(:) = rnfac*xyzWork(:,ijkN1)
1703 ! xyz(XCOORD,ijkN1) = xyzc(XCOORD) - xyzOrig(XCOORD,ijkN1)
1704 ! xyz(YCOORD,ijkN1) = xyzc(YCOORD) - xyzOrig(YCOORD,ijkN1)
1705 ! xyz(ZCOORD,ijkN1) = xyzc(ZCOORD) - xyzOrig(ZCOORD,ijkN1)
1706 !endif
1707  ENDDO ! i
1708  ENDDO ! j
1709  ENDDO ! k
1710 
1711 ! average accumulated grid coordinates
1712 
1713  DO k=kpnbeg+ka,kpnend-ka
1714  DO j=jpnbeg+ja,jpnend-ja
1715  DO i=ipnbeg+ia,ipnend-ia
1716  ijk = indijk(i ,j ,k ,inoff,ijnoff)
1717  xyzorth(:,ijk) = rnfac*xyzwork(:,ijk) - xyzorig(:,ijk)
1718  xyz(xcoord,ijk) = xyzorth(xcoord,ijk)
1719  xyz(ycoord,ijk) = xyzorth(ycoord,ijk)
1720  xyz(zcoord,ijk) = xyzorth(zcoord,ijk)
1721  ENDDO ! i
1722  ENDDO ! j
1723  ENDDO ! k
1724 
1725 !if (region%iRegionGlobal==1) write(*,*)'xyzOrth5',xyzOrth(:,IndIJK(2,4,2 ,iNOff,ijNOff))+xyzOrig(:,IndIJK(2,4,2 ,iNOff,ijNOff))
1726 ! deallocate temporary memory -------------------------------------------------
1727 
1728  DEALLOCATE( xyzwork,stat=errorflag )
1729  global%error = errorflag
1730  IF (global%error /= 0) CALL errorstop( global,err_deallocate,&
1731  __line__ )
1732 
1733 ! finalize --------------------------------------------------------------------
1734 
1735  CALL deregisterfunction( global )
1736 
1737 CONTAINS
1738 
1739 ! ==============================================================================
1740 ! subroutine for accumulation of projected cell corners coordinates
1741 ! ==============================================================================
1742 
1743  SUBROUTINE accumulatexyzwork( cfac )
1744 
1745  IMPLICIT NONE
1746  REAL(RFREAL), POINTER :: cfac(:,:)
1747 
1748  cofgp(:) = cofg(:,ijkc1)
1749  cofgm(:) = cofg(:,ijkc1m)
1750  cfc(:) = cfac(:,ijkn1)
1751 ! cfc(:) = 0.5_RFREAL*(cofgp(:)+cofgm(:))
1752  xyzc(:) = xyzorth(:,ijkn1)
1753  CALL rflo_projectquadcorner( ocell,eps,cofgp,cofgm,cfc,xyzc,xyzp )
1754  xyzwork(:,ijkn1) = xyzwork(:,ijkn1) + xyzp(:)
1755 !write(*,*)'xyzOrth1',xyzp(:),xyzc(:)
1756 !write(*,*)'xyzOrth1',xyzWork(:,ijkN1)
1757  cofgp(:) = cofg(:,ijkc2)
1758  cofgm(:) = cofg(:,ijkc2m)
1759  cfc(:) = cfac(:,ijkn2)
1760 ! cfc(:) = 0.5_RFREAL*(cofgp(:)+cofgm(:))
1761  xyzc(:) = xyzorth(:,ijkn1)
1762  CALL rflo_projectquadcorner( ocell,eps,cofgp,cofgm,cfc,xyzc,xyzp )
1763  xyzwork(:,ijkn1) = xyzwork(:,ijkn1) + xyzp(:)
1764 !write(*,*)'xyzOrth2',xyzp(:),xyzc(:)
1765 !write(*,*)'xyzOrth2',xyzWork(:,ijkN1)
1766  cofgp(:) = cofg(:,ijkc3)
1767  cofgm(:) = cofg(:,ijkc3m)
1768  cfc(:) = cfac(:,ijkn3)
1769 ! cfc(:) = 0.5_RFREAL*(cofgp(:)+cofgm(:))
1770  xyzc(:) = xyzorth(:,ijkn1)
1771  CALL rflo_projectquadcorner( ocell,eps,cofgp,cofgm,cfc,xyzc,xyzp )
1772  xyzwork(:,ijkn1) = xyzwork(:,ijkn1) + xyzp(:)
1773 !write(*,*)'xyzOrth3',xyzp(:),xyzc(:)
1774 !write(*,*)'xyzOrth3',xyzWork(:,ijkN1)
1775  cofgp(:) = cofg(:,ijkc4)
1776  cofgm(:) = cofg(:,ijkc4m)
1777  cfc(:) = cfac(:,ijkn4)
1778 ! cfc(:) = 0.5_RFREAL*(cofgp(:)+cofgm(:))
1779  xyzc(:) = xyzorth(:,ijkn1)
1780  CALL rflo_projectquadcorner( ocell,eps,cofgp,cofgm,cfc,xyzc,xyzp )
1781  xyzwork(:,ijkn1) = xyzwork(:,ijkn1) + xyzp(:)
1782 !write(*,*)'xyzOrth4',xyzp(:),xyzc(:)
1783 !write(*,*)'xyzOrth4',xyzWork(:,ijkN1)
1784  END SUBROUTINE accumulatexyzwork
1785 
1786 END SUBROUTINE rflo_laplacegridortho1
1787 
1788 !******************************************************************************
1789 !
1790 ! Purpose: Enforce orthogonality to Laplacian smoothing
1791 !
1792 ! Description: none.
1793 !
1794 ! Input: region = data of current region, old grid coordinates, and
1795 ! previous grid movements
1796 !
1797 ! Output: region%levels%grid%xyz = new grid movements.
1798 !
1799 ! Notes: on entry xyz holds the actual grid, on exit xyz holds the grid motion
1800 !
1801 !******************************************************************************
1802 
1803 SUBROUTINE rflo_laplacegridortho2( region )
1804 
1805  USE modinterfaces, ONLY : rflo_getdimensdummynodes, &
1806  rflo_getdimensphysnodes, &
1807  rflo_getcelloffset, rflo_getnodeoffset
1808  IMPLICIT NONE
1809 
1810 #include "Indexing.h"
1811 
1812 ! ... parameters
1813  TYPE(t_region) :: region
1814 
1815 ! ... loop variables
1816  INTEGER :: i, j, k
1817 
1818 ! ... local variables
1819  INTEGER :: ilev, ibn, ien, icoff, ijcoff, inoff, ijnoff
1820  INTEGER :: idnbeg, idnend, jdnbeg, jdnend, kdnbeg, kdnend
1821  INTEGER :: ipnbeg, ipnend, jpnbeg, jpnend, kpnbeg, kpnend
1822  INTEGER :: ijk, ijkn, ijkn1, ijkn2, ijkn3, ijkn4, errorflag
1823 
1824  REAL(RFREAL) :: rnfac, ocell, eps, a11, a12, a21, a22, rh1, rh2, rja, c1, c2
1825  REAL(RFREAL) :: delt(xcoord:zcoord), f1(xcoord:zcoord), f2(xcoord:zcoord)
1826  REAL(RFREAL) :: xyzp(xcoord:zcoord)
1827  REAL(RFREAL), POINTER :: xyz(:,:), xyzorig(:,:), xyzold(:,:), xyzorth(:,:)
1828  REAL(RFREAL), POINTER :: xyztemp(:,:)
1829  REAL(RFREAL), ALLOCATABLE :: xyzwork(:,:)
1830 
1831  TYPE(t_global), POINTER :: global
1832 
1833 !******************************************************************************
1834 
1835  global => region%global
1836  CALL registerfunction( global,'RFLO_LaplaceGridOrtho2',&
1837  'RFLO_ModLaplaceSmoothing.F90' )
1838 
1839 ! get local parameters --------------------------------------------------------
1840 
1841  rnfac = 1._rfreal/12._rfreal
1842  ocell = global%moveGridOrthCell
1843  eps = 100._rfreal*epsilon(1._rfreal)
1844 
1845 ! get dimensions and pointers -------------------------------------------------
1846 
1847  ilev = 1
1848 
1849  CALL rflo_getdimensdummynodes( region,ilev,idnbeg,idnend, &
1850  jdnbeg,jdnend,kdnbeg,kdnend )
1851  CALL rflo_getdimensphysnodes( region,ilev,ipnbeg,ipnend, &
1852  jpnbeg,jpnend,kpnbeg,kpnend )
1853  CALL rflo_getcelloffset( region,ilev,icoff,ijcoff )
1854  CALL rflo_getnodeoffset( region,ilev,inoff,ijnoff )
1855  ibn = indijk(idnbeg,jdnbeg,kdnbeg,inoff,ijnoff)
1856  ien = indijk(idnend,jdnend,kdnend,inoff,ijnoff)
1857 
1858  xyz => region%levels(ilev)%grid%xyz
1859  xyzorig => region%levels(ilev)%gridOld%xyz
1860  xyzold => region%levels(ilev)%grid%xyzOld
1861  xyzorth => region%levels(ilev)%grid%xyzOrth
1862  xyztemp => region%levels(ilev)%grid%xyzTemp
1863 
1864  ALLOCATE( xyzwork(3,ibn:ien),stat=errorflag )
1865  global%error = errorflag
1866  IF (global%error /= 0) CALL errorstop( global,err_allocate,&
1867  __line__ )
1868 
1869  xyzwork = 0._rfreal
1870 
1871 ! reset motion vectors --------------------------------------------------------
1872 
1873  xyztemp = xyzold
1874  DO ijk=ibn,ien
1875  xyzold(:,ijk) = xyzold(:,ijk) + xyzorig(:,ijk)
1876  xyzorth(:,ijk) = xyz( :,ijk) + xyzorig(:,ijk)
1877  ENDDO
1878 
1879 !if (region%iRegionGlobal==1) write(*,*)'xyzOrth0',xyzOrth(:,IndIJK(2,2,2 ,iNOff,ijNOff))
1880 
1881 ! reset motion vectors --------------------------------------------------------
1882 
1883  DO k=kpnbeg+1,kpnend-1
1884  DO j=jpnbeg+1,jpnend-1
1885  DO i=ipnbeg+1,ipnend-1
1886 
1887 ! ----- I-Face
1888 
1889  ijkn = indijk(i ,j ,k ,inoff,ijnoff)
1890  ijkn1 = indijk(i ,j+1,k ,inoff,ijnoff)
1891  ijkn2 = indijk(i ,j ,k+1,inoff,ijnoff)
1892  ijkn3 = indijk(i ,j-1,k ,inoff,ijnoff)
1893  ijkn4 = indijk(i ,j ,k-1,inoff,ijnoff)
1894 
1895  CALL accumulatexyzwork
1896 ! if (i==2 .AND. j==4 .AND. k==2) &
1897 ! write(*,*)'xyzOrthI',xyzWork(:,IndIJK(2,4,2 ,iNOff,ijNOff))/4 &
1898 ! -xyzOrig(:,IndIJK(2,4,2 ,iNOff,ijNOff))
1899 
1900 ! ----- J-Face
1901 
1902  ijkn = indijk(i ,j ,k ,inoff,ijnoff)
1903  ijkn1 = indijk(i ,j ,k+1,inoff,ijnoff)
1904  ijkn2 = indijk(i+1,j ,k ,inoff,ijnoff)
1905  ijkn3 = indijk(i ,j ,k-1,inoff,ijnoff)
1906  ijkn4 = indijk(i-1,j ,k ,inoff,ijnoff)
1907 
1908  CALL accumulatexyzwork
1909 ! if (i==2 .AND. j==4 .AND. k==2) &
1910 ! write(*,*)'xyzOrthJ',xyzWork(:,IndIJK(2,4,2 ,iNOff,ijNOff))/8 &
1911 ! -xyzOrig(:,IndIJK(2,4,2 ,iNOff,ijNOff))
1912 
1913 ! ----- K-Face
1914 
1915  ijkn = indijk(i ,j ,k ,inoff,ijnoff)
1916  ijkn1 = indijk(i+1,j ,k ,inoff,ijnoff)
1917  ijkn2 = indijk(i ,j+1,k ,inoff,ijnoff)
1918  ijkn3 = indijk(i-1,j ,k ,inoff,ijnoff)
1919  ijkn4 = indijk(i ,j-1,k ,inoff,ijnoff)
1920 
1921  CALL accumulatexyzwork
1922 ! if (i==2 .AND. j==4 .AND. k==2) &
1923 ! write(*,*)'xyzOrthK',xyzWork(:,IndIJK(2,4,2 ,iNOff,ijNOff))/12 &
1924 ! -xyzOrig(:,IndIJK(2,4,2 ,iNOff,ijNOff))
1925 
1926  ENDDO ! i
1927  ENDDO ! j
1928  ENDDO ! k
1929 
1930 ! average accumulated grid coordinates
1931 
1932  DO k=kpnbeg+1,kpnend-1
1933  DO j=jpnbeg+1,jpnend-1
1934  DO i=ipnbeg+1,ipnend-1
1935  ijk = indijk(i ,j ,k ,inoff,ijnoff)
1936  xyzorth(:,ijk) = rnfac*xyzwork(:,ijk) - xyzorig(:,ijk)
1937  xyz( :,ijk) = xyzorth(:,ijk)
1938  ENDDO ! i
1939  ENDDO ! j
1940  ENDDO ! k
1941 
1942  xyzold = xyztemp
1943 
1944 !if (region%iRegionGlobal==1) write(*,*)'xyzOrth1',xyz(:,IndIJK(2,2,2 ,iNOff,ijNOff)) + xyzOrig(:,IndIJK(2,2,2 ,iNOff,ijNOff))
1945 !if (region%iRegionGlobal==1) write(*,*)'xyzOrth2',xyzOld(:,IndIJK(2,2,2 ,iNOff,ijNOff)) + xyzOrig(:,IndIJK(2,2,2 ,iNOff,ijNOff))
1946 
1947 ! deallocate temporary memory -------------------------------------------------
1948 
1949  DEALLOCATE( xyzwork,stat=errorflag )
1950  global%error = errorflag
1951  IF (global%error /= 0) CALL errorstop( global,err_deallocate,&
1952  __line__ )
1953 
1954 ! finalize --------------------------------------------------------------------
1955 
1956  CALL deregisterfunction( global )
1957 
1958 CONTAINS
1959 
1960 ! ==============================================================================
1961 ! subroutine for accumulation of projected cell corners coordinates
1962 ! ==============================================================================
1963 
1964  SUBROUTINE accumulatexyzwork
1965 
1966  IMPLICIT NONE
1967 
1968  delt(:) = xyzorth(:,ijkn)-xyzold(:,ijkn)
1969 
1970  f1(:) = xyzold(:,ijkn1)-xyzold(:,ijkn)
1971  f2(:) = xyzold(:,ijkn2)-xyzold(:,ijkn)
1972  a11 = dot_product( f1(:),xyzold(:,ijkn1) )
1973  a12 = dot_product( f2(:),xyzold(:,ijkn1) )
1974  a21 = dot_product( f1(:),xyzold(:,ijkn2) )
1975  a22 = dot_product( f2(:),xyzold(:,ijkn2) )
1976  rh1 = dot_product( delt(:),xyzold(:,ijkn1) )
1977  rh2 = dot_product( delt(:),xyzold(:,ijkn2) )
1978  rja = 1._rfreal/(a11*a22-a12*a21+eps)
1979  c1 = (a22*rh1 - a12*rh2)*rja
1980  c2 = (a11*rh2 - a21*rh1)*rja
1981  xyzp(:) = ocell*(xyzold(:,ijkn) + c1*f1(:) + c2*f2(:)) + &
1982  (1._rfreal-ocell)*xyzorth(:,ijkn)
1983  xyzwork(:,ijkn) = xyzwork(:,ijkn) + xyzp(:)
1984 
1985  f1(:) = xyzold(:,ijkn2)-xyzold(:,ijkn)
1986  f2(:) = xyzold(:,ijkn3)-xyzold(:,ijkn)
1987  a11 = dot_product( f1(:),xyzold(:,ijkn2) )
1988  a12 = dot_product( f2(:),xyzold(:,ijkn2) )
1989  a21 = dot_product( f1(:),xyzold(:,ijkn3) )
1990  a22 = dot_product( f2(:),xyzold(:,ijkn3) )
1991  rh1 = dot_product( delt(:),xyzold(:,ijkn2) )
1992  rh2 = dot_product( delt(:),xyzold(:,ijkn3) )
1993  rja = 1._rfreal/(a11*a22-a12*a21+eps)
1994  c1 = (a22*rh1 - a12*rh2)*rja
1995  c2 = (a11*rh2 - a21*rh1)*rja
1996  xyzp(:) = ocell*(xyzold(:,ijkn) + c1*f1(:) + c2*f2(:)) + &
1997  (1._rfreal-ocell)*xyzorth(:,ijkn)
1998  xyzwork(:,ijkn) = xyzwork(:,ijkn) + xyzp(:)
1999 
2000  f1(:) = xyzold(:,ijkn3)-xyzold(:,ijkn)
2001  f2(:) = xyzold(:,ijkn4)-xyzold(:,ijkn)
2002  a11 = dot_product( f1(:),xyzold(:,ijkn3) )
2003  a12 = dot_product( f2(:),xyzold(:,ijkn3) )
2004  a21 = dot_product( f1(:),xyzold(:,ijkn4) )
2005  a22 = dot_product( f2(:),xyzold(:,ijkn4) )
2006  rh1 = dot_product( delt(:),xyzold(:,ijkn3) )
2007  rh2 = dot_product( delt(:),xyzold(:,ijkn4) )
2008  rja = 1._rfreal/(a11*a22-a12*a21+eps)
2009  c1 = (a22*rh1 - a12*rh2)*rja
2010  c2 = (a11*rh2 - a21*rh1)*rja
2011  xyzp(:) = ocell*(xyzold(:,ijkn) + c1*f1(:) + c2*f2(:)) + &
2012  (1._rfreal-ocell)*xyzorth(:,ijkn)
2013  xyzwork(:,ijkn) = xyzwork(:,ijkn) + xyzp(:)
2014 
2015  f1(:) = xyzold(:,ijkn4)-xyzold(:,ijkn)
2016  f2(:) = xyzold(:,ijkn1)-xyzold(:,ijkn)
2017  a11 = dot_product( f1(:),xyzold(:,ijkn4) )
2018  a12 = dot_product( f2(:),xyzold(:,ijkn4) )
2019  a21 = dot_product( f1(:),xyzold(:,ijkn1) )
2020  a22 = dot_product( f2(:),xyzold(:,ijkn1) )
2021  rh1 = dot_product( delt(:),xyzold(:,ijkn4) )
2022  rh2 = dot_product( delt(:),xyzold(:,ijkn1) )
2023  rja = 1._rfreal/(a11*a22-a12*a21+eps)
2024  c1 = (a22*rh1 - a12*rh2)*rja
2025  c2 = (a11*rh2 - a21*rh1)*rja
2026  xyzp(:) = ocell*(xyzold(:,ijkn) + c1*f1(:) + c2*f2(:)) + &
2027  (1._rfreal-ocell)*xyzorth(:,ijkn)
2028  xyzwork(:,ijkn) = xyzwork(:,ijkn) + xyzp(:)
2029 
2030  END SUBROUTINE accumulatexyzwork
2031 
2032 END SUBROUTINE rflo_laplacegridortho2
2033 
2034 !******************************************************************************
2035 !
2036 ! Purpose:
2037 !
2038 ! Description: none.
2039 !
2040 ! Input: region = data of current region, old grid coordinates.
2041 !
2042 ! Output: region%levels%grid%xyz = grid movements.
2043 !
2044 ! Notes: on entry and exit, xyz holds the grid motion
2045 !
2046 !******************************************************************************
2047 
2048 SUBROUTINE rflo_projectquadcorner( ocell,eps,cofgp,cofgm,cfc,xyzc,xyzp )
2049 
2050  IMPLICIT NONE
2051 
2052 ! ... parameters
2053  REAL(RFREAL) :: ocell, eps
2054  REAL(RFREAL) :: cofgp(:), cofgm(:), cfc(:), xyzc(:), xyzp(:)
2055 
2056 ! ... loop variables
2057 
2058 ! ... local variables
2059  REAL(RFREAL) :: cfcgp(xcoord:zcoord), cfcgm(xcoord:zcoord)
2060  REAL(RFREAL) :: sqrp, sqrm, cp, cm, fac
2061 
2062 !******************************************************************************
2063 
2064  cfcgp(:) = cofgp(:) - cfc(:)
2065  cfcgm(:) = cofgm(:) - cfc(:)
2066 
2067  sqrp = cfcgp(xcoord)*cfcgp(xcoord) + &
2068  cfcgp(ycoord)*cfcgp(ycoord) + &
2069  cfcgp(zcoord)*cfcgp(zcoord)
2070 
2071  sqrm = cfcgm(xcoord)*cfcgm(xcoord) + &
2072  cfcgm(ycoord)*cfcgm(ycoord) + &
2073  cfcgm(zcoord)*cfcgm(zcoord)
2074 
2075  cp = dot_product( (cfc(:)-xyzc(:)),cfcgp(:) )/(sqrp + epsilon( 1._rfreal ))
2076  cm = dot_product( (cfc(:)-xyzc(:)),cfcgm(:) )/(sqrm + epsilon( 1._rfreal ))
2077 
2078  fac = 0.5_rfreal
2079 ! IF (cp < eps .OR. cm < eps ) fac = 1._RFREAL
2080  xyzp(:) = xyzc(:) + ocell*fac*(cp*cfcgp(:) + cm*cfcgm(:))
2081 ! xyzp(:) = xyzc(:) + ocell*(ABS(cp)*cp*cfcgp(:) + ABS(cm)*cm*cfcgm(:))/ &
2082 ! (ABS(cp)+ABS(cm))
2083 
2084 ! finalize --------------------------------------------------------------------
2085 
2086 END SUBROUTINE rflo_projectquadcorner
2087 
2088 ! ******************************************************************************
2089 ! End
2090 ! ******************************************************************************
2091 
2092 END MODULE rflo_modlaplacesmoothing
2093 
2094 ! ******************************************************************************
2095 !
2096 ! RCS Revision history:
2097 !
2098 ! $Log: RFLO_ModLaplaceSmoothing.F90,v $
2099 ! Revision 1.13 2009/08/27 14:04:50 mtcampbe
2100 ! Updated to enable burning motion with symmetry boundaries and enhanced
2101 ! burnout code.
2102 !
2103 ! Revision 1.12 2008/12/06 08:44:16 mtcampbe
2104 ! Updated license.
2105 !
2106 ! Revision 1.11 2008/11/19 22:17:27 mtcampbe
2107 ! Added Illinois Open Source License/Copyright
2108 !
2109 ! Revision 1.10 2006/03/05 21:51:34 wasistho
2110 ! changed computational space coordinates to be based on initial grid
2111 !
2112 ! Revision 1.9 2005/11/28 20:39:16 wasistho
2113 ! removed file-Id in rflo_projectQuadCorner
2114 !
2115 ! Revision 1.8 2005/11/28 18:47:19 wasistho
2116 ! outcommented print statement in gridOrtho2
2117 !
2118 ! Revision 1.7 2005/11/28 18:29:56 gzheng
2119 ! removed duplicated RFLO_CalcFaceVectors in USE statement, which broke gnu f95 compiler, also fixed 3 single line statement which is too long for some compilers.
2120 !
2121 ! Revision 1.6 2005/11/17 00:46:26 wasistho
2122 ! dont touch degenerate points in moveGridSolve
2123 !
2124 ! Revision 1.5 2005/11/16 20:45:30 wasistho
2125 ! update on MoveGridFoms (3), still less robust than (2)
2126 !
2127 ! Revision 1.4 2005/11/01 09:01:58 wasistho
2128 ! added volume TFI for movement jump
2129 !
2130 ! Revision 1.3 2005/10/28 07:36:33 wasistho
2131 ! gridOrtho is private
2132 !
2133 ! Revision 1.2 2005/10/27 06:29:46 wasistho
2134 ! increased denominator by epsilon in ProjectQuadCorner
2135 !
2136 ! Revision 1.1 2005/10/27 06:02:46 wasistho
2137 ! initial import
2138 !
2139 !
2140 !
2141 ! ******************************************************************************
2142 
2143 
2144 
2145 
2146 
2147 
2148 
2149 
2150 
2151 
2152 
2153 
2154 
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_copygeometrydummy
subroutine rflo_copygeometrydummy(region)
Definition: RFLO_CopyGeometryDummy.F90:45
idcend
**********************************************************************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 idcend
Definition: RFLO_GetDimensDummy.h:27
jpnbeg
**********************************************************************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 jpnbeg
Definition: RFLO_GetDimensPhysNodes.h:27
rflo_calccellcentroids
subroutine rflo_calccellcentroids(region)
Definition: RFLO_CalcCellCentroids.F90:43
d
const NT & d
Definition: rational_rotation.h:73
k
j indices k indices k
Definition: Indexing.h:6
dx
NT dx
Definition: rational_rotation.h:81
rflo_modlaplacesmoothing::rflo_laplacegridortho2
subroutine rflo_laplacegridortho2(region)
Definition: RFLO_ModLaplaceSmoothing.F90:1803
kpnbeg
**********************************************************************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 kpnbeg
Definition: RFLO_GetDimensPhysNodes.h:27
rflo_modlaplacesmoothing::rflo_laplacegridortho1
subroutine rflo_laplacegridortho1(region)
Definition: RFLO_ModLaplaceSmoothing.F90:1548
moderror::registerfunction
subroutine registerfunction(global, funName, fileName)
Definition: ModError.F90:449
rflo_extrapolategeometry
subroutine rflo_extrapolategeometry(region)
Definition: RFLO_ExtrapolateGeometry.F90:43
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
rflo_exchangednoderecv
subroutine rflo_exchangednoderecv(region, regionSrc, patch, patchSrc, average, dNode)
Definition: RFLO_ExchangeDnodeRecv.F90:51
modgrid::t_grid
Definition: ModGrid.F90:56
sqrt
double sqrt(double d)
Definition: double.h:73
rflo_modlaplacesmoothing::rflo_laplacegridortho
subroutine rflo_laplacegridortho(region)
Definition: RFLO_ModLaplaceSmoothing.F90:1135
rflo_changeinteriorgrid
subroutine rflo_changeinteriorgrid(region, boundMoved, edgeMoved, arcLen12, arcLen34, arcLen56, xyzOld, xyz)
Definition: RFLO_ChangeInteriorGrid.F90:51
rflo_modlaplacesmoothing::rflo_projectquadcorner
subroutine rflo_projectquadcorner(ocell, eps, cofgp, cofgm, cfc, xyzc, xyzp)
Definition: RFLO_ModLaplaceSmoothing.F90:2048
TestBase::f2
virtual int f2()
Definition: meshtest1.C:4115
TestBase::f1
virtual int f1()
Definition: meshtest1.C:4114
jdnbeg
**********************************************************************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 jdnbeg
Definition: RFLO_GetDimensDummyNodes.h:27
rflo_modlaplacesmoothing::rflo_laplacegridjump
subroutine, public rflo_laplacegridjump(region, patch)
Definition: RFLO_ModLaplaceSmoothing.F90:592
jpnend
**********************************************************************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 jpnend
Definition: RFLO_GetDimensPhysNodes.h:27
idnend
**********************************************************************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 idnend
Definition: RFLO_GetDimensDummyNodes.h:27
rflo_getnodeoffset
subroutine rflo_getnodeoffset(region, iLev, iNodeOffset, ijNodeOffset)
Definition: RFLO_GetDimensions.F90:267
jdnend
**********************************************************************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 jdnend
Definition: RFLO_GetDimensDummyNodes.h:27
rflo_getdimensdummy
subroutine rflo_getdimensdummy(region, iLev, idcbeg, idcend, jdcbeg, jdcend, kdcbeg, kdcend)
Definition: RFLO_GetDimensions.F90:67
idnbeg
**********************************************************************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 idnbeg
Definition: RFLO_GetDimensDummyNodes.h:27
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
kdcbeg
**********************************************************************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 kdcbeg
Definition: RFLO_GetDimensDummy.h:27
rflo_modlaplacesmoothing
Definition: RFLO_ModLaplaceSmoothing.F90:37
rflo_exchangegeometry
subroutine rflo_exchangegeometry(regions)
Definition: RFLO_ExchangeGeometry.F90:44
rflo_generatecoarsegrids
subroutine rflo_generatecoarsegrids(region)
Definition: RFLO_GenerateCoarseGrids.F90:44
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
ipnbeg
**********************************************************************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 ipnbeg
Definition: RFLO_GetDimensPhysNodes.h:27
rflo_getcelloffset
subroutine rflo_getcelloffset(region, iLev, iCellOffset, ijCellOffset)
Definition: RFLO_GetDimensions.F90:235
idcbeg
**********************************************************************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 idcbeg
Definition: RFLO_GetDimensDummy.h:27
rflo_clearsendrequests
subroutine rflo_clearsendrequests(regions, iReg, geometry)
Definition: RFLO_ClearSendRequests.F90:46
Direction_3::dz
RT dz() const
Definition: Direction_3.h:133
rflo_laplacegridpatch
subroutine rflo_laplacegridpatch(region, patch)
Definition: RFLO_LaplaceGridPatch.F90:44
jdcend
**********************************************************************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 jdcend
Definition: RFLO_GetDimensDummy.h:27
rflo_getdimensphysnodes
subroutine rflo_getdimensphysnodes(region, iLev, ipnbeg, ipnend, jpnbeg, jpnend, kpnbeg, kpnend)
Definition: RFLO_GetDimensions.F90:193
j
j indices j
Definition: Indexing.h:6
dy
NT dy
Definition: rational_rotation.h:82
rflo_laplacegridsolve
subroutine rflo_laplacegridsolve(region)
Definition: RFLO_LaplaceGridSolve.F90:45
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
rflo_exchangednodecopy
subroutine rflo_exchangednodecopy(region, regionSrc, patch, patchSrc, average, dNode, dNodeSrc)
Definition: RFLO_ExchangeDnodeCopy.F90:54
jdcbeg
**********************************************************************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 jdcbeg
Definition: RFLO_GetDimensDummy.h:27
rflo_getdimensdummynodes
subroutine rflo_getdimensdummynodes(region, iLev, idnbeg, idnend, jdnbeg, jdnend, kdnbeg, kdnend)
Definition: RFLO_GetDimensions.F90:151
modglobal::t_global
Definition: ModGlobal.F90:50
moderror::errorstop
subroutine errorstop(global, errorCode, errorLine, addMessage)
Definition: ModError.F90:483
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
dist
long double dist(long double *coord1, long double *coord2, int size)
Definition: utilities/rc_cmp/util.cpp:305
dot_product
long double dot_product(pnt vec1, pnt vec2)
Definition: utilities/rc_cmp/util.cpp:387
setup_py::grid
subroutine grid(bp)
Definition: setup_py.f90:257
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
ipnend
**********************************************************************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 ipnend
Definition: RFLO_GetDimensPhysNodes.h:27
rflo_exchangednodesend
subroutine rflo_exchangednodesend(region, regionSrc, patch, dNode)
Definition: RFLO_ExchangeDnodeSend.F90:48
denom
CGAL_BEGIN_NAMESPACE void const NT NT NT NT & denom
Definition: rational_rotation.h:61
kdnbeg
**********************************************************************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 kdnbeg
Definition: RFLO_GetDimensDummyNodes.h:27
rflo_laplacegridsmoo
subroutine rflo_laplacegridsmoo(regions, resid)
Definition: RFLO_LaplaceGridSmoo.F90:45
rflo_calcfacecentroids
subroutine rflo_calcfacecentroids(region)
Definition: RFLO_CalcFaceCentroids.F90:42
modbndpatch::t_patch
Definition: ModBndPatch.F90:91
accumulatexyzwork
subroutine accumulatexyzwork(sf, fc1, fc2, fc3, fc4)
Definition: RFLO_ModLaplaceSmoothing.F90:1486