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

WIP: adopting to PETSc 3.9.x and modifications in development branch

This commit is contained in:
Martin Diehl 2018-08-17 11:23:24 +02:00
parent 821860987c
commit d4bcfae82b
8 changed files with 614 additions and 346 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

9
src/CMakeLists.txt
View File

@ -57,7 +57,7 @@ if (PROJECT_NAME STREQUAL "DAMASK_spectral")
add_dependencies(MESH DAMASK_MATH)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:MESH>)
elseif (PROJECT_NAME STREQUAL "DAMASK_FEM")
add_library(FEZoo OBJECT "FEZoo.f90")
add_library(FEZoo OBJECT "FEM_zoo.f90")
add_dependencies(FEZoo DAMASK_MATH)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:FEZoo>)
add_library(MESH OBJECT "meshFEM.f90")
@ -186,14 +186,9 @@ elseif (PROJECT_NAME STREQUAL "DAMASK_FEM")
add_dependencies(FEM_UTILITIES DAMASK_CPFE)
add_library(FEM_SOLVER OBJECT
"FEM_hydrogenflux.f90"
"FEM_porosity.f90"
"FEM_vacancyflux.f90"
"FEM_damage.f90"
"FEM_thermal.f90"
"FEM_mech.f90")
add_dependencies(FEM_SOLVER FEM_UTILITIES)
add_executable(DAMASK_FEM "DAMASK_FEM_driver.f90")
add_executable(DAMASK_FEM "DAMASK_FEM.f90")
add_dependencies(DAMASK_FEM FEM_SOLVER)
endif()

6
src/CPFEM2.f90
View File

@ -50,8 +50,8 @@ subroutine CPFEM_initAll(el,ip)
IO_init
use DAMASK_interface
#ifdef FEM
use FEZoo, only: &
FEZoo_init
use FEM_Zoo, only: &
FEM_Zoo_init
#endif
implicit none
@ -62,7 +62,7 @@ subroutine CPFEM_initAll(el,ip)
call prec_init
call IO_init
#ifdef FEM
call FEZoo_init
call FEM_Zoo_init
#endif
call numerics_init
call debug_init

2
src/FEM_interface.f90
View File

@ -210,7 +210,7 @@ subroutine DAMASK_interface_init()
write(6,'(a,a)') ' Geometry file: ', trim(geometryFile)
write(6,'(a,a)') ' Loadcase file: ', trim(loadCaseFile)
write(6,'(a,a)') ' Solver job name: ', trim(getSolverJobName())
if (SpectralRestartInc > 0_pInt) &
if (FEMRestartInc > 0_pInt) &
write(6,'(a,i6.6)') ' Restart from increment: ', FEMRestartInc
write(6,'(a,l1,/)') ' Append to result file: ', appendToOutFile

2
src/FEM_mech.f90 Executable file → Normal file
View File

@ -23,7 +23,7 @@ module FEM_mech
implicit none
private
#include <petsc/finclude/petsc.h90>
#include <petsc/finclude/petsc.h>
!--------------------------------------------------------------------------------------------------
! derived types

361
src/FEM_utilities.f90
View File

@ -3,14 +3,16 @@
!> @brief Utilities used by the FEM solver
!--------------------------------------------------------------------------------------------------
module FEM_utilities
use, intrinsic :: iso_c_binding
use prec, only: &
pReal, &
pInt
#include <petsc/finclude/petscis.h>
#include <petsc/finclude/petscdmda.h>
use prec, only: pReal, pInt
use PETScdmda
use PETScis
implicit none
private
#include <petsc/finclude/petsc.h90>
#include <petsc/finclude/petsc.h>
!--------------------------------------------------------------------------------------------------
!
logical, public :: cutBack = .false. !< cut back of BVP solver in case convergence is not achieved or a material point is terminally ill
@ -141,36 +143,13 @@ module FEM_utilities
COMPONENT_MGTWIN_PHI_ID
external :: &
MPI_abort, &
MPI_Allreduce, &
PetscOptionsClear, &
PetscOptionsInsertString, &
PetscObjectSetName, &
VecCreateMPI, &
VecSetFromOptions, &
VecGetSize, &
VecAssemblyBegin, &
VecAssemblyEnd, &
VecView, &
VecDestroy, &
ISCreateGeneral, &
ISDuplicate, &
ISDifference, &
ISGetSize, &
ISLocalToGlobalMappingApplyIS, &
ISDestroy, &
DMGetDimension, &
DMGetLocalToGlobalMapping, &
DMGetLabel, &
DMGetStratumSize, &
DMGetStratumIS, &
DMPlexGetHeightStratum, &
DMGetLabelIdIS, &
DMPlexGetChart, &
DMPlexLabelComplete, &
PetscSectionGetStorageSize, &
PetscSectionGetFieldDof, &
PetscSectionGetFieldOffset, &
PetscViewerHDF5Open, &
PetscViewerHDF5PushGroup, &
PetscViewerHDF5PopGroup, &
@ -195,12 +174,7 @@ subroutine utilities_init()
worldsize, &
worldrank, &
petsc_defaultOptions, &
petsc_options, &
structOrder, &
thermalOrder, &
damageOrder, &
soluteOrder, &
mgtwinOrder
petsc_options
use debug, only: &
debug_level, &
debug_SPECTRAL, &
@ -215,22 +189,12 @@ subroutine utilities_init()
mesh_maxNips, &
geomMesh, &
mesh_element
use homogenization, only: &
homogOutput, &
crystalliteOutput, &
phaseOutput
use material, only: &
material_Nhomogenization, &
material_Ncrystallite, &
material_Nphase, &
homogenization_Ngrains, &
homogenization_maxNgrains, &
material_homog, &
material_phase, &
microstructure_crystallite, &
homogenization_name, &
crystallite_name, &
phase_name
microstructure_crystallite
implicit none
@ -262,27 +226,15 @@ subroutine utilities_init()
trim(PETScDebug), &
' add more using the PETSc_Options keyword in numerics.config '
flush(6)
call PetscOptionsClear(PETSC_NULL_OBJECT,ierr)
call PetscOptionsClear(PETSC_NULL_OPTIONS,ierr)
CHKERRQ(ierr)
if(debugPETSc) call PetscOptionsInsertString(PETSC_NULL_OBJECT,trim(PETSCDEBUG),ierr)
if(debugPETSc) call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(PETSCDEBUG),ierr)
CHKERRQ(ierr)
call PetscOptionsInsertString(PETSC_NULL_OBJECT,trim(petsc_defaultOptions),ierr)
call PetscOptionsInsertString(PETSC_NULL_OBJECT,trim(petsc_options),ierr)
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
call PetscOptionsInsertString(PETSC_NULL_OBJECT,trim(petsc_optionsPhysics),ierr)
CHKERRQ(ierr)
write(petsc_optionsPhysics,'(a,i0)') '-thermalFE_petscspace_order ', thermalOrder
call PetscOptionsInsertString(PETSC_NULL_OBJECT,trim(petsc_optionsPhysics),ierr)
CHKERRQ(ierr)
write(petsc_optionsPhysics,'(a,i0)') '-damageFE_petscspace_order ' , damageOrder
call PetscOptionsInsertString(PETSC_NULL_OBJECT,trim(petsc_optionsPhysics),ierr)
CHKERRQ(ierr)
write(petsc_optionsPhysics,'(a,i0)') '-soluteFE_petscspace_order ', soluteOrder
call PetscOptionsInsertString(PETSC_NULL_OBJECT,trim(petsc_optionsPhysics),ierr)
CHKERRQ(ierr)
write(petsc_optionsPhysics,'(a,i0)') '-mgtwinFE_petscspace_order ', mgtwinOrder
call PetscOptionsInsertString(PETSC_NULL_OBJECT,trim(petsc_optionsPhysics),ierr)
!write(petsc_optionsPhysics,'(a,i0)') '-mechFE_petscspace_order ' , structOrder
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_optionsPhysics),ierr)
CHKERRQ(ierr)
wgt = 1.0/real(mesh_maxNips*mesh_NcpElemsGlobal,pReal)
@ -368,129 +320,126 @@ subroutine utilities_init()
call PetscObjectSetName(coordinatesVec, 'NodalCoordinates',ierr)
call VecSetFromOptions(coordinatesVec, ierr); CHKERRQ(ierr)
allocate(mappingCells(worldsize), source = 0)
allocate(homogenizationResultsVec(material_Nhomogenization ))
allocate(crystalliteResultsVec (material_Ncrystallite, homogenization_maxNgrains))
allocate(phaseResultsVec (material_Nphase, homogenization_maxNgrains))
do homog = 1, material_Nhomogenization
mappingCells = 0_pInt; mappingCells(worldrank+1) = homogOutput(homog)%sizeIpCells
call MPI_Allreduce(MPI_IN_PLACE,mappingCells,worldsize,MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr)
call VecCreateMPI(PETSC_COMM_WORLD,mappingCells(worldrank+1),sum(mappingCells), &
homogenizationResultsVec(homog),ierr);CHKERRQ(ierr)
if (sum(mappingCells) > 0) then
call VecCreateMPI(PETSC_COMM_WORLD,mappingCells(worldrank+1)*2**dimPlex,sum(mappingCells)*2**dimPlex, &
connectivityVec,ierr);CHKERRQ(ierr)
call PetscObjectSetName(connectivityVec,'mapping_'//trim(homogenization_name(homog)),ierr)
CHKERRQ(ierr)
call VecGetArrayF90(connectivityVec,results,ierr); CHKERRQ(ierr)
results = 0.0_pReal; ctr = 1_pInt
do cell = cellStart, cellEnd-1; do qPt = 1, mesh_maxNips
if (material_homog(qPt,cell+1) == homog) then
results(ctr:ctr+2**dimPlex-1) = real(reshape(connectivity(1:2**dimPlex,mesh_maxNips*cell+qPt), &
shape=[2**dimPlex]))
ctr = ctr + 2**dimPlex
endif
enddo; enddo
call VecRestoreArrayF90(connectivityVec, results, ierr); CHKERRQ(ierr)
call VecAssemblyBegin(connectivityVec, ierr); CHKERRQ(ierr)
call VecAssemblyEnd (connectivityVec, ierr); CHKERRQ(ierr)
call VecView(connectivityVec, resUnit, ierr); CHKERRQ(ierr)
call VecDestroy(connectivityVec, ierr); CHKERRQ(ierr)
endif
enddo
do cryst = 1, material_Ncrystallite; do grain = 1, homogenization_maxNgrains
mappingCells = 0_pInt
mappingCells(worldrank+1) = crystalliteOutput(cryst,grain)%sizeIpCells
call MPI_Allreduce(MPI_IN_PLACE,mappingCells,worldsize,MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr)
call VecCreateMPI(PETSC_COMM_WORLD,mappingCells(worldrank+1),sum(mappingCells), &
crystalliteResultsVec(cryst,grain),ierr);CHKERRQ(ierr)
if (sum(mappingCells) > 0) then
call VecCreateMPI(PETSC_COMM_WORLD,mappingCells(worldrank+1)*2**dimPlex,sum(mappingCells)*2**dimPlex, &
connectivityVec,ierr);CHKERRQ(ierr)
write(grainStr,'(a,i0)') 'Grain',grain
call PetscObjectSetName(connectivityVec,'mapping_'// &
trim(crystallite_name(cryst))//'_'// &
trim(grainStr),ierr)
CHKERRQ(ierr)
call VecGetArrayF90(connectivityVec, results, ierr); CHKERRQ(ierr)
results = 0.0_pReal; ctr = 1_pInt
do cell = cellStart, cellEnd-1; do qPt = 1, mesh_maxNips
if (homogenization_Ngrains (mesh_element(3,cell+1)) >= grain .and. &
microstructure_crystallite(mesh_element(4,cell+1)) == cryst) then
results(ctr:ctr+2**dimPlex-1) = real(reshape(connectivity(1:2**dimPlex,mesh_maxNips*cell+qPt), &
shape=[2**dimPlex]))
ctr = ctr + 2**dimPlex
endif
enddo; enddo
call VecRestoreArrayF90(connectivityVec, results, ierr); CHKERRQ(ierr)
call VecAssemblyBegin(connectivityVec, ierr); CHKERRQ(ierr)
call VecAssemblyEnd (connectivityVec, ierr); CHKERRQ(ierr)
call VecView(connectivityVec, resUnit, ierr); CHKERRQ(ierr)
call VecDestroy(connectivityVec, ierr); CHKERRQ(ierr)
endif
enddo; enddo
do phase = 1, material_Nphase; do grain = 1, homogenization_maxNgrains
mappingCells = 0_pInt
mappingCells(worldrank+1) = phaseOutput(phase,grain)%sizeIpCells
call MPI_Allreduce(MPI_IN_PLACE,mappingCells,worldsize,MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr)
call VecCreateMPI(PETSC_COMM_WORLD,mappingCells(worldrank+1),sum(mappingCells), &
phaseResultsVec(phase,grain),ierr);CHKERRQ(ierr)
if (sum(mappingCells) > 0) then
call VecCreateMPI(PETSC_COMM_WORLD,mappingCells(worldrank+1)*2**dimPlex,sum(mappingCells)*2**dimPlex, &
connectivityVec,ierr);CHKERRQ(ierr)
write(grainStr,'(a,i0)') 'Grain',grain
call PetscObjectSetName(connectivityVec,&
'mapping_'//trim(phase_name(phase))//'_'// &
trim(grainStr),ierr)
CHKERRQ(ierr)
call VecGetArrayF90(connectivityVec, results, ierr)
CHKERRQ(ierr)
results = 0.0_pReal; ctr = 1_pInt
do cell = cellStart, cellEnd-1; do qPt = 1, mesh_maxNips
if (material_phase(grain,qPt,cell+1) == phase) then
results(ctr:ctr+2**dimPlex-1) = real(reshape(connectivity(1:2**dimPlex,mesh_maxNips*cell+qPt), &
shape=[2**dimPlex]))
ctr = ctr + 2**dimPlex
endif
enddo; enddo
call VecRestoreArrayF90(connectivityVec, results, ierr)
CHKERRQ(ierr)
call VecAssemblyBegin(connectivityVec, ierr);CHKERRQ(ierr)
call VecAssemblyEnd (connectivityVec, ierr);CHKERRQ(ierr)
call VecView(connectivityVec, resUnit, ierr);CHKERRQ(ierr)
call VecDestroy(connectivityVec, ierr); CHKERRQ(ierr)
endif
enddo; enddo
if (worldrank == 0_pInt) then
do homog = 1, material_Nhomogenization
call VecGetSize(homogenizationResultsVec(homog),mappingCells(1),ierr)
CHKERRQ(ierr)
if (mappingCells(1) > 0) &
write(headerID, '(a,i0)') 'number of homog_'// &
trim(homogenization_name(homog))//'_'// &
'cells : ', mappingCells(1)
enddo
do cryst = 1, material_Ncrystallite; do grain = 1, homogenization_maxNgrains
call VecGetSize(crystalliteResultsVec(cryst,grain),mappingCells(1),ierr)
CHKERRQ(ierr)
write(grainStr,'(a,i0)') 'Grain',grain
if (mappingCells(1) > 0) &
write(headerID, '(a,i0)') 'number of cryst_'// &
trim(crystallite_name(cryst))//'_'// &
trim(grainStr)//'_'// &
'cells : ', mappingCells(1)
enddo; enddo
do phase = 1, material_Nphase; do grain = 1, homogenization_maxNgrains
call VecGetSize(phaseResultsVec(phase,grain),mappingCells(1),ierr)
CHKERRQ(ierr)
write(grainStr,'(a,i0)') 'Grain',grain
if (mappingCells(1) > 0) &
write(headerID, '(a,i0)') 'number of phase_'// &
trim(phase_name(phase))//'_'//trim(grainStr)//'_'// &
'cells : ', mappingCells(1)
enddo; enddo
close(headerID)
endif
!allocate(mappingCells(worldsize), source = 0)
!do homog = 1, material_Nhomogenization
! mappingCells = 0_pInt; mappingCells(worldrank+1) = homogOutput(homog)%sizeIpCells
! call MPI_Allreduce(MPI_IN_PLACE,mappingCells,worldsize,MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr)
! call VecCreateMPI(PETSC_COMM_WORLD,mappingCells(worldrank+1),sum(mappingCells), &
! homogenizationResultsVec(homog),ierr);CHKERRQ(ierr)
! if (sum(mappingCells) > 0) then
! call VecCreateMPI(PETSC_COMM_WORLD,mappingCells(worldrank+1)*2**dimPlex,sum(mappingCells)*2**dimPlex, &
! connectivityVec,ierr);CHKERRQ(ierr)
! call PetscObjectSetName(connectivityVec,'mapping_'//trim(homogenization_name(homog)),ierr)
! CHKERRQ(ierr)
! call VecGetArrayF90(connectivityVec,results,ierr); CHKERRQ(ierr)
! results = 0.0_pReal; ctr = 1_pInt
! do cell = cellStart, cellEnd-1; do qPt = 1, mesh_maxNips
! if (material_homog(qPt,cell+1) == homog) then
! results(ctr:ctr+2**dimPlex-1) = real(reshape(connectivity(1:2**dimPlex,mesh_maxNips*cell+qPt), &
! shape=[2**dimPlex]))
! ctr = ctr + 2**dimPlex
! endif
! enddo; enddo
! call VecRestoreArrayF90(connectivityVec, results, ierr); CHKERRQ(ierr)
! call VecAssemblyBegin(connectivityVec, ierr); CHKERRQ(ierr)
! call VecAssemblyEnd (connectivityVec, ierr); CHKERRQ(ierr)
! call VecView(connectivityVec, resUnit, ierr); CHKERRQ(ierr)
! call VecDestroy(connectivityVec, ierr); CHKERRQ(ierr)
! endif
!enddo
!do cryst = 1, material_Ncrystallite; do grain = 1, homogenization_maxNgrains
! mappingCells = 0_pInt
! mappingCells(worldrank+1) = crystalliteOutput(cryst,grain)%sizeIpCells
! call MPI_Allreduce(MPI_IN_PLACE,mappingCells,worldsize,MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr)
! call VecCreateMPI(PETSC_COMM_WORLD,mappingCells(worldrank+1),sum(mappingCells), &
! crystalliteResultsVec(cryst,grain),ierr);CHKERRQ(ierr)
! if (sum(mappingCells) > 0) then
! call VecCreateMPI(PETSC_COMM_WORLD,mappingCells(worldrank+1)*2**dimPlex,sum(mappingCells)*2**dimPlex, &
! connectivityVec,ierr);CHKERRQ(ierr)
! write(grainStr,'(a,i0)') 'Grain',grain
! call PetscObjectSetName(connectivityVec,'mapping_'// &
! trim(crystallite_name(cryst))//'_'// &
! trim(grainStr),ierr)
! CHKERRQ(ierr)
! call VecGetArrayF90(connectivityVec, results, ierr); CHKERRQ(ierr)
! results = 0.0_pReal; ctr = 1_pInt
! do cell = cellStart, cellEnd-1; do qPt = 1, mesh_maxNips
! if (homogenization_Ngrains (mesh_element(3,cell+1)) >= grain .and. &
! microstructure_crystallite(mesh_element(4,cell+1)) == cryst) then
! results(ctr:ctr+2**dimPlex-1) = real(reshape(connectivity(1:2**dimPlex,mesh_maxNips*cell+qPt), &
! shape=[2**dimPlex]))
! ctr = ctr + 2**dimPlex
! endif
! enddo; enddo
! call VecRestoreArrayF90(connectivityVec, results, ierr); CHKERRQ(ierr)
! call VecAssemblyBegin(connectivityVec, ierr); CHKERRQ(ierr)
! call VecAssemblyEnd (connectivityVec, ierr); CHKERRQ(ierr)
! call VecView(connectivityVec, resUnit, ierr); CHKERRQ(ierr)
! call VecDestroy(connectivityVec, ierr); CHKERRQ(ierr)
! endif
!enddo; enddo
!do phase = 1, material_Nphase; do grain = 1, homogenization_maxNgrains
! mappingCells = 0_pInt
! mappingCells(worldrank+1) = phaseOutput(phase,grain)%sizeIpCells
! call MPI_Allreduce(MPI_IN_PLACE,mappingCells,worldsize,MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr)
! call VecCreateMPI(PETSC_COMM_WORLD,mappingCells(worldrank+1),sum(mappingCells), &
! phaseResultsVec(phase,grain),ierr);CHKERRQ(ierr)
! if (sum(mappingCells) > 0) then
! call VecCreateMPI(PETSC_COMM_WORLD,mappingCells(worldrank+1)*2**dimPlex,sum(mappingCells)*2**dimPlex, &
! connectivityVec,ierr);CHKERRQ(ierr)
! write(grainStr,'(a,i0)') 'Grain',grain
! call PetscObjectSetName(connectivityVec,&
! 'mapping_'//trim(phase_name(phase))//'_'// &
! trim(grainStr),ierr)
! CHKERRQ(ierr)
! call VecGetArrayF90(connectivityVec, results, ierr)
! CHKERRQ(ierr)
! results = 0.0_pReal; ctr = 1_pInt
! do cell = cellStart, cellEnd-1; do qPt = 1, mesh_maxNips
! if (material_phase(grain,qPt,cell+1) == phase) then
! results(ctr:ctr+2**dimPlex-1) = real(reshape(connectivity(1:2**dimPlex,mesh_maxNips*cell+qPt), &
! shape=[2**dimPlex]))
! ctr = ctr + 2**dimPlex
! endif
! enddo; enddo
! call VecRestoreArrayF90(connectivityVec, results, ierr)
! CHKERRQ(ierr)
! call VecAssemblyBegin(connectivityVec, ierr);CHKERRQ(ierr)
! call VecAssemblyEnd (connectivityVec, ierr);CHKERRQ(ierr)
! call VecView(connectivityVec, resUnit, ierr);CHKERRQ(ierr)
! call VecDestroy(connectivityVec, ierr); CHKERRQ(ierr)
! endif
!enddo; enddo
!if (worldrank == 0_pInt) then
! do homog = 1, material_Nhomogenization
! call VecGetSize(homogenizationResultsVec(homog),mappingCells(1),ierr)
! CHKERRQ(ierr)
! if (mappingCells(1) > 0) &
! write(headerID, '(a,i0)') 'number of homog_'// &
! trim(homogenization_name(homog))//'_'// &
! 'cells : ', mappingCells(1)
! enddo
! do cryst = 1, material_Ncrystallite; do grain = 1, homogenization_maxNgrains
! call VecGetSize(crystalliteResultsVec(cryst,grain),mappingCells(1),ierr)
! CHKERRQ(ierr)
! write(grainStr,'(a,i0)') 'Grain',grain
! if (mappingCells(1) > 0) &
! write(headerID, '(a,i0)') 'number of cryst_'// &
! trim(crystallite_name(cryst))//'_'// &
! trim(grainStr)//'_'// &
! 'cells : ', mappingCells(1)
! enddo; enddo
! do phase = 1, material_Nphase; do grain = 1, homogenization_maxNgrains
! call VecGetSize(phaseResultsVec(phase,grain),mappingCells(1),ierr)
! CHKERRQ(ierr)
! write(grainStr,'(a,i0)') 'Grain',grain
! if (mappingCells(1) > 0) &
! write(headerID, '(a,i0)') 'number of phase_'// &
! trim(phase_name(phase))//'_'//trim(grainStr)//'_'// &
! 'cells : ', mappingCells(1)
! enddo; enddo
! close(headerID)
!endif
end subroutine utilities_init
@ -509,13 +458,12 @@ subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData)
math_det33
use FEsolving, only: &
restartWrite
use CPFEM2, only: &
CPFEM_general
use homogenization, only: &
materialpoint_F0, &
materialpoint_F, &
materialpoint_P, &
materialpoint_dPdF
materialpoint_dPdF, &
materialpoint_stressAndItsTangent
use mesh, only: &
mesh_NcpElems
@ -560,7 +508,7 @@ subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData)
flush(6)
endif
call CPFEM_general(age,timeinc)
call materialpoint_stressAndItsTangent(.true.,timeinc) ! calculate P field
call debug_info()
@ -791,27 +739,24 @@ end subroutine utilities_indexActiveSet
!--------------------------------------------------------------------------------------------------
subroutine utilities_destroy()
use material, only: &
material_Nhomogenization, &
material_Ncrystallite, &
material_Nphase, &
homogenization_Ngrains
implicit none
PetscInt :: homog, cryst, grain, phase
PetscErrorCode :: ierr
!implicit none
!PetscInt :: homog, cryst, grain, phase
!PetscErrorCode :: ierr
call PetscViewerHDF5PopGroup(resUnit, ierr); CHKERRQ(ierr)
call VecDestroy(coordinatesVec,ierr); CHKERRQ(ierr)
do homog = 1, material_Nhomogenization
call VecDestroy(homogenizationResultsVec(homog),ierr);CHKERRQ(ierr)
do cryst = 1, material_Ncrystallite; do grain = 1, homogenization_Ngrains(homog)
call VecDestroy(crystalliteResultsVec(cryst,grain),ierr);CHKERRQ(ierr)
enddo; enddo
do phase = 1, material_Nphase; do grain = 1, homogenization_Ngrains(homog)
call VecDestroy(phaseResultsVec(phase,grain),ierr);CHKERRQ(ierr)
enddo; enddo
enddo
call PetscViewerDestroy(resUnit, ierr); CHKERRQ(ierr)
!call PetscViewerHDF5PopGroup(resUnit, ierr); CHKERRQ(ierr)
!call VecDestroy(coordinatesVec,ierr); CHKERRQ(ierr)
!do homog = 1, material_Nhomogenization
! call VecDestroy(homogenizationResultsVec(homog),ierr);CHKERRQ(ierr)
! do cryst = 1, material_Ncrystallite; do grain = 1, homogenization_Ngrains(homog)
! call VecDestroy(crystalliteResultsVec(cryst,grain),ierr);CHKERRQ(ierr)
! enddo; enddo
! do phase = 1, material_Nphase; do grain = 1, homogenization_Ngrains(homog)
! call VecDestroy(phaseResultsVec(phase,grain),ierr);CHKERRQ(ierr)
! enddo; enddo
!enddo
!call PetscViewerDestroy(resUnit, ierr); CHKERRQ(ierr)
end subroutine utilities_destroy

16
src/FEM_zoo.f90
View File

@ -6,7 +6,6 @@ module FEM_Zoo
use prec, only: pReal, pInt, p_vec
implicit none
#include <petsc/finclude/petsc.h90>
private
integer(pInt), parameter, public:: &
maxOrder = 5 !< current max interpolation set at cubic (intended to be arbitrary)
@ -35,26 +34,23 @@ contains
!> @brief initializes FEM interpolation data
!--------------------------------------------------------------------------------------------------
subroutine FEM_Zoo_init
use, intrinsic :: iso_fortran_env
#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800
use, intrinsic :: iso_fortran_env, only: &
compiler_version, &
compiler_options
#endif
use IO, only: &
IO_timeStamp
use math, only: &
math_binomial
implicit none
PetscInt :: worldrank
PetscErrorCode :: ierr
external :: &
MPI_Comm_rank, &
MPI_abort
call MPI_Comm_rank(PETSC_COMM_WORLD,worldrank,ierr);CHKERRQ(ierr)
if (worldrank == 0) then
write(6,'(/,a)') ' <<<+- FEM_Zoo init -+>>>'
write(6,'(a)') ' $Id: FEM_Zoo.f90 4354 2015-08-04 15:04:53Z MPIE\p.shanthraj $'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
endif
!--------------------------------------------------------------------------------------------------
! 2D linear
FEM_Zoo_nQuadrature(2,1) = 1

112
src/homogenization.f90
View File

@ -6,9 +6,6 @@
!--------------------------------------------------------------------------------------------------
module homogenization
use prec, only: &
#ifdef FEM
tOutputData, &
#endif
pInt, &
pReal
@ -22,16 +19,8 @@ module homogenization
materialpoint_P !< first P--K stress of IP
real(pReal), dimension(:,:,:,:,:,:), allocatable, public :: &
materialpoint_dPdF !< tangent of first P--K stress at IP
#ifdef FEM
type(tOutputData), dimension(:), allocatable, public :: &
homogOutput
type(tOutputData), dimension(:,:), allocatable, public :: &
crystalliteOutput, &
phaseOutput
#else
real(pReal), dimension(:,:,:), allocatable, public :: &
materialpoint_results !< results array of material point
#endif
integer(pInt), public, protected :: &
materialpoint_sizeResults, &
homogenization_maxSizePostResults, &
@ -90,16 +79,11 @@ subroutine homogenization_init
mesh_element, &
FE_Nips, &
FE_geomtype
#ifdef FEM
use crystallite, only: &
crystallite_sizePostResults
#else
use constitutive, only: &
constitutive_plasticity_maxSizePostResults, &
constitutive_source_maxSizePostResults
use crystallite, only: &
crystallite_maxSizePostResults
#endif
use config, only: &
config_deallocate, &
material_configFile, &
@ -411,33 +395,6 @@ subroutine homogenization_init
hydrogenflux_maxSizePostResults = max(hydrogenflux_maxSizePostResults ,hydrogenfluxState(p)%sizePostResults)
enddo
#ifdef FEM
allocate(homogOutput (material_Nhomogenization ))
allocate(crystalliteOutput(material_Ncrystallite, homogenization_maxNgrains))
allocate(phaseOutput (material_Nphase, homogenization_maxNgrains))
do p = 1, material_Nhomogenization
homogOutput(p)%sizeResults = homogState (p)%sizePostResults + &
thermalState (p)%sizePostResults + &
damageState (p)%sizePostResults + &
vacancyfluxState (p)%sizePostResults + &
porosityState (p)%sizePostResults + &
hydrogenfluxState(p)%sizePostResults
homogOutput(p)%sizeIpCells = count(material_homog==p)
allocate(homogOutput(p)%output(homogOutput(p)%sizeResults,homogOutput(p)%sizeIpCells))
enddo
do p = 1, material_Ncrystallite; do e = 1, homogenization_maxNgrains
crystalliteOutput(p,e)%sizeResults = crystallite_sizePostResults(p)
crystalliteOutput(p,e)%sizeIpCells = count(microstructure_crystallite(mesh_element(4,:)) == p .and. &
homogenization_Ngrains (mesh_element(3,:)) >= e)*mesh_maxNips
allocate(crystalliteOutput(p,e)%output(crystalliteOutput(p,e)%sizeResults,crystalliteOutput(p,e)%sizeIpCells))
enddo; enddo
do p = 1, material_Nphase; do e = 1, homogenization_maxNgrains
phaseOutput(p,e)%sizeResults = plasticState (p)%sizePostResults + &
sum(sourceState (p)%p(:)%sizePostResults)
phaseOutput(p,e)%sizeIpCells = count(material_phase(e,:,:) == p)
allocate(phaseOutput(p,e)%output(phaseOutput(p,e)%sizeResults,phaseOutput(p,e)%sizeIpCells))
enddo; enddo
#else
materialpoint_sizeResults = 1 & ! grain count
+ 1 + homogenization_maxSizePostResults & ! homogSize & homogResult
+ thermal_maxSizePostResults &
@ -449,7 +406,6 @@ subroutine homogenization_init
+ 1 + constitutive_plasticity_maxSizePostResults & ! constitutive size & constitutive results
+ constitutive_source_maxSizePostResults)
allocate(materialpoint_results(materialpoint_sizeResults,mesh_maxNips,mesh_NcpElems))
#endif
write(6,'(/,a)') ' <<<+- homogenization init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
@ -473,9 +429,6 @@ subroutine homogenization_init
write(6,'(a32,1x,7(i8,1x))') 'materialpoint_requested: ', shape(materialpoint_requested)
write(6,'(a32,1x,7(i8,1x))') 'materialpoint_converged: ', shape(materialpoint_converged)
write(6,'(a32,1x,7(i8,1x),/)') 'materialpoint_doneAndHappy: ', shape(materialpoint_doneAndHappy)
#ifndef FEM
write(6,'(a32,1x,7(i8,1x),/)') 'materialpoint_results: ', shape(materialpoint_results)
#endif
write(6,'(a32,1x,7(i8,1x))') 'maxSizePostResults: ', homogenization_maxSizePostResults
endif
flush(6)
@ -904,33 +857,18 @@ subroutine materialpoint_postResults
mesh_element
use material, only: &
mappingHomogenization, &
#ifdef FEM
phaseAt, phasememberAt, &
homogenization_maxNgrains, &
material_Ncrystallite, &
material_Nphase, &
#else
homogState, &
thermalState, &
damageState, &
vacancyfluxState, &
porosityState, &
hydrogenfluxState, &
#endif
plasticState, &
sourceState, &
material_phase, &
homogenization_Ngrains, &
microstructure_crystallite
#ifdef FEM
use constitutive, only: &
constitutive_plasticity_maxSizePostResults, &
constitutive_source_maxSizePostResults
#endif
use crystallite, only: &
#ifdef FEM
crystallite_maxSizePostResults, &
#endif
crystallite_sizePostResults, &
crystallite_postResults
@ -943,55 +881,6 @@ subroutine materialpoint_postResults
g, & !< grain number
i, & !< integration point number
e !< element number
#ifdef FEM
integer(pInt) :: &
myHomog, &
myPhase, &
crystalliteCtr(material_Ncrystallite, homogenization_maxNgrains), &
phaseCtr (material_Nphase, homogenization_maxNgrains)
real(pReal), dimension(1+crystallite_maxSizePostResults + &
1+constitutive_plasticity_maxSizePostResults + &
constitutive_source_maxSizePostResults) :: &
crystalliteResults
crystalliteCtr = 0_pInt; phaseCtr = 0_pInt
elementLooping: do e = FEsolving_execElem(1),FEsolving_execElem(2)
myNgrains = homogenization_Ngrains(mesh_element(3,e))
myCrystallite = microstructure_crystallite(mesh_element(4,e))
IpLooping: do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
myHomog = mappingHomogenization(2,i,e)
thePos = mappingHomogenization(1,i,e)
homogOutput(myHomog)%output(1: &
homogOutput(myHomog)%sizeResults, &
thePos) = homogenization_postResults(i,e)
grainLooping :do g = 1,myNgrains
myPhase = phaseAt(g,i,e)
crystalliteResults(1:1+crystallite_sizePostResults(myCrystallite) + &
1+plasticState(myPhase)%sizePostResults + &
sum(sourceState(myPhase)%p(:)%sizePostResults)) = crystallite_postResults(g,i,e)
if (microstructure_crystallite(mesh_element(4,e)) == myCrystallite .and. &
homogenization_Ngrains (mesh_element(3,e)) >= g) then
crystalliteCtr(myCrystallite,g) = crystalliteCtr(myCrystallite,g) + 1_pInt
crystalliteOutput(myCrystallite,g)% &
output(1:crystalliteOutput(myCrystallite,g)%sizeResults,crystalliteCtr(myCrystallite,g)) = &
crystalliteResults(2:1+crystalliteOutput(myCrystallite,g)%sizeResults)
endif
if (material_phase(g,i,e) == myPhase) then
phaseCtr(myPhase,g) = phaseCtr(myPhase,g) + 1_pInt
phaseOutput(myPhase,g)% &
output(1:phaseOutput(myPhase,g)%sizeResults,phaseCtr(myPhase,g)) = &
crystalliteResults(3 + crystalliteOutput(myCrystallite,g)%sizeResults: &
1 + crystalliteOutput(myCrystallite,g)%sizeResults + &
1 + plasticState (myphase)%sizePostResults + &
sum(sourceState(myphase)%p(:)%sizePostResults))
endif
enddo grainLooping
enddo IpLooping
enddo elementLooping
#else
!$OMP PARALLEL DO PRIVATE(myNgrains,myCrystallite,thePos,theSize)
elementLooping: do e = FEsolving_execElem(1),FEsolving_execElem(2)
@ -1027,7 +916,6 @@ subroutine materialpoint_postResults
enddo IpLooping
enddo elementLooping
!$OMP END PARALLEL DO
#endif
end subroutine materialpoint_postResults

444
src/meshFEM.f90 Normal file
View File

@ -0,0 +1,444 @@
!--------------------------------------------------------------------------------------------------
!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Driver controlling inner and outer load case looping of the FEM solver
!> @details doing cutbacking, forwarding in case of restart, reporting statistics, writing
!> results
!--------------------------------------------------------------------------------------------------
module mesh
#include <petsc/finclude/petscis.h>
#include <petsc/finclude/petscdmda.h>
use prec, only: pReal, pInt
use PETScdmda
use PETScis
implicit none
private
integer(pInt), public, protected :: &
mesh_Nboundaries, &
mesh_NcpElems, & !< total number of CP elements in mesh
mesh_NcpElemsGlobal, &
mesh_Nnodes, & !< total number of nodes in mesh
mesh_maxNnodes, & !< max number of nodes in any CP element
mesh_maxNips, & !< max number of IPs in any CP element
mesh_maxNipNeighbors, &
mesh_Nelems !< total number of elements in mesh
real(pReal), public, protected :: charLength
integer(pInt), dimension(:,:), allocatable, public, protected :: &
mesh_element !< FEid, type(internal representation), material, texture, node indices as CP IDs
real(pReal), dimension(:,:), allocatable, public :: &
mesh_node !< node x,y,z coordinates (after deformation! ONLY FOR MARC!!!)
real(pReal), dimension(:,:), allocatable, public, protected :: &
mesh_ipVolume, & !< volume associated with IP (initially!)
mesh_node0 !< node x,y,z coordinates (initially!)
real(pReal), dimension(:,:,:), allocatable, public :: &
mesh_ipCoordinates !< IP x,y,z coordinates (after deformation!)
real(pReal), dimension(:,:,:), allocatable, public, protected :: &
mesh_ipArea !< area of interface to neighboring IP (initially!)
real(pReal),dimension(:,:,:,:), allocatable, public, protected :: &
mesh_ipAreaNormal !< area normal of interface to neighboring IP (initially!)
integer(pInt), dimension(:,:,:,:), allocatable, public, protected :: &
mesh_ipNeighborhood !< 6 or less neighboring IPs as [element_num, IP_index, neighbor_index that points to me]
logical, dimension(3), public, protected :: mesh_periodicSurface !< flag indicating periodic outer surfaces (used for fluxes)
integer(pInt), dimension(:,:), allocatable, target, private :: &
mesh_mapFEtoCPelem, & !< [sorted FEid, corresponding CPid]
mesh_mapFEtoCPnode !< [sorted FEid, corresponding CPid]
DM, public :: geomMesh
integer(pInt), dimension(:), allocatable, public, protected :: &
mesh_boundaries
! These definitions should actually reside in the FE-solver specific part (different for MARC/ABAQUS)
! Hence, I suggest to prefix with "FE_"
integer(pInt), parameter, public :: &
FE_Nelemtypes = 1_pInt, &
FE_Ngeomtypes = 1_pInt, &
FE_Ncelltypes = 1_pInt, &
FE_maxNnodes = 1_pInt, &
FE_maxNips = 14_pInt
integer(pInt), dimension(FE_Nelemtypes), parameter, public :: FE_geomtype = & !< geometry type of particular element type
int([1],pInt)
integer(pInt), dimension(FE_Ngeomtypes), parameter, public :: FE_celltype = & !< cell type that is used by each geometry type
int([1],pInt)
integer(pInt), dimension(FE_Nelemtypes), parameter, public :: FE_Nnodes = & !< number of nodes that constitute a specific type of element
int([0],pInt)
integer(pInt), dimension(FE_Ngeomtypes), public :: FE_Nips = & !< number of IPs in a specific type of element
int([0],pInt)
integer(pInt), dimension(FE_Ncelltypes), parameter, public :: FE_NipNeighbors = & !< number of ip neighbors / cell faces in a specific cell type
int([6],pInt)
public :: &
mesh_init, &
mesh_FEasCP, &
mesh_FEM_build_ipVolumes, &
mesh_FEM_build_ipCoordinates, &
mesh_cellCenterCoordinates
external :: &
MPI_Bcast, &
DMPlexCreateFromFile, &
DMPlexDistribute, &
DMPlexCopyCoordinates, &
DMGetStratumSize, &
DMPlexGetHeightStratum, &
DMPlexGetLabelValue, &
DMPlexSetLabelValue
contains
!--------------------------------------------------------------------------------------------------
!> @brief initializes the mesh by calling all necessary private routines the mesh module
!! Order and routines strongly depend on type of solver
!--------------------------------------------------------------------------------------------------
subroutine mesh_init(ip,el)
use DAMASK_interface
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use IO, only: &
IO_timeStamp, &
IO_error, &
IO_open_file, &
IO_stringPos, &
IO_intValue, &
IO_EOF, &
IO_read, &
IO_isBlank
use debug, only: &
debug_e, &
debug_i
use numerics, only: &
usePingPong, &
integrationOrder, &
worldrank, &
worldsize
use FEsolving, only: &
FEsolving_execElem, &
FEsolving_execIP, &
calcMode
use FEM_Zoo, only: &
FEM_Zoo_nQuadrature, &
FEM_Zoo_QuadraturePoints
implicit none
integer(pInt), parameter :: FILEUNIT = 222_pInt
integer(pInt), intent(in) :: el, ip
integer(pInt) :: j
integer(pInt), allocatable, dimension(:) :: chunkPos
integer :: dimPlex
character(len=512) :: &
line
logical :: flag
PetscSF :: sf
DM :: globalMesh
PetscInt :: face, nFaceSets
PetscInt, pointer :: pFaceSets(:)
IS :: faceSetIS
PetscErrorCode :: ierr
write(6,'(/,a)') ' <<<+- mesh init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
if (allocated(mesh_mapFEtoCPelem)) deallocate(mesh_mapFEtoCPelem)
if (allocated(mesh_mapFEtoCPnode)) deallocate(mesh_mapFEtoCPnode)
if (allocated(mesh_node0)) deallocate(mesh_node0)
if (allocated(mesh_node)) deallocate(mesh_node)
if (allocated(mesh_element)) deallocate(mesh_element)
if (allocated(mesh_ipCoordinates)) deallocate(mesh_ipCoordinates)
if (allocated(mesh_ipVolume)) deallocate(mesh_ipVolume)
call DMPlexCreateFromFile(PETSC_COMM_WORLD,geometryFile,PETSC_TRUE,globalMesh,ierr)
CHKERRQ(ierr)
call DMGetDimension(globalMesh,dimPlex,ierr)
CHKERRQ(ierr)
call DMGetStratumSize(globalMesh,'depth',dimPlex,mesh_NcpElemsGlobal,ierr)
CHKERRQ(ierr)
call DMGetLabelSize(globalMesh,'Face Sets',mesh_Nboundaries,ierr)
CHKERRQ(ierr)
call MPI_Bcast(mesh_Nboundaries,1,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr)
call MPI_Bcast(mesh_NcpElemsGlobal,1,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr)
call MPI_Bcast(dimPlex,1,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr)
allocate(mesh_boundaries(mesh_Nboundaries), source = 0_pInt)
call DMGetLabelSize(globalMesh,'Face Sets',nFaceSets,ierr)
CHKERRQ(ierr)
call DMGetLabelIdIS(globalMesh,'Face Sets',faceSetIS,ierr)
CHKERRQ(ierr)
if (nFaceSets > 0) call ISGetIndicesF90(faceSetIS,pFaceSets,ierr)
do face = 1, nFaceSets
mesh_boundaries(face) = pFaceSets(face)
enddo
if (nFaceSets > 0) call ISRestoreIndicesF90(faceSetIS,pFaceSets,ierr)
call MPI_Bcast(mesh_boundaries,mesh_Nboundaries,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr)
if (worldrank == 0) then
j = 0
flag = .false.
call IO_open_file(FILEUNIT,trim(geometryFile))
do
read(FILEUNIT,'(a512)') line
if (trim(line) == IO_EOF) exit ! skip empty lines
if (trim(line) == '$Elements') then
read(FILEUNIT,'(a512)') line
read(FILEUNIT,'(a512)') line
flag = .true.
endif
if (trim(line) == '$EndElements') exit
if (flag) then
chunkPos = IO_stringPos(line)
if (chunkPos(1) == 3+IO_intValue(line,chunkPos,3)+dimPlex+1) then
call DMSetLabelValue(globalMesh,'material',j,IO_intValue(line,chunkPos,4),ierr)
CHKERRQ(ierr)
j = j + 1
endif ! count all identifiers to allocate memory and do sanity check
endif
enddo
close (FILEUNIT)
endif
if (worldsize > 1) then
call DMPlexDistribute(globalMesh,0,sf,geomMesh,ierr)
CHKERRQ(ierr)
else
call DMClone(globalMesh,geomMesh,ierr)
CHKERRQ(ierr)
endif
call DMDestroy(globalMesh,ierr); CHKERRQ(ierr)
call DMGetStratumSize(geomMesh,'depth',dimPlex,mesh_Nelems,ierr)
CHKERRQ(ierr)
call DMGetStratumSize(geomMesh,'depth',0,mesh_Nnodes,ierr)
CHKERRQ(ierr)
mesh_NcpElems = mesh_Nelems
call mesh_FEM_mapNodesAndElems
FE_Nips(FE_geomtype(1_pInt)) = FEM_Zoo_nQuadrature(dimPlex,integrationOrder)
mesh_maxNnodes = FE_Nnodes(1_pInt)
mesh_maxNips = FE_Nips(1_pInt)
call mesh_FEM_build_ipCoordinates(dimPlex,FEM_Zoo_QuadraturePoints(dimPlex,integrationOrder)%p)
call mesh_FEM_build_ipVolumes(dimPlex)
allocate (mesh_element (4_pInt+mesh_maxNnodes,mesh_NcpElems)); mesh_element = 0_pInt
do j = 1, mesh_NcpElems
mesh_element( 1,j) = j
mesh_element( 2,j) = 1_pInt ! elem type
mesh_element( 3,j) = 1_pInt ! homogenization
call DMGetLabelValue(geomMesh,'material',j-1,mesh_element(4,j),ierr)
CHKERRQ(ierr)
end do
if (usePingPong .and. (mesh_Nelems /= mesh_NcpElems)) &
call IO_error(600_pInt) ! ping-pong must be disabled when having non-DAMASK elements
if (debug_e < 1 .or. debug_e > mesh_NcpElems) &
call IO_error(602_pInt,ext_msg='element') ! selected element does not exist
if (debug_i < 1 .or. debug_i > FE_Nips(FE_geomtype(mesh_element(2_pInt,debug_e)))) &
call IO_error(602_pInt,ext_msg='IP') ! selected element does not have requested IP
FEsolving_execElem = [ 1_pInt,mesh_NcpElems ] ! parallel loop bounds set to comprise all DAMASK elements
if (allocated(FEsolving_execIP)) deallocate(FEsolving_execIP)
allocate(FEsolving_execIP(2_pInt,mesh_NcpElems)); FEsolving_execIP = 1_pInt ! parallel loop bounds set to comprise from first IP...
forall (j = 1_pInt:mesh_NcpElems) FEsolving_execIP(2,j) = FE_Nips(FE_geomtype(mesh_element(2,j))) ! ...up to own IP count for each element
if (allocated(calcMode)) deallocate(calcMode)
allocate(calcMode(mesh_maxNips,mesh_NcpElems))
calcMode = .false. ! pretend to have collected what first call is asking (F = I)
calcMode(ip,mesh_FEasCP('elem',el)) = .true. ! first ip,el needs to be already pingponged to "calc"
end subroutine mesh_init
!--------------------------------------------------------------------------------------------------
!> @brief Gives the FE to CP ID mapping by binary search through lookup array
!! valid questions (what) are 'elem', 'node'
!--------------------------------------------------------------------------------------------------
integer(pInt) function mesh_FEasCP(what,myID)
use IO, only: &
IO_lc
implicit none
character(len=*), intent(in) :: what
integer(pInt), intent(in) :: myID
integer(pInt), dimension(:,:), pointer :: lookupMap
integer(pInt) :: lower,upper,center
mesh_FEasCP = 0_pInt
select case(IO_lc(what(1:4)))
case('elem')
lookupMap => mesh_mapFEtoCPelem
case('node')
lookupMap => mesh_mapFEtoCPnode
case default
return
endselect
lower = 1_pInt
upper = int(size(lookupMap,2_pInt),pInt)
if (lookupMap(1_pInt,lower) == myID) then ! check at bounds QUESTION is it valid to extend bounds by 1 and just do binary search w/o init check at bounds?
mesh_FEasCP = lookupMap(2_pInt,lower)
return
elseif (lookupMap(1_pInt,upper) == myID) then
mesh_FEasCP = lookupMap(2_pInt,upper)
return
endif
binarySearch: do while (upper-lower > 1_pInt)
center = (lower+upper)/2_pInt
if (lookupMap(1_pInt,center) < myID) then
lower = center
elseif (lookupMap(1_pInt,center) > myID) then
upper = center
else
mesh_FEasCP = lookupMap(2_pInt,center)
exit
endif
enddo binarySearch
end function mesh_FEasCP
!--------------------------------------------------------------------------------------------------
!> @brief Calculates cell center coordinates.
!--------------------------------------------------------------------------------------------------
pure function mesh_cellCenterCoordinates(ip,el)
implicit none
integer(pInt), intent(in) :: el, & !< element number
ip !< integration point number
real(pReal), dimension(3) :: mesh_cellCenterCoordinates !< x,y,z coordinates of the cell center of the requested IP cell
end function mesh_cellCenterCoordinates
!--------------------------------------------------------------------------------------------------
!> @brief Calculates IP volume. Allocates global array 'mesh_ipVolume'
!> @details The IP volume is calculated differently depending on the cell type.
!> 2D cells assume an element depth of one in order to calculate the volume.
!> For the hexahedral cell we subdivide the cell into subvolumes of pyramidal
!> shape with a cell face as basis and the central ip at the tip. This subvolume is
!> calculated as an average of four tetrahedals with three corners on the cell face
!> and one corner at the central ip.
!--------------------------------------------------------------------------------------------------
subroutine mesh_FEM_build_ipVolumes(dimPlex)
use math, only: &
math_I3, &
math_det33
implicit none
PetscInt :: dimPlex
PetscReal :: vol
PetscReal, target :: cent(dimPlex), norm(dimPlex)
PetscReal, pointer :: pCent(:), pNorm(:)
PetscInt :: cellStart, cellEnd, cell
PetscErrorCode :: ierr
if (.not. allocated(mesh_ipVolume)) then
allocate(mesh_ipVolume(mesh_maxNips,mesh_NcpElems))
mesh_ipVolume = 0.0_pReal
endif
call DMPlexGetHeightStratum(geomMesh,0,cellStart,cellEnd,ierr); CHKERRQ(ierr)
pCent => cent
pNorm => norm
do cell = cellStart, cellEnd-1
call DMPlexComputeCellGeometryFVM(geomMesh,cell,vol,pCent,pNorm,ierr)
CHKERRQ(ierr)
mesh_ipVolume(:,cell+1) = vol/real(mesh_maxNips,pReal)
enddo
end subroutine mesh_FEM_build_ipVolumes
!--------------------------------------------------------------------------------------------------
!> @brief Calculates IP Coordinates. Allocates global array 'mesh_ipCoordinates'
! Called by all solvers in mesh_init in order to initialize the ip coordinates.
! Later on the current ip coordinates are directly prvided by the spectral solver and by Abaqus,
! so no need to use this subroutine anymore; Marc however only provides nodal displacements,
! so in this case the ip coordinates are always calculated on the basis of this subroutine.
! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! FOR THE MOMENT THIS SUBROUTINE ACTUALLY CALCULATES THE CELL CENTER AND NOT THE IP COORDINATES,
! AS THE IP IS NOT (ALWAYS) LOCATED IN THE CENTER OF THE IP VOLUME.
! HAS TO BE CHANGED IN A LATER VERSION.
! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!--------------------------------------------------------------------------------------------------
subroutine mesh_FEM_build_ipCoordinates(dimPlex,qPoints)
implicit none
PetscInt, intent(in) :: dimPlex
PetscReal, intent(in) :: qPoints(mesh_maxNips*dimPlex)
PetscReal, target :: v0(dimPlex), cellJ(dimPlex*dimPlex), invcellJ(dimPlex*dimPlex)
PetscReal, pointer :: pV0(:), pCellJ(:), pInvcellJ(:)
PetscReal :: detJ
PetscInt :: cellStart, cellEnd, cell, qPt, dirI, dirJ, qOffset
PetscErrorCode :: ierr
if (.not. allocated(mesh_ipCoordinates)) then
allocate(mesh_ipCoordinates(3,mesh_maxNips,mesh_NcpElems))
mesh_ipCoordinates = 0.0_pReal
endif
pV0 => v0
pCellJ => cellJ
pInvcellJ => invcellJ
call DMPlexGetHeightStratum(geomMesh,0,cellStart,cellEnd,ierr); CHKERRQ(ierr)
do cell = cellStart, cellEnd-1 !< loop over all elements
call DMPlexComputeCellGeometryAffineFEM(geomMesh,cell,pV0,pCellJ,pInvcellJ,detJ,ierr)
CHKERRQ(ierr)
qOffset = 0
do qPt = 1, mesh_maxNips
do dirI = 1, dimPlex
mesh_ipCoordinates(dirI,qPt,cell+1) = pV0(dirI)
do dirJ = 1, dimPlex
mesh_ipCoordinates(dirI,qPt,cell+1) = mesh_ipCoordinates(dirI,qPt,cell+1) + &
pCellJ((dirI-1)*dimPlex+dirJ)*(qPoints(qOffset+dirJ) + 1.0)
enddo
enddo
qOffset = qOffset + dimPlex
enddo
enddo
end subroutine mesh_FEM_build_ipCoordinates
!--------------------------------------------------------------------------------------------------
!> @brief fake map node from FE ID to internal (consecutive) representation for node and element
!! Allocates global array 'mesh_mapFEtoCPnode' and 'mesh_mapFEtoCPelem'
!--------------------------------------------------------------------------------------------------
subroutine mesh_FEM_mapNodesAndElems
use math, only: &
math_range
implicit none
allocate (mesh_mapFEtoCPnode(2_pInt,mesh_Nnodes), source = 0_pInt)
allocate (mesh_mapFEtoCPelem(2_pInt,mesh_NcpElems), source = 0_pInt)
mesh_mapFEtoCPnode = spread(math_range(mesh_Nnodes),1,2)
mesh_mapFEtoCPelem = spread(math_range(mesh_NcpElems),1,2)
end subroutine mesh_FEM_mapNodesAndElems
end module mesh