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

Merge branch 'MiscImprovements' into development

This commit is contained in:
Vitesh Shah 2020-01-27 14:30:50 +01:00
parent c8835947db 9a0d8db23e
commit c433e2448a
30 changed files with 685 additions and 785 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
.gitignore vendored
View File

@ -1,5 +1,6 @@
*.pyc
*.hdf5
*.xdmf
*.bak
*~
.DS_Store

2
PRIVATE

@ -1 +1 @@
Subproject commit 83cb0b6ceea2e26c554ad4c6fbe6aabaee6fa27f
Subproject commit 2077cffcf72afd273751c85dfb77bf7c7f372575

67
processing/post/addIndexed.py
View File

@ -2,6 +2,7 @@
import os
import sys
from io import StringIO
from optparse import OptionParser
import numpy as np
@ -43,71 +44,21 @@ parser.set_defaults(offset = 0,
)
(options,filenames) = parser.parse_args()
if filenames == []: filenames = [None]
if options.label is None:
parser.error('no data columns specified.')
if options.index is None:
parser.error('no index column given.')
# ------------------------------------------ process indexed ASCIItable ---------------------------
if options.asciitable is not None and os.path.isfile(options.asciitable):
indexedTable = damask.ASCIItable(name = options.asciitable,
buffered = False,
readonly = True)
indexedTable.head_read() # read ASCII header info of indexed table
missing_labels = indexedTable.data_readArray(options.label)
indexedTable.close() # close input ASCII table
if len(missing_labels) > 0:
damask.util.croak('column{} {} not found...'.format('s' if len(missing_labels) > 1 else '',', '.join(missing_labels)))
else:
parser.error('no indexed ASCIItable given.')
# --- loop over input files -------------------------------------------------------------------------
if filenames == []: filenames = [None]
for name in filenames:
try: table = damask.ASCIItable(name = name,
buffered = False)
except: continue
damask.util.report(scriptName,name)
damask.util.report(scriptName,name)
# ------------------------------------------ read header ------------------------------------------
table = damask.Table.from_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name)
indexedTable = damask.Table.from_ASCII(options.asciitable)
idx = np.reshape(table.get(options.index).astype(int) + options.offset,(-1))-1
table.head_read()
for data in options.label:
table.add(data+'addIndexed',indexedTable.get(data)[idx],scriptID+' '+' '.join(sys.argv[1:]))
# ------------------------------------------ sanity checks ----------------------------------------
errors = []
indexColumn = table.label_index(options.index)
if indexColumn < 0: errors.append('index column {} not found.'.format(options.index))
if errors != []:
damask.util.croak(errors)
table.close(dismiss = True)
continue
# ------------------------------------------ assemble header --------------------------------------
table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:]))
table.labels_append(indexedTable.labels(raw = True)) # extend ASCII header with new labels
table.head_write()
# ------------------------------------------ process data ------------------------------------------
outputAlive = True
while outputAlive and table.data_read(): # read next data line of ASCII table
try:
table.data_append(indexedTable.data[int(round(float(table.data[indexColumn])))+options.offset-1]) # add all mapped data types
except IndexError:
table.data_append(np.nan*np.ones_like(indexedTable.data[0]))
outputAlive = table.data_write() # output processed line
# ------------------------------------------ output finalization -----------------------------------
table.close() # close ASCII tables
table.to_ASCII(sys.stdout if name is None else name)

95
processing/pre/geom_fromTable.py
View File

@ -2,6 +2,7 @@
import os
import sys
from io import StringIO
from optparse import OptionParser
import numpy as np
@ -54,80 +55,54 @@ parser.set_defaults(homogenization = 1,
)
(options,filenames) = parser.parse_args()
if filenames == []: filenames = [None]
input = [options.quaternion is not None,
options.microstructure is not None,
]
if np.sum(input) != 1:
if np.sum([options.quaternion is not None,
options.microstructure is not None]) != 1:
parser.error('need either microstructure or quaternion (and optionally phase) as input.')
if options.microstructure is not None and options.phase is not None:
parser.error('need either microstructure or phase (and mandatory quaternion) as input.')
if options.axes is not None and not set(options.axes).issubset(set(['x','+x','-x','y','+y','-y','z','+z','-z'])):
parser.error('invalid axes {} {} {}.'.format(*options.axes))
(label,dim,inputtype) = [(options.quaternion,4,'quaternion'),
(options.microstructure,1,'microstructure'),
][np.where(input)[0][0]] # select input label that was requested
if filenames == []: filenames = [None]
for name in filenames:
damask.util.report(scriptName,name)
table = damask.ASCIItable(name = name,readonly=True)
table.head_read() # read ASCII header info
damask.util.report(scriptName,name)
table = damask.Table.from_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name)
table.sort_by(['{}_{}'.format(i,options.pos) for i in range(3,0,-1)]) # x fast, y slow
grid,size,origin = damask.grid_filters.cell_coord0_gridSizeOrigin(table.get(options.pos))
table.data_readArray([options.pos] \
+ (label if isinstance(label, list) else [label]) \
+ ([options.phase] if options.phase else []))
config_header = table.comments
if options.phase is None:
table.data = np.column_stack((table.data,np.ones(len(table.data)))) # add single phase if no phase column given
if options.microstructure:
microstructure = table.get(options.microstructure).reshape(grid,order='F')
grid,size,origin = damask.grid_filters.cell_coord0_gridSizeOrigin(table.data[:,0:3])
elif options.quaternion:
q = table.get(options.quaternion)
phase = table.get(options.phase).astype(int) if options.phase else \
np.ones((table.data.shape[0],1),dtype=int)
indices = np.lexsort((table.data[:,0],table.data[:,1],table.data[:,2])) # indices of position when sorting x fast, z slow
microstructure = np.empty(grid,dtype = int) # initialize empty microstructure
i = 0
unique,unique_inverse = np.unique(np.hstack((q,phase)),return_inverse=True,axis=0)
microstructure = unique_inverse.reshape(grid,order='F') + 1
if inputtype == 'microstructure':
for z in range(grid[2]):
for y in range(grid[1]):
for x in range(grid[0]):
microstructure[x,y,z] = table.data[indices[i],3]
i+=1
config_header = ['<texture>']
for i,data in enumerate(unique):
ori = damask.Rotation(data[0:4])
config_header += ['[Grain{}]'.format(i+1),
'(gauss)\tphi1 {:.2f}\tPhi {:.2f}\tphi2 {:.2f}'.format(*ori.asEulers(degrees = True)),
]
if options.axes is not None: config_header += ['axes\t{} {} {}'.format(*options.axes)]
config_header = []
config_header += ['<microstructure>']
for i,data in enumerate(unique):
config_header += ['[Grain{}]'.format(i+1),
'(constituent)\tphase {}\ttexture {}\tfraction 1.0'.format(int(data[4]),i+1),
]
elif inputtype == 'quaternion':
unique,unique_inverse = np.unique(table.data[:,3:8],return_inverse=True,axis=0)
header = [scriptID + ' ' + ' '.join(sys.argv[1:])]\
+ config_header
geom = damask.Geom(microstructure,size,origin,
homogenization=options.homogenization,comments=header)
damask.util.croak(geom)
for z in range(grid[2]):
for y in range(grid[1]):
for x in range(grid[0]):
microstructure[x,y,z] = unique_inverse[indices[i]]+1
i+=1
config_header = ['<texture>']
for i,data in enumerate(unique):
ori = damask.Rotation(data[0:4])
config_header += ['[Grain{}]'.format(i+1),
'(gauss)\tphi1 {:.2f}\tPhi {:.2f}\tphi2 {:.2f}'.format(*ori.asEulers(degrees = True)),
]
if options.axes is not None: config_header += ['axes\t{} {} {}'.format(*options.axes)]
config_header += ['<microstructure>']
for i,data in enumerate(unique):
config_header += ['[Grain{}]'.format(i+1),
'(constituent)\tphase {}\ttexture {}\tfraction 1.0'.format(int(data[4]),i+1),
]
header = [scriptID + ' ' + ' '.join(sys.argv[1:])]\
+ config_header
geom = damask.Geom(microstructure,size,origin,
homogenization=options.homogenization,comments=header)
damask.util.croak(geom)
geom.to_file(sys.stdout if name is None else os.path.splitext(name)[0]+'.geom',pack=False)
geom.to_file(sys.stdout if name is None else os.path.splitext(name)[0]+'.geom',pack=False)

2
python/damask/dadf5.py
View File

@ -906,7 +906,7 @@ class DADF5():
n_nodes = 8
for i in f['/geometry/T_c']:
vtk_geom.InsertNextCell(n_nodes,vtk_type,i-1)
vtk_geom.InsertNextCell(vtk_type,n_nodes,i-1)
elif mode == 'Point':
Points = vtk.vtkPoints()

4
python/damask/orientation.py
View File

@ -1078,12 +1078,12 @@ class Orientation:
if isinstance(lattice, Lattice):
self.lattice = lattice
else:
self.lattice = Lattice(lattice) # assume string
self.lattice = Lattice(lattice) # assume string
if isinstance(rotation, Rotation):
self.rotation = rotation
else:
self.rotation = Rotation(rotation) # assume quaternion
self.rotation = Rotation.fromQuaternion(rotation) # assume quaternion
def disorientation(self,
other,

46
python/damask/table.py
View File

@ -14,8 +14,8 @@ class Table():
Parameters
----------
data : numpy.ndarray
Data.
data : numpy.ndarray or pandas.DataFrame
Data. Column labels from a pandas.DataFrame will be replaced.
shapes : dict with str:tuple pairs
Shapes of the columns. Example 'F':(3,3) for a deformation gradient.
comments : iterable of str, optional
@ -48,7 +48,7 @@ class Table():
def __add_comment(self,label,shape,info):
if info is not None:
self.comments.append('{}{}: {}'.format(label,
' '+str(shape) if np.prod(shape,dtype=int) > 1 else '',
' '+str(shape) if np.prod(shape) > 1 else '',
info))
@ -73,14 +73,21 @@ class Table():
f = fname
f.seek(0)
header,keyword = f.readline().split()
if keyword == 'header':
header = int(header)
else:
raise TypeError
comments = [f.readline()[:-1] for i in range(1,header)]
labels = f.readline().split()
try:
N_comment_lines,keyword = f.readline().strip().split(maxsplit=1)
if keyword != 'header':
raise ValueError
else:
comments = [f.readline().strip() for i in range(1,int(N_comment_lines))]
labels = f.readline().split()
except ValueError:
f.seek(0)
comments = []
line = f.readline().strip()
while line.startswith('#'):
comments.append(line.lstrip('#').strip())
line = f.readline().strip()
labels = line.split()
shapes = {}
for label in labels:
@ -95,7 +102,7 @@ class Table():
else:
shapes[label] = (1,)
data = pd.read_csv(f,names=list(range(len(labels))),sep=r'\s+').to_numpy()
data = pd.read_csv(f,names=list(range(len(labels))),sep=r'\s+')
return Table(data,shapes,comments)
@ -309,7 +316,7 @@ class Table():
self.shapes[key] = other.shapes[key]
def to_ASCII(self,fname):
def to_ASCII(self,fname,new=False):
"""
Store as plain text file.
@ -329,15 +336,18 @@ class Table():
for i in range(self.shapes[l][0])]
else:
labels += ['{}:{}_{}'.format('x'.join([str(d) for d in self.shapes[l]]),i+1,l) \
for i in range(np.prod(self.shapes[l],dtype=int))]
for i in range(np.prod(self.shapes[l]))]
header = ['{} header'.format(len(self.comments)+1)] \
+ self.comments \
+ [' '.join(labels)]
if new:
header = ['# {}'.format(comment) for comment in self.comments]
else:
header = ['{} header'.format(len(self.comments)+1)] \
+ self.comments \
try:
f = open(fname,'w')
except TypeError:
f = fname
for line in header: f.write(line+'\n')
for line in header + [' '.join(labels)]: f.write(line+'\n')
self.data.to_csv(f,sep=' ',index=False,header=False)

11
python/tests/test_Table.py
View File

@ -38,14 +38,21 @@ class TestTable:
def test_write_read_str(self,default,tmpdir):
default.to_ASCII(str(tmpdir.join('default.txt')))
new = Table.from_ASCII(str(tmpdir.join('default.txt')))
assert all(default.data==new.data)
assert all(default.data==new.data) and default.shapes == new.shapes
def test_write_read_file(self,default,tmpdir):
with open(tmpdir.join('default.txt'),'w') as f:
default.to_ASCII(f)
with open(tmpdir.join('default.txt')) as f:
new = Table.from_ASCII(f)
assert all(default.data==new.data)
assert all(default.data==new.data) and default.shapes == new.shapes
def test_write_read_new_style(self,default,tmpdir):
with open(tmpdir.join('new_style.txt'),'w') as f:
default.to_ASCII(f,new=True)
with open(tmpdir.join('new_style.txt')) as f:
new = Table.from_ASCII(f)
assert all(default.data==new.data) and default.shapes == new.shapes
def test_read_ang_str(self,reference_dir):
new = Table.from_ang(os.path.join(reference_dir,'simple.ang'))

8
src/CPFEM.f90
View File

@ -86,7 +86,6 @@ subroutine CPFEM_initAll(el,ip)
call config_init
call math_init
call rotations_init
call FE_init
call HDF5_utilities_init
call results_init
call mesh_init(ip, el)
@ -119,7 +118,6 @@ subroutine CPFEM_init
write(6,'(a32,1x,6(i8,1x))') 'CPFEM_cs: ', shape(CPFEM_cs)
write(6,'(a32,1x,6(i8,1x))') 'CPFEM_dcsdE: ', shape(CPFEM_dcsdE)
write(6,'(a32,1x,6(i8,1x),/)') 'CPFEM_dcsdE_knownGood: ', shape(CPFEM_dcsdE_knownGood)
write(6,'(a32,l1)') 'symmetricSolver: ', symmetricSolver
flush(6)
endif
@ -266,10 +264,8 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
!* no parallel computation, so we use just one single elFE and ip for computation
if (.not. parallelExecution) then
FEsolving_execElem(1) = elCP
FEsolving_execElem(2) = elCP
FEsolving_execIP(1,elCP) = ip
FEsolving_execIP(2,elCP) = ip
FEsolving_execElem = elCP
FEsolving_execIP = ip
if (iand(debug_level(debug_CPFEM), debug_levelExtensive) /= 0_pInt) &
write(6,'(a,i8,1x,i2)') '<< CPFEM >> calculation for elFE ip ',elFE,ip
call materialpoint_stressAndItsTangent(updateJaco, dt) ! calculate stress and its tangent

27
src/DAMASK_interface.f90
View File

@ -17,6 +17,7 @@
module DAMASK_interface
use, intrinsic :: iso_fortran_env
use PETScSys
use prec
@ -24,15 +25,15 @@ module DAMASK_interface
implicit none
private
logical, public, protected :: &
logical, public, protected :: &
SIGTERM, & !< termination signal
SIGUSR1, & !< 1. user-defined signal
SIGUSR2 !< 2. user-defined signal
integer, public, protected :: &
integer, public, protected :: &
interface_restartInc = 0 !< Increment at which calculation starts
character(len=1024), public, protected :: &
geometryFile = '', & !< parameter given for geometry file
loadCaseFile = '' !< parameter given for load case file
character(len=:), allocatable, public, protected :: &
geometryFile, & !< parameter given for geometry file
loadCaseFile !< parameter given for load case file
public :: &
getSolverJobName, &
@ -146,7 +147,7 @@ subroutine DAMASK_interface_init
write(6,'(/,a)') ' <<<+- DAMASK_interface init -+>>>'
! http://patorjk.com/software/taag/#p=display&f=Lean&t=DAMASK
! http://patorjk.com/software/taag/#p=display&f=Lean&t=DAMASK
write(6,*) achar(27)//'[94m'
write(6,*) ' _/_/_/ _/_/ _/ _/ _/_/ _/_/_/ _/ _/'
write(6,*) ' _/ _/ _/ _/ _/_/ _/_/ _/ _/ _/ _/ _/'
@ -160,7 +161,7 @@ subroutine DAMASK_interface_init
write(6,'(/,a)') ' Version: '//DAMASKVERSION
! https://github.com/jeffhammond/HPCInfo/blob/master/docs/Preprocessor-Macros.md
! https://github.com/jeffhammond/HPCInfo/blob/master/docs/Preprocessor-Macros.md
#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800
write(6,'(/,a)') ' Compiled with: '//compiler_version()
write(6,'(a)') ' Compiler options: '//compiler_options()
@ -270,8 +271,8 @@ subroutine DAMASK_interface_init
write(6,'(a,a)') ' Geometry argument: ', trim(geometryArg)
write(6,'(a,a)') ' Load case argument: ', trim(loadcaseArg)
write(6,'(a,a)') ' Working directory: ', getCWD()
write(6,'(a,a)') ' Geometry file: ', trim(geometryFile)
write(6,'(a,a)') ' Loadcase file: ', trim(loadCaseFile)
write(6,'(a,a)') ' Geometry file: ', geometryFile
write(6,'(a,a)') ' Loadcase file: ', loadCaseFile
write(6,'(a,a)') ' Solver job name: ', getSolverJobName()
if (interface_restartInc > 0) &
write(6,'(a,i6.6)') ' Restart from increment: ', interface_restartInc
@ -348,9 +349,9 @@ function getGeometryFile(geometryParameter)
getGeometryFile = trim(geometryParameter)
if (scan(getGeometryFile,'/') /= 1) getGeometryFile = getCWD()//'/'//trim(getGeometryFile)
getGeometryFile = makeRelativePath(getCWD(), getGeometryFile)
getGeometryFile = trim(makeRelativePath(getCWD(), getGeometryFile))
inquire(file=trim(getGeometryFile), exist=file_exists)
inquire(file=getGeometryFile, exist=file_exists)
if (.not. file_exists) then
write(6,'(/,a)') ' ERROR: Geometry file does not exists ('//trim(getGeometryFile)//')'
call quit(1)
@ -371,9 +372,9 @@ function getLoadCaseFile(loadCaseParameter)
getLoadCaseFile = trim(loadCaseParameter)
if (scan(getLoadCaseFile,'/') /= 1) getLoadCaseFile = getCWD()//'/'//trim(getLoadCaseFile)
getLoadCaseFile = makeRelativePath(getCWD(), getLoadCaseFile)
getLoadCaseFile = trim(makeRelativePath(getCWD(), getLoadCaseFile))
inquire(file=trim(getLoadCaseFile), exist=file_exists)
inquire(file=getLoadCaseFile, exist=file_exists)
if (.not. file_exists) then
write(6,'(/,a)') ' ERROR: Load case file does not exists ('//trim(getLoadCaseFile)//')'
call quit(1)

163
src/DAMASK_marc.f90
View File

@ -29,23 +29,24 @@
#include "prec.f90"
module DAMASK_interface
use prec
use prec
#if __INTEL_COMPILER >= 1800
use, intrinsic :: iso_fortran_env, only: &
compiler_version, &
compiler_options
use, intrinsic :: iso_fortran_env, only: &
compiler_version, &
compiler_options
#endif
use ifport, only: &
CHDIR
use ifport, only: &
CHDIR
implicit none
private
implicit none
private
character(len=4), parameter, public :: INPUTFILEEXTENSION = '.dat'
logical, public :: symmetricSolver
character(len=*), parameter, public :: INPUTFILEEXTENSION = '.dat'
public :: &
DAMASK_interface_init, &
getSolverJobName
public :: &
DAMASK_interface_init, &
getSolverJobName
contains
@ -54,40 +55,41 @@ contains
!--------------------------------------------------------------------------------------------------
subroutine DAMASK_interface_init
integer, dimension(8) :: &
dateAndTime
integer :: ierr
character(len=1024) :: wd
integer, dimension(8) :: &
dateAndTime
integer :: ierr
character(len=1024) :: wd
write(6,'(/,a)') ' <<<+- DAMASK_marc init -+>>>'
write(6,'(/,a)') ' <<<+- DAMASK_marc init -+>>>'
write(6,'(/,a)') ' Roters et al., Computational Materials Science 158:420–478, 2019'
write(6,'(a)') ' https://doi.org/10.1016/j.commatsci.2018.04.030'
write(6,'(/,a)') ' Roters et al., Computational Materials Science 158:420–478, 2019'
write(6,'(a)') ' https://doi.org/10.1016/j.commatsci.2018.04.030'
write(6,'(/,a)') ' Version: '//DAMASKVERSION
write(6,'(/,a)') ' Version: '//DAMASKVERSION
! https://github.com/jeffhammond/HPCInfo/blob/master/docs/Preprocessor-Macros.md
! https://github.com/jeffhammond/HPCInfo/blob/master/docs/Preprocessor-Macros.md
#if __INTEL_COMPILER >= 1800
write(6,'(/,a)') ' Compiled with: '//compiler_version()
write(6,'(a)') ' Compiler options: '//compiler_options()
write(6,'(/,a)') ' Compiled with: '//compiler_version()
write(6,'(a)') ' Compiler options: '//compiler_options()
#else
write(6,'(/,a,i4.4,a,i8.8)') ' Compiled with Intel fortran version :', __INTEL_COMPILER,&
', build date :', __INTEL_COMPILER_BUILD_DATE
write(6,'(/,a,i4.4,a,i8.8)') ' Compiled with Intel fortran version :', __INTEL_COMPILER,&
', build date :', __INTEL_COMPILER_BUILD_DATE
#endif
write(6,'(/,a)') ' Compiled on: '//__DATE__//' at '//__TIME__
write(6,'(/,a)') ' Compiled on: '//__DATE__//' at '//__TIME__
call date_and_time(values = dateAndTime)
write(6,'(/,a,2(i2.2,a),i4.4)') ' Date: ',dateAndTime(3),'/',dateAndTime(2),'/', dateAndTime(1)
write(6,'(a,2(i2.2,a),i2.2)') ' Time: ',dateAndTime(5),':', dateAndTime(6),':', dateAndTime(7)
call date_and_time(values = dateAndTime)
write(6,'(/,a,2(i2.2,a),i4.4)') ' Date: ',dateAndTime(3),'/',dateAndTime(2),'/', dateAndTime(1)
write(6,'(a,2(i2.2,a),i2.2)') ' Time: ',dateAndTime(5),':', dateAndTime(6),':', dateAndTime(7)
inquire(5, name=wd)
wd = wd(1:scan(wd,'/',back=.true.))
ierr = CHDIR(wd)
if (ierr /= 0) then
write(6,'(a20,a,a16)') ' working directory "',trim(wd),'" does not exist'
call quit(1)
endif
inquire(5, name=wd)
wd = wd(1:scan(wd,'/',back=.true.))
ierr = CHDIR(wd)
if (ierr /= 0) then
write(6,'(a20,a,a16)') ' working directory "',trim(wd),'" does not exist'
call quit(1)
endif
symmetricSolver = solverIsSymmetric()
end subroutine DAMASK_interface_init
@ -97,19 +99,65 @@ end subroutine DAMASK_interface_init
!--------------------------------------------------------------------------------------------------
function getSolverJobName()
character(1024) :: getSolverJobName, inputName
character(len=*), parameter :: pathSep = achar(47)//achar(92) ! forward and backward slash
integer :: extPos
character(len=:), allocatable :: getSolverJobName
character(1024) :: inputName
character(len=*), parameter :: pathSep = achar(47)//achar(92) ! forward and backward slash
integer :: extPos
getSolverJobName=''
inputName=''
inquire(5, name=inputName) ! determine inputfile
extPos = len_trim(inputName)-4
getSolverJobName=inputName(scan(inputName,pathSep,back=.true.)+1:extPos)
inputName=''
inquire(5, name=inputName) ! determine inputfile
extPos = len_trim(inputName)-4
getSolverJobName=inputName(scan(inputName,pathSep,back=.true.)+1:extPos)
end function getSolverJobName
!--------------------------------------------------------------------------------------------------
!> @brief determines whether a symmetric solver is used
!--------------------------------------------------------------------------------------------------
logical function solverIsSymmetric()
character(len=pStringLen) :: line
integer :: myStat,fileUnit,s,e
open(newunit=fileUnit, file=getSolverJobName()//INPUTFILEEXTENSION, &
status='old', position='rewind', action='read',iostat=myStat)
do
read (fileUnit,'(A)',END=100) line
if(index(trim(lc(line)),'solver') == 1) then
read (fileUnit,'(A)',END=100) line ! next line
s = verify(line, ' ') ! start of first chunk
s = s + verify(line(s+1:),' ') ! start of second chunk
e = s + scan (line(s+1:),' ') ! end of second chunk
solverIsSymmetric = line(s:e) /= '1'
endif
enddo
100 close(fileUnit)
contains
!--------------------------------------------------------------------------------------------------
!> @brief changes characters in string to lower case
!> @details copied from IO_lc
!--------------------------------------------------------------------------------------------------
function lc(string)
character(len=*), intent(in) :: string !< string to convert
character(len=len(string)) :: lc
character(26), parameter :: LOWER = 'abcdefghijklmnopqrstuvwxyz'
character(26), parameter :: UPPER = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ'
integer :: i,n
do i=1,len(string)
lc(i:i) = string(i:i)
n = index(UPPER,lc(i:i))
if (n/=0) lc(i:i) = LOWER(n:n)
enddo
end function lc
end function solverIsSymmetric
end module DAMASK_interface
@ -128,6 +176,7 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, &
strechn1,eigvn1,ncrd,itel,ndeg,ndm,nnode, &
jtype,lclass,ifr,ifu)
use prec
use DAMASK_interface
use numerics
use FEsolving
use debug
@ -329,21 +378,21 @@ end subroutine hypela2
!> @brief calculate internal heat generated due to inelastic energy dissipation
!--------------------------------------------------------------------------------------------------
subroutine flux(f,ts,n,time)
use prec
use thermal_conduction
use mesh
use prec
use thermal_conduction
use mesh
implicit none
real(pReal), dimension(6), intent(in) :: &
ts
integer, dimension(10), intent(in) :: &
n
real(pReal), intent(in) :: &
time
real(pReal), dimension(2), intent(out) :: &
f
implicit none
real(pReal), dimension(6), intent(in) :: &
ts
integer, dimension(10), intent(in) :: &
n
real(pReal), intent(in) :: &
time
real(pReal), dimension(2), intent(out) :: &
f
call thermal_conduction_getSourceAndItsTangent(f(1), f(2), ts(3), n(3),mesh_FEasCP('elem',n(1)))
call thermal_conduction_getSourceAndItsTangent(f(1), f(2), ts(3), n(3),mesh_FEasCP('elem',n(1)))
end subroutine flux

56
src/FEsolving.f90
View File

@ -1,65 +1,23 @@
!--------------------------------------------------------------------------------------------------
!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH
!> Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @brief holds some global variables and gets extra information for commercial FEM
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @brief global variables for flow control
!--------------------------------------------------------------------------------------------------
module FEsolving
use prec
use IO
use DAMASK_interface
implicit none
private
logical, public :: &
terminallyIll = .false. !< at least one material point is terminally ill
logical :: &
terminallyIll = .false. !< at least one material point is terminally ill
integer, dimension(:,:), allocatable, public :: &
integer, dimension(2) :: &
FEsolving_execElem, & !< for ping-pong scheme always whole range, otherwise one specific element
FEsolving_execIP !< for ping-pong scheme always range to max IP, otherwise one specific IP
integer, dimension(2), public :: &
FEsolving_execElem !< for ping-pong scheme always whole range, otherwise one specific element
#if defined(Marc4DAMASK) || defined(Abaqus)
logical, public, protected :: &
symmetricSolver = .false. !< use a symmetric FEM solver
logical, dimension(:,:), allocatable, public :: &
logical, dimension(:,:), allocatable :: &
calcMode !< do calculation or simply collect when using ping pong scheme
public :: FE_init
#endif
contains
#if defined(Marc4DAMASK) || defined(Abaqus)
!--------------------------------------------------------------------------------------------------
!> @brief determine whether a symmetric solver is used
!--------------------------------------------------------------------------------------------------
subroutine FE_init
write(6,'(/,a)') ' <<<+- FEsolving init -+>>>'
#if defined(Marc4DAMASK)
block
integer, parameter :: FILEUNIT = 222
character(len=pStringLen) :: line
integer, allocatable, dimension(:) :: chunkPos
call IO_open_inputFile(FILEUNIT)
rewind(FILEUNIT)
do
read (FILEUNIT,'(A)',END=100) line
chunkPos = IO_stringPos(line)
if(IO_lc(IO_stringValue(line,chunkPos,1)) == 'solver') then
read (FILEUNIT,'(A)',END=100) line ! next line
chunkPos = IO_stringPos(line)
symmetricSolver = (IO_intValue(line,chunkPos,2) /= 1)
endif
enddo
100 close(FILEUNIT)
end block
#endif
end subroutine FE_init
#endif
end module FEsolving

87
src/IO.f90
View File

@ -36,14 +36,14 @@ module IO
IO_warning
#if defined(Marc4DAMASK) || defined(Abaqus)
public :: &
IO_open_inputFile, &
IO_continuousIntValues
#endif
IO_open_inputFile
#if defined(Abaqus)
public :: &
IO_continuousIntValues, &
IO_extractValue, &
IO_countDataLines
#endif
#endif
contains
@ -120,7 +120,6 @@ function IO_read_ASCII(fileName) result(fileContent)
fileContent(l) = line
l = l + 1
enddo
end function IO_read_ASCII
@ -398,7 +397,7 @@ end function IO_stringValue
!--------------------------------------------------------------------------------------------------
!> @brief reads float value at myChunk from string
!--------------------------------------------------------------------------------------------------
real(pReal) function IO_floatValue (string,chunkPos,myChunk)
real(pReal) function IO_floatValue(string,chunkPos,myChunk)
integer, dimension(:), intent(in) :: chunkPos !< positions of start and end of each tag/chunk in given string
integer, intent(in) :: myChunk !< position number of desired chunk
@ -452,10 +451,10 @@ pure function IO_lc(string)
character(26), parameter :: LOWER = 'abcdefghijklmnopqrstuvwxyz'
character(26), parameter :: UPPER = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ'
integer :: i,n
integer :: i,n
IO_lc = string
do i=1,len(string)
IO_lc(i:i) = string(i:i)
n = index(UPPER,IO_lc(i:i))
if (n/=0) IO_lc(i:i) = LOWER(n:n)
enddo
@ -792,24 +791,6 @@ subroutine IO_warning(warning_ID,el,ip,g,ext_msg)
end subroutine IO_warning
#if defined(Abaqus) || defined(Marc4DAMASK)
!--------------------------------------------------------------------------------------------------
!> @brief reads a line from a text file.
!--------------------------------------------------------------------------------------------------
function IO_read(fileUnit) result(line)
integer, intent(in) :: fileUnit !< file unit
character(len=pStringLen) :: line
read(fileUnit,'(a256)',END=100) line
100 end function IO_read
#ifdef Abaqus
!--------------------------------------------------------------------------------------------------
!> @brief extracts string value from key=value pair and check whether key matches
@ -847,7 +828,7 @@ integer function IO_countDataLines(fileUnit)
line = ''
do while (trim(line) /= IO_EOF)
line = IO_read(fileUnit)
read(fileUnit,'(A)') line
chunkPos = IO_stringPos(line)
tmp = IO_lc(IO_stringValue(line,chunkPos,1))
if (tmp(1:1) == '*' .and. tmp(2:2) /= '*') then ! found keyword
@ -859,14 +840,12 @@ integer function IO_countDataLines(fileUnit)
backspace(fileUnit)
end function IO_countDataLines
#endif
!--------------------------------------------------------------------------------------------------
!> @brief return integer list corresponding to items in consecutive lines.
!! First integer in array is counter
!> @details Marc: ints concatenated by "c" as last char, range of a "to" b, or named set
!! Abaqus: triplet of start,stop,inc or named set
!> @details Abaqus: triplet of start,stop,inc or named set
!--------------------------------------------------------------------------------------------------
function IO_continuousIntValues(fileUnit,maxN,lookupName,lookupMap,lookupMaxN)
@ -878,9 +857,7 @@ function IO_continuousIntValues(fileUnit,maxN,lookupName,lookupMap,lookupMaxN)
integer, dimension(:,:), intent(in) :: lookupMap
character(len=*), dimension(:), intent(in) :: lookupName
integer :: i,first,last
#ifdef Abaqus
integer :: j,l,c
#endif
integer, allocatable, dimension(:) :: chunkPos
character(len=pStringLen) :: line
logical :: rangeGeneration
@ -888,41 +865,6 @@ function IO_continuousIntValues(fileUnit,maxN,lookupName,lookupMap,lookupMaxN)
IO_continuousIntValues = 0
rangeGeneration = .false.
#if defined(Marc4DAMASK)
do
read(fileUnit,'(A)',end=100) line
chunkPos = IO_stringPos(line)
if (chunkPos(1) < 1) then ! empty line
exit
elseif (verify(IO_stringValue(line,chunkPos,1),'0123456789') > 0) then ! a non-int, i.e. set name
do i = 1, lookupMaxN ! loop over known set names
if (IO_stringValue(line,chunkPos,1) == lookupName(i)) then ! found matching name
IO_continuousIntValues = lookupMap(:,i) ! return resp. entity list
exit
endif
enddo
exit
else if (chunkPos(1) > 2 .and. IO_lc(IO_stringValue(line,chunkPos,2)) == 'to' ) then ! found range indicator
first = IO_intValue(line,chunkPos,1)
last = IO_intValue(line,chunkPos,3)
do i = first, last, sign(1,last-first)
IO_continuousIntValues(1) = IO_continuousIntValues(1) + 1
IO_continuousIntValues(1+IO_continuousIntValues(1)) = i
enddo
exit
else
do i = 1,chunkPos(1)-1 ! interpret up to second to last value
IO_continuousIntValues(1) = IO_continuousIntValues(1) + 1
IO_continuousIntValues(1+IO_continuousIntValues(1)) = IO_intValue(line,chunkPos,i)
enddo
if ( IO_lc(IO_stringValue(line,chunkPos,chunkPos(1))) /= 'c' ) then ! line finished, read last value
IO_continuousIntValues(1) = IO_continuousIntValues(1) + 1
IO_continuousIntValues(1+IO_continuousIntValues(1)) = IO_intValue(line,chunkPos,chunkPos(1))
exit
endif
endif
enddo
#elif defined(Abaqus)
c = IO_countDataLines(fileUnit)
do l = 1,c
backspace(fileUnit)
@ -965,7 +907,6 @@ function IO_continuousIntValues(fileUnit,maxN,lookupName,lookupMap,lookupMaxN)
enddo
endif
enddo
#endif
100 end function IO_continuousIntValues
#endif
@ -976,7 +917,7 @@ function IO_continuousIntValues(fileUnit,maxN,lookupName,lookupMap,lookupMaxN)
!--------------------------------------------------------------------------------------------------
!> @brief returns verified integer value in given string
!--------------------------------------------------------------------------------------------------
integer function verifyIntValue (string,validChars,myName)
integer function verifyIntValue(string,validChars,myName)
character(len=*), intent(in) :: string, & !< string for conversion to int value. Must not contain spaces!
validChars, & !< valid characters in string
@ -987,12 +928,12 @@ integer function verifyIntValue (string,validChars,myName)
invalidWhere = verify(string,validChars)
if (invalidWhere == 0) then
read(UNIT=string,iostat=readStatus,FMT=*) verifyIntValue ! no offending chars found
read(string,*,iostat=readStatus) verifyIntValue ! no offending chars found
if (readStatus /= 0) & ! error during string to integer conversion
call IO_warning(203,ext_msg=myName//'"'//string//'"')
else
call IO_warning(202,ext_msg=myName//'"'//string//'"') ! complain about offending characters
read(UNIT=string(1:invalidWhere-1),iostat=readStatus,FMT=*) verifyIntValue ! interpret remaining string
read(string(1:invalidWhere-1),*,iostat=readStatus) verifyIntValue ! interpret remaining string
if (readStatus /= 0) & ! error during string to integer conversion
call IO_warning(203,ext_msg=myName//'"'//string(1:invalidWhere-1)//'"')
endif
@ -1003,7 +944,7 @@ end function verifyIntValue
!--------------------------------------------------------------------------------------------------
!> @brief returns verified float value in given string
!--------------------------------------------------------------------------------------------------
real(pReal) function verifyFloatValue (string,validChars,myName)
real(pReal) function verifyFloatValue(string,validChars,myName)
character(len=*), intent(in) :: string, & !< string for conversion to int value. Must not contain spaces!
validChars, & !< valid characters in string
@ -1015,12 +956,12 @@ real(pReal) function verifyFloatValue (string,validChars,myName)
invalidWhere = verify(string,validChars)
if (invalidWhere == 0) then
read(UNIT=string,iostat=readStatus,FMT=*) verifyFloatValue ! no offending chars found
read(string,*,iostat=readStatus) verifyFloatValue ! no offending chars found
if (readStatus /= 0) & ! error during string to float conversion
call IO_warning(203,ext_msg=myName//'"'//string//'"')
else
call IO_warning(202,ext_msg=myName//'"'//string//'"') ! complain about offending characters
read(UNIT=string(1:invalidWhere-1),iostat=readStatus,FMT=*) verifyFloatValue ! interpret remaining string
read(string(1:invalidWhere-1),*,iostat=readStatus) verifyFloatValue ! interpret remaining string
if (readStatus /= 0) & ! error during string to float conversion
call IO_warning(203,ext_msg=myName//'"'//string(1:invalidWhere-1)//'"')
endif

56
src/crystallite.f90
View File

@ -238,7 +238,7 @@ subroutine crystallite_init
!$OMP PARALLEL DO PRIVATE(myNcomponents,i,c)
do e = FEsolving_execElem(1),FEsolving_execElem(2)
myNcomponents = homogenization_Ngrains(material_homogenizationAt(e))
do i = FEsolving_execIP(1,e), FEsolving_execIP(2,e); do c = 1, myNcomponents
do i = FEsolving_execIP(1), FEsolving_execIP(2); do c = 1, myNcomponents
crystallite_Fp0(1:3,1:3,c,i,e) = math_EulerToR(material_Eulers(1:3,c,i,e)) ! plastic def gradient reflects init orientation
crystallite_Fi0(1:3,1:3,c,i,e) = constitutive_initialFi(c,i,e)
crystallite_F0(1:3,1:3,c,i,e) = math_I3
@ -263,7 +263,7 @@ subroutine crystallite_init
!$OMP PARALLEL DO
do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
do i = FEsolving_execIP(1),FEsolving_execIP(2)
do c = 1,homogenization_Ngrains(material_homogenizationAt(e))
call constitutive_microstructure(crystallite_Fe(1:3,1:3,c,i,e), &
crystallite_Fp(1:3,1:3,c,i,e), &
@ -336,7 +336,7 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC)
crystallite_subStep = 0.0_pReal
!$OMP PARALLEL DO
elementLooping1: do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e); do c = 1,homogenization_Ngrains(material_homogenizationAt(e))
do i = FEsolving_execIP(1),FEsolving_execIP(2); do c = 1,homogenization_Ngrains(material_homogenizationAt(e))
homogenizationRequestsCalculation: if (crystallite_requested(c,i,e)) then
plasticState (material_phaseAt(c,e))%subState0( :,material_phaseMemberAt(c,i,e)) = &
plasticState (material_phaseAt(c,e))%partionedState0(:,material_phaseMemberAt(c,i,e))
@ -361,12 +361,12 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC)
!$OMP END PARALLEL DO
singleRun: if (FEsolving_execELem(1) == FEsolving_execElem(2) .and. &
FEsolving_execIP(1,FEsolving_execELem(1))==FEsolving_execIP(2,FEsolving_execELem(1))) then
startIP = FEsolving_execIP(1,FEsolving_execELem(1))
FEsolving_execIP (1) == FEsolving_execIP (2)) then
startIP = FEsolving_execIP(1)
endIP = startIP
else singleRun
startIP = 1
endIP = discretization_nIP
endIP = discretization_nIP
endif singleRun
NiterationCrystallite = 0
@ -379,7 +379,7 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC)
#endif
!$OMP PARALLEL DO PRIVATE(formerSubStep)
elementLooping3: do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
do i = FEsolving_execIP(1),FEsolving_execIP(2)
do c = 1,homogenization_Ngrains(material_homogenizationAt(e))
!--------------------------------------------------------------------------------------------------
! wind forward
@ -494,14 +494,14 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC)
! return whether converged or not
crystallite_stress = .false.
elementLooping5: do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
do i = FEsolving_execIP(1),FEsolving_execIP(2)
crystallite_stress(i,e) = all(crystallite_converged(:,i,e))
enddo
enddo elementLooping5
#ifdef DEBUG
elementLooping6: do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
do i = FEsolving_execIP(1),FEsolving_execIP(2)
do c = 1,homogenization_Ngrains(material_homogenizationAt(e))
if (.not. crystallite_converged(c,i,e)) then
if(iand(debug_level(debug_crystallite), debug_levelBasic) /= 0) &
@ -563,7 +563,7 @@ subroutine crystallite_stressTangent
!$OMP PARALLEL DO PRIVATE(dSdF,dSdFe,dSdFi,dLpdS,dLpdFi,dFpinvdF,dLidS,dLidFi,dFidS,invSubFi0,o,p, &
!$OMP rhs_3333,lhs_3333,temp_99,temp_33_1,temp_33_2,temp_33_3,temp_33_4,temp_3333,error)
elementLooping: do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
do i = FEsolving_execIP(1),FEsolving_execIP(2)
do c = 1,homogenization_Ngrains(material_homogenizationAt(e))
call constitutive_SandItsTangents(devNull,dSdFe,dSdFi, &
@ -684,7 +684,7 @@ subroutine crystallite_orientations
!$OMP PARALLEL DO
do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
do i = FEsolving_execIP(1),FEsolving_execIP(2)
do c = 1,homogenization_Ngrains(material_homogenizationAt(e))
call crystallite_orientation(c,i,e)%fromMatrix(transpose(math_rotationalPart33(crystallite_Fe(1:3,1:3,c,i,e))))
enddo; enddo; enddo
@ -693,7 +693,7 @@ subroutine crystallite_orientations
nonlocalPresent: if (any(plasticState%nonLocal)) then
!$OMP PARALLEL DO
do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
do i = FEsolving_execIP(1),FEsolving_execIP(2)
if (plasticState(material_phaseAt(1,e))%nonLocal) & ! if nonlocal model
call plastic_nonlocal_updateCompatibility(crystallite_orientation,i,e)
enddo; enddo
@ -1306,7 +1306,7 @@ subroutine integrateStateFPI
!$OMP PARALLEL DO PRIVATE(p,c)
do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
do i = FEsolving_execIP(1),FEsolving_execIP(2)
do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then
p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e)
@ -1334,7 +1334,7 @@ subroutine integrateStateFPI
!$OMP PARALLEL
!$OMP DO PRIVATE(sizeDotState,residuum_plastic,residuum_source,zeta,p,c)
do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
do i = FEsolving_execIP(1),FEsolving_execIP(2)
do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then
p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e)
@ -1389,7 +1389,7 @@ subroutine integrateStateFPI
!$OMP DO
do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
do i = FEsolving_execIP(1),FEsolving_execIP(2)
do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
!$OMP FLUSH(crystallite_todo)
if (crystallite_todo(g,i,e) .and. crystallite_converged(g,i,e)) then ! converged and still alive...
@ -1415,7 +1415,7 @@ subroutine integrateStateFPI
! --- CHECK IF DONE WITH INTEGRATION ---
doneWithIntegration = .true.
do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
do i = FEsolving_execIP(1),FEsolving_execIP(2)
do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then
doneWithIntegration = .false.
@ -1497,7 +1497,7 @@ subroutine integrateStateAdaptiveEuler
!$OMP PARALLEL DO PRIVATE(sizeDotState,p,c)
do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
do i = FEsolving_execIP(1),FEsolving_execIP(2)
do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
if (crystallite_todo(g,i,e)) then
p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e)
@ -1526,7 +1526,7 @@ subroutine integrateStateAdaptiveEuler
!$OMP PARALLEL DO PRIVATE(sizeDotState,p,c)
do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
do i = FEsolving_execIP(1),FEsolving_execIP(2)
do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
if (crystallite_todo(g,i,e)) then
p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e)
@ -1586,7 +1586,7 @@ subroutine integrateStateRK4
!$OMP PARALLEL DO PRIVATE(p,c)
do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
do i = FEsolving_execIP(1),FEsolving_execIP(2)
do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
if (crystallite_todo(g,i,e)) then
p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e)
@ -1677,7 +1677,7 @@ subroutine integrateStateRKCK45
!$OMP PARALLEL DO PRIVATE(p,cc)
do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
do i = FEsolving_execIP(1),FEsolving_execIP(2)
do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
if (crystallite_todo(g,i,e)) then
p = material_phaseAt(g,e); cc = material_phaseMemberAt(g,i,e)
@ -1717,7 +1717,7 @@ subroutine integrateStateRKCK45
!$OMP PARALLEL DO PRIVATE(sizeDotState,p,cc)
do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
do i = FEsolving_execIP(1),FEsolving_execIP(2)
do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
if (crystallite_todo(g,i,e)) then
p = material_phaseAt(g,e); cc = material_phaseMemberAt(g,i,e)
@ -1756,7 +1756,7 @@ subroutine integrateStateRKCK45
!$OMP PARALLEL DO PRIVATE(sizeDotState,p,cc)
do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
do i = FEsolving_execIP(1),FEsolving_execIP(2)
do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
if (crystallite_todo(g,i,e)) then
p = material_phaseAt(g,e); cc = material_phaseMemberAt(g,i,e)
@ -1814,7 +1814,7 @@ subroutine setConvergenceFlag
!OMP DO PARALLEL PRIVATE
do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
do i = FEsolving_execIP(1),FEsolving_execIP(2)
do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
crystallite_converged(g,i,e) = crystallite_todo(g,i,e) .or. crystallite_converged(g,i,e) ! if still "to do" then converged per definition
enddo; enddo; enddo
@ -1854,7 +1854,7 @@ subroutine update_stress(timeFraction)
!$OMP PARALLEL DO
do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
do i = FEsolving_execIP(1),FEsolving_execIP(2)
do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
!$OMP FLUSH(crystallite_todo)
if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then
@ -1884,7 +1884,7 @@ subroutine update_dependentState
!$OMP PARALLEL DO
do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
do i = FEsolving_execIP(1),FEsolving_execIP(2)
do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) &
call constitutive_dependentState(crystallite_Fe(1:3,1:3,g,i,e), &
@ -1914,7 +1914,7 @@ subroutine update_state(timeFraction)
!$OMP PARALLEL DO PRIVATE(mySize,p,c)
do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
do i = FEsolving_execIP(1),FEsolving_execIP(2)
do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then
p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e)
@ -1959,7 +1959,7 @@ subroutine update_dotState(timeFraction)
!$OMP PARALLEL DO PRIVATE (p,c,NaN)
do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
do i = FEsolving_execIP(1),FEsolving_execIP(2)
do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
!$OMP FLUSH(nonlocalStop)
if ((crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) .and. .not. nonlocalStop) then
@ -2005,7 +2005,7 @@ subroutine update_deltaState
!$OMP PARALLEL DO PRIVATE(p,c,myOffset,mySize,NaN)
do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
do i = FEsolving_execIP(1),FEsolving_execIP(2)
do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
!$OMP FLUSH(nonlocalStop)
if ((crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) .and. .not. nonlocalStop) then

4
src/damage_local.f90
View File

@ -74,7 +74,7 @@ subroutine damage_local_init
allocate(damageState(h)%state (1,NofMyHomog), source=damage_initialPhi(h))
nullify(damageMapping(h)%p)
damageMapping(h)%p => mappingHomogenization(1,:,:)
damageMapping(h)%p => material_homogenizationMemberAt
deallocate(damage(h)%p)
damage(h)%p => damageState(h)%state(1,:)
@ -103,7 +103,7 @@ function damage_local_updateState(subdt, ip, el)
phi, phiDot, dPhiDot_dPhi
homog = material_homogenizationAt(el)
offset = mappingHomogenization(1,ip,el)
offset = material_homogenizationMemberAt(ip,el)
phi = damageState(homog)%subState0(1,offset)
call damage_local_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el)
phi = max(residualStiffness,min(1.0_pReal,phi + subdt*phiDot))

2
src/damage_nonlocal.f90
View File

@ -79,7 +79,7 @@ subroutine damage_nonlocal_init
allocate(damageState(h)%state (1,NofMyHomog), source=damage_initialPhi(h))
nullify(damageMapping(h)%p)
damageMapping(h)%p => mappingHomogenization(1,:,:)
damageMapping(h)%p => material_homogenizationMemberAt
deallocate(damage(h)%p)
damage(h)%p => damageState(h)%state(1,:)

20
src/grid/DAMASK_grid.f90
View File

@ -65,12 +65,11 @@ program DAMASK_spectral
currentLoadcase = 0, & !< current load case
inc, & !< current increment in current load case
totalIncsCounter = 0, & !< total # of increments
fileUnit = 0, & !< file unit for reading load case and writing results
myStat, &
statUnit = 0, & !< file unit for statistics output
stagIter, &
nActiveFields = 0
character(len=6) :: loadcase_string
character(len=pStringLen), dimension(:), allocatable :: fileContent
character(len=1024) :: &
incInfo
type(tLoadCase), allocatable, dimension(:) :: loadCases !< array of all load cases
@ -141,17 +140,14 @@ program DAMASK_spectral
!--------------------------------------------------------------------------------------------------
! reading information from load case file and to sanity checks
fileContent = IO_read_ASCII(trim(loadCaseFile))
allocate (loadCases(0)) ! array of load cases
open(newunit=fileunit,iostat=myStat,file=trim(loadCaseFile),action='read')
if (myStat /= 0) call IO_error(100,el=myStat,ext_msg=trim(loadCaseFile))
do
read(fileUnit, '(A)', iostat=myStat) line
if ( myStat /= 0) exit
if (IO_isBlank(line)) cycle ! skip empty lines
currentLoadCase = currentLoadCase + 1
do currentLoadCase = 1, size(fileContent)
line = fileContent(currentLoadCase)
if (IO_isBlank(line)) cycle
chunkPos = IO_stringPos(line)
do i = 1, chunkPos(1) ! reading compulsory parameters for loadcase
select case (IO_lc(IO_stringValue(line,chunkPos,i)))
case('l','fdot','dotf','f')
@ -300,7 +296,7 @@ program DAMASK_spectral
endif reportAndCheck
loadCases = [loadCases,newLoadCase] ! load case is ok, append it
enddo
close(fileUnit)
!--------------------------------------------------------------------------------------------------
! doing initialization depending on active solvers

46
src/homogenization.f90
View File

@ -215,7 +215,7 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
! initialize restoration points of ...
do e = FEsolving_execElem(1),FEsolving_execElem(2)
myNgrains = homogenization_Ngrains(material_homogenizationAt(e))
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e);
do i = FEsolving_execIP(1),FEsolving_execIP(2);
do g = 1,myNgrains
plasticState (material_phaseAt(g,e))%partionedState0(:,material_phasememberAt(g,i,e)) = &
@ -242,16 +242,16 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
materialpoint_requested(i,e) = .true. ! everybody requires calculation
if (homogState(material_homogenizationAt(e))%sizeState > 0) &
homogState(material_homogenizationAt(e))%subState0(:,mappingHomogenization(1,i,e)) = &
homogState(material_homogenizationAt(e))%State0( :,mappingHomogenization(1,i,e)) ! ...internal homogenization state
homogState(material_homogenizationAt(e))%subState0(:,material_homogenizationMemberAt(i,e)) = &
homogState(material_homogenizationAt(e))%State0( :,material_homogenizationMemberAt(i,e)) ! ...internal homogenization state
if (thermalState(material_homogenizationAt(e))%sizeState > 0) &
thermalState(material_homogenizationAt(e))%subState0(:,mappingHomogenization(1,i,e)) = &
thermalState(material_homogenizationAt(e))%State0( :,mappingHomogenization(1,i,e)) ! ...internal thermal state
thermalState(material_homogenizationAt(e))%subState0(:,material_homogenizationMemberAt(i,e)) = &
thermalState(material_homogenizationAt(e))%State0( :,material_homogenizationMemberAt(i,e)) ! ...internal thermal state
if (damageState(material_homogenizationAt(e))%sizeState > 0) &
damageState(material_homogenizationAt(e))%subState0(:,mappingHomogenization(1,i,e)) = &
damageState(material_homogenizationAt(e))%State0( :,mappingHomogenization(1,i,e)) ! ...internal damage state
damageState(material_homogenizationAt(e))%subState0(:,material_homogenizationMemberAt(i,e)) = &
damageState(material_homogenizationAt(e))%State0( :,material_homogenizationMemberAt(i,e)) ! ...internal damage state
enddo
enddo
@ -263,7 +263,7 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
!$OMP PARALLEL DO PRIVATE(myNgrains)
elementLooping1: do e = FEsolving_execElem(1),FEsolving_execElem(2)
myNgrains = homogenization_Ngrains(material_homogenizationAt(e))
IpLooping1: do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
IpLooping1: do i = FEsolving_execIP(1),FEsolving_execIP(2)
converged: if (materialpoint_converged(i,e)) then
#ifdef DEBUG
@ -313,14 +313,14 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
enddo
if(homogState(material_homogenizationAt(e))%sizeState > 0) &
homogState(material_homogenizationAt(e))%subState0(:,mappingHomogenization(1,i,e)) = &
homogState(material_homogenizationAt(e))%State (:,mappingHomogenization(1,i,e))
homogState(material_homogenizationAt(e))%subState0(:,material_homogenizationMemberAt(i,e)) = &
homogState(material_homogenizationAt(e))%State (:,material_homogenizationMemberAt(i,e))
if(thermalState(material_homogenizationAt(e))%sizeState > 0) &
thermalState(material_homogenizationAt(e))%subState0(:,mappingHomogenization(1,i,e)) = &
thermalState(material_homogenizationAt(e))%State (:,mappingHomogenization(1,i,e))
thermalState(material_homogenizationAt(e))%subState0(:,material_homogenizationMemberAt(i,e)) = &
thermalState(material_homogenizationAt(e))%State (:,material_homogenizationMemberAt(i,e))
if(damageState(material_homogenizationAt(e))%sizeState > 0) &
damageState(material_homogenizationAt(e))%subState0(:,mappingHomogenization(1,i,e)) = &
damageState(material_homogenizationAt(e))%State (:,mappingHomogenization(1,i,e))
damageState(material_homogenizationAt(e))%subState0(:,material_homogenizationMemberAt(i,e)) = &
damageState(material_homogenizationAt(e))%State (:,material_homogenizationMemberAt(i,e))
materialpoint_subF0(1:3,1:3,i,e) = materialpoint_subF(1:3,1:3,i,e)
@ -375,14 +375,14 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
enddo
enddo
if(homogState(material_homogenizationAt(e))%sizeState > 0) &
homogState(material_homogenizationAt(e))%State( :,mappingHomogenization(1,i,e)) = &
homogState(material_homogenizationAt(e))%subState0(:,mappingHomogenization(1,i,e))
homogState(material_homogenizationAt(e))%State( :,material_homogenizationMemberAt(i,e)) = &
homogState(material_homogenizationAt(e))%subState0(:,material_homogenizationMemberAt(i,e))
if(thermalState(material_homogenizationAt(e))%sizeState > 0) &
thermalState(material_homogenizationAt(e))%State( :,mappingHomogenization(1,i,e)) = &
thermalState(material_homogenizationAt(e))%subState0(:,mappingHomogenization(1,i,e))
thermalState(material_homogenizationAt(e))%State( :,material_homogenizationMemberAt(i,e)) = &
thermalState(material_homogenizationAt(e))%subState0(:,material_homogenizationMemberAt(i,e))
if(damageState(material_homogenizationAt(e))%sizeState > 0) &
damageState(material_homogenizationAt(e))%State( :,mappingHomogenization(1,i,e)) = &
damageState(material_homogenizationAt(e))%subState0(:,mappingHomogenization(1,i,e))
damageState(material_homogenizationAt(e))%State( :,material_homogenizationMemberAt(i,e)) = &
damageState(material_homogenizationAt(e))%subState0(:,material_homogenizationMemberAt(i,e))
endif
endif converged
@ -414,7 +414,7 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
!$OMP PARALLEL DO PRIVATE(myNgrains)
elementLooping2: do e = FEsolving_execElem(1),FEsolving_execElem(2)
myNgrains = homogenization_Ngrains(material_homogenizationAt(e))
IpLooping2: do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
IpLooping2: do i = FEsolving_execIP(1),FEsolving_execIP(2)
if ( materialpoint_requested(i,e) .and. & ! process requested but...
.not. materialpoint_doneAndHappy(1,i,e)) then ! ...not yet done material points
call partitionDeformation(i,e) ! partition deformation onto constituents
@ -438,7 +438,7 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
! state update
!$OMP PARALLEL DO
elementLooping3: do e = FEsolving_execElem(1),FEsolving_execElem(2)
IpLooping3: do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
IpLooping3: do i = FEsolving_execIP(1),FEsolving_execIP(2)
if ( materialpoint_requested(i,e) .and. &
.not. materialpoint_doneAndHappy(1,i,e)) then
if (.not. materialpoint_converged(i,e)) then
@ -464,7 +464,7 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
call crystallite_orientations() ! calculate crystal orientations
!$OMP PARALLEL DO
elementLooping4: do e = FEsolving_execElem(1),FEsolving_execElem(2)
IpLooping4: do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
IpLooping4: do i = FEsolving_execIP(1),FEsolving_execIP(2)
call averageStressAndItsTangent(i,e)
enddo IpLooping4
enddo elementLooping4

4
src/homogenization_mech_RGC.f90
View File

@ -108,7 +108,7 @@ module subroutine mech_RGC_init
#ifdef DEBUG
if (h==material_homogenizationAt(debug_e)) then
prm%of_debug = mappingHomogenization(1,debug_i,debug_e)
prm%of_debug = material_homogenizationMemberAt(debug_i,debug_e)
endif
#endif
@ -261,7 +261,7 @@ module procedure mech_RGC_updateState
endif zeroTimeStep
instance = homogenization_typeInstance(material_homogenizationAt(el))
of = mappingHomogenization(1,ip,el)
of = material_homogenizationMemberAt(ip,el)
associate(stt => state(instance), st0 => state0(instance), dst => dependentState(instance), prm => param(instance))

52
src/material.f90
View File

@ -18,7 +18,7 @@ module material
implicit none
private
character(len=*), parameter, public :: &
character(len=*), parameter, public :: &
ELASTICITY_hooke_label = 'hooke', &
PLASTICITY_none_label = 'none', &
PLASTICITY_isotropic_label = 'isotropic', &
@ -93,6 +93,8 @@ module material
thermal_type !< thermal transport model
integer(kind(DAMAGE_none_ID)), dimension(:), allocatable, public, protected :: &
damage_type !< nonlocal damage model
integer(kind(HOMOGENIZATION_undefined_ID)), dimension(:), allocatable, public, protected :: &
homogenization_type !< type of each homogenization
integer, public, protected :: &
material_Nphase, & !< number of phases
@ -103,9 +105,6 @@ module material
phase_kinematics, & !< active kinematic mechanisms of each phase
phase_stiffnessDegradation !< active stiffness degradation mechanisms of each phase
integer(kind(HOMOGENIZATION_undefined_ID)), dimension(:), allocatable, public, protected :: &
homogenization_type !< type of each homogenization
integer, public, protected :: &
homogenization_maxNgrains !< max number of grains in any USED homogenization
@ -113,11 +112,9 @@ module material
phase_Nsources, & !< number of source mechanisms active in each phase
phase_Nkinematics, & !< number of kinematic mechanisms active in each phase
phase_NstiffnessDegradations, & !< number of stiffness degradation mechanisms active in each phase
phase_Noutput, & !< number of '(output)' items per phase
phase_elasticityInstance, & !< instance of particular elasticity of each phase
phase_plasticityInstance, & !< instance of particular plasticity of each phase
homogenization_Ngrains, & !< number of grains in each homogenization
homogenization_Noutput, & !< number of '(output)' items per homogenization
homogenization_typeInstance, & !< instance of particular type of each homogenization
thermal_typeInstance, & !< instance of particular type of each thermal transport
damage_typeInstance !< instance of particular type of each nonlocal damage
@ -129,7 +126,7 @@ module material
! NEW MAPPINGS
integer, dimension(:), allocatable, public, protected :: & ! (elem)
material_homogenizationAt !< homogenization ID of each element (copy of discretization_homogenizationAt)
integer, dimension(:,:), allocatable, public, protected :: & ! (ip,elem)
integer, dimension(:,:), allocatable, public, target :: & ! (ip,elem) ToDo: ugly target for mapping hack
material_homogenizationMemberAt !< position of the element within its homogenization instance
integer, dimension(:,:), allocatable, public, protected :: & ! (constituent,elem)
material_phaseAt !< phase ID of each element
@ -153,12 +150,8 @@ module material
material_orientation0 !< initial orientation of each grain,IP,element
logical, dimension(:), allocatable, public, protected :: &
microstructure_active, &
phase_localPlasticity !< flags phases with local constitutive law
integer, private :: &
microstructure_maxNconstituents !< max number of constituents in any phase
integer, dimension(:), allocatable, private :: &
microstructure_Nconstituents !< number of constituents in each microstructure
@ -170,14 +163,9 @@ module material
material_Eulers
type(Rotation), dimension(:), allocatable, private :: &
texture_orientation !< Euler angles in material.config (possibly rotated for alignment)
real(pReal), dimension(:,:), allocatable, private :: &
microstructure_fraction !< vol fraction of each constituent in microstructure
logical, dimension(:), allocatable, private :: &
homogenization_active
! BEGIN DEPRECATED
integer, dimension(:,:,:), allocatable, public, target :: mappingHomogenization !< mapping from material points to offset in heterogenous state/field
integer, dimension(:,:), allocatable, private, target :: mappingHomogenizationConst !< mapping from material points to offset in constant state/field
! END DEPRECATED
@ -223,6 +211,7 @@ module material
HOMOGENIZATION_RGC_ID
contains
!--------------------------------------------------------------------------------------------------
!> @brief parses material configuration file
!--------------------------------------------------------------------------------------------------
@ -235,7 +224,7 @@ subroutine material_init
myDebug = debug_level(debug_material)
write(6,'(/,a)') ' <<<+- material init -+>>>'
write(6,'(/,a)') ' <<<+- material init -+>>>'; flush(6)
call material_parsePhase()
if (iand(myDebug,debug_levelBasic) /= 0) write(6,'(a)') ' Phase parsed'; flush(6)
@ -294,9 +283,8 @@ subroutine material_init
microstructure_Nconstituents(m)
if (microstructure_Nconstituents(m) > 0) then
do c = 1,microstructure_Nconstituents(m)
write(6,'(a1,1x,a32,1x,a32,1x,f7.4)') '>',config_name_phase(microstructure_phase(c,m)),&
config_name_texture(microstructure_texture(c,m)),&
microstructure_fraction(c,m)
write(6,'(a1,1x,a32,1x,a32)') '>',config_name_phase(microstructure_phase(c,m)),&
config_name_texture(microstructure_texture(c,m))
enddo
write(6,*)
endif
@ -362,18 +350,8 @@ subroutine material_init
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! BEGIN DEPRECATED
allocate(mappingHomogenization (2,discretization_nIP,discretization_nElem),source=0)
allocate(mappingHomogenizationConst( discretization_nIP,discretization_nElem),source=1)
CounterHomogenization=0
do e = 1,discretization_nElem
myHomog = discretization_homogenizationAt(e)
do i = 1, discretization_nIP
CounterHomogenization(myHomog) = CounterHomogenization(myHomog) + 1
mappingHomogenization(1:2,i,e) = [CounterHomogenization(myHomog),huge(1)]
enddo
enddo
! hack needed to initialize field values used during constitutive and crystallite initializations
do myHomog = 1,size(config_homogenization)
thermalMapping (myHomog)%p => mappingHomogenizationConst
@ -394,6 +372,8 @@ subroutine material_parseHomogenization
integer :: h
character(len=pStringLen) :: tag
logical, dimension(:), allocatable :: homogenization_active
allocate(homogenization_type(size(config_homogenization)), source=HOMOGENIZATION_undefined_ID)
allocate(thermal_type(size(config_homogenization)), source=THERMAL_isothermal_ID)
allocate(damage_type (size(config_homogenization)), source=DAMAGE_none_ID)
@ -401,7 +381,6 @@ subroutine material_parseHomogenization
allocate(thermal_typeInstance(size(config_homogenization)), source=0)
allocate(damage_typeInstance(size(config_homogenization)), source=0)
allocate(homogenization_Ngrains(size(config_homogenization)), source=0)
allocate(homogenization_Noutput(size(config_homogenization)), source=0)
allocate(homogenization_active(size(config_homogenization)), source=.false.) !!!!!!!!!!!!!!!
allocate(thermal_initialT(size(config_homogenization)), source=300.0_pReal)
allocate(damage_initialPhi(size(config_homogenization)), source=1.0_pReal)
@ -411,7 +390,6 @@ subroutine material_parseHomogenization
do h=1, size(config_homogenization)
homogenization_Noutput(h) = config_homogenization(h)%countKeys('(output)')
tag = config_homogenization(h)%getString('mech')
select case (trim(tag))
@ -488,16 +466,16 @@ subroutine material_parseMicrostructure
integer :: e, m, c, i
character(len=pStringLen) :: &
tag
real(pReal), dimension(:,:), allocatable :: &
microstructure_fraction !< vol fraction of each constituent in microstructure
integer :: &
microstructure_maxNconstituents !< max number of constituents in any phase
allocate(microstructure_Nconstituents(size(config_microstructure)), source=0)
allocate(microstructure_active(size(config_microstructure)), source=.false.)
if(any(discretization_microstructureAt > size(config_microstructure))) &
call IO_error(155,ext_msg='More microstructures in geometry than sections in material.config')
forall (e = 1:discretization_nElem) &
microstructure_active(discretization_microstructureAt(e)) = .true. ! current microstructure used in model? Elementwise view, maximum N operations for N elements
do m=1, size(config_microstructure)
microstructure_Nconstituents(m) = config_microstructure(m)%countKeys('(constituent)')
enddo
@ -548,11 +526,9 @@ subroutine material_parsePhase
allocate(phase_Nsources(size(config_phase)), source=0)
allocate(phase_Nkinematics(size(config_phase)), source=0)
allocate(phase_NstiffnessDegradations(size(config_phase)),source=0)
allocate(phase_Noutput(size(config_phase)), source=0)
allocate(phase_localPlasticity(size(config_phase)), source=.false.)
do p=1, size(config_phase)
phase_Noutput(p) = config_phase(p)%countKeys('(output)')
phase_Nsources(p) = config_phase(p)%countKeys('(source)')
phase_Nkinematics(p) = config_phase(p)%countKeys('(kinematics)')
phase_NstiffnessDegradations(p) = config_phase(p)%countKeys('(stiffness_degradation)')

63
src/mesh/DAMASK_FEM.f90
View File

@ -66,7 +66,7 @@ program DAMASK_FEM
PetscInt :: faceSet, currentFaceSet
PetscInt :: field, dimPlex
PetscErrorCode :: ierr
integer(kind(COMPONENT_UNDEFINED_ID)) :: ID
external :: &
quit
@ -166,45 +166,28 @@ program DAMASK_FEM
!--------------------------------------------------------------------------------------------------
! boundary condition information
case('x') ! X displacement field
do field = 1, nActiveFields
if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_MECH_ID) then
do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents
if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_MECH_X_ID) then
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = &
.true.
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
IO_floatValue(line,chunkPos,i+1)
endif
enddo
endif
enddo
case('y') ! Y displacement field
do field = 1, nActiveFields
if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_MECH_ID) then
do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents
if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_MECH_Y_ID) then
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = &
.true.
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
IO_floatValue(line,chunkPos,i+1)
endif
enddo
endif
enddo
case('z') ! Z displacement field
do field = 1, nActiveFields
if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_MECH_ID) then
do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents
if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_MECH_Z_ID) then
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = &
.true.
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
IO_floatValue(line,chunkPos,i+1)
endif
enddo
endif
enddo
case('x','y','z')
select case(IO_lc(IO_stringValue(line,chunkPos,i)))
case('x')
ID = COMPONENT_MECH_X_ID
case('y')
ID = COMPONENT_MECH_Y_ID
case('z')
ID = COMPONENT_MECH_Z_ID
end select
do field = 1, nActiveFields
if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_MECH_ID) then
do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents
if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == ID) then
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = &
.true.
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
IO_floatValue(line,chunkPos,i+1)
endif
enddo
endif
enddo
end select
enddo; enddo
close(fileUnit)

139
src/mesh/FEM_mech.f90
View File

@ -67,28 +67,33 @@ contains
subroutine FEM_mech_init(fieldBC)
type(tFieldBC), intent(in) :: fieldBC
DM :: mech_mesh
PetscFE :: mechFE
PetscQuadrature :: mechQuad, functional
PetscDS :: mechDS
PetscDualSpace :: mechDualSpace
DMLabel, dimension(:),pointer :: pBCLabel
DMLabel :: BCLabel
PetscInt, dimension(:), pointer :: pNumComp, pNumDof, pBcField, pBcPoint
PetscInt :: numBC, bcSize, nc
PetscInt :: numBC, bcSize, nc, &
field, faceSet, topologDim, nNodalPoints, &
cellStart, cellEnd, cell, basis
IS :: bcPoint
IS, pointer :: pBcComps(:), pBcPoints(:)
IS, dimension(:), pointer :: pBcComps, pBcPoints
PetscSection :: section
PetscInt :: field, faceSet, topologDim, nNodalPoints
PetscReal, dimension(:), pointer :: qPointsP, qWeightsP, &
nodalPointsP, nodalWeightsP
PetscReal, allocatable, target :: nodalPoints(:), nodalWeights(:)
PetscScalar, pointer :: px_scal(:)
PetscScalar, allocatable, target :: x_scal(:)
nodalPointsP, nodalWeightsP,pV0, pCellJ, pInvcellJ
PetscReal :: detJ
PetscReal, allocatable, target :: cellJMat(:,:)
PetscReal, pointer :: pV0(:), pCellJ(:), pInvcellJ(:)
PetscInt :: cellStart, cellEnd, cell, basis
character(len=7), parameter :: prefix = 'mechFE_'
PetscScalar, pointer :: px_scal(:)
PetscScalar, allocatable, target :: x_scal(:)
character(len=*), parameter :: prefix = 'mechFE_'
PetscErrorCode :: ierr
write(6,'(/,a)') ' <<<+- FEM_mech init -+>>>'
@ -125,13 +130,19 @@ subroutine FEM_mech_init(fieldBC)
! Setup FEM mech boundary conditions
call DMGetLabel(mech_mesh,'Face Sets',BCLabel,ierr); CHKERRQ(ierr)
call DMPlexLabelComplete(mech_mesh,BCLabel,ierr); CHKERRQ(ierr)
#if (PETSC_VERSION_MINOR < 11)
call DMGetSection(mech_mesh,section,ierr); CHKERRQ(ierr)
#else
call DMGetLocalSection(mech_mesh,section,ierr); CHKERRQ(ierr)
#endif
allocate(pnumComp(1), source=dimPlex)
allocate(pnumDof(dimPlex+1), source = 0)
allocate(pnumDof(0:dimPlex), source = 0)
do topologDim = 0, dimPlex
call DMPlexGetDepthStratum(mech_mesh,topologDim,cellStart,cellEnd,ierr)
CHKERRQ(ierr)
call PetscSectionGetDof(section,cellStart,pnumDof(topologDim+1),ierr)
call PetscSectionGetDof(section,cellStart,pnumDof(topologDim),ierr)
write(6,*) 'start',cellStart,'end',cellEnd
write(6,*) 'topologDim',topologDim,'numDOF',pNumDOF(topologDim)
CHKERRQ(ierr)
enddo
numBC = 0
@ -163,9 +174,15 @@ subroutine FEM_mech_init(fieldBC)
endif
endif
enddo; enddo
#if (PETSC_VERSION_MINOR < 11)
call DMPlexCreateSection(mech_mesh,dimPlex,1,pNumComp,pNumDof, &
numBC,pBcField,pBcComps,pBcPoints,PETSC_NULL_IS, &
section,ierr)
numBC,pBcField,pBcComps,pBcPoints,PETSC_NULL_IS,section,ierr)
#else
allocate(pBClabel(1),source=BClabel)
call DMPlexCreateSection(mech_mesh,pBClabel,pNumComp,pNumDof, &
numBC,pBcField,pBcComps,pBcPoints,PETSC_NULL_IS,section,ierr)
#endif
CHKERRQ(ierr)
call DMSetSection(mech_mesh,section,ierr); CHKERRQ(ierr)
do faceSet = 1, numBC
@ -196,13 +213,11 @@ subroutine FEM_mech_init(fieldBC)
call VecSet(solution ,0.0,ierr); CHKERRQ(ierr)
call VecSet(solution_rate ,0.0,ierr); CHKERRQ(ierr)
allocate(x_scal(cellDof))
allocate(nodalPoints (dimPlex))
allocate(nodalWeights(1))
nodalPointsP => nodalPoints
nodalWeightsP => nodalWeights
allocate(nodalWeightsP(1))
allocate(nodalPointsP(dimPlex))
allocate(pv0(dimPlex))
allocate(pcellJ(dimPlex*dimPlex))
allocate(pinvcellJ(dimPlex*dimPlex))
allocate(pcellJ(dimPlex**2))
allocate(pinvcellJ(dimPlex**2))
allocate(cellJMat(dimPlex,dimPlex))
call DMGetSection(mech_mesh,section,ierr); CHKERRQ(ierr)
call DMGetDS(mech_mesh,mechDS,ierr); CHKERRQ(ierr)
@ -212,7 +227,7 @@ subroutine FEM_mech_init(fieldBC)
call DMPlexGetHeightStratum(mech_mesh,0,cellStart,cellEnd,ierr)
CHKERRQ(ierr)
do cell = cellStart, cellEnd-1 !< loop over all elements
x_scal = 0.0
x_scal = 0.0_pReal
call DMPlexComputeCellGeometryAffineFEM(mech_mesh,cell,pV0,pCellJ,pInvcellJ,detJ,ierr)
CHKERRQ(ierr)
cellJMat = reshape(pCellJ,shape=[dimPlex,dimPlex])
@ -221,7 +236,7 @@ subroutine FEM_mech_init(fieldBC)
CHKERRQ(ierr)
call PetscQuadratureGetData(functional,dimPlex,nc,nNodalPoints,nodalPointsP,nodalWeightsP,ierr)
CHKERRQ(ierr)
x_scal(basis+1:basis+dimPlex) = pV0 + matmul(transpose(cellJMat),nodalPointsP + 1.0)
x_scal(basis+1:basis+dimPlex) = pV0 + matmul(transpose(cellJMat),nodalPointsP + 1.0_pReal)
enddo
px_scal => x_scal
call DMPlexVecSetClosure(mech_mesh,section,solution_local,cell,px_scal,INSERT_ALL_VALUES,ierr)
@ -283,6 +298,9 @@ end function FEM_mech_solution
subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr)
DM :: dm_local
PetscObject,intent(in) :: dummy
PetscErrorCode :: ierr
PetscDS :: prob
Vec :: x_local, f_local, xx_local
PetscSection :: section
@ -294,10 +312,10 @@ subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr)
qPt, basis, comp, cidx
PetscReal :: detFAvg
PetscReal :: BMat(dimPlex*dimPlex,cellDof)
PetscObject,intent(in) :: dummy
PetscInt :: bcSize
IS :: bcPoints
PetscErrorCode :: ierr
allocate(pV0(dimPlex))
allocate(pcellJ(dimPlex**2))
@ -316,7 +334,7 @@ subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr)
call DMGetStratumIS(dm_local,'Face Sets',mesh_boundaries(face),bcPoints,ierr)
CHKERRQ(ierr)
call utilities_projectBCValues(x_local,section,0,field-1,bcPoints, &
0.0,params%fieldBC%componentBC(field)%Value(face),params%timeinc)
0.0_pReal,params%fieldBC%componentBC(field)%Value(face),params%timeinc)
call ISDestroy(bcPoints,ierr); CHKERRQ(ierr)
endif
endif
@ -403,29 +421,35 @@ end subroutine FEM_mech_formResidual
!--------------------------------------------------------------------------------------------------
subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr)
DM :: dm_local
DM :: dm_local
Mat :: Jac_pre, Jac
PetscObject, intent(in) :: dummy
PetscErrorCode :: ierr
PetscDS :: prob
Vec :: x_local, xx_local
Mat :: Jac_pre, Jac
PetscSection :: section, gSection
PetscReal, dimension(1, cellDof) :: MatB
PetscReal, dimension(dimPlex**2,cellDof) :: BMat, BMatAvg, MatA
PetscReal, dimension(3,3) :: F, FAvg, FInv
PetscReal :: detJ
PetscReal, dimension(:), pointer :: basisField, basisFieldDer, &
pV0, pCellJ, pInvcellJ
PetscScalar, dimension(:), pointer :: pK_e, x_scal
PetscScalar,dimension(cellDOF,cellDOF), target :: K_e
PetscScalar,dimension(cellDOF,cellDOF) :: K_eA , &
K_eB
PetscInt :: cellStart, cellEnd, cell, field, face, &
qPt, basis, comp, cidx,bcSize
PetscScalar,dimension(cellDOF,cellDOF), target :: K_e, &
K_eA , &
K_eB
PetscScalar, target :: K_eVec(cellDof*cellDof)
PetscReal :: BMat (dimPlex*dimPlex,cellDof), &
BMatAvg(dimPlex*dimPlex,cellDof), &
MatA (dimPlex*dimPlex,cellDof), &
MatB (1 ,cellDof)
PetscScalar, dimension(:), pointer :: pK_e, x_scal
PetscReal, dimension(3,3) :: F, FAvg, FInv
PetscObject, intent(in) :: dummy
IS :: bcPoints
PetscErrorCode :: ierr
allocate(pV0(dimPlex))
allocate(pcellJ(dimPlex**2))
@ -440,7 +464,7 @@ subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr)
call DMGetGlobalSection(dm_local,gSection,ierr); CHKERRQ(ierr)
call DMGetLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr)
call VecWAXPY(x_local,1.0,xx_local,solution_local,ierr); CHKERRQ(ierr)
call VecWAXPY(x_local,1.0_pReal,xx_local,solution_local,ierr); CHKERRQ(ierr)
do field = 1, dimPlex; do face = 1, mesh_Nboundaries
if (params%fieldBC%componentBC(field)%Mask(face)) then
call DMGetStratumSize(dm_local,'Face Sets',mesh_boundaries(face),bcSize,ierr)
@ -448,7 +472,7 @@ subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr)
call DMGetStratumIS(dm_local,'Face Sets',mesh_boundaries(face),bcPoints,ierr)
CHKERRQ(ierr)
call utilities_projectBCValues(x_local,section,0,field-1,bcPoints, &
0.0,params%fieldBC%componentBC(field)%Value(face),params%timeinc)
0.0_pReal,params%fieldBC%componentBC(field)%Value(face),params%timeinc)
call ISDestroy(bcPoints,ierr); CHKERRQ(ierr)
endif
endif
@ -501,13 +525,16 @@ subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr)
(matmul(matmul(transpose(BMatAvg), &
reshape(FInv(1:dimPlex,1:dimPlex),shape=[dimPlex*dimPlex,1],order=[2,1])),MatB) + &
K_eB)/real(dimPlex)
else
K_e = K_eA
endif
K_e = K_e + eps*math_identity2nd(cellDof)
K_eVec = reshape(K_e, [cellDof*cellDof])*abs(detJ)
pK_e => K_eVec
K_e = (K_e + eps*math_identity2nd(cellDof)) * abs(detJ)
#ifndef __INTEL_COMPILER
pK_e(1:cellDOF**2) => K_e
#else
! https://software.intel.com/en-us/forums/intel-fortran-compiler/topic/782230 (bug)
allocate(pK_e(cellDOF**2),source = reshape(K_e,[cellDOF**2]))
#endif
call DMPlexMatSetClosure(dm_local,section,gSection,Jac,cell,pK_e,ADD_VALUES,ierr)
CHKERRQ(ierr)
call DMPlexVecRestoreClosure(dm_local,section,x_local,cell,x_scal,ierr)
@ -536,18 +563,18 @@ subroutine FEM_mech_forward(guess,timeinc,timeinc_old,fieldBC)
type(tFieldBC), intent(in) :: &
fieldBC
real(pReal), intent(in) :: &
real(pReal), intent(in) :: &
timeinc_old, &
timeinc
logical, intent(in) :: &
logical, intent(in) :: &
guess
PetscInt :: field, face
DM :: dm_local
Vec :: x_local
PetscSection :: section
PetscInt :: bcSize
IS :: bcPoints
PetscErrorCode :: ierr
PetscInt :: field, face, bcSize
DM :: dm_local
Vec :: x_local
PetscSection :: section
IS :: bcPoints
PetscErrorCode :: ierr
!--------------------------------------------------------------------------------------------------
! forward last inc
@ -557,7 +584,7 @@ subroutine FEM_mech_forward(guess,timeinc,timeinc_old,fieldBC)
call SNESGetDM(mech_snes,dm_local,ierr); CHKERRQ(ierr) !< retrieve mesh info from mech_snes into dm_local
call DMGetSection(dm_local,section,ierr); CHKERRQ(ierr)
call DMGetLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr)
call VecSet(x_local,0.0,ierr); CHKERRQ(ierr)
call VecSet(x_local,0.0_pReal,ierr); CHKERRQ(ierr)
call DMGlobalToLocalBegin(dm_local,solution,INSERT_VALUES,x_local,ierr) !< retrieve my partition of global solution vector
CHKERRQ(ierr)
call DMGlobalToLocalEnd(dm_local,solution,INSERT_VALUES,x_local,ierr)
@ -570,7 +597,7 @@ subroutine FEM_mech_forward(guess,timeinc,timeinc_old,fieldBC)
call DMGetStratumIS(dm_local,'Face Sets',mesh_boundaries(face),bcPoints,ierr)
CHKERRQ(ierr)
call utilities_projectBCValues(solution_local,section,0,field-1,bcPoints, &
0.0,fieldBC%componentBC(field)%Value(face),timeinc_old)
0.0_pReal,fieldBC%componentBC(field)%Value(face),timeinc_old)
call ISDestroy(bcPoints,ierr); CHKERRQ(ierr)
endif
endif

17
src/mesh/FEM_utilities.f90
View File

@ -23,7 +23,7 @@ module FEM_utilities
private
!--------------------------------------------------------------------------------------------------
logical, public :: cutBack = .false. !< cut back of BVP solver in case convergence is not achieved or a material point is terminally ill
logical, public :: cutBack = .false. !< cut back of BVP solver in case convergence is not achieved or a material point is terminally ill
integer, public, parameter :: maxFields = 6
integer, public :: nActiveFields = 0
@ -56,15 +56,15 @@ module FEM_utilities
!--------------------------------------------------------------------------------------------------
! derived types
type, public :: tSolutionState !< return type of solution from FEM solver variants
logical :: converged = .true.
logical :: stagConverged = .true.
integer :: iterationsNeeded = 0
logical :: converged = .true.
logical :: stagConverged = .true.
integer :: iterationsNeeded = 0
end type tSolutionState
type, public :: tComponentBC
integer(kind(COMPONENT_UNDEFINED_ID)) :: ID
real(pReal), allocatable :: Value(:)
logical, allocatable :: Mask(:)
real(pReal), allocatable, dimension(:) :: Value
logical, allocatable, dimension(:) :: Mask
end type tComponentBC
type, public :: tFieldBC
@ -79,8 +79,8 @@ module FEM_utilities
outputfrequency = 1, & !< frequency of result writes
logscale = 0 !< linear/logarithmic time inc flag
logical :: followFormerTrajectory = .true. !< follow trajectory of former loadcase
integer, allocatable :: faceID(:)
type(tFieldBC), allocatable :: fieldBC(:)
integer, allocatable, dimension(:) :: faceID
type(tFieldBC), allocatable, dimension(:) :: fieldBC
end type tLoadCase
public :: &
@ -88,6 +88,7 @@ module FEM_utilities
utilities_constitutiveResponse, &
utilities_projectBCValues, &
FIELD_MECH_ID, &
COMPONENT_UNDEFINED_ID, &
COMPONENT_MECH_X_ID, &
COMPONENT_MECH_Y_ID, &
COMPONENT_MECH_Z_ID

4
src/mesh/FEM_zoo.f90
View File

@ -37,7 +37,7 @@ contains
!--------------------------------------------------------------------------------------------------
subroutine FEM_Zoo_init
write(6,'(/,a)') ' <<<+- FEM_Zoo init -+>>>'
write(6,'(/,a)') ' <<<+- FEM_Zoo init -+>>>'; flush(6)
!--------------------------------------------------------------------------------------------------
! 2D linear
@ -53,7 +53,7 @@ subroutine FEM_Zoo_init
! 2D quadratic
FEM_Zoo_nQuadrature(2,2) = 3
allocate(FEM_Zoo_QuadratureWeights(2,2)%p(3))
allocate(FEM_Zoo_QuadratureWeights(2,2)%p(3))
FEM_Zoo_QuadratureWeights(2,2)%p(1:3) = 1.0_pReal/3.0_pReal
allocate(FEM_Zoo_QuadraturePoints (2,2)%p(6))

241
src/mesh_FEM.f90
View File

@ -26,9 +26,11 @@ module mesh
integer, public, protected :: &
mesh_Nboundaries, &
mesh_NcpElems, & !< total number of CP elements in mesh
mesh_NcpElemsGlobal
integer :: &
mesh_NcpElems !< total number of CP elements in mesh
!!!! BEGIN DEPRECATED !!!!!
integer, public, protected :: &
mesh_maxNips !< max number of IPs in any CP element
@ -44,7 +46,7 @@ module mesh
mesh_ipVolume, & !< volume associated with IP (initially!)
mesh_node0 !< node x,y,z coordinates (initially!)
real(pReal), dimension(:,:,:), allocatable, public :: &
real(pReal), dimension(:,:,:), allocatable :: &
mesh_ipCoordinates !< IP x,y,z coordinates (after deformation!)
DM, public :: geomMesh
@ -66,128 +68,126 @@ contains
!--------------------------------------------------------------------------------------------------
subroutine mesh_init
integer, dimension(1), parameter:: FE_geomtype = [1] !< geometry type of particular element type
integer, dimension(1), parameter:: FE_geomtype = [1] !< geometry type of particular element type
integer, dimension(1) :: FE_Nips !< number of IPs in a specific type of element
integer, dimension(1) :: FE_Nips !< number of IPs in a specific type of element
integer, parameter :: FILEUNIT = 222
integer :: j
integer, allocatable, dimension(:) :: chunkPos
integer :: dimPlex, &
mesh_Nnodes !< total number of nodes in mesh
integer, parameter :: &
mesh_ElemType=1 !< Element type of the mesh (only support homogeneous meshes)
character(len=512) :: &
line
logical :: flag
PetscSF :: sf
DM :: globalMesh
PetscInt :: nFaceSets
PetscInt, pointer :: pFaceSets(:)
IS :: faceSetIS
PetscErrorCode :: ierr
integer, parameter :: FILEUNIT = 222
integer :: j
integer, allocatable, dimension(:) :: chunkPos
integer :: dimPlex, &
mesh_Nnodes !< total number of nodes in mesh
integer, parameter :: &
mesh_ElemType=1 !< Element type of the mesh (only support homogeneous meshes)
character(len=512) :: &
line
logical :: flag
PetscSF :: sf
DM :: globalMesh
PetscInt :: face, nFaceSets
PetscInt, pointer :: pFaceSets(:)
IS :: faceSetIS
PetscErrorCode :: ierr
write(6,'(/,a)') ' <<<+- mesh init -+>>>'
! read in file
call DMPlexCreateFromFile(PETSC_COMM_WORLD,geometryFile,PETSC_TRUE,globalMesh,ierr)
CHKERRQ(ierr)
! get spatial dimension (2 or 3?)
call DMGetDimension(globalMesh,dimPlex,ierr)
CHKERRQ(ierr)
write(6,*) 'dimension',dimPlex;flush(6)
call DMGetStratumSize(globalMesh,'depth',dimPlex,mesh_NcpElemsGlobal,ierr)
CHKERRQ(ierr)
! get number of IDs in face sets (for boundary conditions?)
call DMGetLabelSize(globalMesh,'Face Sets',mesh_Nboundaries,ierr)
CHKERRQ(ierr)
write(6,*) 'number of "Face Sets"',mesh_Nboundaries;flush(6)
call MPI_Bcast(mesh_Nboundaries,1,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr)
call MPI_Bcast(mesh_NcpElemsGlobal,1,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr)
call MPI_Bcast(dimPlex,1,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr)
write(6,'(/,a)') ' <<<+- mesh init -+>>>'
allocate(mesh_boundaries(mesh_Nboundaries), source = 0)
call DMGetLabelSize(globalMesh,'Face Sets',nFaceSets,ierr)
CHKERRQ(ierr)
call DMGetLabelIdIS(globalMesh,'Face Sets',faceSetIS,ierr)
CHKERRQ(ierr)
if (nFaceSets > 0) then
call ISGetIndicesF90(faceSetIS,pFaceSets,ierr)
CHKERRQ(ierr)
mesh_boundaries(1:nFaceSets) = pFaceSets
CHKERRQ(ierr)
call ISRestoreIndicesF90(faceSetIS,pFaceSets,ierr)
endif
call MPI_Bcast(mesh_boundaries,mesh_Nboundaries,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr)
! read in file
call DMPlexCreateFromFile(PETSC_COMM_WORLD,geometryFile,PETSC_TRUE,globalMesh,ierr)
CHKERRQ(ierr)
! get spatial dimension (2 or 3?)
call DMGetDimension(globalMesh,dimPlex,ierr)
CHKERRQ(ierr)
write(6,*) 'dimension',dimPlex;flush(6)
call DMGetStratumSize(globalMesh,'depth',dimPlex,mesh_NcpElemsGlobal,ierr)
CHKERRQ(ierr)
! get number of IDs in face sets (for boundary conditions?)
call DMGetLabelSize(globalMesh,'Face Sets',mesh_Nboundaries,ierr)
CHKERRQ(ierr)
write(6,*) 'number of "Face Sets"',mesh_Nboundaries;flush(6)
call MPI_Bcast(mesh_Nboundaries,1,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr)
call MPI_Bcast(mesh_NcpElemsGlobal,1,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr)
call MPI_Bcast(dimPlex,1,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr)
! this read in function should ignore C and C++ style comments
! it is used for BC only?
if (worldrank == 0) then
j = 0
flag = .false.
call IO_open_file(FILEUNIT,trim(geometryFile))
do
read(FILEUNIT,'(A)') line
if (trim(line) == IO_EOF) exit ! skip empty lines
if (trim(line) == '$Elements') then
read(FILEUNIT,'(A)') line ! number of elements (ignore)
read(FILEUNIT,'(A)') line
flag = .true.
endif
if (trim(line) == '$EndElements') exit
if (flag) then
chunkPos = IO_stringPos(line)
if (chunkPos(1) == 3+IO_intValue(line,chunkPos,3)+dimPlex+1) then
call DMSetLabelValue(globalMesh,'material',j,IO_intValue(line,chunkPos,4),ierr)
CHKERRQ(ierr)
j = j + 1
endif ! count all identifiers to allocate memory and do sanity check
endif
enddo
close (FILEUNIT)
call DMClone(globalMesh,geomMesh,ierr)
CHKERRQ(ierr)
else
call DMPlexDistribute(globalMesh,0,sf,geomMesh,ierr)
CHKERRQ(ierr)
endif
allocate(mesh_boundaries(mesh_Nboundaries), source = 0)
call DMGetLabelSize(globalMesh,'Face Sets',nFaceSets,ierr)
CHKERRQ(ierr)
call DMGetLabelIdIS(globalMesh,'Face Sets',faceSetIS,ierr)
CHKERRQ(ierr)
if (nFaceSets > 0) call ISGetIndicesF90(faceSetIS,pFaceSets,ierr)
do face = 1, nFaceSets
mesh_boundaries(face) = pFaceSets(face)
enddo
if (nFaceSets > 0) call ISRestoreIndicesF90(faceSetIS,pFaceSets,ierr)
call MPI_Bcast(mesh_boundaries,mesh_Nboundaries,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr)
call DMDestroy(globalMesh,ierr); CHKERRQ(ierr)
! this read in function should ignore C and C++ style comments
! it is used for BC only?
if (worldrank == 0) then
j = 0
flag = .false.
call IO_open_file(FILEUNIT,trim(geometryFile))
do
read(FILEUNIT,'(a512)') line
if (trim(line) == IO_EOF) exit ! skip empty lines
if (trim(line) == '$Elements') then
read(FILEUNIT,'(a512)') line ! number of elements (ignore)
read(FILEUNIT,'(a512)') line
flag = .true.
endif
if (trim(line) == '$EndElements') exit
if (flag) then
chunkPos = IO_stringPos(line)
if (chunkPos(1) == 3+IO_intValue(line,chunkPos,3)+dimPlex+1) then
call DMSetLabelValue(globalMesh,'material',j,IO_intValue(line,chunkPos,4),ierr)
CHKERRQ(ierr)
j = j + 1
endif ! count all identifiers to allocate memory and do sanity check
endif
enddo
close (FILEUNIT)
call DMClone(globalMesh,geomMesh,ierr)
CHKERRQ(ierr)
else
call DMPlexDistribute(globalMesh,0,sf,geomMesh,ierr)
CHKERRQ(ierr)
endif
call DMGetStratumSize(geomMesh,'depth',dimPlex,mesh_NcpElems,ierr)
CHKERRQ(ierr)
call DMGetStratumSize(geomMesh,'depth',0,mesh_Nnodes,ierr)
CHKERRQ(ierr)
call DMDestroy(globalMesh,ierr); CHKERRQ(ierr)
FE_Nips(FE_geomtype(1)) = FEM_Zoo_nQuadrature(dimPlex,integrationOrder)
mesh_maxNips = FE_Nips(1)
call DMGetStratumSize(geomMesh,'depth',dimPlex,mesh_NcpElems,ierr)
CHKERRQ(ierr)
call DMGetStratumSize(geomMesh,'depth',0,mesh_Nnodes,ierr)
CHKERRQ(ierr)
write(6,*) 'mesh_maxNips',mesh_maxNips
call mesh_FEM_build_ipCoordinates(dimPlex,FEM_Zoo_QuadraturePoints(dimPlex,integrationOrder)%p)
call mesh_FEM_build_ipVolumes(dimPlex)
FE_Nips(FE_geomtype(1)) = FEM_Zoo_nQuadrature(dimPlex,integrationOrder)
mesh_maxNips = FE_Nips(1)
allocate (mesh_element (4,mesh_NcpElems)); mesh_element = 0
do j = 1, mesh_NcpElems
mesh_element( 1,j) = -1 ! DEPRECATED
mesh_element( 2,j) = mesh_elemType ! elem type
mesh_element( 3,j) = 1 ! homogenization
call DMGetLabelValue(geomMesh,'material',j-1,mesh_element(4,j),ierr)
CHKERRQ(ierr)
end do
write(6,*) 'mesh_maxNips',mesh_maxNips
call mesh_FEM_build_ipCoordinates(dimPlex,FEM_Zoo_QuadraturePoints(dimPlex,integrationOrder)%p)
call mesh_FEM_build_ipVolumes(dimPlex)
if (debug_e < 1 .or. debug_e > mesh_NcpElems) &
call IO_error(602,ext_msg='element') ! selected element does not exist
if (debug_i < 1 .or. debug_i > FE_Nips(FE_geomtype(mesh_element(2,debug_e)))) &
call IO_error(602,ext_msg='IP') ! selected element does not have requested IP
allocate (mesh_element (4,mesh_NcpElems)); mesh_element = 0
do j = 1, mesh_NcpElems
mesh_element( 1,j) = -1 ! DEPRECATED
mesh_element( 2,j) = mesh_elemType ! elem type
mesh_element( 3,j) = 1 ! homogenization
call DMGetLabelValue(geomMesh,'material',j-1,mesh_element(4,j),ierr)
CHKERRQ(ierr)
end do
FEsolving_execElem = [ 1,mesh_NcpElems ] ! parallel loop bounds set to comprise all DAMASK elements
FEsolving_execIP = [1,FE_Nips(FE_geomtype(mesh_element(2,1)))]
if (debug_e < 1 .or. debug_e > mesh_NcpElems) &
call IO_error(602,ext_msg='element') ! selected element does not exist
if (debug_i < 1 .or. debug_i > FE_Nips(FE_geomtype(mesh_element(2,debug_e)))) &
call IO_error(602,ext_msg='IP') ! selected element does not have requested IP
allocate(mesh_node0(3,mesh_Nnodes),source=0.0_pReal)
FEsolving_execElem = [ 1,mesh_NcpElems ] ! parallel loop bounds set to comprise all DAMASK elements
if (allocated(FEsolving_execIP)) deallocate(FEsolving_execIP)
allocate(FEsolving_execIP(2,mesh_NcpElems)); FEsolving_execIP = 1 ! parallel loop bounds set to comprise from first IP...
forall (j = 1:mesh_NcpElems) FEsolving_execIP(2,j) = FE_Nips(FE_geomtype(mesh_element(2,j))) ! ...up to own IP count for each element
allocate(mesh_node0(3,mesh_Nnodes),source=0.0_pReal)
call discretization_init(mesh_element(3,:),mesh_element(4,:),&
call discretization_init(mesh_element(3,:),mesh_element(4,:),&
reshape(mesh_ipCoordinates,[3,mesh_maxNips*mesh_NcpElems]), &
mesh_node0)
@ -199,21 +199,17 @@ end subroutine mesh_init
!--------------------------------------------------------------------------------------------------
subroutine mesh_FEM_build_ipVolumes(dimPlex)
PetscInt :: dimPlex
PetscInt,intent(in):: dimPlex
PetscReal :: vol
PetscReal, target :: cent(dimPlex), norm(dimPlex)
PetscReal, pointer :: pCent(:), pNorm(:)
PetscReal, pointer,dimension(:) :: pCent, pNorm
PetscInt :: cellStart, cellEnd, cell
PetscErrorCode :: ierr
if (.not. allocated(mesh_ipVolume)) then
allocate(mesh_ipVolume(mesh_maxNips,mesh_NcpElems))
mesh_ipVolume = 0.0_pReal
endif
allocate(mesh_ipVolume(mesh_maxNips,mesh_NcpElems),source=0.0_pReal)
call DMPlexGetHeightStratum(geomMesh,0,cellStart,cellEnd,ierr); CHKERRQ(ierr)
pCent => cent
pNorm => norm
allocate(pCent(dimPlex))
allocate(pNorm(dimPlex))
do cell = cellStart, cellEnd-1
call DMPlexComputeCellGeometryFVM(geomMesh,cell,vol,pCent,pNorm,ierr)
CHKERRQ(ierr)
@ -231,8 +227,7 @@ subroutine mesh_FEM_build_ipCoordinates(dimPlex,qPoints)
PetscInt, intent(in) :: dimPlex
PetscReal, intent(in) :: qPoints(mesh_maxNips*dimPlex)
PetscReal, target :: v0(dimPlex), cellJ(dimPlex*dimPlex), invcellJ(dimPlex*dimPlex)
PetscReal, pointer :: pV0(:), pCellJ(:), pInvcellJ(:)
PetscReal, pointer,dimension(:) :: pV0, pCellJ, pInvcellJ
PetscReal :: detJ
PetscInt :: cellStart, cellEnd, cell, qPt, dirI, dirJ, qOffset
PetscErrorCode :: ierr
@ -240,9 +235,9 @@ subroutine mesh_FEM_build_ipCoordinates(dimPlex,qPoints)
allocate(mesh_ipCoordinates(3,mesh_maxNips,mesh_NcpElems),source=0.0_pReal)
pV0 => v0
pCellJ => cellJ
pInvcellJ => invcellJ
allocate(pV0(dimPlex))
allocatE(pCellJ(dimPlex**2))
allocatE(pinvCellJ(dimPlex**2))
call DMPlexGetHeightStratum(geomMesh,0,cellStart,cellEnd,ierr); CHKERRQ(ierr)
do cell = cellStart, cellEnd-1 !< loop over all elements
call DMPlexComputeCellGeometryAffineFEM(geomMesh,cell,pV0,pCellJ,pInvcellJ,detJ,ierr)

2
src/mesh_grid.f90
View File

@ -98,7 +98,7 @@ subroutine mesh_init(ip,el)
worldrank<1))
FEsolving_execElem = [1,product(myGrid)] ! parallel loop bounds set to comprise all elements
allocate(FEsolving_execIP(2,product(myGrid)),source=1) ! parallel loop bounds set to comprise the only IP
FEsolving_execIP = [1,1] ! parallel loop bounds set to comprise the only IP
!--------------------------------------------------------------------------------------------------
! store geometry information for post processing

195
src/mesh_marc.f90
View File

@ -40,7 +40,6 @@ module mesh
mesh_init, &
mesh_FEasCP
contains
!--------------------------------------------------------------------------------------------------
@ -82,8 +81,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]
FEsolving_execIP = spread([1,elem%nIPs],2,nElems)
FEsolving_execElem = [1,nElems]
FEsolving_execIP = [1,elem%nIPs]
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"
@ -185,8 +184,6 @@ subroutine inputRead(elem,node0_elem,connectivity_elem,microstructureAt,homogeni
initialcondTableStyle, &
nNodes, &
nElems
integer, parameter :: &
FILEUNIT = 222
integer, dimension(:), allocatable :: &
matNumber !< material numbers for hypoelastic material
character(len=pStringLen), dimension(:), allocatable :: inputFile !< file content, separated per lines
@ -207,10 +204,9 @@ subroutine inputRead(elem,node0_elem,connectivity_elem,microstructureAt,homogeni
call inputRead_NnodesAndElements(nNodes,nElems,&
inputFile)
call IO_open_inputFile(FILEUNIT) ! ToDo: It would be better to use fileContent
call inputRead_mapElemSets(nameElemSet,mapElemSet,&
FILEUNIT,inputFile)
inputFile)
call inputRead_elemType(elem, &
nElems,inputFile)
@ -229,9 +225,7 @@ subroutine inputRead(elem,node0_elem,connectivity_elem,microstructureAt,homogeni
call inputRead_microstructureAndHomogenization(microstructureAt,homogenizationAt, &
nElems,elem%nNodes,nameElemSet,mapElemSet,&
initialcondTableStyle,FILEUNIT)
close(FILEUNIT)
initialcondTableStyle,inputFile)
end subroutine inputRead
@ -241,8 +235,8 @@ end subroutine inputRead
!--------------------------------------------------------------------------------------------------
subroutine inputRead_fileFormat(fileFormat,fileContent)
integer, intent(out) :: fileFormat
character(len=pStringLen), dimension(:), intent(in) :: fileContent !< file content, separated per lines
integer, intent(out) :: fileFormat
character(len=*), dimension(:), intent(in) :: fileContent !< file content, separated per lines
integer, allocatable, dimension(:) :: chunkPos
integer :: l
@ -263,8 +257,8 @@ end subroutine inputRead_fileFormat
!--------------------------------------------------------------------------------------------------
subroutine inputRead_tableStyles(initialcond,hypoelastic,fileContent)
integer, intent(out) :: initialcond, hypoelastic
character(len=pStringLen), dimension(:), intent(in) :: fileContent !< file content, separated per lines
integer, intent(out) :: initialcond, hypoelastic
character(len=*), dimension(:), intent(in) :: fileContent !< file content, separated per lines
integer, allocatable, dimension(:) :: chunkPos
integer :: l
@ -290,10 +284,9 @@ end subroutine inputRead_tableStyles
subroutine inputRead_matNumber(matNumber, &
tableStyle,fileContent)
integer, allocatable, dimension(:), intent(out) :: matNumber
integer, intent(in) :: tableStyle
character(len=pStringLen), dimension(:), intent(in) :: fileContent !< file content, separated per lines
integer, allocatable, dimension(:), intent(out) :: matNumber
integer, intent(in) :: tableStyle
character(len=*), dimension(:), intent(in) :: fileContent !< file content, separated per lines
integer, allocatable, dimension(:) :: chunkPos
integer :: i, j, data_blocks, l
@ -326,8 +319,8 @@ end subroutine inputRead_matNumber
subroutine inputRead_NnodesAndElements(nNodes,nElems,&
fileContent)
integer, intent(out) :: nNodes, nElems
character(len=pStringLen), dimension(:), intent(in) :: fileContent !< file content, separated per lines
integer, intent(out) :: nNodes, nElems
character(len=*), dimension(:), intent(in) :: fileContent !< file content, separated per lines
integer, allocatable, dimension(:) :: chunkPos
integer :: l
@ -354,11 +347,11 @@ end subroutine inputRead_NnodesAndElements
subroutine inputRead_NelemSets(nElemSets,maxNelemInSet,&
fileContent)
integer, intent(out) :: nElemSets, maxNelemInSet
character(len=pStringLen), dimension(:), intent(in) :: fileContent !< file content, separated per lines
integer, intent(out) :: nElemSets, maxNelemInSet
character(len=*), dimension(:), intent(in) :: fileContent !< file content, separated per lines
integer, allocatable, dimension(:) :: chunkPos
integer :: i,l,elemInCurrentSet
integer :: i,l,elemInCurrentSet
nElemSets = 0
maxNelemInSet = 0
@ -396,16 +389,15 @@ end subroutine inputRead_NelemSets
!--------------------------------------------------------------------------------------------------
!> @brief map element sets
!--------------------------------------------------------------------------------------------------
subroutine inputRead_mapElemSets(nameElemSet,mapElemSet,fileUnit,fileContent)
subroutine inputRead_mapElemSets(nameElemSet,mapElemSet,&
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
character(len=*), dimension(:), intent(in) :: fileContent !< file content, separated per lines
integer, allocatable, dimension(:) :: chunkPos
character(len=300) :: line
integer :: elemSet, NelemSets, maxNelemInSet
integer :: elemSet, NelemSets, maxNelemInSet,l
call inputRead_NelemSets(NelemSets,maxNelemInSet,fileContent)
@ -413,20 +405,17 @@ subroutine inputRead_mapElemSets(nameElemSet,mapElemSet,fileUnit,fileContent)
allocate(mapElemSet(1+maxNelemInSet,NelemSets),source=0)
elemSet = 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
elemSet = elemSet+1
nameElemSet(elemSet) = trim(IO_stringValue(line,chunkPos,4))
mapElemSet(:,elemSet) = IO_continuousIntValues(fileUnit,size(mapElemSet,1)-1,nameElemSet,mapElemSet,size(nameElemSet))
nameElemSet(elemSet) = trim(IO_stringValue(fileContent(l),chunkPos,4))
mapElemSet(:,elemSet) = continuousIntValues(fileContent(l+1:),size(mapElemSet,1)-1,nameElemSet,mapElemSet,size(nameElemSet))
endif
enddo
620 end subroutine inputRead_mapElemSets
end subroutine inputRead_mapElemSets
!--------------------------------------------------------------------------------------------------
@ -437,7 +426,7 @@ 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
character(len=*), dimension(:), intent(in) :: fileContent !< file content, separated per lines
integer, allocatable, dimension(:) :: chunkPos
integer :: i,j,l,nNodesAlreadyRead
@ -470,7 +459,7 @@ end subroutine inputRead_mapElems
!--------------------------------------------------------------------------------------------------
subroutine inputRead_mapNodes(fileContent)
character(len=pStringLen), dimension(:), intent(in) :: fileContent !< file content, separated per lines
character(len=*), dimension(:), intent(in) :: fileContent !< file content, separated per lines
integer, allocatable, dimension(:) :: chunkPos
integer :: i, l
@ -499,7 +488,7 @@ subroutine inputRead_elemNodes(nodes, &
real(pReal), allocatable, dimension(:,:), intent(out) :: nodes
integer, intent(in) :: nNode
character(len=pStringLen), dimension(:), intent(in) :: fileContent !< file content, separated per lines
character(len=*), dimension(:), intent(in) :: fileContent !< file content, separated per lines
integer, allocatable, dimension(:) :: chunkPos
integer :: i,j,m,l
@ -529,9 +518,9 @@ end subroutine inputRead_elemNodes
subroutine inputRead_elemType(elem, &
nElem,fileContent)
type(tElement), intent(out) :: elem
integer, intent(in) :: nElem
character(len=pStringLen), dimension(:), intent(in) :: fileContent !< file content, separated per lines
type(tElement), intent(out) :: elem
integer, intent(in) :: nElem
character(len=*), dimension(:), intent(in) :: fileContent !< file content, separated per lines
integer, allocatable, dimension(:) :: chunkPos
integer :: i,j,t,l,remainingChunks
@ -617,7 +606,7 @@ function inputRead_connectivityElem(nElem,nNodes,fileContent)
integer, intent(in) :: &
nElem, &
nNodes !< number of nodes per element
character(len=pStringLen), dimension(:), intent(in) :: fileContent !< file content, separated per lines
character(len=*), dimension(:), intent(in) :: fileContent !< file content, separated per lines
integer, dimension(nNodes,nElem) :: &
inputRead_connectivityElem
@ -662,7 +651,7 @@ end function inputRead_connectivityElem
!> @brief Stores homogenization and microstructure ID
!--------------------------------------------------------------------------------------------------
subroutine inputRead_microstructureAndHomogenization(microstructureAt,homogenizationAt, &
nElem,nNodes,nameElemSet,mapElemSet,initialcondTableStyle,fileUnit)
nElem,nNodes,nameElemSet,mapElemSet,initialcondTableStyle,fileContent)
integer, dimension(:), allocatable, intent(out) :: &
microstructureAt, &
@ -670,60 +659,46 @@ subroutine inputRead_microstructureAndHomogenization(microstructureAt,homogeniza
integer, intent(in) :: &
nElem, &
nNodes, & !< number of nodes per element
initialcondTableStyle, &
fileUnit
character(len=64), dimension(:), intent(in) :: &
nameElemSet
integer, dimension(:,:), intent(in) :: &
mapElemSet !< list of elements in elementSet
initialcondTableStyle
character(len=*), dimension(:), intent(in) :: nameElemSet
integer, dimension(:,:), intent(in) :: mapElemSet !< list of elements in elementSet
character(len=*), dimension(:), intent(in) :: fileContent !< file content, separated per lines
integer, allocatable, dimension(:) :: chunkPos
character(len=300) line
integer, dimension(1+nElem) :: contInts
integer :: i,j,t,sv,myVal,e,nNodesAlreadyRead
integer :: i,j,t,sv,myVal,e,nNodesAlreadyRead,l,k,m
allocate(microstructureAt(nElem),source=0)
allocate(homogenizationAt(nElem),source=0)
rewind(fileUnit)
read (fileUnit,'(A300)',END=630) line
do
chunkPos = IO_stringPos(line)
if( (IO_lc(IO_stringValue(line,chunkPos,1)) == 'initial') .and. &
(IO_lc(IO_stringValue(line,chunkPos,2)) == 'state') ) then
if (initialcondTableStyle == 2) read (fileUnit,'(A300)',END=630) line ! read extra line for new style
read (fileUnit,'(A300)',END=630) line ! read line with index of state var
chunkPos = IO_stringPos(line)
sv = IO_IntValue(line,chunkPos,1) ! figure state variable index
if( (sv == 2).or.(sv == 3) ) then ! only state vars 2 and 3 of interest
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_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
endif
contInts = IO_continuousIntValues& ! get affected elements
(fileUnit,nElem,nameElemSet,mapElemSet,size(nameElemSet))
do l = 1, size(fileContent)
chunkPos = IO_stringPos(fileContent(l))
if( (IO_lc(IO_stringValue(fileContent(l),chunkPos,1)) == 'initial') .and. &
(IO_lc(IO_stringValue(fileContent(l),chunkPos,2)) == 'state') ) then
k = merge(2,1,initialcondTableStyle == 2)
chunkPos = IO_stringPos(fileContent(l+k))
sv = IO_IntValue(fileContent(l+k),chunkPos,1) ! figure state variable index
if( (sv == 2) .or. (sv == 3) ) then ! only state vars 2 and 3 of interest
m = 1
chunkPos = IO_stringPos(fileContent(l+k+m))
do while (scan(IO_stringValue(fileContent(l+k+m),chunkPos,1),'+-',back=.true.)>1) ! is noEfloat value?
myVal = nint(IO_floatValue(fileContent(l+k+m),chunkPos,1))
if (initialcondTableStyle == 2) m = m + 2
contInts = continuousIntValues(fileContent(l+k+m+1:),nElem,nameElemSet,mapElemSet,size(nameElemSet)) ! get affected elements
do i = 1,contInts(1)
e = mesh_FEasCP('elem',contInts(1+i))
if (sv == 2) microstructureAt(e) = myVal
if (sv == 3) homogenizationAt(e) = myVal
enddo
if (initialcondTableStyle == 0) read (fileUnit,'(A300)',END=630) line ! ignore IP range for old table style
read (fileUnit,'(A300)',END=630) line
chunkPos = IO_stringPos(line)
if (initialcondTableStyle == 0) m = m + 1
enddo
endif
else
read (fileUnit,'(A300)',END=630) line
endif
enddo
630 end subroutine inputRead_microstructureAndHomogenization
end subroutine inputRead_microstructureAndHomogenization
!--------------------------------------------------------------------------------------------------
@ -1092,4 +1067,62 @@ integer function mesh_FEasCP(what,myID)
end function mesh_FEasCP
!--------------------------------------------------------------------------------------------------
!> @brief return integer list corresponding to items in consecutive lines.
!! First integer in array is counter
!> @details ints concatenated by "c" as last char, range of a "to" b, or named set
!--------------------------------------------------------------------------------------------------
function continuousIntValues(fileContent,maxN,lookupName,lookupMap,lookupMaxN)
character(len=*), dimension(:), intent(in) :: fileContent !< file content, separated per lines
integer, intent(in) :: maxN
integer, intent(in) :: lookupMaxN
integer, dimension(:,:), intent(in) :: lookupMap
character(len=*), dimension(:), intent(in) :: lookupName
integer, dimension(1+maxN) :: continuousIntValues
integer :: l,i,first,last
integer, allocatable, dimension(:) :: chunkPos
logical :: rangeGeneration
continuousIntValues = 0
rangeGeneration = .false.
do l = 1, size(fileContent)
chunkPos = IO_stringPos(fileContent(l))
if (chunkPos(1) < 1) then ! empty line
exit
elseif (verify(IO_stringValue(fileContent(l),chunkPos,1),'0123456789') > 0) then ! a non-int, i.e. set name
do i = 1, lookupMaxN ! loop over known set names
if (IO_stringValue(fileContent(l),chunkPos,1) == lookupName(i)) then ! found matching name
continuousIntValues = lookupMap(:,i) ! return resp. entity list
exit
endif
enddo
exit
else if (chunkPos(1) > 2 .and. IO_lc(IO_stringValue(fileContent(l),chunkPos,2)) == 'to' ) then ! found range indicator
first = IO_intValue(fileContent(l),chunkPos,1)
last = IO_intValue(fileContent(l),chunkPos,3)
do i = first, last, sign(1,last-first)
continuousIntValues(1) = continuousIntValues(1) + 1
continuousIntValues(1+continuousIntValues(1)) = i
enddo
exit
else
do i = 1,chunkPos(1)-1 ! interpret up to second to last value
continuousIntValues(1) = continuousIntValues(1) + 1
continuousIntValues(1+continuousIntValues(1)) = IO_intValue(fileContent(l),chunkPos,i)
enddo
if ( IO_lc(IO_stringValue(fileContent(l),chunkPos,chunkPos(1))) /= 'c' ) then ! line finished, read last value
continuousIntValues(1) = continuousIntValues(1) + 1
continuousIntValues(1+continuousIntValues(1)) = IO_intValue(fileContent(l),chunkPos,chunkPos(1))
exit
endif
endif
enddo
end function continuousIntValues
end module mesh

12
src/results.f90
View File

@ -487,9 +487,9 @@ end subroutine results_writeScalarDataset_rotation
!--------------------------------------------------------------------------------------------------
subroutine results_mapping_constituent(phaseAt,memberAt,label)
integer, dimension(:,:), intent(in) :: phaseAt !< phase section at (constituent,element)
integer, dimension(:,:,:), intent(in) :: memberAt !< phase member at (constituent,IP,element)
character(len=pStringLen), dimension(:), intent(in) :: label !< label of each phase section
integer, dimension(:,:), intent(in) :: phaseAt !< phase section at (constituent,element)
integer, dimension(:,:,:), intent(in) :: memberAt !< phase member at (constituent,IP,element)
character(len=*), dimension(:), intent(in) :: label !< label of each phase section
integer, dimension(size(memberAt,1),size(memberAt,2),size(memberAt,3)) :: &
phaseAt_perIP, &
@ -622,9 +622,9 @@ end subroutine results_mapping_constituent
!--------------------------------------------------------------------------------------------------
subroutine results_mapping_materialpoint(homogenizationAt,memberAt,label)
integer, dimension(:), intent(in) :: homogenizationAt !< homogenization section at (element)
integer, dimension(:,:), intent(in) :: memberAt !< homogenization member at (IP,element)
character(len=pStringLen), dimension(:), intent(in) :: label !< label of each homogenization section
integer, dimension(:), intent(in) :: homogenizationAt !< homogenization section at (element)
integer, dimension(:,:), intent(in) :: memberAt !< homogenization member at (IP,element)
character(len=*), dimension(:), intent(in) :: label !< label of each homogenization section
integer, dimension(size(memberAt,1),size(memberAt,2)) :: &
homogenizationAt_perIP, &

4
src/thermal_adiabatic.f90
View File

@ -77,7 +77,7 @@ subroutine thermal_adiabatic_init
allocate(thermalState(h)%state (1,NofMyHomog), source=thermal_initialT(h))
nullify(thermalMapping(h)%p)
thermalMapping(h)%p => mappingHomogenization(1,:,:)
thermalMapping(h)%p => material_homogenizationMemberAt
deallocate(temperature(h)%p)
temperature(h)%p => thermalState(h)%state(1,:)
deallocate(temperatureRate(h)%p)
@ -109,7 +109,7 @@ function thermal_adiabatic_updateState(subdt, ip, el)
T, Tdot, dTdot_dT
homog = material_homogenizationAt(el)
offset = mappingHomogenization(1,ip,el)
offset = material_homogenizationMemberAt(ip,el)
T = thermalState(homog)%subState0(1,offset)
call thermal_adiabatic_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el)

4
src/thermal_conduction.f90
View File

@ -79,7 +79,7 @@ subroutine thermal_conduction_init
allocate(thermalState(h)%state (0,NofMyHomog))
nullify(thermalMapping(h)%p)
thermalMapping(h)%p => mappingHomogenization(1,:,:)
thermalMapping(h)%p => material_homogenizationMemberAt
deallocate(temperature (h)%p)
allocate (temperature (h)%p(NofMyHomog), source=thermal_initialT(h))
deallocate(temperatureRate(h)%p)
@ -114,7 +114,7 @@ subroutine thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el)
constituent
homog = material_homogenizationAt(el)
offset = mappingHomogenization(1,ip,el)
offset = material_homogenizationMemberAt(ip,el)
instance = thermal_typeInstance(homog)
Tdot = 0.0_pReal