Actual source code: ex3f90.F90

      1: ! setting up 3-D DMPlex using DMPlexCreateFromDAG(), DMPlexInterpolate() and
      2: ! DMPlexComputeCellGeometryFVM()
      3: ! Contributed by Adrian Croucher <a.croucher@auckland.ac.nz>
      4:       program main
      5:       use petscsys
      6:       use petscdmplex
      7: #include <petsc/finclude/petscsys.h>
      8: #include <petsc/finclude/petscdmplex.h>
      9:       implicit none
     10:       DM :: dm, dmi
     11:       PetscFV :: fvm
     12:       PetscInt, parameter :: dim = 3

     14:       PetscInt :: depth = 1
     15:       PetscErrorCode :: ierr
     16:       PetscInt, dimension(2) :: numPoints
     17:       PetscInt, dimension(14) :: coneSize
     18:       PetscInt, dimension(16) :: cones
     19:       PetscInt, dimension(16) :: coneOrientations
     20:       PetscScalar, dimension(36) :: vertexCoords
     21:       PetscReal ::  vol = 0.
     22:       PetscReal, target, dimension(dim)  :: centroid
     23:       PetscReal, target, dimension(dim)  :: normal
     24:       PetscReal, pointer :: pcentroid(:)
     25:       PetscReal, pointer :: pnormal(:)

     27:       PetscReal, target, dimension(dim)  :: v0
     28:       PetscReal, target, dimension(dim*dim)  :: J
     29:       PetscReal, target, dimension(dim*dim)  :: invJ
     30:       PetscReal, pointer :: pv0(:)
     31:       PetscReal, pointer :: pJ(:)
     32:       PetscReal, pointer :: pinvJ(:)
     33:       PetscReal :: detJ

     35:       PetscInt :: i

     37:       numPoints = [12, 2]
     38:       coneSize  = [8, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
     39:       cones = [2,5,4,3,6,7,8,9,  3,4,11,10,7,12,13,8]
     40:       coneOrientations = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,0]
     41:       vertexCoords = [                                                  &
     42:      &  -0.5,0.0,0.0, 0.0,0.0,0.0, 0.0,1.0,0.0, -0.5,1.0,0.0,           &
     43:      &  -0.5,0.0,1.0, 0.0,0.0,1.0, 0.0,1.0,1.0, -0.5,1.0,1.0,           &
     44:      &   0.5,0.0,0.0, 0.5,1.0,0.0, 0.5,0.0,1.0,  0.5,1.0,1.0 ]
     45:       pcentroid => centroid
     46:       pnormal => normal

     48:       pv0 => v0
     49:       pJ => J
     50:       pinvJ => invJ

     52:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
     53:       if (ierr .ne. 0) then
     54:         print*,'Unable to initialize PETSc'
     55:         stop
     56:       endif

     58:       call DMPlexCreate(PETSC_COMM_WORLD, dm, ierr);CHKERRA(ierr)
     59:       call PetscObjectSetName(dm, 'testplex', ierr);CHKERRA(ierr)
     60:       call DMSetDimension(dm, dim, ierr);CHKERRA(ierr)

     62:       call DMPlexCreateFromDAG(dm, depth, numPoints, coneSize, cones,coneOrientations, vertexCoords, ierr);CHKERRA(ierr)

     64:       call DMPlexInterpolate(dm, dmi, ierr);CHKERRA(ierr)
     65:       call DMPlexCopyCoordinates(dm, dmi, ierr);CHKERRA(ierr)
     66:       call DMDestroy(dm, ierr);CHKERRA(ierr)
     67:       dm = dmi

     69:       call DMView(dm, PETSC_VIEWER_STDOUT_WORLD, ierr);CHKERRA(ierr)

     71:       do i = 0, 1
     72:         call DMPlexComputeCellGeometryFVM(dm, i, vol, pcentroid, pnormal, ierr);CHKERRA(ierr)
     73:         write(*, '(a, i2, a, f8.4, a, 3(f8.4, 1x))') 'cell: ', i, ' volume: ', vol, ' centroid: ',pcentroid(1), pcentroid(2), pcentroid(3)
     74:         call DMPlexComputeCellGeometryAffineFEM(dm, i, pv0, pJ, pinvJ,detJ, ierr);CHKERRA(ierr)
     75:       end do

     77:       call PetscFVCreate(PETSC_COMM_WORLD, fvm, ierr);CHKERRA(ierr)
     78:       call PetscFVSetUp(fvm, ierr);CHKERRA(ierr)
     79:       call PetscFVDestroy(fvm, ierr);CHKERRA(ierr)

     81:       call DMDestroy(dm, ierr);CHKERRA(ierr)
     82:       call PetscFinalize(ierr)
     83:       end program main

     85: !/*TEST
     86: !
     87: !   test:
     88: !      suffix: 0
     89: !
     90: !TEST*/