1! setting up 3-D DMPlex using DMPlexCreateFromDAG(), DMPlexInterpolate() and 2! DMPlexComputeCellGeometryFVM() 3! Contributed by Adrian Croucher <a.croucher@auckland.ac.nz> 4#include <petsc/finclude/petscdmplex.h> 5program main 6 use petscdmplex 7 implicit none 8 DM :: dm, dmi 9 PetscFV :: fvm 10 PetscInt, parameter :: dim = 3 11 12 PetscInt :: depth = 1 13 PetscErrorCode :: ierr 14 PetscInt, dimension(2) :: numPoints 15 PetscInt, dimension(14) :: coneSize 16 PetscInt, dimension(16) :: cones 17 PetscInt, dimension(16) :: coneOrientations 18 PetscScalar, dimension(36) :: vertexCoords 19 PetscReal :: vol = 0. 20 PetscReal, target, dimension(dim) :: centroid 21 PetscReal, target, dimension(dim) :: normal 22 PetscReal, pointer :: pcentroid(:) 23 PetscReal, pointer :: pnormal(:) 24 25 PetscReal, target, dimension(dim) :: v0 26 PetscReal, target, dimension(dim*dim) :: J 27 PetscReal, target, dimension(dim*dim) :: invJ 28 PetscReal, pointer :: pv0(:) 29 PetscReal, pointer :: pJ(:) 30 PetscReal, pointer :: pinvJ(:) 31 PetscReal :: detJ 32 33 PetscInt :: i 34 35 numPoints = [12, 2] 36 coneSize = [8, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0] 37 cones = [2, 5, 4, 3, 6, 7, 8, 9, 3, 4, 11, 10, 7, 12, 13, 8] 38 coneOrientations = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0] 39 vertexCoords = [ & 40 & -0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, -0.5, 1.0, 0.0, & 41 & -0.5, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, 1.0, -0.5, 1.0, 1.0, & 42 & 0.5, 0.0, 0.0, 0.5, 1.0, 0.0, 0.5, 0.0, 1.0, 0.5, 1.0, 1.0] 43 pcentroid => centroid 44 pnormal => normal 45 46 pv0 => v0 47 pJ => J 48 pinvJ => invJ 49 50 PetscCallA(PetscInitialize(ierr)) 51 PetscCallA(DMPlexCreate(PETSC_COMM_WORLD, dm, ierr)) 52 PetscCallA(PetscObjectSetName(dm, 'testplex', ierr)) 53 PetscCallA(DMSetDimension(dm, dim, ierr)) 54 55 PetscCallA(DMPlexCreateFromDAG(dm, depth, numPoints, coneSize, cones, coneOrientations, vertexCoords, ierr)) 56 57 PetscCallA(DMPlexInterpolate(dm, dmi, ierr)) 58 PetscCallA(DMPlexCopyCoordinates(dm, dmi, ierr)) 59 PetscCallA(DMDestroy(dm, ierr)) 60 dm = dmi 61 62 PetscCallA(DMView(dm, PETSC_VIEWER_STDOUT_WORLD, ierr)) 63 64 do i = 0, 1 65 PetscCallA(DMPlexComputeCellGeometryFVM(dm, i, vol, pcentroid, pnormal, ierr)) 66 write (*, '(a, i2, a, f8.4, a, 3(f8.4, 1x))') 'cell: ', i, ' volume: ', vol, ' centroid: ', pcentroid(1), pcentroid(2), pcentroid(3) 67 PetscCallA(DMPlexComputeCellGeometryAffineFEM(dm, i, pv0, pJ, pinvJ, detJ, ierr)) 68 end do 69 70 PetscCallA(PetscFVCreate(PETSC_COMM_WORLD, fvm, ierr)) 71 PetscCallA(PetscFVSetUp(fvm, ierr)) 72 PetscCallA(PetscFVDestroy(fvm, ierr)) 73 74 PetscCallA(DMDestroy(dm, ierr)) 75 PetscCallA(PetscFinalize(ierr)) 76end program main 77 78!/*TEST 79! 80! test: 81! suffix: 0 82! 83!TEST*/ 84