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 12 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(:) 25 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 33 34 PetscInt :: i 35 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 46 47 pv0 => v0 48 pJ => J 49 pinvJ => invJ 50 51 PetscCallA(PetscInitialize(ierr)) 52 PetscCallA(DMPlexCreate(PETSC_COMM_WORLD, dm, ierr)) 53 PetscCallA(PetscObjectSetName(dm, 'testplex', ierr)) 54 PetscCallA(DMSetDimension(dm, dim, ierr)) 55 56 PetscCallA(DMPlexCreateFromDAG(dm, depth, numPoints, coneSize, cones,coneOrientations, vertexCoords, ierr)) 57 58 PetscCallA(DMPlexInterpolate(dm, dmi, ierr)) 59 PetscCallA(DMPlexCopyCoordinates(dm, dmi, ierr)) 60 PetscCallA(DMDestroy(dm, ierr)) 61 dm = dmi 62 63 PetscCallA(DMView(dm, PETSC_VIEWER_STDOUT_WORLD, ierr)) 64 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 70 71 PetscCallA(PetscFVCreate(PETSC_COMM_WORLD, fvm, ierr)) 72 PetscCallA(PetscFVSetUp(fvm, ierr)) 73 PetscCallA(PetscFVDestroy(fvm, ierr)) 74 75 PetscCallA(DMDestroy(dm, ierr)) 76 PetscCallA(PetscFinalize(ierr)) 77 end program main 78 79!/*TEST 80! 81! test: 82! suffix: 0 83! 84!TEST*/ 85