now compiles with gfortran
This commit is contained in:
Martin Diehl
parent
f7c20d74af
commit
8fb780ab42
|
|
@ -63,7 +63,6 @@ use PETScDMplex
|
||||||
FEM_mech_destroy
|
FEM_mech_destroy
|
||||||
|
|
||||||
external :: &
|
external :: &
|
||||||
MPI_Allreduce, &
|
|
||||||
MatZeroRowsColumnsLocalIS, &
|
MatZeroRowsColumnsLocalIS, &
|
||||||
PetscQuadratureCreate, &
|
PetscQuadratureCreate, &
|
||||||
PetscFECreateDefault, &
|
PetscFECreateDefault, &
|
||||||
|
|
@ -189,7 +188,7 @@ subroutine FEM_mech_init(fieldBC)
|
||||||
do field = 1, dimPlex; do faceSet = 1, mesh_Nboundaries
|
do field = 1, dimPlex; do faceSet = 1, mesh_Nboundaries
|
||||||
if (fieldBC%componentBC(field)%Mask(faceSet)) then
|
if (fieldBC%componentBC(field)%Mask(faceSet)) then
|
||||||
numBC = numBC + 1
|
numBC = numBC + 1
|
||||||
call ISCreateGeneral(PETSC_COMM_WORLD,1,field-1,PETSC_COPY_VALUES,bcComps(numBC),ierr)
|
call ISCreateGeneral(PETSC_COMM_WORLD,1,[field-1],PETSC_COPY_VALUES,bcComps(numBC),ierr)
|
||||||
CHKERRQ(ierr)
|
CHKERRQ(ierr)
|
||||||
call DMGetStratumSize(mech_mesh,'Face Sets',mesh_boundaries(faceSet),bcSize,ierr)
|
call DMGetStratumSize(mech_mesh,'Face Sets',mesh_boundaries(faceSet),bcSize,ierr)
|
||||||
CHKERRQ(ierr)
|
CHKERRQ(ierr)
|
||||||
|
|
|
||||||
446
src/FEM_mesh.f90
446
src/FEM_mesh.f90
|
|
@ -1,446 +0,0 @@
|
||||||
!--------------------------------------------------------------------------------------------------
|
|
||||||
!> @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
|
|
||||||
use, intrinsic :: iso_c_binding
|
|
||||||
use prec, only: pReal, pInt
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
#include <petsc/finclude/petsc.h90>
|
|
||||||
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_abort, &
|
|
||||||
MPI_Bcast, &
|
|
||||||
DMClone, &
|
|
||||||
DMGetDimension, &
|
|
||||||
DMPlexCreateFromFile, &
|
|
||||||
DMPlexDistribute, &
|
|
||||||
DMPlexCopyCoordinates, &
|
|
||||||
DMGetStratumSize, &
|
|
||||||
DMPlexGetHeightStratum, &
|
|
||||||
DMPlexGetLabelValue, &
|
|
||||||
DMPlexSetLabelValue, &
|
|
||||||
DMDestroy
|
|
||||||
|
|
||||||
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
|
|
||||||
|
|
||||||
|
|
||||||
if (worldrank == 0) then
|
|
||||||
write(6,'(/,a)') ' <<<+- mesh init -+>>>'
|
|
||||||
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
|
|
||||||
#include "compilation_info.f90"
|
|
||||||
endif
|
|
||||||
|
|
||||||
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
|
|
||||||
|
|
@ -141,7 +141,6 @@ use PETScis
|
||||||
COMPONENT_MGTWIN_PHI_ID
|
COMPONENT_MGTWIN_PHI_ID
|
||||||
|
|
||||||
external :: &
|
external :: &
|
||||||
MPI_Allreduce, &
|
|
||||||
PetscOptionsInsertString, &
|
PetscOptionsInsertString, &
|
||||||
PetscObjectSetName, &
|
PetscObjectSetName, &
|
||||||
DMPlexGetHeightStratum, &
|
DMPlexGetHeightStratum, &
|
||||||
|
|
|
||||||
|
|
@ -97,7 +97,6 @@ use PETScis
|
||||||
mesh_cellCenterCoordinates
|
mesh_cellCenterCoordinates
|
||||||
|
|
||||||
external :: &
|
external :: &
|
||||||
MPI_Bcast, &
|
|
||||||
DMPlexCreateFromFile, &
|
DMPlexCreateFromFile, &
|
||||||
DMPlexDistribute, &
|
DMPlexDistribute, &
|
||||||
DMPlexCopyCoordinates, &
|
DMPlexCopyCoordinates, &
|
||||||
|
|
|
||||||
Loading…
Reference in New Issue