Merge branch '411_use_petsc_quadrature' into 'development'

Use PETSc quadrature rules

See merge request damask/DAMASK!982

src/homogenization.f90

@@ -185,14 +185,10 @@ end subroutine homogenization_set_phi
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief module initialization
+!> @brief Module initialization.
 !--------------------------------------------------------------------------------------------------
 subroutine homogenization_init()
 
- type(tDict) , pointer :: &
- num_homog, &
- num_homogGeneric
-
  print'(/,1x,a)', '<<<+- homogenization init -+>>>'; flush(IO_STDOUT)
 
 

src/materialpoint.f90

@@ -27,7 +27,10 @@ module materialpoint
  use homogenization
  use discretization
 #if defined(MESH)
+#include "petscversion.h"
+#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR<18)
  use FEM_quadrature
+#endif
  use discretization_mesh
 #elif defined(GRID)
  use base64
@@ -52,7 +55,7 @@ subroutine materialpoint_initAll()
  call prec_init()
  call misc_init()
  call IO_init()
-#if defined(MESH)
+#if defined(MESH) && (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR<18)
  call FEM_quadrature_init()
 #elif defined(GRID)
  call base64_init()

src/mesh/FEM_utilities.f90

@@ -17,7 +17,9 @@ module FEM_utilities
  use parallelization
  use discretization_mesh
  use homogenization
+#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR<18)
  use FEM_quadrature
+#endif
  use constants
 
 #if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY)
@@ -89,7 +91,11 @@ subroutine FEM_utilities_init(num_mesh)
  p_s = num_mesh%get_asInt('p_s',defaultVal = 2)
  p_i = num_mesh%get_asInt('p_i',defaultVal = p_s)
 
+#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR<18)
  if (p_s < 1 .or. p_s > size(FEM_nQuadrature,2)) &
+#else
+ if (p_s < 1) &
+#endif
  call IO_error(821,ext_msg='shape function order (p_s) out of bounds')
  if (p_i < max(1,p_s-1) .or. p_i > p_s) &
  call IO_error(821,ext_msg='integration order (p_i) out of bounds')

src/mesh/discretization_mesh.f90

@@ -11,6 +11,9 @@ module discretization_mesh
  use PETScDMplex
  use PETScDMDA
  use PETScIS
+#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>=18)
+ use PETScDT
+#endif
 #if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY)
  use MPI_f08
 #endif
@@ -20,7 +23,9 @@ module discretization_mesh
  use config
  use discretization
  use result
+#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR<18)
  use FEM_quadrature
+#endif
  use types
  use prec
 
@@ -80,6 +85,9 @@ subroutine discretization_mesh_init()
  PetscInt :: nFaceSets, Nboundaries, NelemsGlobal, Nelems
  PetscInt, pointer, dimension(:) :: pFaceSets
  IS :: faceSetIS
+#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>=18)
+ PetscQuadrature :: quadrature
+#endif
  PetscErrorCode :: err_PETSc
  integer(MPI_INTEGER_KIND) :: err_MPI
  PetscInt, dimension(:), allocatable :: &
@@ -91,6 +99,11 @@ subroutine discretization_mesh_init()
  type(tvec) :: coords_node0
  real(pREAL), pointer, dimension(:) :: &
  mesh_node0_temp
+#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>=18)
+ real(pREAL), pointer, dimension(:) :: &
+ qPointsP, &
+ PETSC_NULL_REAL_PTR => null()
+#endif
 
 
  print'(/,1x,a)', '<<<+- discretization_mesh init -+>>>'
@@ -158,10 +171,27 @@ subroutine discretization_mesh_init()
  CHKERRQ(err_PETSc)
 
 ! Get initial nodal coordinates
+#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR<18)
  mesh_maxNips = FEM_nQuadrature(dimPlex,p_i)
 
  call mesh_FEM_build_ipCoordinates(dimPlex,FEM_quadrature_points(dimPlex,p_i)%p)
  call mesh_FEM_build_ipVolumes(dimPlex)
+#else
+ call PetscDTSimplexQuadrature(dimplex, p_i, -1, quadrature, err_PETSc)
+ CHKERRQ(err_PETSc)
+ call PetscQuadratureGetData(quadrature,PETSC_NULL_INTEGER(1),PETSC_NULL_INTEGER(1), &
+ mesh_maxNips,qPointsP,PETSC_NULL_REAL_PTR,err_PETSc)
+ CHKERRQ(err_PETSc)
+
+ call mesh_FEM_build_ipCoordinates(dimPlex,qPointsP)
+ call mesh_FEM_build_ipVolumes(dimPlex)
+
+ call PetscQuadratureRestoreData(quadrature,PETSC_NULL_INTEGER(1),PETSC_NULL_INTEGER(1), &
+ PETSC_NULL_INTEGER(1),qPointsP,PETSC_NULL_REAL_PTR,err_PETSc)
+ CHKERRQ(err_PETSc)
+ call PetscQuadratureDestroy(quadrature, err_PETSc)
+ CHKERRQ(err_PETSc)
+#endif
 
  allocate(materialAt(mesh_NcpElems))
  do j = 1, mesh_NcpElems

src/mesh/mesh_mech_FEM.f90

@@ -20,7 +20,9 @@ module mesh_mechanical_FEM
  use FEM_utilities
  use discretization
  use discretization_mesh
+#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR<18)
  use FEM_quadrature
+#endif
  use homogenization
  use math
  use constants
@@ -75,7 +77,7 @@ module mesh_mechanical_FEM
 #if defined(PETSC_USE_64BIT_INDICES) || PETSC_VERSION_MINOR < 16
  ISDestroy, &
 #endif
-#if PETSC_VERSION_MINOR > 18
+#if PETSC_VERSION_MINOR > 18 && PETSC_VERSION_MINOR < 22
  DMAddField, &
 #endif
  PetscSectionGetNumFields, &
@@ -120,6 +122,9 @@ subroutine FEM_mechanical_init(mechBC,num_mesh)
  PetscSection :: section
 
  PetscReal, dimension(:), pointer :: qPointsP, qWeightsP, &
+#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>=18)
+ PETSC_NULL_REAL_PTR => null(), &
+#endif
  nodalPointsP, nodalWeightsP,pV0, pCellJ, pInvcellJ
  PetscReal :: detJ
  PetscReal, allocatable, target :: cellJMat(:,:)
@@ -160,6 +165,7 @@ subroutine FEM_mechanical_init(mechBC,num_mesh)
 
 !--------------------------------------------------------------------------------------------------
 ! Setup FEM mech discretization
+#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR<18)
  qPoints = FEM_quadrature_points( dimPlex,num%p_i)%p
  qWeights = FEM_quadrature_weights(dimPlex,num%p_i)%p
  nQuadrature = FEM_nQuadrature( dimPlex,num%p_i)
@@ -170,6 +176,19 @@ subroutine FEM_mechanical_init(mechBC,num_mesh)
  nc = dimPlex
  call PetscQuadratureSetData(mechQuad,dimPlex,nc,int(nQuadrature,pPETSCINT),qPointsP,qWeightsP,err_PETSc)
  CHKERRQ(err_PETSc)
+#else
+ call PetscDTSimplexQuadrature(dimplex,num%p_i,-1,mechQuad,err_PETSc)
+ CHKERRQ(err_PETSc)
+ call PetscQuadratureGetData(mechQuad,PETSC_NULL_INTEGER(1),PETSC_NULL_INTEGER(1), &
+ nQuadrature,PETSC_NULL_REAL_PTR,qWeightsP,err_PETSc)
+ CHKERRQ(err_PETSc)
+ qWeights = qWeightsP
+ call PetscQuadratureRestoreData(mechQuad,PETSC_NULL_INTEGER(1),PETSC_NULL_INTEGER(1), &
+ PETSC_NULL_INTEGER(1),PETSC_NULL_REAL_PTR,qWeightsP, &
+ err_PETSc)
+ CHKERRQ(err_PETSc)
+ nc = dimPlex
+#endif
  call PetscFECreateDefault(PETSC_COMM_SELF,dimPlex,nc,PETSC_TRUE,prefix, &
  num%p_i,mechFE,err_PETSc)
  CHKERRQ(err_PETSc)
@@ -787,14 +806,14 @@ subroutine FEM_mechanical_updateCoords()
  PetscInt :: pStart, pEnd, p, s, e, q, &
  cellStart, cellEnd, c, n
  PetscSection :: section
- PetscQuadrature :: mechQuad
+ PetscDS :: mechDS
  PetscReal, dimension(:), pointer :: basisField, dev_null, &
  nodeCoords_linear !< nodal coordinates (dimPlex*Nnodes)
  real(pREAL), dimension(:), pointer :: x_scal
 
  call SNESGetDM(mechanical_snes,dm_local,err_PETSc)
  CHKERRQ(err_PETSc)
- call DMGetDS(dm_local,mechQuad,err_PETSc)
+ call DMGetDS(dm_local,mechDS,err_PETSc)
  CHKERRQ(err_PETSc)
  call DMGetLocalSection(dm_local,section,err_PETSc)
  CHKERRQ(err_PETSc)
@@ -822,7 +841,7 @@ subroutine FEM_mechanical_updateCoords()
  ! write ip displacements
  call DMPlexGetHeightStratum(dm_local,0_pPETSCINT,cellStart,cellEnd,err_PETSc)
  CHKERRQ(err_PETSc)
- call PetscDSGetTabulation(mechQuad,0_pPETSCINT,basisField,dev_null,err_PETSc)
+ call PetscDSGetTabulation(mechDS,0_pPETSCINT,basisField,dev_null,err_PETSc)
  CHKERRQ(err_PETSc)
  allocate(ipCoords(3,nQuadrature,mesh_NcpElems),source=0.0_pREAL)
  do c=cellStart,cellEnd-1_pPETSCINT
@@ -849,7 +868,7 @@ subroutine FEM_mechanical_updateCoords()
  call discretization_setIPcoords(reshape(ipCoords,[3,mesh_NcpElems*nQuadrature]))
  call DMRestoreLocalVector(dm_local,x_local,err_PETSc)
  CHKERRQ(err_PETSc)
- call PetscDSRestoreTabulation(mechQuad,0_pPETSCINT,basisField,dev_null,err_PETSc)
+ call PetscDSRestoreTabulation(mechDS,0_pPETSCINT,basisField,dev_null,err_PETSc)
  CHKERRQ(err_PETSc)
 
 end subroutine FEM_mechanical_updateCoords