!-------------------------------------------------------------------------
! Subroutine : ppm_mg_init
!-------------------------------------------------------------------------
!
! Purpose : This routine initializes the solver for
! 2D and 3D problems
!
! Input : equation (I) : KIND OF EQUATION TO BE SOLVED
! FOR THE MOMENT ONLY POISSON
! order (I) : ORDER OF FINITE DIFFERENCES
! NOW SECOND. THE GHOSTSIZE IS
! AUTOMATICALLY ADJUSTED
! smoother (I) : NOW GAUSS-SEIDEL
!
! [lda] (I) : LEADING DIMENSION, ONLY TO BE
! GIVEN FOR VECTOR CASES
!
! ibcdef (I) : ARRAY OF BOUNDARY CONDITION
!
!
! bcvalue (F) : ARRAY WHERE THE VALUES OF THE BC
! ARE STORED.IN CASE OF PERIODIC
! JUST GIVE ANY KIND OF VALUE
!
! EPSU (F) : STOPPING CRITERIUM. DETAIL:SHOULD
! BE SCALED WITH THE MAXIMUM VALUE ! OF THE RHS.
!
! limlev (I) :Number of levels that the user
! wants to coarse.
!
! wcycle (L) : TRUE if the user wants W-cycle.
! OTHERWISE FALSE
! lprint (L) : TRUE IF YOU WANT TO DUMP OUT
! INFORMATION
!
! omega (F) : relaxation parameter for SOR
!
!
! Input/output :
!
! Output : info (I) return status. 0 upon success.
!
! Remarks : PLEASE PAY ATTENTION THAT IN ORDER TO DIVIDE
! FURTHER A MESH IT SHOULD BE DIVISIBLE WITH 2.
! IF YOU WANT TO SOLVE DIFFERENT EQUATIONS
! THE WHOLE MACHINERY SHOULD BE CALLED TWICE.
! ALSO THE SOLVER IS NOW PROGRAMMED FOR THE POISSON
! PROBLEM. A FUTURE IMPROVEMENT WOULD BE
! TO USE A GENERAL STENCIL.
!
! References :
!
! Revisions :
!-------------------------------------------------------------------------
! $Log: ppm_mg_init.f,v $
! Revision 1.2 2006/08/22 15:54:37 pchatela
! Added a hopefully appropriate scaling factor in the comparisons against
! lmyeps
!
! Revision 1.1.1.1 2006/07/25 15:18:20 menahel
! initial import
!
! Revision 1.7 2006/02/03 09:34:03 ivos
! Fixed bug 00015: ppm_subs_bc was only allocated and stored for the
! local subs in topo_store. Several mapping routines however need the
! info about all (global) subs.
! Changed subs_bc to hold now the GLOBAL subid and adjusted all
! occurrences.
!
! Revision 1.6 2005/12/08 12:43:16 kotsalie
! commiting dirichlet
!
! Revision 1.5 2005/01/04 09:47:45 kotsalie
! ghostsize=2 for scalar case
!
! Revision 1.4 2004/10/29 15:59:14 kotsalie
! RED BLACK SOR FOR 3d vec case. 2d will soon follow.
!
! Revision 1.3 2004/09/28 14:04:49 kotsalie
! Changes concerning 4th order finite differences
!
! Revision 1.2 2004/09/23 09:38:30 kotsalie
! added details in the header
!
! Revision 1.1 2004/09/22 18:27:09 kotsalie
! MG new version
!
!-------------------------------------------------------------------------
! Parallel Particle Mesh Library (PPM)
! Institute of Computational Science
! ETH Zentrum, Hirschengraben 84
! CH-8092 Zurich, Switzerland
!-------------------------------------------------------------------------
#if __DIM == __SFIELD
#if __MESH_DIM == __2D
#if __KIND == __SINGLE_PRECISION
SUBROUTINE ppm_mg_init_2d_sca_s(topo_id,equation,iorder,smoother,ibcdef,&
& bcvalue,EPSU,mesh_id,limlev,wcycle,lprint,omega,info)
#elif __KIND == __DOUBLE_PRECISION
SUBROUTINE ppm_mg_init_2d_sca_d(topo_id,equation,iorder,smoother,ibcdef,&
& bcvalue,EPSU,mesh_id,limlev,wcycle,lprint,omega,info)
#endif
#elif __MESH_DIM == __3D
#if __KIND == __SINGLE_PRECISION
SUBROUTINE ppm_mg_init_3d_sca_s(topo_id,equation,iorder,smoother,ibcdef,&
& bcvalue,EPSU,mesh_id,limlev,wcycle,lprint,omega,info)
#elif __KIND == __DOUBLE_PRECISION
SUBROUTINE ppm_mg_init_3d_sca_d(topo_id,equation,iorder,smoother,ibcdef,&
& bcvalue,EPSU,mesh_id,limlev,wcycle,lprint,omega,info)
#endif
#endif
#elif __DIM == __VFIELD
#if __MESH_DIM == __2D
#if __KIND == __SINGLE_PRECISION
SUBROUTINE ppm_mg_init_2d_vec_s(topo_id,equation,iorder,smoother,lda,ibcdef,&
& bcvalue,EPSU,mesh_id,limlev,wcycle,lprint,omega,info)
#elif __KIND == __DOUBLE_PRECISION
SUBROUTINE ppm_mg_init_2d_vec_d(topo_id,equation,iorder,smoother,lda,ibcdef,&
& bcvalue,EPSU,mesh_id,limlev,wcycle,lprint,omega,info)
#endif
#elif __MESH_DIM == __3D
#if __KIND == __SINGLE_PRECISION
SUBROUTINE ppm_mg_init_3d_vec_s(topo_id,equation,iorder,smoother,lda,ibcdef,&
& bcvalue,EPSU,mesh_id,limlev,wcycle,lprint,omega,info)
#elif __KIND == __DOUBLE_PRECISION
SUBROUTINE ppm_mg_init_3d_vec_d(topo_id,equation,iorder,smoother,lda,ibcdef,&
& bcvalue,EPSU,mesh_id,limlev,wcycle,lprint,omega,info)
#endif
#endif
#endif
!----------------------------------------------------------------------
! Includes
!----------------------------------------------------------------------
#include "ppm_define.h"
!-----------------------------------------------------------------------
! Modules
!-----------------------------------------------------------------------
USE ppm_module_data
USE ppm_module_data_mesh
USE ppm_module_data_mg
USE ppm_module_alloc
USE ppm_module_mg_alloc
USE ppm_module_error
USE ppm_module_mesh_derive
USE ppm_module_substart
USE ppm_module_substop
USE ppm_module_typedef
IMPLICIT NONE
#if __KIND == __SINGLE_PRECISION
INTEGER, PARAMETER :: MK = ppm_kind_single
#else
INTEGER, PARAMETER :: MK = ppm_kind_double
#endif
!-----------------------------------------------------------------------
! Arguments
!-----------------------------------------------------------------------
INTEGER, INTENT(IN) :: equation
INTEGER, INTENT(IN) :: iorder
INTEGER, INTENT(IN) :: smoother
#if __DIM == __VFIELD
INTEGER, INTENT(IN) :: lda
#endif
#if __DIM == __SFIELD
#if __MESH_DIM == __2D
INTEGER,DIMENSION(:) :: ibcdef
REAL(MK),DIMENSION(:,:) :: bcvalue
#elif __MESH_DIM == __3D
INTEGER,DIMENSION(:) :: ibcdef
REAL(MK),DIMENSION(:,:,:) :: bcvalue
#endif
#elif __DIM == __VFIELD
#if __MESH_DIM == __2D
INTEGER,DIMENSION(:,:) :: ibcdef
REAL(MK),DIMENSION(:,:,:) :: bcvalue
#elif __MESH_DIM == __3D
INTEGER,DIMENSION(:,:) :: ibcdef
REAL(MK),DIMENSION(:,:,:,:) :: bcvalue
#endif
#endif
INTEGER, INTENT(IN) :: mesh_id,topo_id
REAL(MK),INTENT(IN) :: EPSU
INTEGER,INTENT(IN) :: limlev
LOGICAL,INTENT(IN) :: wcycle
LOGICAL,INTENT(IN) :: lprint
REAL(MK),INTENT(IN) :: omega
INTEGER, INTENT(OUT) :: info
!--------------------------------------------------------------------
! Local variables
!-----------------------------------------------------------------------
REAL(MK) :: t0
REAL(MK) :: lmyeps
INTEGER :: meshid,mlev,isub
INTEGER :: idom
INTEGER :: count,ilda,iface
INTEGER :: i,j,k
INTEGER :: kk
TYPE(ppm_t_topo), POINTER :: topo
TYPE(ppm_t_equi_mesh), POINTER :: mesh
#if __MESH_DIM == __2D
INTEGER :: dir
#endif
INTEGER :: iter1,iter2,ix,iy
INTEGER :: ipoint,jpoint
INTEGER :: newmeshid,lmesh_id
INTEGER , DIMENSION(1) :: ldu1
INTEGER , DIMENSION(2) :: ldu2,ldl2 ,direc
INTEGER , DIMENSION(3) :: ldu3,ldl3
#if __MESH_DIM == __3D
INTEGER :: dir1,dir2,jj,iz
INTEGER , DIMENSION(4) :: ldu4,ldl4
#endif
INTEGER , DIMENSION(ppm_dim) :: Nml
REAL(MK), DIMENSION(ppm_dim) :: min_phys,max_phys
REAL(MK), DIMENSION(ppm_dim,ppm_topo(topo_id)%t%nsubs) &
& :: min_sub,max_sub
INTEGER :: iopt,topoid
#if __DIM == __SFIELD
#if __MESH_DIM == __2D
#if __KIND == __SINGLE_PRECISION
TYPE(mg_field_2d_sca_s),DIMENSION(:,:),POINTER :: mgfield
#elif __KIND == __DOUBLE_PRECISION
TYPE(mg_field_2d_sca_d),DIMENSION(:,:),POINTER :: mgfield
#endif
#elif __MESH_DIM == __3D
#if __KIND == __SINGLE_PRECISION
TYPE(mg_field_3d_sca_s),DIMENSION(:,:),POINTER :: mgfield
#elif __KIND == __DOUBLE_PRECISION
TYPE(mg_field_3d_sca_d),DIMENSION(:,:),POINTER :: mgfield
#endif
#endif
#elif __DIM == __VFIELD
#if __MESH_DIM == __2D
#if __KIND == __SINGLE_PRECISION
TYPE(mg_field_2d_vec_s),DIMENSION(:,:),POINTER :: mgfield
#elif __KIND == __DOUBLE_PRECISION
TYPE(mg_field_2d_vec_d),DIMENSION(:,:),POINTER :: mgfield
#endif
#elif __MESH_DIM == __3D
#if __KIND == __SINGLE_PRECISION
TYPE(mg_field_3d_vec_s),DIMENSION(:,:),POINTER :: mgfield
#elif __KIND == __DOUBLE_PRECISION
TYPE(mg_field_3d_vec_d),DIMENSION(:,:),POINTER :: mgfield
#endif
#endif
#endif
#if __DIM == __SFIELD
#if __MESH_DIM == __2D
REAL(MK),DIMENSION(:,:),POINTER :: tuc
REAL(MK),DIMENSION(:,:),POINTER :: terr
#elif __MESH_DIM == __3D
REAL(MK),DIMENSION(:,:,:),POINTER :: tuc
REAL(MK),DIMENSION(:,:,:),POINTER :: terr
#endif
#elif __DIM == __VFIELD
#if __MESH_DIM == __2D
REAL(MK),DIMENSION(:,:,:),POINTER :: tuc
REAL(MK),DIMENSION(:,:,:),POINTER :: terr
#elif __MESH_DIM == __3D
REAL(MK),DIMENSION(:,:,:,:),POINTER :: tuc
REAL(MK),DIMENSION(:,:,:,:),POINTER :: terr
#endif
#endif
!-----------------------------------------------------------------------
! Externals
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
! Initialize
!-----------------------------------------------------------------------
CALL substart('ppm_mg_init',t0,info)
!-----------------------------------------------------------------------
! Check arguments
!-----------------------------------------------------------------------
IF (ppm_debug.GT.0) THEN
#if __DIM == __VFIELD
IF (lda.LE.0) THEN
info = ppm_error_error
CALL ppm_error(ppm_err_argument,'ppm_poiss_mg_init', &
& 'lda must be >0',__LINE__,info)
GOTO 9999
ENDIF
#endif
IF (EPSU.LE.0.0_MK) THEN
info = ppm_error_error
CALL ppm_error(ppm_err_argument,'ppm_poiss_mg_init', &
& 'EPSU must be >0',__LINE__,info)
GOTO 9999
ENDIF
ENDIF
!---------------------------------------------------------------------
! Definition of necessary variables and allocation of arrays
!---------------------------------------------------------------------
#if __DIM == __SFIELD
vecdim = 1
#elif __DIM == __VFIELD
vecdim = lda
#endif
w_cycle=wcycle
l_print=lprint
topoid = topo_id
topo => ppm_topo(topo_id)%t
mesh => topo%mesh(mesh_id)
nsubs = topo%nsublist
!PRINT *,'nsub:',nsubs
meshid = mesh%ID
lmesh_id = mesh_id
#if __KIND == __SINGLE_PRECISION
min_phys(:)=topo%min_physs(:)
max_phys(:)=topo%max_physs(:)
min_sub(:,:)=topo%min_subs(:,:)
max_sub(:,:)=topo%max_subs(:,:)
EPSU_s = EPSU
omega_s=omega
lmyeps=ppm_myepss
#elif __KIND == __DOUBLE_PRECISION
min_phys(:)=topo%min_physd(:)
max_phys(:)=topo%max_physd(:)
min_sub(:,:)=topo%min_subd(:,:)
max_sub(:,:)=topo%max_subd(:,:)
EPSU_d = EPSU
omega_d=omega
lmyeps=ppm_myepsd
#endif
#if __MESH_DIM == __2D
Nml(1) = mesh%Nm(1)
Nml(2) = mesh%Nm(2)
maxlev = INT(log10(Nml(1)*Nml(2)*REAL(ppm_nproc,MK))/log10(2.0_MK))
IF (maxlev.GT.limlev) THEN
maxlev=limlev
ENDIF
#if __KIND == __SINGLE_PRECISION
dx_s = (max_phys(1)-min_phys(1))/(Nml(1)-1)
dy_s = (max_phys(2)-min_phys(2))/(Nml(2)-1)
rdx2_s = 1/(dx_s*dx_s)
rdy2_s = 1/(dy_s*dy_s)
#elif __KIND == __DOUBLE_PRECISION
dx_d = (max_phys(1)-min_phys(1))/(Nml(1)-1)
dy_d = (max_phys(2)-min_phys(2))/(Nml(2)-1)
rdx2_d = 1/(dx_d*dx_d)
rdy2_d = 1/(dy_d*dy_d)
#endif
#elif __MESH_DIM == __3D
Nml(1) = mesh%Nm(1)
Nml(2) = mesh%Nm(2)
Nml(3) = mesh%Nm(3)
maxlev = INT(log10(Nml(1)*Nml(2)*Nml(3)* &
& REAL(ppm_nproc,MK))/log10(2.0_MK))
IF (maxlev.GT.limlev) THEN
maxlev=limlev
ENDIF
#if __KIND == __SINGLE_PRECISION
dx_s = (max_phys(1)-min_phys(1))/(Nml(1)-1)
dy_s = (max_phys(2)-min_phys(2))/(Nml(2)-1)
dz_s = (max_phys(3)-min_phys(3))/(Nml(3)-1)
rdx2_s = 1/(dx_s*dx_s)
rdy2_s = 1/(dy_s*dy_s)
rdz2_s = 1/(dz_s*dz_s)
#elif __KIND == __DOUBLE_PRECISION
dx_d = (max_phys(1)-min_phys(1))/(Nml(1)-1)
dy_d = (max_phys(2)-min_phys(2))/(Nml(2)-1)
dz_d = (max_phys(3)-min_phys(3))/(Nml(3)-1)
rdx2_d = 1/(dx_d*dx_d)
rdy2_d = 1/(dy_d*dy_d)
rdz2_d = 1/(dz_d*dz_d)
#endif
#endif
#if __DIM == __SFIELD
!Print *,'SField'
iopt = ppm_param_alloc_fit
ldu2(1) = nsubs
ldu2(2) = 2*ppm_dim
CALL ppm_alloc(bcdef_sca,ldu2,iopt,info)
IF (info .NE. 0) THEN
info = ppm_error_fatal
CALL ppm_error(ppm_err_alloc,'ppm_poiss_mg_init', &
& 'Boundary condiotions',__LINE__,info)
GOTO 9999
ENDIF
bcdef_sca(:,:)=0
!---------------------------------------------------------
!MICHAEL
!--------------------------------------------------------
DO isub=1,nsubs
idom=topo%isublist(isub)
!----------------------------------------------------------------------
! compare the west boundary
!----------------------------------------------------------------------
IF (ABS(min_sub(1,idom)-min_phys(1)) .LT. lmyeps*(max_sub(1,i)-min_sub(1,i))) THEN
bcdef_sca(isub,1)=ibcdef(1)
ENDIF
!----------------------------------------------------------------------
! compare the east boundary
!----------------------------------------------------------------------
IF (ABS(max_sub(1,idom)-max_phys(1)) .LT. lmyeps*(max_sub(1,i)-min_sub(1,i))) THEN
bcdef_sca(isub,2)=ibcdef(2)
ENDIF
!----------------------------------------------------------------------
! compare the south boundary
!----------------------------------------------------------------------
IF (ABS(min_sub(2,idom)-min_phys(2)) .LT. lmyeps*(max_sub(2,i)-min_sub(2,i))) THEN
bcdef_sca(isub,3)=ibcdef(3)
ENDIF
!----------------------------------------------------------------------
! compare the north boundary
!----------------------------------------------------------------------
IF (ABS(max_sub(2,idom)-max_phys(2)) .LT. lmyeps*(max_sub(2,i)-min_sub(2,i))) THEN
bcdef_sca(isub,4)=ibcdef(4)
ENDIF
!----------------------------------------------------------------------
! compare the south boundary
!----------------------------------------------------------------------
IF (ABS(min_sub(3,idom)-min_phys(3)) .LT. lmyeps*(max_sub(3,i)-min_sub(3,i))) THEN
bcdef_sca(isub,5)=ibcdef(5)
ENDIF
!----------------------------------------------------------------------
! compare the north boundary
!----------------------------------------------------------------------
IF (ABS(max_sub(3,idom)-max_phys(3)) .LT. lmyeps*(max_sub(3,i)-min_sub(3,i))) THEN
bcdef_sca(isub,6)=ibcdef(6)
ENDIF
ENDDO
lperiodic=.TRUE.
DO isub=1,nsubs
DO i=1,2*ppm_dim
IF (bcdef_sca(isub,i).NE.ppm_param_bcdef_periodic) THEN
lperiodic=.FALSE.
EXIT
ENDIF
ENDDO
ENDDO
#elif __DIM == __VFIELD
!print *,'Vfiedl'
iopt = ppm_param_alloc_fit
ldu3(1) = vecdim
ldu3(2) = nsubs
ldu3(3) = 2*ppm_dim
CALL ppm_alloc(bcdef_vec,ldu3,iopt,info)
IF (info .NE. 0) THEN
info = ppm_error_fatal
CALL ppm_error(ppm_err_alloc,'ppm_poiss_mg_init', &
& 'Boundary condiotions',__LINE__,info)
GOTO 9999
ENDIF
bcdef_vec(:,:,:)=0
Do isub=1,nsubs
idom=topo%isublist(isub)
Do ilda=1,vecdim
!----------------------------------------------------------------------
! compare the west boundary
!----------------------------------------------------------------------
IF (ABS(min_sub(1,idom)-min_phys(1)) .LT. lmyeps*(max_sub(1,i)-min_sub(1,i))) THEN
bcdef_vec(ilda,isub,1)=ibcdef(ilda,1)
ENDIF
!----------------------------------------------------------------------
! compare the east boundary
!----------------------------------------------------------------------
IF (ABS(max_sub(1,idom)-max_phys(1)) .LT. lmyeps*(max_sub(1,i)-min_sub(1,i))) THEN
bcdef_vec(ilda,isub,2)=ibcdef(ilda,2)
ENDIF
!----------------------------------------------------------------------
! compare the south boundary
!----------------------------------------------------------------------
IF (ABS(min_sub(2,idom)-min_phys(2)) .LT. lmyeps*(max_sub(2,i)-min_sub(2,i))) THEN
bcdef_vec(ilda,isub,3)=ibcdef(ilda,3)
ENDIF
!----------------------------------------------------------------------
! compare the north boundary
!----------------------------------------------------------------------
IF (ABS(max_sub(2,idom)-max_phys(2)) .LT. lmyeps*(max_sub(2,i)-min_sub(2,i))) THEN
bcdef_vec(ilda,isub,4)=ibcdef(ilda,4)
ENDIF
!----------------------------------------------------------------------
! compare the south boundary
!----------------------------------------------------------------------
IF (ABS(min_sub(3,idom)-min_phys(3)) .LT. lmyeps*(max_sub(3,i)-min_sub(3,i))) THEN
bcdef_vec(ilda,isub,5)=ibcdef(ilda,5)
ENDIF
!----------------------------------------------------------------------
! compare the north boundary
!----------------------------------------------------------------------
IF (ABS(max_sub(3,idom)-max_phys(3)) .LT. lmyeps*(max_sub(3,i)-min_sub(3,i))) THEN
bcdef_vec(ilda,isub,6)=ibcdef(ilda,6)
ENDIF
enddo
enddo
lperiodic=.TRUE.
Do isub=1,nsubs
DO i=1,2*ppm_dim
DO ilda=1,vecdim
IF (bcdef_vec(ilda,isub,i).NE.ppm_param_bcdef_periodic) Then
lperiodic=.FALSE.
EXIT
ENDIF
ENDDO
ENDDO
ENDDO
#endif
iopt = ppm_param_alloc_fit
ldu1(1) = ppm_dim
CALL ppm_alloc(ghostsize,ldu1,iopt,info)
IF (info .NE. 0) THEN
info = ppm_error_fatal
CALL ppm_error(ppm_err_alloc,'ppm_poiss_mg_init', &
& 'ghostsize',__LINE__,info)
GOTO 9999
ENDIF
IF (iorder.EQ.ppm_param_order_2) THEN
ghostsize(:)=1
order=iorder
ELSEIF (iorder.EQ.ppm_param_order_4) THEN
ghostsize(:)=2
order=ppm_param_order_4
ENDIF
iopt = ppm_param_alloc_fit
ldu1(1) = ppm_dim
CALL ppm_alloc(factor,ldu1,iopt,info)
IF (info .NE. 0) THEN
info = ppm_error_fatal
CALL ppm_error(ppm_err_alloc,'ppm_poiss_mg_init', &
& 'factor',__LINE__,info)
GOTO 9999
ENDIF
iopt = ppm_param_alloc_fit
ldu1(1) = maxlev
CALL ppm_alloc(meshid_g,ldu1,iopt,info)
IF (info .NE. 0) THEN
info = ppm_error_fatal
CALL ppm_error(ppm_err_alloc,'ppm_poiss_mg_init', &
& 'meshid_g',__LINE__,info)
GOTO 9999
ENDIF
iopt = ppm_param_alloc_fit
ldu1(1) = maxlev
CALL ppm_alloc(mesh_id_g,ldu1,iopt,info)
IF (info .NE. 0) THEN
info = ppm_error_fatal
CALL ppm_error(ppm_err_alloc,'ppm_poiss_mg_init', &
& 'mesh_id_g',__LINE__,info)
GOTO 9999
ENDIF
iopt = ppm_param_alloc_fit
ldu3(1) = ppm_dim
ldu3(2) = nsubs
ldu3(3) = maxlev
CALL ppm_alloc(start,ldu3,iopt,info)
IF (info .NE. 0) THEN
info = ppm_error_fatal
CALL ppm_error(ppm_err_alloc,'ppm_poiss_mg_init', &
& 'starting indices when updating the field',__LINE__,info)
GOTO 9999
ENDIF
iopt = ppm_param_alloc_fit
ldu3(1) = ppm_dim
ldu3(2) = nsubs
ldu3(3) = maxlev
CALL ppm_alloc(stop,ldu3,iopt,info)
IF (info .NE. 0) THEN
info = ppm_error_fatal
CALL ppm_error(ppm_err_alloc,'ppm_poiss_mg_init', &
& 'stopping indices when updating the field',__LINE__,info)
GOTO 9999
ENDIF
#if __DIM == __SFIELD
#if __MESH_DIM == __2D
#if __KIND == __SINGLE_PRECISION
iopt = ppm_param_alloc_fit
ldu2(1) = nsubs
ldu2(2) = maxlev
CALL ppm_mg_alloc(mgfield_2d_sca_s,ldu2,iopt,info)
IF (info .NE. 0) THEN
info = ppm_error_fatal
CALL ppm_error(ppm_err_alloc,'ppm_poiss_mg_init', &
& 'Multigrid fields used on the different levels',__LINE__,info)
GOTO 9999
ENDIF
mgfield => mgfield_2d_sca_s
#elif __KIND == __DOUBLE_PRECISION
iopt = ppm_param_alloc_fit
ldu2(1) = nsubs
ldu2(2) = maxlev
CALL ppm_mg_alloc(mgfield_2d_sca_d,ldu2,iopt,info)
IF (info .NE. 0) THEN
info = ppm_error_fatal
CALL ppm_error(ppm_err_alloc,'ppm_poiss_mg_init', &
& 'Multigrid fields used on the different levels',__LINE__,info)
GOTO 9999
ENDIF
mgfield => mgfield_2d_sca_d
#endif
#elif __MESH_DIM == __3D
#if __KIND == __SINGLE_PRECISION
iopt = ppm_param_alloc_fit
ldu2(1) = nsubs
ldu2(2) = maxlev
CALL ppm_mg_alloc(mgfield_3d_sca_s,ldu2,iopt,info)
IF (info .NE. 0) THEN
info = ppm_error_fatal
CALL ppm_error(ppm_err_alloc,'ppm_poiss_mg_init', &
& 'Multigrid fields used on the different levels',__LINE__,info)
GOTO 9999
ENDIF
mgfield => mgfield_3d_sca_s
#elif __KIND == __DOUBLE_PRECISION
iopt = ppm_param_alloc_fit
ldu2(1) = nsubs
ldu2(2) = maxlev
CALL ppm_mg_alloc(mgfield_3d_sca_d,ldu2,iopt,info)
IF (info .NE. 0) THEN
info = ppm_error_fatal
CALL ppm_error(ppm_err_alloc,'ppm_poiss_mg_init', &
& 'Multigrid fields used on the different levels',__LINE__,info)
GOTO 9999
ENDIF
mgfield => mgfield_3d_sca_d
#endif
#endif
#elif __DIM == __VFIELD
#if __MESH_DIM == __2D
#if __KIND == __SINGLE_PRECISION
iopt = ppm_param_alloc_fit
ldu2(1) = nsubs
ldu2(2) = maxlev
CALL ppm_mg_alloc(mgfield_2d_vec_s,ldu2,iopt,info)
IF (info .NE. 0) THEN
info = ppm_error_fatal
CALL ppm_error(ppm_err_alloc,'ppm_poiss_mg_init', &
& 'Multigrid fields used on the different levels',__LINE__,info)
GOTO 9999
ENDIF
mgfield => mgfield_2d_vec_s
#elif __KIND == __DOUBLE_PRECISION
iopt = ppm_param_alloc_fit
ldu2(1) = nsubs
ldu2(2) = maxlev
CALL ppm_mg_alloc(mgfield_2d_vec_d,ldu2,iopt,info)
IF (info .NE. 0) THEN
info = ppm_error_fatal
CALL ppm_error(ppm_err_alloc,'ppm_poiss_mg_init', &
& 'Multigrid fields used on the different levels',__LINE__,info)
GOTO 9999
ENDIF
mgfield => mgfield_2d_vec_d
#endif
#elif __MESH_DIM == __3D
#if __KIND == __SINGLE_PRECISION
iopt = ppm_param_alloc_fit
ldu2(1) = nsubs
ldu2(2) = maxlev
CALL ppm_mg_alloc(mgfield_3d_vec_s,ldu2,iopt,info)
IF (info .NE. 0) THEN
info = ppm_error_fatal
CALL ppm_error(ppm_err_alloc,'ppm_poiss_mg_init', &
& 'Multigrid fields used on the different levels',__LINE__,info)
GOTO 9999
ENDIF
mgfield => mgfield_3d_vec_s
#elif __KIND == __DOUBLE_PRECISION
iopt = ppm_param_alloc_fit
ldu2(1) = nsubs
ldu2(2) = maxlev
CALL ppm_mg_alloc(mgfield_3d_vec_d,ldu2,iopt,info)
IF (info .NE. 0) THEN
info = ppm_error_fatal
CALL ppm_error(ppm_err_alloc,'ppm_poiss_mg_init', &
& 'Multigrid fields used on the different levels',__LINE__,info)
GOTO 9999
ENDIF
mgfield => mgfield_3d_vec_d
#endif
#endif
#endif
iopt = ppm_param_alloc_fit
ldu2(1) = 2*ppm_dim
ldu2(2) = nsubs
CALL ppm_alloc(lboundary,ldu2,iopt,info)
IF (info .NE. 0) THEN
info = ppm_error_fatal
CALL ppm_error(ppm_err_alloc,'ppm_poiss_mg_init', &
& 'Problem with the boundary alloc.',__LINE__,info)
GOTO 9999
ENDIF
iopt = ppm_param_alloc_fit
ldu2(1) = ppm_dim
ldu2(2) = maxlev
CALL ppm_alloc(max_node,ldu2,iopt,info)
IF (info .NE. 0) THEN
info = ppm_error_fatal
CALL ppm_error(ppm_err_alloc,'ppm_poiss_mg_init', &
& 'Problem with a maximum number alloc.',__LINE__,info)
GOTO 9999
ENDIF
max_node(:,:)=0
lboundary(:,:)=.FALSE.
start(:,:,:)=1
!-----------------------------------------------------------------------
! Derive coarser meshes
!-----------------------------------------------------------------------
DO mlev=1,maxlev
#if __MESH_DIM == __2D
!--------------------------------------------------------------------
! Go through the subs, define the stopping indices on each mesh,
! check and store if it is on the boundary, allocate the
! multigrid fields, pass the boundary values.
!--------------------------------------------------------------------
DO i=1,nsubs
idom=topo%isublist(i)
stop(:,i,mlev)= mesh%nnodes(:,idom)
DO j=1,ppm_dim
IF (max_node(j,mlev).LT.stop(j,i,mlev)) THEN
max_node(j,mlev)=stop(j,i,mlev)
ENDIF
ENDDO
!-----------------------------------------------------------------
! Allocate the function correction, the restricted errors,
! the residuals and the values on the boundary on each level.
!----------------------------------------------------------------
#if __DIM == __SFIELD
iopt = ppm_param_alloc_fit
ldl2(1) = 1-ghostsize(1)
ldl2(2) = 1-ghostsize(2)
ldu2(1) = mesh%nnodes(1,idom)+ghostsize(1)
ldu2(2) = mesh%nnodes(2,idom)+ghostsize(2)
CALL ppm_alloc(mgfield(i,mlev)%uc,ldl2,ldu2,iopt,info)
IF (info .NE. 0) THEN
info = ppm_error_fatal
CALL ppm_error(ppm_err_alloc,'ppm_poiss_mg_init', &
& 'Problem with the function corr. alloc.',__LINE__,info)
GOTO 9999
ENDIF
tuc=>mgfield(i,mlev)%uc
tuc(:,:)=0.0_MK
iopt = ppm_param_alloc_fit
ldu2(1) = mesh%nnodes(1,idom)
ldu2(2) = mesh%nnodes(2,idom)
CALL ppm_alloc(mgfield(i,mlev)%fc,ldu2,iopt,info)
IF (info .NE. 0) THEN
info = ppm_error_fatal
CALL ppm_error(ppm_err_alloc,'ppm_poiss_mg_init', &
& 'Problem with the restricted err. alloc.',__LINE__,info)
GOTO 9999
ENDIF
mgfield(i,mlev)%fc(:,:)=0.0_MK
iopt = ppm_param_alloc_fit
ldl2(1) = 1-ghostsize(1)
ldl2(2) = 1-ghostsize(2)
ldu2(1) = mesh%nnodes(1,idom)+ghostsize(1)
ldu2(2) = mesh%nnodes(2,idom)+ghostsize(2)
CALL ppm_alloc(mgfield(i,mlev)%err,ldl2,ldu2,iopt,info)
IF (info .NE. 0) THEN
info = ppm_error_fatal
CALL ppm_error(ppm_err_alloc,'ppm_poiss_mg_init', &
& 'Problem with the residual alloc.',__LINE__,info)
GOTO 9999
ENDIF
terr=>mgfield(i,mlev)%err
terr(:,:)=0.0_MK
!------------------------------------------------------------------
!MICHAEL
!------------------------------------------------------------------
!ALLOCATE THE BCVALUE(IT IS A TYPE!!)
!PRINT *,'LPERIODIC:',lperiodic
IF (.NOT.lperiodic) THEN
iopt = ppm_param_alloc_fit
ldu1(1) = 2*ppm_dim
CALL ppm_mg_alloc(mgfield(i,mlev)%bcvalue,ldu1,iopt,info)
IF (info .NE. 0) THEN
info = ppm_error_fatal
CALL ppm_error(ppm_err_alloc,'ppm_poiss_mg_init', &
& 'Problem with the BOUNDARY alloc.',__LINE__,info)
GOTO 9999
ENDIF
!ALLOCATE THE PBCVALUE
DO iface=1,2*ppm_dim
iopt = ppm_param_alloc_fit
IF (iface.EQ.1.OR.iface.EQ.2) THEN
ldu1(1) = max_node(2,mlev)
ELSE
ldu1(1) = max_node(1,mlev)
ENDIF
CALL ppm_alloc(mgfield(i,mlev)%bcvalue(iface)%pbcvalue,ldu1,iopt,info)
IF (info .NE. 0) THEN
info = ppm_error_fatal
CALL ppm_error(ppm_err_alloc,'ppm_poiss_mg_init', &
& 'Problem with the BOUNDARY alloc.',__LINE__,info)
GOTO 9999
ENDIF
ENDDO
DO iface=1,2*ppm_dim
IF (iface.EQ.1.OR.iface.EQ.2) THEN
direc(1)=2
elseif (iface.EQ.3.OR.iface.EQ.4) then
direc(1)=1
endif
DO ipoint=1,max_node(direc(1),mlev)
IF (mlev.EQ.1) THEN
mgfield(i,mlev)%bcvalue(iface)%pbcvalue(ipoint)=bcvalue(iface,ipoint)
ELSE
IF(bcdef_sca(i,iface).EQ.ppm_param_bcdef_neumann) THEN
mgfield(i,mlev)%bcvalue(iface)%pbcvalue(ipoint)=&
& mgfield(i,mlev-1)%bcvalue(iface)%pbcvalue(2*ipoint-1)
ELSE
!NO CORRECTIONS FOR THE DIRICHLET
! mgfield(i,mlev)%bcvalue(iface)%pbcvalue(ipoint)=0.0_MK
mgfield(i,mlev)%bcvalue(iface)%pbcvalue(ipoint)=0.0_MK
ENDIF
ENDIF
ENDDO
\
ENDDO!faces
ENDIF!lperiodic
#elif __DIM == __VFIELD
iopt = ppm_param_alloc_fit
ldl3(1) = 1
ldl3(2) = 1-ghostsize(1)
ldl3(3) = 1-ghostsize(2)
ldu3(1) = vecdim
ldu3(2) = mesh%nnodes(1,idom)+ghostsize(1)
ldu3(3) = mesh%nnodes(2,idom)+ghostsize(2)
CALL ppm_alloc(mgfield(i,mlev)%uc,ldl3,ldu3,iopt,info)
IF (info .NE. 0) THEN
info = ppm_error_fatal
CALL ppm_error(ppm_err_alloc,'ppm_poiss_mg_init', &
& 'Problem with the function corr. alloc.',__LINE__,info)
GOTO 9999
ENDIF
tuc=>mgfield(i,mlev)%uc
tuc(:,:,:)=0.0_MK
iopt = ppm_param_alloc_fit
ldu3(1) = vecdim
ldu3(2) = mesh%nnodes(1,idom)
ldu3(3) = mesh%nnodes(2,idom)
CALL ppm_alloc(mgfield(i,mlev)%fc,ldu3,iopt,info)
IF (info .NE. 0) THEN
info = ppm_error_fatal
CALL ppm_error(ppm_err_alloc,'ppm_poiss_mg_init', &
& 'Problem with the restricted err. alloc.',__LINE__,info)
GOTO 9999
ENDIF
mgfield(i,mlev)%fc(:,:,:)=0.0_MK
iopt = ppm_param_alloc_fit
ldl3(1) = 1
ldl3(2) = 1-ghostsize(1)
ldl3(3) = 1-ghostsize(2)
ldu3(1) = vecdim
ldu3(2) = mesh%nnodes(1,idom)+ghostsize(1)
ldu3(3) = mesh%nnodes(2,idom)+ghostsize(2)
CALL ppm_alloc(mgfield(i,mlev)%err,ldl3,ldu3,iopt,info)
IF (info .NE. 0) THEN
info = ppm_error_fatal
CALL ppm_error(ppm_err_alloc,'ppm_poiss_mg_init', &
& 'Problem with the residual alloc.',__LINE__,info)
GOTO 9999
ENDIF
terr=>mgfield(i,mlev)%err
terr(:,:,:)=0.0_MK
#endif
iopt = ppm_param_alloc_fit
ldl2(1) = 1-ghostsize(1)
ldl2(2) = 1-ghostsize(2)
ldu2(1) = mesh%nnodes(1,idom)+ghostsize(1)
ldu2(2) = mesh%nnodes(2,idom)+ghostsize(2)
CALL ppm_alloc(mgfield(i,mlev)%mask_red,ldl2,ldu2,iopt,info)
IF (info .NE. 0) THEN
info = ppm_error_fatal
CALL ppm_error(ppm_err_alloc,'ppm_poiss_mg_init', &
& 'Problem with the mask alloc.',__LINE__,info)
GOTO 9999
ENDIF
iopt = ppm_param_alloc_fit
ldl2(1) = 1-ghostsize(1)
ldl2(2) = 1-ghostsize(2)
ldu2(1) = mesh%nnodes(1,idom)+ghostsize(1)
ldu2(2) = mesh%nnodes(2,idom)+ghostsize(2)
CALL ppm_alloc(mgfield(i,mlev)%mask_black,ldl2,ldu2,iopt,info)
IF (info .NE. 0) THEN
info = ppm_error_fatal
CALL ppm_error(ppm_err_alloc,'ppm_poiss_mg_init', &
& 'Problem with mask alloc.',__LINE__,info)
GOTO 9999
ENDIF
!----------------------------------------------------------------
!Filling the mask for communication (red black)
!----------------------------------------------------------------
DO iy=1-ghostsize(2),mesh%nnodes(2,idom)+ghostsize(2)
DO ix=1-ghostsize(1),mesh%nnodes(1,idom)+ghostsize(1)
IF (MOD(ix+iy,2).EQ.0) THEN
mgfield(i,mlev)%mask_red(ix,iy)=.TRUE.
mgfield(i,mlev)%mask_black(ix,iy)=.FALSE.
ELSE
mgfield(i,mlev)%mask_red(ix,iy) = .FALSE.
mgfield(i,mlev)%mask_black(ix,iy) = .TRUE.
ENDIF
ENDDO
ENDDO
ENDDO!DO 1,nsubs
#elif __MESH_DIM == __3D
DO i=1,nsubs
idom=topo%isublist(i)
stop(:,i,mlev) = mesh%nnodes(:,idom)
DO j=1,ppm_dim
IF (max_node(j,mlev).LT.stop(j,i,mlev)) THEN
max_node(j,mlev)=stop(j,i,mlev)
ENDIF
ENDDO
IF (topo%subs_bc(1,idom).EQ.1) THEN
lboundary(1,i)=.TRUE.
ELSEIF (topo%subs_bc(3,idom).EQ.1) THEN
lboundary(3,i)=.TRUE.
ELSEIF (topo%subs_bc(2,idom).EQ.1) THEN
lboundary(2,i)=.TRUE.
ELSEIF (topo%subs_bc(4,idom).EQ.1) THEN
lboundary(4,i)=.TRUE.
ELSEIF (topo%subs_bc(5,idom).EQ.1) THEN
lboundary(5,i)=.TRUE.
ELSEIF (topo%subs_bc(6,idom).EQ.1) THEN
lboundary(6,i)=.TRUE.
ENDIF
!----------------------------------------------------------------
! Allocate the function correction, the restricted errors and the
!residuals on each level.
!----------------------------------------------------------------
#if __DIM == __SFIELD
iopt = ppm_param_alloc_fit
ldl3(1) = 1-ghostsize(1)
ldl3(2) = 1-ghostsize(2)
ldl3(3) = 1-ghostsize(3)
ldu3(1) = mesh%nnodes(1,idom)+ghostsize(1)
ldu3(2) = mesh%nnodes(2,idom)+ghostsize(2)
ldu3(3) = mesh%nnodes(3,idom)+ghostsize(3)
CALL ppm_alloc(mgfield(i,mlev)%uc,ldl3,ldu3,iopt,info)
IF (info .NE. 0) THEN
info = ppm_error_fatal
CALL ppm_error(ppm_err_alloc,'ppm_poiss_mg_init', &
& 'Problem with the function corr. alloc.',__LINE__,info)
GOTO 9999
ENDIF
tuc=>mgfield(i,mlev)%uc
tuc(:,:,:)=0.0_MK
iopt = ppm_param_alloc_fit
ldu3(1) = mesh%nnodes(1,idom)
ldu3(2) = mesh%nnodes(2,idom)
ldu3(3) = mesh%nnodes(3,idom)
CALL ppm_alloc(mgfield(i,mlev)%fc,ldu3,iopt,info)
IF (info .NE. 0) THEN
info = ppm_error_fatal
CALL ppm_error(ppm_err_alloc,'ppm_poiss_mg_init', &
& 'Problem with the restricted err. alloc.',__LINE__,info)
GOTO 9999
ENDIF
mgfield(i,mlev)%fc(:,:,:)=0.0_MK
iopt = ppm_param_alloc_fit
ldl3(1) = 1-ghostsize(1)
ldl3(2) = 1-ghostsize(2)
ldl3(3) = 1-ghostsize(3)
ldu3(1) = mesh%nnodes(1,idom)+ghostsize(1)
ldu3(2) = mesh%nnodes(2,idom)+ghostsize(2)
ldu3(3) = mesh%nnodes(3,idom)+ghostsize(3)
CALL ppm_alloc(mgfield(i,mlev)%err,ldl3,ldu3,iopt,info)
IF (info .NE. 0) THEN
info = ppm_error_fatal
CALL ppm_error(ppm_err_alloc,'ppm_poiss_mg_init', &
& 'Problem with the residual alloc.',__LINE__,info)
GOTO 9999
ENDIF
terr=>mgfield(i,mlev)%err
terr(:,:,:)=0.0_MK
!ALLOCATE THE BCVALUE(IT IS A TYPE!!)
!PRINT *,'LPERIODIC:',lperiodic
IF (.NOT.lperiodic) THEN
iopt = ppm_param_alloc_fit
ldu1(1) = 2*ppm_dim
CALL ppm_mg_alloc(mgfield(i,mlev)%bcvalue,ldu1,iopt,info)
IF (info .NE. 0) THEN
info = ppm_error_fatal
CALL ppm_error(ppm_err_alloc,'ppm_poiss_mg_init', &
& 'Problem with the BOUNDARY alloc.',__LINE__,info)
GOTO 9999
ENDIF
!ALLOCATE THE PBCVALUE
DO iface=1,2*ppm_dim
iopt = ppm_param_alloc_fit
IF (iface.EQ.1.OR.iface.EQ.2) THEN
ldu2(1) = max_node(2,mlev)
ldu2(2)= max_node(3,mlev)
ELSEif (iface.EQ.3.OR. iface.EQ.4) then
ldu2(1) = max_node(1,mlev)
ldu2(2)=max_node(3,mlev)
else
ldu2(1)=max_node(1,mlev)
ldu2(2)=max_node(2,mlev)
ENDIF
CALL ppm_alloc(mgfield(i,mlev)%bcvalue(iface)%pbcvalue,ldu2,iopt,info)
!Print *,size(mgfield(i,mlev)%bcvalue(iface)%pbcvalue,1)
IF (info .NE. 0) THEN
info = ppm_error_fatal
CALL ppm_error(ppm_err_alloc,'ppm_poiss_mg_init', &
& 'Problem with the BOUNDARY alloc.',__LINE__,info)
GOTO 9999
ENDIF
ENDDO
DO iface=1,2*ppm_dim
IF (iface.EQ.1.OR.iface.EQ.2) THEN
direc(1)=2
direc(2)=3
elseif (iface.EQ.3.OR.iface.EQ.4) THEN
direc(1)=1
direc(2)=3
else
direc(1)=1
direc(2)=2
endif
DO ipoint=1,max_node(direc(1),mlev)
DO jpoint=1,max_node(direc(2),mlev)
IF (mlev.EQ.1) THEN
mgfield(i,mlev)%bcvalue(iface)%pbcvalue(ipoint,jpoint)=bcvalue(iface,ipoint,jpoint)
ELSE
IF(bcdef_sca(i,iface).EQ.ppm_param_bcdef_neumann) THEN
mgfield(i,mlev)%bcvalue(iface)%pbcvalue(ipoint,jpoint)=&
& mgfield(i,mlev-1)%bcvalue(iface)%pbcvalue(2*ipoint-1,2*jpoint-1)
ELSE
!NO CORRECTIONS FOR THE DIRICHLET
mgfield(i,mlev)%bcvalue(iface)%pbcvalue(ipoint,jpoint)=0.0_MK
! If (mlev.EQ.5) then
! Print *,ipoint,jpoint,mgfield(i,mlev)%bcvalue(iface)%pbcvalue(ipoint,jpoint)
! endif
ENDIF
ENDIF
ENDDO
enddo
ENDDO!faces
endif !lperiodic
#elif __DIM == __VFIELD
iopt = ppm_param_alloc_fit
ldl4(1) = 1
ldl4(2) = 1-ghostsize(1)
ldl4(3) = 1-ghostsize(2)
ldl4(4) = 1-ghostsize(3)
ldu4(1) = vecdim
ldu4(2) = mesh%nnodes(1,idom)+ghostsize(1)
ldu4(3) = mesh%nnodes(2,idom)+ghostsize(2)
ldu4(4) = mesh%nnodes(3,idom)+ghostsize(3)
CALL ppm_alloc(mgfield(i,mlev)%uc,ldl4,ldu4,iopt,info)
IF (info .NE. 0) THEN
info = ppm_error_fatal
CALL ppm_error(ppm_err_alloc,'ppm_poiss_mg_init', &
& 'Problem with the function corr. alloc.',__LINE__,info)
GOTO 9999
ENDIF
tuc=>mgfield(i,mlev)%uc
tuc(:,:,:,:)=0.0_MK
iopt = ppm_param_alloc_fit
ldu4(1) = vecdim
ldu4(2) = mesh%nnodes(1,idom)
ldu4(3) = mesh%nnodes(2,idom)
ldu4(4) = mesh%nnodes(3,idom)
CALL ppm_alloc(mgfield(i,mlev)%fc,ldu4,iopt,info)
IF (info .NE. 0) THEN
info = ppm_error_fatal
CALL ppm_error(ppm_err_alloc,'ppm_poiss_mg_init', &
& 'Problem with the restricted err. alloc.',__LINE__,info)
GOTO 9999
ENDIF
mgfield(i,mlev)%fc(:,:,:,:)=0.0_MK
iopt = ppm_param_alloc_fit
ldl4(1) = 1
ldl4(2) = 1-ghostsize(1)
ldl4(3) = 1-ghostsize(2)
ldl4(4) = 1-ghostsize(3)
ldu4(1) = vecdim
ldu4(2) = mesh%nnodes(1,idom)+ghostsize(1)
ldu4(3) = mesh%nnodes(2,idom)+ghostsize(2)
ldu4(4) = mesh%nnodes(3,idom)+ghostsize(3)
CALL ppm_alloc(mgfield(i,mlev)%err,ldl4,ldu4,iopt,info)
IF (info .NE. 0) THEN
info = ppm_error_fatal
CALL ppm_error(ppm_err_alloc,'ppm_poiss_mg_init', &
& 'Problem with the residual alloc.',__LINE__,info)
GOTO 9999
ENDIF
terr=>mgfield(i,mlev)%err
terr(:,:,:,:)=0.0_MK
!ALLOCATE THE BCVALUE(IT IS A TYPE!!)
!PRINT *,'LPERIODIC:',lperiodic
IF (.NOT.lperiodic) THEN
iopt = ppm_param_alloc_fit
ldu1=2*ppm_dim
!ldu2(1)=vecdim
!ldu2(2) = 2*ppm_dim
!allocate(mgfield(i,mlev)%bcvalue(3,6))
CALL ppm_mg_alloc(mgfield(i,mlev)%bcvalue,ldu1,iopt,info)
IF (info .NE. 0) THEN
info = ppm_error_fatal
CALL ppm_error(ppm_err_alloc,'ppm_poiss_mg_init', &
& 'Problem with the BOUNDARY alloc.',__LINE__,info)
GOTO 9999
ENDIF
!ALLOCATE THE PBCVALUE
DO iface=1,2*ppm_dim
iopt = ppm_param_alloc_fit
ldu3(1)=vecdim
IF (iface.EQ.1.OR.iface.EQ.2) THEN
ldu3(2) = max_node(2,mlev)
ldu3(3)= max_node(3,mlev)
ELSEif (iface.EQ.3.OR. iface.EQ.4) then
ldu3(2) = max_node(1,mlev)
ldu3(3)=max_node(3,mlev)
else
ldu3(2)=max_node(1,mlev)
ldu3(3)=max_node(2,mlev)
ENDIF
CALL ppm_alloc(mgfield(i,mlev)%bcvalue(iface)%pbcvalue,ldu3,iopt,info)
IF (info .NE. 0) THEN
info = ppm_error_fatal
CALL ppm_error(ppm_err_alloc,'ppm_poiss_mg_init', &
& 'Problem with the BOUNDARY alloc.',__LINE__,info)
GOTO 9999
ENDIF
ENDDO
DO iface=1,2*ppm_dim
IF (iface.EQ.1.OR.iface.EQ.2) THEN
direc(1)=2
direc(2)=3
elseif (iface.EQ.3.OR.iface.EQ.4) THEN
direc(1)=1
direc(2)=3
else
direc(1)=1
direc(2)=2
endif
DO ipoint=1,max_node(direc(1),mlev)
DO jpoint=1,max_node(direc(2),mlev)
DO ilda=1,vecdim
IF (mlev.EQ.1) THEN
!PRINT *,'bcef',ilda,iface,bcvalue(ilda,iface,ipoint,jpoint)
mgfield(i,mlev)%bcvalue(iface)%pbcvalue(ilda,ipoint,jpoint)=bcvalue(ilda,iface,ipoint,jpoint)
!PRINT *,ipoint,jpoint,iface,ilda,bcvalue(ilda,iface,ipoint,jpoint),mgfield(i,mlev)%bcvalue(iface)%pbcvalue(ilda,ipoint,jpoint)
! Print*,'Boundariu',bcvalue(ilda,iface,ipoint,jpoint)
ELSE
IF(bcdef_vec(ilda,i,iface).EQ.ppm_param_bcdef_neumann) THEN
mgfield(i,mlev)%bcvalue(iface)%pbcvalue(ilda,ipoint,jpoint)=&
& mgfield(i,mlev-1)%bcvalue(iface)%pbcvalue(ilda,2*ipoint-1,2*jpoint-1)
ELSE
!NO CORRECTIONS FOR THE DIRICHLET
mgfield(i,mlev)%bcvalue(iface)%pbcvalue(ilda,ipoint,jpoint)=0.0_MK
! If (mlev.EQ.5) then
! Print *,ipoint,jpoint,mgfield(i,mlev)%bcvalue(iface)%pbcvalue(ipoint,jpoint)
! endif
ENDIF
ENDIF
ENDDO
ENDDO
enddo
ENDDO!faces
ENDIF !lperiodic
#endif
iopt = ppm_param_alloc_fit
ldl3(1) = 1-ghostsize(1)
ldl3(2) = 1-ghostsize(2)
ldl3(3) = 1-ghostsize(3)
ldu3(1) = mesh%nnodes(1,idom)+ghostsize(1)
ldu3(2) = mesh%nnodes(2,idom)+ghostsize(2)
ldu3(3) = mesh%nnodes(3,idom)+ghostsize(3)
CALL ppm_alloc(mgfield(i,mlev)%mask_red,ldl3,ldu3,iopt,info)
IF (info .NE. 0) THEN
info = ppm_error_fatal
CALL ppm_error(ppm_err_alloc,'ppm_poiss_mg_init', &
& 'Problem with the mask alloc.',__LINE__,info)
GOTO 9999
ENDIF
iopt = ppm_param_alloc_fit
ldl3(1) = 1-ghostsize(1)
ldl3(2) = 1-ghostsize(2)
ldl3(3) = 1-ghostsize(3)
ldu3(1) = mesh%nnodes(1,idom)+ghostsize(1)
ldu3(2) = mesh%nnodes(2,idom)+ghostsize(2)
ldu3(3) = mesh%nnodes(3,idom)+ghostsize(3)
CALL ppm_alloc(mgfield(i,mlev)%mask_black,ldl3,ldu3,iopt,info)
IF (info .NE. 0) THEN
info = ppm_error_fatal
CALL ppm_error(ppm_err_alloc,'ppm_poiss_mg_init', &
& 'Problem with mask alloc.',__LINE__,info)
GOTO 9999
ENDIF
!----------------------------------------------------------------
!Filling the mask for communication (red black)
!-----------------------------------------------------------------
DO iz=1-ghostsize(3),&
mesh%nnodes(3,idom)+ghostsize(3)
DO iy=1-ghostsize(2),&
& mesh%nnodes(2,idom)+ghostsize(2)
DO ix=1-ghostsize(1),&
& mesh%nnodes(1,idom)+ghostsize(1)
IF (MOD(ix+iy+iz,2).EQ.0) THEN
mgfield(i,mlev)%mask_red(ix,iy,iz)=.TRUE.
mgfield(i,mlev)%mask_black(ix,iy,iz)=.FALSE.
!mgfield(i,mlev)%mask_black(ix,iy,iz)=.TRUE.
ELSE
mgfield(i,mlev)%mask_red(ix,iy,iz) = .FALSE.
mgfield(i,mlev)%mask_black(ix,iy,iz) = .TRUE.
ENDIF
ENDDO
ENDDO
ENDDO
ENDDO!DO i=1,nsubs
#endif
factor(:)=2
mesh_id_g(mlev)=lmesh_id
meshid_g(mlev)=meshid
newmeshid=-1
IF (mlev.LT.maxlev) THEN
!Print *,'dfj',meshid,ppm_param_mesh_coarsen,factor,newmeshid
CALL ppm_mesh_derive(topoid,meshid,newmeshid,&
& ppm_param_mesh_coarsen,factor,info)
lmesh_id = newmeshid
meshid = topo%mesh(lmesh_id)%ID
ENDIF
ENDDO!DO mlev=1,maxlev
!----------------------------------------------------------------------
! Return
!----------------------------------------------------------------------
9999 CONTINUE
CALL substop('ppm_mg_init',t0,info)
RETURN
#if __DIM == __SFIELD
#if __MESH_DIM == __2D
#if __KIND == __SINGLE_PRECISION
END SUBROUTINE ppm_mg_init_2d_sca_s
#elif __KIND == __DOUBLE_PRECISION
END SUBROUTINE ppm_mg_init_2d_sca_d
#endif
#elif __MESH_DIM == __3D
#if __KIND == __SINGLE_PRECISION
END SUBROUTINE ppm_mg_init_3d_sca_s
#elif __KIND == __DOUBLE_PRECISION
END SUBROUTINE ppm_mg_init_3d_sca_d
#endif
#endif
#elif __DIM == __VFIELD
#if __MESH_DIM == __2D
#if __KIND == __SINGLE_PRECISION
END SUBROUTINE ppm_mg_init_2d_vec_s
#elif __KIND == __DOUBLE_PRECISION
END SUBROUTINE ppm_mg_init_2d_vec_d
#endif
#elif __MESH_DIM == __3D
#if __KIND == __SINGLE_PRECISION
END SUBROUTINE ppm_mg_init_3d_vec_s
#elif __KIND == __DOUBLE_PRECISION
END SUBROUTINE ppm_mg_init_3d_vec_d
#endif
#endif
#endif