223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

mesh: separate kind for DAMASK and PETSc integers

This commit is contained in:
Martin Diehl 2022-01-13 16:01:06 +01:00
parent 29530da579
commit 82dabe29c1
4 changed files with 62 additions and 47 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
.gitlab-ci.yml
View File

@ -100,6 +100,13 @@ test_grid_GNU-64bit:
- cd PRIVATE/testing/pytest - cd PRIVATE/testing/pytest
- pytest -k 'compile and grid' --basetemp ${TESTROOT}/compile_grid_GNU-64bit - pytest -k 'compile and grid' --basetemp ${TESTROOT}/compile_grid_GNU-64bit
test_mesh_GNU-64bit:
stage: compile
script:
- module load Compiler/GNU/10 Libraries/PETSc/3.16.2/64bit
- cd PRIVATE/testing/pytest
- pytest -k 'compile and mesh' --basetemp ${TESTROOT}/compile_mesh_GNU-64bit
test_grid_IntelLLVM: test_grid_IntelLLVM:
stage: compile stage: compile
script: script:

2
src/mesh/FEM_utilities.f90
View File

@ -50,7 +50,7 @@ module FEM_utilities
type, public :: tSolutionState !< return type of solution from FEM solver variants type, public :: tSolutionState !< return type of solution from FEM solver variants
logical :: converged = .true. logical :: converged = .true.
logical :: stagConverged = .true. logical :: stagConverged = .true.
integer :: iterationsNeeded = 0 PetscInt :: iterationsNeeded = 0_pPETSCINT
end type tSolutionState end type tSolutionState
type, public :: tComponentBC type, public :: tComponentBC

28
src/mesh/discretization_mesh.f90
View File

@ -71,22 +71,22 @@ subroutine discretization_mesh_init(restart)
logical, intent(in) :: restart logical, intent(in) :: restart
integer :: dimPlex, & PetscInt :: dimPlex, &
mesh_Nnodes, & !< total number of nodes in mesh mesh_Nnodes, & !< total number of nodes in mesh
j, & j, &
debug_element, debug_ip debug_element, debug_ip
PetscSF :: sf PetscSF :: sf
DM :: globalMesh DM :: globalMesh
PetscInt :: nFaceSets PetscInt :: nFaceSets, Nboundaries, NelemsGlobal, Nelems
PetscInt, pointer, dimension(:) :: pFaceSets PetscInt, pointer, dimension(:) :: pFaceSets
IS :: faceSetIS IS :: faceSetIS
PetscErrorCode :: err_PETSc PetscErrorCode :: err_PETSc
integer(MPI_INTEGER_KIND) :: err_MPI integer(MPI_INTEGER_KIND) :: err_MPI
integer, dimension(:), allocatable :: & PetscInt, dimension(:), allocatable :: &
materialAt materialAt
class(tNode), pointer :: & class(tNode), pointer :: &
num_mesh num_mesh
integer :: p_i !< integration order (quadrature rule) integer :: p_i, dim !< integration order (quadrature rule)
type(tvec) :: coords_node0 type(tvec) :: coords_node0
print'(/,1x,a)', '<<<+- discretization_mesh init -+>>>' print'(/,1x,a)', '<<<+- discretization_mesh init -+>>>'
@ -105,20 +105,23 @@ subroutine discretization_mesh_init(restart)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call DMGetDimension(globalMesh,dimPlex,err_PETSc) call DMGetDimension(globalMesh,dimPlex,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call DMGetStratumSize(globalMesh,'depth',dimPlex,mesh_NcpElemsGlobal,err_PETSc) call DMGetStratumSize(globalMesh,'depth',dimPlex,NelemsGlobal,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
print'()' mesh_NcpElemsGlobal = int(NelemsGlobal)
call DMView(globalMesh, PETSC_VIEWER_STDOUT_WORLD,err_PETSc) call DMView(globalMesh, PETSC_VIEWER_STDOUT_WORLD,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
! get number of IDs in face sets (for boundary conditions?) ! get number of IDs in face sets (for boundary conditions?)
call DMGetLabelSize(globalMesh,'Face Sets',mesh_Nboundaries,err_PETSc) call DMGetLabelSize(globalMesh,'Face Sets',Nboundaries,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
mesh_Nboundaries = int(Nboundaries)
call MPI_Bcast(mesh_Nboundaries,1_MPI_INTEGER_KIND,MPI_INTEGER,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI) call MPI_Bcast(mesh_Nboundaries,1_MPI_INTEGER_KIND,MPI_INTEGER,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
call MPI_Bcast(mesh_NcpElemsGlobal,1_MPI_INTEGER_KIND,MPI_INTEGER,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI) call MPI_Bcast(mesh_NcpElemsGlobal,1_MPI_INTEGER_KIND,MPI_INTEGER,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
call MPI_Bcast(dimPlex,1_MPI_INTEGER_KIND,MPI_INTEGER,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI) dim = int(dimPlex)
call MPI_Bcast(dim,1_MPI_INTEGER_KIND,MPI_INTEGER,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI)
dimPlex = int(dim,pPETSCINT)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
if (worldrank == 0) then if (worldrank == 0) then
@ -128,7 +131,7 @@ subroutine discretization_mesh_init(restart)
endif endif
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
allocate(mesh_boundaries(mesh_Nboundaries), source = 0) allocate(mesh_boundaries(mesh_Nboundaries), source = 0_pPETSCINT)
call DMGetLabelSize(globalMesh,'Face Sets',nFaceSets,err_PETSc) call DMGetLabelSize(globalMesh,'Face Sets',nFaceSets,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call DMGetLabelIdIS(globalMesh,'Face Sets',faceSetIS,err_PETSc) call DMGetLabelIdIS(globalMesh,'Face Sets',faceSetIS,err_PETSc)
@ -145,8 +148,9 @@ subroutine discretization_mesh_init(restart)
call DMDestroy(globalMesh,err_PETSc); CHKERRQ(err_PETSc) call DMDestroy(globalMesh,err_PETSc); CHKERRQ(err_PETSc)
call DMGetStratumSize(geomMesh,'depth',dimPlex,mesh_NcpElems,err_PETSc) call DMGetStratumSize(geomMesh,'depth',dimPlex,Nelems,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
mesh_NcpElems = int(Nelems)
call DMGetStratumSize(geomMesh,'depth',0_pPETSCINT,mesh_Nnodes,err_PETSc) call DMGetStratumSize(geomMesh,'depth',0_pPETSCINT,mesh_Nnodes,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
@ -166,7 +170,7 @@ subroutine discretization_mesh_init(restart)
call DMGetLabelValue(geomMesh,'Cell Sets',j-1,materialAt(j),err_PETSc) call DMGetLabelValue(geomMesh,'Cell Sets',j-1,materialAt(j),err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
enddo enddo
materialAt = materialAt + 1 materialAt = materialAt + 1_pPETSCINT
if (debug_element < 1 .or. debug_element > mesh_NcpElems) call IO_error(602,ext_msg='element') if (debug_element < 1 .or. debug_element > mesh_NcpElems) call IO_error(602,ext_msg='element')
if (debug_ip < 1 .or. debug_ip > mesh_maxNips) call IO_error(602,ext_msg='IP') if (debug_ip < 1 .or. debug_ip > mesh_maxNips) call IO_error(602,ext_msg='IP')
@ -175,7 +179,7 @@ subroutine discretization_mesh_init(restart)
mesh_node0(1:dimPlex,:) = reshape(mesh_node0_temp,[dimPlex,mesh_Nnodes]) mesh_node0(1:dimPlex,:) = reshape(mesh_node0_temp,[dimPlex,mesh_Nnodes])
call discretization_init(materialAt,& call discretization_init(int(materialAt),&
reshape(mesh_ipCoordinates,[3,mesh_maxNips*mesh_NcpElems]), & reshape(mesh_ipCoordinates,[3,mesh_maxNips*mesh_NcpElems]), &
mesh_node0) mesh_node0)

72
src/mesh/mesh_mech_FEM.f90
View File

@ -40,7 +40,7 @@ module mesh_mechanical_FEM
type(tSolutionParams) :: params type(tSolutionParams) :: params
type, private :: tNumerics type, private :: tNumerics
integer :: & PetscInt :: &
p_i, & !< integration order (quadrature rule) p_i, & !< integration order (quadrature rule)
itmax itmax
logical :: & logical :: &
@ -55,7 +55,8 @@ module mesh_mechanical_FEM
! PETSc data ! PETSc data
SNES :: mechanical_snes SNES :: mechanical_snes
Vec :: solution, solution_rate, solution_local Vec :: solution, solution_rate, solution_local
PetscInt :: dimPlex, cellDof, nQuadrature, nBasis PetscInt :: dimPlex, cellDof, nBasis
integer :: nQuadrature
PetscReal, allocatable, target :: qPoints(:), qWeights(:) PetscReal, allocatable, target :: qPoints(:), qWeights(:)
MatNullSpace :: matnull MatNullSpace :: matnull
@ -104,8 +105,8 @@ subroutine FEM_mechanical_init(fieldBC)
PetscReal :: detJ PetscReal :: detJ
PetscReal, allocatable, target :: cellJMat(:,:) PetscReal, allocatable, target :: cellJMat(:,:)
PetscScalar, pointer :: px_scal(:) PetscScalar, pointer, dimension(:) :: px_scal
PetscScalar, allocatable, target :: x_scal(:) PetscScalar, allocatable, target, dimension(:) :: x_scal
character(len=*), parameter :: prefix = 'mechFE_' character(len=*), parameter :: prefix = 'mechFE_'
PetscErrorCode :: err_PETSc PetscErrorCode :: err_PETSc
@ -118,8 +119,8 @@ subroutine FEM_mechanical_init(fieldBC)
!----------------------------------------------------------------------------- !-----------------------------------------------------------------------------
! read numerical parametes and do sanity checks ! read numerical parametes and do sanity checks
num_mesh => config_numerics%get('mesh',defaultVal=emptyDict) num_mesh => config_numerics%get('mesh',defaultVal=emptyDict)
num%p_i = num_mesh%get_asInt('p_i',defaultVal = 2) num%p_i = int(num_mesh%get_asInt('p_i',defaultVal = 2),pPETSCINT)
num%itmax = num_mesh%get_asInt('itmax',defaultVal=250) num%itmax = int(num_mesh%get_asInt('itmax',defaultVal=250),pPETSCINT)
num%BBarStabilisation = num_mesh%get_asBool('bbarstabilisation',defaultVal = .false.) num%BBarStabilisation = num_mesh%get_asBool('bbarstabilisation',defaultVal = .false.)
num%eps_struct_atol = num_mesh%get_asFloat('eps_struct_atol', defaultVal = 1.0e-10_pReal) num%eps_struct_atol = num_mesh%get_asFloat('eps_struct_atol', defaultVal = 1.0e-10_pReal)
num%eps_struct_rtol = num_mesh%get_asFloat('eps_struct_rtol', defaultVal = 1.0e-4_pReal) num%eps_struct_rtol = num_mesh%get_asFloat('eps_struct_rtol', defaultVal = 1.0e-4_pReal)
@ -143,7 +144,7 @@ subroutine FEM_mechanical_init(fieldBC)
call PetscQuadratureCreate(PETSC_COMM_SELF,mechQuad,err_PETSc) call PetscQuadratureCreate(PETSC_COMM_SELF,mechQuad,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
nc = dimPlex nc = dimPlex
call PetscQuadratureSetData(mechQuad,dimPlex,nc,nQuadrature,qPointsP,qWeightsP,err_PETSc) call PetscQuadratureSetData(mechQuad,dimPlex,nc,int(nQuadrature,pPETSCINT),qPointsP,qWeightsP,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call PetscFECreateDefault(PETSC_COMM_SELF,dimPlex,nc,PETSC_TRUE,prefix, & call PetscFECreateDefault(PETSC_COMM_SELF,dimPlex,nc,PETSC_TRUE,prefix, &
num%p_i,mechFE,err_PETSc); CHKERRQ(err_PETSc) num%p_i,mechFE,err_PETSc); CHKERRQ(err_PETSc)
@ -165,7 +166,7 @@ subroutine FEM_mechanical_init(fieldBC)
call DMPlexLabelComplete(mechanical_mesh,BCLabel,err_PETSc); CHKERRQ(err_PETSc) call DMPlexLabelComplete(mechanical_mesh,BCLabel,err_PETSc); CHKERRQ(err_PETSc)
call DMGetLocalSection(mechanical_mesh,section,err_PETSc); CHKERRQ(err_PETSc) call DMGetLocalSection(mechanical_mesh,section,err_PETSc); CHKERRQ(err_PETSc)
allocate(pnumComp(1), source=dimPlex) allocate(pnumComp(1), source=dimPlex)
allocate(pnumDof(0:dimPlex), source = 0) allocate(pnumDof(0:dimPlex), source = 0_pPETSCINT)
do topologDim = 0, dimPlex do topologDim = 0, dimPlex
call DMPlexGetDepthStratum(mechanical_mesh,topologDim,cellStart,cellEnd,err_PETSc) call DMPlexGetDepthStratum(mechanical_mesh,topologDim,cellStart,cellEnd,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
@ -176,14 +177,14 @@ subroutine FEM_mechanical_init(fieldBC)
do field = 1, dimPlex; do faceSet = 1, mesh_Nboundaries do field = 1, dimPlex; do faceSet = 1, mesh_Nboundaries
if (fieldBC%componentBC(field)%Mask(faceSet)) numBC = numBC + 1 if (fieldBC%componentBC(field)%Mask(faceSet)) numBC = numBC + 1
enddo; enddo enddo; enddo
allocate(pbcField(numBC), source=0) allocate(pbcField(numBC), source=0_pPETSCINT)
allocate(pbcComps(numBC)) allocate(pbcComps(numBC))
allocate(pbcPoints(numBC)) allocate(pbcPoints(numBC))
numBC = 0 numBC = 0
do field = 1, dimPlex; do faceSet = 1, mesh_Nboundaries do field = 1, dimPlex; do faceSet = 1, mesh_Nboundaries
if (fieldBC%componentBC(field)%Mask(faceSet)) then if (fieldBC%componentBC(field)%Mask(faceSet)) then
numBC = numBC + 1 numBC = numBC + 1
call ISCreateGeneral(PETSC_COMM_WORLD,1,[field-1],PETSC_COPY_VALUES,pbcComps(numBC),err_PETSc) call ISCreateGeneral(PETSC_COMM_WORLD,1_pPETSCINT,[field-1],PETSC_COPY_VALUES,pbcComps(numBC),err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call DMGetStratumSize(mechanical_mesh,'Face Sets',mesh_boundaries(faceSet),bcSize,err_PETSc) call DMGetStratumSize(mechanical_mesh,'Face Sets',mesh_boundaries(faceSet),bcSize,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
@ -196,7 +197,7 @@ subroutine FEM_mechanical_init(fieldBC)
call ISRestoreIndicesF90(bcPoint,pBcPoint,err_PETSc); CHKERRQ(err_PETSc) call ISRestoreIndicesF90(bcPoint,pBcPoint,err_PETSc); CHKERRQ(err_PETSc)
call ISDestroy(bcPoint,err_PETSc); CHKERRQ(err_PETSc) call ISDestroy(bcPoint,err_PETSc); CHKERRQ(err_PETSc)
else else
call ISCreateGeneral(PETSC_COMM_WORLD,0,[0],PETSC_COPY_VALUES,pbcPoints(numBC),err_PETSc) call ISCreateGeneral(PETSC_COMM_WORLD,0_pPETSCINT,[0_pPETSCINT],PETSC_COPY_VALUES,pbcPoints(numBC),err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
endif endif
endif endif
@ -226,7 +227,7 @@ subroutine FEM_mechanical_init(fieldBC)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call DMSNESSetJacobianLocal(mechanical_mesh,FEM_mechanical_formJacobian,PETSC_NULL_VEC,err_PETSc) ! function to evaluate stiffness matrix call DMSNESSetJacobianLocal(mechanical_mesh,FEM_mechanical_formJacobian,PETSC_NULL_VEC,err_PETSc) ! function to evaluate stiffness matrix
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call SNESSetMaxLinearSolveFailures(mechanical_snes, huge(1), err_PETSc) ! ignore linear solve failures CHKERRQ(err_PETSc) !< ignore linear solve failures call SNESSetMaxLinearSolveFailures(mechanical_snes, huge(1_pPETSCINT), err_PETSc) ! ignore linear solve failures
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call SNESSetConvergenceTest(mechanical_snes,FEM_mechanical_converged,PETSC_NULL_VEC,PETSC_NULL_FUNCTION,err_PETSc) call SNESSetConvergenceTest(mechanical_snes,FEM_mechanical_converged,PETSC_NULL_VEC,PETSC_NULL_FUNCTION,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
@ -245,10 +246,10 @@ subroutine FEM_mechanical_init(fieldBC)
allocate(pcellJ(dimPlex**2)) allocate(pcellJ(dimPlex**2))
allocate(pinvcellJ(dimPlex**2)) allocate(pinvcellJ(dimPlex**2))
allocate(cellJMat(dimPlex,dimPlex)) allocate(cellJMat(dimPlex,dimPlex))
call PetscDSGetDiscretization(mechDS,0,mechFE,err_PETSc) call PetscDSGetDiscretization(mechDS,0_pPETSCINT,mechFE,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call PetscFEGetDualSpace(mechFE,mechDualSpace,err_PETSc); CHKERRQ(err_PETSc) call PetscFEGetDualSpace(mechFE,mechDualSpace,err_PETSc); CHKERRQ(err_PETSc)
call DMPlexGetHeightStratum(mechanical_mesh,0,cellStart,cellEnd,err_PETSc) call DMPlexGetHeightStratum(mechanical_mesh,0_pPETSCINT,cellStart,cellEnd,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
do cell = cellStart, cellEnd-1 !< loop over all elements do cell = cellStart, cellEnd-1 !< loop over all elements
x_scal = 0.0_pReal x_scal = 0.0_pReal
@ -352,19 +353,21 @@ subroutine FEM_mechanical_formResidual(dm_local,xx_local,f_local,dummy,err_PETSc
call DMGetLocalSection(dm_local,section,err_PETSc); CHKERRQ(err_PETSc) call DMGetLocalSection(dm_local,section,err_PETSc); CHKERRQ(err_PETSc)
call DMGetDS(dm_local,prob,err_PETSc); CHKERRQ(err_PETSc) call DMGetDS(dm_local,prob,err_PETSc); CHKERRQ(err_PETSc)
call PetscDSGetTabulation(prob,0,basisField,basisFieldDer,err_PETSc) call PetscDSGetTabulation(prob,0_pPETSCINT,basisField,basisFieldDer,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call DMPlexGetHeightStratum(dm_local,0,cellStart,cellEnd,err_PETSc) call DMPlexGetHeightStratum(dm_local,0_pPETSCINT,cellStart,cellEnd,err_PETSc)
CHKERRQ(err_PETSc)
call DMGetLocalVector(dm_local,x_local,err_PETSc)
CHKERRQ(err_PETSc)
call VecWAXPY(x_local,1.0,xx_local,solution_local,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call DMGetLocalVector(dm_local,x_local,err_PETSc); CHKERRQ(err_PETSc)
call VecWAXPY(x_local,1.0,xx_local,solution_local,err_PETSc); CHKERRQ(err_PETSc)
do field = 1, dimPlex; do face = 1, mesh_Nboundaries do field = 1, dimPlex; do face = 1, mesh_Nboundaries
if (params%fieldBC%componentBC(field)%Mask(face)) then if (params%fieldBC%componentBC(field)%Mask(face)) then
call DMGetStratumSize(dm_local,'Face Sets',mesh_boundaries(face),bcSize,err_PETSc) call DMGetStratumSize(dm_local,'Face Sets',mesh_boundaries(face),bcSize,err_PETSc)
if (bcSize > 0) then if (bcSize > 0) then
call DMGetStratumIS(dm_local,'Face Sets',mesh_boundaries(face),bcPoints,err_PETSc) call DMGetStratumIS(dm_local,'Face Sets',mesh_boundaries(face),bcPoints,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call utilities_projectBCValues(x_local,section,0,field-1,bcPoints, & call utilities_projectBCValues(x_local,section,0_pPETSCINT,field-1,bcPoints, &
0.0_pReal,params%fieldBC%componentBC(field)%Value(face),params%timeinc) 0.0_pReal,params%fieldBC%componentBC(field)%Value(face),params%timeinc)
call ISDestroy(bcPoints,err_PETSc); CHKERRQ(err_PETSc) call ISDestroy(bcPoints,err_PETSc); CHKERRQ(err_PETSc)
endif endif
@ -515,7 +518,7 @@ subroutine FEM_mechanical_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,err_P
if (bcSize > 0) then if (bcSize > 0) then
call DMGetStratumIS(dm_local,'Face Sets',mesh_boundaries(face),bcPoints,err_PETSc) call DMGetStratumIS(dm_local,'Face Sets',mesh_boundaries(face),bcPoints,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call utilities_projectBCValues(x_local,section,0,field-1,bcPoints, & call utilities_projectBCValues(x_local,section,0_pPETSCINT,field-1,bcPoints, &
0.0_pReal,params%fieldBC%componentBC(field)%Value(face),params%timeinc) 0.0_pReal,params%fieldBC%componentBC(field)%Value(face),params%timeinc)
call ISDestroy(bcPoints,err_PETSc); CHKERRQ(err_PETSc) call ISDestroy(bcPoints,err_PETSc); CHKERRQ(err_PETSc)
endif endif
@ -553,11 +556,11 @@ subroutine FEM_mechanical_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,err_P
FInv = math_inv33(F) FInv = math_inv33(F)
K_eA = K_eA + matmul(transpose(BMat),MatA)*math_det33(FInv)**(1.0/real(dimPlex)) K_eA = K_eA + matmul(transpose(BMat),MatA)*math_det33(FInv)**(1.0/real(dimPlex))
K_eB = K_eB - & K_eB = K_eB - &
matmul(transpose(matmul(reshape(homogenization_F(1:dimPlex,1:dimPlex,m),shape=[dimPlex*dimPlex,1]), & matmul(transpose(matmul(reshape(homogenization_F(1:dimPlex,1:dimPlex,m),shape=[dimPlex**2,1_pPETSCINT]), &
matmul(reshape(FInv(1:dimPlex,1:dimPlex), & matmul(reshape(FInv(1:dimPlex,1:dimPlex), &
shape=[1,dimPlex*dimPlex],order=[2,1]),BMat))),MatA) shape=[1_pPETSCINT,dimPlex**2],order=[2,1]),BMat))),MatA)
MatB = MatB & MatB = MatB &
+ matmul(reshape(homogenization_F(1:dimPlex,1:dimPlex,m),shape=[1,dimPlex*dimPlex]),MatA) + matmul(reshape(homogenization_F(1:dimPlex,1:dimPlex,m),shape=[1_pPETSCINT,dimPlex**2]),MatA)
FAvg = FAvg + F FAvg = FAvg + F
BMatAvg = BMatAvg + BMat BMatAvg = BMatAvg + BMat
else else
@ -568,12 +571,12 @@ subroutine FEM_mechanical_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,err_P
FInv = math_inv33(FAvg) FInv = math_inv33(FAvg)
K_e = K_eA*math_det33(FAvg/real(nQuadrature))**(1.0/real(dimPlex)) + & K_e = K_eA*math_det33(FAvg/real(nQuadrature))**(1.0/real(dimPlex)) + &
(matmul(matmul(transpose(BMatAvg), & (matmul(matmul(transpose(BMatAvg), &
reshape(FInv(1:dimPlex,1:dimPlex),shape=[dimPlex*dimPlex,1],order=[2,1])),MatB) + & reshape(FInv(1:dimPlex,1:dimPlex),shape=[dimPlex**2,1_pPETSCINT],order=[2,1])),MatB) + &
K_eB)/real(dimPlex) K_eB)/real(dimPlex)
else else
K_e = K_eA K_e = K_eA
endif endif
K_e = (K_e + eps*math_eye(cellDof)) * abs(detJ) K_e = (K_e + eps*math_eye(int(cellDof))) * abs(detJ)
#ifndef __INTEL_COMPILER #ifndef __INTEL_COMPILER
pK_e(1:cellDOF**2) => K_e pK_e(1:cellDOF**2) => K_e
#else #else
@ -596,7 +599,8 @@ subroutine FEM_mechanical_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,err_P
#if (PETSC_VERSION_MINOR < 14) #if (PETSC_VERSION_MINOR < 14)
call DMPlexCreateRigidBody(dm_local,matnull,err_PETSc); CHKERRQ(err_PETSc) call DMPlexCreateRigidBody(dm_local,matnull,err_PETSc); CHKERRQ(err_PETSc)
#else #else
call DMPlexCreateRigidBody(dm_local,0,matnull,err_PETSc); CHKERRQ(err_PETSc) call DMPlexCreateRigidBody(dm_local,0_pPETSCINT,matnull,err_PETSc)
CHKERRQ(err_PETSc)
#endif #endif
call MatSetNullSpace(Jac,matnull,err_PETSc); CHKERRQ(err_PETSc) call MatSetNullSpace(Jac,matnull,err_PETSc); CHKERRQ(err_PETSc)
call MatSetNearNullSpace(Jac,matnull,err_PETSc); CHKERRQ(err_PETSc) call MatSetNearNullSpace(Jac,matnull,err_PETSc); CHKERRQ(err_PETSc)
@ -645,7 +649,7 @@ subroutine FEM_mechanical_forward(guess,timeinc,timeinc_old,fieldBC)
if (bcSize > 0) then if (bcSize > 0) then
call DMGetStratumIS(dm_local,'Face Sets',mesh_boundaries(face),bcPoints,err_PETSc) call DMGetStratumIS(dm_local,'Face Sets',mesh_boundaries(face),bcPoints,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call utilities_projectBCValues(solution_local,section,0,field-1,bcPoints, & call utilities_projectBCValues(solution_local,section,0_pPETSCINT,field-1,bcPoints, &
0.0_pReal,fieldBC%componentBC(field)%Value(face),timeinc_old) 0.0_pReal,fieldBC%componentBC(field)%Value(face),timeinc_old)
call ISDestroy(bcPoints,err_PETSc); CHKERRQ(err_PETSc) call ISDestroy(bcPoints,err_PETSc); CHKERRQ(err_PETSc)
endif endif
@ -684,7 +688,7 @@ subroutine FEM_mechanical_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reaso
if (terminallyIll) reason = SNES_DIVERGED_FUNCTION_DOMAIN if (terminallyIll) reason = SNES_DIVERGED_FUNCTION_DOMAIN
print'(/,1x,a,a,i0,a,i0,f0.3)', trim(incInfo), & print'(/,1x,a,a,i0,a,i0,f0.3)', trim(incInfo), &
' @ Iteration ',PETScIter,' mechanical residual norm = ', & ' @ Iteration ',PETScIter,' mechanical residual norm = ', &
int(fnorm/divTol),fnorm/divTol-int(fnorm/divTol) int(fnorm/divTol),fnorm/divTol-int(fnorm/divTol) ! ToDo: int casting?
print'(/,1x,a,/,2(3(2x,f12.4,1x)/),3(2x,f12.4,1x))', & print'(/,1x,a,/,2(3(2x,f12.4,1x)/),3(2x,f12.4,1x))', &
'Piola--Kirchhoff stress / MPa =',transpose(P_av)*1.e-6_pReal 'Piola--Kirchhoff stress / MPa =',transpose(P_av)*1.e-6_pReal
flush(IO_STDOUT) flush(IO_STDOUT)
@ -697,9 +701,7 @@ end subroutine FEM_mechanical_converged
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine FEM_mechanical_updateCoords() subroutine FEM_mechanical_updateCoords()
real(pReal), pointer, dimension(:) :: & PetscReal, pointer, dimension(:,:) :: &
nodeCoords_linear !< nodal coordinates (dimPlex*Nnodes)
real(pReal), pointer, dimension(:,:) :: &
nodeCoords !< nodal coordinates (3,Nnodes) nodeCoords !< nodal coordinates (3,Nnodes)
real(pReal), pointer, dimension(:,:,:) :: & real(pReal), pointer, dimension(:,:,:) :: &
ipCoords !< ip coordinates (3,nQuadrature,mesh_NcpElems) ipCoords !< ip coordinates (3,nQuadrature,mesh_NcpElems)
@ -717,8 +719,9 @@ subroutine FEM_mechanical_updateCoords()
cellStart, cellEnd, c, n cellStart, cellEnd, c, n
PetscSection :: section PetscSection :: section
PetscQuadrature :: mechQuad PetscQuadrature :: mechQuad
PetscReal, dimension(:), pointer :: basisField, basisFieldDer PetscReal, dimension(:), pointer :: basisField, basisFieldDer, &
PetscScalar, dimension(:), pointer :: x_scal nodeCoords_linear !< nodal coordinates (dimPlex*Nnodes)
PetscScalar, dimension(:), pointer :: x_scal
call SNESGetDM(mechanical_snes,dm_local,err_PETSc); CHKERRQ(err_PETSc) call SNESGetDM(mechanical_snes,dm_local,err_PETSc); CHKERRQ(err_PETSc)
call DMGetDS(dm_local,mechQuad,err_PETSc); CHKERRQ(err_PETSc) call DMGetDS(dm_local,mechQuad,err_PETSc); CHKERRQ(err_PETSc)
@ -727,7 +730,8 @@ subroutine FEM_mechanical_updateCoords()
call DMGetDimension(dm_local,dimPlex,err_PETSc); CHKERRQ(err_PETSc) call DMGetDimension(dm_local,dimPlex,err_PETSc); CHKERRQ(err_PETSc)
! write cell vertex displacements ! write cell vertex displacements
call DMPlexGetDepthStratum(dm_local,0,pStart,pEnd,err_PETSc); CHKERRQ(err_PETSc) call DMPlexGetDepthStratum(dm_local,0_pPETSCINT,pStart,pEnd,err_PETSc)
CHKERRQ(err_PETSc)
allocate(nodeCoords(3,pStart:pEnd-1),source=0.0_pReal) allocate(nodeCoords(3,pStart:pEnd-1),source=0.0_pReal)
call VecGetArrayF90(x_local,nodeCoords_linear,err_PETSc); CHKERRQ(err_PETSc) call VecGetArrayF90(x_local,nodeCoords_linear,err_PETSc); CHKERRQ(err_PETSc)
do p=pStart, pEnd-1 do p=pStart, pEnd-1