xref: /petsc/src/dm/impls/plex/tutorials/ex1f90.F90 (revision dcb3e68992f1c4897946af7e8406e2b4165e50f2)
1c4762a1bSJed Brown      program DMPlexTestField
2c4762a1bSJed Brown#include "petsc/finclude/petscdmplex.h"
3c4762a1bSJed Brown#include "petsc/finclude/petscdmlabel.h"
4c4762a1bSJed Brown      use petscdmplex
5c4762a1bSJed Brown      use petscsys
6c4762a1bSJed Brown      implicit none
7c4762a1bSJed Brown
8c4762a1bSJed Brown      DM :: dm
9c4762a1bSJed Brown      DMLabel :: label
10c4762a1bSJed Brown      Vec :: u
11c4762a1bSJed Brown      PetscViewer :: viewer
12c4762a1bSJed Brown      PetscSection :: section
13c4762a1bSJed Brown      PetscInt :: dim,numFields,numBC
14c4762a1bSJed Brown      PetscInt :: i,val
15c4762a1bSJed Brown      DMLabel, pointer :: nolabel(:) => NULL()
16c4762a1bSJed Brown      PetscInt, target, dimension(3) ::  numComp
17c4762a1bSJed Brown      PetscInt, pointer :: pNumComp(:)
18c4762a1bSJed Brown      PetscInt, target, dimension(12) ::  numDof
19c4762a1bSJed Brown      PetscInt, pointer :: pNumDof(:)
20c4762a1bSJed Brown      PetscInt, target, dimension(1) ::  bcField
21c4762a1bSJed Brown      PetscInt, pointer :: pBcField(:)
22c4762a1bSJed Brown      PetscInt, parameter :: zero = 0, one = 1, two = 2, eight = 8
232e3d3ef9SMatthew G. Knepley      PetscMPIInt :: size
24c4762a1bSJed Brown      IS, target, dimension(1) ::   bcCompIS
25c4762a1bSJed Brown      IS, target, dimension(1) ::   bcPointIS
26c4762a1bSJed Brown      IS, pointer :: pBcCompIS(:)
27c4762a1bSJed Brown      IS, pointer :: pBcPointIS(:)
28c4762a1bSJed Brown      PetscErrorCode :: ierr
29c4762a1bSJed Brown
30d8606c27SBarry Smith      PetscCallA(PetscInitialize(ierr))
31d8606c27SBarry Smith      PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD, size, ierr))
32c4762a1bSJed Brown!     Create a mesh
33d8606c27SBarry Smith      PetscCallA(DMCreate(PETSC_COMM_WORLD, dm, ierr))
34d8606c27SBarry Smith      PetscCallA(DMSetType(dm, DMPLEX, ierr))
35d8606c27SBarry Smith      PetscCallA(DMSetFromOptions(dm, ierr))
36d8606c27SBarry Smith      PetscCallA(DMViewFromOptions(dm, PETSC_NULL_VEC, '-dm_view', ierr))
37d8606c27SBarry Smith      PetscCallA(DMGetDimension(dm, dim, ierr))
38c4762a1bSJed Brown!     Create a scalar field u, a vector field v, and a surface vector field w
39c4762a1bSJed Brown      numFields  = 3
40c4762a1bSJed Brown      numComp(1) = 1
41c4762a1bSJed Brown      numComp(2) = dim
42c4762a1bSJed Brown      numComp(3) = dim-1
43c4762a1bSJed Brown      pNumComp => numComp
44c4762a1bSJed Brown      do i = 1, numFields*(dim+1)
45c4762a1bSJed Brown         numDof(i) = 0
46c4762a1bSJed Brown      end do
47c4762a1bSJed Brown!     Let u be defined on vertices
48c4762a1bSJed Brown      numDof(0*(dim+1)+1)     = 1
49c4762a1bSJed Brown!     Let v be defined on cells
50c4762a1bSJed Brown      numDof(1*(dim+1)+dim+1) = dim
51c4762a1bSJed Brown!     Let v be defined on faces
52c4762a1bSJed Brown      numDof(2*(dim+1)+dim)   = dim-1
53c4762a1bSJed Brown      pNumDof => numDof
54c4762a1bSJed Brown!     Setup boundary conditions
55c4762a1bSJed Brown      numBC = 1
56c4762a1bSJed Brown!     Test label retrieval
57d8606c27SBarry Smith      PetscCallA(DMGetLabel(dm, 'marker', label, ierr))
58d8606c27SBarry Smith      PetscCallA(DMLabelGetValue(label, zero, val, ierr))
59*dcb3e689SBarry Smith      PetscCheckA(size .ne. 1 .or. val .eq. -1,PETSC_COMM_SELF,PETSC_ERR_PLIB,'Error in library')
60d8606c27SBarry Smith      PetscCallA(DMLabelGetValue(label, eight, val, ierr))
61*dcb3e689SBarry Smith      PetscCheckA(size .ne. 1 .or. val .eq. 1,PETSC_COMM_SELF,PETSC_ERR_PLIB,'Error in library')
62c4762a1bSJed Brown!     Prescribe a Dirichlet condition on u on the boundary
63c4762a1bSJed Brown!       Label "marker" is made by the mesh creation routine
64c4762a1bSJed Brown      bcField(1) = 0
65c4762a1bSJed Brown      pBcField => bcField
66d8606c27SBarry Smith      PetscCallA(ISCreateStride(PETSC_COMM_WORLD, one, zero, one, bcCompIS(1), ierr))
67c4762a1bSJed Brown      pBcCompIS => bcCompIS
68d8606c27SBarry Smith      PetscCallA(DMGetStratumIS(dm, 'marker', one, bcPointIS(1),ierr))
69c4762a1bSJed Brown      pBcPointIS => bcPointIS
70c4762a1bSJed Brown!     Create a PetscSection with this data layout
71d8606c27SBarry Smith      PetscCallA(DMSetNumFields(dm, numFields,ierr))
72d8606c27SBarry Smith      PetscCallA(DMPlexCreateSection(dm,nolabel,pNumComp,pNumDof,numBC,pBcField,pBcCompIS,pBcPointIS,PETSC_NULL_IS,section,ierr))
73d8606c27SBarry Smith      PetscCallA(ISDestroy(bcCompIS(1), ierr))
74d8606c27SBarry Smith      PetscCallA(ISDestroy(bcPointIS(1), ierr))
75c4762a1bSJed Brown!     Name the Field variables
76d8606c27SBarry Smith      PetscCallA(PetscSectionSetFieldName(section, zero, 'u', ierr))
77d8606c27SBarry Smith      PetscCallA(PetscSectionSetFieldName(section, one,  'v', ierr))
78d8606c27SBarry Smith      PetscCallA(PetscSectionSetFieldName(section, two,  'w', ierr))
792e3d3ef9SMatthew G. Knepley      if (size .eq. 1) then
80d8606c27SBarry Smith        PetscCallA(PetscSectionView(section, PETSC_VIEWER_STDOUT_WORLD, ierr))
812e3d3ef9SMatthew G. Knepley      endif
82c4762a1bSJed Brown!     Tell the DM to use this data layout
83d8606c27SBarry Smith      PetscCallA(DMSetLocalSection(dm, section, ierr))
84c4762a1bSJed Brown!     Create a Vec with this layout and view it
85d8606c27SBarry Smith      PetscCallA(DMGetGlobalVector(dm, u, ierr))
86d8606c27SBarry Smith      PetscCallA(PetscViewerCreate(PETSC_COMM_WORLD, viewer, ierr))
87d8606c27SBarry Smith      PetscCallA(PetscViewerSetType(viewer, PETSCVIEWERVTK, ierr))
88d8606c27SBarry Smith      PetscCallA(PetscViewerFileSetName(viewer, 'sol.vtu', ierr))
89d8606c27SBarry Smith      PetscCallA(VecView(u, viewer, ierr))
90d8606c27SBarry Smith      PetscCallA(PetscViewerDestroy(viewer, ierr))
91d8606c27SBarry Smith      PetscCallA(DMRestoreGlobalVector(dm, u, ierr))
92c4762a1bSJed Brown!     Cleanup
93d8606c27SBarry Smith      PetscCallA(PetscSectionDestroy(section, ierr))
94d8606c27SBarry Smith      PetscCallA(DMDestroy(dm, ierr))
95c4762a1bSJed Brown
96d8606c27SBarry Smith      PetscCallA(PetscFinalize(ierr))
97c4762a1bSJed Brown      end program DMPlexTestField
98c4762a1bSJed Brown
99c4762a1bSJed Brown!/*TEST
100c4762a1bSJed Brown!  build:
101dfd57a17SPierre Jolivet!    requires: defined(PETSC_USING_F90FREEFORM)
102c4762a1bSJed Brown!
103c4762a1bSJed Brown!  test:
104c4762a1bSJed Brown!    suffix: 0
105c4762a1bSJed Brown!    requires: triangle
10630602db0SMatthew G. Knepley!    args: -info :~sys,mat:
107c4762a1bSJed Brown!
108c4762a1bSJed Brown!  test:
1092e3d3ef9SMatthew G. Knepley!    suffix: 0_2
1102e3d3ef9SMatthew G. Knepley!    requires: triangle
1112e3d3ef9SMatthew G. Knepley!    nsize: 2
1122e3d3ef9SMatthew G. Knepley!    args: -info :~sys,mat,partitioner:
1132e3d3ef9SMatthew G. Knepley!
1142e3d3ef9SMatthew G. Knepley!  test:
115c4762a1bSJed Brown!    suffix: 1
116c4762a1bSJed Brown!    requires: ctetgen
11730602db0SMatthew G. Knepley!    args: -dm_plex_dim 3 -info :~sys,mat:
118c4762a1bSJed Brown!
119c4762a1bSJed Brown!TEST*/
120