From 0274a29acec2868f777ec758ea4ea72e0c6296fa Mon Sep 17 00:00:00 2001
From: odemirel <odemirel@7c7fe9aa-52eb-4d9e-b0a8-ba7d787348e9>
Date: Mon, 8 Mar 2010 13:38:59 +0000
Subject: [PATCH] doc is updated for ppm_gmm_* inside ppmnumerics. all the
other files inside ppmnumerics should follow the rules described in
[wiki:Coding_Guide] and take ppm_gmm_* as example.
git-svn-id: https://ppm.inf.ethz.ch/svn/ppmnumerics/branches/ngtopo/libppmnumerics@539 7c7fe9aa-52eb-4d9e-b0a8-ba7d787348e9
---
src/ppm_gmm_cpt.f | 179 ++++++++++++-------------------------
src/ppm_gmm_extend.f | 105 ++++++++--------------
src/ppm_gmm_extend_fwd.f | 69 ++++++--------
src/ppm_gmm_finalize.f | 10 +--
src/ppm_gmm_init.f | 49 +++++-----
src/ppm_gmm_kickoff.f | 5 ++
src/ppm_gmm_reinitialize.f | 161 ++++++++++-----------------------
7 files changed, 201 insertions(+), 377 deletions(-)
diff --git a/src/ppm_gmm_cpt.f b/src/ppm_gmm_cpt.f
index e3f0ddb..18c7f8f 100644
--- a/src/ppm_gmm_cpt.f
+++ b/src/ppm_gmm_cpt.f
@@ -1,94 +1,8 @@
!-------------------------------------------------------------------------
! Subroutine : ppm_gmm_cpt
!-------------------------------------------------------------------------
- !
- ! Purpose : This routine performs a closest point transform.
- ! For each grid point adjacent to the interface,
- ! the closest point ON the interface is returned.
- ! ppm_gmm_init must be called BEFORE this routine is
- ! invoked.
- !
- ! Input : fdata(...) (F) Level data. Either rank 3
- ! (for 2D scalar fields), or rank
- ! 4 (for 3D scalar fields).
- ! Indices: (i,j,[k],isub).
- ! Every zero-crossing is
- ! interpreted as an interface.
- ! A ghostsize of 1 is needed
- ! on all sides which must be
- ! filled with the old level
- ! function value on input!! Only
- ! scalar fdata supported.
- ! tol (F) Relative tolerance for the
- ! determined distance to the
- ! interface. 1E-3 is a good
- ! choice. The tolerance is in
- ! multiples of grid spacings.
- ! chi([:],:,:,:,:)(F) rank 5 (3d) or rank 4 (2d)
- ! field specifying the positions
- ! of the grid nodes. 1st index:
- ! 1..ppm_dim, then i,j,[k],isub.
- ! OPTIONAL. Uniform grid is
- ! assumed if absent. Ghostlayers
- ! of size >=1 must be pre-filled.
- !
- ! Input/output :
- !
- ! Output : npts (I) Number of unique grid points
- ! adjacent to the interface.
- ! ipts(:,:) (I) Mesh indices of these
- ! points. 1st index:
- ! i,j,(k),isub (local sub ID);
- ! 2nd: 1...npts. Will be
- ! allocated by this routine.
- ! closest(:,:) (F) Locations of the closest
- ! points ON the interface. 1st
- ! index: x,y,(z), 2nd:
- ! 1...npts. Will be allocated
- ! by this routine.
- ! info (I) return status. 0 on success.
- !
- ! Remarks :
- !
- ! References :
- !
- ! Revisions :
- !-------------------------------------------------------------------------
- ! $Log: ppm_gmm_cpt.f,v $
- ! Revision 1.1.1.1 2007/07/13 10:18:55 ivos
- ! CBL version of the PPM library
- !
- ! Revision 1.8 2005/07/14 19:58:11 ivos
- ! Added OPTIONAL argument chi for mesh node positions in distorted
- ! (mapped) meshes. For use with AGM for example.
- !
- ! Revision 1.7 2005/06/06 20:33:29 ivos
- ! Added pointer NULLIFication at the end to restore alloccation status.
- !
- ! Revision 1.6 2005/06/01 05:19:05 ivos
- ! Updated to new interface of gmm_kickoff.
- !
- ! Revision 1.5 2005/04/27 01:06:10 ivos
- ! Convergence tests completed, cleaned up code, optmized code (Shark),
- ! and changed structure to allow faster compilation.
- !
- ! Revision 1.4 2005/04/21 04:49:56 ivos
- ! bugfix: pointers to array slabs were incorrectly used.
- !
- ! Revision 1.3 2005/03/16 06:20:07 ivos
- ! Several bugfixes. 1st order version is now tested. Moved all large
- ! data to the module.
- !
- ! Revision 1.2 2005/03/12 04:08:34 ivos
- ! Misc bug fixes.
- !
- ! Revision 1.1 2005/03/11 21:09:10 ivos
- ! Initial implementation.
- !
- !-------------------------------------------------------------------------
! Parallel Particle Mesh Library (PPM)
- ! Institute of Computational Science
- ! ETH Zentrum, Hirschengraben 84
+ ! ETH Zurich
! CH-8092 Zurich, Switzerland
!-------------------------------------------------------------------------
#if __DIM == __2D
@@ -109,6 +23,9 @@
& info,chi)
#endif
#endif
+ !!! This routine performs a closest point transform. For each grid point
+ !!! adjacent to the interface, the closest point ON the interface is
+ !!! returned. ppm_gmm_init must be called BEFORE this routine is invoked.
!-------------------------------------------------------------------------
! Includes
!-------------------------------------------------------------------------
@@ -120,11 +37,6 @@
USE ppm_module_data_mesh
USE ppm_module_data_gmm
USE ppm_module_gmm_kickoff
- USE ppm_module_substart
- USE ppm_module_substop
- USE ppm_module_error
- USE ppm_module_alloc
- USE ppm_module_typedef
IMPLICIT NONE
#if __KIND == __SINGLE_PRECISION | __KIND == __SINGLE_PRECISION_COMPLEX
INTEGER, PARAMETER :: MK = ppm_kind_single
@@ -136,15 +48,40 @@
!-------------------------------------------------------------------------
#if __DIM == __2D
REAL(MK), DIMENSION(:,:,:) , POINTER :: fdata
+ !!! Level data. Rank 3 (for 2D scalar fields),
+ !!! Indices: (i,j,[k],isub). Every zero-crossing is interpreted as an
+ !!! interface. A ghostsize of 1 is needed on all sides which must be
+ !!! filled with the old level function value on input!! Only scalar fdata
+ !!! supported.
REAL(MK), DIMENSION(:,:,:,:) , INTENT(IN), OPTIONAL:: chi
+ !!! Rank 4 (2d) field specifying the positions of the grid nodes.
+ !!! 1st index: 1..ppm_dim, then i,j,[k],isub. OPTIONAL. Uniform grid is
+ !!! assumed if absent. Ghostlayers of size >=1 must be pre-filled.
#elif __DIM == __3D
REAL(MK), DIMENSION(:,:,:,:) , POINTER :: fdata
+ !!! Level data. Rank 4 (for 3D scalar fields).
+ !!! Indices: (i,j,[k],isub). Every zero-crossing is interpreted as an
+ !!! interface. A ghostsize of 1 is needed on all sides which must be
+ !!! filled with the old level function value on input!! Only scalar fdata
+ !!! supported.
REAL(MK), DIMENSION(:,:,:,:,:) , INTENT(IN), OPTIONAL:: chi
+ !!! Rank 5 (3d) field specifying the positions of the grid nodes.
+ !!! 1st index: 1..ppm_dim, then i,j,[k],isub. OPTIONAL. Uniform grid is
+ !!! assumed if absent. Ghostlayers of size >=1 must be pre-filled.
#endif
INTEGER , DIMENSION(:,:) , POINTER :: ipts
+ !!! Mesh indices of these points. 1st index: i,j,(k),isub (local sub ID);
+ !!! 2nd: 1...npts. Will be allocated by this routine.
REAL(MK), DIMENSION(:,:) , POINTER :: closest
+ !!! Locations of the closest points ON the interface. 1st index: x,y,(z),
+ !!! 2nd: 1...npts. Will be allocated by this routine.
REAL(MK) , INTENT(IN ) :: tol
- INTEGER , INTENT( OUT) :: info,npts
+ !!! Relative tolerance for the determined distance to the interface.
+ !!! 1E-3 is a good choice. The tolerance is in multiples of grid spacings.
+ INTEGER , INTENT( OUT) :: info
+ !!! Return status. 0 upon success
+ INTEGER , INTENT( OUT) :: npts
+ !!! Number of unique grid points adjacent to the interface.
!-------------------------------------------------------------------------
! Local variables
!-------------------------------------------------------------------------
@@ -155,14 +92,10 @@
REAL(MK) :: t0,x,y,z,xx,yy,zz,dx,dy,dz
REAL(MK) :: s,sprev,thresh
LOGICAL :: lok
- TYPE(ppm_t_topo), POINTER :: topo
- TYPE(ppm_t_equi_mesh), POINTER :: mesh
!-------------------------------------------------------------------------
! Initialise
!-------------------------------------------------------------------------
CALL substart('ppm_gmm_cpt',t0,info)
- topo => ppm_topo(gmm_topoid)%t
- mesh => topo%mesh(gmm_meshid)
#if __KIND == __SINGLE_PRECISION
clotmp => gmm_clos
#elif __KIND == __DOUBLE_PRECISION
@@ -185,7 +118,7 @@
GOTO 9999
ENDIF
#if __DIM == __3D
- IF (SIZE(fdata,4) .LT. topo%nsublist) THEN
+ IF (SIZE(fdata,4) .LT. ppm_nsublist(gmm_topoid)) THEN
info = ppm_error_error
CALL ppm_error(ppm_err_argument,'ppm_gmm_cpt', &
& 'field data for some subs is missing',__LINE__,info)
@@ -210,7 +143,7 @@
GOTO 9999
ENDIF
#elif __DIM == __2D
- IF (SIZE(fdata,3) .LT. topo%nsublist) THEN
+ IF (SIZE(fdata,3) .LT. ppm_nsublist(gmm_topoid)) THEN
info = ppm_error_error
CALL ppm_error(ppm_err_argument,'ppm_gmm_cpt', &
& 'field data for some subs is missing',__LINE__,info)
@@ -234,23 +167,23 @@
! Find mesh spacing
!-------------------------------------------------------------------------
IF (ppm_kind .EQ. ppm_kind_single) THEN
- dx = (topo%max_physs(1)-topo%min_physs(1))/ &
- & REAL(mesh%Nm(1)-1,ppm_kind_single)
- dy = (topo%max_physs(2)-topo%min_physs(2))/ &
- & REAL(mesh%Nm(2)-1,ppm_kind_single)
+ dx = (ppm_max_physs(1,gmm_topoid)-ppm_min_physs(1,gmm_topoid))/ &
+ & REAL(ppm_cart_mesh(gmm_meshid,gmm_topoid)%Nm(1)-1,ppm_kind_single)
+ dy = (ppm_max_physs(2,gmm_topoid)-ppm_min_physs(2,gmm_topoid))/ &
+ & REAL(ppm_cart_mesh(gmm_meshid,gmm_topoid)%Nm(2)-1,ppm_kind_single)
IF (ppm_dim .GT. 2) THEN
- dz = (topo%max_physs(3)-topo%min_physs(3))/ &
- & REAL(mesh%Nm(3)-1, &
+ dz = (ppm_max_physs(3,gmm_topoid)-ppm_min_physs(3,gmm_topoid))/ &
+ & REAL(ppm_cart_mesh(gmm_meshid,gmm_topoid)%Nm(3)-1, &
& ppm_kind_single)
ENDIF
ELSE
- dx = (topo%max_physs(1)-topo%min_physs(1))/ &
- & REAL(mesh%Nm(1)-1,ppm_kind_double)
- dy = (topo%max_physs(2)-topo%min_physs(2))/ &
- & REAL(mesh%Nm(2)-1,ppm_kind_double)
+ dx = (ppm_max_physd(1,gmm_topoid)-ppm_min_physd(1,gmm_topoid))/ &
+ & REAL(ppm_cart_mesh(gmm_meshid,gmm_topoid)%Nm(1)-1,ppm_kind_double)
+ dy = (ppm_max_physd(2,gmm_topoid)-ppm_min_physd(2,gmm_topoid))/ &
+ & REAL(ppm_cart_mesh(gmm_meshid,gmm_topoid)%Nm(2)-1,ppm_kind_double)
IF (ppm_dim .GT. 2) THEN
- dz = (topo%max_physs(3)-topo%min_physs(3))/ &
- & REAL(mesh%Nm(3)-1, &
+ dz = (ppm_max_physd(3,gmm_topoid)-ppm_min_physd(3,gmm_topoid))/ &
+ & REAL(ppm_cart_mesh(gmm_meshid,gmm_topoid)%Nm(3)-1, &
& ppm_kind_double)
ENDIF
ENDIF
@@ -352,20 +285,20 @@
sprev = HUGE(sprev)
ENDIF
isub = ipts(4,npts)
- jsub = topo%isublist(isub)
+ jsub = ppm_isublist(isub,gmm_topoid)
IF (PRESENT(chi)) THEN
x = chi(1,ipts(1,npts),ipts(2,npts),ipts(3,npts),isub)
y = chi(2,ipts(1,npts),ipts(2,npts),ipts(3,npts),isub)
z = chi(3,ipts(1,npts),ipts(2,npts),ipts(3,npts),isub)
ELSE
IF (ppm_kind .EQ. ppm_kind_single) THEN
- x = topo%min_subs(1,jsub) + (ipts(1,npts)-1)*dx
- y = topo%min_subs(2,jsub)+ (ipts(2,npts)-1)*dy
- z = topo%min_subs(3,jsub)+ (ipts(3,npts)-1)*dz
+ x = ppm_min_subs(1,jsub,gmm_topoid) + (ipts(1,npts)-1)*dx
+ y = ppm_min_subs(2,jsub,gmm_topoid) + (ipts(2,npts)-1)*dy
+ z = ppm_min_subs(3,jsub,gmm_topoid) + (ipts(3,npts)-1)*dz
ELSE
- x = topo%min_subd(1,jsub) + (ipts(1,npts)-1)*dx
- y = topo%min_subd(2,jsub) + (ipts(2,npts)-1)*dy
- z = topo%min_subd(3,jsub) + (ipts(3,npts)-1)*dz
+ x = ppm_min_subd(1,jsub,gmm_topoid) + (ipts(1,npts)-1)*dx
+ y = ppm_min_subd(2,jsub,gmm_topoid) + (ipts(2,npts)-1)*dy
+ z = ppm_min_subd(3,jsub,gmm_topoid) + (ipts(3,npts)-1)*dz
ENDIF
ENDIF
xx = clotmp(1,idx(i))
@@ -390,17 +323,17 @@
sprev = HUGE(sprev)
ENDIF
isub = ipts(3,npts)
- jsub = topo%isublist(isub)
+ jsub = ppm_isublist(isub,gmm_topoid)
IF (PRESENT(chi)) THEN
x = chi(1,ipts(1,npts),ipts(2,npts),isub)
y = chi(2,ipts(1,npts),ipts(2,npts),isub)
ELSE
IF (ppm_kind .EQ. ppm_kind_single) THEN
- x = topo%min_subs(1,jsub) + (ipts(1,npts)-1)*dx
- y = topo%min_subs(2,jsub)+ (ipts(2,npts)-1)*dy
+ x = ppm_min_subs(1,jsub,gmm_topoid) + (ipts(1,npts)-1)*dx
+ y = ppm_min_subs(2,jsub,gmm_topoid) + (ipts(2,npts)-1)*dy
ELSE
- x = topo%min_subd(1,jsub) + (ipts(1,npts)-1)*dx
- y = topo%min_subd(2,jsub) + (ipts(2,npts)-1)*dy
+ x = ppm_min_subd(1,jsub,gmm_topoid) + (ipts(1,npts)-1)*dx
+ y = ppm_min_subd(2,jsub,gmm_topoid) + (ipts(2,npts)-1)*dy
ENDIF
ENDIF
xx = clotmp(1,idx(i))
diff --git a/src/ppm_gmm_extend.f b/src/ppm_gmm_extend.f
index 50fec31..900a896 100644
--- a/src/ppm_gmm_extend.f
+++ b/src/ppm_gmm_extend.f
@@ -84,7 +84,7 @@
#endif
!!! This routine extends a function defined on the interface to the whole
!!! band on which the level function is defined. The extension is done
- !!! such that the gradient of the function is perpendicular to the
+ !!! such that the graident of the function is perpendicular to the
!!! gradient of the level function. ppm_gmm_init must be called BEFORE
!!! this routine is invoked.
!-------------------------------------------------------------------------
@@ -101,11 +101,6 @@
USE ppm_module_gmm_cpt
USE ppm_module_gmm_march
USE ppm_module_gmm_finalize
- USE ppm_module_substart
- USE ppm_module_substop
- USE ppm_module_error
- USE ppm_module_alloc
- USE ppm_module_typedef
IMPLICIT NONE
#if __KIND == __SINGLE_PRECISION | __KIND == __SINGLE_PRECISION_COMPLEX
INTEGER, PARAMETER :: MK = ppm_kind_single
@@ -134,33 +129,16 @@
END FUNCTION ivalue
END INTERFACE
#endif
- !!! Function pointer to the function computing the value of the function
- !!! on the interface: (F) ivalue(x,y[,z])
- !!! The function may assume that the point (x,y[,z]) is on the interface.
#else
REAL(MK) , INTENT(IN ) :: ivalue
!!! A (F) scalar value defining a cutoff. Points closer to the interface
- !!! than this cutoff will be kept to initialize the marching. The cpt
- !!! is skipped in this case.
+ !!! than this cutoff will be kept to initialize the marching. The cpt is
+ !!! skipped in this case.
#endif
#if __TYPE == __SFIELD
REAL(MK), DIMENSION(:,:,:) , POINTER :: udata
- !!! Field of the function to be extended, defined in the same narrow
- !!! band as the level function. Values outside this band are set to
- !!! HUGE. This has to be allocated to proper size incl. ghost layers of
- !!! size order! If ivalue is a function pointer, udata will be
- !!! completely replaced. If ivalue is a scalar, the points closer to
- !!! the interface than this scalar are kept unchanged. Can be a vector
- !!! or a scalar field. If vector, lda must be given.
#elif __TYPE == __VFIELD
REAL(MK), DIMENSION(:,:,:,:) , POINTER :: udata
- !!! Field of the function to be extended, defined in the same narrow
- !!! band as the level function. Values outside this band are set to
- !!! HUGE. This has to be allocated to proper size incl. ghost layers of
- !!! size order! If ivalue is a function pointer, udata will be
- !!! completely replaced. If ivalue is a scalar, the points closer to
- !!! the interface than this scalar are kept unchanged. Can be a vector
- !!! or a scalar field. If vector, lda must be given.
INTEGER , INTENT(IN ) :: lda
!!! Only present for vector udata. Gives the length of the leading
!!! dimension of udata. All elements will be equally extended.
@@ -171,10 +149,10 @@
!!! to be the interface. Always a scalar field. THIS NEEDS TO BE
!!! PROPERLY ALLOCATED ON INPUT, INCLUDING GHOST LAYERS OF SIZE order.
REAL(MK), DIMENSION(:,:,:,:) , INTENT(IN),OPTIONAL:: chi
- !!! rank 5 (3d) or rank 4 (2d) field specifying the positions
- !!! of the grid nodes. 1st index: 1..ppm_dim, then i,j,[k],isub.
- !!! OPTIONAL. Uniform grid is assumed if absent. Ghostlayers of size
- !!! >=1 must be pre-filled.
+ !!! rank 4 (2d) field specifying the positions of the
+ !!! grid nodes. 1st index: 1..ppm_dim, then i,j,[k],isub. OPTIONAL.
+ !!! Uniform grid is assumed if absent. Ghostlayers of size >=1 must
+ !!! be pre-filled.
#elif __DIM == __3D
#if __KICKOFF == __YES
#if __KIND == __SINGLE_PRECISION
@@ -195,20 +173,20 @@
END INTERFACE
#endif
!!! Function pointer to the function computing the value of the function
- !!! on the interface: (F) ivalue(x,y[,z])
- !!! The function may assume that the point (x,y[,z]) is on the interface.
+ !!! on the interface: (F) ivalue(x,y[,z]) The function may assume that
+ !!! the point (x,y[,z]) is on the interface.
#else
REAL(MK) , INTENT(IN ) :: ivalue
!!! A (F) scalar value defining a cutoff. Points closer to the interface
- !!! than this cutoff will be kept to initialize the marching. The cpt
- !!! is skipped in this case.
+ !!! than this cutoff will be kept to initialize the marching. The cpt is
+ !!! skipped in this case.
#endif
#if __TYPE == __SFIELD
REAL(MK), DIMENSION(:,:,:,:) , POINTER :: udata
!!! Field of the function to be extended, defined in the same narrow
!!! band as the level function. Values outside this band are set to
- !!! HUGE. This has to be allocated to proper size incl. ghost layers of
- !!! size order! If ivalue is a function pointer, udata will be
+ !!! HUGE. This has to be allocated to proper size incl. ghost layers
+ !!! of size order! If ivalue is a function pointer, udata will be
!!! completely replaced. If ivalue is a scalar, the points closer to
!!! the interface than this scalar are kept unchanged. Can be a vector
!!! or a scalar field. If vector, lda must be given.
@@ -216,8 +194,8 @@
REAL(MK), DIMENSION(:,:,:,:,:) , POINTER :: udata
!!! Field of the function to be extended, defined in the same narrow
!!! band as the level function. Values outside this band are set to
- !!! HUGE. This has to be allocated to proper size incl. ghost layers of
- !!! size order! If ivalue is a function pointer, udata will be
+ !!! HUGE. This has to be allocated to proper size incl. ghost layers
+ !!! of size order! If ivalue is a function pointer, udata will be
!!! completely replaced. If ivalue is a scalar, the points closer to
!!! the interface than this scalar are kept unchanged. Can be a vector
!!! or a scalar field. If vector, lda must be given.
@@ -231,30 +209,29 @@
!!! to be the interface. Always a scalar field. THIS NEEDS TO BE
!!! PROPERLY ALLOCATED ON INPUT, INCLUDING GHOST LAYERS OF SIZE order.
REAL(MK), DIMENSION(:,:,:,:,:) , INTENT(IN),OPTIONAL:: chi
- !!! rank 5 (3d) or rank 4 (2d) field specifying the positions
- !!! of the grid nodes. 1st index: 1..ppm_dim, then i,j,[k],isub.
- !!! OPTIONAL. Uniform grid is assumed if absent. Ghostlayers of size
- !!! >=1 must be pre-filled.
+ !!! rank 5 (3d) field specifying the positions of the
+ !!! grid nodes. 1st index: 1..ppm_dim, then i,j,[k],isub. OPTIONAL.
+ !!! Uniform grid is assumed if absent. Ghostlayers of size >=1 must
+ !!! be pre-filled.
#endif
INTEGER , INTENT(IN ) :: order
- !!! Order of the method to be used. One of
+ !!! Desired order of the method. One of:
!!!
!!! *ppm_param_order_1
!!! *ppm_param_order_2
!!! *ppm_param_order_3
REAL(MK) , INTENT(IN ) :: tol
- !!! Relative tolerance for the determined distance to the interface. 1E-3
- !!! is a good choice. The tolerance is in multiples of grid spacings.
+ !!! Relative tolerance for the determined distance to the interface.
+ !!! 1E-3 is a good choice. The tolerance is in multiples of grid spacings
REAL(MK) , INTENT(IN ) :: width
- !!! Width of the narrow band to be produced on each side of
- !!! the interface.
+ !!! Width of the narrow band to be produced on each side of the interface.
INTEGER , INTENT( OUT) :: info
- !!! Return status, 0 upon success
+ !!! Return status. 0 upon success
INTEGER , OPTIONAL , INTENT(IN ) :: MaxIter
!!! OPTIONAL argument specifying the maximum number of allowed
- !!! iterations. This can be useful since a cyclic dependency in the GMM
- !!! algorithms could cause infinite loops. In each iteration at least one
- !!! point is computed.
+ !!! iterations. This can be useful since a cyclic dependency in the
+ !!! GMM algorithms could cause infinite loops. In each iteration at least
+ !!! one point is computed.
!-------------------------------------------------------------------------
! Local variables
!-------------------------------------------------------------------------
@@ -265,8 +242,6 @@
REAL(MK) :: t0,x,y,z,big
LOGICAL :: lok
REAL(MK), DIMENSION(:,:), POINTER :: closest
- TYPE(ppm_t_topo), POINTER :: topo
- TYPE(ppm_t_equi_mesh), POINTER :: mesh
#if __TYPE == __VFIELD
#if __DIM == __2D
REAL(MK), DIMENSION(:,:,: ), POINTER :: ext_wrk
@@ -279,8 +254,6 @@
!-------------------------------------------------------------------------
CALL substart('ppm_gmm_extend',t0,info)
big = HUGE(big)
- topo => ppm_topo(gmm_topoid)%t
- mesh => topo%mesh(gmm_meshid)
!-------------------------------------------------------------------------
! Set pointers
!-------------------------------------------------------------------------
@@ -362,7 +335,7 @@
! ldu(1) = maxxhi + ghostsize(1)
! ldu(2) = maxyhi + ghostsize(2)
! ldu(3) = maxzhi + ghostsize(3)
-! ldu(4) = topo%nsublist
+! ldu(4) = ppm_nsublist(gmm_topoid)
! CALL ppm_alloc(udata,ldl,ldu,iopt,info)
! IF (info .NE. ppm_param_success) THEN
! info = ppm_error_fatal
@@ -443,7 +416,7 @@
! ldu(1) = maxxhi + ghostsize(1)
! ldu(2) = maxyhi + ghostsize(2)
! ldu(3) = maxzhi + ghostsize(3)
-! ldu(4) = topo%nsublist
+! ldu(4) = ppm_nsublist(gmm_topoid)
! CALL ppm_alloc(udata,ldu,iopt,info)
! IF (info .NE. ppm_param_success) THEN
! info = ppm_error_fatal
@@ -455,11 +428,11 @@
! Nuke points farther from the interface than ivalue
!-------------------------------------------------------------------------
#if __DIM == __3D
- DO isub=1,topo%nsublist
- jsub = topo%isublist(isub)
- DO k=1,mesh%nnodes(3,jsub)
- DO j=1,mesh%nnodes(2,jsub)
- DO i=1,mesh%nnodes(1,jsub)
+ DO isub=1,ppm_nsublist(gmm_topoid)
+ jsub = ppm_isublist(isub,gmm_topoid)
+ DO k=1,ppm_cart_mesh(gmm_meshid,gmm_topoid)%nnodes(3,jsub)
+ DO j=1,ppm_cart_mesh(gmm_meshid,gmm_topoid)%nnodes(2,jsub)
+ DO i=1,ppm_cart_mesh(gmm_meshid,gmm_topoid)%nnodes(1,jsub)
IF (ABS(fdata(i,j,k,isub)) .GT. ivalue) THEN
#if __TYPE == __VFIELD
DO ida=1,lda
@@ -474,10 +447,10 @@
ENDDO
ENDDO
#elif __DIM == __2D
- DO isub=1,topo%nsublist
- jsub = topo%isublist(isub)
- DO j=1,mesh%nnodes(2,jsub)
- DO i=1,mesh%nnodes(1,jsub)
+ DO isub=1,ppm_nsublist(gmm_topoid)
+ jsub = ppm_isublist(isub,gmm_topoid)
+ DO j=1,ppm_cart_mesh(gmm_meshid,gmm_topoid)%nnodes(2,jsub)
+ DO i=1,ppm_cart_mesh(gmm_meshid,gmm_topoid)%nnodes(1,jsub)
IF (ABS(fdata(i,j,isub)) .GT. ivalue) THEN
#if __TYPE == __VFIELD
DO ida=1,lda
@@ -519,7 +492,7 @@
#if __DIM == __3D
ldu(3) = UBOUND(udata,4)
#endif
- ldu(4) = topo%nsublist
+ ldu(4) = ppm_nsublist(gmm_topoid)
CALL ppm_alloc(ext_wrk,ldl,ldu,iopt,info)
IF (info .NE. ppm_param_success) THEN
info = ppm_error_fatal
diff --git a/src/ppm_gmm_extend_fwd.f b/src/ppm_gmm_extend_fwd.f
index 9ec0a24..6c1dba5 100644
--- a/src/ppm_gmm_extend_fwd.f
+++ b/src/ppm_gmm_extend_fwd.f
@@ -22,7 +22,7 @@
& rhscst,dxinv,dyinv,dzinv,ghostsize,info,speed,chi)
#endif
#endif
- !!! This routine performs the forward marching step of the GMM. See
+ !!! This routine performs the forward extension step of the GMM. See
!!! ppm_gmm_march for details.
!!!
!!! === References ===
@@ -34,11 +34,6 @@
USE ppm_module_data
USE ppm_module_data_mesh
USE ppm_module_data_gmm
- USE ppm_module_substart
- USE ppm_module_substop
- USE ppm_module_error
- USE ppm_module_alloc
- USE ppm_module_typedef
IMPLICIT NONE
#if __KIND == __SINGLE_PRECISION
INTEGER, PARAMETER :: MK = ppm_kind_single
@@ -57,37 +52,33 @@
!-------------------------------------------------------------------------
#if __DIM == __2D
REAL(MK), DIMENSION(:,:,:) , POINTER :: fdta
-#elif __DIM == __3D
- REAL(MK), DIMENSION(:,:,:,:) , POINTER :: fdta
-#endif
!!! pointer to level function. Needs to be defined in a band
- !!! (width+order*dx).
-#if __DIM == __2D
+ !!! (width+order*dx)
REAL(MK), DIMENSION(:,:,:) , POINTER :: dta
-#elif __DIM == __3D
- REAL(MK), DIMENSION(:,:,:,:) , POINTER :: dta
-#endif
!!! pointer to value function.
-#if __DIM == __2D
REAL(MK), DIMENSION(:,:,:) , INTENT(IN), OPTIONAL :: speed
-#elif __DIM == __3D
- REAL(MK), DIMENSION(:,:,:,:) , INTENT(IN), OPTIONAL :: speed
-#endif
- !!! rank 4 (3d) or rank 3 (2d) field of front speeds.
- !!! OPTIONAL to override rhscst.
-#if __DIM == __2D
+ !!! rank 3 (2d) field of front speeds. OPTIONAL to override rhscst.
REAL(MK), DIMENSION(:,:,:,:) , INTENT(IN), OPTIONAL :: chi
+ !!! rank 4 (2d) field specifying the positions of the grid nodes.
+ !!! 1st index: 1..ppm_dim, then i,j,[k],isub. OPTIONAL. Uniform grid is
+ !!! assumed if absent.
#elif __DIM == __3D
+ REAL(MK), DIMENSION(:,:,:,:) , POINTER :: fdta
+ !!! pointer to level function. Needs to be defined in a band
+ !!! (width+order*dx)
+ REAL(MK), DIMENSION(:,:,:,:) , POINTER :: dta
+ !!! pointer to value function.
+ REAL(MK), DIMENSION(:,:,:,:) , INTENT(IN), OPTIONAL :: speed
+ !!! rank 4 (3d) field of front speeds. OPTIONAL to override rhscst.
REAL(MK), DIMENSION(:,:,:,:,:) , INTENT(IN), OPTIONAL :: chi
+ !!! rank 5 (3d) field specifying the positions of the grid nodes.
+ !!! 1st index: 1..ppm_dim, then i,j,[k],isub. OPTIONAL. Uniform grid is
+ !!! assumed if absent.
#endif
- !!! rank 5 (3d) or rank 4 (2d) field specifying the positions
- !!! of the grid nodes. 1st index: 1..ppm_dim, then i,j,[k],isub.
- !!! OPTIONAL. Uniform grid is assumed if absent.
REAL(MK) , INTENT(IN ) :: width
- !!! Width of the narrow band to be produced on each side of
- !!! the interface.
+ !!! Width of the narrow band to be produced on each side of the interface
REAL(MK) , INTENT(IN ) :: rhscst
- !!! constant value for the right hand side of grad u * grad f = c.
+ !!! constant value for the right hand side of grad u * grad f = c.
!!! If speed is present, this argument will be ignored.
REAL(MK) , INTENT(IN ) :: TM
!!! Current threshold for wave front location.
@@ -97,18 +88,18 @@
!!! inverse of the y grid spacing.
REAL(MK) , INTENT(IN ) :: dzinv
!!! inverse of the z grid spacing. (Not used in 2D version).
+ INTEGER, DIMENSION(3) , INTENT(IN ) :: ghostsize
+ !!! Size of the ghostlayer on all sides.
INTEGER , INTENT(IN ) :: order
!!! Order of the method to be used. One of
!!!
!!! *ppm_param_order_1
!!! *ppm_param_order_2
!!! *ppm_param_order_3
- INTEGER, DIMENSION(3) , INTENT(IN ) :: ghostsize
- !!! Size of the ghostlayer on all sides.
INTEGER , INTENT(INOUT) :: npos
!!! Current number of points in the close set.
INTEGER , INTENT( OUT) :: info
- !!! Return status, 0 upon success
+ !!! Return status. 0 upon success
!-------------------------------------------------------------------------
! Local variables
!-------------------------------------------------------------------------
@@ -128,8 +119,6 @@
REAL(MK), DIMENSION(-3:3,ppm_dim):: phi,psi
REAL(MK), DIMENSION(ppm_dim) :: alpha,beta
REAL(MK), DIMENSION(2) :: roots
- TYPE(ppm_t_topo), POINTER :: topo
- TYPE(ppm_t_equi_mesh), POINTER :: mesh
!-------------------------------------------------------------------------
! Externals
!-------------------------------------------------------------------------
@@ -138,8 +127,6 @@
! Initialise
!-------------------------------------------------------------------------
CALL substart('ppm_gmm_extend_fwd',t0,info)
- topo => ppm_topo(gmm_topoid)%t
- mesh => topo%mesh(gmm_meshid)
phi = 0.0_MK
psi = 0.0_MK
big = HUGE(big)
@@ -189,10 +176,10 @@
jj = gmm_ipos(2,p)
kk = gmm_ipos(3,p)
jsub = gmm_ipos(4,p)
- isub = topo%isublist(jsub)
- xhi = mesh%nnodes(1,isub)
- yhi = mesh%nnodes(2,isub)
- zhi = mesh%nnodes(3,isub)
+ isub = ppm_isublist(jsub,gmm_topoid)
+ xhi = ppm_cart_mesh(gmm_meshid,gmm_topoid)%nnodes(1,isub)
+ yhi = ppm_cart_mesh(gmm_meshid,gmm_topoid)%nnodes(2,isub)
+ zhi = ppm_cart_mesh(gmm_meshid,gmm_topoid)%nnodes(3,isub)
fdta0= fdta(ii,jj,kk,jsub)
absfdta0 = fdta0
IF (absfdta0 .LT. 0.0_MK) absfdta0 = -absfdta0
@@ -442,9 +429,9 @@
ii = gmm_ipos(1,p)
jj = gmm_ipos(2,p)
jsub = gmm_ipos(3,p)
- isub = topo%isublist(jsub)
- xhi = mesh%nnodes(1,isub)
- yhi = mesh%nnodes(2,isub)
+ isub = ppm_isublist(jsub,gmm_topoid)
+ xhi = ppm_cart_mesh(gmm_meshid,gmm_topoid)%nnodes(1,isub)
+ yhi = ppm_cart_mesh(gmm_meshid,gmm_topoid)%nnodes(2,isub)
fdta0= fdta(ii,jj,jsub)
!---------------------------------------------------------------------
! GMM update condition (see Kim:2001a)
diff --git a/src/ppm_gmm_finalize.f b/src/ppm_gmm_finalize.f
index cbb5893..b5c0ff9 100644
--- a/src/ppm_gmm_finalize.f
+++ b/src/ppm_gmm_finalize.f
@@ -7,8 +7,8 @@
!-------------------------------------------------------------------------
SUBROUTINE ppm_gmm_finalize(info)
- !!! This routine finalizes the ppm_gmm module and deallocates all data
- !!! structures.
+ !!! This routine finalizes the ppm_gmm module and deallocates all
+ !!! data structures.
!-------------------------------------------------------------------------
! Includes
!-------------------------------------------------------------------------
@@ -18,16 +18,12 @@
!-------------------------------------------------------------------------
USE ppm_module_data
USE ppm_module_data_gmm
- USE ppm_module_substart
- USE ppm_module_substop
- USE ppm_module_error
- USE ppm_module_alloc
IMPLICIT NONE
!-------------------------------------------------------------------------
! Arguments
!-------------------------------------------------------------------------
INTEGER, INTENT( OUT) :: info
- !!! Return status, 0 upon success
+ !!! Return status. 0 upon success
!-------------------------------------------------------------------------
! Local variables
!-------------------------------------------------------------------------
diff --git a/src/ppm_gmm_init.f b/src/ppm_gmm_init.f
index 5a6edc0..d0e169a 100644
--- a/src/ppm_gmm_init.f
+++ b/src/ppm_gmm_init.f
@@ -5,8 +5,8 @@
! ETH Zurich
! CH-8092 Zurich, Switzerland
!-------------------------------------------------------------------------
- SUBROUTINE ppm_gmm_init(field_topoid,meshid,Nest,prec,info)
- !!! This routine initializes the ppm_gmm module and allocates all data
+ SUBROUTINE ppm_gmm_init(meshid,Nest,prec,info)
+ !!! This routine initializes the ppm_gmm module and allocates all data
!!! structures.
!-------------------------------------------------------------------------
! Includes
@@ -18,8 +18,6 @@
USE ppm_module_data
USE ppm_module_data_gmm
USE ppm_module_data_mesh
- USE ppm_module_error
- USE ppm_module_typedef
IMPLICIT NONE
!-------------------------------------------------------------------------
! Arguments
@@ -35,19 +33,15 @@
!!! *ppm_kind_double
INTEGER, INTENT(IN ) :: meshid
!!! Mesh ID (user numbering) for which a GMM should be initialized.
- INTEGER, INTENT(IN ) :: field_topoid
- !!! Topo ID of the field
INTEGER, INTENT( OUT) :: info
- !!! Returns status, 0 upon success
+ !!! Return status. 0 upon success
!-------------------------------------------------------------------------
! Local variables
!-------------------------------------------------------------------------
- INTEGER, DIMENSION(3) :: ldu
- INTEGER :: iopt,i,isub
- LOGICAL :: lok
- REAL(ppm_kind_double) :: t0
- TYPE(ppm_t_topo), POINTER :: topo
- TYPE(ppm_t_equi_mesh), POINTER :: mesh
+ INTEGER, DIMENSION(3) :: ldu
+ INTEGER :: iopt,i,isub
+ LOGICAL :: lok
+ REAL(ppm_kind_double) :: t0
!-------------------------------------------------------------------------
! Externals
!-------------------------------------------------------------------------
@@ -56,8 +50,6 @@
! Initialise
!-------------------------------------------------------------------------
CALL substart('ppm_gmm_init',t0,info)
- topo => ppm_topo(field_topoid)%t
- mesh => topo%mesh(gmm_meshid)
!-------------------------------------------------------------------------
! Check arguments
!-------------------------------------------------------------------------
@@ -68,7 +60,8 @@
& 'Please call ppm_init first!',__LINE__,info)
GOTO 9999
ENDIF
- CALL ppm_check_meshid(field_topoid,meshid,lok,info)
+ CALL ppm_check_meshid(ppm_param_id_user,meshid,ppm_field_topoid, &
+ & lok,info)
IF (.NOT. lok) THEN
info = ppm_error_error
CALL ppm_error(ppm_err_argument,'ppm_gmm_init', &
@@ -122,27 +115,29 @@
& 'sparse data locations GMM_IPOS',__LINE__,info)
GOTO 9999
ENDIF
-
+ !-------------------------------------------------------------------------
+ ! Set topoid for all GMM operations
+ !-------------------------------------------------------------------------
+ gmm_topoid = ppm_field_topoid
!-------------------------------------------------------------------------
! Translate meshid to internal numbering and store it
!-------------------------------------------------------------------------
-
- gmm_meshid = mesh%ID
+ gmm_meshid = ppm_meshid(gmm_topoid)%internal(meshid)
!-------------------------------------------------------------------------
! Determine max extent of mesh in any sub
!-------------------------------------------------------------------------
maxxhi = 0
maxyhi = 0
maxzhi = 0
- DO i=1,topo%nsublist
- isub = topo%isublist(i)
- IF (mesh%nnodes(1,isub).GT.maxxhi) &
- & maxxhi = mesh%nnodes(1,isub)
- IF (mesh%nnodes(2,isub).GT.maxyhi) &
- & maxyhi = mesh%nnodes(2,isub)
+ DO i=1,ppm_nsublist(gmm_topoid)
+ isub = ppm_isublist(i,gmm_topoid)
+ IF (ppm_cart_mesh(gmm_meshid,gmm_topoid)%nnodes(1,isub).GT.maxxhi) &
+ & maxxhi = ppm_cart_mesh(gmm_meshid,gmm_topoid)%nnodes(1,isub)
+ IF (ppm_cart_mesh(gmm_meshid,gmm_topoid)%nnodes(2,isub).GT.maxyhi) &
+ & maxyhi = ppm_cart_mesh(gmm_meshid,gmm_topoid)%nnodes(2,isub)
IF (ppm_dim .GT. 2) THEN
- IF (mesh%nnodes(3,isub).GT.maxzhi)&
- & maxzhi = mesh%nnodes(3,isub)
+ IF (ppm_cart_mesh(gmm_meshid,gmm_topoid)%nnodes(3,isub).GT.maxzhi)&
+ & maxzhi = ppm_cart_mesh(gmm_meshid,gmm_topoid)%nnodes(3,isub)
ENDIF
ENDDO
!-------------------------------------------------------------------------
diff --git a/src/ppm_gmm_kickoff.f b/src/ppm_gmm_kickoff.f
index b278312..3e91839 100644
--- a/src/ppm_gmm_kickoff.f
+++ b/src/ppm_gmm_kickoff.f
@@ -29,18 +29,23 @@
!!! level set, apply shifting as pre-processing.
!!!
!!! [NOTE]
+ !!! ======================================================================
!!! This routine creates second-order accurate data near the interface
!!! and is needed to kick off a reinitialize a level set or extend a
!!! function to a narrow band.
+ !!!
!!! The order of the method is limited by the order of the finite
!!! differences used in computing the rhs of the interpolation system.
!!! Use higher order FD to get higher order initialization.
+ !!!
!!! Tests have shown that storing shifted indices (i.e. ip1 = i+1) is
!!! faster than using i+1 in the array index directly. We thus use
!!! this technique here.
+ !!!
!!! Maybe we shoud actually allocate ipos and copy the stuff so it will
!!! survive ppm_gmm_finalize?? This will be easy: just change the pointer
!!! assignment at the end of the routine to a physical copy operation.
+ !!! ======================================================================
!!!
!!! === References ====
!!!
diff --git a/src/ppm_gmm_reinitialize.f b/src/ppm_gmm_reinitialize.f
index cc7e7bf..0da3ba9 100644
--- a/src/ppm_gmm_reinitialize.f
+++ b/src/ppm_gmm_reinitialize.f
@@ -1,111 +1,8 @@
!-------------------------------------------------------------------------
! Subroutine : ppm_gmm_reinitialize
!-------------------------------------------------------------------------
- !
- ! Purpose : This routine re-initializes a signed distance
- ! level function with the zero level representing
- ! the interface. Shift accordingly if other level
- ! is to be used. The group marching method is
- ! used.
- !
- ! Input : tol (F) Relative tolerance for the
- ! determined distance to the
- ! interface. 1E-3 is a good
- ! choice. The tolerance is in
- ! multiples of grid spacings.
- ! width (F) Width of the narrow band to
- ! be produced on each side of
- ! the interface.
- ! order (I) Desired order of the method.
- ! One of:
- ! ppm_param_order_1
- ! ppm_param_order_2
- ! ppm_param_order_3
- ! thresh (F) OPTIONAL. Threshold for
- ! interface detection. If this is
- ! not specified, it is set to
- ! MAXVAL(ABS(fdata)).
- ! chi([:],:,:,:,:)(F) rank 5 (3d) or rank 4 (2d)
- ! field specifying the positions
- ! of the grid nodes. 1st index:
- ! 1..ppm_dim, then i,j,[k],isub.
- ! OPTIONAL. Uniform grid is
- ! assumed if absent. Ghostlayers
- ! of size >=1 must be pre-filled.
- ! MaxIter (I) OPTIONAL argument specifying the
- ! maximum number of allowed
- ! iterations. This can be useful
- ! since a cyclic dependency in the
- ! GMM algorithms could cause
- ! infinite loops. In each iteration
- ! at least one point is computed.
- !
- ! Input/output : fdata(...) (F) Field data. Either rank 3
- ! (for 2D scalar fields), or rank
- ! 4 (for 3D scalar fields).
- ! Indices: (i,j,[k],isub).
- ! On input: old level function
- ! values. The interface is at
- ! level zero.
- ! A ghostsize of 1 is needed
- ! on all sides which must be
- ! filled with the old level
- ! function value on input!!
- ! On output: reinitialized
- ! signed distance function
- ! using the interpolation
- ! method of Chopp. Points far
- ! from the interface will have
- ! the value HUGE.
- !
- ! Output : info (I) return status. 0 on success.
- !
- ! Remarks :
- !
- ! References : S. Kim. An O(N) level set method for Eikonal equations.
- ! SIAM J. Sci. Comput. 22(6):2178-2193, 2001.
- !
- ! Revisions :
- !-------------------------------------------------------------------------
- ! $Log: ppm_gmm_reinitialize.f,v $
- ! Revision 1.1.1.1 2007/07/13 10:18:55 ivos
- ! CBL version of the PPM library
- !
- ! Revision 1.9 2006/07/03 12:57:36 ivos
- ! Added literature references to header comments.
- !
- ! Revision 1.8 2006/04/06 14:40:29 ivos
- ! Added the MaxIter argument to specify the maximum number of allowed
- ! iterations for the GMM marching.
- !
- ! Revision 1.7 2005/07/14 19:58:15 ivos
- ! Added OPTIONAL argument chi for mesh node positions in distorted
- ! (mapped) meshes. For use with AGM for example.
- !
- ! Revision 1.6 2005/06/01 05:17:48 ivos
- ! Added optional argument thresh for kickoff threshold.
- !
- ! Revision 1.5 2005/04/27 01:06:13 ivos
- ! Convergence tests completed, cleaned up code, optmized code (Shark),
- ! and changed structure to allow faster compilation.
- !
- ! Revision 1.4 2005/04/21 04:49:57 ivos
- ! bugfix: pointers to array slabs were incorrectly used.
- !
- ! Revision 1.3 2005/03/16 06:20:09 ivos
- ! Several bugfixes. 1st order version is now tested. Moved all large
- ! data to the module.
- !
- ! Revision 1.2 2005/03/12 04:08:35 ivos
- ! Misc bug fixes.
- !
- ! Revision 1.1 2005/03/11 04:15:58 ivos
- ! Initial implementation.
- !
- !-------------------------------------------------------------------------
! Parallel Particle Mesh Library (PPM)
- ! Institute of Computational Science
- ! ETH Zentrum, Hirschengraben 84
+ ! ETH Zurich
! CH-8092 Zurich, Switzerland
!-------------------------------------------------------------------------
#if __DIM == __2D
@@ -125,6 +22,14 @@
& order,info,thresh,chi,MaxIter)
#endif
#endif
+ !!! This routine re-initializes a signed distance level function with the
+ !!! zero level representing the interface. Shift accordingly if other
+ !!! level is to be used. The group marching method is used.
+ !!!
+ !!! === References ===
+ !!!
+ !!! S. Kim. An O(N) level set method for Eikonal equations. SIAM J. Sci.
+ !!! Comput. 22(6):2178-2193, 2001.
!-------------------------------------------------------------------------
! Includes
!-------------------------------------------------------------------------
@@ -139,10 +44,6 @@
USE ppm_module_gmm_kickoff
USE ppm_module_gmm_march
USE ppm_module_gmm_finalize
- USE ppm_module_substart
- USE ppm_module_substop
- USE ppm_module_error
- USE ppm_module_typedef
IMPLICIT NONE
#if __KIND == __SINGLE_PRECISION | __KIND == __SINGLE_PRECISION_COMPLEX
INTEGER, PARAMETER :: MK = ppm_kind_single
@@ -154,28 +55,62 @@
!-------------------------------------------------------------------------
#if __DIM == __2D
REAL(MK), DIMENSION(:,:,: ) , POINTER :: fdata
+ !!! Field data. Rank 3 (for 2D scalar fields).
+ !!! Indices: (i,j,[k],isub). On input: old level function values. The
+ !!! interface is at level zero. A ghostsize of 1 is needed on all sides
+ !!! which must be filled with the old level function value on input!!
+ !!! On output: reinitialized signed distance function using the
+ !!! interpolation method of Chopp. Points far from the interface will have
+ !!! the value HUGE.
REAL(MK), DIMENSION(:,:,:,:) , INTENT(IN), OPTIONAL:: chi
+ !!! Rank 4 (2d) field specifying the positions of the grid nodes.
+ !!! 1st index: 1..ppm_dim, then i,j,[k],isub. OPTIONAL. Uniform grid is
+ !!! assumed if absent. Ghostlayers of size >=1 must be pre-filled.
#elif __DIM == __3D
REAL(MK), DIMENSION(:,:,:,:) , POINTER :: fdata
+ !!! Field data. Rank 4 (for 3D scalar fields).
+ !!! Indices: (i,j,[k],isub). On input: old level function values. The
+ !!! interface is at level zero. A ghostsize of 1 is needed on all sides
+ !!! which must be filled with the old level function value on input!!
+ !!! On output: reinitialized signed distance function using the
+ !!! interpolation method of Chopp. Points far from the interface will have
+ !!! the value HUGE.
REAL(MK), DIMENSION(:,:,:,:,:) , INTENT(IN), OPTIONAL:: chi
+ !!! Rank 5 (3d) field specifying the positions
+ !!! of the grid nodes. 1st index: 1..ppm_dim, then i,j,[k],isub.
+ !!! OPTIONAL. Uniform grid is assumed if absent. Ghostlayers of size >=1
+ !!! must be pre-filled.
#endif
INTEGER , INTENT(IN ) :: order
- REAL(MK) , INTENT(IN ) :: tol,width
+ !!! Order of the method to be used. One of
+ !!!
+ !!! *ppm_param_order_1
+ !!! *ppm_param_order_2
+ !!! *ppm_param_order_3
+ REAL(MK) , INTENT(IN ) :: tol
+ !!! Relative tolerance for the determined distance to the interface.
+ !!! 1E-3 is a good choice. The tolerance is in multiples of grid spacings.
+ REAL(MK) , INTENT(IN ) :: width
+ !!! Width of the narrow band to be produced on each side of the interface.
INTEGER , INTENT( OUT) :: info
REAL(MK) , OPTIONAL , INTENT(IN ) :: thresh
+ !!! OPTIONAL. Threshold for interface detection. If this is not specified,
+ !!! it is set to MAXVAL(ABS(fdata)).
INTEGER , OPTIONAL , INTENT(IN ) :: MaxIter
+ !!! OPTIONAL argument specifying the maximum number of allowed iterations.
+ !!! This can be useful since a cyclic dependency in the GMM algorithms
+ !!! could cause infinite loops. In each iteration at least one point is
+ !!! computed.
!-------------------------------------------------------------------------
! Local variables
!-------------------------------------------------------------------------
INTEGER :: xhi,i,isub,Nminit,MaxIt
REAL(MK) :: t0,th
LOGICAL :: lok
- TYPE(ppm_t_topo), POINTER :: topo
!-------------------------------------------------------------------------
! Initialise
!-------------------------------------------------------------------------
CALL substart('ppm_gmm_reinitialize',t0,info)
- topo => ppm_topo(gmm_topoid)%t
!-------------------------------------------------------------------------
! Check arguments
!-------------------------------------------------------------------------
@@ -199,7 +134,7 @@
GOTO 9999
ENDIF
#if __DIM == __3D
- IF (SIZE(fdata,4) .LT. topo%nsublist) THEN
+ IF (SIZE(fdata,4) .LT. ppm_nsublist(gmm_topoid)) THEN
info = ppm_error_error
CALL ppm_error(ppm_err_argument,'ppm_gmm_reinitialize', &
& 'field data for some subs is missing',__LINE__,info)
@@ -224,7 +159,7 @@
GOTO 9999
ENDIF
#elif __DIM == __2D
- IF (SIZE(fdata,3) .LT. topo%nsublist) THEN
+ IF (SIZE(fdata,3) .LT. ppm_nsublist(gmm_topoid)) THEN
info = ppm_error_error
CALL ppm_error(ppm_err_argument,'ppm_gmm_reinitialize', &
& 'field data for some subs is missing',__LINE__,info)
--
GitLab