From d4bcfae82b575e54b457bccf423a5caf02983ede Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 17 Aug 2018 11:23:24 +0200 Subject: [PATCH] WIP: adopting to PETSc 3.9.x and modifications in development branch --- src/CMakeLists.txt | 9 +- src/CPFEM2.f90 | 6 +- src/FEM_interface.f90 | 2 +- src/FEM_mech.f90 | 2 +- src/FEM_utilities.f90 | 363 ++++++++++++++------------------- src/FEM_zoo.f90 | 22 +- src/homogenization.f90 | 112 ----------- src/meshFEM.f90 | 444 +++++++++++++++++++++++++++++++++++++++++ 8 files changed, 614 insertions(+), 346 deletions(-) mode change 100755 => 100644 src/FEM_mech.f90 create mode 100644 src/meshFEM.f90 diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 9418cd56d..caaf0b893 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -57,7 +57,7 @@ if (PROJECT_NAME STREQUAL "DAMASK_spectral") add_dependencies(MESH DAMASK_MATH) list(APPEND OBJECTFILES $) 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 $) 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() diff --git a/src/CPFEM2.f90 b/src/CPFEM2.f90 index c66aa4089..9f75bf8c6 100644 --- a/src/CPFEM2.f90 +++ b/src/CPFEM2.f90 @@ -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 diff --git a/src/FEM_interface.f90 b/src/FEM_interface.f90 index 4a369dd9c..0363ffdaa 100644 --- a/src/FEM_interface.f90 +++ b/src/FEM_interface.f90 @@ -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 diff --git a/src/FEM_mech.f90 b/src/FEM_mech.f90 old mode 100755 new mode 100644 index aa967bec5..6cf47980e --- a/src/FEM_mech.f90 +++ b/src/FEM_mech.f90 @@ -23,7 +23,7 @@ module FEM_mech implicit none private -#include +#include !-------------------------------------------------------------------------------------------------- ! derived types diff --git a/src/FEM_utilities.f90 b/src/FEM_utilities.f90 index 621a32508..e16047da6 100644 --- a/src/FEM_utilities.f90 +++ b/src/FEM_utilities.f90 @@ -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 +#include + use prec, only: pReal, pInt + +use PETScdmda +use PETScis implicit none private -#include +#include !-------------------------------------------------------------------------------------------------- ! 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,8 +508,8 @@ 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() restartWrite = .false. ! reset restartWrite status @@ -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 diff --git a/src/FEM_zoo.f90 b/src/FEM_zoo.f90 index 2c4250098..c34dfb449 100644 --- a/src/FEM_zoo.f90 +++ b/src/FEM_zoo.f90 @@ -6,7 +6,6 @@ module FEM_Zoo use prec, only: pReal, pInt, p_vec implicit none -#include 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() + 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 diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 3565999a8..951527b19 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -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 diff --git a/src/meshFEM.f90 b/src/meshFEM.f90 new file mode 100644 index 000000000..7dc5c93af --- /dev/null +++ b/src/meshFEM.f90 @@ -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 +#include + 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