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

Merge branch 'MiscImprovements' of magit1.mpie.de:/damask/DAMASK into MiscImprovements

This commit is contained in:
Martin Diehl 2020-01-21 19:24:30 +01:00
parent b98ecccf73 91d25ec897
commit 0ef0db586f
10 changed files with 144 additions and 143 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
.gitignore vendored
View File

@ -1,5 +1,6 @@
*.pyc *.pyc
*.hdf5 *.hdf5
*.xdmf
*.bak *.bak
*~ *~
bin bin

2
python/damask/dadf5.py
View File

@ -906,7 +906,7 @@ class DADF5():
n_nodes = 8 n_nodes = 8
for i in f['/geometry/T_c']: for i in f['/geometry/T_c']:
vtk_geom.InsertNextCell(n_nodes,vtk_type,i-1) vtk_geom.InsertNextCell(vtk_type,n_nodes,i-1)
elif mode == 'Point': elif mode == 'Point':
Points = vtk.vtkPoints() Points = vtk.vtkPoints()

4
python/damask/orientation.py
View File

@ -1078,12 +1078,12 @@ class Orientation:
if isinstance(lattice, Lattice): if isinstance(lattice, Lattice):
self.lattice = lattice self.lattice = lattice
else: else:
self.lattice = Lattice(lattice) # assume string self.lattice = Lattice(lattice) # assume string
if isinstance(rotation, Rotation): if isinstance(rotation, Rotation):
self.rotation = rotation self.rotation = rotation
else: else:
self.rotation = Rotation(rotation) # assume quaternion self.rotation = Rotation.fromQuaternion(rotation) # assume quaternion
def disorientation(self, def disorientation(self,
other, other,

6
python/damask/table.py
View File

@ -23,11 +23,7 @@ class Table():
""" """
self.comments = [] if comments is None else [c for c in comments] self.comments = [] if comments is None else [c for c in comments]
if hasattr(data,'columns'): self.data = pd.DataFrame(data=data)
self.data = data
self.data.columns = [''] * len(self.data.columns)
else:
self.data = pd.DataFrame(data=data)
self.shapes = shapes self.shapes = shapes
self.__label_condensed() self.__label_condensed()

20
src/grid/DAMASK_grid.f90
View File

@ -65,12 +65,11 @@ program DAMASK_spectral
currentLoadcase = 0, & !< current load case currentLoadcase = 0, & !< current load case
inc, & !< current increment in current load case inc, & !< current increment in current load case
totalIncsCounter = 0, & !< total # of increments totalIncsCounter = 0, & !< total # of increments
fileUnit = 0, & !< file unit for reading load case and writing results
myStat, &
statUnit = 0, & !< file unit for statistics output statUnit = 0, & !< file unit for statistics output
stagIter, & stagIter, &
nActiveFields = 0 nActiveFields = 0
character(len=6) :: loadcase_string character(len=6) :: loadcase_string
character(len=pStringLen), dimension(:), allocatable :: fileContent
character(len=1024) :: & character(len=1024) :: &
incInfo incInfo
type(tLoadCase), allocatable, dimension(:) :: loadCases !< array of all load cases type(tLoadCase), allocatable, dimension(:) :: loadCases !< array of all load cases
@ -141,17 +140,14 @@ program DAMASK_spectral
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! reading information from load case file and to sanity checks ! reading information from load case file and to sanity checks
fileContent = IO_read_ASCII(trim(loadCaseFile))
allocate (loadCases(0)) ! array of load cases allocate (loadCases(0)) ! array of load cases
open(newunit=fileunit,iostat=myStat,file=trim(loadCaseFile),action='read') do currentLoadCase = 1, size(fileContent)
if (myStat /= 0) call IO_error(100,el=myStat,ext_msg=trim(loadCaseFile)) line = fileContent(currentLoadCase)
do if (IO_isBlank(line)) cycle
read(fileUnit, '(A)', iostat=myStat) line
if ( myStat /= 0) exit
if (IO_isBlank(line)) cycle ! skip empty lines
currentLoadCase = currentLoadCase + 1
chunkPos = IO_stringPos(line) chunkPos = IO_stringPos(line)
do i = 1, chunkPos(1) ! reading compulsory parameters for loadcase do i = 1, chunkPos(1) ! reading compulsory parameters for loadcase
select case (IO_lc(IO_stringValue(line,chunkPos,i))) select case (IO_lc(IO_stringValue(line,chunkPos,i)))
case('l','fdot','dotf','f') case('l','fdot','dotf','f')
@ -300,7 +296,7 @@ program DAMASK_spectral
endif reportAndCheck endif reportAndCheck
loadCases = [loadCases,newLoadCase] ! load case is ok, append it loadCases = [loadCases,newLoadCase] ! load case is ok, append it
enddo enddo
close(fileUnit)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! doing initialization depending on active solvers ! doing initialization depending on active solvers

63
src/mesh/DAMASK_FEM.f90
View File

@ -66,7 +66,7 @@ program DAMASK_FEM
PetscInt :: faceSet, currentFaceSet PetscInt :: faceSet, currentFaceSet
PetscInt :: field, dimPlex PetscInt :: field, dimPlex
PetscErrorCode :: ierr PetscErrorCode :: ierr
integer(kind(COMPONENT_UNDEFINED_ID)) :: ID
external :: & external :: &
quit quit
@ -166,45 +166,28 @@ program DAMASK_FEM
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! boundary condition information ! boundary condition information
case('x') ! X displacement field case('x','y','z')
do field = 1, nActiveFields select case(IO_lc(IO_stringValue(line,chunkPos,i)))
if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_MECH_ID) then case('x')
do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents ID = COMPONENT_MECH_X_ID
if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_MECH_X_ID) then case('y')
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = & ID = COMPONENT_MECH_Y_ID
.true. case('z')
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = & ID = COMPONENT_MECH_Z_ID
IO_floatValue(line,chunkPos,i+1) end select
endif
enddo do field = 1, nActiveFields
endif if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_MECH_ID) then
enddo do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents
case('y') ! Y displacement field if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == ID) then
do field = 1, nActiveFields loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = &
if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_MECH_ID) then .true.
do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_MECH_Y_ID) then IO_floatValue(line,chunkPos,i+1)
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = & endif
.true. enddo
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = & endif
IO_floatValue(line,chunkPos,i+1) enddo
endif
enddo
endif
enddo
case('z') ! Z displacement field
do field = 1, nActiveFields
if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_MECH_ID) then
do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents
if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_MECH_Z_ID) then
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = &
.true.
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
IO_floatValue(line,chunkPos,i+1)
endif
enddo
endif
enddo
end select end select
enddo; enddo enddo; enddo
close(fileUnit) close(fileUnit)

139
src/mesh/FEM_mech.f90
View File

@ -67,28 +67,33 @@ contains
subroutine FEM_mech_init(fieldBC) subroutine FEM_mech_init(fieldBC)
type(tFieldBC), intent(in) :: fieldBC type(tFieldBC), intent(in) :: fieldBC
DM :: mech_mesh DM :: mech_mesh
PetscFE :: mechFE PetscFE :: mechFE
PetscQuadrature :: mechQuad, functional PetscQuadrature :: mechQuad, functional
PetscDS :: mechDS PetscDS :: mechDS
PetscDualSpace :: mechDualSpace PetscDualSpace :: mechDualSpace
DMLabel, dimension(:),pointer :: pBCLabel
DMLabel :: BCLabel DMLabel :: BCLabel
PetscInt, dimension(:), pointer :: pNumComp, pNumDof, pBcField, pBcPoint PetscInt, dimension(:), pointer :: pNumComp, pNumDof, pBcField, pBcPoint
PetscInt :: numBC, bcSize, nc PetscInt :: numBC, bcSize, nc, &
field, faceSet, topologDim, nNodalPoints, &
cellStart, cellEnd, cell, basis
IS :: bcPoint IS :: bcPoint
IS, pointer :: pBcComps(:), pBcPoints(:) IS, dimension(:), pointer :: pBcComps, pBcPoints
PetscSection :: section PetscSection :: section
PetscInt :: field, faceSet, topologDim, nNodalPoints
PetscReal, dimension(:), pointer :: qPointsP, qWeightsP, & PetscReal, dimension(:), pointer :: qPointsP, qWeightsP, &
nodalPointsP, nodalWeightsP nodalPointsP, nodalWeightsP,pV0, pCellJ, pInvcellJ
PetscReal, allocatable, target :: nodalPoints(:), nodalWeights(:)
PetscScalar, pointer :: px_scal(:)
PetscScalar, allocatable, target :: x_scal(:)
PetscReal :: detJ PetscReal :: detJ
PetscReal, allocatable, target :: cellJMat(:,:) PetscReal, allocatable, target :: cellJMat(:,:)
PetscReal, pointer :: pV0(:), pCellJ(:), pInvcellJ(:)
PetscInt :: cellStart, cellEnd, cell, basis PetscScalar, pointer :: px_scal(:)
character(len=7), parameter :: prefix = 'mechFE_' PetscScalar, allocatable, target :: x_scal(:)
character(len=*), parameter :: prefix = 'mechFE_'
PetscErrorCode :: ierr PetscErrorCode :: ierr
write(6,'(/,a)') ' <<<+- FEM_mech init -+>>>' write(6,'(/,a)') ' <<<+- FEM_mech init -+>>>'
@ -125,13 +130,19 @@ subroutine FEM_mech_init(fieldBC)
! Setup FEM mech boundary conditions ! Setup FEM mech boundary conditions
call DMGetLabel(mech_mesh,'Face Sets',BCLabel,ierr); CHKERRQ(ierr) call DMGetLabel(mech_mesh,'Face Sets',BCLabel,ierr); CHKERRQ(ierr)
call DMPlexLabelComplete(mech_mesh,BCLabel,ierr); CHKERRQ(ierr) call DMPlexLabelComplete(mech_mesh,BCLabel,ierr); CHKERRQ(ierr)
#if (PETSC_VERSION_MINOR < 11)
call DMGetSection(mech_mesh,section,ierr); CHKERRQ(ierr) call DMGetSection(mech_mesh,section,ierr); CHKERRQ(ierr)
#else
call DMGetLocalSection(mech_mesh,section,ierr); CHKERRQ(ierr)
#endif
allocate(pnumComp(1), source=dimPlex) allocate(pnumComp(1), source=dimPlex)
allocate(pnumDof(dimPlex+1), source = 0) allocate(pnumDof(0:dimPlex), source = 0)
do topologDim = 0, dimPlex do topologDim = 0, dimPlex
call DMPlexGetDepthStratum(mech_mesh,topologDim,cellStart,cellEnd,ierr) call DMPlexGetDepthStratum(mech_mesh,topologDim,cellStart,cellEnd,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call PetscSectionGetDof(section,cellStart,pnumDof(topologDim+1),ierr) call PetscSectionGetDof(section,cellStart,pnumDof(topologDim),ierr)
write(6,*) 'start',cellStart,'end',cellEnd
write(6,*) 'topologDim',topologDim,'numDOF',pNumDOF(topologDim)
CHKERRQ(ierr) CHKERRQ(ierr)
enddo enddo
numBC = 0 numBC = 0
@ -163,9 +174,15 @@ subroutine FEM_mech_init(fieldBC)
endif endif
endif endif
enddo; enddo enddo; enddo
#if (PETSC_VERSION_MINOR < 11)
call DMPlexCreateSection(mech_mesh,dimPlex,1,pNumComp,pNumDof, & call DMPlexCreateSection(mech_mesh,dimPlex,1,pNumComp,pNumDof, &
numBC,pBcField,pBcComps,pBcPoints,PETSC_NULL_IS, & numBC,pBcField,pBcComps,pBcPoints,PETSC_NULL_IS,section,ierr)
section,ierr) #else
allocate(pBClabel(1),source=BClabel)
call DMPlexCreateSection(mech_mesh,pBClabel,pNumComp,pNumDof, &
numBC,pBcField,pBcComps,pBcPoints,PETSC_NULL_IS,section,ierr)
#endif
CHKERRQ(ierr) CHKERRQ(ierr)
call DMSetSection(mech_mesh,section,ierr); CHKERRQ(ierr) call DMSetSection(mech_mesh,section,ierr); CHKERRQ(ierr)
do faceSet = 1, numBC do faceSet = 1, numBC
@ -196,13 +213,11 @@ subroutine FEM_mech_init(fieldBC)
call VecSet(solution ,0.0,ierr); CHKERRQ(ierr) call VecSet(solution ,0.0,ierr); CHKERRQ(ierr)
call VecSet(solution_rate ,0.0,ierr); CHKERRQ(ierr) call VecSet(solution_rate ,0.0,ierr); CHKERRQ(ierr)
allocate(x_scal(cellDof)) allocate(x_scal(cellDof))
allocate(nodalPoints (dimPlex)) allocate(nodalWeightsP(1))
allocate(nodalWeights(1)) allocate(nodalPointsP(dimPlex))
nodalPointsP => nodalPoints
nodalWeightsP => nodalWeights
allocate(pv0(dimPlex)) allocate(pv0(dimPlex))
allocate(pcellJ(dimPlex*dimPlex)) allocate(pcellJ(dimPlex**2))
allocate(pinvcellJ(dimPlex*dimPlex)) allocate(pinvcellJ(dimPlex**2))
allocate(cellJMat(dimPlex,dimPlex)) allocate(cellJMat(dimPlex,dimPlex))
call DMGetSection(mech_mesh,section,ierr); CHKERRQ(ierr) call DMGetSection(mech_mesh,section,ierr); CHKERRQ(ierr)
call DMGetDS(mech_mesh,mechDS,ierr); CHKERRQ(ierr) call DMGetDS(mech_mesh,mechDS,ierr); CHKERRQ(ierr)
@ -212,7 +227,7 @@ subroutine FEM_mech_init(fieldBC)
call DMPlexGetHeightStratum(mech_mesh,0,cellStart,cellEnd,ierr) call DMPlexGetHeightStratum(mech_mesh,0,cellStart,cellEnd,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
do cell = cellStart, cellEnd-1 !< loop over all elements do cell = cellStart, cellEnd-1 !< loop over all elements
x_scal = 0.0 x_scal = 0.0_pReal
call DMPlexComputeCellGeometryAffineFEM(mech_mesh,cell,pV0,pCellJ,pInvcellJ,detJ,ierr) call DMPlexComputeCellGeometryAffineFEM(mech_mesh,cell,pV0,pCellJ,pInvcellJ,detJ,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
cellJMat = reshape(pCellJ,shape=[dimPlex,dimPlex]) cellJMat = reshape(pCellJ,shape=[dimPlex,dimPlex])
@ -221,7 +236,7 @@ subroutine FEM_mech_init(fieldBC)
CHKERRQ(ierr) CHKERRQ(ierr)
call PetscQuadratureGetData(functional,dimPlex,nc,nNodalPoints,nodalPointsP,nodalWeightsP,ierr) call PetscQuadratureGetData(functional,dimPlex,nc,nNodalPoints,nodalPointsP,nodalWeightsP,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
x_scal(basis+1:basis+dimPlex) = pV0 + matmul(transpose(cellJMat),nodalPointsP + 1.0) x_scal(basis+1:basis+dimPlex) = pV0 + matmul(transpose(cellJMat),nodalPointsP + 1.0_pReal)
enddo enddo
px_scal => x_scal px_scal => x_scal
call DMPlexVecSetClosure(mech_mesh,section,solution_local,cell,px_scal,INSERT_ALL_VALUES,ierr) call DMPlexVecSetClosure(mech_mesh,section,solution_local,cell,px_scal,INSERT_ALL_VALUES,ierr)
@ -283,6 +298,9 @@ end function FEM_mech_solution
subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr) subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr)
DM :: dm_local DM :: dm_local
PetscObject,intent(in) :: dummy
PetscErrorCode :: ierr
PetscDS :: prob PetscDS :: prob
Vec :: x_local, f_local, xx_local Vec :: x_local, f_local, xx_local
PetscSection :: section PetscSection :: section
@ -294,10 +312,10 @@ subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr)
qPt, basis, comp, cidx qPt, basis, comp, cidx
PetscReal :: detFAvg PetscReal :: detFAvg
PetscReal :: BMat(dimPlex*dimPlex,cellDof) PetscReal :: BMat(dimPlex*dimPlex,cellDof)
PetscObject,intent(in) :: dummy
PetscInt :: bcSize PetscInt :: bcSize
IS :: bcPoints IS :: bcPoints
PetscErrorCode :: ierr
allocate(pV0(dimPlex)) allocate(pV0(dimPlex))
allocate(pcellJ(dimPlex**2)) allocate(pcellJ(dimPlex**2))
@ -316,7 +334,7 @@ subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr)
call DMGetStratumIS(dm_local,'Face Sets',mesh_boundaries(face),bcPoints,ierr) call DMGetStratumIS(dm_local,'Face Sets',mesh_boundaries(face),bcPoints,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call utilities_projectBCValues(x_local,section,0,field-1,bcPoints, & call utilities_projectBCValues(x_local,section,0,field-1,bcPoints, &
0.0,params%fieldBC%componentBC(field)%Value(face),params%timeinc) 0.0_pReal,params%fieldBC%componentBC(field)%Value(face),params%timeinc)
call ISDestroy(bcPoints,ierr); CHKERRQ(ierr) call ISDestroy(bcPoints,ierr); CHKERRQ(ierr)
endif endif
endif endif
@ -403,29 +421,35 @@ end subroutine FEM_mech_formResidual
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr) subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr)
DM :: dm_local
DM :: dm_local
Mat :: Jac_pre, Jac
PetscObject, intent(in) :: dummy
PetscErrorCode :: ierr
PetscDS :: prob PetscDS :: prob
Vec :: x_local, xx_local Vec :: x_local, xx_local
Mat :: Jac_pre, Jac
PetscSection :: section, gSection PetscSection :: section, gSection
PetscReal, dimension(1, cellDof) :: MatB
PetscReal, dimension(dimPlex**2,cellDof) :: BMat, BMatAvg, MatA
PetscReal, dimension(3,3) :: F, FAvg, FInv
PetscReal :: detJ PetscReal :: detJ
PetscReal, dimension(:), pointer :: basisField, basisFieldDer, & PetscReal, dimension(:), pointer :: basisField, basisFieldDer, &
pV0, pCellJ, pInvcellJ pV0, pCellJ, pInvcellJ
PetscScalar, dimension(:), pointer :: pK_e, x_scal
PetscScalar,dimension(cellDOF,cellDOF), target :: K_e
PetscScalar,dimension(cellDOF,cellDOF) :: K_eA , &
K_eB
PetscInt :: cellStart, cellEnd, cell, field, face, & PetscInt :: cellStart, cellEnd, cell, field, face, &
qPt, basis, comp, cidx,bcSize qPt, basis, comp, cidx,bcSize
PetscScalar,dimension(cellDOF,cellDOF), target :: K_e, &
K_eA , &
K_eB
PetscScalar, target :: K_eVec(cellDof*cellDof)
PetscReal :: BMat (dimPlex*dimPlex,cellDof), &
BMatAvg(dimPlex*dimPlex,cellDof), &
MatA (dimPlex*dimPlex,cellDof), &
MatB (1 ,cellDof)
PetscScalar, dimension(:), pointer :: pK_e, x_scal
PetscReal, dimension(3,3) :: F, FAvg, FInv
PetscObject, intent(in) :: dummy
IS :: bcPoints IS :: bcPoints
PetscErrorCode :: ierr
allocate(pV0(dimPlex)) allocate(pV0(dimPlex))
allocate(pcellJ(dimPlex**2)) allocate(pcellJ(dimPlex**2))
@ -440,7 +464,7 @@ subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr)
call DMGetGlobalSection(dm_local,gSection,ierr); CHKERRQ(ierr) call DMGetGlobalSection(dm_local,gSection,ierr); CHKERRQ(ierr)
call DMGetLocalVector(dm_local,x_local,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_pReal,xx_local,solution_local,ierr); CHKERRQ(ierr)
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,ierr) call DMGetStratumSize(dm_local,'Face Sets',mesh_boundaries(face),bcSize,ierr)
@ -448,7 +472,7 @@ subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr)
call DMGetStratumIS(dm_local,'Face Sets',mesh_boundaries(face),bcPoints,ierr) call DMGetStratumIS(dm_local,'Face Sets',mesh_boundaries(face),bcPoints,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call utilities_projectBCValues(x_local,section,0,field-1,bcPoints, & call utilities_projectBCValues(x_local,section,0,field-1,bcPoints, &
0.0,params%fieldBC%componentBC(field)%Value(face),params%timeinc) 0.0_pReal,params%fieldBC%componentBC(field)%Value(face),params%timeinc)
call ISDestroy(bcPoints,ierr); CHKERRQ(ierr) call ISDestroy(bcPoints,ierr); CHKERRQ(ierr)
endif endif
endif endif
@ -501,13 +525,16 @@ subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr)
(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*dimPlex,1],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_identity2nd(cellDof) K_e = (K_e + eps*math_identity2nd(cellDof)) * abs(detJ)
K_eVec = reshape(K_e, [cellDof*cellDof])*abs(detJ) #ifndef __INTEL_COMPILER
pK_e => K_eVec pK_e(1:cellDOF**2) => K_e
#else
! https://software.intel.com/en-us/forums/intel-fortran-compiler/topic/782230 (bug)
allocate(pK_e(cellDOF**2),source = reshape(K_e,[cellDOF**2]))
#endif
call DMPlexMatSetClosure(dm_local,section,gSection,Jac,cell,pK_e,ADD_VALUES,ierr) call DMPlexMatSetClosure(dm_local,section,gSection,Jac,cell,pK_e,ADD_VALUES,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call DMPlexVecRestoreClosure(dm_local,section,x_local,cell,x_scal,ierr) call DMPlexVecRestoreClosure(dm_local,section,x_local,cell,x_scal,ierr)
@ -536,18 +563,18 @@ subroutine FEM_mech_forward(guess,timeinc,timeinc_old,fieldBC)
type(tFieldBC), intent(in) :: & type(tFieldBC), intent(in) :: &
fieldBC fieldBC
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
timeinc_old, & timeinc_old, &
timeinc timeinc
logical, intent(in) :: & logical, intent(in) :: &
guess guess
PetscInt :: field, face
DM :: dm_local PetscInt :: field, face, bcSize
Vec :: x_local DM :: dm_local
PetscSection :: section Vec :: x_local
PetscInt :: bcSize PetscSection :: section
IS :: bcPoints IS :: bcPoints
PetscErrorCode :: ierr PetscErrorCode :: ierr
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! forward last inc ! forward last inc
@ -557,7 +584,7 @@ subroutine FEM_mech_forward(guess,timeinc,timeinc_old,fieldBC)
call SNESGetDM(mech_snes,dm_local,ierr); CHKERRQ(ierr) !< retrieve mesh info from mech_snes into dm_local call SNESGetDM(mech_snes,dm_local,ierr); CHKERRQ(ierr) !< retrieve mesh info from mech_snes into dm_local
call DMGetSection(dm_local,section,ierr); CHKERRQ(ierr) call DMGetSection(dm_local,section,ierr); CHKERRQ(ierr)
call DMGetLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr) call DMGetLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr)
call VecSet(x_local,0.0,ierr); CHKERRQ(ierr) call VecSet(x_local,0.0_pReal,ierr); CHKERRQ(ierr)
call DMGlobalToLocalBegin(dm_local,solution,INSERT_VALUES,x_local,ierr) !< retrieve my partition of global solution vector call DMGlobalToLocalBegin(dm_local,solution,INSERT_VALUES,x_local,ierr) !< retrieve my partition of global solution vector
CHKERRQ(ierr) CHKERRQ(ierr)
call DMGlobalToLocalEnd(dm_local,solution,INSERT_VALUES,x_local,ierr) call DMGlobalToLocalEnd(dm_local,solution,INSERT_VALUES,x_local,ierr)
@ -570,7 +597,7 @@ subroutine FEM_mech_forward(guess,timeinc,timeinc_old,fieldBC)
call DMGetStratumIS(dm_local,'Face Sets',mesh_boundaries(face),bcPoints,ierr) call DMGetStratumIS(dm_local,'Face Sets',mesh_boundaries(face),bcPoints,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call utilities_projectBCValues(solution_local,section,0,field-1,bcPoints, & call utilities_projectBCValues(solution_local,section,0,field-1,bcPoints, &
0.0,fieldBC%componentBC(field)%Value(face),timeinc_old) 0.0_pReal,fieldBC%componentBC(field)%Value(face),timeinc_old)
call ISDestroy(bcPoints,ierr); CHKERRQ(ierr) call ISDestroy(bcPoints,ierr); CHKERRQ(ierr)
endif endif
endif endif

17
src/mesh/FEM_utilities.f90
View File

@ -23,7 +23,7 @@ module FEM_utilities
private private
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
logical, public :: cutBack = .false. !< cut back of BVP solver in case convergence is not achieved or a material point is terminally ill logical, public :: cutBack = .false. !< cut back of BVP solver in case convergence is not achieved or a material point is terminally ill
integer, public, parameter :: maxFields = 6 integer, public, parameter :: maxFields = 6
integer, public :: nActiveFields = 0 integer, public :: nActiveFields = 0
@ -56,15 +56,15 @@ module FEM_utilities
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! derived types ! derived types
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 integer :: iterationsNeeded = 0
end type tSolutionState end type tSolutionState
type, public :: tComponentBC type, public :: tComponentBC
integer(kind(COMPONENT_UNDEFINED_ID)) :: ID integer(kind(COMPONENT_UNDEFINED_ID)) :: ID
real(pReal), allocatable :: Value(:) real(pReal), allocatable, dimension(:) :: Value
logical, allocatable :: Mask(:) logical, allocatable, dimension(:) :: Mask
end type tComponentBC end type tComponentBC
type, public :: tFieldBC type, public :: tFieldBC
@ -79,8 +79,8 @@ module FEM_utilities
outputfrequency = 1, & !< frequency of result writes outputfrequency = 1, & !< frequency of result writes
logscale = 0 !< linear/logarithmic time inc flag logscale = 0 !< linear/logarithmic time inc flag
logical :: followFormerTrajectory = .true. !< follow trajectory of former loadcase logical :: followFormerTrajectory = .true. !< follow trajectory of former loadcase
integer, allocatable :: faceID(:) integer, allocatable, dimension(:) :: faceID
type(tFieldBC), allocatable :: fieldBC(:) type(tFieldBC), allocatable, dimension(:) :: fieldBC
end type tLoadCase end type tLoadCase
public :: & public :: &
@ -88,6 +88,7 @@ module FEM_utilities
utilities_constitutiveResponse, & utilities_constitutiveResponse, &
utilities_projectBCValues, & utilities_projectBCValues, &
FIELD_MECH_ID, & FIELD_MECH_ID, &
COMPONENT_UNDEFINED_ID, &
COMPONENT_MECH_X_ID, & COMPONENT_MECH_X_ID, &
COMPONENT_MECH_Y_ID, & COMPONENT_MECH_Y_ID, &
COMPONENT_MECH_Z_ID COMPONENT_MECH_Z_ID

4
src/mesh/FEM_zoo.f90
View File

@ -37,7 +37,7 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine FEM_Zoo_init subroutine FEM_Zoo_init
write(6,'(/,a)') ' <<<+- FEM_Zoo init -+>>>' write(6,'(/,a)') ' <<<+- FEM_Zoo init -+>>>'; flush(6)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! 2D linear ! 2D linear
@ -53,7 +53,7 @@ subroutine FEM_Zoo_init
! 2D quadratic ! 2D quadratic
FEM_Zoo_nQuadrature(2,2) = 3 FEM_Zoo_nQuadrature(2,2) = 3
allocate(FEM_Zoo_QuadratureWeights(2,2)%p(3)) allocate(FEM_Zoo_QuadratureWeights(2,2)%p(3))
FEM_Zoo_QuadratureWeights(2,2)%p(1:3) = 1.0_pReal/3.0_pReal FEM_Zoo_QuadratureWeights(2,2)%p(1:3) = 1.0_pReal/3.0_pReal
allocate(FEM_Zoo_QuadraturePoints (2,2)%p(6)) allocate(FEM_Zoo_QuadraturePoints (2,2)%p(6))

31
src/mesh_FEM.f90
View File

@ -26,9 +26,11 @@ module mesh
integer, public, protected :: & integer, public, protected :: &
mesh_Nboundaries, & mesh_Nboundaries, &
mesh_NcpElems, & !< total number of CP elements in mesh
mesh_NcpElemsGlobal mesh_NcpElemsGlobal
integer :: &
mesh_NcpElems !< total number of CP elements in mesh
!!!! BEGIN DEPRECATED !!!!! !!!! BEGIN DEPRECATED !!!!!
integer, public, protected :: & integer, public, protected :: &
mesh_maxNips !< max number of IPs in any CP element mesh_maxNips !< max number of IPs in any CP element
@ -44,7 +46,7 @@ module mesh
mesh_ipVolume, & !< volume associated with IP (initially!) mesh_ipVolume, & !< volume associated with IP (initially!)
mesh_node0 !< node x,y,z coordinates (initially!) mesh_node0 !< node x,y,z coordinates (initially!)
real(pReal), dimension(:,:,:), allocatable, public :: & real(pReal), dimension(:,:,:), allocatable :: &
mesh_ipCoordinates !< IP x,y,z coordinates (after deformation!) mesh_ipCoordinates !< IP x,y,z coordinates (after deformation!)
DM, public :: geomMesh DM, public :: geomMesh
@ -177,7 +179,7 @@ subroutine mesh_init
call IO_error(602,ext_msg='IP') ! selected element does not have requested IP call IO_error(602,ext_msg='IP') ! selected element does not have requested IP
FEsolving_execElem = [ 1,mesh_NcpElems ] ! parallel loop bounds set to comprise all DAMASK elements FEsolving_execElem = [ 1,mesh_NcpElems ] ! parallel loop bounds set to comprise all DAMASK elements
FEsolving_execIP = spread([1,FE_Nips(FE_geomtype(mesh_element(2,1))),2,nElems) FEsolving_execIP = spread([1,FE_Nips(FE_geomtype(mesh_element(2,1)))],2,mesh_NcpElems)
allocate(mesh_node0(3,mesh_Nnodes),source=0.0_pReal) allocate(mesh_node0(3,mesh_Nnodes),source=0.0_pReal)
@ -193,21 +195,17 @@ end subroutine mesh_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine mesh_FEM_build_ipVolumes(dimPlex) subroutine mesh_FEM_build_ipVolumes(dimPlex)
PetscInt :: dimPlex PetscInt,intent(in):: dimPlex
PetscReal :: vol PetscReal :: vol
PetscReal, target :: cent(dimPlex), norm(dimPlex) PetscReal, pointer,dimension(:) :: pCent, pNorm
PetscReal, pointer :: pCent(:), pNorm(:)
PetscInt :: cellStart, cellEnd, cell PetscInt :: cellStart, cellEnd, cell
PetscErrorCode :: ierr PetscErrorCode :: ierr
if (.not. allocated(mesh_ipVolume)) then allocate(mesh_ipVolume(mesh_maxNips,mesh_NcpElems),source=0.0_pReal)
allocate(mesh_ipVolume(mesh_maxNips,mesh_NcpElems))
mesh_ipVolume = 0.0_pReal
endif
call DMPlexGetHeightStratum(geomMesh,0,cellStart,cellEnd,ierr); CHKERRQ(ierr) call DMPlexGetHeightStratum(geomMesh,0,cellStart,cellEnd,ierr); CHKERRQ(ierr)
pCent => cent allocate(pCent(dimPlex))
pNorm => norm allocate(pNorm(dimPlex))
do cell = cellStart, cellEnd-1 do cell = cellStart, cellEnd-1
call DMPlexComputeCellGeometryFVM(geomMesh,cell,vol,pCent,pNorm,ierr) call DMPlexComputeCellGeometryFVM(geomMesh,cell,vol,pCent,pNorm,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
@ -225,8 +223,7 @@ subroutine mesh_FEM_build_ipCoordinates(dimPlex,qPoints)
PetscInt, intent(in) :: dimPlex PetscInt, intent(in) :: dimPlex
PetscReal, intent(in) :: qPoints(mesh_maxNips*dimPlex) PetscReal, intent(in) :: qPoints(mesh_maxNips*dimPlex)
PetscReal, target :: v0(dimPlex), cellJ(dimPlex*dimPlex), invcellJ(dimPlex*dimPlex) PetscReal, pointer,dimension(:) :: pV0, pCellJ, pInvcellJ
PetscReal, pointer :: pV0(:), pCellJ(:), pInvcellJ(:)
PetscReal :: detJ PetscReal :: detJ
PetscInt :: cellStart, cellEnd, cell, qPt, dirI, dirJ, qOffset PetscInt :: cellStart, cellEnd, cell, qPt, dirI, dirJ, qOffset
PetscErrorCode :: ierr PetscErrorCode :: ierr
@ -234,9 +231,9 @@ subroutine mesh_FEM_build_ipCoordinates(dimPlex,qPoints)
allocate(mesh_ipCoordinates(3,mesh_maxNips,mesh_NcpElems),source=0.0_pReal) allocate(mesh_ipCoordinates(3,mesh_maxNips,mesh_NcpElems),source=0.0_pReal)
pV0 => v0 allocate(pV0(dimPlex))
pCellJ => cellJ allocatE(pCellJ(dimPlex**2))
pInvcellJ => invcellJ allocatE(pinvCellJ(dimPlex**2))
call DMPlexGetHeightStratum(geomMesh,0,cellStart,cellEnd,ierr); CHKERRQ(ierr) call DMPlexGetHeightStratum(geomMesh,0,cellStart,cellEnd,ierr); CHKERRQ(ierr)
do cell = cellStart, cellEnd-1 !< loop over all elements do cell = cellStart, cellEnd-1 !< loop over all elements
call DMPlexComputeCellGeometryAffineFEM(geomMesh,cell,pV0,pCellJ,pInvcellJ,detJ,ierr) call DMPlexComputeCellGeometryAffineFEM(geomMesh,cell,pV0,pCellJ,pInvcellJ,detJ,ierr)