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

clearer logic, no outdated comments

This commit is contained in:
Martin Diehl 2019-06-15 20:49:48 +02:00
parent e30478127d
commit 43a17a17a2
1 changed files with 60 additions and 70 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

130
src/mesh_FEM.f90
View File

@ -8,11 +8,10 @@ module mesh
#include <petsc/finclude/petscdmplex.h> #include <petsc/finclude/petscdmplex.h>
#include <petsc/finclude/petscis.h> #include <petsc/finclude/petscis.h>
#include <petsc/finclude/petscdmda.h> #include <petsc/finclude/petscdmda.h>
use prec
use mesh_base
use PETScdmplex use PETScdmplex
use PETScdmda use PETScdmda
use PETScis use PETScis
use DAMASK_interface use DAMASK_interface
use IO use IO
use debug use debug
@ -20,6 +19,8 @@ module mesh
use numerics use numerics
use FEsolving use FEsolving
use FEM_Zoo use FEM_Zoo
use prec
use mesh_base
implicit none implicit none
private private
@ -35,13 +36,13 @@ module mesh
mesh_maxNips !< max number of IPs in any CP element mesh_maxNips !< max number of IPs in any CP element
!!!! BEGIN DEPRECATED !!!!! !!!! BEGIN DEPRECATED !!!!!
integer, dimension(:,:), allocatable, public, protected :: & integer, dimension(:,:), allocatable :: &
mesh_element !DEPRECATED mesh_element !DEPRECATED
real(pReal), dimension(:,:), allocatable, public :: & real(pReal), dimension(:,:), allocatable :: &
mesh_node !< node x,y,z coordinates (after deformation! ONLY FOR MARC!!!) mesh_node !< node x,y,z coordinates (after deformation! ONLY FOR MARC!!!)
real(pReal), dimension(:,:), allocatable, public, protected :: & real(pReal), dimension(:,:), allocatable :: &
mesh_ipVolume, & !< volume associated with IP (initially!) mesh_ipVolume, & !< volume associated with IP (initially!)
mesh_node0 !< node x,y,z coordinates (initially!) mesh_node0 !< node x,y,z coordinates (initially!)
@ -176,15 +177,13 @@ subroutine mesh_init
endif endif
enddo enddo
close (FILEUNIT) close (FILEUNIT)
endif
if (worldsize > 1) then
call DMPlexDistribute(globalMesh,0,sf,geomMesh,ierr)
CHKERRQ(ierr)
else
call DMClone(globalMesh,geomMesh,ierr) call DMClone(globalMesh,geomMesh,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
else
call DMPlexDistribute(globalMesh,0,sf,geomMesh,ierr)
CHKERRQ(ierr)
endif endif
call DMDestroy(globalMesh,ierr); CHKERRQ(ierr) call DMDestroy(globalMesh,ierr); CHKERRQ(ierr)
call DMGetStratumSize(geomMesh,'depth',dimPlex,mesh_NcpElems,ierr) call DMGetStratumSize(geomMesh,'depth',dimPlex,mesh_NcpElems,ierr)
@ -255,75 +254,66 @@ end function mesh_cellCenterCoordinates
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine mesh_FEM_build_ipVolumes(dimPlex) subroutine mesh_FEM_build_ipVolumes(dimPlex)
PetscInt :: dimPlex PetscInt :: dimPlex
PetscReal :: vol PetscReal :: vol
PetscReal, target :: cent(dimPlex), norm(dimPlex) PetscReal, target :: cent(dimPlex), norm(dimPlex)
PetscReal, pointer :: pCent(:), pNorm(:) PetscReal, pointer :: pCent(:), pNorm(:)
PetscInt :: cellStart, cellEnd, cell PetscInt :: cellStart, cellEnd, cell
PetscErrorCode :: ierr PetscErrorCode :: ierr
if (.not. allocated(mesh_ipVolume)) then if (.not. allocated(mesh_ipVolume)) then
allocate(mesh_ipVolume(mesh_maxNips,mesh_NcpElems)) allocate(mesh_ipVolume(mesh_maxNips,mesh_NcpElems))
mesh_ipVolume = 0.0_pReal mesh_ipVolume = 0.0_pReal
endif endif
call DMPlexGetHeightStratum(geomMesh,0,cellStart,cellEnd,ierr); CHKERRQ(ierr) call DMPlexGetHeightStratum(geomMesh,0,cellStart,cellEnd,ierr); CHKERRQ(ierr)
pCent => cent pCent => cent
pNorm => norm pNorm => norm
do cell = cellStart, cellEnd-1 do cell = cellStart, cellEnd-1
call DMPlexComputeCellGeometryFVM(geomMesh,cell,vol,pCent,pNorm,ierr) call DMPlexComputeCellGeometryFVM(geomMesh,cell,vol,pCent,pNorm,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
mesh_ipVolume(:,cell+1) = vol/real(mesh_maxNips,pReal) mesh_ipVolume(:,cell+1) = vol/real(mesh_maxNips,pReal)
enddo enddo
end subroutine mesh_FEM_build_ipVolumes end subroutine mesh_FEM_build_ipVolumes
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Calculates IP Coordinates. Allocates global array 'mesh_ipCoordinates' !> @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) subroutine mesh_FEM_build_ipCoordinates(dimPlex,qPoints)
PetscInt, intent(in) :: dimPlex PetscInt, intent(in) :: dimPlex
PetscReal, intent(in) :: qPoints(mesh_maxNips*dimPlex) PetscReal, intent(in) :: qPoints(mesh_maxNips*dimPlex)
PetscReal, target :: v0(dimPlex), cellJ(dimPlex*dimPlex), invcellJ(dimPlex*dimPlex)
PetscReal, pointer :: pV0(:), pCellJ(:), pInvcellJ(:)
PetscReal :: detJ
PetscInt :: cellStart, cellEnd, cell, qPt, dirI, dirJ, qOffset
PetscErrorCode :: ierr
PetscReal, target :: v0(dimPlex), cellJ(dimPlex*dimPlex), invcellJ(dimPlex*dimPlex)
PetscReal, pointer :: pV0(:), pCellJ(:), pInvcellJ(:) allocate(mesh_ipCoordinates(3,mesh_maxNips,mesh_NcpElems),source=0.0_pReal)
PetscReal :: detJ
PetscInt :: cellStart, cellEnd, cell, qPt, dirI, dirJ, qOffset pV0 => v0
PetscErrorCode :: ierr pCellJ => cellJ
pInvcellJ => invcellJ
call DMPlexGetHeightStratum(geomMesh,0,cellStart,cellEnd,ierr); CHKERRQ(ierr)
allocate(mesh_ipCoordinates(3,mesh_maxNips,mesh_NcpElems),source=0.0_pReal) do cell = cellStart, cellEnd-1 !< loop over all elements
call DMPlexComputeCellGeometryAffineFEM(geomMesh,cell,pV0,pCellJ,pInvcellJ,detJ,ierr)
pV0 => v0 CHKERRQ(ierr)
pCellJ => cellJ qOffset = 0
pInvcellJ => invcellJ do qPt = 1, mesh_maxNips
call DMPlexGetHeightStratum(geomMesh,0,cellStart,cellEnd,ierr); CHKERRQ(ierr) do dirI = 1, dimPlex
do cell = cellStart, cellEnd-1 !< loop over all elements mesh_ipCoordinates(dirI,qPt,cell+1) = pV0(dirI)
call DMPlexComputeCellGeometryAffineFEM(geomMesh,cell,pV0,pCellJ,pInvcellJ,detJ,ierr) do dirJ = 1, dimPlex
CHKERRQ(ierr) mesh_ipCoordinates(dirI,qPt,cell+1) = mesh_ipCoordinates(dirI,qPt,cell+1) + &
qOffset = 0 pCellJ((dirI-1)*dimPlex+dirJ)*(qPoints(qOffset+dirJ) + 1.0)
do qPt = 1, mesh_maxNips enddo
do dirI = 1, dimPlex enddo
mesh_ipCoordinates(dirI,qPt,cell+1) = pV0(dirI) qOffset = qOffset + dimPlex
do dirJ = 1, dimPlex enddo
mesh_ipCoordinates(dirI,qPt,cell+1) = mesh_ipCoordinates(dirI,qPt,cell+1) + & enddo
pCellJ((dirI-1)*dimPlex+dirJ)*(qPoints(qOffset+dirJ) + 1.0)
enddo
enddo
qOffset = qOffset + dimPlex
enddo
enddo
end subroutine mesh_FEM_build_ipCoordinates end subroutine mesh_FEM_build_ipCoordinates