diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 937a327bd..e06742282 100644 --- a/src/homogenization.f90 +++ b/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) diff --git a/src/materialpoint.f90 b/src/materialpoint.f90 index bbc970e95..fd9053aae 100644 --- a/src/materialpoint.f90 +++ b/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() diff --git a/src/mesh/FEM_utilities.f90 b/src/mesh/FEM_utilities.f90 index 3f7d0b4ad..beb0a276d 100644 --- a/src/mesh/FEM_utilities.f90 +++ b/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') diff --git a/src/mesh/discretization_mesh.f90 b/src/mesh/discretization_mesh.f90 index 67656fd1b..ed48e6f76 100644 --- a/src/mesh/discretization_mesh.f90 +++ b/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 diff --git a/src/mesh/mesh_mech_FEM.f90 b/src/mesh/mesh_mech_FEM.f90 index 1220c9196..aee11d1e5 100644 --- a/src/mesh/mesh_mech_FEM.f90 +++ b/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