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

Merge remote-tracking branch 'remotes/origin/marc-improvements' into development

This commit is contained in:
Franz Roters 2020-01-18 12:21:24 +01:00
parent 58c3235820 f7d0b75f06
commit 701a9b18cc
11 changed files with 468 additions and 389 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

9
.gitlab-ci.yml
View File

@ -433,6 +433,15 @@ J2_plasticBehavior:
- master
- release
Marc_elementLib:
stage: marc
script:
- module load $IntelMarc $HDF5Marc $MSC
- Marc_elementLib/test.py
except:
- master
- release
###################################################################################################
grid_all_example:
stage: example

2
PRIVATE

@ -1 +1 @@
Subproject commit 99b076706a186ec7deb051ae181c2d9697c799fc
Subproject commit 4a8f3ccc2f9241b4afceece57ce71cae8025faa3

28
python/damask/dadf5.py
View File

@ -40,8 +40,9 @@ class DADF5():
self.version_major = f.attrs['DADF5-major']
self.version_minor = f.attrs['DADF5-minor']
if self.version_major != 0 or not 2 <= self.version_minor <= 5:
raise TypeError('Unsupported DADF5 version {} '.format(f.attrs['DADF5-version']))
if self.version_major != 0 or not 2 <= self.version_minor <= 6:
raise TypeError('Unsupported DADF5 version {}.{} '.format(f.attrs['DADF5_version_major'],
f.attrs['DADF5_version_minor']))
self.structured = 'grid' in f['geometry'].attrs.keys()
@ -886,8 +887,27 @@ class DADF5():
vtk_geom = vtk.vtkUnstructuredGrid()
vtk_geom.SetPoints(nodes)
vtk_geom.Allocate(f['/geometry/T_c'].shape[0])
if self.version_major == 0 and self.version_minor <= 5:
vtk_type = vtk.VTK_HEXAHEDRON
n_nodes = 8
else:
if f['/geometry/T_c'].attrs['VTK_TYPE'] == b'TRIANGLE':
vtk_type = vtk.VTK_TRIANGLE
n_nodes = 3
elif f['/geometry/T_c'].attrs['VTK_TYPE'] == b'QUAD':
vtk_type = vtk.VTK_QUAD
n_nodes = 4
elif f['/geometry/T_c'].attrs['VTK_TYPE'] == b'TETRA': # not tested
vtk_type = vtk.VTK_TETRA
n_nodes = 4
elif f['/geometry/T_c'].attrs['VTK_TYPE'] == b'HEXAHEDRON':
vtk_type = vtk.VTK_HEXAHEDRON
n_nodes = 8
for i in f['/geometry/T_c']:
vtk_geom.InsertNextCell(vtk.VTK_HEXAHEDRON,8,i-1) # not for all elements!
vtk_geom.InsertNextCell(n_nodes,vtk_type,i-1)
elif mode == 'Point':
Points = vtk.vtkPoints()
Vertices = vtk.vtkCellArray()
@ -936,7 +956,7 @@ class DADF5():
vtk_geom.GetCellData().AddArray(vtk_data[-1])
self.set_visible('materialpoints',materialpoints_backup)
constituents_backup = self.visible['constituents'].copy()
self.set_visible('constituents',False)
for label in (labels if isinstance(labels,list) else [labels]):

1
src/CPFEM.f90
View File

@ -377,6 +377,7 @@ subroutine CPFEM_results(inc,time)
call constitutive_results
call crystallite_results
call homogenization_results
call discretization_results
call results_finalizeIncrement
call results_closeJobFile

9
src/DAMASK_marc.f90
View File

@ -350,6 +350,9 @@ subroutine flux(f,ts,n,time)
!--------------------------------------------------------------------------------------------------
!> @brief trigger writing of results
!> @details uedinc is called before each new increment, not at the end of a converged one.
!> Therefore, storing the last written inc with an 'save' variable is required to avoid writing the
! same increment multiple times.
!--------------------------------------------------------------------------------------------------
subroutine uedinc(inc,incsub)
use prec
@ -357,8 +360,12 @@ subroutine uedinc(inc,incsub)
implicit none
integer, intent(in) :: inc, incsub
integer, save :: inc_written
#include QUOTE(PASTE(./MarcInclude/creeps,Marc4DAMASK)) ! creeps is needed for timinc (time increment)
call CPFEM_results(inc,cptim)
if (inc > inc_written) then
call CPFEM_results(inc,cptim)
inc_written = inc
endif
end subroutine uedinc

118
src/IO.f90
View File

@ -37,15 +37,12 @@ module IO
#if defined(Marc4DAMASK) || defined(Abaqus)
public :: &
IO_open_inputFile, &
IO_countContinuousIntValues, &
IO_continuousIntValues, &
IO_continuousIntValues
#endif
#if defined(Abaqus)
public :: &
IO_extractValue, &
IO_countDataLines
#elif defined(Marc4DAMASK)
IO_fixedNoEFloatValue, &
IO_fixedIntValue
#endif
#endif
contains
@ -57,7 +54,7 @@ contains
!--------------------------------------------------------------------------------------------------
subroutine IO_init
write(6,'(/,a)') ' <<<+- IO init -+>>>'
write(6,'(/,a)') ' <<<+- IO init -+>>>'; flush(6)
end subroutine IO_init
@ -444,57 +441,6 @@ integer function IO_intValue(string,chunkPos,myChunk)
end function IO_intValue
#ifdef Marc4DAMASK
!--------------------------------------------------------------------------------------------------
!> @brief reads float x.y+z value at myChunk from format string
!--------------------------------------------------------------------------------------------------
real(pReal) function IO_fixedNoEFloatValue (string,ends,myChunk)
character(len=*), intent(in) :: string !< raw input with known ends of each chunk
integer, intent(in) :: myChunk !< position number of desired chunk
integer, dimension(:), intent(in) :: ends !< positions of end of each tag/chunk in given string
character(len=*), parameter :: MYNAME = 'IO_fixedNoEFloatValue '
character(len=*), parameter :: VALIDBASE = '0123456789.+-'
character(len=*), parameter :: VALIDEXP = '0123456789+-'
real(pReal) :: base
integer :: expon
integer :: pos_exp
pos_exp = scan(string(ends(myChunk)+1:ends(myChunk+1)),'+-',back=.true.)
hasExponent: if (pos_exp > 1) then
base = verifyFloatValue(trim(adjustl(string(ends(myChunk)+1:ends(myChunk)+pos_exp-1))),&
VALIDBASE,MYNAME//'(base): ')
expon = verifyIntValue(trim(adjustl(string(ends(myChunk)+pos_exp:ends(myChunk+1)))),&
VALIDEXP,MYNAME//'(exp): ')
else hasExponent
base = verifyFloatValue(trim(adjustl(string(ends(myChunk)+1:ends(myChunk+1)))),&
VALIDBASE,MYNAME//'(base): ')
expon = 0
endif hasExponent
IO_fixedNoEFloatValue = base*10.0_pReal**real(expon,pReal)
end function IO_fixedNoEFloatValue
!--------------------------------------------------------------------------------------------------
!> @brief reads integer value at myChunk from fixed format string
!--------------------------------------------------------------------------------------------------
integer function IO_fixedIntValue(string,ends,myChunk)
character(len=*), intent(in) :: string !< raw input with known ends of each chunk
integer, intent(in) :: myChunk !< position number of desired chunk
integer, dimension(:), intent(in) :: ends !< positions of end of each tag/chunk in given string
character(len=*), parameter :: MYNAME = 'IO_fixedIntValue: '
character(len=*), parameter :: VALIDCHARACTERS = '0123456789+-'
IO_fixedIntValue = verifyIntValue(trim(adjustl(string(ends(myChunk)+1:ends(myChunk+1)))),&
VALIDCHARACTERS,MYNAME)
end function IO_fixedIntValue
#endif
!--------------------------------------------------------------------------------------------------
!> @brief changes characters in string to lower case
!--------------------------------------------------------------------------------------------------
@ -916,62 +862,6 @@ end function IO_countDataLines
#endif
!--------------------------------------------------------------------------------------------------
!> @brief count items in consecutive lines depending on lines
!> @details Marc: ints concatenated by "c" as last char or range of values a "to" b
!> Abaqus: triplet of start,stop,inc
!--------------------------------------------------------------------------------------------------
integer function IO_countContinuousIntValues(fileUnit)
integer, intent(in) :: fileUnit
#ifdef Abaqus
integer :: l,c
#endif
integer, allocatable, dimension(:) :: chunkPos
character(len=pStringLen) :: line
IO_countContinuousIntValues = 0
line = ''
#if defined(Marc4DAMASK)
do while (trim(line) /= IO_EOF)
line = IO_read(fileUnit)
chunkPos = IO_stringPos(line)
if (chunkPos(1) < 1) then ! empty line
exit
elseif (IO_lc(IO_stringValue(line,chunkPos,2)) == 'to' ) then ! found range indicator
IO_countContinuousIntValues = 1 + abs( IO_intValue(line,chunkPos,3) &
-IO_intValue(line,chunkPos,1))
exit ! only one single range indicator allowed
else
IO_countContinuousIntValues = IO_countContinuousIntValues+chunkPos(1)-1 ! add line's count when assuming 'c'
if ( IO_lc(IO_stringValue(line,chunkPos,chunkPos(1))) /= 'c' ) then ! line finished, read last value
IO_countContinuousIntValues = IO_countContinuousIntValues+1
exit ! data ended
endif
endif
enddo
#elif defined(Abaqus)
c = IO_countDataLines(fileUnit)
do l = 1,c
backspace(fileUnit)
enddo
l = 1
do while (trim(line) /= IO_EOF .and. l <= c) ! ToDo: is this correct?
l = l + 1
line = IO_read(fileUnit)
chunkPos = IO_stringPos(line)
IO_countContinuousIntValues = IO_countContinuousIntValues + 1 + & ! assuming range generation
(IO_intValue(line,chunkPos,2)-IO_intValue(line,chunkPos,1))/&
max(1,IO_intValue(line,chunkPos,3))
enddo
#endif
end function IO_countContinuousIntValues
!--------------------------------------------------------------------------------------------------
!> @brief return integer list corresponding to items in consecutive lines.
!! First integer in array is counter

10
src/discretization.f90
View File

@ -48,10 +48,10 @@ subroutine discretization_init(homogenizationAt,microstructureAt,&
real(pReal), dimension(:,:), intent(in) :: &
IPcoords0, &
NodeCoords0
integer, optional, intent(in) :: &
integer, optional, intent(in) :: &
sharedNodesBeginn
write(6,'(/,a)') ' <<<+- discretization init -+>>>'
write(6,'(/,a)') ' <<<+- discretization init -+>>>'; flush(6)
discretization_nElem = size(microstructureAt,1)
discretization_nIP = size(IPcoords0,2)/discretization_nElem
@ -81,15 +81,15 @@ subroutine discretization_results
real(pReal), dimension(:,:), allocatable :: u
call results_closeGroup(results_addGroup(trim('current/geometry')))
call results_closeGroup(results_addGroup('current/geometry'))
u = discretization_NodeCoords (1:3,:discretization_sharedNodesBeginn) &
- discretization_NodeCoords0(1:3,:discretization_sharedNodesBeginn)
call results_writeDataset('current/geometry',u,'u_n','nodal displacements','m')
call results_writeDataset('current/geometry',u,'u_n','displacements of the nodes','m')
u = discretization_IPcoords &
- discretization_IPcoords0
call results_writeDataset('current/geometry',u,'u_c','cell center displacements','m')
call results_writeDataset('current/geometry',u,'u_p','displacements of the materialpoints','m')
end subroutine discretization_results

498
src/element.f90
View File

@ -22,6 +22,8 @@ module element
NcellnodesPerCell, &
nIPs, &
nIPneighbors
character(len=:), allocatable :: &
vtkType
integer, dimension(:,:), allocatable :: &
Cell, & !< intra-element (cell) nodes that constitute a cell
IPneighbor, &
@ -36,10 +38,10 @@ module element
end type tElement
integer, parameter, private :: &
integer, parameter :: &
NELEMTYPE = 13
integer, dimension(NELEMTYPE), parameter, private :: NNODE = &
integer, dimension(NELEMTYPE), parameter :: NNODE = &
[ &
3, & ! 2D 3node 1ip
6, & ! 2D 6node 3ip
@ -57,7 +59,7 @@ module element
20 & ! 3D 20node 27ip
] !< number of nodes that constitute a specific type of element
integer, dimension(NELEMTYPE), parameter, public :: GEOMTYPE = &
integer, dimension(NELEMTYPE), parameter :: GEOMTYPE = &
[ &
1, &
2, &
@ -74,7 +76,7 @@ module element
10 &
] !< geometry type of particular element type
integer, dimension(maxval(GEOMTYPE)), parameter, private :: NCELLNODE = &
integer, dimension(maxval(GEOMTYPE)), parameter :: NCELLNODE = &
[ &
3, &
7, &
@ -88,7 +90,7 @@ module element
64 &
] !< number of cell nodes in a specific geometry type
integer, dimension(maxval(GEOMTYPE)), parameter, private :: NIP = &
integer, dimension(maxval(GEOMTYPE)), parameter :: NIP = &
[ &
1, &
3, &
@ -102,7 +104,7 @@ module element
27 &
] !< number of IPs in a specific geometry type
integer, dimension(maxval(GEOMTYPE)), parameter, private :: CELLTYPE = &
integer, dimension(maxval(GEOMTYPE)), parameter :: CELLTYPE = &
[ &
1, & ! 2D 3node
2, & ! 2D 4node
@ -116,7 +118,7 @@ module element
4 & ! 3D 8node
] !< cell type that is used by each geometry type
integer, dimension(maxval(CELLTYPE)), parameter, private :: NIPNEIGHBOR = &
integer, dimension(maxval(CELLTYPE)), parameter :: NIPNEIGHBOR = &
[ &
3, & ! 2D 3node
4, & ! 2D 4node
@ -124,7 +126,7 @@ module element
6 & ! 3D 8node
] !< number of ip neighbors / cell faces in a specific cell type
integer, dimension(maxval(CELLTYPE)), parameter, private :: NCELLNODEPERCELLFACE = &
integer, dimension(maxval(CELLTYPE)), parameter :: NCELLNODEPERCELLFACE = &
[ &
2, & ! 2D 3node
2, & ! 2D 4node
@ -132,7 +134,7 @@ module element
4 & ! 3D 8node
] !< number of cell nodes in a specific cell type
integer, dimension(maxval(CELLTYPE)), parameter, private :: NCELLNODEPERCELL = &
integer, dimension(maxval(CELLTYPE)), parameter :: NCELLNODEPERCELL = &
[ &
3, & ! 2D 3node
4, & ! 2D 4node
@ -146,27 +148,39 @@ module element
! Positive integers denote an intra-element IP identifier.
! Negative integers denote the interface behind which the neighboring (extra-element) IP will be located.
integer, dimension(nIPneighbor(cellType(1)),nIP(1)), parameter, private :: IPneighbor1 = &
integer, dimension(NIPNEIGHBOR(CELLTYPE(1)),NIP(1)), parameter :: IPNEIGHBOR1 = &
reshape([&
-2,-3,-1 &
],[nIPneighbor(cellType(1)),nIP(1)])
#if !defined(__GFORTRAN__)
],shape(IPNEIGHBOR1))
#else
],[NIPNEIGHBOR(CELLTYPE(1)),NIP(1)])
#endif
integer, dimension(nIPneighbor(cellType(2)),nIP(2)), parameter, private :: IPneighbor2 = &
integer, dimension(NIPNEIGHBOR(CELLTYPE(2)),NIP(2)), parameter :: IPNEIGHBOR2 = &
reshape([&
2,-3, 3,-1, &
-2, 1, 3,-1, &
2,-3,-2, 1 &
],[nIPneighbor(cellType(2)),nIP(2)])
#if !defined(__GFORTRAN__)
],shape(IPNEIGHBOR2))
#else
],[NIPNEIGHBOR(CELLTYPE(2)),NIP(2)])
#endif
integer, dimension(nIPneighbor(cellType(3)),nIP(3)), parameter, private :: IPneighbor3 = &
integer, dimension(NIPNEIGHBOR(CELLTYPE(3)),NIP(3)), parameter :: IPNEIGHBOR3 = &
reshape([&
2,-4, 3,-1, &
-2, 1, 4,-1, &
4,-4,-3, 1, &
-2, 3,-3, 2 &
],[nIPneighbor(cellType(3)),nIP(3)])
#if !defined(__GFORTRAN__)
],shape(IPNEIGHBOR3))
#else
],[NIPNEIGHBOR(CELLTYPE(3)),NIP(3)])
#endif
integer, dimension(nIPneighbor(cellType(4)),nIP(4)), parameter, private :: IPneighbor4 = &
integer, dimension(NIPNEIGHBOR(CELLTYPE(4)),NIP(4)), parameter :: IPNEIGHBOR4 = &
reshape([&
2,-4, 4,-1, &
3, 1, 5,-1, &
@ -177,22 +191,34 @@ module element
8,-4,-3, 4, &
9, 7,-3, 5, &
-2, 8,-3, 6 &
],[nIPneighbor(cellType(4)),nIP(4)])
#if !defined(__GFORTRAN__)
],shape(IPNEIGHBOR4))
#else
],[NIPNEIGHBOR(CELLTYPE(4)),NIP(4)])
#endif
integer, dimension(nIPneighbor(cellType(5)),nIP(5)), parameter, private :: IPneighbor5 = &
integer, dimension(NIPNEIGHBOR(CELLTYPE(5)),NIP(5)), parameter :: IPNEIGHBOR5 = &
reshape([&
-1,-2,-3,-4 &
],[nIPneighbor(cellType(5)),nIP(5)])
#if !defined(__GFORTRAN__)
],shape(IPNEIGHBOR5))
#else
],[NIPNEIGHBOR(CELLTYPE(5)),NIP(5)])
#endif
integer, dimension(nIPneighbor(cellType(6)),nIP(6)), parameter, private :: IPneighbor6 = &
integer, dimension(NIPNEIGHBOR(CELLTYPE(6)),NIP(6)), parameter :: IPNEIGHBOR6 = &
reshape([&
2,-4, 3,-2, 4,-1, &
-2, 1, 3,-2, 4,-1, &
2,-4,-3, 1, 4,-1, &
2,-4, 3,-2,-3, 1 &
],[nIPneighbor(cellType(6)),nIP(6)])
#if !defined(__GFORTRAN__)
],shape(IPNEIGHBOR6))
#else
],[NIPNEIGHBOR(CELLTYPE(6)),NIP(6)])
#endif
integer, dimension(nIPneighbor(cellType(7)),nIP(7)), parameter, private :: IPneighbor7 = &
integer, dimension(NIPNEIGHBOR(CELLTYPE(7)),NIP(7)), parameter :: IPNEIGHBOR7 = &
reshape([&
2,-4, 3,-2, 4,-1, &
-3, 1, 3,-2, 5,-1, &
@ -200,14 +226,22 @@ module element
5,-4, 6,-2,-5, 1, &
-3, 4, 6,-2,-5, 2, &
5,-4,-3, 4,-5, 3 &
],[nIPneighbor(cellType(7)),nIP(7)])
#if !defined(__GFORTRAN__)
],shape(IPNEIGHBOR7))
#else
],[NIPNEIGHBOR(CELLTYPE(7)),NIP(7)])
#endif
integer, dimension(nIPneighbor(cellType(8)),nIP(8)), parameter, private :: IPneighbor8 = &
integer, dimension(NIPNEIGHBOR(CELLTYPE(8)),NIP(8)), parameter :: IPNEIGHBOR8 = &
reshape([&
-3,-5,-4,-2,-6,-1 &
],[nIPneighbor(cellType(8)),nIP(8)])
-3,-5,-4,-2,-6,-1 &
#if !defined(__GFORTRAN__)
],shape(IPNEIGHBOR8))
#else
],[NIPNEIGHBOR(CELLTYPE(8)),NIP(8)])
#endif
integer, dimension(nIPneighbor(cellType(9)),nIP(9)), parameter, private :: IPneighbor9 = &
integer, dimension(NIPNEIGHBOR(CELLTYPE(9)),NIP(9)), parameter :: IPNEIGHBOR9 = &
reshape([&
2,-5, 3,-2, 5,-1, &
-3, 1, 4,-2, 6,-1, &
@ -217,9 +251,13 @@ module element
-3, 5, 8,-2,-6, 2, &
8,-5,-4, 5,-6, 3, &
-3, 7,-4, 6,-6, 4 &
],[nIPneighbor(cellType(9)),nIP(9)])
#if !defined(__GFORTRAN__)
],shape(IPNEIGHBOR9))
#else
],[NIPNEIGHBOR(CELLTYPE(9)),NIP(9)])
#endif
integer, dimension(nIPneighbor(cellType(10)),nIP(10)), parameter, private :: IPneighbor10 = &
integer, dimension(NIPNEIGHBOR(CELLTYPE(10)),NIP(10)), parameter :: IPNEIGHBOR10 = &
reshape([&
2,-5, 4,-2,10,-1, &
3, 1, 5,-2,11,-1, &
@ -248,17 +286,25 @@ module element
26,-5,-4,22,-6,16, &
27,25,-4,23,-6,17, &
-3,26,-4,24,-6,18 &
],[nIPneighbor(cellType(10)),nIP(10)])
#if !defined(__GFORTRAN__)
],shape(IPNEIGHBOR10))
#else
],[NIPNEIGHBOR(CELLTYPE(10)),NIP(10)])
#endif
integer, dimension(nNode(1),NcellNode(geomType(1))), parameter :: cellNodeParentNodeWeights1 = &
integer, dimension(NNODE(1),NCELLNODE(GEOMTYPE(1))), parameter :: CELLNODEPARENTNODEWEIGHTS1 = &
reshape([&
1, 0, 0, &
0, 1, 0, &
0, 0, 1 &
],[nNode(1),NcellNode(geomType(1))]) !< 2D 3node 1ip
#if !defined(__GFORTRAN__)
],shape(CELLNODEPARENTNODEWEIGHTS1)) !< 2D 3node 1ip
#else
],[NNODE(1),NCELLNODE(GEOMTYPE(1))])
#endif
integer, dimension(nNode(2),NcellNode(geomType(2))), parameter :: cellNodeParentNodeWeights2 = &
integer, dimension(NNODE(2),NCELLNODE(GEOMTYPE(2))), parameter :: CELLNODEPARENTNODEWEIGHTS2 = &
reshape([&
1, 0, 0, 0, 0, 0, &
0, 1, 0, 0, 0, 0, &
@ -267,9 +313,13 @@ module element
0, 0, 0, 0, 1, 0, &
0, 0, 0, 0, 0, 1, &
1, 1, 1, 2, 2, 2 &
],[nNode(2),NcellNode(geomType(2))]) !< 2D 6node 3ip
#if !defined(__GFORTRAN__)
],shape(CELLNODEPARENTNODEWEIGHTS2)) !< 2D 6node 3ip
#else
],[NNODE(2),NCELLNODE(GEOMTYPE(2))])
#endif
integer, dimension(nNode(3),NcellNode(geomType(3))), parameter :: cellNodeParentNodeWeights3 = &
integer, dimension(NNODE(3),NCELLNODE(GEOMTYPE(3))), parameter :: CELLNODEPARENTNODEWEIGHTS3 = &
reshape([&
1, 0, 0, 0, &
0, 1, 0, 0, &
@ -280,9 +330,13 @@ module element
0, 0, 1, 1, &
1, 0, 0, 1, &
1, 1, 1, 1 &
],[nNode(3),NcellNode(geomType(3))]) !< 2D 6node 3ip
#if !defined(__GFORTRAN__)
],shape(CELLNODEPARENTNODEWEIGHTS3)) !< 2D 6node 3ip
#else
],[NNODE(3),NCELLNODE(GEOMTYPE(3))])
#endif
integer, dimension(nNode(4),NcellNode(geomType(4))), parameter :: cellNodeParentNodeWeights4 = &
integer, dimension(NNODE(4),NCELLNODE(GEOMTYPE(4))), parameter :: CELLNODEPARENTNODEWEIGHTS4 = &
reshape([&
1, 0, 0, 0, 0, 0, 0, 0, &
0, 1, 0, 0, 0, 0, 0, 0, &
@ -300,9 +354,13 @@ module element
1, 4, 1, 1, 8, 8, 2, 2, &
1, 1, 4, 1, 2, 8, 8, 2, &
1, 1, 1, 4, 2, 2, 8, 8 &
],[nNode(4),NcellNode(geomType(4))]) !< 2D 8node 9ip
#if !defined(__GFORTRAN__)
],shape(CELLNODEPARENTNODEWEIGHTS4)) !< 2D 8node 9ip
#else
],[NNODE(4),NCELLNODE(GEOMTYPE(4))])
#endif
integer, dimension(nNode(5),NcellNode(geomType(5))), parameter :: cellNodeParentNodeWeights5 = &
integer, dimension(NNODE(5),NCELLNODE(GEOMTYPE(5))), parameter :: CELLNODEPARENTNODEWEIGHTS5 = &
reshape([&
1, 0, 0, 0, 0, 0, 0, 0, &
0, 1, 0, 0, 0, 0, 0, 0, &
@ -313,17 +371,25 @@ module element
0, 0, 0, 0, 0, 0, 1, 0, &
0, 0, 0, 0, 0, 0, 0, 1, &
1, 1, 1, 1, 2, 2, 2, 2 &
],[nNode(5),NcellNode(geomType(5))]) !< 2D 8node 4ip
#if !defined(__GFORTRAN__)
],shape(CELLNODEPARENTNODEWEIGHTS5)) !< 2D 8node 4ip
#else
],[NNODE(5),NCELLNODE(GEOMTYPE(5))])
#endif
integer, dimension(nNode(6),NcellNode(geomType(6))), parameter :: cellNodeParentNodeWeights6 = &
integer, dimension(NNODE(6),NcellNode(GEOMTYPE(6))), parameter :: CELLNODEPARENTNODEWEIGHTS6 = &
reshape([&
1, 0, 0, 0, &
0, 1, 0, 0, &
0, 0, 1, 0, &
0, 0, 0, 1 &
],[nNode(6),NcellNode(geomType(6))]) !< 3D 4node 1ip
#if !defined(__GFORTRAN__)
],shape(CELLNODEPARENTNODEWEIGHTS6)) !< 3D 4node 1ip
#else
],[NNODE(6),NcellNode(GEOMTYPE(6))])
#endif
integer, dimension(nNode(7),NcellNode(geomType(7))), parameter :: cellNodeParentNodeWeights7 = &
integer, dimension(NNODE(7),NCELLNODE(GEOMTYPE(7))), parameter :: CELLNODEPARENTNODEWEIGHTS7 = &
reshape([&
1, 0, 0, 0, 0, &
0, 1, 0, 0, 0, &
@ -340,9 +406,13 @@ module element
0, 1, 1, 1, 0, &
1, 0, 1, 1, 0, &
0, 0, 0, 0, 1 &
],[nNode(7),NcellNode(geomType(7))]) !< 3D 5node 4ip
#if !defined(__GFORTRAN__)
],shape(CELLNODEPARENTNODEWEIGHTS7)) !< 3D 5node 4ip
#else
],[NNODE(7),NCELLNODE(GEOMTYPE(7))])
#endif
integer, dimension(nNode(8),NcellNode(geomType(8))), parameter :: cellNodeParentNodeWeights8 = &
integer, dimension(NNODE(8),NCELLNODE(GEOMTYPE(8))), parameter :: CELLNODEPARENTNODEWEIGHTS8 = &
reshape([&
1, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
0, 1, 0, 0, 0, 0, 0, 0, 0, 0, &
@ -359,9 +429,13 @@ module element
0, 1, 1, 1, 0, 2, 0, 0, 2, 2, &
1, 0, 1, 1, 0, 0, 2, 2, 0, 2, &
3, 3, 3, 3, 4, 4, 4, 4, 4, 4 &
],[nNode(8),NcellNode(geomType(8))]) !< 3D 10node 4ip
#if !defined(__GFORTRAN__)
],shape(CELLNODEPARENTNODEWEIGHTS8)) !< 3D 10node 4ip
#else
],[NNODE(8),NCELLNODE(GEOMTYPE(8))])
#endif
integer, dimension(nNode(9),NcellNode(geomType(9))), parameter :: cellNodeParentNodeWeights9 = &
integer, dimension(NNODE(9),NCELLNODE(GEOMTYPE(9))), parameter :: CELLNODEPARENTNODEWEIGHTS9 = &
reshape([&
1, 0, 0, 0, 0, 0, &
0, 1, 0, 0, 0, 0, &
@ -384,9 +458,13 @@ module element
1, 0, 1, 1, 0, 1, &
0, 0, 0, 1, 1, 1, &
1, 1, 1, 1, 1, 1 &
],[nNode(9),NcellNode(geomType(9))]) !< 3D 6node 6ip
#if !defined(__GFORTRAN__)
],shape(CELLNODEPARENTNODEWEIGHTS9)) !< 3D 6node 6ip
#else
],[NNODE(9),NCELLNODE(GEOMTYPE(9))])
#endif
integer, dimension(nNode(10),NcellNode(geomType(10))), parameter :: cellNodeParentNodeWeights10 = &
integer, dimension(NNODE(10),NCELLNODE(GEOMTYPE(10))), parameter :: CELLNODEPARENTNODEWEIGHTS10 = &
reshape([&
1, 0, 0, 0, 0, 0, 0, 0, &
0, 1, 0, 0, 0, 0, 0, 0, &
@ -396,9 +474,13 @@ module element
0, 0, 0, 0, 0, 1, 0, 0, &
0, 0, 0, 0, 0, 0, 1, 0, &
0, 0, 0, 0, 0, 0, 0, 1 &
],[nNode(10),NcellNode(geomType(10))]) !< 3D 8node 1ip
#if !defined(__GFORTRAN__)
],shape(CELLNODEPARENTNODEWEIGHTS10)) !< 3D 8node 1ip
#else
],[NNODE(10),NCELLNODE(GEOMTYPE(10))])
#endif
integer, dimension(nNode(11),NcellNode(geomType(11))), parameter :: cellNodeParentNodeWeights11 = &
integer, dimension(NNODE(11),NCELLNODE(GEOMTYPE(11))), parameter :: CELLNODEPARENTNODEWEIGHTS11 = &
reshape([&
1, 0, 0, 0, 0, 0, 0, 0, & !
0, 1, 0, 0, 0, 0, 0, 0, & !
@ -427,9 +509,13 @@ module element
1, 0, 0, 1, 1, 0, 0, 1, & ! 25
0, 0, 0, 0, 1, 1, 1, 1, & !
1, 1, 1, 1, 1, 1, 1, 1 & !
],[nNode(11),NcellNode(geomType(11))]) !< 3D 8node 8ip
#if !defined(__GFORTRAN__)
],shape(CELLNODEPARENTNODEWEIGHTS11)) !< 3D 8node 8ip
#else
],[NNODE(11),NCELLNODE(GEOMTYPE(11))])
#endif
integer, dimension(nNode(12),NcellNode(geomType(12))), parameter :: cellNodeParentNodeWeights12 = &
integer, dimension(NNODE(12),NCELLNODE(GEOMTYPE(12))), parameter :: CELLNODEPARENTNODEWEIGHTS12 = &
reshape([&
1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & !
0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & !
@ -458,9 +544,13 @@ module element
1, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 2, 0, 0, 0, 2, 2, 0, 0, 2, & ! 25
0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 2, 2, 2, 2, 0, 0, 0, 0, & !
3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4 & !
],[nNode(12),NcellNode(geomType(12))]) !< 3D 20node 8ip
#if !defined(__GFORTRAN__)
],shape(CELLNODEPARENTNODEWEIGHTS12)) !< 3D 20node 8ip
#else
],[NNODE(12),NCELLNODE(GEOMTYPE(12))])
#endif
integer, dimension(nNode(13),NcellNode(geomType(13))), parameter :: cellNodeParentNodeWeights13 = &
integer, dimension(NNODE(13),NCELLNODE(GEOMTYPE(13))), parameter :: CELLNODEPARENTNODEWEIGHTS13 = &
reshape([&
1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & !
0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & !
@ -526,30 +616,46 @@ module element
4, 8, 4, 3, 8,24, 8, 4, 12,12, 4, 4, 32,32,12,12, 12,32,12, 4, & !
3, 4, 8, 4, 4, 8,24, 8, 4,12,12, 4, 12,32,32,12, 4,12,32,12, & !
4, 3, 4, 8, 8, 4, 8,24, 4, 4,12,12, 12,12,32,32, 12, 4,12,32 & !
],[nNode(13),NcellNode(geomType(13))]) !< 3D 20node 27ip
#if !defined(__GFORTRAN__)
],shape(CELLNODEPARENTNODEWEIGHTS13)) !< 3D 20node 27ip
#else
],[NNODE(13),NCELLNODE(GEOMTYPE(13))])
#endif
integer, dimension(NCELLNODEPERCELL(CELLTYPE(1)),NIP(1)), parameter :: CELL1 = &
integer, dimension(NCELLNODEPERCELL(CELLTYPE(1)),NIP(1)), parameter :: CELL1 = &
reshape([&
1,2,3 &
#if !defined(__GFORTRAN__)
],shape(CELL1))
#else
],[NCELLNODEPERCELL(CELLTYPE(1)),NIP(1)])
#endif
integer, dimension(NCELLNODEPERCELL(CELLTYPE(2)),NIP(2)), parameter :: CELL2 = &
integer, dimension(NCELLNODEPERCELL(CELLTYPE(2)),NIP(2)), parameter :: CELL2 = &
reshape([&
1, 4, 7, 6, &
2, 5, 7, 4, &
3, 6, 7, 5 &
#if !defined(__GFORTRAN__)
],shape(CELL2))
#else
],[NCELLNODEPERCELL(CELLTYPE(2)),NIP(2)])
#endif
integer, dimension(NCELLNODEPERCELL(CELLTYPE(3)),NIP(3)), parameter :: CELL3 = &
integer, dimension(NCELLNODEPERCELL(CELLTYPE(3)),NIP(3)), parameter :: CELL3 = &
reshape([&
1, 5, 9, 8, &
5, 2, 6, 9, &
8, 9, 7, 4, &
9, 6, 3, 7 &
#if !defined(__GFORTRAN__)
],shape(CELL3))
#else
],[NCELLNODEPERCELL(CELLTYPE(3)),NIP(3)])
#endif
integer, dimension(NCELLNODEPERCELL(CELLTYPE(4)),NIP(4)), parameter :: CELL4 = &
integer, dimension(NCELLNODEPERCELL(CELLTYPE(4)),NIP(4)), parameter :: CELL4 = &
reshape([&
1, 5,13,12, &
5, 6,14,13, &
@ -560,22 +666,34 @@ module element
11,16,10, 4, &
16,15, 9,10, &
15, 8, 3, 9 &
#if !defined(__GFORTRAN__)
],shape(CELL4))
#else
],[NCELLNODEPERCELL(CELLTYPE(4)),NIP(4)])
#endif
integer, dimension(NCELLNODEPERCELL(CELLTYPE(5)),NIP(5)), parameter :: CELL5 = &
integer, dimension(NCELLNODEPERCELL(CELLTYPE(5)),NIP(5)), parameter :: CELL5 = &
reshape([&
1, 2, 3, 4 &
#if !defined(__GFORTRAN__)
],shape(CELL5))
#else
],[NCELLNODEPERCELL(CELLTYPE(5)),NIP(5)])
#endif
integer, dimension(NCELLNODEPERCELL(CELLTYPE(6)),NIP(6)), parameter :: CELL6 = &
integer, dimension(NCELLNODEPERCELL(CELLTYPE(6)),NIP(6)), parameter :: CELL6 = &
reshape([&
1, 5,11, 7, 8,12,15,14, &
5, 2, 6,11,12, 9,13,15, &
7,11, 6, 3,14,15,13,10, &
8,12,15, 4, 4, 9,13,10 &
#if !defined(__GFORTRAN__)
],shape(CELL6))
#else
],[NCELLNODEPERCELL(CELLTYPE(6)),NIP(6)])
#endif
integer, dimension(NCELLNODEPERCELL(CELLTYPE(7)),NIP(7)), parameter :: CELL7 = &
integer, dimension(NCELLNODEPERCELL(CELLTYPE(7)),NIP(7)), parameter :: CELL7 = &
reshape([&
1, 7,16, 9,10,17,21,19, &
7, 2, 8,16,17,11,18,21, &
@ -583,14 +701,22 @@ module element
10,17,21,19, 4,13,20,15, &
17,11,18,21,13, 5,14,20, &
19,21,18,12,15,20,14, 6 &
#if !defined(__GFORTRAN__)
],shape(CELL7))
#else
],[NCELLNODEPERCELL(CELLTYPE(7)),NIP(7)])
#endif
integer, dimension(NCELLNODEPERCELL(CELLTYPE(8)),NIP(8)), parameter :: CELL8 = &
integer, dimension(NCELLNODEPERCELL(CELLTYPE(8)),NIP(8)), parameter :: CELL8 = &
reshape([&
1, 2, 3, 4, 5, 6, 7, 8 &
#if !defined(__GFORTRAN__)
],shape(CELL8))
#else
],[NCELLNODEPERCELL(CELLTYPE(8)),NIP(8)])
#endif
integer, dimension(NCELLNODEPERCELL(CELLTYPE(9)),NIP(9)), parameter :: CELL9 = &
integer, dimension(NCELLNODEPERCELL(CELLTYPE(9)),NIP(9)), parameter :: CELL9 = &
reshape([&
1, 9,21,12,13,22,27,25, &
9, 2,10,21,22,14,23,27, &
@ -600,7 +726,11 @@ module element
22,14,23,27,17, 6,18,26, &
25,27,24,16,20,26,19, 8, &
27,23,15,24,26,18, 7,19 &
#if !defined(__GFORTRAN__)
],shape(CELL9))
#else
],[NCELLNODEPERCELL(CELLTYPE(9)),NIP(9)])
#endif
integer, dimension(NCELLNODEPERCELL(CELLTYPE(10)),NIP(10)), parameter :: CELL10 = &
reshape([&
@ -631,7 +761,11 @@ module element
51,64,50,24,31,56,30, 8, &
64,63,49,50,56,55,29,30, &
63,48,23,49,55,28, 7,29 &
#if !defined(__GFORTRAN__)
],shape(CELL10))
#else
],[NCELLNODEPERCELL(CELLTYPE(10)),NIP(10)])
#endif
integer, dimension(NCELLNODEPERCELLFACE(1),NIPNEIGHBOR(1)), parameter :: CELLFACE1 = &
@ -639,7 +773,11 @@ module element
2,3, &
3,1, &
1,2 &
],[NCELLNODEPERCELLFACE(1),NIPNEIGHBOR(1)]) !< 2D 3node, VTK_TRIANGLE (5)
#if !defined(__GFORTRAN__)
],shape(CELLFACE1)) !< 2D 3node, VTK_TRIANGLE (5)
#else
],[NCELLNODEPERCELLFACE(1),NIPNEIGHBOR(1)])
#endif
integer, dimension(NCELLNODEPERCELLFACE(2),NIPNEIGHBOR(2)), parameter :: CELLFACE2 = &
reshape([&
@ -647,7 +785,11 @@ module element
4,1, &
3,4, &
1,2 &
],[NCELLNODEPERCELLFACE(2),NIPNEIGHBOR(2)]) !< 2D 4node, VTK_QUAD (9)
#if !defined(__GFORTRAN__)
],shape(CELLFACE2)) !< 2D 4node, VTK_QUAD (9)
#else
],[NCELLNODEPERCELLFACE(2),NIPNEIGHBOR(2)])
#endif
integer, dimension(NCELLNODEPERCELLFACE(3),NIPNEIGHBOR(3)), parameter :: CELLFACE3 = &
reshape([&
@ -655,7 +797,11 @@ module element
1,2,4, &
2,3,4, &
1,4,3 &
],[NCELLNODEPERCELLFACE(3),NIPNEIGHBOR(3)]) !< 3D 4node, VTK_TETRA (10)
#if !defined(__GFORTRAN__)
],shape(CELLFACE3)) !< 3D 4node, VTK_TETRA (10)
#else
],[NCELLNODEPERCELLFACE(3),NIPNEIGHBOR(3)])
#endif
integer, dimension(NCELLNODEPERCELLFACE(4),NIPNEIGHBOR(4)), parameter :: CELLFACE4 = &
reshape([&
@ -665,114 +811,128 @@ module element
1,2,6,5, &
5,6,7,8, &
1,4,3,2 &
],[NCELLNODEPERCELLFACE(4),NIPNEIGHBOR(4)]) !< 3D 8node, VTK_HEXAHEDRON (12)
#if !defined(__GFORTRAN__)
],shape(CELLFACE4)) !< 3D 8node, VTK_HEXAHEDRON (12)
#else
],[NCELLNODEPERCELLFACE(4),NIPNEIGHBOR(4)])
#endif
contains
!---------------------------------------------------------------------------------------------------
!> define properties of an element
!---------------------------------------------------------------------------------------------------
subroutine tElement_init(self,elemType)
class(tElement) :: self
integer, intent(in) :: elemType
self%elemType = elemType
self%Nnodes = Nnode (self%elemType)
self%geomType = geomType (self%elemType)
select case (self%elemType)
case(1)
self%cellNodeParentNodeWeights = cellNodeParentNodeWeights1
case(2)
self%cellNodeParentNodeWeights = cellNodeParentNodeWeights2
case(3)
self%cellNodeParentNodeWeights = cellNodeParentNodeWeights3
case(4)
self%cellNodeParentNodeWeights = cellNodeParentNodeWeights4
case(5)
self%cellNodeParentNodeWeights = cellNodeParentNodeWeights5
case(6)
self%cellNodeParentNodeWeights = cellNodeParentNodeWeights6
case(7)
self%cellNodeParentNodeWeights = cellNodeParentNodeWeights7
case(8)
self%cellNodeParentNodeWeights = cellNodeParentNodeWeights8
case(9)
self%cellNodeParentNodeWeights = cellNodeParentNodeWeights9
case(10)
self%cellNodeParentNodeWeights = cellNodeParentNodeWeights10
case(11)
self%cellNodeParentNodeWeights = cellNodeParentNodeWeights11
case(12)
self%cellNodeParentNodeWeights = cellNodeParentNodeWeights12
case(13)
self%cellNodeParentNodeWeights = cellNodeParentNodeWeights13
case default
call IO_error(0,ext_msg='invalid element type')
end select
self%NcellNodes = NcellNode (self%geomType)
self%nIPs = nIP (self%geomType)
self%cellType = cellType (self%geomType)
select case (self%geomType)
case(1)
self%IPneighbor = IPneighbor1
self%cell = CELL1
case(2)
self%IPneighbor = IPneighbor2
self%cell = CELL2
case(3)
self%IPneighbor = IPneighbor3
self%cell = CELL3
case(4)
self%IPneighbor = IPneighbor4
self%cell = CELL4
case(5)
self%IPneighbor = IPneighbor5
self%cell = CELL5
case(6)
self%IPneighbor = IPneighbor6
self%cell = CELL6
case(7)
self%IPneighbor = IPneighbor7
self%cell = CELL7
case(8)
self%IPneighbor = IPneighbor8
self%cell = CELL8
case(9)
self%IPneighbor = IPneighbor9
self%cell = CELL9
case(10)
self%IPneighbor = IPneighbor10
self%cell = CELL10
end select
self%NcellnodesPerCell = NCELLNODEPERCELL(self%cellType)
class(tElement) :: self
integer, intent(in) :: elemType
select case(self%cellType)
case(1)
self%cellFace = CELLFACE1
case(2)
self%cellFace = CELLFACE2
case(3)
self%cellFace = CELLFACE3
case(4)
self%cellFace = CELLFACE4
end select
self%nIPneighbors = size(self%IPneighbor,1)
write(6,'(/,a)') ' <<<+- element_init -+>>>'
write(6,*) ' element type: ',self%elemType
write(6,*) ' geom type: ',self%geomType
write(6,*) ' cell type: ',self%cellType
write(6,*) ' # node: ',self%Nnodes
write(6,*) ' # IP: ',self%nIPs
write(6,*) ' # cellnode: ',self%Ncellnodes
write(6,*) ' # cellnode/cell: ',self%NcellnodesPerCell
write(6,*) ' # IP neighbor: ',self%nIPneighbors
self%elemType = elemType
self%Nnodes = NNODE (self%elemType)
self%geomType = GEOMTYPE(self%elemType)
select case (self%elemType)
case(1)
self%cellNodeParentNodeWeights = CELLNODEPARENTNODEWEIGHTS1
case(2)
self%cellNodeParentNodeWeights = CELLNODEPARENTNODEWEIGHTS2
case(3)
self%cellNodeParentNodeWeights = CELLNODEPARENTNODEWEIGHTS3
case(4)
self%cellNodeParentNodeWeights = CELLNODEPARENTNODEWEIGHTS4
case(5)
self%cellNodeParentNodeWeights = CELLNODEPARENTNODEWEIGHTS5
case(6)
self%cellNodeParentNodeWeights = CELLNODEPARENTNODEWEIGHTS6
case(7)
self%cellNodeParentNodeWeights = CELLNODEPARENTNODEWEIGHTS7
case(8)
self%cellNodeParentNodeWeights = CELLNODEPARENTNODEWEIGHTS8
case(9)
self%cellNodeParentNodeWeights = CELLNODEPARENTNODEWEIGHTS9
case(10)
self%cellNodeParentNodeWeights = CELLNODEPARENTNODEWEIGHTS10
case(11)
self%cellNodeParentNodeWeights = CELLNODEPARENTNODEWEIGHTS11
case(12)
self%cellNodeParentNodeWeights = CELLNODEPARENTNODEWEIGHTS12
case(13)
self%cellNodeParentNodeWeights = CELLNODEPARENTNODEWEIGHTS13
case default
call IO_error(0,ext_msg='invalid element type')
end select
self%NcellNodes = NCELLNODE(self%geomType)
self%nIPs = NIP (self%geomType)
self%cellType = CELLTYPE (self%geomType)
select case (self%geomType)
case(1)
self%IPneighbor = IPNEIGHBOR1
self%cell = CELL1
case(2)
self%IPneighbor = IPNEIGHBOR2
self%cell = CELL2
case(3)
self%IPneighbor = IPNEIGHBOR3
self%cell = CELL3
case(4)
self%IPneighbor = IPNEIGHBOR4
self%cell = CELL4
case(5)
self%IPneighbor = IPNEIGHBOR5
self%cell = CELL5
case(6)
self%IPneighbor = IPNEIGHBOR6
self%cell = CELL6
case(7)
self%IPneighbor = IPNEIGHBOR7
self%cell = CELL7
case(8)
self%IPneighbor = IPNEIGHBOR8
self%cell = CELL8
case(9)
self%IPneighbor = IPNEIGHBOR9
self%cell = CELL9
case(10)
self%IPneighbor = IPNEIGHBOR10
self%cell = CELL10
end select
self%NcellnodesPerCell = NCELLNODEPERCELL(self%cellType)
select case(self%cellType)
case(1)
self%cellFace = CELLFACE1
self%vtkType = 'TRIANGLE'
case(2)
self%cellFace = CELLFACE2
self%vtkType = 'QUAD'
case(3)
self%cellFace = CELLFACE3
self%vtkType = 'TETRA'
case(4)
self%cellFace = CELLFACE4
self%vtkType = 'HEXAHEDRON'
end select
self%nIPneighbors = size(self%IPneighbor,1)
write(6,'(/,a)') ' <<<+- element_init -+>>>'; flush(6)
write(6,*) ' element type: ',self%elemType
write(6,*) ' geom type: ',self%geomType
write(6,*) ' cell type: ',self%cellType
write(6,*) ' # node: ',self%Nnodes
write(6,*) ' # IP: ',self%nIPs
write(6,*) ' # cellnode: ',self%Ncellnodes
write(6,*) ' # cellnode/cell: ',self%NcellnodesPerCell
write(6,*) ' # IP neighbor: ',self%nIPneighbors
end subroutine tElement_init
end module element

2
src/mesh_grid.f90
View File

@ -61,7 +61,7 @@ subroutine mesh_init(ip,el)
integer(C_INTPTR_T) :: &
devNull, z, z_offset
write(6,'(/,a)') ' <<<+- mesh_grid init -+>>>'; flush(6)
write(6,'(/,a)') ' <<<+- mesh_grid init -+>>>'; flush(6)
call readGeom(grid,geomSize,origin,microstructureAt,homogenizationAt)

178
src/mesh_marc.f90
View File

@ -72,7 +72,7 @@ subroutine mesh_init(ip,el)
real(pReal), dimension(:,:,:,:),allocatable :: &
unscaledNormals
write(6,'(/,a)') ' <<<+- mesh init -+>>>'
write(6,'(/,a)') ' <<<+- mesh init -+>>>'; flush(6)
mesh_unitlength = numerics_unitlength ! set physical extent of a length unit in mesh
@ -82,9 +82,8 @@ subroutine mesh_init(ip,el)
if (debug_e < 1 .or. debug_e > nElems) call IO_error(602,ext_msg='element')
if (debug_i < 1 .or. debug_i > elem%nIPs) call IO_error(602,ext_msg='IP')
FEsolving_execElem = [ 1,nElems ] ! parallel loop bounds set to comprise all DAMASK elements
allocate(FEsolving_execIP(2,nElems), source=1) ! parallel loop bounds set to comprise from first IP...
FEsolving_execIP(2,:) = elem%nIPs
FEsolving_execElem = [1,nElems]
FEsolving_execIP = spread([1,elem%nIPs],2,nElems)
allocate(calcMode(elem%nIPs,nElems),source=.false.) ! pretend to have collected what first call is asking (F = I)
calcMode(ip,mesh_FEasCP('elem',el)) = .true. ! first ip,el needs to be already pingponged to "calc"
@ -103,9 +102,9 @@ subroutine mesh_init(ip,el)
call discretization_init(microstructureAt,homogenizationAt,&
ip_reshaped,&
node0_elem)
node0_cell)
call writeGeometry(0,connectivity_elem,&
call writeGeometry(elem,connectivity_elem,&
reshape(connectivity_cell,[elem%NcellNodesPerCell,elem%nIPs*nElems]),&
node0_cell,ip_reshaped)
@ -121,13 +120,14 @@ end subroutine mesh_init
!--------------------------------------------------------------------------------------------------
!> @brief Writes all information needed for the DADF5 geometry
!> @brief Write all information needed for the DADF5 geometry
!--------------------------------------------------------------------------------------------------
subroutine writeGeometry(elemType, &
subroutine writeGeometry(elem, &
connectivity_elem,connectivity_cell, &
coordinates_nodes,coordinates_points)
integer, intent(in) :: elemType
type(tElement), intent(in) :: &
elem
integer, dimension(:,:), intent(in) :: &
connectivity_elem, &
connectivity_cell
@ -150,28 +150,31 @@ subroutine writeGeometry(elemType, &
connectivity_temp = connectivity_cell
call results_writeDataset('geometry',connectivity_temp,'T_c', &
'connectivity of the cells','-')
call results_addAttribute('VTK_TYPE',elem%vtkType,'geometry/T_c')
coordinates_temp = coordinates_nodes
call results_writeDataset('geometry',coordinates_temp,'x_n', &
'coordinates of the nodes','m')
'initial coordinates of the nodes','m')
coordinates_temp = coordinates_points
call results_writeDataset('geometry',coordinates_temp,'x_p', &
'coordinates of the material points','m')
'initial coordinates of the materialpoints','m')
call results_closeJobFile
end subroutine writeGeometry
!--------------------------------------------------------------------------------------------------
!> @brief Read mesh from marc input file
!--------------------------------------------------------------------------------------------------
subroutine inputRead(elem,node0_elem,connectivity_elem,microstructureAt,homogenizationAt)
type(tElement), intent(out) :: elem
real(pReal), dimension(:,:), allocatable, intent(out) :: &
node0_elem !< node x,y,z coordinates (initially!)
integer, dimension(:,:), allocatable, intent(out) :: &
integer, dimension(:,:), allocatable, intent(out) :: &
connectivity_elem
integer, dimension(:), allocatable, intent(out) :: &
microstructureAt, &
homogenizationAt
@ -207,16 +210,18 @@ subroutine inputRead(elem,node0_elem,connectivity_elem,microstructureAt,homogeni
call IO_open_inputFile(FILEUNIT) ! ToDo: It would be better to use fileContent
call inputRead_mapElemSets(nameElemSet,mapElemSet,&
FILEUNIT)
FILEUNIT,inputFile)
call inputRead_elemType(elem, &
nElems,inputFile)
allocate (mesh_mapFEtoCPelem(2,nElems), source=0)
call inputRead_mapElems(hypoelasticTableStyle,nameElemSet,mapElemSet,fileFormatVersion,matNumber,FILEUNIT)
call inputRead_mapElems(elem%nNodes,nElems,&
inputFile)
allocate (mesh_mapFEtoCPnode(2,Nnodes), source=0)
call inputRead_mapNodes(inputFile)
call inputRead_elemType(elem, &
nElems,inputFile)
call inputRead_elemNodes(node0_elem, &
Nnodes,inputFile)
@ -347,47 +352,63 @@ end subroutine inputRead_NnodesAndElements
!> @brief Count overall number of element sets in mesh.
!--------------------------------------------------------------------------------------------------
subroutine inputRead_NelemSets(nElemSets,maxNelemInSet,&
fileUnit)
fileContent)
integer, intent(out) :: nElemSets, maxNelemInSet
integer, intent(in) :: fileUnit
integer, intent(out) :: nElemSets, maxNelemInSet
character(len=pStringLen), dimension(:), intent(in) :: fileContent !< file content, separated per lines
integer, allocatable, dimension(:) :: chunkPos
character(len=300) :: line
integer :: i,l,elemInCurrentSet
nElemSets = 0
maxNelemInSet = 0
rewind(fileUnit)
do
read (fileUnit,'(A300)',END=620) line
chunkPos = IO_stringPos(line)
if ( IO_lc(IO_StringValue(line,chunkPos,1)) == 'define' .and. &
IO_lc(IO_StringValue(line,chunkPos,2)) == 'element' ) then
do l = 1, size(fileContent)
chunkPos = IO_stringPos(fileContent(l))
if ( IO_lc(IO_StringValue(fileContent(l),chunkPos,1)) == 'define' .and. &
IO_lc(IO_StringValue(fileContent(l),chunkPos,2)) == 'element' ) then
nElemSets = nElemSets + 1
maxNelemInSet = max(maxNelemInSet, IO_countContinuousIntValues(fileUnit))
chunkPos = IO_stringPos(fileContent(l+1))
if(IO_lc(IO_StringValue(fileContent(l+1),chunkPos,2)) == 'to' ) then
elemInCurrentSet = 1 + abs( IO_intValue(fileContent(l+1),chunkPos,3) &
-IO_intValue(fileContent(l+1),chunkPos,1))
else
elemInCurrentSet = 0
i = 0
do while (.true.)
i = i + 1
chunkPos = IO_stringPos(fileContent(l+i))
elemInCurrentSet = elemInCurrentSet + chunkPos(1) - 1 ! add line's count when assuming 'c'
if(IO_lc(IO_stringValue(fileContent(l+i),chunkPos,chunkPos(1))) /= 'c') then ! line finished, read last value
elemInCurrentSet = elemInCurrentSet + 1 ! data ended
exit
endif
enddo
endif
maxNelemInSet = max(maxNelemInSet, elemInCurrentSet)
endif
enddo
620 end subroutine inputRead_NelemSets
end subroutine inputRead_NelemSets
!--------------------------------------------------------------------------------------------------
!> @brief map element sets
!--------------------------------------------------------------------------------------------------
subroutine inputRead_mapElemSets(nameElemSet,mapElemSet,fileUnit)
subroutine inputRead_mapElemSets(nameElemSet,mapElemSet,fileUnit,fileContent)
character(len=64), dimension(:), allocatable, intent(out) :: nameElemSet
integer, dimension(:,:), allocatable, intent(out) :: mapElemSet
integer, intent(in) :: fileUnit
character(len=pStringLen), dimension(:), intent(in) :: fileContent !< file content, separated per lines
integer, allocatable, dimension(:) :: chunkPos
character(len=300) :: line
integer :: elemSet, NelemSets, maxNelemInSet
call inputRead_NelemSets(NelemSets,maxNelemInSet,fileUnit)
call inputRead_NelemSets(NelemSets,maxNelemInSet,fileContent)
allocate(nameElemSet(NelemSets)); nameElemSet = 'n/a'
allocate(mapElemSet(1+maxNelemInSet,NelemSets),source=0)
elemSet = 0
@ -411,65 +432,35 @@ subroutine inputRead_mapElemSets(nameElemSet,mapElemSet,fileUnit)
!--------------------------------------------------------------------------------------------------
!> @brief Maps elements from FE ID to internal (consecutive) representation.
!--------------------------------------------------------------------------------------------------
subroutine inputRead_mapElems(tableStyle,nameElemSet,mapElemSet,fileFormatVersion,matNumber,fileUnit)
integer, intent(in) :: fileUnit,tableStyle,fileFormatVersion
integer, dimension(:), intent(in) :: matNumber
character(len=64), intent(in), dimension(:) :: nameElemSet
integer, dimension(:,:), intent(in) :: &
mapElemSet
subroutine inputRead_mapElems(nNodes,nElem,fileContent)
integer, intent(in) :: &
nElem, &
nNodes !< number of nodes per element
character(len=pStringLen), dimension(:), intent(in) :: fileContent !< file content, separated per lines
integer, allocatable, dimension(:) :: chunkPos
character(len=300) :: line, &
tmp
integer, dimension(:), allocatable :: contInts
integer :: i,cpElem
allocate(contInts(size(mesh_mapFEtoCPelem,2)+1))
cpElem = 0
contInts = 0
rewind(fileUnit)
do
read (fileUnit,'(A300)',END=620) line
chunkPos = IO_stringPos(line)
Marc2016andEarlier: if (fileFormatVersion < 13) then
if( IO_lc(IO_stringValue(line,chunkPos,1)) == 'hypoelastic' ) then
skipLines: do i=1,3+TableStyle
read (fileUnit,'(A300)') line
enddo skipLines
contInts = IO_continuousIntValues(fileUnit,size(mesh_mapFEtoCPelem,2),nameElemSet,&
mapElemSet,size(nameElemSet))
exit
endif
else Marc2016andEarlier
if ( IO_lc(IO_stringValue(line,chunkPos,1)) == 'connectivity') then
read (fileUnit,'(A300)',END=620) line
chunkPos = IO_stringPos(line)
if(any(matNumber==IO_intValue(line,chunkPos,6))) then
do
read (fileUnit,'(A300)',END=620) line
chunkPos = IO_stringPos(line)
tmp = IO_lc(IO_stringValue(line,chunkPos,1))
if (verify(trim(tmp),"0123456789")/=0) then ! found keyword
exit
else
contInts(1) = contInts(1) + 1
read (tmp,*) contInts(contInts(1)+1)
endif
enddo
endif
endif
endif Marc2016andEarlier
enddo
620 do i = 1,contInts(1)
cpElem = cpElem+1
mesh_mapFEtoCPelem(1,cpElem) = contInts(1+i)
mesh_mapFEtoCPelem(2,cpElem) = cpElem
integer :: i,j,l,nNodesAlreadyRead
do l = 1, size(fileContent)
chunkPos = IO_stringPos(fileContent(l))
if( IO_lc(IO_stringValue(fileContent(l),chunkPos,1)) == 'connectivity' ) then
j = 0
do i = 1,nElem
chunkPos = IO_stringPos(fileContent(l+1+i+j))
mesh_mapFEtoCPelem(:,i) = [IO_intValue(fileContent(l+1+i+j),chunkPos,1),i]
nNodesAlreadyRead = chunkPos(1) - 2
do while(nNodesAlreadyRead < nNodes) ! read on if not all nodes in one line
j = j + 1
chunkPos = IO_stringPos(fileContent(l+1+i+j))
nNodesAlreadyRead = nNodesAlreadyRead + chunkPos(1)
enddo
enddo
exit
endif
enddo
call math_sort(mesh_mapFEtoCPelem)
call math_sort(mesh_mapFEtoCPelem)
end subroutine inputRead_mapElems
@ -488,7 +479,8 @@ subroutine inputRead_mapNodes(fileContent)
chunkPos = IO_stringPos(fileContent(l))
if( IO_lc(IO_stringValue(fileContent(l),chunkPos,1)) == 'coordinates' ) then
do i = 1,size(mesh_mapFEtoCPnode,2)
mesh_mapFEtoCPnode(1:2,i) = [IO_fixedIntValue (fileContent(l+1+i),[0,10],1),i] ! ToDo: use IO_intValue
chunkPos = IO_stringPos(fileContent(l+1+i))
mesh_mapFEtoCPnode(:,i) = [IO_intValue(fileContent(l+1+i),chunkPos,1),i]
enddo
exit
endif
@ -509,7 +501,6 @@ subroutine inputRead_elemNodes(nodes, &
integer, intent(in) :: nNode
character(len=pStringLen), dimension(:), intent(in) :: fileContent !< file content, separated per lines
integer, dimension(5), parameter :: node_ends = [0,10,30,50,70]
integer, allocatable, dimension(:) :: chunkPos
integer :: i,j,m,l
@ -519,9 +510,10 @@ subroutine inputRead_elemNodes(nodes, &
chunkPos = IO_stringPos(fileContent(l))
if( IO_lc(IO_stringValue(fileContent(l),chunkPos,1)) == 'coordinates' ) then
do i=1,nNode
m = mesh_FEasCP('node',IO_fixedIntValue(fileContent(l+1+i),node_ends,1)) !ToDo: use IO_intValue
chunkPos = IO_stringPos(fileContent(l+1+i))
m = mesh_FEasCP('node',IO_intValue(fileContent(l+1+i),chunkPos,1))
do j = 1,3
nodes(j,m) = mesh_unitlength * IO_fixedNoEFloatValue(fileContent(l+1+i),node_ends,j+1) !ToDo: use IO_floatValue
nodes(j,m) = mesh_unitlength * IO_floatValue(fileContent(l+1+i),chunkPos,j+1)
enddo
enddo
exit
@ -709,7 +701,7 @@ subroutine inputRead_microstructureAndHomogenization(microstructureAt,homogeniza
read (fileUnit,'(A300)',END=630) line ! read line with value of state var
chunkPos = IO_stringPos(line)
do while (scan(IO_stringValue(line,chunkPos,1),'+-',back=.true.)>1) ! is noEfloat value?
myVal = nint(IO_fixedNoEFloatValue(line,[0,20],1),pInt) ! state var's value
myVal = nint(IO_floatValue(line,chunkPos,1))
if (initialcondTableStyle == 2) then
read (fileUnit,'(A300)',END=630) line ! read extra line
read (fileUnit,'(A300)',END=630) line ! read extra line

2
src/results.f90
View File

@ -70,7 +70,7 @@ subroutine results_init
resultsFile = HDF5_openFile(trim(getSolverJobName())//'.hdf5','w',.true.)
call results_addAttribute('DADF5_version_major',0)
call results_addAttribute('DADF5_version_minor',5)
call results_addAttribute('DADF5_version_minor',6)
call results_addAttribute('DAMASK_version',DAMASKVERSION)
call get_command(commandLine)
call results_addAttribute('call',trim(commandLine))