Actual source code: ex11f90.F90

      1: !     Tests DMDAGetVecGetArray()

      3: program main
      4: #include <petsc/finclude/petscdm.h>
      5: #include <petsc/finclude/petscdmda.h>
      6:   use petscdmda
      7:   use petsc
      8:   implicit none

     10:   Type(tVec) g
     11:   Type(tDM) ada

     13:   PetscScalar, pointer :: x1(:), x2(:, :)
     14:   PetscScalar, pointer :: x3(:, :, :), x4(:, :, :, :)
     15:   PetscErrorCode ierr
     16:   PetscInt m, n, p, dof, s, i, j, k, xs, xl
     17:   PetscInt ys, yl
     18:   PetscInt zs, zl, sw

     20:   PetscInt nen, nel
     21:   PetscInt, pointer :: elements(:)

     23:   PetscInt nfields
     24:   character(80), pointer :: namefields(:)
     25:   IS, pointer :: isfields(:)
     26:   DM, pointer :: dmfields(:)
     27:   PetscInt zero, one

     29:   m = 5
     30:   n = 6
     31:   p = 4
     32:   s = 1
     33:   dof = 1
     34:   sw = 1
     35:   zero = 0
     36:   one = 1
     37:   PetscCallA(PetscInitialize(ierr))
     38:   PetscCallA(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, m, dof, sw, PETSC_NULL_INTEGER_ARRAY, ada, ierr))
     39:   PetscCallA(DMSetUp(ada, ierr))
     40:   PetscCallA(DMGetGlobalVector(ada, g, ierr))
     41:   PetscCallA(DMDAGetCorners(ada, xs, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, xl, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, ierr))
     42:   PetscCallA(DMDAVecGetArray(ada, g, x1, ierr))
     43:   do i = xs, xs + xl - 1
     44: !         CHKMEMQ
     45:     x1(i) = i
     46: !         CHKMEMQ
     47:   end do
     48:   PetscCallA(DMDAVecRestoreArray(ada, g, x1, ierr))
     49:   PetscCallA(VecView(g, PETSC_VIEWER_STDOUT_WORLD, ierr))
     50:   PetscCallA(DMRestoreGlobalVector(ada, g, ierr))
     51:   PetscCallA(DMDestroy(ada, ierr))

     53:   PetscCallA(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, m, n, PETSC_DECIDE, PETSC_DECIDE, dof, s, PETSC_NULL_INTEGER_ARRAY, PETSC_NULL_INTEGER_ARRAY, ada, ierr))
     54:   PetscCallA(DMSetUp(ada, ierr))
     55:   PetscCallA(DMGetGlobalVector(ada, g, ierr))
     56:   PetscCallA(DMDAGetCorners(ada, xs, ys, PETSC_NULL_INTEGER, xl, yl, PETSC_NULL_INTEGER, ierr))
     57:   PetscCallA(DMDAVecGetArray(ada, g, x2, ierr))
     58:   do i = xs, xs + xl - 1
     59:     do j = ys, ys + yl - 1
     60: !           CHKMEMQ
     61:       x2(i, j) = i + j
     62: !           CHKMEMQ
     63:     end do
     64:   end do
     65:   PetscCallA(DMDAVecRestoreArray(ada, g, x2, ierr))
     66:   PetscCallA(VecView(g, PETSC_VIEWER_STDOUT_WORLD, ierr))
     67:   PetscCallA(DMRestoreGlobalVector(ada, g, ierr))

     69:   PetscCallA(DMDAGetElements(ada, nen, nel, elements, ierr))
     70:   do i = 1, nen*nel
     71:     PetscCheckA(elements(i) >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, 'Error getting DMDA elements')
     72:   end do
     73:   PetscCallA(DMDARestoreElements(ada, nen, nel, elements, ierr))
     74:   PetscCallA(DMDestroy(ada, ierr))

     76:   PetscCallA(DMDACreate3d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, m, n, p, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, dof, s, PETSC_NULL_INTEGER_ARRAY, PETSC_NULL_INTEGER_ARRAY, PETSC_NULL_INTEGER_ARRAY, ada, ierr))
     77:   PetscCallA(DMSetUp(ada, ierr))
     78:   PetscCallA(DMGetGlobalVector(ada, g, ierr))
     79:   PetscCallA(DMDAGetCorners(ada, xs, ys, zs, xl, yl, zl, ierr))
     80:   PetscCallA(DMDAVecGetArray(ada, g, x3, ierr))
     81:   do i = xs, xs + xl - 1
     82:     do j = ys, ys + yl - 1
     83:       do k = zs, zs + zl - 1
     84: !            CHKMEMQ
     85:         x3(i, j, k) = i + j + k
     86: !            CHKMEMQ
     87:       end do
     88:     end do
     89:   end do
     90:   PetscCallA(DMDAVecRestoreArray(ada, g, x3, ierr))
     91:   PetscCallA(VecView(g, PETSC_VIEWER_STDOUT_WORLD, ierr))
     92:   PetscCallA(DMRestoreGlobalVector(ada, g, ierr))
     93:   PetscCallA(DMDestroy(ada, ierr))

     95: !
     96: !  Same tests but now with DOF > 1, so dimensions of array are one higher
     97: !
     98:   dof = 2
     99:   PetscCallA(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, m, dof, sw, PETSC_NULL_INTEGER_ARRAY, ada, ierr))
    100:   PetscCallA(DMSetUp(ada, ierr))
    101:   PetscCallA(DMGetGlobalVector(ada, g, ierr))
    102:   PetscCallA(DMDAGetCorners(ada, xs, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, xl, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, ierr))
    103:   PetscCallA(DMDAVecGetArray(ada, g, x2, ierr))
    104:   do i = xs, xs + xl - 1
    105: !         CHKMEMQ
    106:     x2(0, i) = i
    107:     x2(1, i) = -i
    108: !         CHKMEMQ
    109:   end do
    110:   PetscCallA(DMDAVecRestoreArray(ada, g, x1, ierr))
    111:   PetscCallA(VecView(g, PETSC_VIEWER_STDOUT_WORLD, ierr))
    112:   PetscCallA(DMRestoreGlobalVector(ada, g, ierr))

    114:   ! some testing unrelated to the example
    115:   PetscCallA(DMDASetFieldName(ada, zero, 'Field 0', ierr))
    116:   PetscCallA(DMDASetFieldName(ada, one, 'Field 1', ierr))
    117:   PetscCallA(DMCreateFieldDecomposition(ada, nfields, namefields, PETSC_NULL_IS_POINTER, PETSC_NULL_DM_POINTER, ierr))
    118:   ! print*,nfields,trim(namefields(1)),trim(namefields(2))
    119:   PetscCallA(DMDestroyFieldDecomposition(ada, nfields, namefields, PETSC_NULL_IS_POINTER, PETSC_NULL_DM_POINTER, ierr))
    120:   PetscCallA(DMCreateFieldDecomposition(ada, nfields, namefields, isfields, dmfields, ierr))
    121:   PetscCallA(DMDestroyFieldDecomposition(ada, nfields, namefields, isfields, dmfields, ierr))

    123:   PetscCallA(DMDestroy(ada, ierr))

    125:   dof = 2
    126:   PetscCallA(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, m, n, PETSC_DECIDE, PETSC_DECIDE, dof, s, PETSC_NULL_INTEGER_ARRAY, PETSC_NULL_INTEGER_ARRAY, ada, ierr))
    127:   PetscCallA(DMSetUp(ada, ierr))
    128:   PetscCallA(DMGetGlobalVector(ada, g, ierr))
    129:   PetscCallA(DMDAGetCorners(ada, xs, ys, PETSC_NULL_INTEGER, xl, yl, PETSC_NULL_INTEGER, ierr))
    130:   PetscCallA(DMDAVecGetArray(ada, g, x3, ierr))
    131:   do i = xs, xs + xl - 1
    132:     do j = ys, ys + yl - 1
    133: !           CHKMEMQ
    134:       x3(0, i, j) = i + j
    135:       x3(1, i, j) = -(i + j)
    136: !           CHKMEMQ
    137:     end do
    138:   end do
    139:   PetscCallA(DMDAVecRestoreArray(ada, g, x3, ierr))
    140:   PetscCallA(VecView(g, PETSC_VIEWER_STDOUT_WORLD, ierr))
    141:   PetscCallA(DMRestoreGlobalVector(ada, g, ierr))
    142:   PetscCallA(DMDestroy(ada, ierr))

    144:   dof = 3
    145:   PetscCallA(DMDACreate3d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, m, n, p, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, dof, s, PETSC_NULL_INTEGER_ARRAY, PETSC_NULL_INTEGER_ARRAY, PETSC_NULL_INTEGER_ARRAY, ada, ierr))
    146:   PetscCallA(DMSetUp(ada, ierr))
    147:   PetscCallA(DMGetGlobalVector(ada, g, ierr))
    148:   PetscCallA(DMDAGetCorners(ada, xs, ys, zs, xl, yl, zl, ierr))
    149:   PetscCallA(DMDAVecGetArray(ada, g, x4, ierr))
    150:   do i = xs, xs + xl - 1
    151:     do j = ys, ys + yl - 1
    152:       do k = zs, zs + zl - 1
    153: !            CHKMEMQ
    154:         x4(0, i, j, k) = i + j + k
    155:         x4(1, i, j, k) = -(i + j + k)
    156:         x4(2, i, j, k) = i + j + k
    157: !            CHKMEMQ
    158:       end do
    159:     end do
    160:   end do
    161:   PetscCallA(DMDAVecRestoreArray(ada, g, x4, ierr))
    162:   PetscCallA(VecView(g, PETSC_VIEWER_STDOUT_WORLD, ierr))
    163:   PetscCallA(DMRestoreGlobalVector(ada, g, ierr))
    164:   PetscCallA(DMDestroy(ada, ierr))

    166:   PetscCallA(PetscFinalize(ierr))
    167: END PROGRAM

    169: !
    170: !/*TEST
    171: !
    172: !   build:
    173: !     requires: !complex
    174: !
    175: !   test:
    176: !     filter: Error: grep -v "Vec Object" | grep -v "Warning: ieee_inexact is signaling"
    177: !
    178: !TEST*/