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

3493 lines
123 KiB
Fortran
Raw Blame History

! Copyright 2011 Max-Planck-Institut für Eisenforschung GmbH
!
! This file is part of DAMASK,
! the Düsseldorf Advanced MAterial Simulation Kit.
!
! DAMASK is free software: you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation, either version 3 of the License, or
! (at your option) any later version.
!
! DAMASK is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License
! along with DAMASK. If not, see <http://www.gnu.org/licenses/>.
!
!##############################################################
!* $Id$
!##############################################################
MODULE mesh
!##############################################################
use prec, only: pReal,pInt
implicit none
! ---------------------------
! _Nelems : total number of elements in mesh
! _NcpElems : total number of CP elements in mesh
! _Nnodes : total number of nodes in mesh
! _maxNnodes : max number of nodes in any CP element
! _maxNips : max number of IPs in any CP element
! _maxNipNeighbors : max number of IP neighbors in any CP element
! _maxNsharedElems : max number of CP elements sharing a node
!
! _element : FEid, type(internal representation), material, texture, node indices
! _node0 : x,y,z coordinates (initially!)
! _node : x,y,z coordinates (after deformation!)
! _sharedElem : entryCount and list of elements containing node
!
! _mapFEtoCPelem : [sorted FEid, corresponding CPid]
! _mapFEtoCPnode : [sorted FEid, corresponding CPid]
!
! MISSING: these definitions should actually reside in the
! FE-solver specific part (different for MARC/ABAQUS)..!
! Hence, I suggest to prefix with "FE_"
!
! _Nnodes : # nodes in a specific type of element (how we use it)
! _NoriginalNodes : # nodes in a specific type of element (how it is originally defined by marc)
! _Nips : # IPs in a specific type of element
! _NipNeighbors : # IP neighbors in a specific type of element
! _ipNeighbor : +x,-x,+y,-y,+z,-z list of intra-element IPs and
! (negative) neighbor faces per own IP in a specific type of element
! _NfaceNodes : # nodes per face in a specific type of element
! _nodeOnFace : list of node indices on each face of a specific type of element
! _maxNnodesAtIP : max number of (equivalent) nodes attached to an IP
! _nodesAtIP : map IP index to two node indices in a specific type of element
! _NsubNodes : # subnodes required to fully define all IP volumes
! order is +x,-x,+y,-y,+z,-z but meaning strongly depends on Elemtype
! ---------------------------
integer(pInt) mesh_Nelems, mesh_NcpElems, mesh_NelemSets, mesh_maxNelemInSet
integer(pInt) mesh_Nmaterials
integer(pInt) mesh_Nnodes, mesh_maxNnodes, mesh_maxNips, mesh_maxNipNeighbors, mesh_maxNsharedElems, mesh_maxNsubNodes
integer(pInt), dimension(2) :: mesh_maxValStateVar = 0_pInt
character(len=64), dimension(:), allocatable :: mesh_nameElemSet, & ! names of elementSet
mesh_nameMaterial, & ! names of material in solid section
mesh_mapMaterial ! name of elementSet for material
integer(pInt), dimension(:,:), allocatable :: mesh_mapElemSet ! list of elements in elementSet
integer(pInt), dimension(:,:), allocatable, target :: mesh_mapFEtoCPelem, mesh_mapFEtoCPnode
integer(pInt), dimension(:,:), allocatable :: mesh_element, & ! FEid, type(internal representation), material, texture, node indices
mesh_sharedElem, & ! entryCount and list of elements containing node
mesh_nodeTwins ! node twins are surface nodes that lie exactly on opposite sides of the mesh (surfaces nodes with equal coordinate values in two dimensions)
integer(pInt), dimension(:,:,:,:), allocatable :: mesh_ipNeighborhood ! 6 or less neighboring IPs as [element_num, IP_index]
real(pReal), dimension(:,:,:), allocatable :: mesh_subNodeCoord ! coordinates of subnodes per element
real(pReal), dimension(:,:), allocatable :: mesh_node0, & ! node coordinates (initially!)
mesh_node, & ! node coordinates (after deformation! ONLY FOR MARC!!!)
mesh_ipVolume ! volume associated with IP (initially!)
real(pReal), dimension(:,:,:), allocatable :: mesh_ipArea, & ! area of interface to neighboring IP (initially!)
mesh_ipCenterOfGravity ! center of gravity of IP (after deformation!)
real(pReal), dimension(:,:,:,:), allocatable :: mesh_ipAreaNormal ! area normal of interface to neighboring IP (initially!)
integer(pInt), dimension(:,:,:), allocatable :: FE_nodesAtIP ! map IP index to two node indices in a specific type of element
integer(pInt), dimension(:,:,:), allocatable :: FE_ipNeighbor
integer(pInt), dimension(:,:,:), allocatable :: FE_subNodeParent
integer(pInt), dimension(:,:,:,:), allocatable :: FE_subNodeOnIPFace
logical :: noPart ! for cases where the ABAQUS input file does not use part/assembly information
logical, dimension(3) :: mesh_periodicSurface ! flag indicating periodic outer surfaces (used for fluxes)
integer(pInt) :: hypoelasticTableStyle
integer(pInt) :: initialcondTableStyle
integer(pInt), parameter :: FE_Nelemtypes = 10
integer(pInt), parameter :: FE_maxNnodes = 8
integer(pInt), parameter :: FE_maxNsubNodes = 56
integer(pInt), parameter :: FE_maxNips = 27
integer(pInt), parameter :: FE_maxNipNeighbors = 6
integer(pInt), parameter :: FE_maxmaxNnodesAtIP = 8
integer(pInt), parameter :: FE_NipFaceNodes = 4
integer(pInt), dimension(FE_Nelemtypes), parameter :: FE_Nnodes = &
(/8, & ! element 7
4, & ! element 134
4, & ! element 11
4, & ! element 27
4, & ! element 157
6, & ! element 136
8, & ! element 21
8, & ! element 117
8, & ! element 57 (c3d20r == c3d8 --> copy of 7)
3 & ! element 155, 125, 128
/)
integer(pInt), dimension(FE_Nelemtypes), parameter :: FE_NoriginalNodes = &
(/8, & ! element 7
4, & ! element 134
4, & ! element 11
8, & ! element 27
4, & ! element 157
6, & ! element 136
20,& ! element 21
8, & ! element 117
20,& ! element 57 (c3d20r == c3d8 --> copy of 7)
6 & ! element 155, 125, 128
/)
integer(pInt), dimension(FE_Nelemtypes), parameter :: FE_Nips = &
(/8, & ! element 7
1, & ! element 134
4, & ! element 11
9, & ! element 27
4, & ! element 157
6, & ! element 136
27,& ! element 21
1, & ! element 117
8, & ! element 57 (c3d20r == c3d8 --> copy of 7)
3 & ! element 155, 125, 128
/)
integer(pInt), dimension(FE_Nelemtypes), parameter :: FE_NipNeighbors = &
(/6, & ! element 7
4, & ! element 134
4, & ! element 11
4, & ! element 27
6, & ! element 157
6, & ! element 136
6, & ! element 21
6, & ! element 117
6, & ! element 57 (c3d20r == c3d8 --> copy of 7)
4 & ! element 155, 125, 128
/)
integer(pInt), dimension(FE_Nelemtypes), parameter :: FE_NsubNodes = &
(/19,& ! element 7
0, & ! element 134
5, & ! element 11
12,& ! element 27
0, & ! element 157
15,& ! element 136
56,& ! element 21
0, & ! element 117
19,& ! element 57 (c3d20r == c3d8 --> copy of 7)
4 & ! element 155, 125, 128
/)
integer(pInt), dimension(FE_maxNipNeighbors,FE_Nelemtypes), parameter :: FE_NfaceNodes = &
reshape((/&
4,4,4,4,4,4, & ! element 7
3,3,3,3,0,0, & ! element 134
2,2,2,2,0,0, & ! element 11
2,2,2,2,0,0, & ! element 27
3,3,3,3,0,0, & ! element 157
3,4,4,4,3,0, & ! element 136
4,4,4,4,4,4, & ! element 21
4,4,4,4,4,4, & ! element 117
4,4,4,4,4,4, & ! element 57 (c3d20r == c3d8 --> copy of 7)
2,2,2,0,0,0 & ! element 155, 125, 128
/),(/FE_maxNipNeighbors,FE_Nelemtypes/))
integer(pInt), dimension(FE_Nelemtypes), parameter :: FE_maxNnodesAtIP = &
(/1, & ! element 7
4, & ! element 134
1, & ! element 11
2, & ! element 27
1, & ! element 157
1, & ! element 136
4, & ! element 21
8, & ! element 117
1, & ! element 57 (c3d20r == c3d8 --> copy of 7)
1 & ! element 155, 125, 128
/)
integer(pInt), dimension(FE_NipFaceNodes,FE_maxNipNeighbors,FE_Nelemtypes), parameter :: FE_nodeOnFace = &
reshape((/&
1,2,3,4 , & ! element 7
2,1,5,6 , &
3,2,6,7 , &
4,3,7,8 , &
4,1,5,8 , &
8,7,6,5 , &
1,2,3,0 , & ! element 134
1,4,2,0 , &
2,3,4,0 , &
1,3,4,0 , &
0,0,0,0 , &
0,0,0,0 , &
1,2,0,0 , & ! element 11
2,3,0,0 , &
3,4,0,0 , &
4,1,0,0 , &
0,0,0,0 , &
0,0,0,0 , &
1,2,0,0 , & ! element 27
2,3,0,0 , &
3,4,0,0 , &
4,1,0,0 , &
0,0,0,0 , &
0,0,0,0 , &
1,2,3,0 , & ! element 157
1,4,2,0 , &
2,3,4,0 , &
1,3,4,0 , &
0,0,0,0 , &
0,0,0,0 , &
1,2,3,0 , & ! element 136
1,4,5,2 , &
2,5,6,3 , &
1,3,6,4 , &
4,6,5,0 , &
0,0,0,0 , &
1,2,3,4 , & ! element 21
2,1,5,6 , &
3,2,6,7 , &
4,3,7,8 , &
4,1,5,8 , &
8,7,6,5 , &
1,2,3,4 , & ! element 117
2,1,5,6 , &
3,2,6,7 , &
4,3,7,8 , &
4,1,5,8 , &
8,7,6,5 , &
1,2,3,4 , & ! element 57 (c3d20r == c3d8 --> copy of 7)
2,1,5,6 , &
3,2,6,7 , &
4,3,7,8 , &
4,1,5,8 , &
8,7,6,5 , &
1,2,0,0 , & ! element 155,125,128
2,3,0,0 , &
3,1,0,0 , &
0,0,0,0 , &
0,0,0,0 , &
0,0,0,0 &
/),(/FE_NipFaceNodes,FE_maxNipNeighbors,FE_Nelemtypes/))
CONTAINS
! ---------------------------
! subroutine mesh_init()
! function mesh_FEtoCPelement(FEid)
! function mesh_build_ipNeighorhood()
! ---------------------------
!***********************************************************
! initialization
!***********************************************************
subroutine mesh_init (ip,element)
use DAMASK_interface
use prec, only: pInt
use IO, only: IO_error,IO_open_InputFile,IO_abaqus_hasNoPart
use FEsolving, only: parallelExecution, FEsolving_execElem, FEsolving_execIP, calcMode, lastMode, FEmodelGeometry
implicit none
integer(pInt), parameter :: fileUnit = 222
integer(pInt) e,element,ip
!$OMP CRITICAL (write2out)
write(6,*)
write(6,*) '<<<+- mesh init -+>>>'
write(6,*) '$Id$'
#include "compilation_info.f90"
!$OMP END CRITICAL (write2out)
call mesh_build_FEdata() ! --- get properties of the different types of elements
if (IO_open_inputFile(fileUnit,FEmodelGeometry)) then ! --- parse info from input file...
select case (FEsolver)
case ('Spectral')
call mesh_spectral_count_nodesAndElements(fileUnit)
call mesh_spectral_count_cpElements()
call mesh_spectral_map_elements()
call mesh_spectral_map_nodes()
call mesh_spectral_count_cpSizes()
call mesh_spectral_build_nodes(fileUnit)
call mesh_spectral_build_elements(fileUnit)
case ('Marc')
call mesh_marc_get_tableStyles(fileUnit)
call mesh_marc_count_nodesAndElements(fileUnit)
call mesh_marc_count_elementSets(fileUnit)
call mesh_marc_map_elementSets(fileUnit)
call mesh_marc_count_cpElements(fileUnit)
call mesh_marc_map_elements(fileUnit)
call mesh_marc_map_nodes(fileUnit)
call mesh_marc_build_nodes(fileUnit)
call mesh_marc_count_cpSizes(fileunit)
call mesh_marc_build_elements(fileUnit)
case ('Abaqus')
noPart = IO_abaqus_hasNoPart(fileUnit)
call mesh_abaqus_count_nodesAndElements(fileUnit)
call mesh_abaqus_count_elementSets(fileUnit)
call mesh_abaqus_count_materials(fileUnit)
call mesh_abaqus_map_elementSets(fileUnit)
call mesh_abaqus_map_materials(fileUnit)
call mesh_abaqus_count_cpElements(fileUnit)
call mesh_abaqus_map_elements(fileUnit)
call mesh_abaqus_map_nodes(fileUnit)
call mesh_abaqus_build_nodes(fileUnit)
call mesh_abaqus_count_cpSizes(fileunit)
call mesh_abaqus_build_elements(fileUnit)
end select
call mesh_get_damaskOptions(fileUnit)
close (fileUnit)
call mesh_build_subNodeCoords()
call mesh_build_ipCoordinates()
call mesh_build_ipVolumes()
call mesh_build_ipAreas()
call mesh_build_nodeTwins()
call mesh_build_sharedElems()
call mesh_build_ipNeighborhood()
call mesh_tell_statistics()
parallelExecution = (parallelExecution .and. (mesh_Nelems == mesh_NcpElems)) ! plus potential killer from non-local constitutive
else
call IO_error(error_ID=101_pInt) ! cannot open input file
endif
FEsolving_execElem = [ 1_pInt,mesh_NcpElems]
allocate(FEsolving_execIP(2,mesh_NcpElems)); FEsolving_execIP = 1_pInt
forall (e = 1_pInt:mesh_NcpElems) FEsolving_execIP(2,e) = FE_Nips(mesh_element(2,e))
allocate(calcMode(mesh_maxNips,mesh_NcpElems))
calcMode = .false. ! pretend to have collected what first call is asking (F = I)
calcMode(ip,mesh_FEasCP('elem',element)) = .true. ! first ip,el needs to be already pingponged to "calc"
lastMode = .true. ! and its mode is already known...
endsubroutine
!***********************************************************
! mapping of FE element types to internal representation
!***********************************************************
function FE_mapElemtype(what)
use IO, only: IO_lc
implicit none
character(len=*), intent(in) :: what
integer(pInt) FE_mapElemtype
select case (IO_lc(what))
case ( '7', &
'c3d8')
FE_mapElemtype = 1_pInt ! Three-dimensional Arbitrarily Distorted Brick
case ('134', &
'c3d4')
FE_mapElemtype = 2_pInt ! Three-dimensional Four-node Tetrahedron
case ( '11', &
'cpe4')
FE_mapElemtype = 3_pInt ! Arbitrary Quadrilateral Plane-strain
case ( '27', &
'cpe8')
FE_mapElemtype = 4_pInt ! Plane Strain, Eight-node Distorted Quadrilateral
case ('157')
FE_mapElemtype = 5_pInt ! Three-dimensional, Low-order, Tetrahedron, Herrmann Formulations
case ('136', &
'c3d6')
FE_mapElemtype = 6_pInt ! Three-dimensional Arbitrarily Distorted Pentahedral
case ( '21', &
'c3d20')
FE_mapElemtype = 7_pInt ! Three-dimensional Arbitrarily Distorted quadratic hexahedral
case ( '117', &
'123', &
'c3d8r')
FE_mapElemtype = 8_pInt ! Three-dimensional Arbitrarily Distorted linear hexahedral with reduced integration
case ( '57', &
'c3d20r')
FE_mapElemtype = 9_pInt ! Three-dimensional Arbitrarily Distorted quad hexahedral with reduced integration
case ( '155', &
'125', &
'128')
FE_mapElemtype = 10_pInt ! Two-dimensional Plane Strain triangle (155: cubic shape function, 125/128: second order isoparametric)
case default
FE_mapElemtype = 0_pInt ! unknown element --> should raise an error upstream..!
endselect
endfunction
!***********************************************************
! FE to CP id mapping by binary search thru lookup array
!
! valid questions are 'elem', 'node'
!***********************************************************
function mesh_FEasCP(what,id)
use prec, only: pInt
use IO, only: IO_lc
implicit none
character(len=*), intent(in) :: what
integer(pInt), intent(in) :: id
integer(pInt), dimension(:,:), pointer :: lookupMap
integer(pInt) mesh_FEasCP, 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)
! check at bounds QUESTION is it valid to extend bounds by 1 and just do binary search w/o init check at bounds?
if (lookupMap(1,lower) == id) then
mesh_FEasCP = lookupMap(2,lower)
return
elseif (lookupMap(1,upper) == id) then
mesh_FEasCP = lookupMap(2,upper)
return
endif
! binary search in between bounds
do while (upper-lower > 1)
center = (lower+upper)/2
if (lookupMap(1,center) < id) then
lower = center
elseif (lookupMap(1,center) > id) then
upper = center
else
mesh_FEasCP = lookupMap(2,center)
exit
endif
enddo
endfunction
!***********************************************************
! find face-matching element of same type
!***********************************************************
subroutine mesh_faceMatch(elem, face ,matchingElem, matchingFace)
use prec, only: pInt
implicit none
!*** output variables
integer(pInt), intent(out) :: matchingElem, & ! matching CP element ID
matchingFace ! matching FE face ID
!*** input variables
integer(pInt), intent(in) :: face, & ! FE face ID
elem ! FE elem ID
!*** local variables
integer(pInt), dimension(FE_NfaceNodes(face,mesh_element(2,elem))) :: &
myFaceNodes ! global node ids on my face
integer(pInt) myType, &
candidateType, &
candidateElem, &
candidateFace, &
candidateFaceNode, &
minNsharedElems, &
NsharedElems, &
lonelyNode, &
i, &
n, &
dir ! periodicity direction
integer(pInt), dimension(:), allocatable :: element_seen
logical checkTwins
minNsharedElems = mesh_maxNsharedElems + 1_pInt ! init to worst case
matchingFace = 0_pInt
matchingElem = 0_pInt ! intialize to "no match found"
myType = mesh_element(2,elem) ! figure elemType
do n = 1,FE_NfaceNodes(face,myType) ! loop over nodes on face
myFaceNodes(n) = mesh_FEasCP('node',mesh_element(4+FE_nodeOnFace(n,face,myType),elem)) ! CP id of face node
NsharedElems = mesh_sharedElem(1,myFaceNodes(n)) ! figure # shared elements for this node
if (NsharedElems < minNsharedElems) then
minNsharedElems = NsharedElems ! remember min # shared elems
lonelyNode = n ! remember most lonely node
endif
enddo
allocate(element_seen(minNsharedElems))
element_seen = 0_pInt
checkCandidate: do i = 1,minNsharedElems ! iterate over lonelyNode's shared elements
candidateElem = mesh_sharedElem(1+i,myFaceNodes(lonelyNode)) ! present candidate elem
if (all(element_seen /= candidateElem)) then ! element seen for the first time?
element_seen(i) = candidateElem
candidateType = mesh_element(2,candidateElem) ! figure elemType of candidate
checkCandidateFace: do candidateFace = 1,FE_maxNipNeighbors ! check each face of candidate
if (FE_NfaceNodes(candidateFace,candidateType) /= FE_NfaceNodes(face,myType) & ! incompatible face
.or. (candidateElem == elem .and. candidateFace == face)) then ! this is my face
cycle checkCandidateFace
endif
checkTwins = .false.
do n = 1,FE_NfaceNodes(candidateFace,candidateType) ! loop through nodes on face
candidateFaceNode = mesh_FEasCP('node', mesh_element(4+FE_nodeOnFace(n,candidateFace,candidateType),candidateElem))
if (all(myFaceNodes /= candidateFaceNode)) then ! candidate node does not match any of my face nodes
checkTwins = .true. ! perhaps the twin nodes do match
exit
endif
enddo
if(checkTwins) then
checkCandidateFaceTwins: do dir = 1,3
do n = 1,FE_NfaceNodes(candidateFace,candidateType) ! loop through nodes on face
candidateFaceNode = mesh_FEasCP('node', mesh_element(4+FE_nodeOnFace(n,candidateFace,candidateType),candidateElem))
if (all(myFaceNodes /= mesh_nodeTwins(dir,candidateFaceNode))) then ! node twin does not match either
if (dir == 3) then
cycle checkCandidateFace
else
cycle checkCandidateFaceTwins ! try twins in next dimension
endif
endif
enddo
exit checkCandidateFaceTwins
enddo checkCandidateFaceTwins
endif
matchingFace = candidateFace
matchingElem = candidateElem
exit checkCandidate ! found my matching candidate
enddo checkCandidateFace
endif
enddo checkCandidate
deallocate(element_seen)
endsubroutine
!********************************************************************
! get properties of different types of finite elements
!
! assign globals:
! FE_nodesAtIP, FE_ipNeighbor, FE_subNodeParent, FE_subNodeOnIPFace
!********************************************************************
subroutine mesh_build_FEdata ()
use prec, only: pInt
implicit none
allocate(FE_nodesAtIP(FE_maxmaxNnodesAtIP,FE_maxNips,FE_Nelemtypes)) ; FE_nodesAtIP = 0_pInt
allocate(FE_ipNeighbor(FE_maxNipNeighbors,FE_maxNips,FE_Nelemtypes)) ; FE_ipNeighbor = 0_pInt
allocate(FE_subNodeParent(FE_maxNips,FE_maxNsubNodes,FE_Nelemtypes)) ; FE_subNodeParent = 0_pInt
allocate(FE_subNodeOnIPFace(FE_NipFaceNodes,FE_maxNipNeighbors,FE_maxNips,FE_Nelemtypes)) ; FE_subNodeOnIPFace = 0_pInt
! fill FE_nodesAtIP with data
FE_nodesAtIP(:,:FE_Nips(1),1) = & ! element 7
reshape((/&
1, &
2, &
4, &
3, &
5, &
6, &
8, &
7 &
/),(/FE_maxNnodesAtIP(1),FE_Nips(1)/))
FE_nodesAtIP(:,:FE_Nips(2),2) = & ! element 134
reshape((/&
1,2,3,4 &
/),(/FE_maxNnodesAtIP(2),FE_Nips(2)/))
FE_nodesAtIP(:,:FE_Nips(3),3) = & ! element 11
reshape((/&
1, &
2, &
4, &
3 &
/),(/FE_maxNnodesAtIP(3),FE_Nips(3)/))
FE_nodesAtIP(:,:FE_Nips(4),4) = & ! element 27
reshape((/&
1,0, &
1,2, &
2,0, &
1,4, &
0,0, &
2,3, &
4,0, &
3,4, &
3,0 &
/),(/FE_maxNnodesAtIP(4),FE_Nips(4)/))
FE_nodesAtIP(:,:FE_Nips(5),5) = & ! element 157
reshape((/&
1, &
2, &
3, &
4 &
/),(/FE_maxNnodesAtIP(5),FE_Nips(5)/))
FE_nodesAtIP(:,:FE_Nips(6),6) = & ! element 136
reshape((/&
1, &
2, &
3, &
4, &
5, &
6 &
/),(/FE_maxNnodesAtIP(6),FE_Nips(6)/))
FE_nodesAtIP(:,:FE_Nips(7),7) = & ! element 21
reshape((/&
1,0, 0,0, &
1,2, 0,0, &
2,0, 0,0, &
1,4, 0,0, &
1,3, 2,4, &
2,3, 0,0, &
4,0, 0,0, &
3,4, 0,0, &
3,0, 0,0, &
1,5, 0,0, &
1,6, 2,5, &
2,6, 0,0, &
1,8, 4,5, &
0,0, 0,0, &
2,7, 3,6, &
4,8, 0,0, &
3,8, 4,7, &
3,7, 0,0, &
5,0, 0,0, &
5,6, 0,0, &
6,0, 0,0, &
5,8, 0,0, &
5,7, 6,8, &
6,7, 0,0, &
8,0, 0,0, &
7,8, 0,0, &
7,0, 0,0 &
/),(/FE_maxNnodesAtIP(7),FE_Nips(7)/))
! FE_nodesAtIP(:,:FE_Nips(8),8) = & ! element 117 (c3d8r --> single IP per element, so no need for this mapping)
! reshape((/&
! 1,2,3,4,5,6,7,8 &
! /),(/FE_maxNnodesAtIP(8),FE_Nips(8)/))
FE_nodesAtIP(:,:FE_Nips(9),9) = & ! element 57 (c3d20r == c3d8 --> copy of 7)
reshape((/&
1, &
2, &
4, &
3, &
5, &
6, &
8, &
7 &
/),(/FE_maxNnodesAtIP(9),FE_Nips(9)/))
FE_nodesAtIP(:,:FE_Nips(10),10) = & ! element 155, 125, 128
reshape((/&
1, &
2, &
3 &
/),(/FE_maxNnodesAtIP(10),FE_Nips(10)/))
! *** FE_ipNeighbor ***
! is a list of the neighborhood of each IP.
! It is sorted in (local) +x,-x, +y,-y, +z,-z direction.
! Positive integers denote an intra-FE IP identifier.
! Negative integers denote the interface behind which the neighboring (extra-FE) IP will be located.
FE_ipNeighbor(:FE_NipNeighbors(1),:FE_Nips(1),1) = & ! element 7
reshape((/&
2,-5, 3,-2, 5,-1, &
-3, 1, 4,-2, 6,-1, &
4,-5,-4, 1, 7,-1, &
-3, 3,-4, 2, 8,-1, &
6,-5, 7,-2,-6, 1, &
-3, 5, 8,-2,-6, 2, &
8,-5,-4, 5,-6, 3, &
-3, 7,-4, 6,-6, 4 &
/),(/FE_NipNeighbors(1),FE_Nips(1)/))
FE_ipNeighbor(:FE_NipNeighbors(2),:FE_Nips(2),2) = & ! element 134
reshape((/&
-1,-2,-3,-4 &
/),(/FE_NipNeighbors(2),FE_Nips(2)/))
FE_ipNeighbor(:FE_NipNeighbors(3),:FE_Nips(3),3) = & ! element 11
reshape((/&
2,-4, 3,-1, &
-2, 1, 4,-1, &
4,-4,-3, 1, &
-2, 3,-3, 2 &
/),(/FE_NipNeighbors(3),FE_Nips(3)/))
FE_ipNeighbor(:FE_NipNeighbors(4),:FE_Nips(4),4) = & ! element 27
reshape((/&
2,-4, 4,-1, &
3, 1, 5,-1, &
-2, 2, 6,-1, &
5,-4, 7, 1, &
6, 4, 8, 2, &
-2, 5, 9, 3, &
8,-4,-3, 4, &
9, 7,-3, 5, &
-2, 8,-3, 6 &
/),(/FE_NipNeighbors(4),FE_Nips(4)/))
FE_ipNeighbor(:FE_NipNeighbors(5),:FE_Nips(5),5) = & ! element 157
reshape((/&
2,-4, 3,-2, 4,-1, &
3,-2, 1,-3, 4,-1, &
1,-3, 2,-4, 4,-1, &
1,-3, 2,-4, 3,-2 &
/),(/FE_NipNeighbors(5),FE_Nips(5)/))
FE_ipNeighbor(:FE_NipNeighbors(6),:FE_Nips(6),6) = & ! element 136
reshape((/&
2,-4, 3,-2, 4,-1, &
-3, 1, 3,-2, 5,-1, &
2,-4,-3, 1, 6,-1, &
5,-4, 6,-2,-5, 1, &
-3, 4, 6,-2,-5, 2, &
5,-4,-3, 4,-5, 3 &
/),(/FE_NipNeighbors(6),FE_Nips(6)/))
FE_ipNeighbor(:FE_NipNeighbors(7),:FE_Nips(7),7) = & ! element 21
reshape((/&
2,-5, 4,-2,10,-1, &
3, 1, 5,-2,11,-1, &
-3, 2, 6,-2,12,-1, &
5,-5, 7, 1,13,-1, &
6, 4, 8, 2,14,-1, &
-3, 5, 9, 3,15,-1, &
8,-5,-4, 4,16,-1, &
9, 7,-4, 5,17,-1, &
-3, 8,-4, 6,18,-1, &
11,-5,13,-2,19, 1, &
12,10,14,-2,20, 2, &
-3,11,15,-2,21, 3, &
14,-5,16,10,22, 4, &
15,13,17,11,23, 5, &
-3,14,18,12,24, 6, &
17,-5,-4,13,25, 7, &
18,16,-4,14,26, 8, &
-3,17,-4,15,27, 9, &
20,-5,22,-2,-6,10, &
21,19,23,-2,-6,11, &
-3,20,24,-2,-6,12, &
23,-5,25,19,-6,13, &
24,22,26,20,-6,14, &
-3,23,27,21,-6,15, &
26,-5,-4,22,-6,16, &
27,25,-4,23,-6,17, &
-3,26,-4,24,-6,18 &
/),(/FE_NipNeighbors(7),FE_Nips(7)/))
FE_ipNeighbor(:FE_NipNeighbors(8),:FE_Nips(8),8) = & ! element 117
reshape((/&
-3,-5,-4,-2,-6,-1 &
/),(/FE_NipNeighbors(8),FE_Nips(8)/))
FE_ipNeighbor(:FE_NipNeighbors(9),:FE_Nips(9),9) = & ! element 57 (c3d20r == c3d8 --> copy of 7)
reshape((/&
2,-5, 3,-2, 5,-1, &
-3, 1, 4,-2, 6,-1, &
4,-5,-4, 1, 7,-1, &
-3, 3,-4, 2, 8,-1, &
6,-5, 7,-2,-6, 1, &
-3, 5, 8,-2,-6, 2, &
8,-5,-4, 5,-6, 3, &
-3, 7,-4, 6,-6, 4 &
/),(/FE_NipNeighbors(9),FE_Nips(9)/))
FE_ipNeighbor(:FE_NipNeighbors(10),:FE_Nips(10),10) = & ! element 155, 125, 128
reshape((/&
2,-3, 3,-1, &
-2, 1, 3,-1, &
2,-3,-2, 1 &
/),(/FE_NipNeighbors(10),FE_Nips(10)/))
! *** FE_subNodeParent ***
! lists the group of nodes for which the center of gravity
! corresponds to the location of a each subnode.
! example: face-centered subnode with faceNodes 1,2,3,4 to be used in,
! e.g., a 8 IP grid, would be encoded:
! 1, 2, 3, 4, 0, 0, 0, 0
FE_subNodeParent(:FE_Nips(1),:FE_NsubNodes(1),1) = & ! element 7
reshape((/&
1, 2, 0, 0, 0, 0, 0, 0, &
2, 3, 0, 0, 0, 0, 0, 0, &
3, 4, 0, 0, 0, 0, 0, 0, &
4, 1, 0, 0, 0, 0, 0, 0, &
1, 5, 0, 0, 0, 0, 0, 0, &
2, 6, 0, 0, 0, 0, 0, 0, &
3, 7, 0, 0, 0, 0, 0, 0, &
4, 8, 0, 0, 0, 0, 0, 0, &
5, 6, 0, 0, 0, 0, 0, 0, &
6, 7, 0, 0, 0, 0, 0, 0, &
7, 8, 0, 0, 0, 0, 0, 0, &
8, 5, 0, 0, 0, 0, 0, 0, &
1, 2, 3, 4, 0, 0, 0, 0, &
1, 2, 6, 5, 0, 0, 0, 0, &
2, 3, 7, 6, 0, 0, 0, 0, &
3, 4, 8, 7, 0, 0, 0, 0, &
1, 4, 8, 5, 0, 0, 0, 0, &
5, 6, 7, 8, 0, 0, 0, 0, &
1, 2, 3, 4, 5, 6, 7, 8 &
/),(/FE_Nips(1),FE_NsubNodes(1)/))
!FE_subNodeParent(:FE_Nips(2),:FE_NsubNodes(2),2) ! element 134 has no subnodes
FE_subNodeParent(:FE_Nips(3),:FE_NsubNodes(3),3) = & ! element 11
reshape((/&
1, 2, 0, 0, &
2, 3, 0, 0, &
3, 4, 0, 0, &
4, 1, 0, 0, &
1, 2, 3, 4 &
/),(/FE_Nips(3),FE_NsubNodes(3)/))
FE_subNodeParent(:FE_Nips(4),:FE_NsubNodes(4),4) = & ! element 27
reshape((/&
1, 1, 2, 0, 0, 0, 0, 0, 0, &
1, 2, 2, 0, 0, 0, 0, 0, 0, &
2, 2, 3, 0, 0, 0, 0, 0, 0, &
2, 3, 3, 0, 0, 0, 0, 0, 0, &
3, 3, 4, 0, 0, 0, 0, 0, 0, &
3, 4, 4, 0, 0, 0, 0, 0, 0, &
4, 4, 1, 0, 0, 0, 0, 0, 0, &
4, 1, 1, 0, 0, 0, 0, 0, 0, &
1, 1, 1, 1, 2, 2, 4, 4, 3, &
2, 2, 2, 2, 1, 1, 3, 3, 4, &
3, 3, 3, 3, 2, 2, 4, 4, 1, &
4, 4, 4, 4, 1, 1, 3, 3, 2 &
/),(/FE_Nips(4),FE_NsubNodes(4)/))
!FE_subNodeParent(:FE_Nips(5),:FE_NsubNodes(5),5) = & ! element 157
! reshape((/&
! *still to be defined*
! /),(/FE_Nips(5),FE_NsubNodes(5)/))
FE_subNodeParent(:FE_Nips(6),:FE_NsubNodes(6),6) = & ! element 136
reshape((/&
1, 2, 0, 0, 0, 0, &
2, 3, 0, 0, 0, 0, &
3, 1, 0, 0, 0, 0, &
1, 4, 0, 0, 0, 0, &
2, 5, 0, 0, 0, 0, &
3, 6, 0, 0, 0, 0, &
4, 5, 0, 0, 0, 0, &
5, 6, 0, 0, 0, 0, &
6, 4, 0, 0, 0, 0, &
1, 2, 3, 0, 0, 0, &
1, 2, 4, 5, 0, 0, &
2, 3, 5, 6, 0, 0, &
1, 3, 4, 6, 0, 0, &
4, 5, 6, 0, 0, 0, &
1, 2, 3, 4, 5, 6 &
/),(/FE_Nips(6),FE_NsubNodes(6)/))
FE_subNodeParent(:FE_Nips(7),:FE_NsubNodes(7),7) = & ! element 21
reshape((/&
1, 1, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
1, 2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
2, 2, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
2, 3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
3, 3, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
3, 4, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
4, 4, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
4, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
1, 1, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
2, 2, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
3, 3, 7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
4, 4, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
1, 5, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
2, 6, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
3, 7, 7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
4, 8, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
5, 5, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
5, 6, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
6, 6, 7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
6, 7, 7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
7, 7, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
7, 8, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
8, 8, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
8, 5, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
1, 1, 1, 1, 2, 2, 4, 4, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
2, 2, 2, 2, 1, 1, 3, 3, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
3, 3, 3, 3, 2, 2, 4, 4, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
4, 4, 4, 4, 1, 1, 3, 3, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
1, 1, 1, 1, 2, 2, 5, 5, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
2, 2, 2, 2, 1, 1, 6, 6, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
2, 2, 2, 2, 3, 3, 6, 6, 7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
3, 3, 3, 3, 2, 2, 7, 7, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
3, 3, 3, 3, 4, 4, 7, 7, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
4, 4, 4, 4, 3, 3, 8, 8, 7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
4, 4, 4, 4, 1, 1, 8, 8, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
1, 1, 1, 1, 4, 4, 5, 5, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
5, 5, 5, 5, 1, 1, 6, 6, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
6, 6, 6, 6, 2, 2, 5, 5, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
6, 6, 6, 6, 2, 2, 7, 7, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
7, 7, 7, 7, 3, 3, 6, 6, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
7, 7, 7, 7, 3, 3, 8, 8, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
8, 8, 8, 8, 4, 4, 7, 7, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
8, 8, 8, 8, 4, 4, 5, 5, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
5, 5, 5, 5, 1, 1, 8, 8, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
5, 5, 5, 5, 6, 6, 8, 8, 7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
6, 6, 6, 6, 5, 5, 7, 7, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
7, 7, 7, 7, 6, 6, 8, 8, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
8, 8, 8, 8, 5, 5, 7, 7, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 4, 4, 4, 4, 5, 5, 5, 5, 3, 3, 6, 6, 8, 8, 7, &
2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 3, 3, 3, 3, 6, 6, 6, 6, 4, 4, 5, 5, 7, 7, 8, &
3, 3, 3, 3, 3, 3, 3, 3, 2, 2, 2, 2, 4, 4, 4, 4, 7, 7, 7, 7, 1, 1, 6, 6, 8, 8, 5, &
4, 4, 4, 4, 4, 4, 4, 4, 1, 1, 1, 1, 3, 3, 3, 3, 8, 8, 8, 8, 2, 2, 5, 5, 7, 7, 6, &
5, 5, 5, 5, 5, 5, 5, 5, 1, 1, 1, 1, 6, 6, 6, 6, 8, 8, 8, 8, 2, 2, 4, 4, 7, 7, 3, &
6, 6, 6, 6, 6, 6, 6, 6, 2, 2, 2, 2, 5, 5, 5, 5, 7, 7, 7, 7, 1, 1, 3, 3, 8, 8, 4, &
7, 7, 7, 7, 7, 7, 7, 7, 3, 3, 3, 3, 6, 6, 6, 6, 8, 8, 8, 8, 2, 2, 4, 4, 5, 5, 1, &
8, 8, 8, 8, 8, 8, 8, 8, 4, 4, 4, 4, 5, 5, 5, 5, 7, 7, 7, 7, 1, 1, 3, 3, 6, 6, 2 &
/),(/FE_Nips(7),FE_NsubNodes(7)/))
!FE_subNodeParent(:FE_Nips(8),:FE_NsubNodes(8),8) ! element 117 has no subnodes
FE_subNodeParent(:FE_Nips(9),:FE_NsubNodes(9),9) = & ! element 57 (c3d20r == c3d8 --> copy of 7)
reshape((/&
1, 2, 0, 0, 0, 0, 0, 0, &
2, 3, 0, 0, 0, 0, 0, 0, &
3, 4, 0, 0, 0, 0, 0, 0, &
4, 1, 0, 0, 0, 0, 0, 0, &
1, 5, 0, 0, 0, 0, 0, 0, &
2, 6, 0, 0, 0, 0, 0, 0, &
3, 7, 0, 0, 0, 0, 0, 0, &
4, 8, 0, 0, 0, 0, 0, 0, &
5, 6, 0, 0, 0, 0, 0, 0, &
6, 7, 0, 0, 0, 0, 0, 0, &
7, 8, 0, 0, 0, 0, 0, 0, &
8, 5, 0, 0, 0, 0, 0, 0, &
1, 2, 3, 4, 0, 0, 0, 0, &
1, 2, 6, 5, 0, 0, 0, 0, &
2, 3, 7, 6, 0, 0, 0, 0, &
3, 4, 8, 7, 0, 0, 0, 0, &
1, 4, 8, 5, 0, 0, 0, 0, &
5, 6, 7, 8, 0, 0, 0, 0, &
1, 2, 3, 4, 5, 6, 7, 8 &
/),(/FE_Nips(9),FE_NsubNodes(9)/))
FE_subNodeParent(:FE_Nips(10),:FE_NsubNodes(10),10) = & ! element 155, 125, 128
reshape((/&
1, 2, 0, &
2, 3, 0, &
3, 1, 0, &
1, 2, 3 &
/),(/FE_Nips(10),FE_NsubNodes(10)/))
! *** FE_subNodeOnIPFace ***
! indicates which subnodes make up the interfaces enclosing the IP volume.
! The sorting convention is such that the outward pointing normal
! follows from a right-handed traversal of the face node list.
! For two-dimensional elements, which only have lines as "interface"
! one nevertheless has to specify each interface by a closed path,
! e.g., 1,2, 2,1, assuming the line connects nodes 1 and 2.
! This will result in zero ipVolume and interfaceArea, but is not
! detrimental at the moment since non-local constitutive laws are
! currently not foreseen in 2D cases.
FE_subNodeOnIPFace(:FE_NipFaceNodes,:FE_NipNeighbors(1),:FE_Nips(1),1) = & ! element 7
reshape((/&
9,21,27,22, & ! 1
1,13,25,12, &
12,25,27,21, &
1, 9,22,13, &
13,22,27,25, &
1,12,21, 9, &
2,10,23,14, & ! 2
9,22,27,21, &
10,21,27,23, &
2,14,22, 9, &
14,23,27,22, &
2, 9,21,10, &
11,24,27,21, & ! 3
4,12,25,16, &
4,16,24,11, &
12,21,27,25, &
16,25,27,24, &
4,11,21,12, &
3,15,23,10, & ! 4
11,21,27,24, &
3,11,24,15, &
10,23,27,21, &
15,24,27,23, &
3,10,21,11, &
17,22,27,26, & ! 5
5,20,25,13, &
20,26,27,25, &
5,13,22,17, &
5,17,26,20, &
13,25,27,22, &
6,14,23,18, & ! 6
17,26,27,22, &
18,23,27,26, &
6,17,22,14, &
6,18,26,17, &
14,22,27,23, &
19,26,27,24, & ! 7
8,16,25,20, &
8,19,24,16, &
20,25,27,26, &
8,20,26,19, &
16,24,27,25, &
7,18,23,15, & ! 8
19,24,27,26, &
7,15,24,19, &
18,26,27,23, &
7,19,26,18, &
15,23,27,24 &
/),(/FE_NipFaceNodes,FE_NipNeighbors(1),FE_Nips(1)/))
FE_subNodeOnIPFace(:FE_NipFaceNodes,:FE_NipNeighbors(2),:FE_Nips(2),2) = & ! element 134
reshape((/&
1, 1, 3, 2, & ! 1
1, 1, 2, 4, &
2, 2, 3, 4, &
1, 1, 4, 3 &
/),(/FE_NipFaceNodes,FE_NipNeighbors(2),FE_Nips(2)/))
FE_subNodeOnIPFace(:FE_NipFaceNodes,:FE_NipNeighbors(3),:FE_Nips(3),3) = & ! element 11
reshape((/&
5, 9, 9, 5 , & ! 1
1, 8, 8, 1 , &
8, 9, 9, 8 , &
1, 5, 5, 1 , &
2, 6, 6, 2 , & ! 2
5, 9, 9, 5 , &
6, 9, 9, 6 , &
2, 5, 5, 2 , &
3, 6, 6, 3 , & ! 3
7, 9, 9, 7 , &
3, 7, 7, 3 , &
6, 9, 9, 6 , &
7, 9, 9, 7 , & ! 4
4, 8, 8, 4 , &
4, 7, 7, 4 , &
8, 9, 9, 8 &
/),(/FE_NipFaceNodes,FE_NipNeighbors(3),FE_Nips(3)/))
FE_subNodeOnIPFace(:FE_NipFaceNodes,:FE_NipNeighbors(4),:FE_Nips(4),4) = & ! element 27
reshape((/&
9,17,17, 9 , & ! 1
1,16,16, 1 , &
16,17,17,16 , &
1, 9, 9, 1 , &
10,18,18,10 , & ! 2
9,17,17, 9 , &
17,18,18,17 , &
9,10,10, 9 , &
2,11,11, 2 , & ! 3
10,18,18,10 , &
11,18,18,11 , &
2,10,10, 2 , &
17,20,20,17 , & ! 4
15,16,16,15 , &
15,20,20,15 , &
16,17,17,16 , &
18,19,19,18 , & ! 5
17,20,20,17 , &
19,20,20,19 , &
17,18,18,17 , &
11,12,12,11 , & ! 6
18,19,19,18 , &
12,19,19,12 , &
11,18,18,11 , &
14,20,20,14 , & ! 7
4,15,15, 4 , &
4,14,14, 4 , &
15,20,20,15 , &
13,19,19,13 , & ! 8
14,20,20,14 , &
13,14,14,13 , &
19,20,20,19 , &
3,12,12, 3 , & ! 9
13,19,19,13 , &
3,13,13, 3 , &
12,19,19,12 &
/),(/FE_NipFaceNodes,FE_NipNeighbors(4),FE_Nips(4)/))
!FE_subNodeOnIPFace(:FE_NipFaceNodes,:FE_NipNeighbors(5),:FE_Nips(5),5) = & ! element 157
! reshape((/&
! *still to be defined*
! /),(/FE_NipFaceNodes,FE_NipNeighbors(5),FE_Nips(5)/))
FE_subNodeOnIPFace(:FE_NipFaceNodes,:FE_NipNeighbors(6),:FE_Nips(6),6) = & ! element 136
reshape((/&
7,16,21,17, & ! 1
1,10,19, 9, &
9,19,21,16, &
1, 7,17,10, &
10,17,21,19, &
1, 9,16, 7, &
2, 8,18,11, & ! 2
7,17,21,16, &
8,16,21,18, &
2,11,17, 7, &
11,18,21,17, &
2, 7,16, 8, &
8,18,21,16, & ! 3
3, 9,19,12, &
3,12,18, 8, &
9,16,21,19, &
12,19,21,18, &
3, 8,16, 9, &
13,17,21,20, & ! 4
4,15,19,10, &
15,20,21,19, &
4,10,17,13, &
4,13,20,15, &
10,19,21,17, &
5,11,18,14, & ! 5
13,20,21,17, &
14,18,21,20, &
5,13,17,11, &
5,14,20,13, &
11,17,21,18, &
14,20,21,18, & ! 6
6,12,19,15, &
6,14,18,12, &
15,19,21,20, &
6,15,20,14, &
12,18,21,19 &
/),(/FE_NipFaceNodes,FE_NipNeighbors(6),FE_Nips(6)/))
FE_subNodeOnIPFace(:FE_NipFaceNodes,:FE_NipNeighbors(7),:FE_Nips(7),7) = & ! element 21
reshape((/&
9,33,57,37, & ! 1
1,17,44,16, &
33,16,44,57, &
1, 9,37,17, &
17,37,57,44, &
1,16,33, 9, &
10,34,58,38, & ! 2
9,37,57,33, &
34,33,57,58, &
9,10,38,37, &
37,38,58,57, &
9,33,34,10, &
2,11,39,18, & ! 3
10,38,58,34, &
11,34,58,39, &
10, 2,18,38, &
38,18,39,58, &
10,34,11, 2, &
33,36,60,57, & ! 4
16,44,43,15, &
36,15,43,60, &
16,33,57,44, &
44,57,60,43, &
16,15,36,33, &
34,35,59,58, & ! 5
33,57,60,36, &
35,36,60,59, &
33,34,58,57, &
57,58,59,60, &
33,36,35,34, &
11,12,40,39, & ! 6
34,58,59,35, &
12,35,59,40, &
34,11,39,58, &
58,39,40,59, &
34,35,12,11, &
36,14,42,60, & ! 7
15,43,20, 4, &
14, 4,20,42, &
15,36,60,43, &
43,60,42,20, &
15, 4,14,36, &
35,13,41,59, & ! 8
36,60,42,14, &
13,14,42,41, &
36,35,59,60, &
60,59,41,42, &
36,14,13,35, &
12, 3,19,40, & ! 9
35,59,41,13, &
3,13,41,19, &
35,12,40,59, &
59,40,19,41, &
35,13, 3,12, &
37,57,61,45, & ! 10
17,21,52,44, &
57,44,52,61, &
17,37,45,21, &
21,45,61,52, &
17,44,57,37, &
38,58,62,46, & ! 11
37,45,61,57, &
58,57,61,62, &
37,38,46,45, &
45,46,62,61, &
37,57,58,38, &
18,39,47,22, & ! 12
38,46,62,58, &
39,58,62,47, &
38,18,22,46, &
46,22,47,62, &
38,58,39,18, &
57,60,64,61, & ! 13
44,52,51,43, &
60,43,51,64, &
44,57,61,52, &
52,61,64,51, &
44,43,60,57, &
58,59,63,62, & ! 14
57,61,64,60, &
59,60,64,63, &
57,58,62,61, &
61,62,63,64, &
57,60,59,58, &
39,40,48,47, & ! 15
58,62,63,59, &
40,59,63,48, &
58,39,47,62, &
62,47,48,63, &
58,59,40,39, &
60,42,50,64, & ! 16
43,51,24,20, &
42,20,24,50, &
43,60,64,51, &
51,64,50,24, &
43,20,42,60, &
59,41,49,63, & ! 17
60,64,50,42, &
41,42,50,49, &
60,59,63,64, &
64,63,49,50, &
60,42,41,59, &
40,19,23,48, & ! 18
59,63,49,41, &
19,41,49,23, &
59,40,48,63, &
63,48,23,49, &
59,41,19,40, &
45,61,53,25, & ! 19
21, 5,32,52, &
61,52,32,53, &
21,45,25, 5, &
5,25,53,32, &
21,52,61,45, &
46,62,54,26, & ! 20
45,25,53,61, &
62,61,53,54, &
45,46,26,25, &
25,26,54,53, &
45,61,62,46, &
22,47,27, 6, & ! 21
46,26,54,62, &
47,62,54,27, &
46,22, 6,26, &
26, 6,27,54, &
46,62,47,22, &
61,64,56,53, & ! 22
52,32,31,51, &
64,51,31,56, &
52,61,53,32, &
32,53,56,31, &
52,51,64,61, &
62,63,55,54, & ! 23
61,53,56,64, &
63,64,56,55, &
61,62,54,53, &
53,54,55,56, &
61,64,63,62, &
47,48,28,27, & ! 24
62,54,55,63, &
48,63,55,28, &
62,47,27,54, &
54,27,28,55, &
62,63,48,47, &
64,50,30,56, & ! 25
51,31, 8,24, &
50,24, 8,30, &
51,64,56,31, &
31,56,30, 8, &
51,24,50,64, &
63,49,29,55, & ! 26
64,56,30,50, &
49,50,30,29, &
64,63,55,56, &
56,55,29,30, &
64,50,49,63, &
48,23, 7,28, & ! 27
63,55,29,49, &
23,49,29, 7, &
63,48,28,55, &
55,28, 7,29, &
63,49,23,48 &
/),(/FE_NipFaceNodes,FE_NipNeighbors(7),FE_Nips(7)/))
FE_subNodeOnIPFace(:FE_NipFaceNodes,:FE_NipNeighbors(8),:FE_Nips(8),8) = & ! element 117
reshape((/&
2, 3, 7, 6, & ! 1
1, 5, 8, 4, &
3, 4, 8, 7, &
1, 2, 6, 5, &
5, 6, 7, 8, &
1, 4, 3, 2 &
/),(/FE_NipFaceNodes,FE_NipNeighbors(8),FE_Nips(8)/))
FE_subNodeOnIPFace(:FE_NipFaceNodes,:FE_NipNeighbors(9),:FE_Nips(9),9) = & ! element 57 (c3d20r == c3d8 --> copy of 7)
reshape((/&
9,21,27,22, & ! 1
1,13,25,12, &
12,25,27,21, &
1, 9,22,13, &
13,22,27,25, &
1,12,21, 9, &
2,10,23,14, & ! 2
9,22,27,21, &
10,21,27,23, &
2,14,22, 9, &
14,23,27,22, &
2, 9,21,10, &
11,24,27,21, & ! 3
4,12,25,16, &
4,16,24,11, &
12,21,27,25, &
16,25,27,24, &
4,11,21,12, &
3,15,23,10, & ! 4
11,21,27,24, &
3,11,24,15, &
10,23,27,21, &
15,24,27,23, &
3,10,21,11, &
17,22,27,26, & ! 5
5,20,25,13, &
20,26,27,25, &
5,13,22,17, &
5,17,26,20, &
13,25,27,22, &
6,14,23,18, & ! 6
17,26,27,22, &
18,23,27,26, &
6,17,22,14, &
6,18,26,17, &
14,22,27,23, &
19,26,27,24, & ! 7
8,16,25,20, &
8,19,24,16, &
20,25,27,26, &
8,20,26,19, &
16,24,27,25, &
7,18,23,15, & ! 8
19,24,27,26, &
7,15,24,19, &
18,26,27,23, &
7,19,26,18, &
15,23,27,24 &
/),(/FE_NipFaceNodes,FE_NipNeighbors(9),FE_Nips(9)/))
FE_subNodeOnIPFace(:FE_NipFaceNodes,:FE_NipNeighbors(10),:FE_Nips(10),10) = & ! element 155, 125, 128
reshape((/&
4, 7, 7, 4 , & ! 1
1, 6, 6, 1 , &
6, 7, 7, 6 , &
1, 4, 4, 1 , &
2, 5, 5, 2 , & ! 2
4, 7, 7, 4 , &
5, 7, 7, 5 , &
2, 4, 4, 2 , &
5, 7, 7, 5 , & ! 3
3, 6, 6, 3 , &
3, 5, 5, 3 , &
6, 7, 7, 6 &
/),(/FE_NipFaceNodes,FE_NipNeighbors(10),FE_Nips(10)/))
endsubroutine
!********************************************************************
! figure out table styles (Marc only)
!
! initialcondTableStyle, hypoelasticTableStyle
!********************************************************************
subroutine mesh_marc_get_tableStyles (myUnit)
use prec, only: pInt
use IO
implicit none
integer(pInt), parameter :: maxNchunks = 6
integer(pInt), dimension (1+2*maxNchunks) :: myPos
integer(pInt) myUnit
character(len=300) line
initialcondTableStyle = 0_pInt
hypoelasticTableStyle = 0_pInt
610 FORMAT(A300)
rewind(myUnit)
do
read (myUnit,610,END=620) line
myPos = IO_stringPos(line,maxNchunks)
if ( IO_lc(IO_stringValue(line,myPos,1_pInt)) == 'table' .and. myPos(1_pInt) .GT. 5) then
initialcondTableStyle = IO_intValue(line,myPos,4_pInt)
hypoelasticTableStyle = IO_intValue(line,myPos,5_pInt)
exit
endif
enddo
620 endsubroutine
!********************************************************************
! get any additional damask options from input file
!
! mesh_periodicSurface
!********************************************************************
subroutine mesh_get_damaskOptions(myUnit)
use DAMASK_interface, only: FEsolver
use prec, only: pInt
use IO
implicit none
integer(pInt), intent(in) :: myUnit
integer(pInt), parameter :: maxNchunks = 5
integer(pInt), dimension (1+2*maxNchunks) :: myPos
integer(pInt) chunk, Nchunks
character(len=300) line, keyword, damaskOption, v
mesh_periodicSurface = (/.false., .false., .false./)
610 FORMAT(A300)
select case (FEsolver)
case ('Spectral') ! no special keyword needed, the damask option directly goes into the header
case ('Marc')
keyword = '$damask'
case ('Abaqus')
keyword = '**damask'
end select
rewind(myUnit)
do
read (myUnit,610,END=620) line
myPos = IO_stringPos(line,maxNchunks)
Nchunks = myPos(1)
select case (FEsolver)
case ('Marc','Abaqus')
if (IO_lc(IO_stringValue(line,myPos,1_pInt)) == keyword .and. Nchunks > 1_pInt) then ! found keyword for damask option and there is at least one more chunk to read
damaskOption = IO_lc(IO_stringValue(line,myPos,2_pInt))
select case(damaskOption)
case('periodic') ! damask Option that allows to specify periodic fluxes
do chunk = 3_pInt,Nchunks ! loop through chunks (skipping the keyword)
v = IO_lc(IO_stringValue(line,myPos,chunk)) ! chunk matches keyvalues x,y, or z?
mesh_periodicSurface(1) = mesh_periodicSurface(1) .or. v == 'x'
mesh_periodicSurface(2) = mesh_periodicSurface(2) .or. v == 'y'
mesh_periodicSurface(3) = mesh_periodicSurface(3) .or. v == 'z'
enddo
endselect
endif
case('Spectral')
damaskOption = IO_lc(IO_stringValue(line,myPos,1_pInt))
select case(damaskOption)
case('periodic') ! damask Option that allows to specify periodic fluxes
do chunk = 2_pInt,Nchunks ! loop through chunks (skipping the keyword)
v = IO_lc(IO_stringValue(line,myPos,chunk)) ! chunk matches keyvalues x,y, or z?
mesh_periodicSurface(1) = mesh_periodicSurface(1) .or. v == 'x'
mesh_periodicSurface(2) = mesh_periodicSurface(2) .or. v == 'y'
mesh_periodicSurface(3) = mesh_periodicSurface(3) .or. v == 'z'
enddo
endselect
endselect
enddo
620 endsubroutine
!********************************************************************
! count overall number of nodes and elements in mesh
!
! mesh_Nelems, mesh_Nnodes
!********************************************************************
subroutine mesh_spectral_count_nodesAndElements (myUnit)
use prec, only: pInt
use IO
implicit none
integer(pInt), parameter :: maxNchunks = 7
integer(pInt), dimension (1+2*maxNchunks) :: myPos
integer(pInt) a,b,c,i,j,headerLength
integer(pInt) myUnit
character(len=1024) line,keyword
mesh_Nnodes = 0_pInt
mesh_Nelems = 0_pInt
rewind(myUnit)
read(myUnit,'(a1024)') line
myPos = IO_stringPos(line,2_pInt)
keyword = IO_lc(IO_StringValue(line,myPos,2_pInt))
if (keyword(1:4) == 'head') then
headerLength = IO_intValue(line,myPos,1_pInt) + 1_pInt
else
call IO_error(error_ID=42_pInt)
endif
rewind(myUnit)
do i = 1, headerLength
read(myUnit,'(a1024)') line
myPos = IO_stringPos(line,maxNchunks)
if ( IO_lc(IO_StringValue(line,myPos,1_pInt)) == 'resolution') then
do j = 2,6,2
select case (IO_lc(IO_stringValue(line,myPos,j)))
case('a')
a = IO_intValue(line,myPos,j+1_pInt)
case('b')
b = IO_intValue(line,myPos,j+1_pInt)
case('c')
c = IO_intValue(line,myPos,j+1_pInt)
end select
enddo
mesh_Nelems = a * b * c
mesh_Nnodes = (1 + a)*(1 + b)*(1 + c)
endif
enddo
endsubroutine
!********************************************************************
! count overall number of nodes and elements in mesh
!
! mesh_Nelems, mesh_Nnodes
!********************************************************************
subroutine mesh_marc_count_nodesAndElements (myUnit)
use prec, only: pInt
use IO
implicit none
integer(pInt), parameter :: maxNchunks = 4
integer(pInt), dimension (1+2*maxNchunks) :: myPos
integer(pInt) myUnit
character(len=300) line
mesh_Nnodes = 0_pInt
mesh_Nelems = 0_pInt
610 FORMAT(A300)
rewind(myUnit)
do
read (myUnit,610,END=620) line
myPos = IO_stringPos(line,maxNchunks)
if ( IO_lc(IO_StringValue(line,myPos,1_pInt)) == 'sizing') then
mesh_Nelems = IO_IntValue (line,myPos,3_pInt)
mesh_Nnodes = IO_IntValue (line,myPos,4_pInt)
exit
endif
enddo
620 endsubroutine
!********************************************************************
! count overall number of nodes and elements in mesh
!
! mesh_Nelems, mesh_Nnodes
!********************************************************************
subroutine mesh_abaqus_count_nodesAndElements (myUnit)
use prec, only: pInt
use IO
implicit none
integer(pInt), parameter :: maxNchunks = 2
integer(pInt), dimension (1+2*maxNchunks) :: myPos
character(len=300) line
integer(pInt) myUnit
logical inPart
mesh_Nnodes = 0
mesh_Nelems = 0
610 FORMAT(A300)
inPart = .false.
rewind(myUnit)
do
read (myUnit,610,END=620) line
myPos = IO_stringPos(line,maxNchunks)
if ( IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*part' ) inPart = .true.
if ( IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*end' .and. &
IO_lc(IO_stringValue(line,myPos,2_pInt)) == 'part' ) inPart = .false.
if (inPart .or. noPart) then
select case ( IO_lc(IO_stringValue(line,myPos,1_pInt)))
case('*node')
if( &
IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'output' .and. &
IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'print' .and. &
IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'file' .and. &
IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'response' &
) &
mesh_Nnodes = mesh_Nnodes + IO_countDataLines(myUnit)
case('*element')
if( &
IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'output' .and. &
IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'matrix' .and. &
IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'response' &
) then
mesh_Nelems = mesh_Nelems + IO_countDataLines(myUnit)
endif
endselect
endif
enddo
620 if (mesh_Nnodes < 2) call IO_error(error_ID=900_pInt)
if (mesh_Nelems == 0) call IO_error(error_ID=901_pInt)
endsubroutine
!********************************************************************
! count overall number of element sets in mesh
!
! mesh_NelemSets, mesh_maxNelemInSet
!********************************************************************
subroutine mesh_marc_count_elementSets (myUnit)
use prec, only: pInt
use IO
implicit none
integer(pInt), parameter :: maxNchunks = 2
integer(pInt), dimension (1+2*maxNchunks) :: myPos
integer(pInt) myUnit
character(len=300) line
mesh_NelemSets = 0_pInt
mesh_maxNelemInSet = 0_pInt
610 FORMAT(A300)
rewind(myUnit)
do
read (myUnit,610,END=620) line
myPos = IO_stringPos(line,maxNchunks)
if ( IO_lc(IO_StringValue(line,myPos,1_pInt)) == 'define' .and. &
IO_lc(IO_StringValue(line,myPos,2_pInt)) == 'element' ) then
mesh_NelemSets = mesh_NelemSets + 1_pInt
mesh_maxNelemInSet = max(mesh_maxNelemInSet, &
IO_countContinousIntValues(myUnit))
endif
enddo
620 endsubroutine
!********************************************************************
! count overall number of element sets in mesh
!
! mesh_NelemSets, mesh_maxNelemInSet
!********************************************************************
subroutine mesh_abaqus_count_elementSets (myUnit)
use prec, only: pInt
use IO
implicit none
integer(pInt), parameter :: maxNchunks = 2
integer(pInt), dimension (1+2*maxNchunks) :: myPos
character(len=300) line
integer(pInt) myUnit
logical inPart
mesh_NelemSets = 0_pInt
mesh_maxNelemInSet = mesh_Nelems ! have to be conservative, since Abaqus allows for recursive definitons
610 FORMAT(A300)
inPart = .false.
rewind(myUnit)
do
read (myUnit,610,END=620) line
myPos = IO_stringPos(line,maxNchunks)
if ( IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*part' ) inPart = .true.
if ( IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*end' .and. &
IO_lc(IO_stringValue(line,myPos,2_pInt)) == 'part' ) inPart = .false.
if ( (inPart .or. noPart) .and. IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*elset' ) &
mesh_NelemSets = mesh_NelemSets + 1_pInt
enddo
620 continue
if (mesh_NelemSets == 0) call IO_error(error_ID=902_pInt)
endsubroutine
!********************************************************************
! count overall number of solid sections sets in mesh (Abaqus only)
!
! mesh_Nmaterials
!********************************************************************
subroutine mesh_abaqus_count_materials (myUnit)
use prec, only: pInt
use IO
implicit none
integer(pInt), parameter :: maxNchunks = 2
integer(pInt), dimension (1+2*maxNchunks) :: myPos
character(len=300) line
integer(pInt) myUnit
logical inPart
mesh_Nmaterials = 0_pInt
610 FORMAT(A300)
inPart = .false.
rewind(myUnit)
do
read (myUnit,610,END=620) line
myPos = IO_stringPos(line,maxNchunks)
if ( IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*part' ) inPart = .true.
if ( IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*end' .and. &
IO_lc(IO_stringValue(line,myPos,2_pInt)) == 'part' ) inPart = .false.
if ( (inPart .or. noPart) .and. &
IO_lc(IO_StringValue(line,myPos,1_pInt)) == '*solid' .and. &
IO_lc(IO_StringValue(line,myPos,2_pInt)) == 'section' ) &
mesh_Nmaterials = mesh_Nmaterials + 1_pInt
enddo
620 if (mesh_Nmaterials == 0) call IO_error(error_ID=903_pInt)
endsubroutine
!********************************************************************
! count overall number of cpElements in mesh
!
! mesh_NcpElems
!********************************************************************
subroutine mesh_spectral_count_cpElements ()
use prec, only: pInt
implicit none
mesh_NcpElems = mesh_Nelems
endsubroutine
!********************************************************************
! count overall number of cpElements in mesh
!
! mesh_NcpElems
!********************************************************************
subroutine mesh_marc_count_cpElements (myUnit)
use prec, only: pInt
use IO
implicit none
integer(pInt), parameter :: maxNchunks = 1
integer(pInt), dimension (1+2*maxNchunks) :: myPos
integer(pInt) myUnit,i
character(len=300) line
mesh_NcpElems = 0_pInt
610 FORMAT(A300)
rewind(myUnit)
do
read (myUnit,610,END=620) line
myPos = IO_stringPos(line,maxNchunks)
if ( IO_lc(IO_stringValue(line,myPos,1_pInt)) == 'hypoelastic') then
do i=1,3+hypoelasticTableStyle ! Skip 3 or 4 lines
read (myUnit,610,END=620) line
enddo
mesh_NcpElems = mesh_NcpElems + IO_countContinousIntValues(myUnit)
exit
endif
enddo
620 endsubroutine
!********************************************************************
! count overall number of cpElements in mesh
!
! mesh_NcpElems
!********************************************************************
subroutine mesh_abaqus_count_cpElements (myUnit)
use prec, only: pInt
use IO
implicit none
integer(pInt), parameter :: maxNchunks = 2
integer(pInt), dimension (1+2*maxNchunks) :: myPos
character(len=300) line
integer(pInt) myUnit,i,k
logical materialFound
character (len=64) materialName,elemSetName
mesh_NcpElems = 0
610 FORMAT(A300)
rewind(myUnit)
do
read (myUnit,610,END=620) line
myPos = IO_stringPos(line,maxNchunks)
select case ( IO_lc(IO_stringValue(line,myPos,1_pInt)) )
case('*material')
materialName = trim(IO_extractValue(IO_lc(IO_stringValue(line,myPos,2_pInt)),'name')) ! extract name=value
materialFound = materialName /= '' ! valid name?
case('*user')
if (IO_lc(IO_StringValue(line,myPos,2_pInt)) == 'material' .and. materialFound) then
do i = 1,mesh_Nmaterials ! look thru material names
if (materialName == mesh_nameMaterial(i)) then ! found one
elemSetName = mesh_mapMaterial(i) ! take corresponding elemSet
do k = 1,mesh_NelemSets ! look thru all elemSet definitions
if (elemSetName == mesh_nameElemSet(k)) & ! matched?
mesh_NcpElems = mesh_NcpElems + mesh_mapElemSet(1,k) ! add those elem count
enddo
endif
enddo
materialFound = .false.
endif
endselect
enddo
620 if (mesh_NcpElems == 0) call IO_error(error_ID=906_pInt)
endsubroutine
!********************************************************************
! map element sets
!
! allocate globals: mesh_nameElemSet, mesh_mapElemSet
!********************************************************************
subroutine mesh_marc_map_elementSets (myUnit)
use prec, only: pInt
use IO
implicit none
integer(pInt), parameter :: maxNchunks = 4
integer(pInt), dimension (1+2*maxNchunks) :: myPos
character(len=300) line
integer(pInt) myUnit,elemSet
allocate (mesh_nameElemSet(mesh_NelemSets)) ; mesh_nameElemSet = ''
allocate (mesh_mapElemSet(1_pInt+mesh_maxNelemInSet,mesh_NelemSets)) ; mesh_mapElemSet = 0_pInt
610 FORMAT(A300)
elemSet = 0_pInt
rewind(myUnit)
do
read (myUnit,610,END=640) line
myPos = IO_stringPos(line,maxNchunks)
if( (IO_lc(IO_stringValue(line,myPos,1_pInt)) == 'define' ) .and. &
(IO_lc(IO_stringValue(line,myPos,2_pInt)) == 'element' ) ) then
elemSet = elemSet+1
mesh_nameElemSet(elemSet) = trim(IO_stringValue(line,myPos,4_pInt))
mesh_mapElemSet(:,elemSet) = IO_continousIntValues(myUnit,mesh_maxNelemInSet,mesh_nameElemSet,mesh_mapElemSet,mesh_NelemSets)
endif
enddo
640 endsubroutine
!********************************************************************
! Build element set mapping
!
! allocate globals: mesh_nameElemSet, mesh_mapElemSet
!********************************************************************
subroutine mesh_abaqus_map_elementSets (myUnit)
use prec, only: pInt
use IO
implicit none
integer(pInt), parameter :: maxNchunks = 4_pInt
integer(pInt), dimension (1_pInt+2_pInt*maxNchunks) :: myPos
character(len=300) line
integer(pInt) myUnit,elemSet,i
logical inPart
allocate (mesh_nameElemSet(mesh_NelemSets)) ; mesh_nameElemSet = ''
allocate (mesh_mapElemSet(1_pInt+mesh_maxNelemInSet,mesh_NelemSets)) ; mesh_mapElemSet = 0_pInt
610 FORMAT(A300)
elemSet = 0_pInt
inPart = .false.
rewind(myUnit)
do
read (myUnit,610,END=640) line
myPos = IO_stringPos(line,maxNchunks)
if ( IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*part' ) inPart = .true.
if ( IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*end' .and. &
IO_lc(IO_stringValue(line,myPos,2_pInt)) == 'part' ) inPart = .false.
if ( (inPart .or. noPart) .and. IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*elset' ) then
elemSet = elemSet + 1_pInt
mesh_nameElemSet(elemSet) = trim(IO_extractValue(IO_lc(IO_stringValue(line,myPos,2_pInt)),'elset'))
mesh_mapElemSet(:,elemSet) = IO_continousIntValues(myUnit,mesh_Nelems,mesh_nameElemSet,&
mesh_mapElemSet,elemSet-1_pInt)
endif
enddo
640 do i = 1_pInt,elemSet
! write(6,*)'elemSetName: ',mesh_nameElemSet(i)
! write(6,*)'elems in Elset',mesh_mapElemSet(:,i)
if (mesh_mapElemSet(1,i) == 0_pInt) call IO_error(error_ID=904_pInt,ext_msg=mesh_nameElemSet(i))
enddo
endsubroutine
!********************************************************************
! map solid section (Abaqus only)
!
! allocate globals: mesh_nameMaterial, mesh_mapMaterial
!********************************************************************
subroutine mesh_abaqus_map_materials (myUnit)
use prec, only: pInt
use IO
implicit none
integer(pInt), parameter :: maxNchunks = 20
integer(pInt), dimension (1+2*maxNchunks) :: myPos
character(len=300) line
integer(pInt) myUnit,i,count
logical inPart
character(len=64) elemSetName,materialName
allocate (mesh_nameMaterial(mesh_Nmaterials)) ; mesh_nameMaterial = ''
allocate (mesh_mapMaterial(mesh_Nmaterials)) ; mesh_mapMaterial = ''
610 FORMAT(A300)
count = 0
inPart = .false.
rewind(myUnit)
do
read (myUnit,610,END=620) line
myPos = IO_stringPos(line,maxNchunks)
if ( IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*part' ) inPart = .true.
if ( IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*end' .and. &
IO_lc(IO_stringValue(line,myPos,2_pInt)) == 'part' ) inPart = .false.
if ( (inPart .or. noPart) .and. &
IO_lc(IO_StringValue(line,myPos,1_pInt)) == '*solid' .and. &
IO_lc(IO_StringValue(line,myPos,2_pInt)) == 'section' ) then
elemSetName = ''
materialName = ''
do i = 3,myPos(1_pInt)
if (IO_extractValue(IO_lc(IO_stringValue(line,myPos,i)),'elset') /= '') &
elemSetName = trim(IO_extractValue(IO_lc(IO_stringValue(line,myPos,i)),'elset'))
if (IO_extractValue(IO_lc(IO_stringValue(line,myPos,i)),'material') /= '') &
materialName = trim(IO_extractValue(IO_lc(IO_stringValue(line,myPos,i)),'material'))
enddo
if (elemSetName /= '' .and. materialName /= '') then
count = count + 1_pInt
mesh_nameMaterial(count) = materialName ! name of material used for this section
mesh_mapMaterial(count) = elemSetName ! mapped to respective element set
endif
endif
enddo
620 if (count==0) call IO_error(error_ID=905_pInt)
do i=1,count
! write(6,*)'name of materials: ',i,mesh_nameMaterial(i)
! write(6,*)'name of elemSets: ',i,mesh_mapMaterial(i)
if (mesh_nameMaterial(i)=='' .or. mesh_mapMaterial(i)=='') call IO_error(error_ID=905_pInt)
enddo
endsubroutine
!********************************************************************
! map nodes from FE id to internal (consecutive) representation
!
! allocate globals: mesh_mapFEtoCPnode
!********************************************************************
subroutine mesh_spectral_map_nodes ()
use prec, only: pInt
implicit none
integer(pInt) i
allocate (mesh_mapFEtoCPnode(2_pInt,mesh_Nnodes)) ; mesh_mapFEtoCPnode = 0_pInt
forall (i = 1:mesh_Nnodes) &
mesh_mapFEtoCPnode(1:2,i) = i
endsubroutine
!********************************************************************
! map nodes from FE id to internal (consecutive) representation
!
! allocate globals: mesh_mapFEtoCPnode
!********************************************************************
subroutine mesh_marc_map_nodes (myUnit)
use prec, only: pInt
use math, only: qsort
use IO
implicit none
integer(pInt), parameter :: maxNchunks = 1_pInt
integer(pInt), dimension (1_pInt+2_pInt*maxNchunks) :: myPos
character(len=300) line
integer(pInt), dimension (mesh_Nnodes) :: node_count
integer(pInt) myUnit,i
allocate (mesh_mapFEtoCPnode(2_pInt,mesh_Nnodes)) ; mesh_mapFEtoCPnode = 0_pInt
610 FORMAT(A300)
node_count = 0_pInt
rewind(myUnit)
do
read (myUnit,610,END=650) line
myPos = IO_stringPos(line,maxNchunks)
if( IO_lc(IO_stringValue(line,myPos,1_pInt)) == 'coordinates' ) then
read (myUnit,610,END=650) line ! skip crap line
do i = 1,mesh_Nnodes
read (myUnit,610,END=650) line
mesh_mapFEtoCPnode(1_pInt,i) = IO_fixedIntValue (line,[ 0_pInt,10_pInt],1_pInt)
mesh_mapFEtoCPnode(2_pInt,i) = i
enddo
exit
endif
enddo
650 call qsort(mesh_mapFEtoCPnode,1_pInt,int(size(mesh_mapFEtoCPnode,2_pInt),pInt))
endsubroutine
!********************************************************************
! map nodes from FE id to internal (consecutive) representation
!
! allocate globals: mesh_mapFEtoCPnode
!********************************************************************
subroutine mesh_abaqus_map_nodes (myUnit)
use prec, only: pInt
use math, only: qsort
use IO
implicit none
integer(pInt), parameter :: maxNchunks = 2
integer(pInt), dimension (1+2*maxNchunks) :: myPos
character(len=300) line
integer(pInt) myUnit,i,count,cpNode
logical inPart
allocate (mesh_mapFEtoCPnode(2_pInt,mesh_Nnodes)) ; mesh_mapFEtoCPnode = 0_pInt
610 FORMAT(A300)
cpNode = 0_pInt
inPart = .false.
rewind(myUnit)
do
read (myUnit,610,END=650) line
myPos = IO_stringPos(line,maxNchunks)
if ( IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*part' ) inPart = .true.
if ( IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*end' .and. &
IO_lc(IO_stringValue(line,myPos,2_pInt)) == 'part' ) inPart = .false.
if( (inPart .or. noPart) .and. &
IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*node' .and. &
( IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'output' .and. &
IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'print' .and. &
IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'file' .and. &
IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'response' ) &
) then
count = IO_countDataLines(myUnit)
do i = 1,count
backspace(myUnit)
enddo
do i = 1,count
read (myUnit,610,END=650) line
myPos = IO_stringPos(line,maxNchunks)
cpNode = cpNode + 1_pInt
mesh_mapFEtoCPnode(1_pInt,cpNode) = IO_intValue(line,myPos,1_pInt)
mesh_mapFEtoCPnode(2_pInt,cpNode) = cpNode
enddo
endif
enddo
650 call qsort(mesh_mapFEtoCPnode,1_pInt,int(size(mesh_mapFEtoCPnode,2_pInt),pInt))
if (int(size(mesh_mapFEtoCPnode),pInt) == 0_pInt) call IO_error(error_ID=908_pInt)
endsubroutine
!********************************************************************
! map elements from FE id to internal (consecutive) representation
!
! allocate globals: mesh_mapFEtoCPelem
!********************************************************************
subroutine mesh_spectral_map_elements ()
use prec, only: pInt
implicit none
integer(pInt) i
allocate (mesh_mapFEtoCPelem(2_pInt,mesh_NcpElems)) ; mesh_mapFEtoCPelem = 0_pInt
forall (i = 1:mesh_NcpElems) &
mesh_mapFEtoCPelem(1:2,i) = i
endsubroutine
!********************************************************************
! map elements from FE id to internal (consecutive) representation
!
! allocate globals: mesh_mapFEtoCPelem
!********************************************************************
subroutine mesh_marc_map_elements (myUnit)
use prec, only: pInt
use math, only: qsort
use IO
implicit none
integer(pInt), parameter :: maxNchunks = 1
integer(pInt), dimension (1+2*maxNchunks) :: myPos
character(len=300) line
integer(pInt), dimension (1+mesh_NcpElems) :: contInts
integer(pInt) myUnit,i,cpElem
allocate (mesh_mapFEtoCPelem(2,mesh_NcpElems)) ; mesh_mapFEtoCPelem = 0_pInt
610 FORMAT(A300)
cpElem = 0_pInt
rewind(myUnit)
do
read (myUnit,610,END=660) line
myPos = IO_stringPos(line,maxNchunks)
if( IO_lc(IO_stringValue(line,myPos,1_pInt)) == 'hypoelastic' ) then
do i=1,3+hypoelasticTableStyle ! skip three (or four if new table style!) lines
read (myUnit,610,END=660) line
enddo
contInts = IO_continousIntValues(myUnit,mesh_NcpElems,mesh_nameElemSet,mesh_mapElemSet,mesh_NelemSets)
do i = 1,contInts(1)
cpElem = cpElem+1
mesh_mapFEtoCPelem(1,cpElem) = contInts(1+i)
mesh_mapFEtoCPelem(2,cpElem) = cpElem
enddo
endif
enddo
660 call qsort(mesh_mapFEtoCPelem,1_pInt,int(size(mesh_mapFEtoCPelem,2_pInt),pInt)) ! should be mesh_NcpElems
endsubroutine
!********************************************************************
! map elements from FE id to internal (consecutive) representation
!
! allocate globals: mesh_mapFEtoCPelem
!********************************************************************
subroutine mesh_abaqus_map_elements (myUnit)
use prec, only: pInt
use math, only: qsort
use IO
implicit none
integer(pInt), parameter :: maxNchunks = 2
integer(pInt), dimension (1+2*maxNchunks) :: myPos
character(len=300) line
integer(pInt) myUnit,i,j,k,cpElem
logical materialFound
character (len=64) materialName,elemSetName ! why limited to 64? ABAQUS?
allocate (mesh_mapFEtoCPelem(2,mesh_NcpElems)) ; mesh_mapFEtoCPelem = 0_pInt
610 FORMAT(A300)
cpElem = 0_pInt
rewind(myUnit)
do
read (myUnit,610,END=660) line
myPos = IO_stringPos(line,maxNchunks)
select case ( IO_lc(IO_stringValue(line,myPos,1_pInt)) )
case('*material')
materialName = trim(IO_extractValue(IO_lc(IO_stringValue(line,myPos,2_pInt)),'name')) ! extract name=value
materialFound = materialName /= '' ! valid name?
case('*user')
if (IO_lc(IO_stringValue(line,myPos,2_pInt)) == 'material' .and. materialFound) then
do i = 1,mesh_Nmaterials ! look thru material names
if (materialName == mesh_nameMaterial(i)) then ! found one
elemSetName = mesh_mapMaterial(i) ! take corresponding elemSet
do k = 1,mesh_NelemSets ! look thru all elemSet definitions
if (elemSetName == mesh_nameElemSet(k)) then ! matched?
do j = 1,mesh_mapElemSet(1,k)
cpElem = cpElem + 1_pInt
mesh_mapFEtoCPelem(1,cpElem) = mesh_mapElemSet(1+j,k) ! store FE id
mesh_mapFEtoCPelem(2,cpElem) = cpElem ! store our id
enddo
endif
enddo
endif
enddo
materialFound = .false.
endif
endselect
enddo
660 call qsort(mesh_mapFEtoCPelem,1_pInt,int(size(mesh_mapFEtoCPelem,2_pInt),pInt)) ! should be mesh_NcpElems
if (int(size(mesh_mapFEtoCPelem),pInt) < 2_pInt) call IO_error(error_ID=907_pInt)
endsubroutine
!********************************************************************
! get maximum count of nodes, IPs, IP neighbors, and subNodes
! among cpElements
!
! _maxNnodes, _maxNips, _maxNipNeighbors, _maxNsubNodes
!********************************************************************
subroutine mesh_spectral_count_cpSizes ()
use prec, only: pInt
implicit none
integer(pInt) t
t = FE_mapElemtype('C3D8R') ! fake 3D hexahedral 8 node 1 IP element
mesh_maxNnodes = FE_Nnodes(t)
mesh_maxNips = FE_Nips(t)
mesh_maxNipNeighbors = FE_NipNeighbors(t)
mesh_maxNsubNodes = FE_NsubNodes(t)
endsubroutine
!********************************************************************
! get maximum count of nodes, IPs, IP neighbors, and subNodes
! among cpElements
!
! _maxNnodes, _maxNips, _maxNipNeighbors, _maxNsubNodes
!********************************************************************
subroutine mesh_marc_count_cpSizes (myUnit)
use prec, only: pInt
use IO
implicit none
integer(pInt), parameter :: maxNchunks = 2
integer(pInt), dimension (1+2*maxNchunks) :: myPos
character(len=300) line
integer(pInt) myUnit,i,t,e
mesh_maxNnodes = 0_pInt
mesh_maxNips = 0_pInt
mesh_maxNipNeighbors = 0_pInt
mesh_maxNsubNodes = 0_pInt
610 FORMAT(A300)
rewind(myUnit)
do
read (myUnit,610,END=630) line
myPos = IO_stringPos(line,maxNchunks)
if( IO_lc(IO_stringValue(line,myPos,1_pInt)) == 'connectivity' ) then
read (myUnit,610,END=630) line ! Garbage line
do i=1,mesh_Nelems ! read all elements
read (myUnit,610,END=630) line
myPos = IO_stringPos(line,maxNchunks) ! limit to id and type
e = mesh_FEasCP('elem',IO_intValue(line,myPos,1_pInt))
if (e /= 0) then
t = FE_mapElemtype(IO_stringValue(line,myPos,2_pInt))
mesh_maxNnodes = max(mesh_maxNnodes,FE_Nnodes(t))
mesh_maxNips = max(mesh_maxNips,FE_Nips(t))
mesh_maxNipNeighbors = max(mesh_maxNipNeighbors,FE_NipNeighbors(t))
mesh_maxNsubNodes = max(mesh_maxNsubNodes,FE_NsubNodes(t))
call IO_skipChunks(myUnit,FE_NoriginalNodes(t)-(myPos(1_pInt)-2_pInt)) ! read on if FE_Nnodes exceeds node count present on current line
endif
enddo
exit
endif
enddo
630 endsubroutine
!********************************************************************
! get maximum count of nodes, IPs, IP neighbors, and subNodes
! among cpElements
!
! _maxNnodes, _maxNips, _maxNipNeighbors, _maxNsubNodes
!********************************************************************
subroutine mesh_abaqus_count_cpSizes (myUnit)
use prec, only: pInt
use IO
implicit none
integer(pInt), parameter :: maxNchunks = 2
integer(pInt), dimension (1+2*maxNchunks) :: myPos
character(len=300) line
integer(pInt) myUnit,i,count,t
logical inPart
mesh_maxNnodes = 0_pInt
mesh_maxNips = 0_pInt
mesh_maxNipNeighbors = 0_pInt
mesh_maxNsubNodes = 0_pInt
610 FORMAT(A300)
inPart = .false.
rewind(myUnit)
do
read (myUnit,610,END=620) line
myPos = IO_stringPos(line,maxNchunks)
if ( IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*part' ) inPart = .true.
if ( IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*end' .and. &
IO_lc(IO_stringValue(line,myPos,2_pInt)) == 'part' ) inPart = .false.
if( (inPart .or. noPart) .and. &
IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*element' .and. &
( IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'output' .and. &
IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'matrix' .and. &
IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'response' ) &
) then
t = FE_mapElemtype(IO_extractValue(IO_lc(IO_stringValue(line,myPos,2_pInt)),'type')) ! remember elem type
if (t==0) call IO_error(error_ID=910_pInt,ext_msg='mesh_abaqus_count_cpSizes')
count = IO_countDataLines(myUnit)
do i = 1,count
backspace(myUnit)
enddo
do i = 1,count
read (myUnit,610,END=620) line
myPos = IO_stringPos(line,maxNchunks) ! limit to 64 nodes max
if (mesh_FEasCP('elem',IO_intValue(line,myPos,1_pInt)) /= 0) then ! disregard non CP elems
mesh_maxNnodes = max(mesh_maxNnodes,FE_Nnodes(t))
mesh_maxNips = max(mesh_maxNips,FE_Nips(t))
mesh_maxNipNeighbors = max(mesh_maxNipNeighbors,FE_NipNeighbors(t))
mesh_maxNsubNodes = max(mesh_maxNsubNodes,FE_NsubNodes(t))
endif
enddo
endif
enddo
620 endsubroutine
!********************************************************************
! store x,y,z coordinates of all nodes in mesh
!
! allocate globals:
! _node
!********************************************************************
subroutine mesh_spectral_build_nodes (myUnit)
use prec, only: pInt
use IO
implicit none
integer(pInt), parameter :: maxNchunks = 7
integer(pInt), dimension (1+2*maxNchunks) :: myPos
integer(pInt) a,b,c,n,i,j,headerLength
real(pReal) x,y,z
logical gotResolution,gotDimension
integer(pInt) myUnit
character(len=1024) line, keyword
allocate ( mesh_node0 (3,mesh_Nnodes) ); mesh_node0 = 0.0_pReal
allocate ( mesh_node (3,mesh_Nnodes) ); mesh_node = 0.0_pReal
a = 1_pInt
b = 1_pInt
c = 1_pInt
x = 1.0_pReal
y = 1.0_pReal
z = 1.0_pReal
gotResolution = .false.
gotDimension = .false.
rewind(myUnit)
read(myUnit,'(a1024)') line
myPos = IO_stringPos(line,2_pInt)
keyword = IO_lc(IO_StringValue(line,myPos,2_pInt))
if (keyword(1:4) == 'head') then
headerLength = IO_intValue(line,myPos,1_pInt) + 1_pInt
else
call IO_error(error_ID=42_pInt)
endif
rewind(myUnit)
do i = 1, headerLength
read(myUnit,'(a1024)') line
myPos = IO_stringPos(line,maxNchunks)
select case ( IO_lc(IO_StringValue(line,myPos,1_pInt)) )
case ('dimension')
gotDimension = .true.
do j = 2,6,2
select case (IO_lc(IO_stringValue(line,myPos,j)))
case('x')
x = IO_floatValue(line,myPos,j+1_pInt)
case('y')
y = IO_floatValue(line,myPos,j+1_pInt)
case('z')
z = IO_floatValue(line,myPos,j+1_pInt)
end select
enddo
case ('resolution')
gotResolution = .true.
do j = 2,6,2
select case (IO_lc(IO_stringValue(line,myPos,j)))
case('a')
a = 1_pInt + IO_intValue(line,myPos,j+1_pInt)
case('b')
b = 1_pInt + IO_intValue(line,myPos,j+1_pInt)
case('c')
c = 1_pInt + IO_intValue(line,myPos,j+1_pInt)
end select
enddo
end select
enddo
! --- sanity checks ---
if ((.not. gotDimension) .or. (.not. gotResolution)) call IO_error(error_ID=42_pInt)
if ((a < 1) .or. (b < 1) .or. (c < 0)) call IO_error(error_ID=43_pInt) ! 1_pInt is already added
if ((x <= 0.0_pReal) .or. (y <= 0.0_pReal) .or. (z <= 0.0_pReal)) call IO_error(error_ID=44_pInt)
forall (n = 0:mesh_Nnodes-1)
mesh_node0(1,n+1) = x * dble(mod(n,a)) / dble(a-1_pInt)
mesh_node0(2,n+1) = y * dble(mod(n/a,b)) / dble(b-1_pInt)
mesh_node0(3,n+1) = z * dble(mod(n/a/b,c)) / dble(c-1_pInt)
end forall
mesh_node = mesh_node0 !why?
endsubroutine
!********************************************************************
! store x,y,z coordinates of all nodes in mesh
!
! allocate globals:
! _node
!********************************************************************
subroutine mesh_marc_build_nodes (myUnit)
use prec, only: pInt
use IO
implicit none
integer(pInt), dimension(5), parameter :: node_ends = (/0,10,30,50,70/)
integer(pInt), parameter :: maxNchunks = 1
integer(pInt), dimension (1+2*maxNchunks) :: myPos
character(len=300) line
integer(pInt) myUnit,i,j,m
allocate ( mesh_node0 (3,mesh_Nnodes) ); mesh_node0 = 0.0_pReal
allocate ( mesh_node (3,mesh_Nnodes) ); mesh_node = 0.0_pReal
610 FORMAT(A300)
rewind(myUnit)
do
read (myUnit,610,END=670) line
myPos = IO_stringPos(line,maxNchunks)
if( IO_lc(IO_stringValue(line,myPos,1_pInt)) == 'coordinates' ) then
read (myUnit,610,END=670) line ! skip crap line
do i=1,mesh_Nnodes
read (myUnit,610,END=670) line
m = mesh_FEasCP('node',IO_fixedIntValue(line,node_ends,1_pInt))
forall (j = 1_pInt:3_pInt) mesh_node0(j,m) = IO_fixedNoEFloatValue(line,node_ends,j+1_pInt)
enddo
exit
endif
enddo
670 mesh_node = mesh_node0
endsubroutine
!********************************************************************
! store x,y,z coordinates of all nodes in mesh
!
! allocate globals:
! _node
!********************************************************************
subroutine mesh_abaqus_build_nodes (myUnit)
use prec, only: pInt
use IO
implicit none
integer(pInt), parameter :: maxNchunks = 4
integer(pInt), dimension (1+2*maxNchunks) :: myPos
character(len=300) line
integer(pInt) myUnit,i,j,m,count
logical inPart
allocate ( mesh_node0 (3,mesh_Nnodes) ); mesh_node0 = 0.0_pReal
allocate ( mesh_node (3,mesh_Nnodes) ); mesh_node = 0.0_pReal
610 FORMAT(A300)
inPart = .false.
rewind(myUnit)
do
read (myUnit,610,END=670) line
myPos = IO_stringPos(line,maxNchunks)
if ( IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*part' ) inPart = .true.
if ( IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*end' .and. &
IO_lc(IO_stringValue(line,myPos,2_pInt)) == 'part' ) inPart = .false.
if( (inPart .or. noPart) .and. &
IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*node' .and. &
( IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'output' .and. &
IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'print' .and. &
IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'file' .and. &
IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'response' ) &
) then
count = IO_countDataLines(myUnit) ! how many nodes are defined here?
do i = 1,count
backspace(myUnit) ! rewind to first entry
enddo
do i = 1,count
read (myUnit,610,END=670) line
myPos = IO_stringPos(line,maxNchunks)
m = mesh_FEasCP('node',IO_intValue(line,myPos,1_pInt))
forall (j=1:3) mesh_node0(j,m) = IO_floatValue(line,myPos,j+1_pInt)
enddo
endif
enddo
670 if (int(size(mesh_node0,2_pInt),pInt) /= mesh_Nnodes) call IO_error(error_ID=909_pInt)
mesh_node = mesh_node0
endsubroutine
!********************************************************************
! store FEid, type, mat, tex, and node list per element
!
! allocate globals:
! _element
!********************************************************************
subroutine mesh_spectral_build_elements (myUnit)
use prec, only: pInt
use IO
implicit none
integer(pInt), parameter :: maxNchunks = 7
integer(pInt), dimension (1+2*maxNchunks) :: myPos
integer(pInt) a,b,c,e,i,j,homog,headerLength,maxIntCount
integer(pInt), dimension(:), allocatable :: microstructures
integer(pInt), dimension(1,1) :: dummySet = 0_pInt
integer(pInt) myUnit
character(len=65536) line,keyword
character(len=64), dimension(1) :: dummyName = ''
a = 1_pInt
b = 1_pInt
c = 1_pInt
rewind(myUnit)
read(myUnit,'(a65536)') line
myPos = IO_stringPos(line,2_pInt)
keyword = IO_lc(IO_StringValue(line,myPos,2_pInt))
if (keyword(1:4) == 'head') then
headerLength = IO_intValue(line,myPos,1_pInt) + 1_pInt
else
call IO_error(error_ID=42_pInt)
endif
rewind(myUnit)
do i = 1, headerLength
read(myUnit,'(a65536)') line
myPos = IO_stringPos(line,maxNchunks)
select case ( IO_lc(IO_StringValue(line,myPos,1_pInt)) )
case ('resolution')
do j = 2,6,2
select case (IO_lc(IO_stringValue(line,myPos,j)))
case('a')
a = 1_pInt + IO_intValue(line,myPos,j+1_pInt)
case('b')
b = 1_pInt + IO_intValue(line,myPos,j+1_pInt)
case('c')
c = 1_pInt + IO_intValue(line,myPos,j+1_pInt)
end select
enddo
case ('homogenization')
homog = IO_intValue(line,myPos,2_pInt)
end select
enddo
maxIntCount = 0_pInt
i = 1_pInt
do while (i > 0_pInt)
i = IO_countContinousIntValues(myUnit)
maxIntCount = max(maxIntCount, i)
enddo
rewind (myUnit)
do i=1,headerLength ! skip header
read(myUnit,'(a65536)') line
enddo
100 allocate (mesh_element (4+mesh_maxNnodes,mesh_NcpElems)) ; mesh_element = 0_pInt
allocate (microstructures (1_pInt+maxIntCount)) ; microstructures = 2_pInt
e = 0_pInt
do while (e < mesh_NcpElems .and. microstructures(1) > 0_pInt) ! fill expected number of elements, stop at end of data (or blank line!)
microstructures = IO_continousIntValues(myUnit,maxIntCount,dummyName,dummySet,0_pInt) ! get affected elements
do i = 1,microstructures(1_pInt)
e = e+1 ! valid element entry
mesh_element( 1,e) = e ! FE id
mesh_element( 2,e) = FE_mapElemtype('C3D8R') ! elem type
mesh_element( 3,e) = homog ! homogenization
mesh_element( 4,e) = microstructures(1_pInt+i) ! microstructure
mesh_element( 5,e) = e + (e-1_pInt)/(a-1_pInt) + ((e-1_pInt)/((a-1_pInt)*(b-1_pInt)))*a ! base node
mesh_element( 6,e) = mesh_element(5,e) + 1_pInt
mesh_element( 7,e) = mesh_element(5,e) + a + 1_pInt
mesh_element( 8,e) = mesh_element(5,e) + a
mesh_element( 9,e) = mesh_element(5,e) + a * b ! second floor base node
mesh_element(10,e) = mesh_element(9,e) + 1_pInt
mesh_element(11,e) = mesh_element(9,e) + a + 1_pInt
mesh_element(12,e) = mesh_element(9,e) + a
mesh_maxValStateVar(1) = max(mesh_maxValStateVar(1),mesh_element(3,e)) !needed for statistics
mesh_maxValStateVar(2) = max(mesh_maxValStateVar(2),mesh_element(4,e))
enddo
enddo
110 deallocate(microstructures)
if (e /= mesh_NcpElems) call IO_error(180_pInt,e)
endsubroutine
!********************************************************************
! store FEid, type, mat, tex, and node list per element
!
! allocate globals:
! _element
!********************************************************************
subroutine mesh_marc_build_elements (myUnit)
use prec, only: pInt
use IO
implicit none
integer(pInt), parameter :: maxNchunks = 66
integer(pInt), dimension (1+2*maxNchunks) :: myPos
character(len=300) line
integer(pInt), dimension(1+mesh_NcpElems) :: contInts
integer(pInt) myUnit,i,j,sv,val,e
allocate (mesh_element (4+mesh_maxNnodes,mesh_NcpElems)) ; mesh_element = 0_pInt
610 FORMAT(A300)
rewind(myUnit)
do
read (myUnit,610,END=620) line
myPos(1:1+2*1) = IO_stringPos(line,1_pInt)
if( IO_lc(IO_stringValue(line,myPos,1_pInt)) == 'connectivity' ) then
read (myUnit,610,END=620) line ! Garbage line
do i = 1,mesh_Nelems
read (myUnit,610,END=620) line
myPos = IO_stringPos(line,maxNchunks) ! limit to 64 nodes max (plus ID, type)
e = mesh_FEasCP('elem',IO_intValue(line,myPos,1_pInt))
if (e /= 0) then ! disregard non CP elems
mesh_element(1,e) = IO_IntValue (line,myPos,1_pInt) ! FE id
mesh_element(2,e) = FE_mapElemtype(IO_StringValue(line,myPos,2_pInt)) ! elem type
forall (j = 1:FE_Nnodes(mesh_element(2,e))) &
mesh_element(j+4,e) = IO_IntValue(line,myPos,j+2_pInt) ! copy FE ids of nodes
call IO_skipChunks(myUnit,FE_NoriginalNodes(mesh_element(2_pInt,e))-(myPos(1_pInt)-2_pInt)) ! read on if FE_Nnodes exceeds node count present on current line
endif
enddo
exit
endif
enddo
620 rewind(myUnit) ! just in case "initial state" apears before "connectivity"
read (myUnit,610,END=620) line
do
myPos(1:1+2*2) = IO_stringPos(line,2_pInt)
if( (IO_lc(IO_stringValue(line,myPos,1_pInt)) == 'initial') .and. &
(IO_lc(IO_stringValue(line,myPos,2_pInt)) == 'state') ) then
if (initialcondTableStyle == 2_pInt) read (myUnit,610,END=620) line ! read extra line for new style
read (myUnit,610,END=630) line ! read line with index of state var
myPos(1:1+2*1) = IO_stringPos(line,1_pInt)
sv = IO_IntValue(line,myPos,1_pInt) ! figure state variable index
if( (sv == 2_pInt).or.(sv == 3_pInt) ) then ! only state vars 2 and 3 of interest
read (myUnit,610,END=620) line ! read line with value of state var
myPos(1:1+2*1) = IO_stringPos(line,1_pInt)
do while (scan(IO_stringValue(line,myPos,1_pInt),'+-',back=.true.)>1_pInt) ! is noEfloat value?
val = NINT(IO_fixedNoEFloatValue(line,(/0_pInt,20_pInt/),1_pInt)) ! state var's value
mesh_maxValStateVar(sv-1) = max(val,mesh_maxValStateVar(sv-1_pInt)) ! remember max val of homogenization and microstructure index
if (initialcondTableStyle == 2) then
read (myUnit,610,END=630) line ! read extra line
read (myUnit,610,END=630) line ! read extra line
endif
contInts = IO_continousIntValues(myUnit,mesh_Nelems,mesh_nameElemSet,mesh_mapElemSet,mesh_NelemSets) ! get affected elements
do i = 1,contInts(1)
e = mesh_FEasCP('elem',contInts(1_pInt+i))
mesh_element(1_pInt+sv,e) = val
enddo
if (initialcondTableStyle == 0) read (myUnit,610,END=620) line ! ignore IP range for old table style
read (myUnit,610,END=630) line
myPos(1:1+2*1) = IO_stringPos(line,1_pInt)
enddo
endif
else
read (myUnit,610,END=630) line
endif
enddo
630 endsubroutine
!********************************************************************
! store FEid, type, mat, tex, and node list per element
!
! allocate globals:
! _element
!********************************************************************
subroutine mesh_abaqus_build_elements (myUnit)
use prec, only: pInt
use IO
implicit none
integer(pInt), parameter :: maxNchunks = 65
integer(pInt), dimension (1+2*maxNchunks) :: myPos
integer(pInt) myUnit,i,j,k,count,e,t,homog,micro
logical inPart,materialFound
character (len=64) materialName,elemSetName
character(len=300) line
allocate (mesh_element (4+mesh_maxNnodes,mesh_NcpElems)) ; mesh_element = 0_pInt
610 FORMAT(A300)
inPart = .false.
rewind(myUnit)
do
read (myUnit,610,END=620) line
myPos(1:1+2*2) = IO_stringPos(line,2_pInt)
if ( IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*part' ) inPart = .true.
if ( IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*end' .and. &
IO_lc(IO_stringValue(line,myPos,2_pInt)) == 'part' ) inPart = .false.
if( (inPart .or. noPart) .and. &
IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*element' .and. &
( IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'output' .and. &
IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'matrix' .and. &
IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'response' ) &
) then
t = FE_mapElemtype(IO_extractValue(IO_lc(IO_stringValue(line,myPos,2_pInt)),'type')) ! remember elem type
if (t==0) call IO_error(error_ID=910_pInt,ext_msg='mesh_abaqus_build_elements')
count = IO_countDataLines(myUnit)
do i = 1,count
backspace(myUnit)
enddo
do i = 1,count
read (myUnit,610,END=620) line
myPos = IO_stringPos(line,maxNchunks) ! limit to 64 nodes max
e = mesh_FEasCP('elem',IO_intValue(line,myPos,1_pInt))
if (e /= 0) then ! disregard non CP elems
mesh_element(1,e) = IO_intValue(line,myPos,1_pInt) ! FE id
mesh_element(2,e) = t ! elem type
forall (j=1:FE_Nnodes(t)) &
mesh_element(4_pInt+j,e) = IO_intValue(line,myPos,1_pInt+j) ! copy FE ids of nodes to position 5:
call IO_skipChunks(myUnit,FE_NoriginalNodes(t)-(myPos(1_pInt)-1_pInt)) ! read on (even multiple lines) if FE_NoriginalNodes exceeds required node count
endif
enddo
endif
enddo
620 rewind(myUnit) ! just in case "*material" definitions apear before "*element"
materialFound = .false.
do
read (myUnit,610,END=630) line
myPos = IO_stringPos(line,maxNchunks)
select case ( IO_lc(IO_StringValue(line,myPos,1_pInt)))
case('*material')
materialName = trim(IO_extractValue(IO_lc(IO_StringValue(line,myPos,2_pInt)),'name')) ! extract name=value
materialFound = materialName /= '' ! valid name?
case('*user')
if ( IO_lc(IO_StringValue(line,myPos,2_pInt)) == 'material' .and. &
materialFound ) then
read (myUnit,610,END=630) line ! read homogenization and microstructure
myPos(1:1+2*2) = IO_stringPos(line,2_pInt)
homog = NINT(IO_floatValue(line,myPos,1_pInt))
micro = NINT(IO_floatValue(line,myPos,2_pInt))
do i = 1,mesh_Nmaterials ! look thru material names
if (materialName == mesh_nameMaterial(i)) then ! found one
elemSetName = mesh_mapMaterial(i) ! take corresponding elemSet
do k = 1,mesh_NelemSets ! look thru all elemSet definitions
if (elemSetName == mesh_nameElemSet(k)) then ! matched?
do j = 1,mesh_mapElemSet(1,k)
e = mesh_FEasCP('elem',mesh_mapElemSet(1+j,k))
mesh_element(3,e) = homog ! store homogenization
mesh_element(4,e) = micro ! store microstructure
mesh_maxValStateVar(1) = max(mesh_maxValStateVar(1),homog)
mesh_maxValStateVar(2) = max(mesh_maxValStateVar(2),micro)
enddo
endif
enddo
endif
enddo
materialFound = .false.
endif
endselect
enddo
630 endsubroutine
!********************************************************************
! get maximum count of shared elements among cpElements and
! build list of elements shared by each node in mesh
!
! _maxNsharedElems
! _sharedElem
!********************************************************************
subroutine mesh_build_sharedElems()
use prec, only: pInt
implicit none
integer(pint) e, & ! element index
t, & ! element type
node, & ! CP node index
j, & ! node index per element
dim, & ! dimension index
nodeTwin ! node twin in the specified dimension
integer(pInt), dimension (mesh_Nnodes) :: node_count
integer(pInt), dimension (:), allocatable :: node_seen
allocate(node_seen(maxval(FE_Nnodes)))
node_count = 0_pInt
do e = 1,mesh_NcpElems
t = mesh_element(2,e) ! get element type
node_seen = 0_pInt ! reset node duplicates
do j = 1,FE_Nnodes(t) ! check each node of element
node = mesh_FEasCP('node',mesh_element(4+j,e)) ! translate to internal (consecutive) numbering
if (all(node_seen /= node)) then
node_count(node) = node_count(node) + 1_pInt ! if FE node not yet encountered -> count it
do dim = 1,3 ! check in each dimension...
nodeTwin = mesh_nodeTwins(dim,node)
if (nodeTwin > 0_pInt) & ! if I am a twin of some node...
node_count(nodeTwin) = node_count(nodeTwin) + 1_pInt ! -> count me again for the twin node
enddo
endif
node_seen(j) = node ! remember this node to be counted already
enddo
enddo
mesh_maxNsharedElems = maxval(node_count) ! most shared node
allocate(mesh_sharedElem(1+mesh_maxNsharedElems,mesh_Nnodes))
mesh_sharedElem = 0_pInt
do e = 1,mesh_NcpElems
t = mesh_element(2,e)
node_seen = 0_pInt
do j = 1,FE_Nnodes(t)
node = mesh_FEasCP('node',mesh_element(4+j,e))
if (all(node_seen /= node)) then
mesh_sharedElem(1,node) = mesh_sharedElem(1,node) + 1_pInt ! count for each node the connected elements
mesh_sharedElem(mesh_sharedElem(1,node)+1,node) = e ! store the respective element id
do dim = 1,3 ! check in each dimension...
nodeTwin = mesh_nodeTwins(dim,node)
if (nodeTwin > 0) then ! if i am a twin of some node...
mesh_sharedElem(1,nodeTwin) = mesh_sharedElem(1,nodeTwin) + 1_pInt ! ...count me again for the twin
mesh_sharedElem(mesh_sharedElem(1,nodeTwin)+1,nodeTwin) = e ! store the respective element id
endif
enddo
endif
node_seen(j) = node
enddo
enddo
deallocate(node_seen)
endsubroutine
!***********************************************************
! build up of IP neighborhood
!
! allocate globals
! _ipNeighborhood
!***********************************************************
subroutine mesh_build_ipNeighborhood()
use prec, only: pInt
implicit none
integer(pInt) myElem, & ! my CP element index
myIP, &
myType, & ! my element type
myFace, &
neighbor, & ! neighor index
neighboringIPkey, & ! positive integer indicating the neighboring IP (for intra-element) and negative integer indicating the face towards neighbor (for neighboring element)
candidateIP, &
neighboringType, & ! element type of neighbor
NlinkedNodes, & ! number of linked nodes
twin_of_linkedNode, & ! node twin of a specific linkedNode
NmatchingNodes, & ! number of matching nodes
dir, & ! direction of periodicity
matchingElem, & ! CP elem number of matching element
matchingFace, & ! face ID of matching element
a, anchor
integer(pInt), dimension(FE_maxmaxNnodesAtIP) :: &
linkedNodes, &
matchingNodes
logical checkTwins
allocate(mesh_ipNeighborhood(2,mesh_maxNipNeighbors,mesh_maxNips,mesh_NcpElems))
mesh_ipNeighborhood = 0_pInt
linkedNodes = 0_pInt
do myElem = 1,mesh_NcpElems ! loop over cpElems
myType = mesh_element(2,myElem) ! get elemType
do myIP = 1,FE_Nips(myType) ! loop over IPs of elem
do neighbor = 1,FE_NipNeighbors(myType) ! loop over neighbors of IP
neighboringIPkey = FE_ipNeighbor(neighbor,myIP,myType)
!*** if the key is positive, the neighbor is inside the element
!*** that means, we have already found our neighboring IP
if (neighboringIPkey > 0) then
mesh_ipNeighborhood(1,neighbor,myIP,myElem) = myElem
mesh_ipNeighborhood(2,neighbor,myIP,myElem) = neighboringIPkey
!*** if the key is negative, the neighbor resides in a neighboring element
!*** that means, we have to look through the face indicated by the key and see which element is behind that face
elseif (neighboringIPkey < 0) then ! neighboring element's IP
myFace = -neighboringIPkey
call mesh_faceMatch(myElem, myFace, matchingElem, matchingFace) ! get face and CP elem id of face match
if (matchingElem > 0_pInt) then ! found match?
neighboringType = mesh_element(2,matchingElem)
!*** trivial solution if neighbor has only one IP
if (FE_Nips(neighboringType) == 1_pInt) then
mesh_ipNeighborhood(1,neighbor,myIP,myElem) = matchingElem
mesh_ipNeighborhood(2,neighbor,myIP,myElem) = 1_pInt
cycle
endif
!*** find those nodes which build the link to the neighbor
NlinkedNodes = 0_pInt
linkedNodes = 0_pInt
do a = 1,FE_maxNnodesAtIP(myType) ! figure my anchor nodes on connecting face
anchor = FE_nodesAtIP(a,myIP,myType)
if (anchor /= 0_pInt) then ! valid anchor node
if (any(FE_nodeOnFace(:,myFace,myType) == anchor)) then ! ip anchor sits on face?
NlinkedNodes = NlinkedNodes + 1_pInt
linkedNodes(NlinkedNodes) = &
mesh_FEasCP('node',mesh_element(4+anchor,myElem)) ! CP id of anchor node
else ! something went wrong with the linkage, since not all anchors sit on my face
NlinkedNodes = 0_pInt
linkedNodes = 0_pInt
exit
endif
endif
enddo
!*** loop through the ips of my neighbor
!*** and try to find an ip with matching nodes
!*** also try to match with node twins
checkCandidateIP: do candidateIP = 1,FE_Nips(neighboringType)
NmatchingNodes = 0_pInt
matchingNodes = 0_pInt
do a = 1,FE_maxNnodesAtIP(neighboringType) ! check each anchor node of that ip
anchor = FE_nodesAtIP(a,candidateIP,neighboringType)
if (anchor /= 0_pInt) then ! valid anchor node
if (any(FE_nodeOnFace(:,matchingFace,neighboringType) == anchor)) then ! sits on matching face?
NmatchingNodes = NmatchingNodes + 1_pInt
matchingNodes(NmatchingNodes) = &
mesh_FEasCP('node',mesh_element(4+anchor,matchingElem)) ! CP id of neighbor's anchor node
else ! no matching, because not all nodes sit on the matching face
NmatchingNodes = 0_pInt
matchingNodes = 0_pInt
exit
endif
endif
enddo
if (NmatchingNodes /= NlinkedNodes) & ! this ip has wrong count of anchors on face
cycle checkCandidateIP
!*** check "normal" nodes whether they match or not
checkTwins = .false.
do a = 1,NlinkedNodes
if (all(matchingNodes /= linkedNodes(a))) then ! this linkedNode does not match any matchingNode
checkTwins = .true.
exit ! no need to search further
endif
enddo
!*** if no match found, then also check node twins
if(checkTwins) then
dir = maxloc(abs(mesh_ipAreaNormal(1:3,neighbor,myIP,myElem)),1) ! check for twins only in direction of the surface normal
do a = 1,NlinkedNodes
twin_of_linkedNode = mesh_nodeTwins(dir,linkedNodes(a))
if (twin_of_linkedNode == 0_pInt .or. & ! twin of linkedNode does not exist...
all(matchingNodes /= twin_of_linkedNode)) then ! ... or it does not match any matchingNode
cycle checkCandidateIP ! ... then check next candidateIP
endif
enddo
endif
!*** we found a match !!!
mesh_ipNeighborhood(1,neighbor,myIP,myElem) = matchingElem
mesh_ipNeighborhood(2,neighbor,myIP,myElem) = candidateIP
exit checkCandidateIP
enddo checkCandidateIP
endif ! end of valid external matching
endif ! end of internal/external matching
enddo
enddo
enddo
endsubroutine
!***********************************************************
! assignment of coordinates for subnodes in each cp element
!
! allocate globals
! _subNodeCoord
!***********************************************************
subroutine mesh_build_subNodeCoords()
use prec, only: pInt,pReal
implicit none
integer(pInt) e,t,n,p
if (.not. allocated(mesh_subNodeCoord)) then
allocate(mesh_subNodeCoord(3,mesh_maxNnodes+mesh_maxNsubNodes,mesh_NcpElems))
endif
mesh_subNodeCoord = 0.0_pReal
do e = 1,mesh_NcpElems ! loop over cpElems
t = mesh_element(2,e) ! get elemType
do n = 1,FE_Nnodes(t)
mesh_subNodeCoord(1:3,n,e) = mesh_node(1:3,mesh_FEasCP('node',mesh_element(4+n,e))) ! loop over nodes of this element type
enddo
do n = 1,FE_NsubNodes(t) ! now for the true subnodes
do p = 1,FE_Nips(t) ! loop through possible parent nodes
if (FE_subNodeParent(p,n,t) > 0) & ! valid parent node
mesh_subNodeCoord(1:3,FE_Nnodes(t)+n,e) = mesh_subNodeCoord(1:3,FE_Nnodes(t)+n,e) &
+ mesh_node(1:3,mesh_FEasCP('node',mesh_element(4+FE_subNodeParent(p,n,t),e))) ! add up parents
enddo
mesh_subNodeCoord(1:3,n+FE_Nnodes(t),e) = mesh_subNodeCoord(1:3,n+FE_Nnodes(t),e) &
/real(count(FE_subNodeParent(:,n,t) > 0),pReal)
enddo
enddo
endsubroutine
!***********************************************************
! calculation of IP coordinates
!
! allocate globals
! _ipCenterOfGravity
!***********************************************************
subroutine mesh_build_ipCoordinates()
use prec, only: pInt, tol_gravityNodePos
implicit none
integer(pInt) e,f,t,i,j,k,n
logical(pInt), dimension(mesh_maxNnodes+mesh_maxNsubNodes) :: gravityNode ! flagList to find subnodes determining center of grav
real(pReal), dimension(3,mesh_maxNnodes+mesh_maxNsubNodes) :: gravityNodePos ! coordinates of subnodes determining center of grav
real(pReal), dimension(3) :: centerOfGravity
if (.not. allocated(mesh_ipCenterOfGravity)) then
allocate(mesh_ipCenterOfGravity(3,mesh_maxNips,mesh_NcpElems))
endif
do e = 1,mesh_NcpElems ! loop over cpElems
t = mesh_element(2,e) ! get elemType
do i = 1,FE_Nips(t) ! loop over IPs of elem
gravityNode = .false. ! reset flagList
gravityNodePos = 0.0_pReal ! reset coordinates
do f = 1,FE_NipNeighbors(t) ! loop over interfaces of IP
do n = 1,FE_NipFaceNodes ! loop over nodes on interface
gravityNode(FE_subNodeOnIPFace(n,f,i,t)) = .true.
gravityNodePos(:,FE_subNodeOnIPFace(n,f,i,t)) = mesh_subNodeCoord(:,FE_subNodeOnIPFace(n,f,i,t),e)
enddo
enddo
do j = 1,mesh_maxNnodes+mesh_maxNsubNodes-1 ! walk through entire flagList except last
if (gravityNode(j)) then ! valid node index
do k = j+1,mesh_maxNnodes+mesh_maxNsubNodes ! walk through remainder of list
if (gravityNode(k) .and. all(abs(gravityNodePos(:,j) - gravityNodePos(:,k)) < tol_gravityNodePos)) then ! found duplicate
gravityNode(j) = .false. ! delete first instance
gravityNodePos(:,j) = 0.0_pReal
exit ! continue with next suspect
endif
enddo
endif
enddo
centerOfGravity = sum(gravityNodePos,2)/real(count(gravityNode),pReal)
mesh_ipCenterOfGravity(:,i,e) = centerOfGravity
enddo
enddo
endsubroutine
!***********************************************************
! calculation of IP volume
!
! allocate globals
! _ipVolume
!***********************************************************
subroutine mesh_build_ipVolumes()
use prec, only: pInt, tol_gravityNodePos
use math, only: math_volTetrahedron
implicit none
integer(pInt) e,f,t,i,j,n
integer(pInt), parameter :: Ntriangles = FE_NipFaceNodes-2 ! each interface is made up of this many triangles
real(pReal), dimension(3,FE_NipFaceNodes) :: nPos ! coordinates of nodes on IP face
real(pReal), dimension(Ntriangles,FE_NipFaceNodes) :: volume ! volumes of possible tetrahedra
if (.not. allocated(mesh_ipVolume)) then
allocate(mesh_ipVolume(mesh_maxNips,mesh_NcpElems))
endif
mesh_ipVolume = 0.0_pReal
do e = 1,mesh_NcpElems ! loop over cpElems
t = mesh_element(2_pInt,e) ! get elemType
do i = 1,FE_Nips(t) ! loop over IPs of elem
do f = 1,FE_NipNeighbors(t) ! loop over interfaces of IP and add tetrahedra which connect to CoG
forall (n = 1:FE_NipFaceNodes) &
nPos(:,n) = mesh_subNodeCoord(:,FE_subNodeOnIPFace(n,f,i,t),e)
forall (n = 1:FE_NipFaceNodes, j = 1:Ntriangles) & ! start at each interface node and build valid triangles to cover interface
volume(j,n) = math_volTetrahedron(nPos(:,n), & ! calc volume of respective tetrahedron to CoG
nPos(:,1_pInt+mod(n-1_pInt +j ,FE_NipFaceNodes)), & ! start at offset j
nPos(:,1+mod(n-1_pInt +j+1_pInt,FE_NipFaceNodes)), & ! and take j's neighbor
mesh_ipCenterOfGravity(:,i,e))
mesh_ipVolume(i,e) = mesh_ipVolume(i,e) + sum(volume) ! add contribution from this interface
enddo
mesh_ipVolume(i,e) = mesh_ipVolume(i,e) / FE_NipFaceNodes ! renormalize with interfaceNodeNum due to loop over them
enddo
enddo
endsubroutine
!***********************************************************
! calculation of IP interface areas
!
! allocate globals
! _ipArea, _ipAreaNormal
!***********************************************************
subroutine mesh_build_ipAreas()
use prec, only: pInt,pReal
use math
implicit none
integer(pInt) e,f,t,i,j,n
integer(pInt), parameter :: Ntriangles = FE_NipFaceNodes-2 ! each interface is made up of this many triangles
real(pReal), dimension (3,FE_NipFaceNodes) :: nPos ! coordinates of nodes on IP face
real(pReal), dimension(3,Ntriangles,FE_NipFaceNodes) :: normal
real(pReal), dimension(Ntriangles,FE_NipFaceNodes) :: area
allocate(mesh_ipArea(mesh_maxNipNeighbors,mesh_maxNips,mesh_NcpElems)) ; mesh_ipArea = 0.0_pReal
allocate(mesh_ipAreaNormal(3_pInt,mesh_maxNipNeighbors,mesh_maxNips,mesh_NcpElems)) ; mesh_ipAreaNormal = 0.0_pReal
do e = 1,mesh_NcpElems ! loop over cpElems
t = mesh_element(2,e) ! get elemType
do i = 1,FE_Nips(t) ! loop over IPs of elem
do f = 1,FE_NipNeighbors(t) ! loop over interfaces of IP
forall (n = 1:FE_NipFaceNodes) nPos(:,n) = mesh_subNodeCoord(:,FE_subNodeOnIPFace(n,f,i,t),e)
forall (n = 1:FE_NipFaceNodes, j = 1:Ntriangles) ! start at each interface node and build valid triangles to cover interface
normal(:,j,n) = math_vectorproduct(nPos(:,1_pInt+mod(n+j-1_pInt,FE_NipFaceNodes)) - nPos(:,n), & ! calc their normal vectors
nPos(:,1_pInt+mod(n+j-0_pInt,FE_NipFaceNodes)) - nPos(:,n))
area(j,n) = sqrt(sum(normal(:,j,n)*normal(:,j,n))) ! and area
end forall
forall (n = 1:FE_NipFaceNodes, j = 1:Ntriangles, area(j,n) > 0.0_pReal) &
normal(1:3,j,n) = normal(1:3,j,n) / area(j,n) ! make myUnit normal
mesh_ipArea(f,i,e) = sum(area) / (FE_NipFaceNodes*2.0_pReal) ! area of parallelograms instead of triangles
mesh_ipAreaNormal(:,f,i,e) = sum(sum(normal,3),2)/& ! average of all valid normals
real(count(area > 0.0_pReal),pReal)
enddo
enddo
enddo
endsubroutine
!***********************************************************
! assignment of twin nodes for each cp node
!
! allocate globals
! _nodeTwins
!***********************************************************
subroutine mesh_build_nodeTwins()
use prec, only: pInt, pReal
implicit none
integer(pInt) dir, & ! direction of periodicity
node, &
minimumNode, &
maximumNode, &
n1, &
n2
integer(pInt), dimension(mesh_Nnodes+1) :: minimumNodes, maximumNodes ! list of surface nodes (minimum and maximum coordinate value) with first entry giving the number of nodes
real(pReal) minCoord, maxCoord, & ! extreme positions in one dimension
tolerance ! tolerance below which positions are assumed identical
real(pReal), dimension(3) :: distance ! distance between two nodes in all three coordinates
logical, dimension(mesh_Nnodes) :: unpaired
allocate(mesh_nodeTwins(3,mesh_Nnodes))
mesh_nodeTwins = 0_pInt
tolerance = 0.001_pReal * minval(mesh_ipVolume) ** 0.333_pReal
do dir = 1,3 ! check periodicity in directions of x,y,z
if (mesh_periodicSurface(dir)) then ! only if periodicity is requested
!*** find out which nodes sit on the surface
!*** and have a minimum or maximum position in this dimension
minimumNodes = 0_pInt
maximumNodes = 0_pInt
minCoord = minval(mesh_node0(dir,:))
maxCoord = maxval(mesh_node0(dir,:))
do node = 1,mesh_Nnodes ! loop through all nodes and find surface nodes
if (abs(mesh_node0(dir,node) - minCoord) <= tolerance) then
minimumNodes(1) = minimumNodes(1) + 1_pInt
minimumNodes(minimumNodes(1)+1) = node
elseif (abs(mesh_node0(dir,node) - maxCoord) <= tolerance) then
maximumNodes(1) = maximumNodes(1) + 1_pInt
maximumNodes(maximumNodes(1)+1) = node
endif
enddo
!*** find the corresponding node on the other side with the same position in this dimension
unpaired = .true.
do n1 = 1,minimumNodes(1)
minimumNode = minimumNodes(n1+1)
if (unpaired(minimumNode)) then
do n2 = 1,maximumNodes(1)
maximumNode = maximumNodes(n2+1)
distance = abs(mesh_node0(:,minimumNode) - mesh_node0(:,maximumNode))
if (sum(distance) - distance(dir) <= tolerance) then ! minimum possible distance (within tolerance)
mesh_nodeTwins(dir,minimumNode) = maximumNode
mesh_nodeTwins(dir,maximumNode) = minimumNode
unpaired(maximumNode) = .false. ! remember this node, we don't have to look for his partner again
exit
endif
enddo
endif
enddo
endif
enddo
endsubroutine
!***********************************************************
! write statistics regarding input file parsing
! to the output file
!
!***********************************************************
subroutine mesh_tell_statistics()
use prec, only: pInt
use math, only: math_range
use IO, only: IO_error
use debug, only: debug_verbosity, &
debug_e, &
debug_i, &
debug_selectiveDebugger
implicit none
integer(pInt), dimension (:,:), allocatable :: mesh_HomogMicro
character(len=64) fmt
integer(pInt) i,e,n,f,t
if (mesh_maxValStateVar(1) < 1_pInt) call IO_error(error_ID=110_pInt) ! no homogenization specified
if (mesh_maxValStateVar(2) < 1_pInt) call IO_error(error_ID=120_pInt) ! no microstructure specified
allocate (mesh_HomogMicro(mesh_maxValStateVar(1),mesh_maxValStateVar(2))); mesh_HomogMicro = 0_pInt
do e = 1,mesh_NcpElems
if (mesh_element(3,e) < 1_pInt) call IO_error(error_ID=110_pInt,e=e) ! no homogenization specified
if (mesh_element(4,e) < 1_pInt) call IO_error(error_ID=120_pInt,e=e) ! no microstructure specified
mesh_HomogMicro(mesh_element(3,e),mesh_element(4,e)) = &
mesh_HomogMicro(mesh_element(3,e),mesh_element(4,e)) + 1 ! count combinations of homogenization and microstructure
enddo
if (debug_verbosity > 0) then
!$OMP CRITICAL (write2out)
write (6,*)
write (6,*) 'Input Parser: STATISTICS'
write (6,*)
write (6,*) mesh_Nelems, ' : total number of elements in mesh'
write (6,*) mesh_NcpElems, ' : total number of CP elements in mesh'
write (6,*) mesh_Nnodes, ' : total number of nodes in mesh'
write (6,*) mesh_maxNnodes, ' : max number of nodes in any CP element'
write (6,*) mesh_maxNips, ' : max number of IPs in any CP element'
write (6,*) mesh_maxNipNeighbors, ' : max number of IP neighbors in any CP element'
write (6,*) mesh_maxNsubNodes, ' : max number of (additional) subnodes in any CP element'
write (6,*) mesh_maxNsharedElems, ' : max number of CP elements sharing a node'
write (6,*)
write (6,*) 'Input Parser: HOMOGENIZATION/MICROSTRUCTURE'
write (6,*)
write (6,*) mesh_maxValStateVar(1), ' : maximum homogenization index'
write (6,*) mesh_maxValStateVar(2), ' : maximum microstructure index'
write (6,*)
write (fmt,'(a,i32.32,a)') '(9x,a2,1x,',mesh_maxValStateVar(2),'(i8))'
write (6,fmt) '+-',math_range(mesh_maxValStateVar(2))
write (fmt,'(a,i32.32,a)') '(i8,1x,a2,1x,',mesh_maxValStateVar(2),'(i8))'
do i=1,mesh_maxValStateVar(1) ! loop over all (possibly assigned) homogenizations
write (6,fmt) i,'| ',mesh_HomogMicro(i,:) ! loop over all (possibly assigned) microstructures
enddo
write(6,*)
write(6,*) 'Input Parser: ADDITIONAL MPIE OPTIONS'
write(6,*)
write(6,*) 'periodic surface : ', mesh_periodicSurface
write(6,*)
call flush(6)
!$OMP END CRITICAL (write2out)
endif
if (debug_verbosity > 1) then
!$OMP CRITICAL (write2out)
write (6,*)
write (6,*) 'Input Parser: SUBNODE COORDINATES'
write (6,*)
write(6,'(a8,1x,a5,1x,a15,1x,a15,1x,a20,3(1x,a12))') 'elem','IP','IP neighbor','IPFaceNodes','subNodeOnIPFace','x','y','z'
do e = 1,mesh_NcpElems ! loop over cpElems
if (debug_selectiveDebugger .and. debug_e /= e) cycle
t = mesh_element(2,e) ! get elemType
do i = 1,FE_Nips(t) ! loop over IPs of elem
if (debug_selectiveDebugger .and. debug_i /= i) cycle
do f = 1,FE_NipNeighbors(t) ! loop over interfaces of IP
do n = 1,FE_NipFaceNodes ! loop over nodes on interface
write(6,'(i8,1x,i5,1x,i15,1x,i15,1x,i20,3(1x,f12.8))') e,i,f,n,FE_subNodeOnIPFace(n,f,i,t),&
mesh_subNodeCoord(1,FE_subNodeOnIPFace(n,f,i,t),e),&
mesh_subNodeCoord(2,FE_subNodeOnIPFace(n,f,i,t),e),&
mesh_subNodeCoord(3,FE_subNodeOnIPFace(n,f,i,t),e)
enddo
enddo
enddo
enddo
write(6,*)
write(6,*) 'Input Parser: IP COORDINATES'
write(6,'(a8,1x,a5,3(1x,a12))') 'elem','IP','x','y','z'
do e = 1,mesh_NcpElems
if (debug_selectiveDebugger .and. debug_e /= e) cycle
do i = 1,FE_Nips(mesh_element(2,e))
if (debug_selectiveDebugger .and. debug_i /= i) cycle
write (6,'(i8,1x,i5,3(1x,f12.8))') e, i, mesh_ipCenterOfGravity(:,i,e)
enddo
enddo
write (6,*)
write (6,*) 'Input Parser: ELEMENT VOLUME'
write (6,*)
write (6,'(a13,1x,e15.8)') 'total volume', sum(mesh_ipVolume)
write (6,*)
write (6,'(a8,1x,a5,1x,a15,1x,a5,1x,a15,1x,a16)') 'elem','IP','volume','face','area','-- normal --'
do e = 1,mesh_NcpElems
if (debug_selectiveDebugger .and. debug_e /= e) cycle
do i = 1,FE_Nips(mesh_element(2,e))
if (debug_selectiveDebugger .and. debug_i /= i) cycle
write (6,'(i8,1x,i5,1x,e15.8)') e,i,mesh_IPvolume(i,e)
do f = 1,FE_NipNeighbors(mesh_element(2,e))
write (6,'(i33,1x,e15.8,1x,3(f6.3,1x))') f,mesh_ipArea(f,i,e),mesh_ipAreaNormal(:,f,i,e)
enddo
enddo
enddo
write (6,*)
write (6,*) 'Input Parser: NODE TWINS'
write (6,*)
write(6,'(a6,3(3x,a6))') ' node','twin_x','twin_y','twin_z'
do n = 1,mesh_Nnodes ! loop over cpNodes
if (debug_e <= mesh_NcpElems) then
if (any(mesh_element(5:,debug_e) == n)) then
write(6,'(i6,3(3x,i6))') n, mesh_nodeTwins(1:3,n)
endif
endif
enddo
write(6,*)
write(6,*) 'Input Parser: IP NEIGHBORHOOD'
write(6,*)
write(6,'(a8,1x,a10,1x,a10,1x,a3,1x,a13,1x,a13)') 'elem','IP','neighbor','','elemNeighbor','ipNeighbor'
do e = 1,mesh_NcpElems ! loop over cpElems
if (debug_selectiveDebugger .and. debug_e /= e) cycle
t = mesh_element(2,e) ! get elemType
do i = 1,FE_Nips(t) ! loop over IPs of elem
if (debug_selectiveDebugger .and. debug_i /= i) cycle
do n = 1,FE_NipNeighbors(t) ! loop over neighbors of IP
write (6,'(i8,1x,i10,1x,i10,1x,a3,1x,i13,1x,i13)') e,i,n,'-->',mesh_ipNeighborhood(1,n,i,e),mesh_ipNeighborhood(2,n,i,e)
enddo
enddo
enddo
!$OMP END CRITICAL (write2out)
endif
deallocate(mesh_HomogMicro)
endsubroutine
END MODULE mesh
Reference in New Issue View Git Blame Copy Permalink