program testDMView implicit none #include DM :: dmMesh PetscErrorCode :: ierr MPI_Comm :: comm integer :: topologyDim integer :: numCell integer :: numNodePerCell integer :: geometryDim integer :: numNode integer,allocatable :: cellVertices(:) PetscReal, allocatable :: coordinates(:) PetscBool :: interpolate integer :: numX,numY PetscReal :: delX,delY integer :: ic,jc,kc,dim integer :: numStrata IS :: cellIS ! cell id index set integer :: globalidStart,globalidEnd PetscReal, target, dimension(2) :: v0 PetscReal, target, dimension(4) :: J PetscReal, target, dimension(4) :: invJ PetscReal, pointer :: pv0(:), pJ(:), pinvJ(:) PetscReal :: detJ integer :: cellId,coneSize PetscSection :: coordSection Vec :: CoordVec PetscScalar,pointer :: coords(:) integer :: locCrdSize ! Begin program call PetscInitialize(PETSC_NULL_CHARACTER, ierr) comm = PETSC_COMM_SELF !------------------------------ ! Create a box mesh with quadrangle elements ! Mesh region = [0 1] x [0 1] ! numX - number of points in x-direction ! numY - number of points in y-direction !------------------------------ numX = 3 numY = 3 ! topology dimension is 2 since mesh is two dimentsional topologyDim=2 ! Number of Cells numCell = (numX-1)*(numY-1) ! Number of Nodes per Cell is 4 for quadrangle element numNodePerCell = 4 ! geometry dimensio is 2, since the mesh is in x-y plane geometryDim = 2 ! Grid Spacing in x and y direction delX = 1.0/(numX-1) delY = 1.0/(numY-1) ! Number of Nodes/vertices in the mesh numNode = numX*numY allocate(cellVertices(numCell*numNodePerCell)) allocate(coordinates(numNode*geometryDim)) interpolate = .true. ! Prepare mesh Coordiantes ! Mesh Nodes are numberd in x-direction first then in y direction kc=0 do jc=1,numY do ic=1,numX kc=kc+1 coordinates(kc) = (ic-1)*delX kc=kc+1 coordinates(kc) = (jc-1)*delY end do end do ! Form Connectivity ! -1 becuase of zero based indexing for nodes ! Cells are counted in X-direction first then in Y-direction ! In each cells the nodes are listed in counter clockwise direction ! For example, for cellid = 0, the node coordiantes are ! node1 = [0,0] node2 = [0,delX], node3 = [delX,delY], node4 = [0,delY] !kc=0 do jc=1,numY-1 do ic=1,numX-1 kc=ic + (numX-1)*(jc-1) cellVertices((kc-1)*numNodePerCell + 1) = ic + (jc-1)*numX -1 cellVertices((kc-1)*numNodePerCell + 2) = ic + (jc-1)*numX + 1 -1 cellVertices((kc-1)*numNodePerCell + 3) = ic + (jc)*numX + 1 -1 cellVertices((kc-1)*numNodePerCell + 4) = ic + (jc)*numX -1 end do end do !-------------------------------- ! Create DM using vertices and Connectivity ! Any entity of a mesh is identified as cell with dimension dim ! or a cell with codimension cdim ? ! Terminology Notes: ! vertex = 0 dimensional cell or dim codimension cell ! element = dim dimensional cell or 0 codimension cell ! facet = dim-1 dimensional cell ! ! node or point = 0 dimension topology ! edge or line = 1 dimension topology ! Triangle or qudrangle = 2 dimension topology ! Terahedron or pyramid or Any 3d element is = 3 dimension topology ! By defualt DM created from cell list has only 0 and dim dimensional entities ! use interpolateflag=.true. to create other entites !-------------------------------- call DMPlexCreateFromCellList(comm,topologydim,numCell,numNode,numNodePerCell, & interpolate,cellVertices,geometryDim,coordinates, & dmMesh, ierr) print*,'Surcessfully created DM From Cell List' !-------------------------------- ! Lets print the dmMesh using DMView to see how the DM information ! is organized !-------------------------------- call DMView(dmMesh,PETSC_VIEWER_STDOUT_WORLD, ierr) !------------------------------- ! You should see an output like this on screen ! Mesh in 2 dimensions: ! 0-cells: 9 ! = numX*numY ! 1-cells: 12 ! 2-cells: 4 ! = (numX-1)*(numY-1) ! Labels: ! depth: 3 strata of sizes (9, 12, 4) !-------------------------------- ! What is a strata ? < I do not clearly understand this > ! ! We can see below using the "depth" label how we can ! access the globalids of cells in a strata ! For example we can get the number of strata of a particular ! value or id ? (here it is dim=1) call DMPlexGetStratumSize(dmMesh,"depth",1,numstrata,ierr) print*,'number of strat of vlaue',1,'is',numstrata ! We can get the index sets of cells in the strata ! Note every entity is a cell and is enumerated continuously ! irrespective of its dimension ! Here we get the cell index set of cells with dimension = 1 call DMPlexGetStratumIS(dmMesh,"depth",1,cellIS,ierr) call ISView(cellIS,PETSC_VIEWER_STDOUT_WORLD, ierr) ! Here we print the start and endids of the different dimensional cells ! get the start and end ids of the 0-cells call DMPlexGetDepthStratum(dmMesh,0,globalidStart,globalidEnd,ierr) print*,'Node id start=',globalidStart,'node id end=',globalidEnd-1 ! get the start and end ids of the 1-cells, 1-cell here is a edge call DMPlexGetDepthStratum(dmMesh,1,globalidStart,globalidEnd,ierr) print*,'edge id start=',globalidStart,'edge id end=',globalidEnd-1 ! get the start and end ids of the 2-cells call DMPlexGetDepthStratum(dmMesh,2,globalidStart,globalidEnd,ierr) print*,'element id start=',globalidStart,'element id end=',globalidEnd-1 ! Once we have the start and end ids of dim dimensional cells ! we are ready for element by element assembly ! During the assembly we will need to compute the cellGeometry ! which can be done using the call below ! compute cell geometry for first cell cellId = 0 !------------------------- !<<<<<< Jacobian data seems weired >>>>>>>>>>>>>>>>>>> !------------------------ pv0 => v0 pJ => J pinvJ => invJ call DMPlexComputeCellGeometry(dmMesh,cellId,pv0,pJ,pinvJ,detJ,ierr) CHKERRQ(ierr) print*,'Translation Vector ::', v0 print*,'Jacobian ::',J print*,'invJ ::',invJ print*,'detJ::',detJ ! Alternatively,if the cell geometry is no implimented for your mesh element ! say a pyramid, you can try the follwoing appraoch. ! The following code is Fortran version of code in ! DMPlexComputeRectangleGeometry_private, defined in src/dm/impls/plex/plex.c print*,'Calling get ConeSize for cell to get the number of nodes per cell' call DMPlexGetConeSize(dmMesh, cellId,coneSize,ierr) print*,'Cone Size i.e., number of nodes per cell = ',conesize ! Get the local coordinates in a vector call DMGetCoordinatesLocal(dmMesh,coordVec,ierr) !call VecView(coordVec,PETSC_VIEWER_STDOUT_WORLD, ierr) call DMPlexGetCoordinateSection(dmMesh, coordSection,ierr) !call PetscSectionView(coordSection, PETSC_VIEWER_STDOUT_WORLD, ierr) ! Get access to data via a pointer using coordinate Vec and coordSection ! For a give cell Id ! in this case we access the vertex coordinates for cellId=0 call DMPlexVecGetClosure(dmMesh, coordSection, coordVec,cellId,coords,ierr) print*,'Cell Node Coordiantes in counter clockwise order' do ic=1,conesize print*,'Node=',ic,'crd=',coords((ic-1)*2+1:(ic)*2) end do ! Do the local calculations with coords dim = topologyDim do ic=1,dim v0(ic) = PetscRealPart(coords(ic)) ! Translation part of Affine vector is First Node end do do ic=0,dim-1 do jc=0,dim-1 J(ic*(dim) + jc + 1) = 0.5*(PetscRealPart(coords((ic*2+1)*dim+jc+1)) & - PetscRealPart(coords(0*dim+jc+1))) end do end do detJ = J(1)*J(4) - J(2)*J(3) invJ(1) = (1.0/detJ)*J(4) invJ(2) = -(1.0/detJ)*J(2) invJ(3) = -(1.0/detJ)*J(3) invJ(4) = (1.0/detJ)*J(1) print*,'Translation Vector ::', v0 print*,'Jacobian ::',J print*,'invJ ::',invJ print*,'detJ::',detJ ! restore the access call DMPlexVecRestoreClosure(dmMesh, coordSection, coordVec,cellId,coords,ierr) call PetscFinalize(ierr) end program testDMView