PROGRAM ParAssembly
	!
	!      Mesh : Q1_4cells.msh
	!
	!
	!      7----------10---------6
	!      |          |          |
	!      |          |          |
	!      |     3    |    2     |
	!      |          |          |
	!      |          |          |
	!     11----------12---------9
	!      |          |          |
	!      |          |          |
	!      |     0    |    1     |
	!      |          |          |
	!      |          |          |
	!      4----------8----------5
	!
	!  Unknowns in u_global = [u9,u11,u12]
	!
	implicit none
#include <petsc/finclude/petsc.h90>
	!#include <petsc/finclude/petscsys.h>

	PetscInt, parameter :: ndime = 2, ndof = 1, nodes = 4
	PetscInt, parameter :: edges = 4, ndof_e=0
	PetscInt, parameter :: faces = 1, ndof_f=0

	PetscErrorCode      :: ierr
	PetscErrorCode      :: ierr_petsc
	PetscMPIInt         :: nranks,rank
	DM                  :: dm,dm_distrib
	DMLabel									:: dm_label
	Mat                 :: A
	!Vec                 :: u_global
	PetscBool           :: interpolate
	PetscSection        :: section
	PetscScalar         :: trace
	PetscViewer :: viewer

	PetscSF :: sf

	PetscInt, pointer :: cone(:),closure(:),indices(:)
	PetscInt, target, dimension(:),allocatable :: numComp ! We define a scalar field
	PetscInt, pointer :: pNumComp(:)
	PetscInt, target, dimension(:), allocatable :: numDof  ! number of dof per point
	PetscInt, pointer :: pNumDof(:)
	PetscInt, target, dimension(1) :: bcField
	PetscInt, pointer :: pBcField(:)
	IS, target, dimension(1) :: bcCompIS
	IS, pointer :: pBcCompIS(:)
	IS, target, dimension(1) :: bcPointIS
	IS, pointer :: pBcPointIS(:)

	PetscInt , dimension(nodes*ndof+edges*ndof_e+faces*ndof_f)    :: local_ldof
	PetscReal, dimension( nodes*ndof+edges*ndof_e+faces*ndof_f, nodes*ndof+edges*ndof_e+faces*ndof_f ) :: elstf

	PetscInt       :: nfields, numBC, dofs, dofsfix, offset
	PetscInt       :: cStart, cEnd, vStart, eStart, eEnd, vEnd, pStart, pEnd
	PetscInt :: fstart, fend, pstart1, pend1
	PetscInt       :: idface,r_sta, r_end
	PetscInt       :: size_cone,size_closure
	PetscInt       :: i,ic,ie,icell,idof,inod,ipoin,indx

	PetscInt       :: ntot, m, n, num_roots, num_leaves
	PetscInt, pointer :: locpoints(:)
	!	PetscSFNode,pointer :: remotepoints(:)

	ISLocalToGlobalMapping :: ltog

	! cone and support stuff
	PetscInt :: conesize, supportsize
	PetscInt, pointer :: conearr(:),supportarr(:) ! though 3 is enough
	PetscInt, pointer :: orientations(:)

	PetscInt :: nproc_nbr, numprocmax
	PetscInt, allocatable :: procs_nbr(:), numprocs(:), indices_nbr(:,:)

	!*************************************
	!--- Beginning of the program
	!*************************************
	call PetscInitialize(PETSC_NULL_CHARACTER,ierr)

	call MPI_Comm_size(PETSC_COMM_WORLD,nranks,ierr)
	call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)


	!*************************************
	!--- Creating Mesh
	!*************************************

	interpolate = PETSC_TRUE
	call DMPlexCreateGmshFromFile(PETSC_COMM_WORLD,"Q1_4cells.msh",interpolate,dm,ierr)

	call DMView(dm,PETSC_VIEWER_STDOUT_WORLD, ierr)

	!	if(rank==0)then
	!
	!	call DMPlexGetChart(dm, pstart1, pend1, ierr); CHKERRQ(ierr);
	!  call DMPlexGetDepthStratum(dm,0,vStart,vEnd,ierr);CHKERRQ(ierr);  ! Start & End vertices
	!  call DMPlexGetDepthStratum(dm,1,eStart,eEnd,ierr);CHKERRQ(ierr);  ! Start & End vertices
	!  call DMPlexGetDepthStratum(dm,2,fStart,fEnd,ierr);CHKERRQ(ierr);  ! Start & End vertices
	!  call DMPlexGetHeightStratum(dm,0,cStart,cEnd,ierr);CHKERRQ(ierr); ! Start & End cells
	!
	!	print *,'rank=', rank, 'all point start and end:', pstart1, pend1
	!	print *,'rank=', rank, 'vertex start and end:', vstart, vend
	!	print *,'rank=', rank, 'edge start and end:', estart, eend
	!	print *,'rank=', rank, 'face start and end:', fstart, fend
	!	print *,'rank=', rank, 'cell start and end:', cstart, cend
	!
	!!	pause
	!  !--- Loop over cells
	!  call PetscPrintf(PETSC_COMM_WORLD,"\n Loop over cells \n\n",ierr)
	!
	!  do icell = cStart, cEnd-1
	!
	!     call DMPlexGetTransitiveClosure(dm,icell,PETSC_TRUE,closure,ierr)
	!     size_closure = size(closure)/2  !closure contains (points,orientation)
	!
	!		 print *,'closure and orientations', closure(1:size(closure))
	!
	!		 !call DMPlexGetConeSize(dm, icell, conesize, ierr)
	!		 !call DMPlexGetCone(dm,icell,conearr, ierr)
	!
	!		 !call DMPlexGetConeOrientation(dm, icell,orientations, ierr)
	!
	!		 !print *, 'conesize',conesize, conearr(1:conesize)
	!		 !print *, 'conesize',conesize, orientations(1:conesize)
	!
	!		 !if(rank==0)		 print *,'size of closure', size_closure
	!
	!			 print *, 'cell index:', icell
	!     inod = 0
	!     do ic = 1 , size_closure
	!        ipoin = closure((ic-1)*2+1)
	!				if(ipoin>=estart .and. ipoin<eend)then
	!					!print *, 'this is a vertex', ipoin, 'orrientation',closure(2*ic)
	!					print *, 'this is a edge', ipoin, 'orrientation',closure(2*ic)
	!		      call DMPlexGetConeSize(dm, ipoin, conesize, ierr)
	!		     call DMPlexGetCone(dm,ipoin,conearr, ierr)
	!		     call DMPlexGetConeOrientation(dm, ipoin,orientations, ierr)
	!		 print *, 'conesize',conesize, conearr(1:conesize)
	!		 print *, 'conesize',conesize, orientations(1:conesize)
	!
	!
	!
	!				endif
	!     enddo
	!     call DMPlexRestoreTransitiveClosure(dm,icell,PETSC_TRUE,closure,ierr)
	!	 enddo
	! endif

	call PetscPrintf(PETSC_COMM_WORLD,"\n Distribute the mesh ***************** \n\n",ierr)
	!--- Distributing the mesh over processes
	if (nranks > 1) then
		!call DMPlexDistribute(dm,0,PETSC_NULL_OBJECT,dm_distrib,ierr);
		call DMPlexDistribute(dm,0,sf,dm_distrib,ierr);
		call DMDestroy(dm,ierr);
		dm = dm_distrib;
	endif

	!	call PetscSFView(sf, PETSC_VIEWER_STDOUT_WORLD, ierr)

	!call PetscSFGetGraph(sf,num_roots, num_leaves,locpoints, remotepoints)
	!call PetscSFDestroy(sf, ierr)

	!*************************************
	!--- Set Petsc Layout
	!*************************************

	! !--- We create additional labels for vertices on labeled faces
	!  call DMCreateLabel(dm,"AllVertBC",ierr)
	!
	!  !--- We get the range Id of Sieve Points which are faces/edges in 2D
	!  call DMPlexGetDepthStratum(dm,ndime-1,eStart,eEnd,ierr)
	!  do ie = eStart, eEnd-1
	!      !- This returns FaceNum if ifac is external face or -1 if not
	!      !- IMPORTANT : Physical Entities must be defined in .geo to have FaceNums
	!      call DMGetLabelValue(dm,"Face Sets",ie,idface,ierr)
	!      if (idface .gt. 0) then
	!         !-- For Dirichlet BCs on vertices we want to extract vertices
	!         call DMPlexGetConeSize(dm,ie,size_cone,ierr)
	!         call DMPlexGetCone(dm,ie,cone,ierr)
	!         do i=1,size_cone
	!            !-- Select only 2 faces for Bcs
	!            if (idface .eq. 1 .or. idface .eq. 3) then
	!               call DMSetLabelValue(dm,"AllVertBC",cone(i),1,ierr)
	!            endif
	!         enddo        
	!      endif
	!  enddo

	call DMView(dm,PETSC_VIEWER_STDOUT_WORLD, ierr)


	!---- Defining a Scalar Field on vertices
	nFields = 2    ! Num of fields
	allocate(numComp(nFields))
	numComp(1) = 1 ! Num of components
	numComp(2) = 1 ! Num of components
	pnumComp => numComp

	allocate(numdof(nFields * (ndime+1)))
	do i = 1, nFields*(ndime+1)
		numDof(i) = 0
	end do

	numDof(0*(ndime+1)+1) = 1  !1 DOFs on vertices
	numDof(0*(ndime+1)+2) = 0  !1 DOFs on edges
	numDof(0*(ndime+1)+3) = 0  !1 DOFs on edges
	numDof(1*(ndime+1)+1) = 1  !1 DOFs on vertices
	numDof(1*(ndime+1)+2) = 0  !1 DOFs on edges
	numDof(1*(ndime+1)+3) = 0  !1 DOFs on edges
	!numDof(1*(ndime+1)+1) = 1  !1 DOFs on vertices
	!numDof(2*(ndime+1)+1) = 1  !1 DOFs on vertices
	pnumDof => numDof

	!--- Defining Boundary Conditions on faces 1 and 3
	! numBC = 1
	numBC = 0
	bcField(1) = 0  ! Field 0 : u
	!--- Data structure defining contrained components of u
	!  call ISCreateStride(PETSC_COMM_WORLD, 1, 0, 1, bcCompIS(1), ierr);
	!--- Data structure with num of SievePoints with val=1 in Label 'Vert Sets'
	!  call DMGetStratumIS(dm, "AllVertBC", 1, bcPointIS(1),ierr);
	!  pBcField => bcField
	!  pBcCompIS => bcCompIS
	!  pBcPointIS => bcPointIS

	!--- Create a PetscSection with this data layout
	!call DMPlexCreateSection(dm, ndime, nFields, pNumComp,        &
	!        pNumDof, numBC, pBcField, pBcCompIS, pBcPointIS,     &
	!          PETSC_NULL_OBJECT, section, ierr)
	!  call DMPlexCreateSection(dm, ndime, nFields, pNumComp,        &
	!           pNumDof, 0, PETSC_NULL_OBJECT, PETSC_NULL_OBJECT, PETSC_NULL_OBJECT,     &
	!           PETSC_NULL_OBJECT, section, ierr)
	call DMPlexCreateSection(dm, ndime, nFields, pNumComp,        &
		pNumDof, numBC, pBcField, pBcCompIS, pBcPointIS,     &
		PETSC_NULL_OBJECT, section, ierr)

	!--- Name the Field variables
	call PetscSectionSetFieldName(section, 0, 'u', ierr)

	!--- Tell DM to use this section
	call DMSetDefaultSection(dm,section,ierr)

	!--- Visualizing Section
	call PetscPrintf(PETSC_COMM_WORLD,"\n Visualizing Section ... \n",ierr)
	call PetscSectionView(section, PETSC_VIEWER_STDOUT_WORLD, ierr)

	call PetscPrintf(PETSC_COMM_WORLD,"\n dm label view ... \n",ierr)
	call DMGetLabel(dm, 'ghost', dm_label, ierr)
	call DMLabelView(dm_label, PETSC_VIEWER_STDOUT_WORLD, ierr)


	!*************************************
	!--- Matrix and RHS Assembly
	!*************************************

	!--- Create a Sparse Matrix with sparsity struct. based on dm
	call DMCreateMatrix(dm,A,ierr)
	!call DMCreateGlobalVector(dm,u_global,ierr)
	!--- Set Mapping Local to Global Indexes
	call DMGetLocalToGlobalMapping(dm,ltog,ierr)
	call MatSetLocalToGlobalMapping(A,ltog,ltog,ierr);

	call ISLocalToGlobalMappingView(ltog, PETSC_VIEWER_STDOUT_WORLD, ierr);

	call MatZeroEntries(A,ierr)

	call MatGetOwnershipRange(A, r_sta, r_end, ierr)
	print *, 'rank-start-end',rank,r_sta, r_end

	call MPI_Barrier(MPI_COMM_WORLD, ierr)

	call ISLocalToGlobalMpngGetInfoSize(ltog, nproc_nbr, numprocmax, ierr)
!	call ISLocalToGlobalMappingGetInfo(ltog, nproc_nbr, ierr)
	!numprocmax = 4
	print *, 'rank-nproc-max', rank, nproc_nbr, numprocmax
	allocate(procs_nbr(nproc_nbr))
	allocate(numprocs(nproc_nbr))
	allocate(indices_nbr(nproc_nbr,numprocmax))
	call ISLocalToGlobalMappingGetInfo(ltog, nproc_nbr,procs_nbr, numprocs,	indices_nbr, ierr)

	call MPI_Barrier(MPI_COMM_WORLD, ierr)
	if(rank==0) then
		print *, ' '
		print *, 'rank', rank
	print *, 'nproc_nbr',nproc_nbr
	print *, 'procs_nbr', procs_nbr
	print *, 'numprocs', numprocs
	print *, 'indices', indices_nbr
endif
	call MPI_Barrier(MPI_COMM_WORLD, ierr)
	if(rank==1) then
		print *, ' '
		print *, 'rank', rank
	print *, 'nproc_nbr',nproc_nbr
	print *, 'procs_nbr', procs_nbr
	print *, 'numprocs', numprocs
	print *, 'indices', indices_nbr
endif

	call MPI_Barrier(MPI_COMM_WORLD, ierr)
	call ISLocalToGlobalMappingDestroy(ltog,ierr);
	!************ Starting assembly
!	call DMPlexGetChart(dm, pstart1, pend1, ierr); CHKERRQ(ierr);
!	call DMPlexGetDepthStratum(dm,0,vStart,vEnd,ierr);CHKERRQ(ierr);  ! Start & End vertices
!	call DMPlexGetDepthStratum(dm,1,eStart,eEnd,ierr);CHKERRQ(ierr);  ! Start & End vertices
!	call DMPlexGetDepthStratum(dm,2,fStart,fEnd,ierr);CHKERRQ(ierr);  ! Start & End vertices
!	call DMPlexGetHeightStratum(dm,0,cStart,cEnd,ierr);CHKERRQ(ierr); ! Start & End cells
!
!	print *,'rank=', rank, 'all point start and end:', pstart1, pend1
!
!	print *,'rank=', rank, 'point start and end:', vstart, vend
!	print *,'rank=', rank, 'edge start and end:', estart, eend
!	print *,'rank=', rank, 'face start and end:', fstart, fend
!	print *,'rank=', rank, 'cell start and end:', cstart, cend
!
!	!	pause
!	!--- Loop over cells
!	call PetscPrintf(PETSC_COMM_WORLD,"\n Loop over cells \n\n",ierr)
!
!	do icell = cStart, cEnd-1
!
!		call DMPlexGetTransitiveClosure(dm,icell,PETSC_TRUE,closure,ierr)
!		size_closure = size(closure)/2  !closure contains (points,orientation)
!
!
!		call DMPlexGetConeSize(dm, icell, conesize, ierr)
!		call DMPlexGetCone(dm,icell,conearr, ierr)
!
!		print *,'rank', rank, 'conesize',conesize, conearr(1:conesize)
!
!		if(rank==0)		 print *,'size of closure', size_closure
!
!		inod = 0
!		do ic = 1 , size_closure
!			ipoin = closure((ic-1)*2+1)
!			call PetscSectionGetDof(section,ipoin,dofs,ierr)
!			call PetscSectionGetOffSet(section,ipoin,offset,ierr)
!			if(rank==0)	print *,'dof and offset', dofs, offset, 'at sieve point',				ipoin
!			if (dofs .ne. 0) then  ! To get Vertices
!				! !--- Get number of vertex constrained dofs
!				! call PetscSectionGetConstraintDof(section,ipoin,dofsfix,ierr)
!				! if (dofsfix /= 0) then
!				!     call PetscSectionGetConstraintIndicesF90(section,ipoin,indices,ierr)
!				! endif
!				!--- Get location of point in local vector
!				call DMPlexGetPointLocal(dm,ipoin,pStart,pEnd,ierr)
!				if(rank==0) print *,'dof',ipoin, 'pstart, and pend', pstart, pend
!				!--- Get location of DOF(s) number(s)
!				do idof = 1 , dofs
!					inod = inod + 1
!					i = inod 
!					local_ldof(i) = offset + idof !idof  !Location of DOF in local_locdof()
!				enddo
!				!if (dofsfix /= 0) then ! If constrained dofs
!				!   do indx = 1 , dofsfix
!				!      idof = indices(indx) + 1
!				!      i = (inod-1)*dofs + idof
!				!      ! Change local_dof(i) to -local_dof(i) if constrained
!				!      ! Will be used to modify RHS for BCs
!				!      local_ldof(i) = -local_ldof(i)
!				!   enddo
!				!endif
!				!if (dofsfix /= 0) then
!				!   call PetscSectionRestoreConstraintIndicesF90(section,ipoin,indices,ierr)
!				!endif
!			endif
!		enddo
!		call DMPlexRestoreTransitiveClosure(dm,icell,PETSC_TRUE,closure,ierr)
!
!
!		elstf = 1.0 ! Elementary Stiffness
!
!		ntot = nodes*ndof+edges*ndof_e + faces*ndof_f
!
!		do i = 1 , ntot
!			if (local_ldof(i) > 0) then !Set only if DOF not constrained ??
!				!print *,'local_global list rank=',rank,  local_ldof
!				call MatSetValuesLocal(A,1,local_ldof(i)-1,ntot,local_ldof(:)-1, &
!					elstf(i,:),ADD_VALUES,ierr)
!			endif
!		enddo
!
!	enddo
!
!	call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)
!	call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)
!
!	!	call MatGetOwnershipRange(A, m, n, ierr)
!
!	!	print *, 'rank', rank, m,n, (n-m)
!
!
!	call MatView(A,PETSC_VIEWER_STDOUT_WORLD,ierr)
!
!	call MatGetTrace(A, trace, ierr)
!
!	if(rank==0)	print *, rank, 'trace=', trace
!
!	call PetscViewerCreate(MPI_COMM_WORLD, viewer, ierr_petsc)
!	!call PetscViewerSetType(viewer,PETSCVIEWERBINARY, ierr_petsc)
!	!call PetscViewerSetFormat(viewer,PETSC_VIEWER_DEFAULT, ierr_petsc)
!	!call PetscViewerFileSetMode(viewer, FILE_MODE_WRITE, ierr_petsc)
!	!call PetscViewerFileSetName(viewer, 'GmatA.m', ierr_petsc)
!
!	!------------------- ascii format
!	call PetscViewerASCIIOpen(Petsc_comm_world, 'xx.m', viewer, ierr_petsc)
!	call PetscViewerSetFOrmat(viewer, PETSC_VIEWER_ASCII_MATLAB, ierr_petsc)
!
!	call Matview(A, viewer, ierr_petsc)
!
!	call PetscViewerDestroy(viewer, ierr_petsc)
!

	!*************************************
	!--- Terminating Program
	!*************************************


	!--- Destroying objects
	call PetscPrintf(PETSC_COMM_WORLD,"\n ... Destroying objects \n",ierr)
	call DMDestroy(dm,ierr)
	call ISDestroy(bcCompIS, ierr)
	call ISDestroy(bcPointIS, ierr)
	!call MatDestroy(A,ierr)
	!call VecDestroy(u_global,ierr)

	call PetscPrintf(PETSC_COMM_WORLD,"\n ---- End the Program \n",ierr)
	call PetscFinalize(ierr)

	!  return

ENDPROGRAM