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

Merge branch 'write_ipdisplacements' into 'development' 

Write ip displacements

See merge request damask/DAMASK!414
This commit is contained in:
Abisheik Panneerselvam 2021-07-26 09:01:04 +00:00
parent 18dbda7843 11cdc63a27
commit 6a9c892d7b
10 changed files with 326 additions and 343 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
PRIVATE

@ -1 +1 @@
Subproject commit 174ecac2d3ab7596bdb60184d6bb9e1a52cb7378
Subproject commit f2e2d6c71ea798bfc63230a756b7cf9748599bec

4
examples/mesh/tensionZ.load
View File

@ -1,11 +1,11 @@
#initial elastic step
$Loadcase 1 time 0.0005 incs 1 frequency 5
$Loadcase 1 t 0.0005 N 1 f_out 5
Face 1 X 0.01
Face 2 X 0.0
Face 2 Y 0.0
Face 2 Z 0.0
$EndLoadcase
$Loadcase 2 time 10.0 incs 200 frequency 5
$Loadcase 2 t 10.0 N 200 f_out 5
Face 1 X 0.01
Face 2 X 0.0
Face 2 Y 0.0

178
src/mesh/DAMASK_mesh.f90
View File

@ -22,6 +22,15 @@ program DAMASK_mesh
implicit none
type :: tLoadCase
real(pReal) :: time = 0.0_pReal !< length of increment
integer :: incs = 0, & !< number of increments
outputfrequency = 1 !< frequency of result writes
logical :: followFormerTrajectory = .true. !< follow trajectory of former loadcase
integer, allocatable, dimension(:) :: faceID
type(tFieldBC), allocatable, dimension(:) :: fieldBC
end type tLoadCase
!--------------------------------------------------------------------------------------------------
! variables related to information from load case and geom file
integer, allocatable, dimension(:) :: chunkPos ! this is longer than needed for geometry parsing
@ -68,7 +77,7 @@ program DAMASK_mesh
type(tLoadCase), allocatable, dimension(:) :: loadCases !< array of all load cases
type(tSolutionState), allocatable, dimension(:) :: solres
PetscInt :: faceSet, currentFaceSet, field, dimPlex
PetscInt :: faceSet, currentFaceSet, dimPlex
PetscErrorCode :: ierr
integer(kind(COMPONENT_UNDEFINED_ID)) :: ID
external :: &
@ -91,8 +100,7 @@ program DAMASK_mesh
! reading basic information from load case file and allocate data structure containing load cases
call DMGetDimension(geomMesh,dimPlex,ierr) !< dimension of mesh (2D or 3D)
CHKERRA(ierr)
nActiveFields = 1
allocate(solres(nActiveFields))
allocate(solres(1))
!--------------------------------------------------------------------------------------------------
! reading basic information from load case file and allocate data structure containing load cases
@ -103,8 +111,8 @@ program DAMASK_mesh
chunkPos = IO_stringPos(line)
do i = 1, chunkPos(1) ! reading compulsory parameters for loadcase
select case (IO_lc(IO_stringValue(line,chunkPos,i)))
case('$loadcase')
select case (IO_stringValue(line,chunkPos,i))
case('$Loadcase')
N_def = N_def + 1
end select
enddo ! count all identifiers to allocate memory and do sanity check
@ -114,32 +122,26 @@ program DAMASK_mesh
allocate(loadCases(N_def))
do i = 1, size(loadCases)
allocate(loadCases(i)%fieldBC(nActiveFields))
field = 1
loadCases(i)%fieldBC(field)%ID = FIELD_MECH_ID
allocate(loadCases(i)%fieldBC(1))
loadCases(i)%fieldBC(1)%ID = FIELD_MECH_ID
enddo
do i = 1, size(loadCases)
do field = 1, nActiveFields
select case (loadCases(i)%fieldBC(field)%ID)
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
loadCases(i)%fieldBC(1)%nComponents = dimPlex !< X, Y (, Z) displacements
allocate(loadCases(i)%fieldBC(1)%componentBC(loadCases(i)%fieldBC(1)%nComponents))
do component = 1, loadCases(i)%fieldBC(1)%nComponents
select case (component)
case (1)
loadCases(i)%fieldBC(1)%componentBC(component)%ID = COMPONENT_MECH_X_ID
case (2)
loadCases(i)%fieldBC(1)%componentBC(component)%ID = COMPONENT_MECH_Y_ID
case (3)
loadCases(i)%fieldBC(1)%componentBC(component)%ID = COMPONENT_MECH_Z_ID
end select
do component = 1, loadCases(i)%fieldBC(field)%nComponents
allocate(loadCases(i)%fieldBC(field)%componentBC(component)%Value(mesh_Nboundaries), source = 0.0_pReal)
allocate(loadCases(i)%fieldBC(field)%componentBC(component)%Mask (mesh_Nboundaries), source = .false.)
enddo
enddo
do component = 1, loadCases(i)%fieldBC(1)%nComponents
allocate(loadCases(i)%fieldBC(1)%componentBC(component)%Value(mesh_Nboundaries), source = 0.0_pReal)
allocate(loadCases(i)%fieldBC(1)%componentBC(component)%Mask (mesh_Nboundaries), source = .false.)
enddo
enddo
@ -151,52 +153,45 @@ program DAMASK_mesh
chunkPos = IO_stringPos(line)
do i = 1, chunkPos(1)
select case (IO_lc(IO_stringValue(line,chunkPos,i)))
select case (IO_stringValue(line,chunkPos,i))
!--------------------------------------------------------------------------------------------------
! loadcase information
case('$loadcase')
case('$Loadcase')
currentLoadCase = IO_intValue(line,chunkPos,i+1)
case('face')
case('Face')
currentFace = IO_intValue(line,chunkPos,i+1)
currentFaceSet = -1
do faceSet = 1, mesh_Nboundaries
if (mesh_boundaries(faceSet) == currentFace) currentFaceSet = faceSet
enddo
if (currentFaceSet < 0) call IO_error(error_ID = 837, ext_msg = 'invalid BC')
case('t','time','delta') ! increment time
case('t')
loadCases(currentLoadCase)%time = IO_floatValue(line,chunkPos,i+1)
case('n','incs','increments','steps') ! number of increments
case('N')
loadCases(currentLoadCase)%incs = IO_intValue(line,chunkPos,i+1)
case('logincs','logincrements','logsteps') ! number of increments (switch to log time scaling)
loadCases(currentLoadCase)%incs = IO_intValue(line,chunkPos,i+1)
loadCases(currentLoadCase)%logscale = 1
case('freq','frequency','outputfreq') ! frequency of result writings
case('f_out')
loadCases(currentLoadCase)%outputfrequency = IO_intValue(line,chunkPos,i+1)
case('guessreset','dropguessing')
case('estimate_rate')
loadCases(currentLoadCase)%followFormerTrajectory = .false. ! do not continue to predict deformation along former trajectory
!--------------------------------------------------------------------------------------------------
! boundary condition information
case('x','y','z')
select case(IO_lc(IO_stringValue(line,chunkPos,i)))
case('x')
case('X','Y','Z')
select case(IO_stringValue(line,chunkPos,i))
case('X')
ID = COMPONENT_MECH_X_ID
case('y')
case('Y')
ID = COMPONENT_MECH_Y_ID
case('z')
case('Z')
ID = COMPONENT_MECH_Z_ID
end select
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 == 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
do component = 1, loadcases(currentLoadCase)%fieldBC(1)%nComponents
if (loadCases(currentLoadCase)%fieldBC(1)%componentBC(component)%ID == ID) then
loadCases(currentLoadCase)%fieldBC(1)%componentBC(component)%Mask (currentFaceSet) = &
.true.
loadCases(currentLoadCase)%fieldBC(1)%componentBC(component)%Value(currentFaceSet) = &
IO_floatValue(line,chunkPos,i+1)
endif
enddo
end select
@ -212,21 +207,16 @@ program DAMASK_mesh
print'(a,i0)', ' load case: ', currentLoadCase
if (.not. loadCases(currentLoadCase)%followFormerTrajectory) &
print'(a)', ' drop guessing along trajectory'
do field = 1, nActiveFields
select case (loadCases(currentLoadCase)%fieldBC(field)%ID)
case(FIELD_MECH_ID)
print'(a)', ' Field '//trim(FIELD_MECH_label)
print'(a)', ' Field '//trim(FIELD_MECH_label)
end select
do faceSet = 1, mesh_Nboundaries
do component = 1, loadCases(currentLoadCase)%fieldBC(field)%nComponents
if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask(faceSet)) &
print'(a,i2,a,i2,a,f12.7)', ' Face ', mesh_boundaries(faceSet), &
' Component ', component, &
' Value ', loadCases(currentLoadCase)%fieldBC(field)% &
componentBC(component)%Value(faceSet)
enddo
enddo
do faceSet = 1, mesh_Nboundaries
do component = 1, loadCases(currentLoadCase)%fieldBC(1)%nComponents
if (loadCases(currentLoadCase)%fieldBC(1)%componentBC(component)%Mask(faceSet)) &
print'(a,i2,a,i2,a,f12.7)', ' Face ', mesh_boundaries(faceSet), &
' Component ', component, &
' Value ', loadCases(currentLoadCase)%fieldBC(1)% &
componentBC(component)%Value(faceSet)
enddo
enddo
print'(a,f12.6)', ' time: ', loadCases(currentLoadCase)%time
if (loadCases(currentLoadCase)%incs < 1) errorID = 835 ! non-positive incs count
@ -240,12 +230,7 @@ program DAMASK_mesh
!--------------------------------------------------------------------------------------------------
! doing initialization depending on active solvers
call FEM_Utilities_init
do field = 1, nActiveFields
select case (loadCases(1)%fieldBC(field)%ID)
case(FIELD_MECH_ID)
call FEM_mechanical_init(loadCases(1)%fieldBC(field))
end select
enddo
call FEM_mechanical_init(loadCases(1)%fieldBC(1))
if (worldrank == 0) then
open(newunit=statUnit,file=trim(getSolverJobName())//'.sta',form='FORMATTED',status='REPLACE')
@ -266,32 +251,14 @@ program DAMASK_mesh
!--------------------------------------------------------------------------------------------------
! forwarding time
timeIncOld = timeinc ! last timeinc that brought former inc to an end
if (loadCases(currentLoadCase)%logscale == 0) then ! linear scale
timeinc = loadCases(currentLoadCase)%time/real(loadCases(currentLoadCase)%incs,pReal)
else
if (currentLoadCase == 1) then ! 1st load case of logarithmic scale
if (inc == 1) then ! 1st inc of 1st load case of logarithmic scale
timeinc = loadCases(1)%time*(2.0_pReal**real( 1-loadCases(1)%incs ,pReal)) ! assume 1st inc is equal to 2nd
else ! not-1st inc of 1st load case of logarithmic scale
timeinc = loadCases(1)%time*(2.0_pReal**real(inc-1-loadCases(1)%incs ,pReal))
endif
else ! not-1st load case of logarithmic scale
timeinc = time0 * &
( (1.0_pReal + loadCases(currentLoadCase)%time/time0 )**(real( inc,pReal)/&
real(loadCases(currentLoadCase)%incs ,pReal))&
-(1.0_pReal + loadCases(currentLoadCase)%time/time0 )**(real( inc-1 ,pReal)/&
real(loadCases(currentLoadCase)%incs ,pReal)))
endif
endif
timeinc = loadCases(currentLoadCase)%time/real(loadCases(currentLoadCase)%incs,pReal)
timeinc = timeinc * real(subStepFactor,pReal)**real(-cutBackLevel,pReal) ! depending on cut back level, decrease time step
stepFraction = 0 ! fraction scaled by stepFactor**cutLevel
stepFraction = 0 ! fraction scaled by stepFactor**cutLevel
subStepLooping: do while (stepFraction < subStepFactor**cutBackLevel)
remainingLoadCaseTime = loadCases(currentLoadCase)%time+time0 - time
time = time + timeinc ! forward target time
stepFraction = stepFraction + 1 ! count step
time = time + timeinc ! forward target time
stepFraction = stepFraction + 1 ! count step
!--------------------------------------------------------------------------------------------------
! report begin of new step
@ -306,33 +273,16 @@ program DAMASK_mesh
'-',stepFraction, '/', subStepFactor**cutBackLevel
flush(IO_STDOUT)
!--------------------------------------------------------------------------------------------------
! forward fields
do field = 1, nActiveFields
select case (loadCases(currentLoadCase)%fieldBC(field)%ID)
case(FIELD_MECH_ID)
call FEM_mechanical_forward (&
guess,timeinc,timeIncOld,loadCases(currentLoadCase)%fieldBC(field))
end select
enddo
call FEM_mechanical_forward(guess,timeinc,timeIncOld,loadCases(currentLoadCase)%fieldBC(1))
!--------------------------------------------------------------------------------------------------
! solve fields
stagIter = 0
stagIterate = .true.
do while (stagIterate)
do field = 1, nActiveFields
select case (loadCases(currentLoadCase)%fieldBC(field)%ID)
case(FIELD_MECH_ID)
solres(field) = FEM_mechanical_solution (&
incInfo,timeinc,timeIncOld,loadCases(currentLoadCase)%fieldBC(field))
solres(1) = FEM_mechanical_solution(incInfo,timeinc,timeIncOld,loadCases(currentLoadCase)%fieldBC(1))
if(.not. solres(1)%converged) exit
end select
if(.not. solres(field)%converged) exit ! no solution found
enddo
stagIter = stagIter + 1
stagIterate = stagIter < stagItMax &
.and. all(solres(:)%converged) &

340
src/mesh/FEM_quadrature.f90
View File

@ -9,7 +9,7 @@ module FEM_quadrature
private
integer, parameter :: &
maxOrder = 5 !< current max interpolation set at cubic (intended to be arbitrary)
maxOrder = 5 !< maximum integration order
real(pReal), dimension(2,3), parameter :: &
triangle = reshape([-1.0_pReal, -1.0_pReal, &
1.0_pReal, -1.0_pReal, &
@ -20,8 +20,12 @@ module FEM_quadrature
-1.0_pReal, 1.0_pReal, -1.0_pReal, &
-1.0_pReal, -1.0_pReal, 1.0_pReal], shape=[3,4])
type :: group_float !< variable length datatype
real(pReal), dimension(:), allocatable :: p
end type group_float
integer, dimension(2:3,maxOrder), public, protected :: &
FEM_nQuadrature !< number of quadrature points for a given spatial dimension(2-3) and interpolation order(1-maxOrder)
FEM_nQuadrature !< number of quadrature points for spatial dimension(2-3) and interpolation order (1-maxOrder)
type(group_float), dimension(2:3,maxOrder), public, protected :: &
FEM_quadrature_weights, & !< quadrature weights for each quadrature rule
FEM_quadrature_points !< quadrature point coordinates (in simplical system) for each quadrature rule
@ -35,145 +39,146 @@ contains
!--------------------------------------------------------------------------------------------------
!> @brief initializes FEM interpolation data
!--------------------------------------------------------------------------------------------------
subroutine FEM_quadrature_init
subroutine FEM_quadrature_init()
print'(/,a)', ' <<<+- FEM_quadrature init -+>>>'; flush(6)
print*, 'L. Zhang et al., Journal of Computational Mathematics 27(1):89-96, 2009'
print*, 'https://www.jstor.org/stable/43693493'
!--------------------------------------------------------------------------------------------------
! 2D linear
FEM_nQuadrature(2,1) = 1
allocate(FEM_quadrature_weights(2,1)%p(1))
FEM_quadrature_weights(2,1)%p(1) = 1.0_pReal
allocate(FEM_quadrature_weights(2,1)%p(FEM_nQuadrature(2,1)))
FEM_quadrature_weights(2,1)%p(1) = 1._pReal
allocate(FEM_quadrature_points (2,1)%p(2))
FEM_quadrature_points (2,1)%p(1:2) = permutationStar3([1.0_pReal/3.0_pReal])
FEM_quadrature_points (2,1)%p = permutationStar3([1._pReal/3._pReal])
!--------------------------------------------------------------------------------------------------
! 2D quadratic
FEM_nQuadrature(2,2) = 3
allocate(FEM_quadrature_weights(2,2)%p(3))
FEM_quadrature_weights(2,2)%p(1:3) = 1.0_pReal/3.0_pReal
allocate(FEM_quadrature_weights(2,2)%p(FEM_nQuadrature(2,2)))
FEM_quadrature_weights(2,2)%p(1:3) = 1._pReal/3._pReal
allocate(FEM_quadrature_points (2,2)%p(6))
FEM_quadrature_points (2,2)%p(1:6) = permutationStar21([1.0_pReal/6.0_pReal])
FEM_quadrature_points (2,2)%p = permutationStar21([1._pReal/6._pReal])
!--------------------------------------------------------------------------------------------------
! 2D cubic
FEM_nQuadrature(2,3) = 6
allocate(FEM_quadrature_weights(2,3)%p(6))
FEM_quadrature_weights(2,3)%p(1:3) = 0.22338158967801146570_pReal
FEM_quadrature_weights(2,3)%p(4:6) = 0.10995174365532186764_pReal
allocate(FEM_quadrature_weights(2,3)%p(FEM_nQuadrature(2,3)))
FEM_quadrature_weights(2,3)%p(1:3) = 2.2338158967801147e-1_pReal
FEM_quadrature_weights(2,3)%p(4:6) = 1.0995174365532187e-1_pReal
allocate(FEM_quadrature_points (2,3)%p(12))
FEM_quadrature_points (2,3)%p(1:6) = permutationStar21([0.44594849091596488632_pReal])
FEM_quadrature_points (2,3)%p(7:12)= permutationStar21([0.091576213509770743460_pReal])
FEM_quadrature_points (2,3)%p = [ &
permutationStar21([4.4594849091596489e-1_pReal]), &
permutationStar21([9.157621350977074e-2_pReal]) ]
!--------------------------------------------------------------------------------------------------
! 2D quartic
FEM_nQuadrature(2,4) = 12
allocate(FEM_quadrature_weights(2,4)%p(12))
FEM_quadrature_weights(2,4)%p(1:3) = 0.11678627572638_pReal
FEM_quadrature_weights(2,4)%p(4:6) = 0.05084490637021_pReal
FEM_quadrature_weights(2,4)%p(7:12) = 0.08285107561837_pReal
allocate(FEM_quadrature_weights(2,4)%p(FEM_nQuadrature(2,4)))
FEM_quadrature_weights(2,4)%p(1:3) = 1.1678627572637937e-1_pReal
FEM_quadrature_weights(2,4)%p(4:6) = 5.0844906370206817e-2_pReal
FEM_quadrature_weights(2,4)%p(7:12) = 8.285107561837358e-2_pReal
allocate(FEM_quadrature_points (2,4)%p(24))
FEM_quadrature_points (2,4)%p(1:6) = permutationStar21([0.24928674517091_pReal])
FEM_quadrature_points (2,4)%p(7:12) = permutationStar21([0.06308901449150_pReal])
FEM_quadrature_points (2,4)%p(13:24)= permutationStar111([0.31035245103378_pReal, 0.63650249912140_pReal])
FEM_quadrature_points (2,4)%p = [ &
permutationStar21([2.4928674517091042e-1_pReal]), &
permutationStar21([6.308901449150223e-2_pReal]), &
permutationStar111([3.1035245103378440e-1_pReal, 5.3145049844816947e-2_pReal]) ]
!--------------------------------------------------------------------------------------------------
! 2D quintic
FEM_nQuadrature(2,5) = 16
allocate(FEM_quadrature_weights(2,5)%p(16))
FEM_quadrature_weights(2,5)%p(1 ) = 0.14431560767779_pReal
FEM_quadrature_weights(2,5)%p(2:4) = 0.09509163426728_pReal
FEM_quadrature_weights(2,5)%p(5:7) = 0.10321737053472_pReal
FEM_quadrature_weights(2,5)%p(8:10) = 0.03245849762320_pReal
FEM_quadrature_weights(2,5)%p(11:16)= 0.02723031417443_pReal
allocate(FEM_quadrature_weights(2,5)%p(FEM_nQuadrature(2,5)))
FEM_quadrature_weights(2,5)%p(1:1) = 1.4431560767778717e-1_pReal
FEM_quadrature_weights(2,5)%p(2:4) = 9.509163426728463e-2_pReal
FEM_quadrature_weights(2,5)%p(5:7) = 1.0321737053471825e-1_pReal
FEM_quadrature_weights(2,5)%p(8:10) = 3.2458497623198080e-2_pReal
FEM_quadrature_weights(2,5)%p(11:16) = 2.7230314174434994e-2_pReal
allocate(FEM_quadrature_points (2,5)%p(32))
FEM_quadrature_points (2,5)%p(1:2) = permutationStar3([0.33333333333333_pReal])
FEM_quadrature_points (2,5)%p(3:8) = permutationStar21([0.45929258829272_pReal])
FEM_quadrature_points (2,5)%p(9:14) = permutationStar21([0.17056930775176_pReal])
FEM_quadrature_points (2,5)%p(15:20)= permutationStar21([0.05054722831703_pReal])
FEM_quadrature_points (2,5)%p(21:32)= permutationStar111([0.26311282963464_pReal, 0.72849239295540_pReal])
FEM_quadrature_points (2,5)%p = [ &
permutationStar3([1._pReal/3._pReal]), &
permutationStar21([4.5929258829272316e-1_pReal]), &
permutationStar21([1.705693077517602e-1_pReal]), &
permutationStar21([5.0547228317030975e-2_pReal]), &
permutationStar111([2.631128296346381e-1_pReal, 8.3947774099576053e-2_pReal]) ]
!--------------------------------------------------------------------------------------------------
! 3D linear
FEM_nQuadrature(3,1) = 1
allocate(FEM_quadrature_weights(3,1)%p(1))
FEM_quadrature_weights(3,1)%p(1) = 1.0_pReal
allocate(FEM_quadrature_weights(3,1)%p(FEM_nQuadrature(3,1)))
FEM_quadrature_weights(3,1)%p(1) = 1.0_pReal
allocate(FEM_quadrature_points (3,1)%p(3))
FEM_quadrature_points (3,1)%p(1:3)= permutationStar4([0.25_pReal])
FEM_quadrature_points (3,1)%p = permutationStar4([0.25_pReal])
!--------------------------------------------------------------------------------------------------
! 3D quadratic
FEM_nQuadrature(3,2) = 4
allocate(FEM_quadrature_weights(3,2)%p(4))
allocate(FEM_quadrature_weights(3,2)%p(FEM_nQuadrature(3,2)))
FEM_quadrature_weights(3,2)%p(1:4) = 0.25_pReal
allocate(FEM_quadrature_points (3,2)%p(12))
FEM_quadrature_points (3,2)%p(1:12)= permutationStar31([0.13819660112501051518_pReal])
FEM_quadrature_points (3,2)%p = permutationStar31([1.3819660112501052e-1_pReal])
!--------------------------------------------------------------------------------------------------
! 3D cubic
FEM_nQuadrature(3,3) = 14
allocate(FEM_quadrature_weights(3,3)%p(14))
FEM_quadrature_weights(3,3)%p(5:8) = 0.11268792571801585080_pReal
FEM_quadrature_weights(3,3)%p(1:4) = 0.073493043116361949544_pReal
FEM_quadrature_weights(3,3)%p(9:14) = 0.042546020777081466438_pReal
allocate(FEM_quadrature_weights(3,3)%p(FEM_nQuadrature(3,3)))
FEM_quadrature_weights(3,3)%p(1:4) = 7.3493043116361949e-2_pReal
FEM_quadrature_weights(3,3)%p(5:8) = 1.1268792571801585e-1_pReal
FEM_quadrature_weights(3,3)%p(9:14) = 4.2546020777081467e-2_pReal
allocate(FEM_quadrature_points (3,3)%p(42))
FEM_quadrature_points (3,3)%p(1:12) = permutationStar31([0.092735250310891226402_pReal])
FEM_quadrature_points (3,3)%p(13:24)= permutationStar31([0.31088591926330060980_pReal])
FEM_quadrature_points (3,3)%p(25:42)= permutationStar22([0.045503704125649649492_pReal])
FEM_quadrature_points (3,3)%p = [ &
permutationStar31([9.273525031089123e-2_pReal]), &
permutationStar31([3.108859192633006e-1_pReal]), &
permutationStar22([4.5503704125649649e-2_pReal]) ]
!--------------------------------------------------------------------------------------------------
! 3D quartic
! 3D quartic (lower precision/unknown source)
FEM_nQuadrature(3,4) = 35
allocate(FEM_quadrature_weights(3,4)%p(35))
FEM_quadrature_weights(3,4)%p(1:4) = 0.0021900463965388_pReal
FEM_quadrature_weights(3,4)%p(5:16) = 0.0143395670177665_pReal
FEM_quadrature_weights(3,4)%p(17:22) = 0.0250305395686746_pReal
FEM_quadrature_weights(3,4)%p(23:34) = 0.0479839333057554_pReal
FEM_quadrature_weights(3,4)%p(35) = 0.0931745731195340_pReal
allocate(FEM_quadrature_weights(3,4)%p(FEM_nQuadrature(3,4)))
FEM_quadrature_weights(3,4)%p(1:4) = 0.0021900463965388_pReal
FEM_quadrature_weights(3,4)%p(5:16) = 0.0143395670177665_pReal
FEM_quadrature_weights(3,4)%p(17:22) = 0.0250305395686746_pReal
FEM_quadrature_weights(3,4)%p(23:34) = 0.0479839333057554_pReal
FEM_quadrature_weights(3,4)%p(35) = 0.0931745731195340_pReal
allocate(FEM_quadrature_points (3,4)%p(105))
FEM_quadrature_points (3,4)%p(1:12) = permutationStar31([0.0267367755543735_pReal])
FEM_quadrature_points (3,4)%p(13:48) = permutationStar211([0.0391022406356488_pReal, 0.7477598884818090_pReal])
FEM_quadrature_points (3,4)%p(49:66) = permutationStar22([0.4547545999844830_pReal])
FEM_quadrature_points (3,4)%p(67:102) = permutationStar211([0.2232010379623150_pReal, 0.0504792790607720_pReal])
FEM_quadrature_points (3,4)%p(103:105)= permutationStar4([0.25_pReal])
FEM_quadrature_points (3,4)%p = [ &
permutationStar31([0.0267367755543735_pReal]), &
permutationStar211([0.0391022406356488_pReal, 0.7477598884818090_pReal]), &
permutationStar22([0.4547545999844830_pReal]), &
permutationStar211([0.2232010379623150_pReal, 0.0504792790607720_pReal]), &
permutationStar4([0.25_pReal]) ]
!--------------------------------------------------------------------------------------------------
! 3D quintic
! 3D quintic (lower precision/unknown source)
FEM_nQuadrature(3,5) = 56
allocate(FEM_quadrature_weights(3,5)%p(56))
FEM_quadrature_weights(3,5)%p(1:4) = 0.0010373112336140_pReal
FEM_quadrature_weights(3,5)%p(5:16) = 0.0096016645399480_pReal
FEM_quadrature_weights(3,5)%p(17:28) = 0.0164493976798232_pReal
FEM_quadrature_weights(3,5)%p(29:40) = 0.0153747766513310_pReal
FEM_quadrature_weights(3,5)%p(41:52) = 0.0293520118375230_pReal
FEM_quadrature_weights(3,5)%p(53:56) = 0.0366291366405108_pReal
allocate(FEM_quadrature_weights(3,5)%p(FEM_nQuadrature(3,5)))
FEM_quadrature_weights(3,5)%p(1:4) = 0.0010373112336140_pReal
FEM_quadrature_weights(3,5)%p(5:16) = 0.0096016645399480_pReal
FEM_quadrature_weights(3,5)%p(17:28) = 0.0164493976798232_pReal
FEM_quadrature_weights(3,5)%p(29:40) = 0.0153747766513310_pReal
FEM_quadrature_weights(3,5)%p(41:52) = 0.0293520118375230_pReal
FEM_quadrature_weights(3,5)%p(53:56) = 0.0366291366405108_pReal
allocate(FEM_quadrature_points (3,5)%p(168))
FEM_quadrature_points (3,5)%p(1:12) = permutationStar31([0.0149520651530592_pReal])
FEM_quadrature_points (3,5)%p(13:48) = permutationStar211([0.0340960211962615_pReal, 0.1518319491659370_pReal])
FEM_quadrature_points (3,5)%p(49:84) = permutationStar211([0.0462051504150017_pReal, 0.3549340560639790_pReal])
FEM_quadrature_points (3,5)%p(85:120) = permutationStar211([0.2281904610687610_pReal, 0.0055147549744775_pReal])
FEM_quadrature_points (3,5)%p(121:156)= permutationStar211([0.3523052600879940_pReal, 0.0992057202494530_pReal])
FEM_quadrature_points (3,5)%p(157:168)= permutationStar31([0.1344783347929940_pReal])
FEM_quadrature_points (3,5)%p = [ &
permutationStar31([0.0149520651530592_pReal]), &
permutationStar211([0.0340960211962615_pReal, 0.1518319491659370_pReal]), &
permutationStar211([0.0462051504150017_pReal, 0.3549340560639790_pReal]), &
permutationStar211([0.2281904610687610_pReal, 0.0055147549744775_pReal]), &
permutationStar211([0.3523052600879940_pReal, 0.0992057202494530_pReal]), &
permutationStar31([0.1344783347929940_pReal]) ]
call selfTest
end subroutine FEM_quadrature_init
@ -186,11 +191,9 @@ pure function permutationStar3(point) result(qPt)
real(pReal), dimension(2) :: qPt
real(pReal), dimension(1), intent(in) :: point
real(pReal), dimension(3,1) :: temp
temp(:,1) = [point(1), point(1), point(1)]
qPt = reshape(matmul(triangle, temp),[2])
qPt = pack(matmul(triangle,reshape([ &
point(1), point(1), point(1)],[3,1])),.true.)
end function permutationStar3
@ -203,13 +206,11 @@ pure function permutationStar21(point) result(qPt)
real(pReal), dimension(6) :: qPt
real(pReal), dimension(1), intent(in) :: point
real(pReal), dimension(3,3) :: temp
temp(:,1) = [point(1), point(1), 1.0_pReal - 2.0_pReal*point(1)]
temp(:,2) = [point(1), 1.0_pReal - 2.0_pReal*point(1), point(1)]
temp(:,3) = [1.0_pReal - 2.0_pReal*point(1), point(1), point(1)]
qPt = reshape(matmul(triangle, temp),[6])
qPt = pack(matmul(triangle,reshape([ &
point(1), point(1), 1.0_pReal - 2.0_pReal*point(1), &
point(1), 1.0_pReal - 2.0_pReal*point(1), point(1), &
1.0_pReal - 2.0_pReal*point(1), point(1), point(1)],[3,3])),.true.)
end function permutationStar21
@ -222,16 +223,14 @@ pure function permutationStar111(point) result(qPt)
real(pReal), dimension(12) :: qPt
real(pReal), dimension(2), intent(in) :: point
real(pReal), dimension(3,6) :: temp
temp(:,1) = [point(1), point(2), 1.0_pReal - point(1) - point(2)]
temp(:,2) = [point(1), 1.0_pReal - point(1) - point(2), point(2)]
temp(:,3) = [point(2), point(1), 1.0_pReal - point(1) - point(2)]
temp(:,4) = [point(2), 1.0_pReal - point(1) - point(2), point(1)]
temp(:,5) = [1.0_pReal - point(1) - point(2), point(2), point(1)]
temp(:,6) = [1.0_pReal - point(1) - point(2), point(1), point(2)]
qPt = reshape(matmul(triangle, temp),[12])
qPt = pack(matmul(triangle,reshape([ &
point(1), point(2), 1.0_pReal - point(1) - point(2), &
point(1), 1.0_pReal - point(1) - point(2), point(2), &
point(2), point(1), 1.0_pReal - point(1) - point(2), &
point(2), 1.0_pReal - point(1) - point(2), point(1), &
1.0_pReal - point(1) - point(2), point(2), point(1), &
1.0_pReal - point(1) - point(2), point(1), point(2)],[3,6])),.true.)
end function permutationStar111
@ -244,11 +243,9 @@ pure function permutationStar4(point) result(qPt)
real(pReal), dimension(3) :: qPt
real(pReal), dimension(1), intent(in) :: point
real(pReal), dimension(4,1) :: temp
temp(:,1) = [point(1), point(1), point(1), point(1)]
qPt = reshape(matmul(tetrahedron, temp),[3])
qPt = pack(matmul(tetrahedron,reshape([ &
point(1), point(1), point(1), point(1)],[4,1])),.true.)
end function permutationStar4
@ -261,14 +258,12 @@ pure function permutationStar31(point) result(qPt)
real(pReal), dimension(12) :: qPt
real(pReal), dimension(1), intent(in) :: point
real(pReal), dimension(4,4) :: temp
temp(:,1) = [point(1), point(1), point(1), 1.0_pReal - 3.0_pReal*point(1)]
temp(:,2) = [point(1), point(1), 1.0_pReal - 3.0_pReal*point(1), point(1)]
temp(:,3) = [point(1), 1.0_pReal - 3.0_pReal*point(1), point(1), point(1)]
temp(:,4) = [1.0_pReal - 3.0_pReal*point(1), point(1), point(1), point(1)]
qPt = reshape(matmul(tetrahedron, temp),[12])
qPt = pack(matmul(tetrahedron,reshape([ &
point(1), point(1), point(1), 1.0_pReal - 3.0_pReal*point(1), &
point(1), point(1), 1.0_pReal - 3.0_pReal*point(1), point(1), &
point(1), 1.0_pReal - 3.0_pReal*point(1), point(1), point(1), &
1.0_pReal - 3.0_pReal*point(1), point(1), point(1), point(1)],[4,4])),.true.)
end function permutationStar31
@ -276,21 +271,19 @@ end function permutationStar31
!--------------------------------------------------------------------------------------------------
!> @brief star 22 permutation of input
!--------------------------------------------------------------------------------------------------
pure function permutationStar22(point) result(qPt)
function permutationStar22(point) result(qPt)
real(pReal), dimension(18) :: qPt
real(pReal), dimension(1), intent(in) :: point
real(pReal), dimension(4,6) :: temp
temp(:,1) = [point(1), point(1), 0.5_pReal - point(1), 0.5_pReal - point(1)]
temp(:,2) = [point(1), 0.5_pReal - point(1), point(1), 0.5_pReal - point(1)]
temp(:,3) = [0.5_pReal - point(1), point(1), point(1), 0.5_pReal - point(1)]
temp(:,4) = [0.5_pReal - point(1), point(1), 0.5_pReal - point(1), point(1)]
temp(:,5) = [0.5_pReal - point(1), 0.5_pReal - point(1), point(1), point(1)]
temp(:,6) = [point(1), 0.5_pReal - point(1), 0.5_pReal - point(1), point(1)]
qPt = reshape(matmul(tetrahedron, temp),[18])
qPt = pack(matmul(tetrahedron,reshape([ &
point(1), point(1), 0.5_pReal - point(1), 0.5_pReal - point(1), &
point(1), 0.5_pReal - point(1), point(1), 0.5_pReal - point(1), &
0.5_pReal - point(1), point(1), point(1), 0.5_pReal - point(1), &
0.5_pReal - point(1), point(1), 0.5_pReal - point(1), point(1), &
0.5_pReal - point(1), 0.5_pReal - point(1), point(1), point(1), &
point(1), 0.5_pReal - point(1), 0.5_pReal - point(1), point(1)],[4,6])),.true.)
end function permutationStar22
@ -303,22 +296,20 @@ pure function permutationStar211(point) result(qPt)
real(pReal), dimension(36) :: qPt
real(pReal), dimension(2), intent(in) :: point
real(pReal), dimension(4,12) :: temp
temp(:,1 ) = [point(1), point(1), point(2), 1.0_pReal - 2.0_pReal*point(1) - point(2)]
temp(:,2 ) = [point(1), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(2)]
temp(:,3 ) = [point(1), point(2), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2)]
temp(:,4 ) = [point(1), point(2), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1)]
temp(:,5 ) = [point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(2)]
temp(:,6 ) = [point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(2), point(1)]
temp(:,7 ) = [point(2), point(1), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2)]
temp(:,8 ) = [point(2), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1)]
temp(:,9 ) = [point(2), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(1)]
temp(:,10) = [1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(1), point(2)]
temp(:,11) = [1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(2), point(1)]
temp(:,12) = [1.0_pReal - 2.0_pReal*point(1) - point(2), point(2), point(1), point(1)]
qPt = reshape(matmul(tetrahedron, temp),[36])
qPt = pack(matmul(tetrahedron,reshape([ &
point(1), point(1), point(2), 1.0_pReal - 2.0_pReal*point(1) - point(2), &
point(1), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(2), &
point(1), point(2), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), &
point(1), point(2), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), &
point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(2), &
point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(2), point(1), &
point(2), point(1), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), &
point(2), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), &
point(2), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(1), &
1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(1), point(2), &
1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(2), point(1), &
1.0_pReal - 2.0_pReal*point(1) - point(2), point(2), point(1), point(1)],[4,12])),.true.)
end function permutationStar211
@ -331,35 +322,60 @@ pure function permutationStar1111(point) result(qPt)
real(pReal), dimension(72) :: qPt
real(pReal), dimension(3), intent(in) :: point
real(pReal), dimension(4,24) :: temp
temp(:,1 ) = [point(1), point(2), point(3), 1.0_pReal - point(1) - point(2)- point(3)]
temp(:,2 ) = [point(1), point(2), 1.0_pReal - point(1) - point(2)- point(3), point(3)]
temp(:,3 ) = [point(1), point(3), point(2), 1.0_pReal - point(1) - point(2)- point(3)]
temp(:,4 ) = [point(1), point(3), 1.0_pReal - point(1) - point(2)- point(3), point(2)]
temp(:,5 ) = [point(1), 1.0_pReal - point(1) - point(2)- point(3), point(2), point(3)]
temp(:,6 ) = [point(1), 1.0_pReal - point(1) - point(2)- point(3), point(3), point(2)]
temp(:,7 ) = [point(2), point(1), point(3), 1.0_pReal - point(1) - point(2)- point(3)]
temp(:,8 ) = [point(2), point(1), 1.0_pReal - point(1) - point(2)- point(3), point(3)]
temp(:,9 ) = [point(2), point(3), point(1), 1.0_pReal - point(1) - point(2)- point(3)]
temp(:,10) = [point(2), point(3), 1.0_pReal - point(1) - point(2)- point(3), point(1)]
temp(:,11) = [point(2), 1.0_pReal - point(1) - point(2)- point(3), point(1), point(3)]
temp(:,12) = [point(2), 1.0_pReal - point(1) - point(2)- point(3), point(3), point(1)]
temp(:,13) = [point(3), point(1), point(2), 1.0_pReal - point(1) - point(2)- point(3)]
temp(:,14) = [point(3), point(1), 1.0_pReal - point(1) - point(2)- point(3), point(2)]
temp(:,15) = [point(3), point(2), point(1), 1.0_pReal - point(1) - point(2)- point(3)]
temp(:,16) = [point(3), point(2), 1.0_pReal - point(1) - point(2)- point(3), point(1)]
temp(:,17) = [point(3), 1.0_pReal - point(1) - point(2)- point(3), point(1), point(2)]
temp(:,18) = [point(3), 1.0_pReal - point(1) - point(2)- point(3), point(2), point(1)]
temp(:,19) = [1.0_pReal - point(1) - point(2)- point(3), point(1), point(2), point(3)]
temp(:,20) = [1.0_pReal - point(1) - point(2)- point(3), point(1), point(3), point(2)]
temp(:,21) = [1.0_pReal - point(1) - point(2)- point(3), point(2), point(1), point(3)]
temp(:,22) = [1.0_pReal - point(1) - point(2)- point(3), point(2), point(3), point(1)]
temp(:,23) = [1.0_pReal - point(1) - point(2)- point(3), point(3), point(1), point(2)]
temp(:,24) = [1.0_pReal - point(1) - point(2)- point(3), point(3), point(2), point(1)]
qPt = reshape(matmul(tetrahedron, temp),[72])
qPt = pack(matmul(tetrahedron,reshape([ &
point(1), point(2), point(3), 1.0_pReal - point(1) - point(2)- point(3), &
point(1), point(2), 1.0_pReal - point(1) - point(2)- point(3), point(3), &
point(1), point(3), point(2), 1.0_pReal - point(1) - point(2)- point(3), &
point(1), point(3), 1.0_pReal - point(1) - point(2)- point(3), point(2), &
point(1), 1.0_pReal - point(1) - point(2)- point(3), point(2), point(3), &
point(1), 1.0_pReal - point(1) - point(2)- point(3), point(3), point(2), &
point(2), point(1), point(3), 1.0_pReal - point(1) - point(2)- point(3), &
point(2), point(1), 1.0_pReal - point(1) - point(2)- point(3), point(3), &
point(2), point(3), point(1), 1.0_pReal - point(1) - point(2)- point(3), &
point(2), point(3), 1.0_pReal - point(1) - point(2)- point(3), point(1), &
point(2), 1.0_pReal - point(1) - point(2)- point(3), point(1), point(3), &
point(2), 1.0_pReal - point(1) - point(2)- point(3), point(3), point(1), &
point(3), point(1), point(2), 1.0_pReal - point(1) - point(2)- point(3), &
point(3), point(1), 1.0_pReal - point(1) - point(2)- point(3), point(2), &
point(3), point(2), point(1), 1.0_pReal - point(1) - point(2)- point(3), &
point(3), point(2), 1.0_pReal - point(1) - point(2)- point(3), point(1), &
point(3), 1.0_pReal - point(1) - point(2)- point(3), point(1), point(2), &
point(3), 1.0_pReal - point(1) - point(2)- point(3), point(2), point(1), &
1.0_pReal - point(1) - point(2)- point(3), point(1), point(2), point(3), &
1.0_pReal - point(1) - point(2)- point(3), point(1), point(3), point(2), &
1.0_pReal - point(1) - point(2)- point(3), point(2), point(1), point(3), &
1.0_pReal - point(1) - point(2)- point(3), point(2), point(3), point(1), &
1.0_pReal - point(1) - point(2)- point(3), point(3), point(1), point(2), &
1.0_pReal - point(1) - point(2)- point(3), point(3), point(2), point(1)],[4,24])),.true.)
end function permutationStar1111
!--------------------------------------------------------------------------------------------------
!> @brief Check correctness of quadrature weights and points.
!--------------------------------------------------------------------------------------------------
subroutine selfTest
integer :: o, d, n
real(pReal), dimension(2:3), parameter :: w = [3.0_pReal,2.0_pReal]
do d = lbound(FEM_quadrature_weights,1), ubound(FEM_quadrature_weights,1)
do o = lbound(FEM_quadrature_weights(d,:),1), ubound(FEM_quadrature_weights(d,:),1)
if (dNeq(sum(FEM_quadrature_weights(d,o)%p),1.0_pReal,5e-15_pReal)) &
error stop 'quadrature weights'
enddo
enddo
do d = lbound(FEM_quadrature_points,1), ubound(FEM_quadrature_points,1)
do o = lbound(FEM_quadrature_points(d,:),1), ubound(FEM_quadrature_points(d,:),1)
n = size(FEM_quadrature_points(d,o)%p,1)/d
if (any(dNeq(sum(reshape(FEM_quadrature_points(d,o)%p,[d,n]),2),-real(n,pReal)/w(d),1.e-14_pReal))) &
error stop 'quadrature points'
enddo
enddo
end subroutine selfTest
end module FEM_quadrature

35
src/mesh/FEM_utilities.f90
View File

@ -23,14 +23,8 @@ module FEM_utilities
implicit none
private
!--------------------------------------------------------------------------------------------------
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 :: nActiveFields = 0
!--------------------------------------------------------------------------------------------------
! grid related information information
real(pReal), public :: wgt !< weighting factor 1/Nelems
logical, public :: cutBack = .false. !< cut back of BVP solver in case convergence is not achieved or a material point is terminally ill
real(pReal), public, protected :: wgt !< weighting factor 1/Nelems
!--------------------------------------------------------------------------------------------------
@ -49,10 +43,6 @@ module FEM_utilities
COMPONENT_MECH_Z_ID
end enum
!--------------------------------------------------------------------------------------------------
! variables controlling debugging
logical :: &
debugPETSc !< use some in debug defined options for more verbose PETSc solution
!--------------------------------------------------------------------------------------------------
! derived types
@ -63,27 +53,17 @@ module FEM_utilities
end type tSolutionState
type, public :: tComponentBC
integer(kind(COMPONENT_UNDEFINED_ID)) :: ID
real(pReal), allocatable, dimension(:) :: Value
logical, allocatable, dimension(:) :: Mask
integer(kind(COMPONENT_UNDEFINED_ID)) :: ID
real(pReal), allocatable, dimension(:) :: Value
logical, allocatable, dimension(:) :: Mask
end type tComponentBC
type, public :: tFieldBC
integer(kind(FIELD_UNDEFINED_ID)) :: ID
integer :: nComponents = 0
type(tComponentBC), allocatable :: componentBC(:)
type(tComponentBC), allocatable, dimension(:) :: componentBC
end type tFieldBC
type, public :: tLoadCase
real(pReal) :: time = 0.0_pReal !< length of increment
integer :: incs = 0, & !< number of increments
outputfrequency = 1, & !< frequency of result writes
logscale = 0 !< linear/logarithmic time inc flag
logical :: followFormerTrajectory = .true. !< follow trajectory of former loadcase
integer, allocatable, dimension(:) :: faceID
type(tFieldBC), allocatable, dimension(:) :: fieldBC
end type tLoadCase
public :: &
FEM_utilities_init, &
utilities_constitutiveResponse, &
@ -109,8 +89,9 @@ subroutine FEM_utilities_init
integer :: structOrder !< order of displacement shape functions
character(len=*), parameter :: &
PETSCDEBUG = ' -snes_view -snes_monitor '
PetscErrorCode :: ierr
logical :: debugPETSc !< use some in debug defined options for more verbose PETSc solution
print'(/,a)', ' <<<+- FEM_utilities init -+>>>'

16
src/mesh/discretization_mesh.f90
View File

@ -40,6 +40,11 @@ module discretization_mesh
mesh_maxNips !< max number of IPs in any CP element
!!!! BEGIN DEPRECATED !!!!!
DM, public :: geomMesh
PetscInt, dimension(:), allocatable, public, protected :: &
mesh_boundaries
real(pReal), dimension(:,:), allocatable :: &
mesh_ipVolume, & !< volume associated with IP (initially!)
mesh_node0 !< node x,y,z coordinates (initially!)
@ -50,11 +55,6 @@ module discretization_mesh
real(pReal), dimension(:,:,:), allocatable :: &
mesh_ipCoordinates !< IP x,y,z coordinates (after deformation!)
DM, public :: geomMesh
PetscInt, dimension(:), allocatable, public, protected :: &
mesh_boundaries
public :: &
discretization_mesh_init, &
mesh_FEM_build_ipVolumes, &
@ -71,16 +71,14 @@ subroutine discretization_mesh_init(restart)
logical, intent(in) :: restart
integer, allocatable, dimension(:) :: chunkPos
integer :: dimPlex, &
mesh_Nnodes, & !< total number of nodes in mesh
j, l, &
j, &
debug_element, debug_ip
PetscSF :: sf
DM :: globalMesh
PetscInt :: nFaceSets
PetscInt, pointer, dimension(:) :: pFaceSets
character(len=pStringLen), dimension(:), allocatable :: fileContent
IS :: faceSetIS
PetscErrorCode :: ierr
integer, dimension(:), allocatable :: &
@ -88,7 +86,7 @@ subroutine discretization_mesh_init(restart)
class(tNode), pointer :: &
num_mesh
integer :: integrationOrder !< order of quadrature rule required
type(tvec) :: coords_node0
type(tvec) :: coords_node0
print'(/,a)', ' <<<+- discretization_mesh init -+>>>'

54
src/mesh/mesh_mech_FEM.f90
View File

@ -109,7 +109,7 @@ subroutine FEM_mechanical_init(fieldBC)
character(len=*), parameter :: prefix = 'mechFE_'
PetscErrorCode :: ierr
real(pReal), dimension(3,3) :: devNull
class(tNode), pointer :: &
num_mesh
@ -258,6 +258,7 @@ subroutine FEM_mechanical_init(fieldBC)
call DMPlexVecSetClosure(mechanical_mesh,section,solution_local,cell,px_scal,5,ierr)
CHKERRQ(ierr)
enddo
call utilities_constitutiveResponse(0.0_pReal,devNull,.true.)
end subroutine FEM_mechanical_init
@ -288,8 +289,8 @@ type(tSolutionState) function FEM_mechanical_solution( &
params%timeinc = timeinc
params%fieldBC = fieldBC
call SNESSolve(mechanical_snes,PETSC_NULL_VEC,solution,ierr); CHKERRQ(ierr) ! solve mechanical_snes based on solution guess (result in solution)
call SNESGetConvergedReason(mechanical_snes,reason,ierr); CHKERRQ(ierr) ! solution converged?
call SNESSolve(mechanical_snes,PETSC_NULL_VEC,solution,ierr); CHKERRQ(ierr) ! solve mechanical_snes based on solution guess (result in solution)
call SNESGetConvergedReason(mechanical_snes,reason,ierr); CHKERRQ(ierr) ! solution converged?
terminallyIll = .false.
if (reason < 1) then ! 0: still iterating (will not occur), negative -> convergence error
@ -397,7 +398,7 @@ subroutine FEM_mechanical_formResidual(dm_local,xx_local,f_local,dummy,ierr)
!--------------------------------------------------------------------------------------------------
! evaluate constitutive response
call Utilities_constitutiveResponse(params%timeinc,P_av,ForwardData)
call utilities_constitutiveResponse(params%timeinc,P_av,ForwardData)
call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,ierr)
ForwardData = .false.
@ -670,7 +671,7 @@ end subroutine FEM_mechanical_converged
!--------------------------------------------------------------------------------------------------
!> @brief Calculate current coordinates (FEM nodal coordinates only at the moment)
!> @brief Calculate current coordinates (both nodal and ip coordinates)
!--------------------------------------------------------------------------------------------------
subroutine FEM_mechanical_updateCoords()
@ -678,21 +679,35 @@ subroutine FEM_mechanical_updateCoords()
nodeCoords_linear !< nodal coordinates (dimPlex*Nnodes)
real(pReal), pointer, dimension(:,:) :: &
nodeCoords !< nodal coordinates (3,Nnodes)
real(pReal), pointer, dimension(:,:,:) :: &
ipCoords !< ip coordinates (3,nQuadrature,mesh_NcpElems)
integer :: &
qPt, &
comp, &
qOffset, &
nOffset
DM :: dm_local
Vec :: x_local
PetscErrorCode :: ierr
PetscInt :: dimPlex, pStart, pEnd, p, s, e
PetscInt :: pStart, pEnd, p, s, e, q, &
cellStart, cellEnd, c, n
PetscSection :: section
PetscQuadrature :: mechQuad
PetscReal, dimension(:), pointer :: basisField, basisFieldDer
PetscScalar, dimension(:), pointer :: x_scal
call SNESGetDM(mechanical_snes,dm_local,ierr); CHKERRQ(ierr)
call DMGetDS(dm_local,mechQuad,ierr); CHKERRQ(ierr)
call DMGetLocalSection(dm_local,section,ierr); CHKERRQ(ierr)
call DMGetLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr)
call DMGetDimension(dm_local,dimPlex,ierr); CHKERRQ(ierr)
! write cell vertex displacements
call DMPlexGetDepthStratum(dm_local,0,pStart,pEnd,ierr); CHKERRQ(ierr)
allocate(nodeCoords(3,pStart:pEnd-1),source=0.0_pReal)
call VecGetArrayF90(x_local,nodeCoords_linear,ierr); CHKERRQ(ierr)
do p=pStart, pEnd-1
call DMPlexGetPointLocal(dm_local, p, s, e, ierr); CHKERRQ(ierr)
nodeCoords(1:dimPlex,p)=nodeCoords_linear(s+1:e)
@ -700,6 +715,31 @@ subroutine FEM_mechanical_updateCoords()
call discretization_setNodeCoords(nodeCoords)
call VecRestoreArrayF90(x_local,nodeCoords_linear,ierr); CHKERRQ(ierr)
! write ip displacements
call DMPlexGetHeightStratum(dm_local,0,cellStart,cellEnd,ierr); CHKERRQ(ierr)
call PetscDSGetTabulation(mechQuad,0,basisField,basisFieldDer,ierr); CHKERRQ(ierr)
allocate(ipCoords(3,nQuadrature,mesh_NcpElems),source=0.0_pReal)
do c=cellStart,cellEnd-1
qOffset=0
call DMPlexVecGetClosure(dm_local,section,x_local,c,x_scal,ierr); CHKERRQ(ierr) !< get nodal coordinates of each element
do qPt=0,nQuadrature-1
qOffset= qPt * (size(basisField)/nQuadrature)
do comp=0,dimPlex-1 !< loop over components
nOffset=0
q = comp
do n=0,nBasis-1
ipCoords(comp+1,qPt+1,c+1)=ipCoords(comp+1,qPt+1,c+1)+&
sum(basisField(qOffset+(q*dimPlex)+1:qOffset+(q*dimPlex)+dimPlex)*&
x_scal(nOffset+1:nOffset+dimPlex))
q = q+dimPlex
nOffset = nOffset+dimPlex
enddo
enddo
enddo
call DMPlexVecRestoreClosure(dm_local,section,x_local,c,x_scal,ierr); CHKERRQ(ierr)
end do
call discretization_setIPcoords(reshape(ipCoords,[3,mesh_NcpElems*nQuadrature]))
call DMRestoreLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr)
end subroutine FEM_mechanical_updateCoords

30
src/phase_mechanical.f90
View File

@ -229,8 +229,8 @@ module subroutine mechanical_init(materials,phases)
allocate(phase_mechanical_F0(phases%length))
allocate(phase_mechanical_Li(phases%length))
allocate(phase_mechanical_Li0(phases%length))
allocate(phase_mechanical_Lp0(phases%length))
allocate(phase_mechanical_Lp(phases%length))
allocate(phase_mechanical_Lp0(phases%length))
allocate(phase_mechanical_S(phases%length))
allocate(phase_mechanical_P(phases%length))
allocate(phase_mechanical_S0(phases%length))
@ -238,20 +238,20 @@ module subroutine mechanical_init(materials,phases)
do ph = 1, phases%length
Nmembers = count(material_phaseID == ph)
allocate(phase_mechanical_Fi(ph)%data(3,3,Nmembers))
allocate(phase_mechanical_Fe(ph)%data(3,3,Nmembers))
allocate(phase_mechanical_Fi(ph)%data(3,3,Nmembers))
allocate(phase_mechanical_Fi0(ph)%data(3,3,Nmembers))
allocate(phase_mechanical_Fp(ph)%data(3,3,Nmembers))
allocate(phase_mechanical_Fp0(ph)%data(3,3,Nmembers))
allocate(phase_mechanical_Li(ph)%data(3,3,Nmembers))
allocate(phase_mechanical_Li0(ph)%data(3,3,Nmembers))
allocate(phase_mechanical_Lp0(ph)%data(3,3,Nmembers))
allocate(phase_mechanical_Lp(ph)%data(3,3,Nmembers))
allocate(phase_mechanical_F(ph)%data(3,3,Nmembers))
allocate(phase_mechanical_F0(ph)%data(3,3,Nmembers))
allocate(phase_mechanical_Li(ph)%data(3,3,Nmembers),source=0.0_pReal)
allocate(phase_mechanical_Li0(ph)%data(3,3,Nmembers),source=0.0_pReal)
allocate(phase_mechanical_Lp(ph)%data(3,3,Nmembers),source=0.0_pReal)
allocate(phase_mechanical_Lp0(ph)%data(3,3,Nmembers),source=0.0_pReal)
allocate(phase_mechanical_S(ph)%data(3,3,Nmembers),source=0.0_pReal)
allocate(phase_mechanical_P(ph)%data(3,3,Nmembers),source=0.0_pReal)
allocate(phase_mechanical_S0(ph)%data(3,3,Nmembers),source=0.0_pReal)
allocate(phase_mechanical_F(ph)%data(3,3,Nmembers))
allocate(phase_mechanical_F0(ph)%data(3,3,Nmembers))
phase => phases%get(ph)
mech => phase%get('mechanical')
@ -508,7 +508,7 @@ function integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) result(broken)
enddo LpLoop
call phase_LiAndItsTangents(Li_constitutive, dLi_dS, dLi_dFi, &
S, Fi_new, ph,en)
S, Fi_new, ph,en)
!* update current residuum and check for convergence of loop
atol_Li = max(num%rtol_crystalliteStress * max(norm2(Liguess),norm2(Li_constitutive)), & ! absolute tolerance from largest acceptable relative error
@ -1021,7 +1021,7 @@ module function phase_mechanical_constitutive(Delta_t,co,ip,el) result(converged
subLi0 = phase_mechanical_Li0(ph)%data(1:3,1:3,en)
subLp0 = phase_mechanical_Lp0(ph)%data(1:3,1:3,en)
subState0 = plasticState(ph)%State0(:,en)
allocate(subState0,source=plasticState(ph)%State0(:,en))
subFp0 = phase_mechanical_Fp0(ph)%data(1:3,1:3,en)
subFi0 = phase_mechanical_Fi0(ph)%data(1:3,1:3,en)
subF0 = phase_mechanical_F0(ph)%data(1:3,1:3,en)
@ -1144,12 +1144,12 @@ module function phase_mechanical_dPdF(Delta_t,co,ce) result(dPdF)
en = material_phaseEntry(co,ce)
call phase_hooke_SandItsTangents(devNull,dSdFe,dSdFi, &
phase_mechanical_Fe(ph)%data(1:3,1:3,en), &
phase_mechanical_Fi(ph)%data(1:3,1:3,en),ph,en)
phase_mechanical_Fe(ph)%data(1:3,1:3,en), &
phase_mechanical_Fi(ph)%data(1:3,1:3,en),ph,en)
call phase_LiAndItsTangents(devNull,dLidS,dLidFi, &
phase_mechanical_S(ph)%data(1:3,1:3,en), &
phase_mechanical_Fi(ph)%data(1:3,1:3,en), &
ph,en)
phase_mechanical_S(ph)%data(1:3,1:3,en), &
phase_mechanical_Fi(ph)%data(1:3,1:3,en), &
ph,en)
invFp = math_inv33(phase_mechanical_Fp(ph)%data(1:3,1:3,en))
invFi = math_inv33(phase_mechanical_Fi(ph)%data(1:3,1:3,en))

5
src/prec.f90
View File

@ -27,11 +27,6 @@ module prec
real(pReal), parameter :: tol_math_check = 1.0e-8_pReal !< tolerance for internal math self-checks (rotation)
type :: group_float !< variable length datatype used for storage of state
real(pReal), dimension(:), pointer :: p
end type group_float
type :: tState
integer :: &
sizeState = 0, & !< size of state

5
src/quit.f90
View File

@ -6,7 +6,10 @@
!--------------------------------------------------------------------------------------------------
subroutine quit(stop_id)
#include <petsc/finclude/petscsys.h>
use PetscSys
use PETScSys
#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY)
use MPI_f08
#endif
use HDF5
implicit none