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: #include <petsc/finclude/petscdmplex.h>
5: program main
6: use petscdmplex
7: implicit none
8: DM :: dm, dmi
9: PetscFV :: fvm
10: PetscInt, parameter :: dim = 3
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(:)
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
33: PetscInt :: i
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
46: pv0 => v0
47: pJ => J
48: pinvJ => invJ
50: PetscCallA(PetscInitialize(ierr))
51: PetscCallA(DMPlexCreate(PETSC_COMM_WORLD, dm, ierr))
52: PetscCallA(PetscObjectSetName(dm, 'testplex', ierr))
53: PetscCallA(DMSetDimension(dm, dim, ierr))
55: PetscCallA(DMPlexCreateFromDAG(dm, depth, numPoints, coneSize, cones, coneOrientations, vertexCoords, ierr))
57: PetscCallA(DMPlexInterpolate(dm, dmi, ierr))
58: PetscCallA(DMPlexCopyCoordinates(dm, dmi, ierr))
59: PetscCallA(DMDestroy(dm, ierr))
60: dm = dmi
62: PetscCallA(DMView(dm, PETSC_VIEWER_STDOUT_WORLD, ierr))
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
70: PetscCallA(PetscFVCreate(PETSC_COMM_WORLD, fvm, ierr))
71: PetscCallA(PetscFVSetUp(fvm, ierr))
72: PetscCallA(PetscFVDestroy(fvm, ierr))
74: PetscCallA(DMDestroy(dm, ierr))
75: PetscCallA(PetscFinalize(ierr))
76: end program main
78: !/*TEST
79: !
80: ! test:
81: ! suffix: 0
82: !
83: !TEST*/