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

Merge branch '38-introduce-rudimentary-PETSc-based-FEM-solver' of magit1.mpie.de:damask/DAMASK into 38-introduce-rudimentary-PETSc-based-FEM-solver

This commit is contained in:
Martin Diehl 2018-09-22 13:41:48 +02:00
parent 42bc1b468a d9bdf53628
commit 9eb318b9bf
4 changed files with 20 additions and 12 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
examples/FEM/polyXtal/numerics.config
View File

@ -1,2 +1,3 @@
residualStiffness 0.001
charLength 0.02
petsc_options -mech_snes_type newtonls -mech_ksp_type fgmres -mech_pc_type ml -mech_ksp_monitor

10
src/DAMASK_FEM.f90
View File

@ -163,6 +163,16 @@ program DAMASK_FEM
case(FIELD_MECH_ID)
loadCases(i)%fieldBC(field)%nComponents = dimPlex !< X, Y (, Z) displacements
allocate(loadCases(i)%fieldBC(field)%componentBC(loadCases(i)%fieldBC(field)%nComponents))
do component = 1, loadCases(i)%fieldBC(field)%nComponents
select case (component)
case (1)
loadCases(i)%fieldBC(field)%componentBC(component)%ID = COMPONENT_MECH_X_ID
case (2)
loadCases(i)%fieldBC(field)%componentBC(component)%ID = COMPONENT_MECH_Y_ID
case (3)
loadCases(i)%fieldBC(field)%componentBC(component)%ID = COMPONENT_MECH_Z_ID
end select
enddo
end select
do component = 1, loadCases(i)%fieldBC(field)%nComponents
allocate(loadCases(i)%fieldBC(field)%componentBC(component)%Value(mesh_Nboundaries), source = 0.0_pReal)

17
src/FEM_mech.f90
View File

@ -7,10 +7,8 @@
module FEM_mech
#include <petsc/finclude/petscdmplex.h>
#include <petsc/finclude/petscdm.h>
#include <petsc/finclude/petscdmda.h>
#include <petsc/finclude/petsc.h>
use PETScdmda
use PETScsnes
use PETScDM
use PETScDMplex
@ -111,7 +109,6 @@ subroutine FEM_mech_init(fieldBC)
PetscInt :: cellStart, cellEnd, cell, basis
character(len=7) :: prefix = 'mechFE_'
PetscErrorCode :: ierr
PetscReal, allocatable, target, dimension(:) :: qWeights
write(6,'(/,a)') ' <<<+- FEM_mech init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
@ -131,17 +128,14 @@ subroutine FEM_mech_init(fieldBC)
qWeightsP => qWeights
call PetscQuadratureCreate(PETSC_COMM_SELF,mechQuad,ierr); CHKERRQ(ierr)
CHKERRQ(ierr)
! what is the number of components, nc?
nc = dimPlex
call PetscQuadratureSetData(mechQuad,dimPlex,nc,nQuadrature,qPointsP,qWeightsP,ierr)
CHKERRQ(ierr)
! call PetscFECreateDefault(mech_mesh,dimPlex,nc,PETSC_TRUE,prefix, &
! integrationOrder,mechFE,ierr); CHKERRQ(ierr)
call PetscFECreateDefault(PETSC_COMM_SELF,dimPlex,nc,PETSC_TRUE,prefix, &
integrationOrder,mechFE,ierr); CHKERRQ(ierr)
! Polar decomposition failed in run time
call PetscFESetQuadrature(mechFE,mechQuad,ierr); CHKERRQ(ierr)
call PetscFEGetDimension(mechFE,nBasis,ierr); CHKERRQ(ierr)
nBasis = nBasis/nc
call DMGetDS(mech_mesh,mechDS,ierr); CHKERRQ(ierr)
call PetscDSAddDiscretization(mechDS,mechFE,ierr); CHKERRQ(ierr)
call PetscDSGetTotalDimension(mechDS,cellDof,ierr); CHKERRQ(ierr)
@ -248,7 +242,7 @@ subroutine FEM_mech_init(fieldBC)
call DMPlexComputeCellGeometryAffineFEM(mech_mesh,cell,pV0,pCellJ,pInvcellJ,detJ,ierr)
CHKERRQ(ierr)
cellJMat = reshape(pCellJ,shape=[dimPlex,dimPlex])
do basis = 0, nBasis/dimPlex-1
do basis = 0, nBasis-1
call PetscDualSpaceGetFunctional(mechDualSpace,basis,functional,ierr)
CHKERRQ(ierr)
call PetscQuadratureGetData(functional,dimPlex,nc,nNodalPoints,nodalPointsP,nodalWeightsP,ierr)
@ -328,6 +322,7 @@ subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr)
materialpoint_F, &
materialpoint_P
use math, only: &
math_I3, &
math_det33, &
math_inv33
use FEsolving, only: &
@ -363,7 +358,8 @@ subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr)
CHKERRQ(ierr)
call DMPlexGetHeightStratum(dm_local,0,cellStart,cellEnd,ierr); CHKERRQ(ierr)
call DMGetLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr)
call VecWAXPY(x_local,1.0,xx_local,solution_local,ierr); CHKERRQ(ierr)
!call VecWAXPY(x_local,1.0,xx_local,solution_local,ierr); CHKERRQ(ierr)
call VecCopy(xx_local,x_local,ierr); CHKERRQ(ierr)
do field = 1, dimPlex; do face = 1, mesh_Nboundaries
if (params%fieldBC%componentBC(field)%Mask(face)) then
call DMGetStratumSize(dm_local,'Face Sets',mesh_boundaries(face),bcSize,ierr)
@ -396,6 +392,7 @@ subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr)
enddo
enddo
materialpoint_F(1:dimPlex,1:dimPlex,qPt+1,cell+1) = &
math_I3 + &
reshape(matmul(BMat,x_scal),shape=[dimPlex,dimPlex], order=[2,1])
enddo
if (BBarStabilisation) then
@ -591,7 +588,7 @@ subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr)
!--------------------------------------------------------------------------------------------------
! apply boundary conditions
call DMPlexCreateRigidBody(dm_local,matnull,ierr); CHKERRQ(ierr) !MD: linker error
call DMPlexCreateRigidBody(dm_local,matnull,ierr); CHKERRQ(ierr)
call MatSetNullSpace(Jac,matnull,ierr); CHKERRQ(ierr)
call MatSetNearNullSpace(Jac,matnull,ierr); CHKERRQ(ierr)
call MatNullSpaceDestroy(matnull,ierr); CHKERRQ(ierr)

4
src/FEM_utilities.f90
View File

@ -188,7 +188,7 @@ subroutine utilities_init()
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_defaultOptions),ierr)
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_options),ierr)
CHKERRQ(ierr)
write(petsc_optionsPhysics,'(a,i0)') '-mechFE_petscspace_order ' , structOrder
write(petsc_optionsPhysics,'(a,i0)') '-mechFE_petscspace_degree ' , structOrder
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_optionsPhysics),ierr)
CHKERRQ(ierr)
@ -411,7 +411,7 @@ subroutine utilities_projectBCValues(localVec,section,field,comp,bcPointsIS,BCVa
call PetscSectionGetFieldOffset(section,bcPoints(point),field,offset,ierr)
CHKERRQ(ierr)
do dof = offset+comp+1, offset+numDof, numComp
localArray(dof) = localArray(dof) + BCValue + BCDotValue*timeinc
localArray(dof) = BCValue + BCDotValue*timeinc
enddo
enddo
call VecRestoreArrayF90(localVec,localArray,ierr); CHKERRQ(ierr)