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: #include <petsc/finclude/petscsys.h>
      6: #include <petsc/finclude/petscdmplex.h>
      7:   use petscdmplex
      8:   implicit none
      9:   DM :: dm, dmi
     10:   PetscFV :: fvm
     11:   PetscInt, parameter :: dim = 3

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

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

     34:   PetscInt :: i

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

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

     51:   PetscCallA(PetscInitialize(ierr))
     52:   PetscCallA(DMPlexCreate(PETSC_COMM_WORLD, dm, ierr))
     53:   PetscCallA(PetscObjectSetName(dm, 'testplex', ierr))
     54:   PetscCallA(DMSetDimension(dm, dim, ierr))

     56:   PetscCallA(DMPlexCreateFromDAG(dm, depth, numPoints, coneSize, cones, coneOrientations, vertexCoords, ierr))

     58:   PetscCallA(DMPlexInterpolate(dm, dmi, ierr))
     59:   PetscCallA(DMPlexCopyCoordinates(dm, dmi, ierr))
     60:   PetscCallA(DMDestroy(dm, ierr))
     61:   dm = dmi

     63:   PetscCallA(DMView(dm, PETSC_VIEWER_STDOUT_WORLD, ierr))

     65:   do i = 0, 1
     66:     PetscCallA(DMPlexComputeCellGeometryFVM(dm, i, vol, pcentroid, pnormal, ierr))
     67:     write (*, '(a, i2, a, f8.4, a, 3(f8.4, 1x))') 'cell: ', i, ' volume: ', vol, ' centroid: ', pcentroid(1), pcentroid(2), pcentroid(3)
     68:     PetscCallA(DMPlexComputeCellGeometryAffineFEM(dm, i, pv0, pJ, pinvJ, detJ, ierr))
     69:   end do

     71:   PetscCallA(PetscFVCreate(PETSC_COMM_WORLD, fvm, ierr))
     72:   PetscCallA(PetscFVSetUp(fvm, ierr))
     73:   PetscCallA(PetscFVDestroy(fvm, ierr))

     75:   PetscCallA(DMDestroy(dm, ierr))
     76:   PetscCallA(PetscFinalize(ierr))
     77: end program main

     79: !/*TEST
     80: !
     81: !   test:
     82: !      suffix: 0
     83: !
     84: !TEST*/