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

Merge branch 'development' into more_crystallite_cleaning

This commit is contained in:
Franz Roters 2020-10-30 16:42:19 +01:00
parent 0c51c3d8d4 0af81e1af7
commit 98cb7be21a
61 changed files with 753 additions and 740 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
PRIVATE

@ -1 +1 @@
Subproject commit 1e8c66897820468ab46958d995005e2b69204d0e Subproject commit 3112a4dbfa1e926c07b7f9443161239b8a7e85ca

2
VERSION
View File

@ -1 +1 @@
v3.0.0-alpha-553-g5eee0d74f v3.0.0-alpha-594-g46e5023f8

3
examples/.gitignore vendored
View File

@ -3,6 +3,3 @@
*.xdmf *.xdmf
*.sta *.sta
*.vt* *.vt*
*.geom
*.seeds
postProc

21
python/damask/_config.py
View File

@ -1,6 +1,7 @@
from io import StringIO from io import StringIO
import abc import abc
import numpy as np
import yaml import yaml
class NiceDumper(yaml.SafeDumper): class NiceDumper(yaml.SafeDumper):
@ -15,6 +16,11 @@ class NiceDumper(yaml.SafeDumper):
def increase_indent(self, flow=False, indentless=False): def increase_indent(self, flow=False, indentless=False):
return super().increase_indent(flow, False) return super().increase_indent(flow, False)
def represent_data(self, data):
"""Cast Config objects and its subclasses to dict."""
return self.represent_data(dict(data)) if isinstance(data, dict) and type(data) != dict else \
super().represent_data(data)
class Config(dict): class Config(dict):
"""YAML-based configuration.""" """YAML-based configuration."""
@ -64,7 +70,20 @@ class Config(dict):
kwargs['width'] = 256 kwargs['width'] = 256
if 'default_flow_style' not in kwargs: if 'default_flow_style' not in kwargs:
kwargs['default_flow_style'] = None kwargs['default_flow_style'] = None
fhandle.write(yaml.dump(dict(self),Dumper=NiceDumper,**kwargs)) if 'sort_keys' not in kwargs:
kwargs['sort_keys'] = False
def array_representer(dumper, data):
"""Convert numpy array to list of native types."""
return dumper.represent_list([d.item() for d in data])
NiceDumper.add_representer(np.ndarray, array_representer)
try:
fhandle.write(yaml.dump(self,Dumper=NiceDumper,**kwargs))
except TypeError: # compatibility with old pyyaml
del kwargs['sort_keys']
fhandle.write(yaml.dump(self,Dumper=NiceDumper,**kwargs))
@property @property

38
python/damask/_configmaterial.py
View File

@ -2,7 +2,6 @@ import copy
import numpy as np import numpy as np
from . import grid_filters
from . import Config from . import Config
from . import Lattice from . import Lattice
from . import Rotation from . import Rotation
@ -26,7 +25,7 @@ class ConfigMaterial(Config):
@staticmethod @staticmethod
def from_table(table,coordinates=None,constituents={},**kwargs): def from_table(table,constituents={},**kwargs):
""" """
Load from an ASCII table. Load from an ASCII table.
@ -34,10 +33,6 @@ class ConfigMaterial(Config):
---------- ----------
table : damask.Table table : damask.Table
Table that contains material information. Table that contains material information.
coordinates : str, optional
Label of spatial coordiates. Used for sorting and performing a
sanity check. Default to None, in which case no sorting or checking is
peformed.
constituents : dict, optional constituents : dict, optional
Entries for 'constituents'. The key is the name and the value specifies Entries for 'constituents'. The key is the name and the value specifies
the label of the data column in the table the label of the data column in the table
@ -54,7 +49,7 @@ class ConfigMaterial(Config):
pos pos pos qu qu qu qu phase homog pos pos pos qu qu qu qu phase homog
0 0 0 0 0.19 0.8 0.24 -0.51 Aluminum SX 0 0 0 0 0.19 0.8 0.24 -0.51 Aluminum SX
1 1 0 0 0.8 0.19 0.24 -0.51 Steel SX 1 1 0 0 0.8 0.19 0.24 -0.51 Steel SX
>>> cm.from_table(t,'pos',{'O':'qu','phase':'phase'},homogenization='homog') >>> cm.from_table(t,{'O':'qu','phase':'phase'},homogenization='homog')
material: material:
- constituents: - constituents:
- O: [0.19, 0.8, 0.24, -0.51] - O: [0.19, 0.8, 0.24, -0.51]
@ -68,17 +63,12 @@ class ConfigMaterial(Config):
homogenization: SX homogenization: SX
""" """
if coordinates is not None: constituents_ = {k:table.get(v) for k,v in constituents.items()}
t = table.sort_by([f'{i}_{coordinates}' for i in range(3,0,-1)]) kwargs_ = {k:table.get(v) for k,v in kwargs.items()}
grid_filters.coord0_check(t.get(coordinates))
else:
t = table
constituents_ = {k:t.get(v) for k,v in constituents.items()}
kwargs_ = {k:t.get(v) for k,v in kwargs.items()}
_,idx = np.unique(np.hstack(list({**constituents_,**kwargs_}.values())),return_index=True,axis=0) _,idx = np.unique(np.hstack(list({**constituents_,**kwargs_}.values())),return_index=True,axis=0)
idx = np.sort(idx)
constituents_ = {k:v[idx].squeeze() for k,v in constituents_.items()} constituents_ = {k:v[idx].squeeze() for k,v in constituents_.items()}
kwargs_ = {k:v[idx].squeeze() for k,v in kwargs_.items()} kwargs_ = {k:v[idx].squeeze() for k,v in kwargs_.items()}
@ -229,19 +219,18 @@ class ConfigMaterial(Config):
Parameters Parameters
---------- ----------
constituents : dict constituents : dict
Entries for 'constituents'. The key is the name and the value specifies Entries for 'constituents' as key-value pair.
the label of the data column in the table
**kwargs **kwargs
Keyword arguments where the key is the name and the value specifies Key-value pairs.
the label of the data column in the table
Examples Examples
-------- --------
>>> import damask >>> import damask
>>> m = damask.ConfigMaterial()
>>> O = damask.Rotation.from_random(3).as_quaternion() >>> O = damask.Rotation.from_random(3).as_quaternion()
>>> phase = ['Aluminum','Steel','Aluminum'] >>> phase = ['Aluminum','Steel','Aluminum']
>>> m.material_add(constituents={'phase':phase,'O':O},homogenization='SX') >>> m = damask.ConfigMaterial().material_add(constituents={'phase':phase,'O':O},
... homogenization='SX')
>>> m
material: material:
- constituents: - constituents:
- O: [0.577764, -0.146299, -0.617669, 0.513010] - O: [0.577764, -0.146299, -0.617669, 0.513010]
@ -264,15 +253,14 @@ class ConfigMaterial(Config):
for k,v in kwargs.items(): for k,v in kwargs.items():
if hasattr(v,'__len__') and not isinstance(v,str): if hasattr(v,'__len__') and not isinstance(v,str):
if len(v) != len(c): if len(v) != len(c):
raise ValueError('len mismatch 1') raise ValueError('Cannot add entries of different length')
for i,vv in enumerate(v): for i,vv in enumerate(v):
c[i][k] = [w.item() for w in vv] if isinstance(vv,np.ndarray) else vv.item() c[i][k] = [w.item() for w in vv] if isinstance(vv,np.ndarray) else vv.item()
else: else:
for i in range(len(c)): for i in range(len(c)):
c[i][k] = v c[i][k] = v
dup = copy.deepcopy(self) dup = copy.deepcopy(self)
if 'material' not in dup: dup['material'] = [] dup['material'] = dup['material'] + c if 'material' in dup else c
dup['material'] +=c
return dup return dup
@ -288,7 +276,7 @@ class ConfigMaterial(Config):
for k,v in kwargs.items(): for k,v in kwargs.items():
if hasattr(v,'__len__') and not isinstance(v,str): if hasattr(v,'__len__') and not isinstance(v,str):
if len(v) != N_material: if len(v) != N_material:
raise ValueError('len mismatch 2') raise ValueError('Cannot add entries of different length')
for i,vv in enumerate(np.array(v)): for i,vv in enumerate(np.array(v)):
m[i][0][k] = [w.item() for w in vv] if isinstance(vv,np.ndarray) else vv.item() m[i][0][k] = [w.item() for w in vv] if isinstance(vv,np.ndarray) else vv.item()
else: else:

44
python/damask/_geom.py
View File

@ -4,6 +4,7 @@ from functools import partial
from os import path from os import path
import numpy as np import numpy as np
import pandas as pd
import h5py import h5py
from scipy import ndimage,spatial from scipy import ndimage,spatial
@ -254,21 +255,26 @@ class Geom:
table : damask.Table table : damask.Table
Table that contains material information. Table that contains material information.
coordinates : str coordinates : str
Label of the column containing the spatial coordinates. Label of the column containing the vector of spatial coordinates.
Need to be ordered (1./x fast, 3./z slow).
labels : str or list of str labels : str or list of str
Label(s) of the columns containing the material definition. Label(s) of the columns containing the material definition.
Each unique combintation of values results in a material. Each unique combintation of values results in a material.
""" """
t = table.sort_by([f'{i}_{coordinates}' for i in range(3,0,-1)]) grid,size,origin = grid_filters.cell_coord0_gridSizeOrigin(table.get(coordinates))
grid,size,origin = grid_filters.cell_coord0_gridSizeOrigin(t.get(coordinates))
labels_ = [labels] if isinstance(labels,str) else labels labels_ = [labels] if isinstance(labels,str) else labels
_,unique_inverse = np.unique(np.hstack([t.get(l) for l in labels_]),return_inverse=True,axis=0) unique,unique_inverse = np.unique(np.hstack([table.get(l) for l in labels_]),return_inverse=True,axis=0)
ma = unique_inverse.reshape(grid,order='F') + 1 if len(unique) == grid.prod():
ma = np.arange(grid.prod())
else:
from_ma = pd.unique(unique_inverse)
sort_idx = np.argsort(from_ma)
idx = np.searchsorted(from_ma,unique_inverse,sorter = sort_idx)
ma = np.arange(from_ma.size)[sort_idx][idx]
return Geom(ma,size,origin,util.execution_stamp('Geom','from_table')) return Geom(ma.reshape(grid,order='F'),size,origin,util.execution_stamp('Geom','from_table'))
@staticmethod @staticmethod
@ -417,7 +423,7 @@ class Geom:
Number of periods per unit cell. Defaults to 1. Number of periods per unit cell. Defaults to 1.
materials : (int, int), optional materials : (int, int), optional
Material IDs. Defaults to (1,2). Material IDs. Defaults to (1,2).
Notes Notes
----- -----
The following triply-periodic minimal surfaces are implemented: The following triply-periodic minimal surfaces are implemented:
@ -473,7 +479,7 @@ class Geom:
Compress with zlib algorithm. Defaults to True. Compress with zlib algorithm. Defaults to True.
""" """
v = VTK.from_rectilinearGrid(self.grid,self.size,self.origin) v = VTK.from_rectilinear_grid(self.grid,self.size,self.origin)
v.add(self.material.flatten(order='F'),'material') v.add(self.material.flatten(order='F'),'material')
v.add_comments(self.comments) v.add_comments(self.comments)
@ -508,7 +514,7 @@ class Geom:
def show(self): def show(self):
"""Show on screen.""" """Show on screen."""
v = VTK.from_rectilinearGrid(self.grid,self.size,self.origin) v = VTK.from_rectilinear_grid(self.grid,self.size,self.origin)
v.show() v.show()
@ -687,12 +693,10 @@ class Geom:
def renumber(self): def renumber(self):
"""Renumber sorted material indices to 1,...,N.""" """Renumber sorted material indices to 0,...,N-1."""
renumbered = np.empty(self.grid,dtype=self.material.dtype) _,renumbered = np.unique(self.material,return_inverse=True)
for i, oldID in enumerate(np.unique(self.material)):
renumbered = np.where(self.material == oldID, i+1, renumbered)
return Geom(material = renumbered, return Geom(material = renumbered.reshape(self.grid),
size = self.size, size = self.size,
origin = self.origin, origin = self.origin,
comments = self.comments+[util.execution_stamp('Geom','renumber')], comments = self.comments+[util.execution_stamp('Geom','renumber')],
@ -783,11 +787,13 @@ class Geom:
New material indices. New material indices.
""" """
substituted = self.material.copy() def mp(entry,mapper):
for from_ms,to_ms in zip(from_material,to_material): return mapper[entry] if entry in mapper else entry
substituted[self.material==from_ms] = to_ms
mp = np.vectorize(mp)
mapper = dict(zip(from_material,to_material))
return Geom(material = substituted, return Geom(material = mp(self.material,mapper).reshape(self.grid),
size = self.size, size = self.size,
origin = self.origin, origin = self.origin,
comments = self.comments+[util.execution_stamp('Geom','substitute')], comments = self.comments+[util.execution_stamp('Geom','substitute')],

10
python/damask/_result.py
View File

@ -1212,14 +1212,14 @@ class Result:
if mode.lower()=='cell': if mode.lower()=='cell':
if self.structured: if self.structured:
v = VTK.from_rectilinearGrid(self.grid,self.size,self.origin) v = VTK.from_rectilinear_grid(self.grid,self.size,self.origin)
else: else:
with h5py.File(self.fname,'r') as f: with h5py.File(self.fname,'r') as f:
v = VTK.from_unstructuredGrid(f['/geometry/x_n'][()], v = VTK.from_unstructured_grid(f['/geometry/x_n'][()],
f['/geometry/T_c'][()]-1, f['/geometry/T_c'][()]-1,
f['/geometry/T_c'].attrs['VTK_TYPE'].decode()) f['/geometry/T_c'].attrs['VTK_TYPE'].decode())
elif mode.lower()=='point': elif mode.lower()=='point':
v = VTK.from_polyData(self.cell_coordinates()) v = VTK.from_poly_data(self.cell_coordinates)
N_digits = int(np.floor(np.log10(max(1,int(self.increments[-1][3:])))))+1 N_digits = int(np.floor(np.log10(max(1,int(self.increments[-1][3:])))))+1

22
python/damask/_vtk.py
View File

@ -36,7 +36,7 @@ class VTK:
@staticmethod @staticmethod
def from_rectilinearGrid(grid,size,origin=np.zeros(3)): def from_rectilinear_grid(grid,size,origin=np.zeros(3)):
""" """
Create VTK of type vtk.vtkRectilinearGrid. Create VTK of type vtk.vtkRectilinearGrid.
@ -64,7 +64,7 @@ class VTK:
@staticmethod @staticmethod
def from_unstructuredGrid(nodes,connectivity,cell_type): def from_unstructured_grid(nodes,connectivity,cell_type):
""" """
Create VTK of type vtk.vtkUnstructuredGrid. Create VTK of type vtk.vtkUnstructuredGrid.
@ -96,7 +96,7 @@ class VTK:
@staticmethod @staticmethod
def from_polyData(points): def from_poly_data(points):
""" """
Create VTK of type vtk.polyData. Create VTK of type vtk.polyData.
@ -108,11 +108,18 @@ class VTK:
Spatial position of the points. Spatial position of the points.
""" """
N = points.shape[0]
vtk_points = vtk.vtkPoints() vtk_points = vtk.vtkPoints()
vtk_points.SetData(np_to_vtk(points)) vtk_points.SetData(np_to_vtk(points))
vtk_cells = vtk.vtkCellArray()
vtk_cells.SetNumberOfCells(N)
vtk_cells.SetCells(N,np_to_vtkIdTypeArray(np.stack((np.ones (N,dtype=np.int64),
np.arange(N,dtype=np.int64)),axis=1).ravel(),deep=True))
vtk_data = vtk.vtkPolyData() vtk_data = vtk.vtkPolyData()
vtk_data.SetPoints(vtk_points) vtk_data.SetPoints(vtk_points)
vtk_data.SetVerts(vtk_cells)
return VTK(vtk_data) return VTK(vtk_data)
@ -164,6 +171,7 @@ class VTK:
return VTK(vtk_data) return VTK(vtk_data)
@staticmethod @staticmethod
def _write(writer): def _write(writer):
"""Wrapper for parallel writing.""" """Wrapper for parallel writing."""
@ -192,7 +200,7 @@ class VTK:
default_ext = writer.GetDefaultFileExtension() default_ext = writer.GetDefaultFileExtension()
ext = Path(fname).suffix ext = Path(fname).suffix
if ext and ext != '.'+default_ext: if ext and ext != '.'+default_ext:
raise ValueError(f'Given extension {ext} does not match default .{default_ext}') raise ValueError(f'Given extension "{ext}" does not match default ".{default_ext}"')
writer.SetFileName(str(Path(fname).with_suffix('.'+default_ext))) writer.SetFileName(str(Path(fname).with_suffix('.'+default_ext)))
if compress: if compress:
writer.SetCompressorTypeToZLib() writer.SetCompressorTypeToZLib()
@ -238,10 +246,10 @@ class VTK:
else data).reshape(N_data,-1),deep=True) # avoid large files else data).reshape(N_data,-1),deep=True) # avoid large files
d.SetName(label) d.SetName(label)
if N_data == N_cells: if N_data == N_points:
self.vtk_data.GetCellData().AddArray(d)
elif N_data == N_points:
self.vtk_data.GetPointData().AddArray(d) self.vtk_data.GetPointData().AddArray(d)
elif N_data == N_cells:
self.vtk_data.GetCellData().AddArray(d)
else: else:
raise ValueError(f'Cell / point count ({N_cells} / {N_points}) differs from data ({N_data}).') raise ValueError(f'Cell / point count ({N_cells} / {N_points}) differs from data ({N_data}).')
elif isinstance(data,pd.DataFrame): elif isinstance(data,pd.DataFrame):

61
python/tests/reference/Result/12grains6x7x8.geom
View File

@ -1,61 +0,0 @@
4 header
grid a 6 b 7 c 8
size x 0.75 y 0.875 z 1.0
origin x 0.0 y 0.0 z 0.0
homogenization 1
9 3 3 10 9 9
9 1 1 1 9 9
9 11 1 1 7 9
7 11 11 7 7 7
7 11 11 7 7 7
12 3 3 10 7 12
12 3 3 10 10 12
12 3 3 1 9 9
9 1 1 1 9 9
9 1 1 1 7 7
7 1 1 7 7 7
12 12 3 7 7 7
12 3 3 3 12 12
12 3 3 3 12 12
12 3 3 1 1 12
9 1 1 1 1 9
6 1 1 1 8 8
7 6 8 8 8 8
12 12 8 8 8 12
12 3 3 3 12 12
12 3 3 3 12 12
5 6 6 6 1 12
6 6 6 6 8 8
6 6 6 8 8 8
8 6 8 8 8 8
12 5 8 8 8 8
12 5 5 8 8 12
5 5 5 3 12 12
5 5 6 6 6 5
6 6 6 6 6 6
6 6 6 6 8 8
4 4 6 8 8 8
4 4 2 2 2 8
5 5 5 2 2 2
5 5 5 5 2 5
5 5 5 10 10 5
6 6 6 6 10 4
4 4 11 11 2 4
4 4 11 2 2 4
4 4 2 2 2 2
5 5 5 2 2 2
5 5 5 10 10 5
5 5 10 10 10 9
4 11 11 11 10 9
4 4 11 11 11 4
4 4 11 11 2 4
4 4 2 2 2 2
5 5 2 2 2 2
5 5 10 10 10 10
9 10 10 10 10 9
9 11 11 10 9 9
4 11 11 11 9 9
4 11 11 11 7 7
4 4 11 2 7 7
12 10 10 10 10 7
9 10 10 10 10 9

30
python/tests/reference/Result/12grains6x7x8.vtr Normal file
View File

@ -0,0 +1,30 @@
<?xml version="1.0"?>
<VTKFile type="RectilinearGrid" version="0.1" byte_order="LittleEndian" header_type="UInt32" compressor="vtkZLibDataCompressor">
<RectilinearGrid WholeExtent="0 6 0 7 0 8">
<FieldData>
<Array type="String" Name="comments" NumberOfTuples="1" format="binary">
AQAAAACAAAARAAAAGQAAAA==eF7LyM/NT0/Ny6xKLMnMz1MwZAAAPsIGPQ==
</Array>
</FieldData>
<Piece Extent="0 6 0 7 0 8">
<PointData>
</PointData>
<CellData>
<DataArray type="Int64" Name="material" format="binary" RangeMin="0" RangeMax="11">
AQAAAACAAACACgAA2wAAAA==eF6tlssOgkAQBNcHiPj//2uM25eOZc8oc6kYagzdISzbeM/ZeJ/cgDTk7x/c16zm6fduXAOr/mOS8rqXfDH5mqP6pKHcVY9yJN/ziv5/R/k+8lI/GnnLV2uMm1G5tefXnZ6j6v3bD/nXyQWokU8e5a96tJc8z9H1aY88MfWZek3Xf6XnuBhTr+6fgPKpH9oj3/eS5+/bap/yPafo54buJ/mek3zqJeXu+tRP8vycp15E8tNe6rPaf7fPrk/9eO6q598/1GO1/67v53V6nul8T3n9Oy758p47Sgdl
</DataArray>
</CellData>
<Coordinates>
<DataArray type="Float64" Name="x" format="binary" RangeMin="0" RangeMax="0.75">
AQAAAACAAAA4AAAAHAAAAA==eF5jYEAGB+wh9AUofQNKP4DST6D0C3sAisIGjw==
</DataArray>
<DataArray type="Float64" Name="y" format="binary" RangeMin="0" RangeMax="0.875">
AQAAAACAAABAAAAAHwAAAA==eF5jYEAGB+wh9AUofQNKP4DST6D0Cyj9xh4AwVEHug==
</DataArray>
<DataArray type="Float64" Name="z" format="binary" RangeMin="0" RangeMax="1">
AQAAAACAAABIAAAAIgAAAA==eF5jYEAGB+wh9AUofQNKP4DST6D0Cyj9Bkp/sAcAAU8I6Q==
</DataArray>
</Coordinates>
</Piece>
</RectilinearGrid>
</VTKFile>

61
python/tests/reference/Result/6grains6x7x8_single_phase.geom
View File

@ -1,61 +0,0 @@
4 header
grid a 6 b 7 c 8
size x 0.75 y 0.875 z 1.0
origin x 0.0 y 0.0 z 0.0
homogenization 1
3 3 3 4 3 3
3 1 1 1 3 3
3 5 1 1 1 3
1 5 5 1 1 1
1 5 5 1 1 1
6 3 3 4 1 6
6 3 3 4 4 6
6 3 3 1 3 3
3 1 1 1 3 3
3 1 1 1 1 1
1 1 1 1 1 1
6 6 3 1 1 1
6 3 3 3 6 6
6 3 3 3 6 6
6 3 3 1 1 6
3 1 1 1 1 3
6 1 1 1 2 2
1 6 2 2 2 2
6 6 2 2 2 6
6 3 3 3 6 6
6 3 3 3 6 6
5 6 6 6 1 6
6 6 6 6 2 2
6 6 6 2 2 2
2 6 2 2 2 2
6 5 2 2 2 2
6 5 5 2 2 6
5 5 5 3 6 6
5 5 6 6 6 5
6 6 6 6 6 6
6 6 6 6 2 2
4 4 6 2 2 2
4 4 2 2 2 2
5 5 5 2 2 2
5 5 5 5 2 5
5 5 5 4 4 5
6 6 6 6 4 4
4 4 5 5 2 4
4 4 5 2 2 4
4 4 2 2 2 2
5 5 5 2 2 2
5 5 5 4 4 5
5 5 4 4 4 3
4 5 5 5 4 3
4 4 5 5 5 4
4 4 5 5 2 4
4 4 2 2 2 2
5 5 2 2 2 2
5 5 4 4 4 4
3 4 4 4 4 3
3 5 5 4 3 3
4 5 5 5 3 3
4 5 5 5 1 1
4 4 5 2 1 1
6 4 4 4 4 1
3 4 4 4 4 3

30
python/tests/reference/Result/6grains6x7x8_single_phase.vtr Normal file
View File

@ -0,0 +1,30 @@
<?xml version="1.0"?>
<VTKFile type="RectilinearGrid" version="0.1" byte_order="LittleEndian" header_type="UInt32" compressor="vtkZLibDataCompressor">
<RectilinearGrid WholeExtent="0 6 0 7 0 8">
<FieldData>
<Array type="String" Name="comments" NumberOfTuples="1" format="binary">
AQAAAACAAAARAAAAGQAAAA==eF7LyM/NT0/Ny6xKLMnMz1MwZAAAPsIGPQ==
</Array>
</FieldData>
<Piece Extent="0 6 0 7 0 8">
<PointData>
</PointData>
<CellData>
<DataArray type="Int64" Name="material" format="binary" RangeMin="0" RangeMax="5">
AQAAAACAAACACgAAwAAAAA==eF69lcsOwjAQA6GU//9lDsUSGjHyBgq+WGpm09jqY7sc2uA3uR43Gb+/YV/FfXd405S/P93ykmt8vPHRWX3+SpbDZHnj3O8snpqeN+L9TFd4lDmu05ljyn3bj/F5P2yfyNZbnilnc1MuOVZ5mzMu3vpsvbb1T5057Ltk/ZBvfVo/qzznGsdzTvvknO3D8zS+9fjvPlsu4/ifn87bf7DNNf7sPlf59rxYbuPircdp/6s81Z5navoeRav9PACxsANv
</DataArray>
</CellData>
<Coordinates>
<DataArray type="Float64" Name="x" format="binary" RangeMin="0" RangeMax="0.75">
AQAAAACAAAA4AAAAHAAAAA==eF5jYEAGB+wh9AUofQNKP4DST6D0C3sAisIGjw==
</DataArray>
<DataArray type="Float64" Name="y" format="binary" RangeMin="0" RangeMax="0.875">
AQAAAACAAABAAAAAHwAAAA==eF5jYEAGB+wh9AUofQNKP4DST6D0Cyj9xh4AwVEHug==
</DataArray>
<DataArray type="Float64" Name="z" format="binary" RangeMin="0" RangeMax="1">
AQAAAACAAABIAAAAIgAAAA==eF5jYEAGB+wh9AUofQNKP4DST6D0Cyj9Bkp/sAcAAU8I6Q==
</DataArray>
</Coordinates>
</Piece>
</RectilinearGrid>
</VTKFile>

10
python/tests/reference/VTK/polyData.vtp
View File

@ -1,7 +1,7 @@
<?xml version="1.0"?> <?xml version="1.0"?>
<VTKFile type="PolyData" version="0.1" byte_order="LittleEndian" header_type="UInt32" compressor="vtkZLibDataCompressor"> <VTKFile type="PolyData" version="0.1" byte_order="LittleEndian" header_type="UInt32" compressor="vtkZLibDataCompressor">
<PolyData> <PolyData>
<Piece NumberOfPoints="10" NumberOfVerts="0" NumberOfLines="0" NumberOfStrips="0" NumberOfPolys="0"> <Piece NumberOfPoints="10" NumberOfVerts="10" NumberOfLines="0" NumberOfStrips="0" NumberOfPolys="0">
<PointData> <PointData>
<DataArray type="Float32" Name="coordinates" NumberOfComponents="3" format="binary" RangeMin="0.7453560147132696" RangeMax="2.449489742783178"> <DataArray type="Float32" Name="coordinates" NumberOfComponents="3" format="binary" RangeMin="0.7453560147132696" RangeMax="2.449489742783178">
AQAAAACAAAB4AAAAVgAAAA==eF5jYICBhv2WfY9tLfuS7Ypk3PeDaCDf7okF3/7Vq1bZrV6lZQ+k94HEgHL2QHovUM7+iUUfiG0LlQdhkH77Ipnj9iB5qFp7kBjQDiBmcADRANsaLXM= AQAAAACAAAB4AAAAVgAAAA==eF5jYICBhv2WfY9tLfuS7Ypk3PeDaCDf7okF3/7Vq1bZrV6lZQ+k94HEgHL2QHovUM7+iUUfiG0LlQdhkH77Ipnj9iB5qFp7kBjQDiBmcADRANsaLXM=
@ -31,11 +31,11 @@
</DataArray> </DataArray>
</Points> </Points>
<Verts> <Verts>
<DataArray type="Int64" Name="connectivity" format="binary" RangeMin="1e+299" RangeMax="-1e+299"> <DataArray type="Int64" Name="connectivity" format="binary" RangeMin="0" RangeMax="9">
AAAAAACAAAAAAAAA AQAAAACAAABQAAAAIgAAAA==eF4txbcBACAIADAsiP7/sAPJkog2PL28nT4uXz9/BXgALg==
</DataArray> </DataArray>
<DataArray type="Int64" Name="offsets" format="binary" RangeMin="1e+299" RangeMax="-1e+299"> <DataArray type="Int64" Name="offsets" format="binary" RangeMin="1" RangeMax="10">
AAAAAACAAAAAAAAA AQAAAACAAABQAAAAIgAAAA==eF4txbcBACAIADA76v8HM5As6a0MTy9vH4evn78TBzAAOA==
</DataArray> </DataArray>
</Verts> </Verts>
<Lines> <Lines>

3
python/tests/test_Config.py
View File

@ -1,4 +1,5 @@
import pytest import pytest
import numpy as np
from damask import Config from damask import Config
@ -29,6 +30,8 @@ class TestConfig:
f.write(config.__repr__()) f.write(config.__repr__())
assert Config.load(tmp_path/'config.yaml') == config assert Config.load(tmp_path/'config.yaml') == config
def test_numpy(self,tmp_path):
assert Config({'A':np.ones(3,'i')}).__repr__() == Config({'A':[1,1,1]}).__repr__()
def test_abstract_is_valid(self): def test_abstract_is_valid(self):
assert Config().is_valid is None assert Config().is_valid is None

13
python/tests/test_ConfigMaterial.py
View File

@ -1,8 +1,10 @@
import os import os
import pytest import pytest
import numpy as np
from damask import ConfigMaterial from damask import ConfigMaterial
from damask import Table
@pytest.fixture @pytest.fixture
def reference_dir(reference_dir_base): def reference_dir(reference_dir_base):
@ -74,3 +76,14 @@ class TestConfigMaterial:
material_config = ConfigMaterial.load(reference_dir/'material.yaml') material_config = ConfigMaterial.load(reference_dir/'material.yaml')
new = material_config.material_rename_homogenization({'Taylor':'isostrain'}) new = material_config.material_rename_homogenization({'Taylor':'isostrain'})
assert not new.is_complete assert not new.is_complete
def test_from_table(self):
N = np.random.randint(3,10)
a = np.vstack((np.hstack((np.arange(N),np.arange(N)[::-1])),np.ones(N*2),np.zeros(N*2),np.ones(N*2))).T
t = Table(a,{'varying':2,'constant':2})
c = ConfigMaterial.from_table(t,constituents={'a':'varying','b':'1_constant'},c='2_constant')
assert len(c['material']) == N
for i,m in enumerate(c['material']):
c = m['constituents'][0]
assert m['c'] == 1 and c['b'] == 0 and c['a'] == [i,1]

23
python/tests/test_Geom.py
View File

@ -3,8 +3,10 @@ import numpy as np
from damask import VTK from damask import VTK
from damask import Geom from damask import Geom
from damask import Table
from damask import Rotation from damask import Rotation
from damask import util from damask import util
from damask import grid_filters
def geom_equal(a,b): def geom_equal(a,b):
@ -49,7 +51,7 @@ class TestGeom:
assert geom_equal(new,default) assert geom_equal(new,default)
def test_invalid_vtr(self,tmp_path): def test_invalid_vtr(self,tmp_path):
v = VTK.from_rectilinearGrid(np.random.randint(5,10,3)*2,np.random.random(3) + 1.0) v = VTK.from_rectilinear_grid(np.random.randint(5,10,3)*2,np.random.random(3) + 1.0)
v.save(tmp_path/'no_materialpoint.vtr') v.save(tmp_path/'no_materialpoint.vtr')
with pytest.raises(ValueError): with pytest.raises(ValueError):
Geom.load(tmp_path/'no_materialpoint.vtr') Geom.load(tmp_path/'no_materialpoint.vtr')
@ -176,6 +178,7 @@ class TestGeom:
material = default.material.copy() material = default.material.copy()
for m in np.unique(material): for m in np.unique(material):
material[material==m] = material.max() + np.random.randint(1,30) material[material==m] = material.max() + np.random.randint(1,30)
default.material -= 1
modified = Geom(material, modified = Geom(material,
default.size, default.size,
default.origin) default.origin)
@ -194,6 +197,13 @@ class TestGeom:
modified.substitute(np.arange(default.material.max())+1+offset, modified.substitute(np.arange(default.material.max())+1+offset,
np.arange(default.material.max())+1)) np.arange(default.material.max())+1))
def test_substitute_invariant(self,default):
f = np.unique(default.material.flatten())[:np.random.randint(1,default.material.max())]
t = np.random.permutation(f)
modified = default.substitute(f,t)
assert np.array_equiv(t,f) or (not geom_equal(modified,default))
assert geom_equal(default, modified.substitute(t,f))
@pytest.mark.parametrize('axis_angle',[np.array([1,0,0,86.7]), np.array([0,1,0,90.4]), np.array([0,0,1,90]), @pytest.mark.parametrize('axis_angle',[np.array([1,0,0,86.7]), np.array([0,1,0,90.4]), np.array([0,0,1,90]),
np.array([1,0,0,175]),np.array([0,-1,0,178]),np.array([0,0,1,180])]) np.array([1,0,0,175]),np.array([0,-1,0,178]),np.array([0,0,1,180])])
@ -363,3 +373,14 @@ class TestGeom:
grid = np.ones(3,dtype=int)*64 grid = np.ones(3,dtype=int)*64
geom = Geom.from_minimal_surface(grid,np.ones(3),surface,threshold) geom = Geom.from_minimal_surface(grid,np.ones(3),surface,threshold)
assert np.isclose(np.count_nonzero(geom.material==1)/np.prod(geom.grid),.5,rtol=1e-3) assert np.isclose(np.count_nonzero(geom.material==1)/np.prod(geom.grid),.5,rtol=1e-3)
def test_from_table(self):
grid = np.random.randint(60,100,3)
size = np.ones(3)+np.random.rand(3)
coords = grid_filters.cell_coord0(grid,size).reshape(-1,3,order='F')
z=np.ones(grid.prod())
z[grid[:2].prod()*int(grid[2]/2):]=0
t = Table(np.column_stack((coords,z)),{'coords':3,'z':1})
g = Geom.from_table(t,'coords',['1_coords','z'])
assert g.N_materials == g.grid[0]*2 and (g.material[:,:,-1]-g.material[:,:,0] == grid[0]).all()

5
python/tests/test_Result.py
View File

@ -341,6 +341,11 @@ class TestResult:
os.chdir(tmp_path) os.chdir(tmp_path)
default.save_vtk(output) default.save_vtk(output)
@pytest.mark.parametrize('mode',['point','cell'])
def test_vtk_mode(self,tmp_path,single_phase,mode):
os.chdir(tmp_path)
single_phase.save_vtk(mode=mode)
def test_XDMF(self,tmp_path,single_phase): def test_XDMF(self,tmp_path,single_phase):
os.chdir(tmp_path) os.chdir(tmp_path)
single_phase.save_XDMF() single_phase.save_XDMF()

14
python/tests/test_VTK.py
View File

@ -18,7 +18,7 @@ def default():
"""Simple VTK.""" """Simple VTK."""
grid = np.array([5,6,7],int) grid = np.array([5,6,7],int)
size = np.array([.6,1.,.5]) size = np.array([.6,1.,.5])
return VTK.from_rectilinearGrid(grid,size) return VTK.from_rectilinear_grid(grid,size)
class TestVTK: class TestVTK:
@ -30,7 +30,7 @@ class TestVTK:
grid = np.random.randint(5,10,3)*2 grid = np.random.randint(5,10,3)*2
size = np.random.random(3) + 1.0 size = np.random.random(3) + 1.0
origin = np.random.random(3) origin = np.random.random(3)
v = VTK.from_rectilinearGrid(grid,size,origin) v = VTK.from_rectilinear_grid(grid,size,origin)
string = v.__repr__() string = v.__repr__()
v.save(tmp_path/'rectilinearGrid',False) v.save(tmp_path/'rectilinearGrid',False)
vtr = VTK.load(tmp_path/'rectilinearGrid.vtr') vtr = VTK.load(tmp_path/'rectilinearGrid.vtr')
@ -41,7 +41,7 @@ class TestVTK:
def test_polyData(self,tmp_path): def test_polyData(self,tmp_path):
points = np.random.rand(100,3) points = np.random.rand(100,3)
v = VTK.from_polyData(points) v = VTK.from_poly_data(points)
string = v.__repr__() string = v.__repr__()
v.save(tmp_path/'polyData',False) v.save(tmp_path/'polyData',False)
vtp = VTK.load(tmp_path/'polyData.vtp') vtp = VTK.load(tmp_path/'polyData.vtp')
@ -60,7 +60,7 @@ class TestVTK:
def test_unstructuredGrid(self,tmp_path,cell_type,n): def test_unstructuredGrid(self,tmp_path,cell_type,n):
nodes = np.random.rand(n,3) nodes = np.random.rand(n,3)
connectivity = np.random.choice(np.arange(n),n,False).reshape(-1,n) connectivity = np.random.choice(np.arange(n),n,False).reshape(-1,n)
v = VTK.from_unstructuredGrid(nodes,connectivity,cell_type) v = VTK.from_unstructured_grid(nodes,connectivity,cell_type)
string = v.__repr__() string = v.__repr__()
v.save(tmp_path/'unstructuredGrid',False) v.save(tmp_path/'unstructuredGrid',False)
vtu = VTK.load(tmp_path/'unstructuredGrid.vtu') vtu = VTK.load(tmp_path/'unstructuredGrid.vtu')
@ -72,7 +72,7 @@ class TestVTK:
def test_parallel_out(self,tmp_path): def test_parallel_out(self,tmp_path):
points = np.random.rand(102,3) points = np.random.rand(102,3)
v = VTK.from_polyData(points) v = VTK.from_poly_data(points)
fname_s = tmp_path/'single.vtp' fname_s = tmp_path/'single.vtp'
fname_p = tmp_path/'parallel.vtp' fname_p = tmp_path/'parallel.vtp'
v.save(fname_s,False) v.save(fname_s,False)
@ -121,7 +121,7 @@ class TestVTK:
def test_compare_reference_polyData(self,update,reference_dir,tmp_path): def test_compare_reference_polyData(self,update,reference_dir,tmp_path):
points=np.dstack((np.linspace(0.,1.,10),np.linspace(0.,2.,10),np.linspace(-1.,1.,10))).squeeze() points=np.dstack((np.linspace(0.,1.,10),np.linspace(0.,2.,10),np.linspace(-1.,1.,10))).squeeze()
polyData = VTK.from_polyData(points) polyData = VTK.from_poly_data(points)
polyData.add(points,'coordinates') polyData.add(points,'coordinates')
if update: if update:
polyData.save(reference_dir/'polyData') polyData.save(reference_dir/'polyData')
@ -133,7 +133,7 @@ class TestVTK:
def test_compare_reference_rectilinearGrid(self,update,reference_dir,tmp_path): def test_compare_reference_rectilinearGrid(self,update,reference_dir,tmp_path):
grid = np.array([5,6,7],int) grid = np.array([5,6,7],int)
size = np.array([.6,1.,.5]) size = np.array([.6,1.,.5])
rectilinearGrid = VTK.from_rectilinearGrid(grid,size) rectilinearGrid = VTK.from_rectilinear_grid(grid,size)
c = grid_filters.cell_coord0(grid,size).reshape(-1,3,order='F') c = grid_filters.cell_coord0(grid,size).reshape(-1,3,order='F')
n = grid_filters.node_coord0(grid,size).reshape(-1,3,order='F') n = grid_filters.node_coord0(grid,size).reshape(-1,3,order='F')
rectilinearGrid.add(c,'cell') rectilinearGrid.add(c,'cell')

20
src/CPFEM.f90
View File

@ -107,9 +107,9 @@ subroutine CPFEM_init
print'(/,a)', ' <<<+- CPFEM init -+>>>'; flush(IO_STDOUT) print'(/,a)', ' <<<+- CPFEM init -+>>>'; flush(IO_STDOUT)
allocate(CPFEM_cs( 6,discretization_nIP,discretization_nElem), source= 0.0_pReal) allocate(CPFEM_cs( 6,discretization_nIPs,discretization_Nelems), source= 0.0_pReal)
allocate(CPFEM_dcsdE( 6,6,discretization_nIP,discretization_nElem), source= 0.0_pReal) allocate(CPFEM_dcsdE( 6,6,discretization_nIPs,discretization_Nelems), source= 0.0_pReal)
allocate(CPFEM_dcsdE_knownGood(6,6,discretization_nIP,discretization_nElem), source= 0.0_pReal) allocate(CPFEM_dcsdE_knownGood(6,6,discretization_nIPs,discretization_Nelems), source= 0.0_pReal)
!------------------------------------------------------------------------------ !------------------------------------------------------------------------------
! read debug options ! read debug options
@ -184,8 +184,8 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS
temperature(material_homogenizationAt(elCP))%p(thermalMapping(material_homogenizationAt(elCP))%p(ip,elCP)) = & temperature(material_homogenizationAt(elCP))%p(thermalMapping(material_homogenizationAt(elCP))%p(ip,elCP)) = &
temperature_inp temperature_inp
end select chosenThermal1 end select chosenThermal1
materialpoint_F0(1:3,1:3,ip,elCP) = ffn homogenization_F0(1:3,1:3,ip,elCP) = ffn
materialpoint_F(1:3,1:3,ip,elCP) = ffn1 homogenization_F(1:3,1:3,ip,elCP) = ffn1
if (iand(mode, CPFEM_CALCRESULTS) /= 0_pInt) then if (iand(mode, CPFEM_CALCRESULTS) /= 0_pInt) then
@ -212,17 +212,17 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS
else terminalIllness else terminalIllness
! translate from P to sigma ! translate from P to sigma
Kirchhoff = matmul(materialpoint_P(1:3,1:3,ip,elCP), transpose(materialpoint_F(1:3,1:3,ip,elCP))) Kirchhoff = matmul(homogenization_P(1:3,1:3,ip,elCP), transpose(homogenization_F(1:3,1:3,ip,elCP)))
J_inverse = 1.0_pReal / math_det33(materialpoint_F(1:3,1:3,ip,elCP)) J_inverse = 1.0_pReal / math_det33(homogenization_F(1:3,1:3,ip,elCP))
CPFEM_cs(1:6,ip,elCP) = math_sym33to6(J_inverse * Kirchhoff,weighted=.false.) CPFEM_cs(1:6,ip,elCP) = math_sym33to6(J_inverse * Kirchhoff,weighted=.false.)
! translate from dP/dF to dCS/dE ! translate from dP/dF to dCS/dE
H = 0.0_pReal H = 0.0_pReal
do i=1,3; do j=1,3; do k=1,3; do l=1,3; do m=1,3; do n=1,3 do i=1,3; do j=1,3; do k=1,3; do l=1,3; do m=1,3; do n=1,3
H(i,j,k,l) = H(i,j,k,l) & H(i,j,k,l) = H(i,j,k,l) &
+ materialpoint_F(j,m,ip,elCP) * materialpoint_F(l,n,ip,elCP) & + homogenization_F(j,m,ip,elCP) * homogenization_F(l,n,ip,elCP) &
* materialpoint_dPdF(i,m,k,n,ip,elCP) & * homogenization_dPdF(i,m,k,n,ip,elCP) &
- math_delta(j,l) * materialpoint_F(i,m,ip,elCP) * materialpoint_P(k,m,ip,elCP) & - math_delta(j,l) * homogenization_F(i,m,ip,elCP) * homogenization_P(k,m,ip,elCP) &
+ 0.5_pReal * ( Kirchhoff(j,l)*math_delta(i,k) + Kirchhoff(i,k)*math_delta(j,l) & + 0.5_pReal * ( Kirchhoff(j,l)*math_delta(i,k) + Kirchhoff(i,k)*math_delta(j,l) &
+ Kirchhoff(j,k)*math_delta(i,l) + Kirchhoff(i,l)*math_delta(j,k)) + Kirchhoff(j,k)*math_delta(i,l) + Kirchhoff(i,l)*math_delta(j,k))
enddo; enddo; enddo; enddo; enddo; enddo enddo; enddo; enddo; enddo; enddo; enddo

2
src/IO.f90
View File

@ -491,7 +491,7 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg)
case (705) case (705)
msg = 'Unsupported feature' msg = 'Unsupported feature'
case (706) case (706)
msg = 'Access by incorrect node type' msg = 'Type mismatch in YAML data node'
case (707) case (707)
msg = 'Abrupt end of file' msg = 'Abrupt end of file'
case (708) case (708)

24
src/YAML_types.f90
View File

@ -311,7 +311,7 @@ function tNode_asScalar(self) result(scalar)
class is(tScalar) class is(tScalar)
scalar => self scalar => self
class default class default
call IO_error(706,ext_msg='tNode_asScalar') call IO_error(706,ext_msg='Expected "scalar"')
end select end select
end function tNode_asScalar end function tNode_asScalar
@ -329,7 +329,7 @@ function tNode_asList(self) result(list)
class is(tList) class is(tList)
list => self list => self
class default class default
call IO_error(706,ext_msg='tNode_asList') call IO_error(706,ext_msg='Expected "list"')
end select end select
end function tNode_asList end function tNode_asList
@ -347,7 +347,7 @@ function tNode_asDict(self) result(dict)
class is(tDict) class is(tDict)
dict => self dict => self
class default class default
call IO_error(706,ext_msg='tNode_asDict') call IO_error(706,ext_msg='Expected "dict"')
end select end select
end function tNode_asDict end function tNode_asDict
@ -583,9 +583,9 @@ function tNode_get_byIndex_asStrings(self,i) result(nodeAsStrings)
end function tNode_get_byIndex_asStrings end function tNode_get_byIndex_asStrings
!------------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Returns the key in a dictionary as a string !> @brief Returns the key in a dictionary as a string
!------------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function tNode_getKey_byIndex(self,i) result(key) function tNode_getKey_byIndex(self,i) result(key)
class(tNode), intent(in), target :: self class(tNode), intent(in), target :: self
@ -641,7 +641,7 @@ function tNode_contains(self,k) result(exists)
endif endif
enddo enddo
else else
call IO_error(706,ext_msg='tNode_contains') call IO_error(706,ext_msg='Expected "list" or "dict"')
endif endif
end function tNode_contains end function tNode_contains
@ -1169,17 +1169,21 @@ function tList_asStrings(self)
type(tScalar), pointer :: scalar type(tScalar), pointer :: scalar
len_max = 0 len_max = 0
allocate(character(len=pStringLen) :: tList_asStrings(self%length)) item => self%first
do i = 1, self%length
scalar => item%node%asScalar()
len_max = max(len_max, len_trim(scalar%asString()))
item => item%next
enddo
allocate(character(len=len_max) :: tList_asStrings(self%length))
item => self%first item => self%first
do i = 1, self%length do i = 1, self%length
scalar => item%node%asScalar() scalar => item%node%asScalar()
tList_asStrings(i) = scalar%asString() tList_asStrings(i) = scalar%asString()
len_max = max(len_max, len_trim(tList_asStrings(i)))
item => item%next item => item%next
enddo enddo
!ToDo: trim to len_max
end function tList_asStrings end function tList_asStrings

22
src/constitutive.f90
View File

@ -127,7 +127,7 @@ module constitutive
instance,of,ip,el) instance,of,ip,el)
real(pReal), dimension(3,3), intent(in) :: & real(pReal), dimension(3,3), intent(in) :: &
Mp !< MandelStress Mp !< MandelStress
real(pReal), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem), intent(in) :: & real(pReal), dimension(3,3,homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems), intent(in) :: &
F, & !< deformation gradient F, & !< deformation gradient
Fp !< plastic deformation gradient Fp !< plastic deformation gradient
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
@ -218,7 +218,7 @@ module constitutive
instance, & instance, &
i, & i, &
e e
type(rotation), dimension(1,discretization_nIP,discretization_nElem), intent(in) :: & type(rotation), dimension(1,discretization_nIPs,discretization_Nelems), intent(in) :: &
orientation !< crystal orientation orientation !< crystal orientation
end subroutine plastic_nonlocal_updateCompatibility end subroutine plastic_nonlocal_updateCompatibility
@ -753,7 +753,7 @@ function constitutive_collectDotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el
of of
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
subdt !< timestep subdt !< timestep
real(pReal), intent(in), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: & real(pReal), intent(in), dimension(3,3,homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems) :: &
FArray, & !< elastic deformation gradient FArray, & !< elastic deformation gradient
FpArray !< plastic deformation gradient FpArray !< plastic deformation gradient
real(pReal), intent(in), dimension(3,3) :: & real(pReal), intent(in), dimension(3,3) :: &
@ -905,12 +905,12 @@ end function constitutive_deltaState
!> @brief Allocate the components of the state structure for a given phase !> @brief Allocate the components of the state structure for a given phase
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_allocateState(state, & subroutine constitutive_allocateState(state, &
NipcMyPhase,sizeState,sizeDotState,sizeDeltaState) Nconstituents,sizeState,sizeDotState,sizeDeltaState)
class(tState), intent(out) :: & class(tState), intent(out) :: &
state state
integer, intent(in) :: & integer, intent(in) :: &
NipcMyPhase, & Nconstituents, &
sizeState, & sizeState, &
sizeDotState, & sizeDotState, &
sizeDeltaState sizeDeltaState
@ -921,14 +921,14 @@ subroutine constitutive_allocateState(state, &
state%offsetDeltaState = sizeState-sizeDeltaState ! deltaState occupies latter part of state by definition state%offsetDeltaState = sizeState-sizeDeltaState ! deltaState occupies latter part of state by definition
allocate(state%atol (sizeState), source=0.0_pReal) allocate(state%atol (sizeState), source=0.0_pReal)
allocate(state%state0 (sizeState,NipcMyPhase), source=0.0_pReal) allocate(state%state0 (sizeState,Nconstituents), source=0.0_pReal)
allocate(state%partitionedState0(sizeState,NipcMyPhase), source=0.0_pReal) allocate(state%partitionedState0(sizeState,Nconstituents), source=0.0_pReal)
allocate(state%subState0 (sizeState,NipcMyPhase), source=0.0_pReal) allocate(state%subState0 (sizeState,Nconstituents), source=0.0_pReal)
allocate(state%state (sizeState,NipcMyPhase), source=0.0_pReal) allocate(state%state (sizeState,Nconstituents), source=0.0_pReal)
allocate(state%dotState (sizeDotState,NipcMyPhase), source=0.0_pReal) allocate(state%dotState (sizeDotState,Nconstituents), source=0.0_pReal)
allocate(state%deltaState(sizeDeltaState,NipcMyPhase), source=0.0_pReal) allocate(state%deltaState(sizeDeltaState,Nconstituents), source=0.0_pReal)
end subroutine constitutive_allocateState end subroutine constitutive_allocateState

2
src/constitutive_damage.f90
View File

@ -184,7 +184,7 @@ module subroutine constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi
phiDot = 0.0_pReal phiDot = 0.0_pReal
dPhiDot_dPhi = 0.0_pReal dPhiDot_dPhi = 0.0_pReal
do grain = 1, homogenization_Ngrains(material_homogenizationAt(el)) do grain = 1, homogenization_Nconstituents(material_homogenizationAt(el))
phase = material_phaseAt(grain,el) phase = material_phaseAt(grain,el)
constituent = material_phasememberAt(grain,ip,el) constituent = material_phasememberAt(grain,ip,el)
do source = 1, phase_Nsources(phase) do source = 1, phase_Nsources(phase)

30
src/constitutive_plastic_disloTungsten.f90
View File

@ -78,9 +78,9 @@ module function plastic_disloTungsten_init() result(myPlasticity)
logical, dimension(:), allocatable :: myPlasticity logical, dimension(:), allocatable :: myPlasticity
integer :: & integer :: &
Ninstance, & Ninstances, &
p, i, & p, i, &
NipcMyPhase, & Nconstituents, &
sizeState, sizeDotState, & sizeState, sizeDotState, &
startIndex, endIndex startIndex, endIndex
integer, dimension(:), allocatable :: & integer, dimension(:), allocatable :: &
@ -99,17 +99,17 @@ module function plastic_disloTungsten_init() result(myPlasticity)
print'(/,a)', ' <<<+- plastic_dislotungsten init -+>>>' print'(/,a)', ' <<<+- plastic_dislotungsten init -+>>>'
myPlasticity = plastic_active('disloTungsten') myPlasticity = plastic_active('disloTungsten')
Ninstance = count(myPlasticity) Ninstances = count(myPlasticity)
print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT) print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
if(Ninstance == 0) return if(Ninstances == 0) return
print*, 'Cereceda et al., International Journal of Plasticity 78:242–256, 2016' print*, 'Cereceda et al., International Journal of Plasticity 78:242–256, 2016'
print*, 'https://dx.doi.org/10.1016/j.ijplas.2015.09.002' print*, 'https://dx.doi.org/10.1016/j.ijplas.2015.09.002'
allocate(param(Ninstance)) allocate(param(Ninstances))
allocate(state(Ninstance)) allocate(state(Ninstances))
allocate(dotState(Ninstance)) allocate(dotState(Ninstances))
allocate(dependentState(Ninstance)) allocate(dependentState(Ninstances))
phases => config_material%get('phase') phases => config_material%get('phase')
i = 0 i = 0
@ -221,18 +221,18 @@ module function plastic_disloTungsten_init() result(myPlasticity)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocate state arrays ! allocate state arrays
NipcMyPhase = count(material_phaseAt == p) * discretization_nIP Nconstituents = count(material_phaseAt == p) * discretization_nIPs
sizeDotState = size(['rho_mob ','rho_dip ','gamma_sl']) * prm%sum_N_sl sizeDotState = size(['rho_mob ','rho_dip ','gamma_sl']) * prm%sum_N_sl
sizeState = sizeDotState sizeState = sizeDotState
call constitutive_allocateState(plasticState(p),NipcMyPhase,sizeState,sizeDotState,0) call constitutive_allocateState(plasticState(p),Nconstituents,sizeState,sizeDotState,0)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! state aliases and initialization ! state aliases and initialization
startIndex = 1 startIndex = 1
endIndex = prm%sum_N_sl endIndex = prm%sum_N_sl
stt%rho_mob => plasticState(p)%state(startIndex:endIndex,:) stt%rho_mob => plasticState(p)%state(startIndex:endIndex,:)
stt%rho_mob = spread(rho_mob_0,2,NipcMyPhase) stt%rho_mob = spread(rho_mob_0,2,Nconstituents)
dot%rho_mob => plasticState(p)%dotState(startIndex:endIndex,:) dot%rho_mob => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal) plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal)
if (any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_rho' if (any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_rho'
@ -240,7 +240,7 @@ module function plastic_disloTungsten_init() result(myPlasticity)
startIndex = endIndex + 1 startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_sl endIndex = endIndex + prm%sum_N_sl
stt%rho_dip => plasticState(p)%state(startIndex:endIndex,:) stt%rho_dip => plasticState(p)%state(startIndex:endIndex,:)
stt%rho_dip = spread(rho_dip_0,2,NipcMyPhase) stt%rho_dip = spread(rho_dip_0,2,Nconstituents)
dot%rho_dip => plasticState(p)%dotState(startIndex:endIndex,:) dot%rho_dip => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal) plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal)
@ -252,8 +252,8 @@ module function plastic_disloTungsten_init() result(myPlasticity)
! global alias ! global alias
plasticState(p)%slipRate => plasticState(p)%dotState(startIndex:endIndex,:) plasticState(p)%slipRate => plasticState(p)%dotState(startIndex:endIndex,:)
allocate(dst%Lambda_sl(prm%sum_N_sl,NipcMyPhase), source=0.0_pReal) allocate(dst%Lambda_sl(prm%sum_N_sl,Nconstituents), source=0.0_pReal)
allocate(dst%threshold_stress(prm%sum_N_sl,NipcMyPhase), source=0.0_pReal) allocate(dst%threshold_stress(prm%sum_N_sl,Nconstituents), source=0.0_pReal)
plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally

46
src/constitutive_plastic_dislotwin.f90
View File

@ -126,9 +126,9 @@ module function plastic_dislotwin_init() result(myPlasticity)
logical, dimension(:), allocatable :: myPlasticity logical, dimension(:), allocatable :: myPlasticity
integer :: & integer :: &
Ninstance, & Ninstances, &
p, i, & p, i, &
NipcMyPhase, & Nconstituents, &
sizeState, sizeDotState, & sizeState, sizeDotState, &
startIndex, endIndex startIndex, endIndex
integer, dimension(:), allocatable :: & integer, dimension(:), allocatable :: &
@ -146,9 +146,9 @@ module function plastic_dislotwin_init() result(myPlasticity)
print'(/,a)', ' <<<+- plastic_dislotwin init -+>>>' print'(/,a)', ' <<<+- plastic_dislotwin init -+>>>'
myPlasticity = plastic_active('dislotwin') myPlasticity = plastic_active('dislotwin')
Ninstance = count(myPlasticity) Ninstances = count(myPlasticity)
print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT) print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
if(Ninstance == 0) return if(Ninstances == 0) return
print*, 'Ma and Roters, Acta Materialia 52(12):3603–3612, 2004' print*, 'Ma and Roters, Acta Materialia 52(12):3603–3612, 2004'
print*, 'https://doi.org/10.1016/j.actamat.2004.04.012'//IO_EOL print*, 'https://doi.org/10.1016/j.actamat.2004.04.012'//IO_EOL
@ -159,10 +159,10 @@ module function plastic_dislotwin_init() result(myPlasticity)
print*, 'Wong et al., Acta Materialia 118:140–151, 2016' print*, 'Wong et al., Acta Materialia 118:140–151, 2016'
print*, 'https://doi.org/10.1016/j.actamat.2016.07.032' print*, 'https://doi.org/10.1016/j.actamat.2016.07.032'
allocate(param(Ninstance)) allocate(param(Ninstances))
allocate(state(Ninstance)) allocate(state(Ninstances))
allocate(dotState(Ninstance)) allocate(dotState(Ninstances))
allocate(dependentState(Ninstance)) allocate(dependentState(Ninstances))
phases => config_material%get('phase') phases => config_material%get('phase')
i = 0 i = 0
@ -407,21 +407,21 @@ module function plastic_dislotwin_init() result(myPlasticity)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocate state arrays ! allocate state arrays
NipcMyPhase = count(material_phaseAt == p) * discretization_nIP Nconstituents = count(material_phaseAt == p) * discretization_nIPs
sizeDotState = size(['rho_mob ','rho_dip ','gamma_sl']) * prm%sum_N_sl & sizeDotState = size(['rho_mob ','rho_dip ','gamma_sl']) * prm%sum_N_sl &
+ size(['f_tw']) * prm%sum_N_tw & + size(['f_tw']) * prm%sum_N_tw &
+ size(['f_tr']) * prm%sum_N_tr + size(['f_tr']) * prm%sum_N_tr
sizeState = sizeDotState sizeState = sizeDotState
call constitutive_allocateState(plasticState(p),NipcMyPhase,sizeState,sizeDotState,0) call constitutive_allocateState(plasticState(p),Nconstituents,sizeState,sizeDotState,0)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! locally defined state aliases and initialization of state0 and atol ! locally defined state aliases and initialization of state0 and atol
startIndex = 1 startIndex = 1
endIndex = prm%sum_N_sl endIndex = prm%sum_N_sl
stt%rho_mob=>plasticState(p)%state(startIndex:endIndex,:) stt%rho_mob=>plasticState(p)%state(startIndex:endIndex,:)
stt%rho_mob= spread(rho_mob_0,2,NipcMyPhase) stt%rho_mob= spread(rho_mob_0,2,Nconstituents)
dot%rho_mob=>plasticState(p)%dotState(startIndex:endIndex,:) dot%rho_mob=>plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal) plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal)
if (any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_rho' if (any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_rho'
@ -429,7 +429,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
startIndex = endIndex + 1 startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_sl endIndex = endIndex + prm%sum_N_sl
stt%rho_dip=>plasticState(p)%state(startIndex:endIndex,:) stt%rho_dip=>plasticState(p)%state(startIndex:endIndex,:)
stt%rho_dip= spread(rho_dip_0,2,NipcMyPhase) stt%rho_dip= spread(rho_dip_0,2,Nconstituents)
dot%rho_dip=>plasticState(p)%dotState(startIndex:endIndex,:) dot%rho_dip=>plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal) plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal)
@ -455,18 +455,18 @@ module function plastic_dislotwin_init() result(myPlasticity)
plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('f_trans',defaultVal=1.0e-6_pReal) plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('f_trans',defaultVal=1.0e-6_pReal)
if (any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' f_trans' if (any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' f_trans'
allocate(dst%Lambda_sl (prm%sum_N_sl,NipcMyPhase),source=0.0_pReal) allocate(dst%Lambda_sl (prm%sum_N_sl,Nconstituents),source=0.0_pReal)
allocate(dst%tau_pass (prm%sum_N_sl,NipcMyPhase),source=0.0_pReal) allocate(dst%tau_pass (prm%sum_N_sl,Nconstituents),source=0.0_pReal)
allocate(dst%Lambda_tw (prm%sum_N_tw,NipcMyPhase),source=0.0_pReal) allocate(dst%Lambda_tw (prm%sum_N_tw,Nconstituents),source=0.0_pReal)
allocate(dst%tau_hat_tw (prm%sum_N_tw,NipcMyPhase),source=0.0_pReal) allocate(dst%tau_hat_tw (prm%sum_N_tw,Nconstituents),source=0.0_pReal)
allocate(dst%tau_r_tw (prm%sum_N_tw,NipcMyPhase),source=0.0_pReal) allocate(dst%tau_r_tw (prm%sum_N_tw,Nconstituents),source=0.0_pReal)
allocate(dst%V_tw (prm%sum_N_tw,NipcMyPhase),source=0.0_pReal) allocate(dst%V_tw (prm%sum_N_tw,Nconstituents),source=0.0_pReal)
allocate(dst%Lambda_tr (prm%sum_N_tr,NipcMyPhase),source=0.0_pReal) allocate(dst%Lambda_tr (prm%sum_N_tr,Nconstituents),source=0.0_pReal)
allocate(dst%tau_hat_tr (prm%sum_N_tr,NipcMyPhase),source=0.0_pReal) allocate(dst%tau_hat_tr (prm%sum_N_tr,Nconstituents),source=0.0_pReal)
allocate(dst%tau_r_tr (prm%sum_N_tr,NipcMyPhase),source=0.0_pReal) allocate(dst%tau_r_tr (prm%sum_N_tr,Nconstituents),source=0.0_pReal)
allocate(dst%V_tr (prm%sum_N_tr,NipcMyPhase),source=0.0_pReal) allocate(dst%V_tr (prm%sum_N_tr,Nconstituents),source=0.0_pReal)
plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally

20
src/constitutive_plastic_isotropic.f90
View File

@ -53,10 +53,10 @@ module function plastic_isotropic_init() result(myPlasticity)
logical, dimension(:), allocatable :: myPlasticity logical, dimension(:), allocatable :: myPlasticity
integer :: & integer :: &
Ninstance, & Ninstances, &
p, & p, &
i, & i, &
NipcMyPhase, & Nconstituents, &
sizeState, sizeDotState sizeState, sizeDotState
real(pReal) :: & real(pReal) :: &
xi_0 !< initial critical stress xi_0 !< initial critical stress
@ -70,16 +70,16 @@ module function plastic_isotropic_init() result(myPlasticity)
print'(/,a)', ' <<<+- plastic_isotropic init -+>>>' print'(/,a)', ' <<<+- plastic_isotropic init -+>>>'
myPlasticity = plastic_active('isotropic') myPlasticity = plastic_active('isotropic')
Ninstance = count(myPlasticity) Ninstances = count(myPlasticity)
print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT) print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
if(Ninstance == 0) return if(Ninstances == 0) return
print*, 'Maiti and Eisenlohr, Scripta Materialia 145:37–40, 2018' print*, 'Maiti and Eisenlohr, Scripta Materialia 145:37–40, 2018'
print*, 'https://doi.org/10.1016/j.scriptamat.2017.09.047' print*, 'https://doi.org/10.1016/j.scriptamat.2017.09.047'
allocate(param(Ninstance)) allocate(param(Ninstances))
allocate(state(Ninstance)) allocate(state(Ninstances))
allocate(dotState(Ninstance)) allocate(dotState(Ninstances))
phases => config_material%get('phase') phases => config_material%get('phase')
i = 0 i = 0
@ -130,11 +130,11 @@ module function plastic_isotropic_init() result(myPlasticity)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocate state arrays ! allocate state arrays
NipcMyPhase = count(material_phaseAt == p) * discretization_nIP Nconstituents = count(material_phaseAt == p) * discretization_nIPs
sizeDotState = size(['xi ','gamma']) sizeDotState = size(['xi ','gamma'])
sizeState = sizeDotState sizeState = sizeDotState
call constitutive_allocateState(plasticState(p),NipcMyPhase,sizeState,sizeDotState,0) call constitutive_allocateState(plasticState(p),Nconstituents,sizeState,sizeDotState,0)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! state aliases and initialization ! state aliases and initialization

24
src/constitutive_plastic_kinehardening.f90
View File

@ -62,9 +62,9 @@ module function plastic_kinehardening_init() result(myPlasticity)
logical, dimension(:), allocatable :: myPlasticity logical, dimension(:), allocatable :: myPlasticity
integer :: & integer :: &
Ninstance, & Ninstances, &
p, i, o, & p, i, o, &
NipcMyPhase, & Nconstituents, &
sizeState, sizeDeltaState, sizeDotState, & sizeState, sizeDeltaState, sizeDotState, &
startIndex, endIndex startIndex, endIndex
integer, dimension(:), allocatable :: & integer, dimension(:), allocatable :: &
@ -82,14 +82,14 @@ module function plastic_kinehardening_init() result(myPlasticity)
print'(/,a)', ' <<<+- plastic_kinehardening init -+>>>' print'(/,a)', ' <<<+- plastic_kinehardening init -+>>>'
myPlasticity = plastic_active('kinehardening') myPlasticity = plastic_active('kinehardening')
Ninstance = count(myPlasticity) Ninstances = count(myPlasticity)
print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT) print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
if(Ninstance == 0) return if(Ninstances == 0) return
allocate(param(Ninstance)) allocate(param(Ninstances))
allocate(state(Ninstance)) allocate(state(Ninstances))
allocate(dotState(Ninstance)) allocate(dotState(Ninstances))
allocate(deltaState(Ninstance)) allocate(deltaState(Ninstances))
phases => config_material%get('phase') phases => config_material%get('phase')
i = 0 i = 0
@ -174,19 +174,19 @@ module function plastic_kinehardening_init() result(myPlasticity)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocate state arrays ! allocate state arrays
NipcMyPhase = count(material_phaseAt == p) * discretization_nIP Nconstituents = count(material_phaseAt == p) * discretization_nIPs
sizeDotState = size(['crss ','crss_back', 'accshear ']) * prm%sum_N_sl!ToDo: adjust names, ask Philip sizeDotState = size(['crss ','crss_back', 'accshear ']) * prm%sum_N_sl!ToDo: adjust names, ask Philip
sizeDeltaState = size(['sense ', 'chi0 ', 'gamma0' ]) * prm%sum_N_sl !ToDo: adjust names sizeDeltaState = size(['sense ', 'chi0 ', 'gamma0' ]) * prm%sum_N_sl !ToDo: adjust names
sizeState = sizeDotState + sizeDeltaState sizeState = sizeDotState + sizeDeltaState
call constitutive_allocateState(plasticState(p),NipcMyPhase,sizeState,sizeDotState,sizeDeltaState) call constitutive_allocateState(plasticState(p),Nconstituents,sizeState,sizeDotState,sizeDeltaState)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! state aliases and initialization ! state aliases and initialization
startIndex = 1 startIndex = 1
endIndex = prm%sum_N_sl endIndex = prm%sum_N_sl
stt%crss => plasticState(p)%state (startIndex:endIndex,:) stt%crss => plasticState(p)%state (startIndex:endIndex,:)
stt%crss = spread(xi_0, 2, NipcMyPhase) stt%crss = spread(xi_0, 2, Nconstituents)
dot%crss => plasticState(p)%dotState(startIndex:endIndex,:) dot%crss => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal) plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal)
if(any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi' if(any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi'

14
src/constitutive_plastic_none.f90
View File

@ -16,9 +16,9 @@ module function plastic_none_init() result(myPlasticity)
logical, dimension(:), allocatable :: myPlasticity logical, dimension(:), allocatable :: myPlasticity
integer :: & integer :: &
Ninstance, & Ninstances, &
p, & p, &
NipcMyPhase Nconstituents
class(tNode), pointer :: & class(tNode), pointer :: &
phases, & phases, &
phase, & phase, &
@ -34,15 +34,15 @@ module function plastic_none_init() result(myPlasticity)
if(pl%get_asString('type') == 'none') myPlasticity(p) = .true. if(pl%get_asString('type') == 'none') myPlasticity(p) = .true.
enddo enddo
Ninstance = count(myPlasticity) Ninstances = count(myPlasticity)
print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT) print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
if(Ninstance == 0) return if(Ninstances == 0) return
do p = 1, phases%length do p = 1, phases%length
phase => phases%get(p) phase => phases%get(p)
if(.not. myPlasticity(p)) cycle if(.not. myPlasticity(p)) cycle
NipcMyPhase = count(material_phaseAt == p) * discretization_nIP Nconstituents = count(material_phaseAt == p) * discretization_nIPs
call constitutive_allocateState(plasticState(p),NipcMyPhase,0,0,0) call constitutive_allocateState(plasticState(p),Nconstituents,0,0,0)
enddo enddo
end function plastic_none_init end function plastic_none_init

84
src/constitutive_plastic_nonlocal.f90
View File

@ -153,7 +153,7 @@ submodule(constitutive:constitutive_plastic) plastic_nonlocal
state, & state, &
state0 state0
type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance) type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstances)
type(tNonlocalMicrostructure), dimension(:), allocatable :: microstructure type(tNonlocalMicrostructure), dimension(:), allocatable :: microstructure
@ -168,9 +168,9 @@ module function plastic_nonlocal_init() result(myPlasticity)
logical, dimension(:), allocatable :: myPlasticity logical, dimension(:), allocatable :: myPlasticity
integer :: & integer :: &
Ninstance, & Ninstances, &
p, i, & p, i, &
NipcMyPhase, & Nconstituents, &
sizeState, sizeDotState, sizeDependentState, sizeDeltaState, & sizeState, sizeDotState, sizeDependentState, sizeDeltaState, &
s1, s2, & s1, s2, &
s, t, l s, t, l
@ -188,9 +188,9 @@ module function plastic_nonlocal_init() result(myPlasticity)
print'(/,a)', ' <<<+- plastic_nonlocal init -+>>>' print'(/,a)', ' <<<+- plastic_nonlocal init -+>>>'
myPlasticity = plastic_active('nonlocal') myPlasticity = plastic_active('nonlocal')
Ninstance = count(myPlasticity) Ninstances = count(myPlasticity)
print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT) print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
if(Ninstance == 0) then if(Ninstances == 0) then
call geometry_plastic_nonlocal_disable call geometry_plastic_nonlocal_disable
return return
endif endif
@ -201,12 +201,12 @@ module function plastic_nonlocal_init() result(myPlasticity)
print*, 'Kords, Dissertation RWTH Aachen, 2014' print*, 'Kords, Dissertation RWTH Aachen, 2014'
print*, 'http://publications.rwth-aachen.de/record/229993' print*, 'http://publications.rwth-aachen.de/record/229993'
allocate(param(Ninstance)) allocate(param(Ninstances))
allocate(state(Ninstance)) allocate(state(Ninstances))
allocate(state0(Ninstance)) allocate(state0(Ninstances))
allocate(dotState(Ninstance)) allocate(dotState(Ninstances))
allocate(deltaState(Ninstance)) allocate(deltaState(Ninstances))
allocate(microstructure(Ninstance)) allocate(microstructure(Ninstances))
phases => config_material%get('phase') phases => config_material%get('phase')
i = 0 i = 0
@ -391,7 +391,7 @@ module function plastic_nonlocal_init() result(myPlasticity)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocate state arrays ! allocate state arrays
NipcMyPhase = count(material_phaseAt==p) * discretization_nIP Nconstituents = count(material_phaseAt==p) * discretization_nIPs
sizeDotState = size([ 'rhoSglEdgePosMobile ','rhoSglEdgeNegMobile ', & sizeDotState = size([ 'rhoSglEdgePosMobile ','rhoSglEdgeNegMobile ', &
'rhoSglScrewPosMobile ','rhoSglScrewNegMobile ', & 'rhoSglScrewPosMobile ','rhoSglScrewNegMobile ', &
'rhoSglEdgePosImmobile ','rhoSglEdgeNegImmobile ', & 'rhoSglEdgePosImmobile ','rhoSglEdgeNegImmobile ', &
@ -405,7 +405,7 @@ module function plastic_nonlocal_init() result(myPlasticity)
'maxDipoleHeightEdge ','maxDipoleHeightScrew' ]) * prm%sum_N_sl !< other dependent state variables that are not updated by microstructure 'maxDipoleHeightEdge ','maxDipoleHeightScrew' ]) * prm%sum_N_sl !< other dependent state variables that are not updated by microstructure
sizeDeltaState = sizeDotState sizeDeltaState = sizeDotState
call constitutive_allocateState(plasticState(p),NipcMyPhase,sizeState,sizeDotState,sizeDeltaState) call constitutive_allocateState(plasticState(p),Nconstituents,sizeState,sizeDotState,sizeDeltaState)
plasticState(p)%nonlocal = pl%get_asBool('nonlocal') plasticState(p)%nonlocal = pl%get_asBool('nonlocal')
if(plasticState(p)%nonlocal .and. .not. allocated(IPneighborhood)) & if(plasticState(p)%nonlocal .and. .not. allocated(IPneighborhood)) &
@ -476,26 +476,26 @@ module function plastic_nonlocal_init() result(myPlasticity)
dot%rho_dip_scr => plasticState(p)%dotState (9*prm%sum_N_sl+1:10*prm%sum_N_sl,:) dot%rho_dip_scr => plasticState(p)%dotState (9*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
del%rho_dip_scr => plasticState(p)%deltaState (9*prm%sum_N_sl+1:10*prm%sum_N_sl,:) del%rho_dip_scr => plasticState(p)%deltaState (9*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
stt%gamma => plasticState(p)%state (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:NipcMyPhase) stt%gamma => plasticState(p)%state (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nconstituents)
dot%gamma => plasticState(p)%dotState (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:NipcMyPhase) dot%gamma => plasticState(p)%dotState (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nconstituents)
del%gamma => plasticState(p)%deltaState (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:NipcMyPhase) del%gamma => plasticState(p)%deltaState (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nconstituents)
plasticState(p)%atol(10*prm%sum_N_sl+1:11*prm%sum_N_sl ) = pl%get_asFloat('atol_gamma', defaultVal = 1.0e-2_pReal) plasticState(p)%atol(10*prm%sum_N_sl+1:11*prm%sum_N_sl ) = pl%get_asFloat('atol_gamma', defaultVal = 1.0e-2_pReal)
if(any(plasticState(p)%atol(10*prm%sum_N_sl+1:11*prm%sum_N_sl) < 0.0_pReal)) & if(any(plasticState(p)%atol(10*prm%sum_N_sl+1:11*prm%sum_N_sl) < 0.0_pReal)) &
extmsg = trim(extmsg)//' atol_gamma' extmsg = trim(extmsg)//' atol_gamma'
plasticState(p)%slipRate => plasticState(p)%dotState (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:NipcMyPhase) plasticState(p)%slipRate => plasticState(p)%dotState (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nconstituents)
stt%rho_forest => plasticState(p)%state (11*prm%sum_N_sl + 1:12*prm%sum_N_sl,1:NipcMyPhase) stt%rho_forest => plasticState(p)%state (11*prm%sum_N_sl + 1:12*prm%sum_N_sl,1:Nconstituents)
stt%v => plasticState(p)%state (12*prm%sum_N_sl + 1:16*prm%sum_N_sl,1:NipcMyPhase) stt%v => plasticState(p)%state (12*prm%sum_N_sl + 1:16*prm%sum_N_sl,1:Nconstituents)
stt%v_edg_pos => plasticState(p)%state (12*prm%sum_N_sl + 1:13*prm%sum_N_sl,1:NipcMyPhase) stt%v_edg_pos => plasticState(p)%state (12*prm%sum_N_sl + 1:13*prm%sum_N_sl,1:Nconstituents)
stt%v_edg_neg => plasticState(p)%state (13*prm%sum_N_sl + 1:14*prm%sum_N_sl,1:NipcMyPhase) stt%v_edg_neg => plasticState(p)%state (13*prm%sum_N_sl + 1:14*prm%sum_N_sl,1:Nconstituents)
stt%v_scr_pos => plasticState(p)%state (14*prm%sum_N_sl + 1:15*prm%sum_N_sl,1:NipcMyPhase) stt%v_scr_pos => plasticState(p)%state (14*prm%sum_N_sl + 1:15*prm%sum_N_sl,1:Nconstituents)
stt%v_scr_neg => plasticState(p)%state (15*prm%sum_N_sl + 1:16*prm%sum_N_sl,1:NipcMyPhase) stt%v_scr_neg => plasticState(p)%state (15*prm%sum_N_sl + 1:16*prm%sum_N_sl,1:Nconstituents)
allocate(dst%tau_pass(prm%sum_N_sl,NipcMyPhase),source=0.0_pReal) allocate(dst%tau_pass(prm%sum_N_sl,Nconstituents),source=0.0_pReal)
allocate(dst%tau_back(prm%sum_N_sl,NipcMyPhase),source=0.0_pReal) allocate(dst%tau_back(prm%sum_N_sl,Nconstituents),source=0.0_pReal)
end associate end associate
if (NipcMyPhase > 0) call stateInit(ini,p,NipcMyPhase,i) if (Nconstituents > 0) call stateInit(ini,p,Nconstituents,i)
plasticState(p)%state0 = plasticState(p)%state plasticState(p)%state0 = plasticState(p)%state
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -505,12 +505,12 @@ module function plastic_nonlocal_init() result(myPlasticity)
enddo enddo
allocate(compatibility(2,maxval(param%sum_N_sl),maxval(param%sum_N_sl),nIPneighbors,& allocate(compatibility(2,maxval(param%sum_N_sl),maxval(param%sum_N_sl),nIPneighbors,&
discretization_nIP,discretization_nElem), source=0.0_pReal) discretization_nIPs,discretization_Nelems), source=0.0_pReal)
! BEGIN DEPRECATED---------------------------------------------------------------------------------- ! BEGIN DEPRECATED----------------------------------------------------------------------------------
allocate(iRhoU(maxval(param%sum_N_sl),4,Ninstance), source=0) allocate(iRhoU(maxval(param%sum_N_sl),4,Ninstances), source=0)
allocate(iV(maxval(param%sum_N_sl),4,Ninstance), source=0) allocate(iV(maxval(param%sum_N_sl),4,Ninstances), source=0)
allocate(iD(maxval(param%sum_N_sl),2,Ninstance), source=0) allocate(iD(maxval(param%sum_N_sl),2,Ninstances), source=0)
i = 0 i = 0
do p = 1, phases%length do p = 1, phases%length
@ -519,7 +519,7 @@ module function plastic_nonlocal_init() result(myPlasticity)
if(.not. myPlasticity(p)) cycle if(.not. myPlasticity(p)) cycle
i = i + 1 i = i + 1
NipcMyPhase = count(material_phaseAt==p) * discretization_nIP Nconstituents = count(material_phaseAt==p) * discretization_nIPs
l = 0 l = 0
do t = 1,4 do t = 1,4
do s = 1,param(i)%sum_N_sl do s = 1,param(i)%sum_N_sl
@ -976,7 +976,7 @@ module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature,timestep, &
real(pReal), dimension(3,3), intent(in) :: & real(pReal), dimension(3,3), intent(in) :: &
Mp !< MandelStress Mp !< MandelStress
real(pReal), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem), intent(in) :: & real(pReal), dimension(3,3,homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems), intent(in) :: &
F, & !< elastic deformation gradient F, & !< elastic deformation gradient
Fp !< plastic deformation gradient Fp !< plastic deformation gradient
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
@ -1176,7 +1176,7 @@ end subroutine plastic_nonlocal_dotState
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
function rhoDotFlux(F,Fp,timestep, instance,of,ip,el) function rhoDotFlux(F,Fp,timestep, instance,of,ip,el)
real(pReal), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem), intent(in) :: & real(pReal), dimension(3,3,homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems), intent(in) :: &
F, & !< elastic deformation gradient F, & !< elastic deformation gradient
Fp !< plastic deformation gradient Fp !< plastic deformation gradient
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
@ -1416,7 +1416,7 @@ end function rhoDotFlux
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine plastic_nonlocal_updateCompatibility(orientation,instance,i,e) module subroutine plastic_nonlocal_updateCompatibility(orientation,instance,i,e)
type(rotation), dimension(1,discretization_nIP,discretization_nElem), intent(in) :: & type(rotation), dimension(1,discretization_nIPs,discretization_Nelems), intent(in) :: &
orientation ! crystal orientation orientation ! crystal orientation
integer, intent(in) :: & integer, intent(in) :: &
instance, & instance, &
@ -1601,13 +1601,13 @@ end subroutine plastic_nonlocal_results
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief populates the initial dislocation density !> @brief populates the initial dislocation density
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine stateInit(ini,phase,NipcMyPhase,instance) subroutine stateInit(ini,phase,Nconstituents,instance)
type(tInitialParameters) :: & type(tInitialParameters) :: &
ini ini
integer,intent(in) :: & integer,intent(in) :: &
phase, & phase, &
NipcMyPhase, & Nconstituents, &
instance instance
integer :: & integer :: &
e, & e, &
@ -1625,15 +1625,15 @@ subroutine stateInit(ini,phase,NipcMyPhase,instance)
totalVolume, & totalVolume, &
densityBinning, & densityBinning, &
minimumIpVolume minimumIpVolume
real(pReal), dimension(NipcMyPhase) :: & real(pReal), dimension(Nconstituents) :: &
volume volume
associate(stt => state(instance)) associate(stt => state(instance))
if (ini%random_rho_u > 0.0_pReal) then ! randomly distribute dislocation segments on random slip system and of random type in the volume if (ini%random_rho_u > 0.0_pReal) then ! randomly distribute dislocation segments on random slip system and of random type in the volume
do e = 1,discretization_nElem do e = 1,discretization_Nelems
do i = 1,discretization_nIP do i = 1,discretization_nIPs
if (material_phaseAt(1,e) == phase) volume(material_phasememberAt(1,i,e)) = IPvolume(i,e) if (material_phaseAt(1,e) == phase) volume(material_phasememberAt(1,i,e)) = IPvolume(i,e)
enddo enddo
enddo enddo
@ -1645,13 +1645,13 @@ subroutine stateInit(ini,phase,NipcMyPhase,instance)
meanDensity = 0.0_pReal meanDensity = 0.0_pReal
do while(meanDensity < ini%random_rho_u) do while(meanDensity < ini%random_rho_u)
call random_number(rnd) call random_number(rnd)
phasemember = nint(rnd(1)*real(NipcMyPhase,pReal) + 0.5_pReal) phasemember = nint(rnd(1)*real(Nconstituents,pReal) + 0.5_pReal)
s = nint(rnd(2)*real(sum(ini%N_sl),pReal)*4.0_pReal + 0.5_pReal) s = nint(rnd(2)*real(sum(ini%N_sl),pReal)*4.0_pReal + 0.5_pReal)
meanDensity = meanDensity + densityBinning * volume(phasemember) / totalVolume meanDensity = meanDensity + densityBinning * volume(phasemember) / totalVolume
stt%rhoSglMobile(s,phasemember) = densityBinning stt%rhoSglMobile(s,phasemember) = densityBinning
enddo enddo
else ! homogeneous distribution with noise else ! homogeneous distribution with noise
do e = 1, NipcMyPhase do e = 1, Nconstituents
do f = 1,size(ini%N_sl,1) do f = 1,size(ini%N_sl,1)
from = 1 + sum(ini%N_sl(1:f-1)) from = 1 + sum(ini%N_sl(1:f-1))
upto = sum(ini%N_sl(1:f)) upto = sum(ini%N_sl(1:f))

24
src/constitutive_plastic_phenopowerlaw.f90
View File

@ -70,9 +70,9 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
logical, dimension(:), allocatable :: myPlasticity logical, dimension(:), allocatable :: myPlasticity
integer :: & integer :: &
Ninstance, & Ninstances, &
p, i, & p, i, &
NipcMyPhase, & Nconstituents, &
sizeState, sizeDotState, & sizeState, sizeDotState, &
startIndex, endIndex startIndex, endIndex
integer, dimension(:), allocatable :: & integer, dimension(:), allocatable :: &
@ -91,13 +91,13 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
print'(/,a)', ' <<<+- plastic_phenopowerlaw init -+>>>' print'(/,a)', ' <<<+- plastic_phenopowerlaw init -+>>>'
myPlasticity = plastic_active('phenopowerlaw') myPlasticity = plastic_active('phenopowerlaw')
Ninstance = count(myPlasticity) Ninstances = count(myPlasticity)
print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT) print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
if(Ninstance == 0) return if(Ninstances == 0) return
allocate(param(Ninstance)) allocate(param(Ninstances))
allocate(state(Ninstance)) allocate(state(Ninstances))
allocate(dotState(Ninstance)) allocate(dotState(Ninstances))
phases => config_material%get('phase') phases => config_material%get('phase')
i = 0 i = 0
@ -224,20 +224,20 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocate state arrays ! allocate state arrays
NipcMyPhase = count(material_phaseAt == p) * discretization_nIP Nconstituents = count(material_phaseAt == p) * discretization_nIPs
sizeDotState = size(['xi_sl ','gamma_sl']) * prm%sum_N_sl & sizeDotState = size(['xi_sl ','gamma_sl']) * prm%sum_N_sl &
+ size(['xi_tw ','gamma_tw']) * prm%sum_N_tw + size(['xi_tw ','gamma_tw']) * prm%sum_N_tw
sizeState = sizeDotState sizeState = sizeDotState
call constitutive_allocateState(plasticState(p),NipcMyPhase,sizeState,sizeDotState,0) call constitutive_allocateState(plasticState(p),Nconstituents,sizeState,sizeDotState,0)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! state aliases and initialization ! state aliases and initialization
startIndex = 1 startIndex = 1
endIndex = prm%sum_N_sl endIndex = prm%sum_N_sl
stt%xi_slip => plasticState(p)%state (startIndex:endIndex,:) stt%xi_slip => plasticState(p)%state (startIndex:endIndex,:)
stt%xi_slip = spread(xi_0_sl, 2, NipcMyPhase) stt%xi_slip = spread(xi_0_sl, 2, Nconstituents)
dot%xi_slip => plasticState(p)%dotState(startIndex:endIndex,:) dot%xi_slip => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal) plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal)
if(any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi' if(any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi'
@ -245,7 +245,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
startIndex = endIndex + 1 startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_tw endIndex = endIndex + prm%sum_N_tw
stt%xi_twin => plasticState(p)%state (startIndex:endIndex,:) stt%xi_twin => plasticState(p)%state (startIndex:endIndex,:)
stt%xi_twin = spread(xi_0_tw, 2, NipcMyPhase) stt%xi_twin = spread(xi_0_tw, 2, Nconstituents)
dot%xi_twin => plasticState(p)%dotState(startIndex:endIndex,:) dot%xi_twin => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal) plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal)
if(any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi' if(any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi'

2
src/constitutive_thermal.f90
View File

@ -95,7 +95,7 @@ module subroutine constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T,
homog = material_homogenizationAt(el) homog = material_homogenizationAt(el)
instance = thermal_typeInstance(homog) instance = thermal_typeInstance(homog)
do grain = 1, homogenization_Ngrains(homog) do grain = 1, homogenization_Nconstituents(homog)
phase = material_phaseAt(grain,el) phase = material_phaseAt(grain,el)
constituent = material_phasememberAt(grain,ip,el) constituent = material_phasememberAt(grain,ip,el)
do source = 1, phase_Nsources(phase) do source = 1, phase_Nsources(phase)

42
src/crystallite.f90
View File

@ -135,7 +135,7 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine crystallite_init subroutine crystallite_init
logical, dimension(discretization_nIP,discretization_nElem) :: devNull logical, dimension(discretization_nIPs,discretization_Nelems) :: devNull
integer :: & integer :: &
c, & !< counter in integration point component loop c, & !< counter in integration point component loop
i, & !< counter in integration point loop i, & !< counter in integration point loop
@ -162,9 +162,9 @@ subroutine crystallite_init
debugCrystallite%ip = config_debug%get_asInt('integrationpoint', defaultVal=1) debugCrystallite%ip = config_debug%get_asInt('integrationpoint', defaultVal=1)
debugCrystallite%grain = config_debug%get_asInt('grain', defaultVal=1) debugCrystallite%grain = config_debug%get_asInt('grain', defaultVal=1)
cMax = homogenization_maxNgrains cMax = homogenization_maxNconstituents
iMax = discretization_nIP iMax = discretization_nIPs
eMax = discretization_nElem eMax = discretization_Nelems
allocate(crystallite_partitionedF(3,3,cMax,iMax,eMax),source=0.0_pReal) allocate(crystallite_partitionedF(3,3,cMax,iMax,eMax),source=0.0_pReal)
@ -253,7 +253,7 @@ subroutine crystallite_init
! initialize ! initialize
!$OMP PARALLEL DO PRIVATE(i,c) !$OMP PARALLEL DO PRIVATE(i,c)
do e = FEsolving_execElem(1),FEsolving_execElem(2) do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1), FEsolving_execIP(2); do c = 1, homogenization_Ngrains(material_homogenizationAt(e)) do i = FEsolving_execIP(1), FEsolving_execIP(2); do c = 1, homogenization_Nconstituents(material_homogenizationAt(e))
crystallite_Fp0(1:3,1:3,c,i,e) = material_orientation0(c,i,e)%asMatrix() ! Fp reflects initial orientation (see 10.1016/j.actamat.2006.01.005) crystallite_Fp0(1:3,1:3,c,i,e) = material_orientation0(c,i,e)%asMatrix() ! Fp reflects initial orientation (see 10.1016/j.actamat.2006.01.005)
crystallite_Fp0(1:3,1:3,c,i,e) = crystallite_Fp0(1:3,1:3,c,i,e) & crystallite_Fp0(1:3,1:3,c,i,e) = crystallite_Fp0(1:3,1:3,c,i,e) &
/ math_det33(crystallite_Fp0(1:3,1:3,c,i,e))**(1.0_pReal/3.0_pReal) / math_det33(crystallite_Fp0(1:3,1:3,c,i,e))**(1.0_pReal/3.0_pReal)
@ -279,7 +279,7 @@ subroutine crystallite_init
!$OMP PARALLEL DO !$OMP PARALLEL DO
do e = FEsolving_execElem(1),FEsolving_execElem(2) do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1),FEsolving_execIP(2) do i = FEsolving_execIP(1),FEsolving_execIP(2)
do c = 1,homogenization_Ngrains(material_homogenizationAt(e)) do c = 1,homogenization_Nconstituents(material_homogenizationAt(e))
call constitutive_dependentState(crystallite_partitionedF0(1:3,1:3,c,i,e), & call constitutive_dependentState(crystallite_partitionedF0(1:3,1:3,c,i,e), &
crystallite_partitionedFp0(1:3,1:3,c,i,e), & crystallite_partitionedFp0(1:3,1:3,c,i,e), &
c,i,e) ! update dependent state variables to be consistent with basic states c,i,e) ! update dependent state variables to be consistent with basic states
@ -307,7 +307,7 @@ end subroutine crystallite_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function crystallite_stress() function crystallite_stress()
logical, dimension(discretization_nIP,discretization_nElem) :: crystallite_stress logical, dimension(discretization_nIPs,discretization_Nelems) :: crystallite_stress
real(pReal) :: & real(pReal) :: &
formerSubStep formerSubStep
integer :: & integer :: &
@ -316,7 +316,7 @@ function crystallite_stress()
i, & !< counter in integration point loop i, & !< counter in integration point loop
e, & !< counter in element loop e, & !< counter in element loop
s s
logical, dimension(homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: todo !ToDo: need to set some values to false for different Ngrains logical, dimension(homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems) :: todo !ToDo: need to set some values to false for different Ngrains
real(pReal), dimension(:,:,:,:,:), allocatable :: & real(pReal), dimension(:,:,:,:,:), allocatable :: &
subLp0,& !< plastic velocity grad at start of crystallite inc subLp0,& !< plastic velocity grad at start of crystallite inc
subLi0 !< intermediate velocity grad at start of crystallite inc subLi0 !< intermediate velocity grad at start of crystallite inc
@ -334,7 +334,7 @@ function crystallite_stress()
crystallite_subStep = 0.0_pReal crystallite_subStep = 0.0_pReal
!$OMP PARALLEL DO !$OMP PARALLEL DO
elementLooping1: do e = FEsolving_execElem(1),FEsolving_execElem(2) elementLooping1: do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1),FEsolving_execIP(2); do c = 1,homogenization_Ngrains(material_homogenizationAt(e)) do i = FEsolving_execIP(1),FEsolving_execIP(2); do c = 1,homogenization_Nconstituents(material_homogenizationAt(e))
homogenizationRequestsCalculation: if (crystallite_requested(c,i,e)) then homogenizationRequestsCalculation: if (crystallite_requested(c,i,e)) then
plasticState (material_phaseAt(c,e))%subState0( :,material_phaseMemberAt(c,i,e)) = & plasticState (material_phaseAt(c,e))%subState0( :,material_phaseMemberAt(c,i,e)) = &
plasticState (material_phaseAt(c,e))%partitionedState0(:,material_phaseMemberAt(c,i,e)) plasticState (material_phaseAt(c,e))%partitionedState0(:,material_phaseMemberAt(c,i,e))
@ -366,7 +366,7 @@ function crystallite_stress()
!$OMP PARALLEL DO PRIVATE(formerSubStep) !$OMP PARALLEL DO PRIVATE(formerSubStep)
elementLooping3: do e = FEsolving_execElem(1),FEsolving_execElem(2) elementLooping3: do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1),FEsolving_execIP(2) do i = FEsolving_execIP(1),FEsolving_execIP(2)
do c = 1,homogenization_Ngrains(material_homogenizationAt(e)) do c = 1,homogenization_Nconstituents(material_homogenizationAt(e))
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! wind forward ! wind forward
if (crystallite_converged(c,i,e)) then if (crystallite_converged(c,i,e)) then
@ -462,7 +462,7 @@ subroutine crystallite_initializeRestorationPoints(i,e)
c, & !< constituent number c, & !< constituent number
s s
do c = 1,homogenization_Ngrains(material_homogenizationAt(e)) do c = 1,homogenization_Nconstituents(material_homogenizationAt(e))
crystallite_partitionedFp0(1:3,1:3,c,i,e) = crystallite_Fp0(1:3,1:3,c,i,e) crystallite_partitionedFp0(1:3,1:3,c,i,e) = crystallite_Fp0(1:3,1:3,c,i,e)
crystallite_partitionedLp0(1:3,1:3,c,i,e) = crystallite_Lp0(1:3,1:3,c,i,e) crystallite_partitionedLp0(1:3,1:3,c,i,e) = crystallite_Lp0(1:3,1:3,c,i,e)
crystallite_partitionedFi0(1:3,1:3,c,i,e) = crystallite_Fi0(1:3,1:3,c,i,e) crystallite_partitionedFi0(1:3,1:3,c,i,e) = crystallite_Fi0(1:3,1:3,c,i,e)
@ -493,7 +493,7 @@ subroutine crystallite_windForward(i,e)
c, & !< constituent number c, & !< constituent number
s s
do c = 1,homogenization_Ngrains(material_homogenizationAt(e)) do c = 1,homogenization_Nconstituents(material_homogenizationAt(e))
crystallite_partitionedF0 (1:3,1:3,c,i,e) = crystallite_partitionedF(1:3,1:3,c,i,e) crystallite_partitionedF0 (1:3,1:3,c,i,e) = crystallite_partitionedF(1:3,1:3,c,i,e)
crystallite_partitionedFp0(1:3,1:3,c,i,e) = crystallite_Fp (1:3,1:3,c,i,e) crystallite_partitionedFp0(1:3,1:3,c,i,e) = crystallite_Fp (1:3,1:3,c,i,e)
crystallite_partitionedLp0(1:3,1:3,c,i,e) = crystallite_Lp (1:3,1:3,c,i,e) crystallite_partitionedLp0(1:3,1:3,c,i,e) = crystallite_Lp (1:3,1:3,c,i,e)
@ -526,7 +526,7 @@ subroutine crystallite_restore(i,e,includeL)
c, & !< constituent number c, & !< constituent number
s s
do c = 1,homogenization_Ngrains(material_homogenizationAt(e)) do c = 1,homogenization_Nconstituents(material_homogenizationAt(e))
if (includeL) then if (includeL) then
crystallite_Lp(1:3,1:3,c,i,e) = crystallite_partitionedLp0(1:3,1:3,c,i,e) crystallite_Lp(1:3,1:3,c,i,e) = crystallite_partitionedLp0(1:3,1:3,c,i,e)
crystallite_Li(1:3,1:3,c,i,e) = crystallite_partitionedLi0(1:3,1:3,c,i,e) crystallite_Li(1:3,1:3,c,i,e) = crystallite_partitionedLi0(1:3,1:3,c,i,e)
@ -687,7 +687,7 @@ subroutine crystallite_orientations
!$OMP PARALLEL DO !$OMP PARALLEL DO
do e = FEsolving_execElem(1),FEsolving_execElem(2) do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1),FEsolving_execIP(2) do i = FEsolving_execIP(1),FEsolving_execIP(2)
do c = 1,homogenization_Ngrains(material_homogenizationAt(e)) do c = 1,homogenization_Nconstituents(material_homogenizationAt(e))
call crystallite_orientation(c,i,e)%fromMatrix(transpose(math_rotationalPart(crystallite_Fe(1:3,1:3,c,i,e)))) call crystallite_orientation(c,i,e)%fromMatrix(transpose(math_rotationalPart(crystallite_Fe(1:3,1:3,c,i,e))))
enddo; enddo; enddo enddo; enddo; enddo
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
@ -811,11 +811,11 @@ subroutine crystallite_results
real(pReal), allocatable, dimension(:,:,:) :: select_tensors real(pReal), allocatable, dimension(:,:,:) :: select_tensors
integer :: e,i,c,j integer :: e,i,c,j
allocate(select_tensors(3,3,count(material_phaseAt==instance)*discretization_nIP)) allocate(select_tensors(3,3,count(material_phaseAt==instance)*discretization_nIPs))
j=0 j=0
do e = 1, size(material_phaseAt,2) do e = 1, size(material_phaseAt,2)
do i = 1, discretization_nIP do i = 1, discretization_nIPs
do c = 1, size(material_phaseAt,1) !ToDo: this needs to be changed for varying Ngrains do c = 1, size(material_phaseAt,1) !ToDo: this needs to be changed for varying Ngrains
if (material_phaseAt(c,e) == instance) then if (material_phaseAt(c,e) == instance) then
j = j + 1 j = j + 1
@ -838,11 +838,11 @@ subroutine crystallite_results
type(rotation), allocatable, dimension(:) :: select_rotations type(rotation), allocatable, dimension(:) :: select_rotations
integer :: e,i,c,j integer :: e,i,c,j
allocate(select_rotations(count(material_phaseAt==instance)*homogenization_maxNgrains*discretization_nIP)) allocate(select_rotations(count(material_phaseAt==instance)*homogenization_maxNconstituents*discretization_nIPs))
j=0 j=0
do e = 1, size(material_phaseAt,2) do e = 1, size(material_phaseAt,2)
do i = 1, discretization_nIP do i = 1, discretization_nIPs
do c = 1, size(material_phaseAt,1) !ToDo: this needs to be changed for varying Ngrains do c = 1, size(material_phaseAt,1) !ToDo: this needs to be changed for varying Ngrains
if (material_phaseAt(c,e) == instance) then if (material_phaseAt(c,e) == instance) then
j = j + 1 j = j + 1
@ -1561,7 +1561,7 @@ subroutine crystallite_restartWrite
call HDF5_closeGroup(groupHandle) call HDF5_closeGroup(groupHandle)
groupHandle = HDF5_addGroup(fileHandle,'materialpoint') groupHandle = HDF5_addGroup(fileHandle,'materialpoint')
do i = 1, material_Nhomogenization do i = 1, size(material_name_homogenization)
write(datasetName,'(i0,a)') i,'_omega_homogenization' write(datasetName,'(i0,a)') i,'_omega_homogenization'
call HDF5_write(groupHandle,homogState(i)%state,datasetName) call HDF5_write(groupHandle,homogState(i)%state,datasetName)
enddo enddo
@ -1602,7 +1602,7 @@ subroutine crystallite_restartRead
call HDF5_closeGroup(groupHandle) call HDF5_closeGroup(groupHandle)
groupHandle = HDF5_openGroup(fileHandle,'materialpoint') groupHandle = HDF5_openGroup(fileHandle,'materialpoint')
do i = 1, material_Nhomogenization do i = 1,size(material_name_homogenization)
write(datasetName,'(i0,a)') i,'_omega_homogenization' write(datasetName,'(i0,a)') i,'_omega_homogenization'
call HDF5_read(groupHandle,homogState(i)%state0,datasetName) call HDF5_read(groupHandle,homogState(i)%state0,datasetName)
enddo enddo
@ -1635,7 +1635,7 @@ subroutine crystallite_forward
do j = 1,phase_Nsources(i) do j = 1,phase_Nsources(i)
sourceState(i)%p(j)%state0 = sourceState(i)%p(j)%state sourceState(i)%p(j)%state0 = sourceState(i)%p(j)%state
enddo; enddo enddo; enddo
do i = 1, material_Nhomogenization do i = 1,size(material_name_homogenization)
homogState (i)%state0 = homogState (i)%state homogState (i)%state0 = homogState (i)%state
thermalState(i)%state0 = thermalState(i)%state thermalState(i)%state0 = thermalState(i)%state
damageState (i)%state0 = damageState (i)%state damageState (i)%state0 = damageState (i)%state

18
src/damage_local.f90
View File

@ -42,7 +42,7 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine damage_local_init subroutine damage_local_init
integer :: Ninstance,NofMyHomog,h integer :: Ninstances,Nmaterialpoints,h
class(tNode), pointer :: & class(tNode), pointer :: &
num_generic, & num_generic, &
material_homogenization, & material_homogenization, &
@ -57,8 +57,8 @@ subroutine damage_local_init
num%residualStiffness = num_generic%get_asFloat('residualStiffness', defaultVal=1.0e-6_pReal) num%residualStiffness = num_generic%get_asFloat('residualStiffness', defaultVal=1.0e-6_pReal)
if (num%residualStiffness < 0.0_pReal) call IO_error(301,ext_msg='residualStiffness') if (num%residualStiffness < 0.0_pReal) call IO_error(301,ext_msg='residualStiffness')
Ninstance = count(damage_type == DAMAGE_local_ID) Ninstances = count(damage_type == DAMAGE_local_ID)
allocate(param(Ninstance)) allocate(param(Ninstances))
material_homogenization => config_material%get('homogenization') material_homogenization => config_material%get('homogenization')
do h = 1, material_homogenization%length do h = 1, material_homogenization%length
@ -73,11 +73,11 @@ subroutine damage_local_init
prm%output = homogDamage%get_asStrings('output',defaultVal=emptyStringArray) prm%output = homogDamage%get_asStrings('output',defaultVal=emptyStringArray)
#endif #endif
NofMyHomog = count(material_homogenizationAt == h) Nmaterialpoints = count(material_homogenizationAt == h)
damageState(h)%sizeState = 1 damageState(h)%sizeState = 1
allocate(damageState(h)%state0 (1,NofMyHomog), source=damage_initialPhi(h)) allocate(damageState(h)%state0 (1,Nmaterialpoints), source=damage_initialPhi(h))
allocate(damageState(h)%subState0(1,NofMyHomog), source=damage_initialPhi(h)) allocate(damageState(h)%subState0(1,Nmaterialpoints), source=damage_initialPhi(h))
allocate(damageState(h)%state (1,NofMyHomog), source=damage_initialPhi(h)) allocate(damageState(h)%state (1,Nmaterialpoints), source=damage_initialPhi(h))
nullify(damageMapping(h)%p) nullify(damageMapping(h)%p)
damageMapping(h)%p => material_homogenizationMemberAt damageMapping(h)%p => material_homogenizationMemberAt
@ -143,8 +143,8 @@ subroutine damage_local_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el
call constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el) call constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el)
phiDot = phiDot/real(homogenization_Ngrains(material_homogenizationAt(el)),pReal) phiDot = phiDot/real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal)
dPhiDot_dPhi = dPhiDot_dPhi/real(homogenization_Ngrains(material_homogenizationAt(el)),pReal) dPhiDot_dPhi = dPhiDot_dPhi/real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal)
end subroutine damage_local_getSourceAndItsTangent end subroutine damage_local_getSourceAndItsTangent

12
src/damage_none.f90
View File

@ -16,18 +16,18 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine damage_none_init subroutine damage_none_init
integer :: h,NofMyHomog integer :: h,Nmaterialpoints
print'(/,a)', ' <<<+- damage_none init -+>>>'; flush(6) print'(/,a)', ' <<<+- damage_none init -+>>>'; flush(6)
do h = 1, material_Nhomogenization do h = 1, size(material_name_homogenization)
if (damage_type(h) /= DAMAGE_NONE_ID) cycle if (damage_type(h) /= DAMAGE_NONE_ID) cycle
NofMyHomog = count(material_homogenizationAt == h) Nmaterialpoints = count(material_homogenizationAt == h)
damageState(h)%sizeState = 0 damageState(h)%sizeState = 0
allocate(damageState(h)%state0 (0,NofMyHomog)) allocate(damageState(h)%state0 (0,Nmaterialpoints))
allocate(damageState(h)%subState0(0,NofMyHomog)) allocate(damageState(h)%subState0(0,Nmaterialpoints))
allocate(damageState(h)%state (0,NofMyHomog)) allocate(damageState(h)%state (0,Nmaterialpoints))
deallocate(damage(h)%p) deallocate(damage(h)%p)
allocate (damage(h)%p(1), source=damage_initialPhi(h)) allocate (damage(h)%p(1), source=damage_initialPhi(h))

26
src/damage_nonlocal.f90
View File

@ -46,7 +46,7 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine damage_nonlocal_init subroutine damage_nonlocal_init
integer :: Ninstance,NofMyHomog,h integer :: Ninstances,Nmaterialpoints,h
class(tNode), pointer :: & class(tNode), pointer :: &
num_generic, & num_generic, &
material_homogenization, & material_homogenization, &
@ -60,8 +60,8 @@ subroutine damage_nonlocal_init
num_generic => config_numerics%get('generic',defaultVal= emptyDict) num_generic => config_numerics%get('generic',defaultVal= emptyDict)
num%charLength = num_generic%get_asFloat('charLength',defaultVal=1.0_pReal) num%charLength = num_generic%get_asFloat('charLength',defaultVal=1.0_pReal)
Ninstance = count(damage_type == DAMAGE_nonlocal_ID) Ninstances = count(damage_type == DAMAGE_nonlocal_ID)
allocate(param(Ninstance)) allocate(param(Ninstances))
material_homogenization => config_material%get('homogenization') material_homogenization => config_material%get('homogenization')
do h = 1, material_homogenization%length do h = 1, material_homogenization%length
@ -76,11 +76,11 @@ subroutine damage_nonlocal_init
prm%output = homogDamage%get_asStrings('output',defaultVal=emptyStringArray) prm%output = homogDamage%get_asStrings('output',defaultVal=emptyStringArray)
#endif #endif
NofMyHomog = count(material_homogenizationAt == h) Nmaterialpoints = count(material_homogenizationAt == h)
damageState(h)%sizeState = 1 damageState(h)%sizeState = 1
allocate(damageState(h)%state0 (1,NofMyHomog), source=damage_initialPhi(h)) allocate(damageState(h)%state0 (1,Nmaterialpoints), source=damage_initialPhi(h))
allocate(damageState(h)%subState0(1,NofMyHomog), source=damage_initialPhi(h)) allocate(damageState(h)%subState0(1,Nmaterialpoints), source=damage_initialPhi(h))
allocate(damageState(h)%state (1,NofMyHomog), source=damage_initialPhi(h)) allocate(damageState(h)%state (1,Nmaterialpoints), source=damage_initialPhi(h))
nullify(damageMapping(h)%p) nullify(damageMapping(h)%p)
damageMapping(h)%p => material_homogenizationMemberAt damageMapping(h)%p => material_homogenizationMemberAt
@ -110,8 +110,8 @@ subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip,
dPhiDot_dPhi = 0.0_pReal dPhiDot_dPhi = 0.0_pReal
call constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el) call constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el)
phiDot = phiDot/real(homogenization_Ngrains(material_homogenizationAt(el)),pReal) phiDot = phiDot/real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal)
dPhiDot_dPhi = dPhiDot_dPhi/real(homogenization_Ngrains(material_homogenizationAt(el)),pReal) dPhiDot_dPhi = dPhiDot_dPhi/real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal)
end subroutine damage_nonlocal_getSourceAndItsTangent end subroutine damage_nonlocal_getSourceAndItsTangent
@ -132,13 +132,13 @@ function damage_nonlocal_getDiffusion(ip,el)
homog = material_homogenizationAt(el) homog = material_homogenizationAt(el)
damage_nonlocal_getDiffusion = 0.0_pReal damage_nonlocal_getDiffusion = 0.0_pReal
do grain = 1, homogenization_Ngrains(homog) do grain = 1, homogenization_Nconstituents(homog)
damage_nonlocal_getDiffusion = damage_nonlocal_getDiffusion + & damage_nonlocal_getDiffusion = damage_nonlocal_getDiffusion + &
crystallite_push33ToRef(grain,ip,el,lattice_D(1:3,1:3,material_phaseAt(grain,el))) crystallite_push33ToRef(grain,ip,el,lattice_D(1:3,1:3,material_phaseAt(grain,el)))
enddo enddo
damage_nonlocal_getDiffusion = & damage_nonlocal_getDiffusion = &
num%charLength**2*damage_nonlocal_getDiffusion/real(homogenization_Ngrains(homog),pReal) num%charLength**2*damage_nonlocal_getDiffusion/real(homogenization_Nconstituents(homog),pReal)
end function damage_nonlocal_getDiffusion end function damage_nonlocal_getDiffusion
@ -156,12 +156,12 @@ real(pReal) function damage_nonlocal_getMobility(ip,el)
damage_nonlocal_getMobility = 0.0_pReal damage_nonlocal_getMobility = 0.0_pReal
do ipc = 1, homogenization_Ngrains(material_homogenizationAt(el)) do ipc = 1, homogenization_Nconstituents(material_homogenizationAt(el))
damage_nonlocal_getMobility = damage_nonlocal_getMobility + lattice_M(material_phaseAt(ipc,el)) damage_nonlocal_getMobility = damage_nonlocal_getMobility + lattice_M(material_phaseAt(ipc,el))
enddo enddo
damage_nonlocal_getMobility = damage_nonlocal_getMobility/& damage_nonlocal_getMobility = damage_nonlocal_getMobility/&
real(homogenization_Ngrains(material_homogenizationAt(el)),pReal) real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal)
end function damage_nonlocal_getMobility end function damage_nonlocal_getMobility

8
src/discretization.f90
View File

@ -11,8 +11,8 @@ module discretization
private private
integer, public, protected :: & integer, public, protected :: &
discretization_nIP, & discretization_nIPs, &
discretization_nElem discretization_Nelems
integer, public, protected, dimension(:), allocatable :: & integer, public, protected, dimension(:), allocatable :: &
discretization_materialAt discretization_materialAt
@ -51,8 +51,8 @@ subroutine discretization_init(materialAt,&
print'(/,a)', ' <<<+- discretization init -+>>>'; flush(6) print'(/,a)', ' <<<+- discretization init -+>>>'; flush(6)
discretization_nElem = size(materialAt,1) discretization_Nelems = size(materialAt,1)
discretization_nIP = size(IPcoords0,2)/discretization_nElem discretization_nIPs = size(IPcoords0,2)/discretization_Nelems
discretization_materialAt = materialAt discretization_materialAt = materialAt

2
src/grid/DAMASK_grid.f90
View File

@ -210,7 +210,7 @@ program DAMASK_grid
if(.not. step_discretization%contains('t')) call IO_error(error_ID=837,ext_msg = 't missing') if(.not. step_discretization%contains('t')) call IO_error(error_ID=837,ext_msg = 't missing')
if(.not. step_discretization%contains('N')) call IO_error(error_ID=837,ext_msg = 'N missing') if(.not. step_discretization%contains('N')) call IO_error(error_ID=837,ext_msg = 'N missing')
loadCases(l)%time = step_discretization%get_asFloat('t') loadCases(l)%time = step_discretization%get_asFloat('t')
loadCases(l)%incs = step_discretization%get_asFloat('N') loadCases(l)%incs = step_discretization%get_asInt ('N')
loadCases(l)%logscale = step_discretization%get_asBool ('log_timestep', defaultVal= .false.) loadCases(l)%logscale = step_discretization%get_asBool ('log_timestep', defaultVal= .false.)
loadCases(l)%outputfrequency = step_discretization%get_asInt ('f_out', defaultVal=1) loadCases(l)%outputfrequency = step_discretization%get_asInt ('f_out', defaultVal=1)
loadCases(l)%restartfrequency = step_discretization%get_asInt ('f_restart', defaultVal=huge(0)) loadCases(l)%restartfrequency = step_discretization%get_asInt ('f_restart', defaultVal=huge(0))

30
src/grid/grid_mech_FEM.f90
View File

@ -236,7 +236,7 @@ subroutine grid_mech_FEM_init
F = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) F = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3)
endif restartRead endif restartRead
materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent homogenization_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent
call utilities_updateCoords(F) call utilities_updateCoords(F)
call utilities_constitutiveResponse(P_current,temp33_Real,C_volAvg,devNull, & ! stress field, stress avg, global average of stiffness and (min+max)/2 call utilities_constitutiveResponse(P_current,temp33_Real,C_volAvg,devNull, & ! stress field, stress avg, global average of stiffness and (min+max)/2
F, & ! target F F, & ! target F
@ -364,7 +364,7 @@ subroutine grid_mech_FEM_forward(cutBack,guess,timeinc,timeinc_old,loadCaseTime,
F_lastInc = F F_lastInc = F
materialpoint_F0 = reshape(F, [3,3,1,product(grid(1:2))*grid3]) homogenization_F0 = reshape(F, [3,3,1,product(grid(1:2))*grid3])
endif endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -557,9 +557,9 @@ subroutine formResidual(da_local,x_local, &
ii = i-xstart+1; jj = j-ystart+1; kk = k-zstart+1 ii = i-xstart+1; jj = j-ystart+1; kk = k-zstart+1
ele = ele + 1 ele = ele + 1
f_elem = matmul(transpose(BMat),transpose(P_current(1:3,1:3,ii,jj,kk)))*detJ + & f_elem = matmul(transpose(BMat),transpose(P_current(1:3,1:3,ii,jj,kk)))*detJ + &
matmul(HGMat,x_elem)*(materialpoint_dPdF(1,1,1,1,1,ele) + & matmul(HGMat,x_elem)*(homogenization_dPdF(1,1,1,1,1,ele) + &
materialpoint_dPdF(2,2,2,2,1,ele) + & homogenization_dPdF(2,2,2,2,1,ele) + &
materialpoint_dPdF(3,3,3,3,1,ele))/3.0_pReal homogenization_dPdF(3,3,3,3,1,ele))/3.0_pReal
ctr = 0 ctr = 0
do kk = 0, 1; do jj = 0, 1; do ii = 0, 1 do kk = 0, 1; do jj = 0, 1; do ii = 0, 1
ctr = ctr + 1 ctr = ctr + 1
@ -636,18 +636,18 @@ subroutine formJacobian(da_local,x_local,Jac_pre,Jac,dummy,ierr)
row = col row = col
ele = ele + 1 ele = ele + 1
K_ele = 0.0 K_ele = 0.0
K_ele(1 :8 ,1 :8 ) = HGMat*(materialpoint_dPdF(1,1,1,1,1,ele) + & K_ele(1 :8 ,1 :8 ) = HGMat*(homogenization_dPdF(1,1,1,1,1,ele) + &
materialpoint_dPdF(2,2,2,2,1,ele) + & homogenization_dPdF(2,2,2,2,1,ele) + &
materialpoint_dPdF(3,3,3,3,1,ele))/3.0_pReal homogenization_dPdF(3,3,3,3,1,ele))/3.0_pReal
K_ele(9 :16,9 :16) = HGMat*(materialpoint_dPdF(1,1,1,1,1,ele) + & K_ele(9 :16,9 :16) = HGMat*(homogenization_dPdF(1,1,1,1,1,ele) + &
materialpoint_dPdF(2,2,2,2,1,ele) + & homogenization_dPdF(2,2,2,2,1,ele) + &
materialpoint_dPdF(3,3,3,3,1,ele))/3.0_pReal homogenization_dPdF(3,3,3,3,1,ele))/3.0_pReal
K_ele(17:24,17:24) = HGMat*(materialpoint_dPdF(1,1,1,1,1,ele) + & K_ele(17:24,17:24) = HGMat*(homogenization_dPdF(1,1,1,1,1,ele) + &
materialpoint_dPdF(2,2,2,2,1,ele) + & homogenization_dPdF(2,2,2,2,1,ele) + &
materialpoint_dPdF(3,3,3,3,1,ele))/3.0_pReal homogenization_dPdF(3,3,3,3,1,ele))/3.0_pReal
K_ele = K_ele + & K_ele = K_ele + &
matmul(transpose(BMatFull), & matmul(transpose(BMatFull), &
matmul(reshape(reshape(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,ele), & matmul(reshape(reshape(homogenization_dPdF(1:3,1:3,1:3,1:3,1,ele), &
shape=[3,3,3,3], order=[2,1,4,3]),shape=[9,9]),BMatFull))*detJ shape=[3,3,3,3], order=[2,1,4,3]),shape=[9,9]),BMatFull))*detJ
call MatSetValuesStencil(Jac,24,row,24,col,K_ele,ADD_VALUES,ierr) call MatSetValuesStencil(Jac,24,row,24,col,K_ele,ADD_VALUES,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)

4
src/grid/grid_mech_spectral_basic.f90
View File

@ -198,7 +198,7 @@ subroutine grid_mech_spectral_basic_init
F = reshape(F_lastInc,[9,grid(1),grid(2),grid3]) F = reshape(F_lastInc,[9,grid(1),grid(2),grid3])
endif restartRead endif restartRead
materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent homogenization_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent
call utilities_updateCoords(reshape(F,shape(F_lastInc))) call utilities_updateCoords(reshape(F,shape(F_lastInc)))
call utilities_constitutiveResponse(P,temp33_Real,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2 call utilities_constitutiveResponse(P,temp33_Real,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2
reshape(F,shape(F_lastInc)), & ! target F reshape(F,shape(F_lastInc)), & ! target F
@ -324,7 +324,7 @@ subroutine grid_mech_spectral_basic_forward(cutBack,guess,timeinc,timeinc_old,lo
rotation_BC%rotate(F_aimDot,active=.true.)) rotation_BC%rotate(F_aimDot,active=.true.))
F_lastInc = reshape(F,[3,3,grid(1),grid(2),grid3]) F_lastInc = reshape(F,[3,3,grid(1),grid(2),grid3])
materialpoint_F0 = reshape(F, [3,3,1,product(grid(1:2))*grid3]) homogenization_F0 = reshape(F, [3,3,1,product(grid(1:2))*grid3])
endif endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------

6
src/grid/grid_mech_spectral_polarisation.f90
View File

@ -224,7 +224,7 @@ subroutine grid_mech_spectral_polarisation_init
F_tau_lastInc = 2.0_pReal*F_lastInc F_tau_lastInc = 2.0_pReal*F_lastInc
endif restartRead endif restartRead
materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent homogenization_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent
call utilities_updateCoords(reshape(F,shape(F_lastInc))) call utilities_updateCoords(reshape(F,shape(F_lastInc)))
call utilities_constitutiveResponse(P,temp33_Real,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2 call utilities_constitutiveResponse(P,temp33_Real,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2
reshape(F,shape(F_lastInc)), & ! target F reshape(F,shape(F_lastInc)), & ! target F
@ -364,7 +364,7 @@ subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,timeinc,timeinc
F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid3]) F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid3])
F_tau_lastInc = reshape(F_tau,[3,3,grid(1),grid(2),grid3]) F_tau_lastInc = reshape(F_tau,[3,3,grid(1),grid(2),grid3])
materialpoint_F0 = reshape(F,[3,3,1,product(grid(1:2))*grid3]) homogenization_F0 = reshape(F,[3,3,1,product(grid(1:2))*grid3])
endif endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -606,7 +606,7 @@ subroutine formResidual(in, FandF_tau, &
do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1)
e = e + 1 e = e + 1
residual_F(1:3,1:3,i,j,k) = & residual_F(1:3,1:3,i,j,k) = &
math_mul3333xx33(math_invSym3333(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,e) + C_scale), & math_mul3333xx33(math_invSym3333(homogenization_dPdF(1:3,1:3,1:3,1:3,1,e) + C_scale), &
residual_F(1:3,1:3,i,j,k) - matmul(F(1:3,1:3,i,j,k), & residual_F(1:3,1:3,i,j,k) - matmul(F(1:3,1:3,i,j,k), &
math_mul3333xx33(C_scale,F_tau(1:3,1:3,i,j,k) - F(1:3,1:3,i,j,k) - math_I3))) & math_mul3333xx33(C_scale,F_tau(1:3,1:3,i,j,k) - F(1:3,1:3,i,j,k) - math_I3))) &
+ residual_F_tau(1:3,1:3,i,j,k) + residual_F_tau(1:3,1:3,i,j,k)

20
src/grid/spectral_utilities.f90
View File

@ -802,7 +802,7 @@ end subroutine utilities_fourierTensorDivergence
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculate constitutive response from materialpoint_F0 to F during timeinc !> @brief calculate constitutive response from homogenization_F0 to F during timeinc
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,& subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
F,timeinc,rotation_BC) F,timeinc,rotation_BC)
@ -824,11 +824,11 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
print'(/,a)', ' ... evaluating constitutive response ......................................' print'(/,a)', ' ... evaluating constitutive response ......................................'
flush(IO_STDOUT) flush(IO_STDOUT)
materialpoint_F = reshape(F,[3,3,1,product(grid(1:2))*grid3]) ! set materialpoint target F to estimated field homogenization_F = reshape(F,[3,3,1,product(grid(1:2))*grid3]) ! set materialpoint target F to estimated field
call materialpoint_stressAndItsTangent(timeinc) ! calculate P field call materialpoint_stressAndItsTangent(timeinc) ! calculate P field
P = reshape(materialpoint_P, [3,3,grid(1),grid(2),grid3]) P = reshape(homogenization_P, [3,3,grid(1),grid(2),grid3])
P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt ! average of P P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt ! average of P
call MPI_Allreduce(MPI_IN_PLACE,P_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,P_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
if (debugRotation) & if (debugRotation) &
@ -845,13 +845,13 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
dPdF_min = huge(1.0_pReal) dPdF_min = huge(1.0_pReal)
dPdF_norm_min = huge(1.0_pReal) dPdF_norm_min = huge(1.0_pReal)
do i = 1, product(grid(1:2))*grid3 do i = 1, product(grid(1:2))*grid3
if (dPdF_norm_max < sum(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,i)**2.0_pReal)) then if (dPdF_norm_max < sum(homogenization_dPdF(1:3,1:3,1:3,1:3,1,i)**2.0_pReal)) then
dPdF_max = materialpoint_dPdF(1:3,1:3,1:3,1:3,1,i) dPdF_max = homogenization_dPdF(1:3,1:3,1:3,1:3,1,i)
dPdF_norm_max = sum(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,i)**2.0_pReal) dPdF_norm_max = sum(homogenization_dPdF(1:3,1:3,1:3,1:3,1,i)**2.0_pReal)
endif endif
if (dPdF_norm_min > sum(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,i)**2.0_pReal)) then if (dPdF_norm_min > sum(homogenization_dPdF(1:3,1:3,1:3,1:3,1,i)**2.0_pReal)) then
dPdF_min = materialpoint_dPdF(1:3,1:3,1:3,1:3,1,i) dPdF_min = homogenization_dPdF(1:3,1:3,1:3,1:3,1,i)
dPdF_norm_min = sum(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,i)**2.0_pReal) dPdF_norm_min = sum(homogenization_dPdF(1:3,1:3,1:3,1:3,1,i)**2.0_pReal)
endif endif
end do end do
@ -869,7 +869,7 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
C_minmaxAvg = 0.5_pReal*(dPdF_max + dPdF_min) C_minmaxAvg = 0.5_pReal*(dPdF_max + dPdF_min)
C_volAvg = sum(sum(materialpoint_dPdF,dim=6),dim=5) C_volAvg = sum(sum(homogenization_dPdF,dim=6),dim=5)
call MPI_Allreduce(MPI_IN_PLACE,C_volAvg,81,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,C_volAvg,81,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
C_volAvg = C_volAvg * wgt C_volAvg = C_volAvg * wgt

76
src/homogenization.f90
View File

@ -31,12 +31,12 @@ module homogenization
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! General variables for the homogenization at a material point ! General variables for the homogenization at a material point
real(pReal), dimension(:,:,:,:), allocatable, public :: & real(pReal), dimension(:,:,:,:), allocatable, public :: &
materialpoint_F0, & !< def grad of IP at start of FE increment homogenization_F0, & !< def grad of IP at start of FE increment
materialpoint_F !< def grad of IP to be reached at end of FE increment homogenization_F !< def grad of IP to be reached at end of FE increment
real(pReal), dimension(:,:,:,:), allocatable, public, protected :: & real(pReal), dimension(:,:,:,:), allocatable, public, protected :: &
materialpoint_P !< first P--K stress of IP homogenization_P !< first P--K stress of IP
real(pReal), dimension(:,:,:,:,:,:), allocatable, public, protected :: & real(pReal), dimension(:,:,:,:,:,:), allocatable, public, protected :: &
materialpoint_dPdF !< tangent of first P--K stress at IP homogenization_dPdF !< tangent of first P--K stress at IP
type :: tNumerics type :: tNumerics
integer :: & integer :: &
@ -158,7 +158,7 @@ subroutine homogenization_init
debugHomog%grain = config_debug%get_asInt('grain',defaultVal = 1) debugHomog%grain = config_debug%get_asInt('grain',defaultVal = 1)
if (debugHomog%grain < 1 & if (debugHomog%grain < 1 &
.or. debugHomog%grain > homogenization_Ngrains(material_homogenizationAt(debugHomog%element))) & .or. debugHomog%grain > homogenization_Nconstituents(material_homogenizationAt(debugHomog%element))) &
call IO_error(602,ext_msg='constituent', el=debugHomog%element, g=debugHomog%grain) call IO_error(602,ext_msg='constituent', el=debugHomog%element, g=debugHomog%grain)
@ -181,10 +181,10 @@ subroutine homogenization_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocate and initialize global variables ! allocate and initialize global variables
allocate(materialpoint_dPdF(3,3,3,3,discretization_nIP,discretization_nElem), source=0.0_pReal) allocate(homogenization_dPdF(3,3,3,3,discretization_nIPs,discretization_Nelems), source=0.0_pReal)
materialpoint_F0 = spread(spread(math_I3,3,discretization_nIP),4,discretization_nElem) ! initialize to identity homogenization_F0 = spread(spread(math_I3,3,discretization_nIPs),4,discretization_Nelems) ! initialize to identity
materialpoint_F = materialpoint_F0 ! initialize to identity homogenization_F = homogenization_F0 ! initialize to identity
allocate(materialpoint_P(3,3,discretization_nIP,discretization_nElem), source=0.0_pReal) allocate(homogenization_P(3,3,discretization_nIPs,discretization_Nelems), source=0.0_pReal)
print'(/,a)', ' <<<+- homogenization init -+>>>'; flush(IO_STDOUT) print'(/,a)', ' <<<+- homogenization init -+>>>'; flush(IO_STDOUT)
@ -213,13 +213,13 @@ subroutine materialpoint_stressAndItsTangent(dt)
i, & !< integration point number i, & !< integration point number
e, & !< element number e, & !< element number
myNgrains myNgrains
real(pReal), dimension(discretization_nIP,discretization_nElem) :: & real(pReal), dimension(discretization_nIPs,discretization_Nelems) :: &
subFrac, & subFrac, &
subStep subStep
logical, dimension(discretization_nIP,discretization_nElem) :: & logical, dimension(discretization_nIPs,discretization_Nelems) :: &
requested, & requested, &
converged converged
logical, dimension(2,discretization_nIP,discretization_nElem) :: & logical, dimension(2,discretization_nIPs,discretization_Nelems) :: &
doneAndHappy doneAndHappy
@ -257,7 +257,7 @@ subroutine materialpoint_stressAndItsTangent(dt)
!$OMP PARALLEL DO !$OMP PARALLEL DO
elementLooping1: do e = FEsolving_execElem(1),FEsolving_execElem(2) elementLooping1: do e = FEsolving_execElem(1),FEsolving_execElem(2)
myNgrains = homogenization_Ngrains(material_homogenizationAt(e)) myNgrains = homogenization_Nconstituents(material_homogenizationAt(e))
IpLooping1: do i = FEsolving_execIP(1),FEsolving_execIP(2) IpLooping1: do i = FEsolving_execIP(1),FEsolving_execIP(2)
if (converged(i,e)) then if (converged(i,e)) then
@ -327,11 +327,11 @@ subroutine materialpoint_stressAndItsTangent(dt)
! deformation partitioning ! deformation partitioning
!$OMP PARALLEL DO PRIVATE(myNgrains) !$OMP PARALLEL DO PRIVATE(myNgrains)
elementLooping2: do e = FEsolving_execElem(1),FEsolving_execElem(2) elementLooping2: do e = FEsolving_execElem(1),FEsolving_execElem(2)
myNgrains = homogenization_Ngrains(material_homogenizationAt(e)) myNgrains = homogenization_Nconstituents(material_homogenizationAt(e))
IpLooping2: do i = FEsolving_execIP(1),FEsolving_execIP(2) IpLooping2: do i = FEsolving_execIP(1),FEsolving_execIP(2)
if(requested(i,e) .and. .not. doneAndHappy(1,i,e)) then ! requested but not yet done if(requested(i,e) .and. .not. doneAndHappy(1,i,e)) then ! requested but not yet done
call partitionDeformation(materialpoint_F0(1:3,1:3,i,e) & call partitionDeformation(homogenization_F0(1:3,1:3,i,e) &
+ (materialpoint_F(1:3,1:3,i,e)-materialpoint_F0(1:3,1:3,i,e))& + (homogenization_F(1:3,1:3,i,e)-homogenization_F0(1:3,1:3,i,e))&
*(subStep(i,e)+subFrac(i,e)), & *(subStep(i,e)+subFrac(i,e)), &
i,e) i,e)
crystallite_dt(1:myNgrains,i,e) = dt*subStep(i,e) ! propagate materialpoint dt to grains crystallite_dt(1:myNgrains,i,e) = dt*subStep(i,e) ! propagate materialpoint dt to grains
@ -357,8 +357,8 @@ subroutine materialpoint_stressAndItsTangent(dt)
doneAndHappy(1:2,i,e) = [.true.,.false.] doneAndHappy(1:2,i,e) = [.true.,.false.]
else else
doneAndHappy(1:2,i,e) = updateState(dt*subStep(i,e), & doneAndHappy(1:2,i,e) = updateState(dt*subStep(i,e), &
materialpoint_F0(1:3,1:3,i,e) & homogenization_F0(1:3,1:3,i,e) &
+ (materialpoint_F(1:3,1:3,i,e)-materialpoint_F0(1:3,1:3,i,e)) & + (homogenization_F(1:3,1:3,i,e)-homogenization_F0(1:3,1:3,i,e)) &
*(subStep(i,e)+subFrac(i,e)), & *(subStep(i,e)+subFrac(i,e)), &
i,e) i,e)
converged(i,e) = all(doneAndHappy(1:2,i,e)) ! converged if done and happy converged(i,e) = all(doneAndHappy(1:2,i,e)) ! converged if done and happy
@ -408,12 +408,12 @@ subroutine partitionDeformation(subF,ip,el)
case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization
call mech_isostrain_partitionDeformation(& call mech_isostrain_partitionDeformation(&
crystallite_partitionedF(1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el), & crystallite_partitionedF(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), &
subF) subF)
case (HOMOGENIZATION_RGC_ID) chosenHomogenization case (HOMOGENIZATION_RGC_ID) chosenHomogenization
call mech_RGC_partitionDeformation(& call mech_RGC_partitionDeformation(&
crystallite_partitionedF(1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el), & crystallite_partitionedF(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), &
subF,& subF,&
ip, & ip, &
el) el)
@ -437,19 +437,19 @@ function updateState(subdt,subF,ip,el)
el !< element number el !< element number
integer :: c integer :: c
logical, dimension(2) :: updateState logical, dimension(2) :: updateState
real(pReal) :: dPdFs(3,3,3,3,homogenization_Ngrains(material_homogenizationAt(el))) real(pReal) :: dPdFs(3,3,3,3,homogenization_Nconstituents(material_homogenizationAt(el)))
updateState = .true. updateState = .true.
chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el))) chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el)))
case (HOMOGENIZATION_RGC_ID) chosenHomogenization case (HOMOGENIZATION_RGC_ID) chosenHomogenization
do c=1,homogenization_Ngrains(material_homogenizationAt(el)) do c=1,homogenization_Nconstituents(material_homogenizationAt(el))
dPdFs(:,:,:,:,c) = crystallite_stressTangent(c,ip,el) dPdFs(:,:,:,:,c) = crystallite_stressTangent(c,ip,el)
enddo enddo
updateState = & updateState = &
updateState .and. & updateState .and. &
mech_RGC_updateState(crystallite_P(1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el), & mech_RGC_updateState(crystallite_P(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), &
crystallite_partitionedF(1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el), & crystallite_partitionedF(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), &
crystallite_partitionedF0(1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el),& crystallite_partitionedF0(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el),&
subF,& subF,&
subdt, & subdt, &
dPdFs, & dPdFs, &
@ -487,33 +487,33 @@ subroutine averageStressAndItsTangent(ip,el)
ip, & !< integration point ip, & !< integration point
el !< element number el !< element number
integer :: c integer :: c
real(pReal) :: dPdFs(3,3,3,3,homogenization_Ngrains(material_homogenizationAt(el))) real(pReal) :: dPdFs(3,3,3,3,homogenization_Nconstituents(material_homogenizationAt(el)))
chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el))) chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el)))
case (HOMOGENIZATION_NONE_ID) chosenHomogenization case (HOMOGENIZATION_NONE_ID) chosenHomogenization
materialpoint_P(1:3,1:3,ip,el) = crystallite_P(1:3,1:3,1,ip,el) homogenization_P(1:3,1:3,ip,el) = crystallite_P(1:3,1:3,1,ip,el)
materialpoint_dPdF(1:3,1:3,1:3,1:3,ip,el) = crystallite_stressTangent(1,ip,el) homogenization_dPdF(1:3,1:3,1:3,1:3,ip,el) = crystallite_stressTangent(1,ip,el)
case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization
do c = 1, homogenization_Ngrains(material_homogenizationAt(el)) do c = 1, homogenization_Nconstituents(material_homogenizationAt(el))
dPdFs(:,:,:,:,c) = crystallite_stressTangent(c,ip,el) dPdFs(:,:,:,:,c) = crystallite_stressTangent(c,ip,el)
enddo enddo
call mech_isostrain_averageStressAndItsTangent(& call mech_isostrain_averageStressAndItsTangent(&
materialpoint_P(1:3,1:3,ip,el), & homogenization_P(1:3,1:3,ip,el), &
materialpoint_dPdF(1:3,1:3,1:3,1:3,ip,el),& homogenization_dPdF(1:3,1:3,1:3,1:3,ip,el),&
crystallite_P(1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el), & crystallite_P(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), &
dPdFs, & dPdFs, &
homogenization_typeInstance(material_homogenizationAt(el))) homogenization_typeInstance(material_homogenizationAt(el)))
case (HOMOGENIZATION_RGC_ID) chosenHomogenization case (HOMOGENIZATION_RGC_ID) chosenHomogenization
do c = 1, homogenization_Ngrains(material_homogenizationAt(el)) do c = 1, homogenization_Nconstituents(material_homogenizationAt(el))
dPdFs(:,:,:,:,c) = crystallite_stressTangent(c,ip,el) dPdFs(:,:,:,:,c) = crystallite_stressTangent(c,ip,el)
enddo enddo
call mech_RGC_averageStressAndItsTangent(& call mech_RGC_averageStressAndItsTangent(&
materialpoint_P(1:3,1:3,ip,el), & homogenization_P(1:3,1:3,ip,el), &
materialpoint_dPdF(1:3,1:3,1:3,1:3,ip,el),& homogenization_dPdF(1:3,1:3,1:3,1:3,ip,el),&
crystallite_P(1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el), & crystallite_P(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), &
dPdFs, & dPdFs, &
homogenization_typeInstance(material_homogenizationAt(el))) homogenization_typeInstance(material_homogenizationAt(el)))
end select chosenHomogenization end select chosenHomogenization
@ -539,10 +539,10 @@ subroutine homogenization_results
group = trim(group_base)//'/generic' group = trim(group_base)//'/generic'
call results_closeGroup(results_addGroup(group)) call results_closeGroup(results_addGroup(group))
!temp = reshape(materialpoint_F,[3,3,discretization_nIP*discretization_nElem]) !temp = reshape(homogenization_F,[3,3,discretization_nIPs*discretization_Nelems])
!call results_writeDataset(group,temp,'F',& !call results_writeDataset(group,temp,'F',&
! 'deformation gradient','1') ! 'deformation gradient','1')
!temp = reshape(materialpoint_P,[3,3,discretization_nIP*discretization_nElem]) !temp = reshape(homogenization_P,[3,3,discretization_nIPs*discretization_Nelems])
!call results_writeDataset(group,temp,'P',& !call results_writeDataset(group,temp,'P',&
! '1st Piola-Kirchhoff stress','Pa') ! '1st Piola-Kirchhoff stress','Pa')

38
src/homogenization_mech_RGC.f90
View File

@ -81,9 +81,9 @@ module subroutine mech_RGC_init(num_homogMech)
num_homogMech !< pointer to mechanical homogenization numerics data num_homogMech !< pointer to mechanical homogenization numerics data
integer :: & integer :: &
Ninstance, & Ninstances, &
h, & h, &
NofMyHomog, & Nmaterialpoints, &
sizeState, nIntFaceTot sizeState, nIntFaceTot
class (tNode), pointer :: & class (tNode), pointer :: &
@ -94,8 +94,8 @@ module subroutine mech_RGC_init(num_homogMech)
print'(/,a)', ' <<<+- homogenization_mech_rgc init -+>>>' print'(/,a)', ' <<<+- homogenization_mech_rgc init -+>>>'
Ninstance = count(homogenization_type == HOMOGENIZATION_RGC_ID) Ninstances = count(homogenization_type == HOMOGENIZATION_RGC_ID)
print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT) print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
print*, 'Tjahjanto et al., International Journal of Material Forming 2(1):939–942, 2009' print*, 'Tjahjanto et al., International Journal of Material Forming 2(1):939–942, 2009'
print*, 'https://doi.org/10.1007/s12289-009-0619-1'//IO_EOL print*, 'https://doi.org/10.1007/s12289-009-0619-1'//IO_EOL
@ -105,10 +105,10 @@ module subroutine mech_RGC_init(num_homogMech)
allocate(param(Ninstance)) allocate(param(Ninstances))
allocate(state(Ninstance)) allocate(state(Ninstances))
allocate(state0(Ninstance)) allocate(state0(Ninstances))
allocate(dependentState(Ninstance)) allocate(dependentState(Ninstances))
num_RGC => num_homogMech%get('RGC',defaultVal=emptyDict) num_RGC => num_homogMech%get('RGC',defaultVal=emptyDict)
@ -164,7 +164,7 @@ module subroutine mech_RGC_init(num_homogMech)
#endif #endif
prm%N_constituents = homogMech%get_asInts('cluster_size',requiredSize=3) prm%N_constituents = homogMech%get_asInts('cluster_size',requiredSize=3)
if (homogenization_Ngrains(h) /= product(prm%N_constituents)) & if (homogenization_Nconstituents(h) /= product(prm%N_constituents)) &
call IO_error(211,ext_msg='N_constituents (mech_RGC)') call IO_error(211,ext_msg='N_constituents (mech_RGC)')
prm%xi_alpha = homogMech%get_asFloat('xi_alpha') prm%xi_alpha = homogMech%get_asFloat('xi_alpha')
@ -173,7 +173,7 @@ module subroutine mech_RGC_init(num_homogMech)
prm%D_alpha = homogMech%get_asFloats('D_alpha', requiredSize=3) prm%D_alpha = homogMech%get_asFloats('D_alpha', requiredSize=3)
prm%a_g = homogMech%get_asFloats('a_g', requiredSize=3) prm%a_g = homogMech%get_asFloats('a_g', requiredSize=3)
NofMyHomog = count(material_homogenizationAt == h) Nmaterialpoints = count(material_homogenizationAt == h)
nIntFaceTot = 3*( (prm%N_constituents(1)-1)*prm%N_constituents(2)*prm%N_constituents(3) & nIntFaceTot = 3*( (prm%N_constituents(1)-1)*prm%N_constituents(2)*prm%N_constituents(3) &
+ prm%N_constituents(1)*(prm%N_constituents(2)-1)*prm%N_constituents(3) & + prm%N_constituents(1)*(prm%N_constituents(2)-1)*prm%N_constituents(3) &
+ prm%N_constituents(1)*prm%N_constituents(2)*(prm%N_constituents(3)-1)) + prm%N_constituents(1)*prm%N_constituents(2)*(prm%N_constituents(3)-1))
@ -181,24 +181,24 @@ module subroutine mech_RGC_init(num_homogMech)
+ size(['avg constitutive work ','average penalty energy']) + size(['avg constitutive work ','average penalty energy'])
homogState(h)%sizeState = sizeState homogState(h)%sizeState = sizeState
allocate(homogState(h)%state0 (sizeState,NofMyHomog), source=0.0_pReal) allocate(homogState(h)%state0 (sizeState,Nmaterialpoints), source=0.0_pReal)
allocate(homogState(h)%subState0(sizeState,NofMyHomog), source=0.0_pReal) allocate(homogState(h)%subState0(sizeState,Nmaterialpoints), source=0.0_pReal)
allocate(homogState(h)%state (sizeState,NofMyHomog), source=0.0_pReal) allocate(homogState(h)%state (sizeState,Nmaterialpoints), source=0.0_pReal)
stt%relaxationVector => homogState(h)%state(1:nIntFaceTot,:) stt%relaxationVector => homogState(h)%state(1:nIntFaceTot,:)
st0%relaxationVector => homogState(h)%state0(1:nIntFaceTot,:) st0%relaxationVector => homogState(h)%state0(1:nIntFaceTot,:)
stt%work => homogState(h)%state(nIntFaceTot+1,:) stt%work => homogState(h)%state(nIntFaceTot+1,:)
stt%penaltyEnergy => homogState(h)%state(nIntFaceTot+2,:) stt%penaltyEnergy => homogState(h)%state(nIntFaceTot+2,:)
allocate(dst%volumeDiscrepancy( NofMyHomog), source=0.0_pReal) allocate(dst%volumeDiscrepancy( Nmaterialpoints), source=0.0_pReal)
allocate(dst%relaxationRate_avg( NofMyHomog), source=0.0_pReal) allocate(dst%relaxationRate_avg( Nmaterialpoints), source=0.0_pReal)
allocate(dst%relaxationRate_max( NofMyHomog), source=0.0_pReal) allocate(dst%relaxationRate_max( Nmaterialpoints), source=0.0_pReal)
allocate(dst%mismatch( 3,NofMyHomog), source=0.0_pReal) allocate(dst%mismatch( 3,Nmaterialpoints), source=0.0_pReal)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! assigning cluster orientations ! assigning cluster orientations
dependentState(homogenization_typeInstance(h))%orientation = spread(eu2om(prm%a_g*inRad),3,NofMyHomog) dependentState(homogenization_typeInstance(h))%orientation = spread(eu2om(prm%a_g*inRad),3,Nmaterialpoints)
!dst%orientation = spread(eu2om(prm%a_g*inRad),3,NofMyHomog) ifort version 18.0.1 crashes (for whatever reason) !dst%orientation = spread(eu2om(prm%a_g*inRad),3,Nmaterialpoints) ifort version 18.0.1 crashes (for whatever reason)
end associate end associate

22
src/homogenization_mech_isostrain.f90
View File

@ -18,7 +18,7 @@ submodule(homogenization) homogenization_mech_isostrain
mapping mapping
end type end type
type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance) type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstances)
contains contains
@ -29,9 +29,9 @@ contains
module subroutine mech_isostrain_init module subroutine mech_isostrain_init
integer :: & integer :: &
Ninstance, & Ninstances, &
h, & h, &
NofMyHomog Nmaterialpoints
class(tNode), pointer :: & class(tNode), pointer :: &
material_homogenization, & material_homogenization, &
homog, & homog, &
@ -39,10 +39,10 @@ module subroutine mech_isostrain_init
print'(/,a)', ' <<<+- homogenization_mech_isostrain init -+>>>' print'(/,a)', ' <<<+- homogenization_mech_isostrain init -+>>>'
Ninstance = count(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID) Ninstances = count(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID)
print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT) print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
allocate(param(Ninstance)) ! one container of parameters per instance allocate(param(Ninstances)) ! one container of parameters per instance
material_homogenization => config_material%get('homogenization') material_homogenization => config_material%get('homogenization')
do h = 1, size(homogenization_type) do h = 1, size(homogenization_type)
@ -51,7 +51,7 @@ module subroutine mech_isostrain_init
homogMech => homog%get('mech') homogMech => homog%get('mech')
associate(prm => param(homogenization_typeInstance(h))) associate(prm => param(homogenization_typeInstance(h)))
prm%N_constituents = homogenization_Ngrains(h) prm%N_constituents = homogenization_Nconstituents(h)
select case(homogMech%get_asString('mapping',defaultVal = 'sum')) select case(homogMech%get_asString('mapping',defaultVal = 'sum'))
case ('sum') case ('sum')
prm%mapping = parallel_ID prm%mapping = parallel_ID
@ -61,11 +61,11 @@ module subroutine mech_isostrain_init
call IO_error(211,ext_msg='sum'//' (mech_isostrain)') call IO_error(211,ext_msg='sum'//' (mech_isostrain)')
end select end select
NofMyHomog = count(material_homogenizationAt == h) Nmaterialpoints = count(material_homogenizationAt == h)
homogState(h)%sizeState = 0 homogState(h)%sizeState = 0
allocate(homogState(h)%state0 (0,NofMyHomog)) allocate(homogState(h)%state0 (0,Nmaterialpoints))
allocate(homogState(h)%subState0(0,NofMyHomog)) allocate(homogState(h)%subState0(0,Nmaterialpoints))
allocate(homogState(h)%state (0,NofMyHomog)) allocate(homogState(h)%state (0,Nmaterialpoints))
end associate end associate

18
src/homogenization_mech_none.f90
View File

@ -14,26 +14,26 @@ contains
module subroutine mech_none_init module subroutine mech_none_init
integer :: & integer :: &
Ninstance, & Ninstances, &
h, & h, &
NofMyHomog Nmaterialpoints
print'(/,a)', ' <<<+- homogenization_mech_none init -+>>>' print'(/,a)', ' <<<+- homogenization_mech_none init -+>>>'
Ninstance = count(homogenization_type == HOMOGENIZATION_NONE_ID) Ninstances = count(homogenization_type == HOMOGENIZATION_NONE_ID)
print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT) print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
do h = 1, size(homogenization_type) do h = 1, size(homogenization_type)
if(homogenization_type(h) /= HOMOGENIZATION_NONE_ID) cycle if(homogenization_type(h) /= HOMOGENIZATION_NONE_ID) cycle
if(homogenization_Ngrains(h) /= 1) & if(homogenization_Nconstituents(h) /= 1) &
call IO_error(211,ext_msg='N_constituents (mech_none)') call IO_error(211,ext_msg='N_constituents (mech_none)')
NofMyHomog = count(material_homogenizationAt == h) Nmaterialpoints = count(material_homogenizationAt == h)
homogState(h)%sizeState = 0 homogState(h)%sizeState = 0
allocate(homogState(h)%state0 (0,NofMyHomog)) allocate(homogState(h)%state0 (0,Nmaterialpoints))
allocate(homogState(h)%subState0(0,NofMyHomog)) allocate(homogState(h)%subState0(0,Nmaterialpoints))
allocate(homogState(h)%state (0,NofMyHomog)) allocate(homogState(h)%state (0,Nmaterialpoints))
enddo enddo

12
src/kinematics_cleavage_opening.f90
View File

@ -20,7 +20,7 @@ submodule(constitutive:constitutive_damage) kinematics_cleavage_opening
cleavage_systems cleavage_systems
end type tParameters end type tParameters
type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance) type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstances)
contains contains
@ -35,7 +35,7 @@ module function kinematics_cleavage_opening_init(kinematics_length) result(myKin
integer, intent(in) :: kinematics_length integer, intent(in) :: kinematics_length
logical, dimension(:,:), allocatable :: myKinematics logical, dimension(:,:), allocatable :: myKinematics
integer :: Ninstance,p,k integer :: Ninstances,p,k
integer, dimension(:), allocatable :: N_cl !< active number of cleavage systems per family integer, dimension(:), allocatable :: N_cl !< active number of cleavage systems per family
character(len=pStringLen) :: extmsg = '' character(len=pStringLen) :: extmsg = ''
class(tNode), pointer :: & class(tNode), pointer :: &
@ -48,12 +48,12 @@ module function kinematics_cleavage_opening_init(kinematics_length) result(myKin
print'(/,a)', ' <<<+- kinematics_cleavage_opening init -+>>>' print'(/,a)', ' <<<+- kinematics_cleavage_opening init -+>>>'
myKinematics = kinematics_active('cleavage_opening',kinematics_length) myKinematics = kinematics_active('cleavage_opening',kinematics_length)
Ninstance = count(myKinematics) Ninstances = count(myKinematics)
print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT) print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
if(Ninstance == 0) return if(Ninstances == 0) return
phases => config_material%get('phase') phases => config_material%get('phase')
allocate(param(Ninstance)) allocate(param(Ninstances))
allocate(kinematics_cleavage_opening_instance(phases%length), source=0) allocate(kinematics_cleavage_opening_instance(phases%length), source=0)
do p = 1, phases%length do p = 1, phases%length

12
src/kinematics_slipplane_opening.f90
View File

@ -22,7 +22,7 @@ submodule(constitutive:constitutive_damage) kinematics_slipplane_opening
P_n P_n
end type tParameters end type tParameters
type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance) type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstances)
contains contains
@ -37,7 +37,7 @@ module function kinematics_slipplane_opening_init(kinematics_length) result(myKi
integer, intent(in) :: kinematics_length integer, intent(in) :: kinematics_length
logical, dimension(:,:), allocatable :: myKinematics logical, dimension(:,:), allocatable :: myKinematics
integer :: Ninstance,p,i,k integer :: Ninstances,p,i,k
character(len=pStringLen) :: extmsg = '' character(len=pStringLen) :: extmsg = ''
integer, dimension(:), allocatable :: N_sl integer, dimension(:), allocatable :: N_sl
real(pReal), dimension(:,:), allocatable :: d,n,t real(pReal), dimension(:,:), allocatable :: d,n,t
@ -51,13 +51,13 @@ module function kinematics_slipplane_opening_init(kinematics_length) result(myKi
print'(/,a)', ' <<<+- kinematics_slipplane init -+>>>' print'(/,a)', ' <<<+- kinematics_slipplane init -+>>>'
myKinematics = kinematics_active('slipplane_opening',kinematics_length) myKinematics = kinematics_active('slipplane_opening',kinematics_length)
Ninstance = count(myKinematics) Ninstances = count(myKinematics)
print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT) print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
if(Ninstance == 0) return if(Ninstances == 0) return
phases => config_material%get('phase') phases => config_material%get('phase')
allocate(kinematics_slipplane_opening_instance(phases%length), source=0) allocate(kinematics_slipplane_opening_instance(phases%length), source=0)
allocate(param(Ninstance)) allocate(param(Ninstances))
do p = 1, phases%length do p = 1, phases%length
if(any(myKinematics(:,p))) kinematics_slipplane_opening_instance(p) = count(myKinematics(:,1:p)) if(any(myKinematics(:,p))) kinematics_slipplane_opening_instance(p) = count(myKinematics(:,1:p))

10
src/kinematics_thermal_expansion.f90
View File

@ -29,7 +29,7 @@ module function kinematics_thermal_expansion_init(kinematics_length) result(myKi
integer, intent(in) :: kinematics_length integer, intent(in) :: kinematics_length
logical, dimension(:,:), allocatable :: myKinematics logical, dimension(:,:), allocatable :: myKinematics
integer :: Ninstance,p,i,k integer :: Ninstances,p,i,k
real(pReal), dimension(:), allocatable :: temp real(pReal), dimension(:), allocatable :: temp
class(tNode), pointer :: & class(tNode), pointer :: &
phases, & phases, &
@ -41,12 +41,12 @@ module function kinematics_thermal_expansion_init(kinematics_length) result(myKi
print'(/,a)', ' <<<+- kinematics_thermal_expansion init -+>>>' print'(/,a)', ' <<<+- kinematics_thermal_expansion init -+>>>'
myKinematics = kinematics_active('thermal_expansion',kinematics_length) myKinematics = kinematics_active('thermal_expansion',kinematics_length)
Ninstance = count(myKinematics) Ninstances = count(myKinematics)
print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT) print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
if(Ninstance == 0) return if(Ninstances == 0) return
phases => config_material%get('phase') phases => config_material%get('phase')
allocate(param(Ninstance)) allocate(param(Ninstances))
allocate(kinematics_thermal_expansion_instance(phases%length), source=0) allocate(kinematics_thermal_expansion_instance(phases%length), source=0)
do p = 1, phases%length do p = 1, phases%length

207
src/material.f90
View File

@ -52,7 +52,7 @@ module material
HOMOGENIZATION_RGC_ID HOMOGENIZATION_RGC_ID
end enum end enum
character(len=pStringLen), public, protected, allocatable, dimension(:) :: & character(len=pStringLen), public, protected, allocatable, dimension(:) :: &
material_name_phase, & !< name of each phase material_name_phase, & !< name of each phase
material_name_homogenization !< name of each homogenization material_name_homogenization !< name of each homogenization
@ -64,13 +64,10 @@ module material
homogenization_type !< type of each homogenization homogenization_type !< type of each homogenization
integer, public, protected :: & integer, public, protected :: &
material_Nhomogenization !< number of homogenizations homogenization_maxNconstituents !< max number of grains in any USED homogenization
integer, public, protected :: &
homogenization_maxNgrains !< max number of grains in any USED homogenization
integer, dimension(:), allocatable, public, protected :: & integer, dimension(:), allocatable, public, protected :: &
homogenization_Ngrains, & !< number of grains in each homogenization homogenization_Nconstituents, & !< number of grains in each homogenization
homogenization_typeInstance, & !< instance of particular type of each homogenization homogenization_typeInstance, & !< instance of particular type of each homogenization
thermal_typeInstance, & !< instance of particular type of each thermal transport thermal_typeInstance, & !< instance of particular type of each thermal transport
damage_typeInstance !< instance of particular type of each nonlocal damage damage_typeInstance !< instance of particular type of each nonlocal damage
@ -83,9 +80,9 @@ module material
material_homogenizationAt !< homogenization ID of each element material_homogenizationAt !< homogenization ID of each element
integer, dimension(:,:), allocatable, public, target :: & ! (ip,elem) ToDo: ugly target for mapping hack integer, dimension(:,:), allocatable, public, target :: & ! (ip,elem) ToDo: ugly target for mapping hack
material_homogenizationMemberAt !< position of the element within its homogenization instance material_homogenizationMemberAt !< position of the element within its homogenization instance
integer, dimension(:,:), allocatable, public, protected :: & ! (constituent,elem) integer, dimension(:,:), allocatable, public, protected :: & ! (constituent,elem)
material_phaseAt !< phase ID of each element material_phaseAt !< phase ID of each element
integer, dimension(:,:,:), allocatable, public, protected :: & ! (constituent,elem) integer, dimension(:,:,:), allocatable, public, protected :: & ! (constituent,IP,elem)
material_phaseMemberAt !< position of the element within its phase instance material_phaseMemberAt !< position of the element within its phase instance
type(tState), allocatable, dimension(:), public :: & type(tState), allocatable, dimension(:), public :: &
@ -96,11 +93,6 @@ module material
type(Rotation), dimension(:,:,:), allocatable, public, protected :: & type(Rotation), dimension(:,:,:), allocatable, public, protected :: &
material_orientation0 !< initial orientation of each grain,IP,element material_orientation0 !< initial orientation of each grain,IP,element
integer, dimension(:), allocatable, private :: &
material_Nconstituents !< number of constituents in each material
! BEGIN DEPRECATED ! BEGIN DEPRECATED
integer, dimension(:,:), allocatable, private, target :: mappingHomogenizationConst !< mapping from material points to offset in constant state/field integer, dimension(:,:), allocatable, private, target :: mappingHomogenizationConst !< mapping from material points to offset in constant state/field
! END DEPRECATED ! END DEPRECATED
@ -157,28 +149,10 @@ contains
subroutine material_init(restart) subroutine material_init(restart)
logical, intent(in) :: restart logical, intent(in) :: restart
integer :: myHomog
integer :: ph, myHomog
class(tNode), pointer :: &
phases, &
material_homogenization
character(len=pStringLen) :: sectionName
print'(/,a)', ' <<<+- material init -+>>>'; flush(IO_STDOUT) print'(/,a)', ' <<<+- material init -+>>>'; flush(IO_STDOUT)
phases => config_material%get('phase')
allocate(material_name_phase(phases%length))
do ph = 1, phases%length
write(sectionName,'(i0,a)') ph,'_'
material_name_phase(ph) = trim(adjustl(sectionName))//phases%getKey(ph) !ToDO: No reason to do. Update damage tests
enddo
material_homogenization => config_material%get('homogenization')
allocate(material_name_homogenization(material_homogenization%length))
do myHomog = 1, material_homogenization%length
write(sectionName,'(i0,a)') myHomog,'_'
material_name_homogenization(myHomog) = trim(adjustl(sectionName))//material_homogenization%getKey(myHomog)
enddo
call material_parseMaterial call material_parseMaterial
print*, 'Material parsed' print*, 'Material parsed'
@ -187,34 +161,32 @@ subroutine material_init(restart)
print*, 'Homogenization parsed' print*, 'Homogenization parsed'
if(homogenization_maxNgrains > size(material_phaseAt,1)) call IO_error(148) allocate(homogState (size(material_name_homogenization)))
allocate(thermalState (size(material_name_homogenization)))
allocate(damageState (size(material_name_homogenization)))
allocate(homogState (material_Nhomogenization)) allocate(thermalMapping (size(material_name_homogenization)))
allocate(thermalState (material_Nhomogenization)) allocate(damageMapping (size(material_name_homogenization)))
allocate(damageState (material_Nhomogenization))
allocate(thermalMapping (material_Nhomogenization)) allocate(temperature (size(material_name_homogenization)))
allocate(damageMapping (material_Nhomogenization)) allocate(damage (size(material_name_homogenization)))
allocate(temperature (material_Nhomogenization)) allocate(temperatureRate (size(material_name_homogenization)))
allocate(damage (material_Nhomogenization))
allocate(temperatureRate (material_Nhomogenization))
if (.not. restart) then if (.not. restart) then
call results_openJobFile call results_openJobFile
call results_mapping_constituent(material_phaseAt,material_phaseMemberAt,material_name_phase) call results_mapping_constituent(material_phaseAt,material_phaseMemberAt,material_name_phase)
call results_mapping_materialpoint(material_homogenizationAt,material_homogenizationMemberAt,material_name_homogenization) call results_mapping_homogenization(material_homogenizationAt,material_homogenizationMemberAt,material_name_homogenization)
call results_closeJobFile call results_closeJobFile
endif endif
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! BEGIN DEPRECATED ! BEGIN DEPRECATED
allocate(mappingHomogenizationConst( discretization_nIP,discretization_nElem),source=1) allocate(mappingHomogenizationConst( discretization_nIPs,discretization_Nelems),source=1)
! hack needed to initialize field values used during constitutive initialization ! hack needed to initialize field values used during constitutive initialization
do myHomog = 1,material_Nhomogenization do myHomog = 1, size(material_name_homogenization)
thermalMapping (myHomog)%p => mappingHomogenizationConst thermalMapping (myHomog)%p => mappingHomogenizationConst
damageMapping (myHomog)%p => mappingHomogenizationConst damageMapping (myHomog)%p => mappingHomogenizationConst
allocate(temperature (myHomog)%p(1), source=thermal_initialT(myHomog)) allocate(temperature (myHomog)%p(1), source=thermal_initialT(myHomog))
@ -225,6 +197,7 @@ subroutine material_init(restart)
end subroutine material_init end subroutine material_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief parses the homogenization part from the material configuration !> @brief parses the homogenization part from the material configuration
! ToDo: This should be done in homogenization ! ToDo: This should be done in homogenization
@ -241,22 +214,19 @@ subroutine material_parseHomogenization
integer :: h integer :: h
material_homogenization => config_material%get('homogenization') material_homogenization => config_material%get('homogenization')
material_Nhomogenization = material_homogenization%length
allocate(homogenization_type(material_Nhomogenization), source=HOMOGENIZATION_undefined_ID) allocate(homogenization_type(size(material_name_homogenization)), source=HOMOGENIZATION_undefined_ID)
allocate(thermal_type(material_Nhomogenization), source=THERMAL_isothermal_ID) allocate(thermal_type(size(material_name_homogenization)), source=THERMAL_isothermal_ID)
allocate(damage_type (material_Nhomogenization), source=DAMAGE_none_ID) allocate(damage_type (size(material_name_homogenization)), source=DAMAGE_none_ID)
allocate(homogenization_typeInstance(material_Nhomogenization), source=0) allocate(homogenization_typeInstance(size(material_name_homogenization)), source=0)
allocate(thermal_typeInstance(material_Nhomogenization), source=0) allocate(thermal_typeInstance(size(material_name_homogenization)), source=0)
allocate(damage_typeInstance(material_Nhomogenization), source=0) allocate(damage_typeInstance(size(material_name_homogenization)), source=0)
allocate(homogenization_Ngrains(material_Nhomogenization), source=0) allocate(thermal_initialT(size(material_name_homogenization)), source=300.0_pReal)
allocate(thermal_initialT(material_Nhomogenization), source=300.0_pReal) allocate(damage_initialPhi(size(material_name_homogenization)), source=1.0_pReal)
allocate(damage_initialPhi(material_Nhomogenization), source=1.0_pReal)
do h=1, material_Nhomogenization do h=1, size(material_name_homogenization)
homog => material_homogenization%get(h) homog => material_homogenization%get(h)
homogMech => homog%get('mech') homogMech => homog%get('mech')
homogenization_Ngrains(h) = homog%get_asInt('N_constituents')
select case (homogMech%get_asString('type')) select case (homogMech%get_asString('type'))
case('none') case('none')
homogenization_type(h) = HOMOGENIZATION_NONE_ID homogenization_type(h) = HOMOGENIZATION_NONE_ID
@ -302,15 +272,12 @@ subroutine material_parseHomogenization
endif endif
enddo enddo
do h=1, material_Nhomogenization do h=1, size(material_name_homogenization)
homogenization_typeInstance(h) = count(homogenization_type(1:h) == homogenization_type(h)) homogenization_typeInstance(h) = count(homogenization_type(1:h) == homogenization_type(h))
thermal_typeInstance(h) = count(thermal_type (1:h) == thermal_type (h)) thermal_typeInstance(h) = count(thermal_type (1:h) == thermal_type (h))
damage_typeInstance(h) = count(damage_type (1:h) == damage_type (h)) damage_typeInstance(h) = count(damage_type (1:h) == damage_type (h))
enddo enddo
homogenization_maxNgrains = maxval(homogenization_Ngrains)
end subroutine material_parseHomogenization end subroutine material_parseHomogenization
@ -324,7 +291,8 @@ subroutine material_parseMaterial
constituents, & !> list of constituents constituents, & !> list of constituents
constituent, & !> constituent definition constituent, & !> constituent definition
phases, & phases, &
homogenizations homogenizations, &
homogenization
integer, dimension(:), allocatable :: & integer, dimension(:), allocatable :: &
counterPhase, & counterPhase, &
@ -333,65 +301,108 @@ subroutine material_parseMaterial
real(pReal) :: & real(pReal) :: &
frac frac
integer :: & integer :: &
e, & e, i, c, &
i, & h
m, &
c, &
maxNconstituents
materials => config_material%get('material') materials => config_material%get('material')
if(any(discretization_materialAt > materials%length)) & phases => config_material%get('phase')
call IO_error(155,ext_msg='More materials requested than found in material.yaml')
allocate(material_Nconstituents(materials%length),source=0)
do m = 1, materials%length
material => materials%get(m)
constituents => material%get('constituents')
material_Nconstituents(m) = constituents%length
enddo
maxNconstituents = maxval(material_Nconstituents)
allocate(material_homogenizationAt(discretization_nElem),source=0)
allocate(material_homogenizationMemberAt(discretization_nIP,discretization_nElem),source=0)
allocate(material_phaseAt(maxNconstituents,discretization_nElem),source=0)
allocate(material_phaseMemberAt(maxNconstituents,discretization_nIP,discretization_nElem),source=0)
allocate(material_orientation0(maxNconstituents,discretization_nIP,discretization_nElem))
phases => config_material%get('phase')
allocate(counterPhase(phases%length),source=0)
homogenizations => config_material%get('homogenization') homogenizations => config_material%get('homogenization')
call sanityCheck(materials, homogenizations)
material_name_phase = getKeys(phases)
material_name_homogenization = getKeys(homogenizations)
allocate(homogenization_Nconstituents(homogenizations%length))
do h=1, homogenizations%length
homogenization => homogenizations%get(h)
homogenization_Nconstituents(h) = homogenization%get_asInt('N_constituents')
enddo
homogenization_maxNconstituents = maxval(homogenization_Nconstituents)
allocate(counterPhase(phases%length),source=0)
allocate(counterHomogenization(homogenizations%length),source=0) allocate(counterHomogenization(homogenizations%length),source=0)
do e = 1, discretization_nElem allocate(material_homogenizationAt(discretization_Nelems),source=0)
allocate(material_homogenizationMemberAt(discretization_nIPs,discretization_Nelems),source=0)
allocate(material_phaseAt(homogenization_maxNconstituents,discretization_Nelems),source=0)
allocate(material_phaseMemberAt(homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems),source=0)
allocate(material_orientation0(homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems))
do e = 1, discretization_Nelems
material => materials%get(discretization_materialAt(e)) material => materials%get(discretization_materialAt(e))
constituents => material%get('constituents') constituents => material%get('constituents')
material_homogenizationAt(e) = homogenizations%getIndex(material%get_asString('homogenization')) material_homogenizationAt(e) = homogenizations%getIndex(material%get_asString('homogenization'))
do i = 1, discretization_nIP do i = 1, discretization_nIPs
counterHomogenization(material_homogenizationAt(e)) = counterHomogenization(material_homogenizationAt(e)) + 1 counterHomogenization(material_homogenizationAt(e)) = counterHomogenization(material_homogenizationAt(e)) + 1
material_homogenizationMemberAt(i,e) = counterHomogenization(material_homogenizationAt(e)) material_homogenizationMemberAt(i,e) = counterHomogenization(material_homogenizationAt(e))
enddo enddo
frac = 0.0_pReal frac = 0.0_pReal
do c = 1, constituents%length do c = 1, constituents%length
constituent => constituents%get(c) constituent => constituents%get(c)
frac = frac + constituent%get_asFloat('fraction') frac = frac + constituent%get_asFloat('fraction')
material_phaseAt(c,e) = phases%getIndex(constituent%get_asString('phase')) material_phaseAt(c,e) = phases%getIndex(constituent%get_asString('phase'))
do i = 1, discretization_nIP do i = 1, discretization_nIPs
counterPhase(material_phaseAt(c,e)) = counterPhase(material_phaseAt(c,e)) + 1 counterPhase(material_phaseAt(c,e)) = counterPhase(material_phaseAt(c,e)) + 1
material_phaseMemberAt(c,i,e) = counterPhase(material_phaseAt(c,e)) material_phaseMemberAt(c,i,e) = counterPhase(material_phaseAt(c,e))
call material_orientation0(c,i,e)%fromQuaternion(constituent%get_asFloats('O',requiredSize=4)) call material_orientation0(c,i,e)%fromQuaternion(constituent%get_asFloats('O',requiredSize=4)) ! should be done in crystallite
enddo enddo
enddo enddo
if (dNeq(frac,1.0_pReal)) call IO_error(153,ext_msg='constituent') if (dNeq(frac,1.0_pReal)) call IO_error(153,ext_msg='constituent')
enddo enddo
end subroutine material_parseMaterial end subroutine material_parseMaterial
!--------------------------------------------------------------------------------------------------
!> @brief Check if material.yaml is consistent and contains sufficient # of materials
!--------------------------------------------------------------------------------------------------
subroutine sanityCheck(materials,homogenizations)
class(tNode), intent(in) :: materials, &
homogenizations
class(tNode), pointer :: material, &
homogenization, &
constituents
integer :: m
if(maxval(discretization_materialAt) > materials%length) &
call IO_error(155,ext_msg='More materials requested than found in material.yaml')
do m = 1, materials%length
material => materials%get(m)
constituents => material%get('constituents')
homogenization => homogenizations%get(material%get_asString('homogenization'))
if(constituents%length /= homogenization%get_asInt('N_constituents')) call IO_error(148)
enddo
end subroutine sanityCheck
!--------------------------------------------------------------------------------------------------
!> @brief Get all keys from a dictionary (currently with #_ prefix)
!--------------------------------------------------------------------------------------------------
function getKeys(dict)
class(tNode), intent(in) :: dict
character(len=pStringLen), dimension(:), allocatable :: getKeys
integer :: i
character(len=pStringLen) :: sectionName
allocate(getKeys(dict%length))
do i=1, dict%length
write(sectionName,'(i0,a)') i,'_'
getKeys(i) = trim(adjustl(sectionName))//dict%getKey(i) !ToDo: remove prefix
enddo
end function getKeys
end module material end module material

2
src/mesh/FEM_utilities.f90
View File

@ -164,7 +164,7 @@ subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData)
cutBack = .false. ! reset cutBack status cutBack = .false. ! reset cutBack status
P_av = sum(sum(materialpoint_P,dim=4),dim=3) * wgt ! average of P P_av = sum(sum(homogenization_P,dim=4),dim=3) * wgt ! average of P
call MPI_Allreduce(MPI_IN_PLACE,P_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,P_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
end subroutine utilities_constitutiveResponse end subroutine utilities_constitutiveResponse

20
src/mesh/mesh_mech_FEM.f90
View File

@ -400,15 +400,15 @@ subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr)
(((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp+1)*dimPlex)) (((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp+1)*dimPlex))
enddo enddo
enddo enddo
materialpoint_F(1:dimPlex,1:dimPlex,qPt+1,cell+1) = & homogenization_F(1:dimPlex,1:dimPlex,qPt+1,cell+1) = &
reshape(matmul(BMat,x_scal),shape=[dimPlex,dimPlex], order=[2,1]) reshape(matmul(BMat,x_scal),shape=[dimPlex,dimPlex], order=[2,1])
enddo enddo
if (num%BBarStabilisation) then if (num%BBarStabilisation) then
detFAvg = math_det33(sum(materialpoint_F(1:3,1:3,1:nQuadrature,cell+1),dim=3)/real(nQuadrature)) detFAvg = math_det33(sum(homogenization_F(1:3,1:3,1:nQuadrature,cell+1),dim=3)/real(nQuadrature))
do qPt = 1, nQuadrature do qPt = 1, nQuadrature
materialpoint_F(1:dimPlex,1:dimPlex,qPt,cell+1) = & homogenization_F(1:dimPlex,1:dimPlex,qPt,cell+1) = &
materialpoint_F(1:dimPlex,1:dimPlex,qPt,cell+1)* & homogenization_F(1:dimPlex,1:dimPlex,qPt,cell+1)* &
(detFAvg/math_det33(materialpoint_F(1:3,1:3,qPt,cell+1)))**(1.0/real(dimPlex)) (detFAvg/math_det33(homogenization_F(1:3,1:3,qPt,cell+1)))**(1.0/real(dimPlex))
enddo enddo
endif endif
@ -443,7 +443,7 @@ subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr)
enddo enddo
f_scal = f_scal + & f_scal = f_scal + &
matmul(transpose(BMat), & matmul(transpose(BMat), &
reshape(transpose(materialpoint_P(1:dimPlex,1:dimPlex,qPt+1,cell+1)), & reshape(transpose(homogenization_P(1:dimPlex,1:dimPlex,qPt+1,cell+1)), &
shape=[dimPlex*dimPlex]))*qWeights(qPt+1) shape=[dimPlex*dimPlex]))*qWeights(qPt+1)
enddo enddo
f_scal = f_scal*abs(detJ) f_scal = f_scal*abs(detJ)
@ -545,7 +545,7 @@ subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr)
(((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp+1)*dimPlex)) (((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp+1)*dimPlex))
enddo enddo
enddo enddo
MatA = matmul(reshape(reshape(materialpoint_dPdF(1:dimPlex,1:dimPlex,1:dimPlex,1:dimPlex,qPt+1,cell+1), & MatA = matmul(reshape(reshape(homogenization_dPdF(1:dimPlex,1:dimPlex,1:dimPlex,1:dimPlex,qPt+1,cell+1), &
shape=[dimPlex,dimPlex,dimPlex,dimPlex], order=[2,1,4,3]), & shape=[dimPlex,dimPlex,dimPlex,dimPlex], order=[2,1,4,3]), &
shape=[dimPlex*dimPlex,dimPlex*dimPlex]),BMat)*qWeights(qPt+1) shape=[dimPlex*dimPlex,dimPlex*dimPlex]),BMat)*qWeights(qPt+1)
if (num%BBarStabilisation) then if (num%BBarStabilisation) then
@ -553,12 +553,12 @@ subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr)
FInv = math_inv33(F) FInv = math_inv33(F)
K_eA = K_eA + matmul(transpose(BMat),MatA)*math_det33(FInv)**(1.0/real(dimPlex)) K_eA = K_eA + matmul(transpose(BMat),MatA)*math_det33(FInv)**(1.0/real(dimPlex))
K_eB = K_eB - & K_eB = K_eB - &
matmul(transpose(matmul(reshape(materialpoint_F(1:dimPlex,1:dimPlex,qPt+1,cell+1), & matmul(transpose(matmul(reshape(homogenization_F(1:dimPlex,1:dimPlex,qPt+1,cell+1), &
shape=[dimPlex*dimPlex,1]), & shape=[dimPlex*dimPlex,1]), &
matmul(reshape(FInv(1:dimPlex,1:dimPlex), & matmul(reshape(FInv(1:dimPlex,1:dimPlex), &
shape=[1,dimPlex*dimPlex],order=[2,1]),BMat))),MatA) shape=[1,dimPlex*dimPlex],order=[2,1]),BMat))),MatA)
MatB = MatB + & MatB = MatB + &
matmul(reshape(materialpoint_F(1:dimPlex,1:dimPlex,qPt+1,cell+1),shape=[1,dimPlex*dimPlex]),MatA) matmul(reshape(homogenization_F(1:dimPlex,1:dimPlex,qPt+1,cell+1),shape=[1,dimPlex*dimPlex]),MatA)
FAvg = FAvg + F FAvg = FAvg + F
BMatAvg = BMatAvg + BMat BMatAvg = BMatAvg + BMat
else else
@ -630,7 +630,7 @@ subroutine FEM_mech_forward(guess,timeinc,timeinc_old,fieldBC)
! forward last inc ! forward last inc
if (guess .and. .not. cutBack) then if (guess .and. .not. cutBack) then
ForwardData = .True. ForwardData = .True.
materialpoint_F0 = materialpoint_F homogenization_F0 = homogenization_F
call SNESGetDM(mech_snes,dm_local,ierr); CHKERRQ(ierr) !< retrieve mesh info from mech_snes into dm_local 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 DMGetSection(dm_local,section,ierr); CHKERRQ(ierr)
call DMGetLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr) call DMGetLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr)

24
src/results.f90
View File

@ -56,7 +56,7 @@ module results
results_addAttribute, & results_addAttribute, &
results_removeLink, & results_removeLink, &
results_mapping_constituent, & results_mapping_constituent, &
results_mapping_materialpoint results_mapping_homogenization
contains contains
subroutine results_init(restart) subroutine results_init(restart)
@ -644,7 +644,7 @@ end subroutine results_mapping_constituent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief adds the unique mapping from spatial position and constituent ID to results !> @brief adds the unique mapping from spatial position and constituent ID to results
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine results_mapping_materialpoint(homogenizationAt,memberAtLocal,label) subroutine results_mapping_homogenization(homogenizationAt,memberAtLocal,label)
integer, dimension(:), intent(in) :: homogenizationAt !< homogenization section at (element) integer, dimension(:), intent(in) :: homogenizationAt !< homogenization section at (element)
integer, dimension(:,:), intent(in) :: memberAtLocal !< homogenization member at (IP,element) integer, dimension(:,:), intent(in) :: memberAtLocal !< homogenization member at (IP,element)
@ -711,13 +711,13 @@ subroutine results_mapping_materialpoint(homogenizationAt,memberAtLocal,label)
! MPI settings and communication ! MPI settings and communication
#ifdef PETSc #ifdef PETSc
call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, ierr) call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, ierr)
if (ierr < 0) call IO_error(1,ext_msg='results_mapping_materialpoint: h5pset_dxpl_mpio_f') if (ierr < 0) call IO_error(1,ext_msg='results_mapping_homogenization: h5pset_dxpl_mpio_f')
call MPI_allreduce(MPI_IN_PLACE,writeSize,worldsize,MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr) ! get output at each process call MPI_allreduce(MPI_IN_PLACE,writeSize,worldsize,MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr) ! get output at each process
if (ierr /= 0) call IO_error(894,ext_msg='results_mapping_materialpoint: MPI_allreduce/writeSize') if (ierr /= 0) call IO_error(894,ext_msg='results_mapping_homogenization: MPI_allreduce/writeSize')
call MPI_allreduce(MPI_IN_PLACE,memberOffset,size(memberOffset),MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr)! get offset at each process call MPI_allreduce(MPI_IN_PLACE,memberOffset,size(memberOffset),MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr)! get offset at each process
if (ierr /= 0) call IO_error(894,ext_msg='results_mapping_materialpoint: MPI_allreduce/memberOffset') if (ierr /= 0) call IO_error(894,ext_msg='results_mapping_homogenization: MPI_allreduce/memberOffset')
#endif #endif
myShape = int([writeSize(worldrank)], HSIZE_T) myShape = int([writeSize(worldrank)], HSIZE_T)
@ -727,13 +727,13 @@ subroutine results_mapping_materialpoint(homogenizationAt,memberAtLocal,label)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! create dataspace in memory (local shape = hyperslab) and in file (global shape) ! create dataspace in memory (local shape = hyperslab) and in file (global shape)
call h5screate_simple_f(1,myShape,memspace_id,ierr,myShape) call h5screate_simple_f(1,myShape,memspace_id,ierr,myShape)
if (ierr < 0) call IO_error(1,ext_msg='results_mapping_materialpoint: h5screate_simple_f/memspace_id') if (ierr < 0) call IO_error(1,ext_msg='results_mapping_homogenization: h5screate_simple_f/memspace_id')
call h5screate_simple_f(1,totalShape,filespace_id,ierr,totalShape) call h5screate_simple_f(1,totalShape,filespace_id,ierr,totalShape)
if (ierr < 0) call IO_error(1,ext_msg='results_mapping_materialpoint: h5screate_simple_f/filespace_id') if (ierr < 0) call IO_error(1,ext_msg='results_mapping_homogenization: h5screate_simple_f/filespace_id')
call h5sselect_hyperslab_f(filespace_id, H5S_SELECT_SET_F, myOffset, myShape, ierr) call h5sselect_hyperslab_f(filespace_id, H5S_SELECT_SET_F, myOffset, myShape, ierr)
if (ierr < 0) call IO_error(1,ext_msg='results_mapping_materialpoint: h5sselect_hyperslab_f') if (ierr < 0) call IO_error(1,ext_msg='results_mapping_homogenization: h5sselect_hyperslab_f')
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
! expand phaseAt to consider IPs (is not stored per IP) ! expand phaseAt to consider IPs (is not stored per IP)
@ -753,14 +753,14 @@ subroutine results_mapping_materialpoint(homogenizationAt,memberAtLocal,label)
loc_id = results_openGroup('/mapping') loc_id = results_openGroup('/mapping')
call h5dcreate_f(loc_id, 'homogenization', dtype_id, filespace_id, dset_id, ierr) call h5dcreate_f(loc_id, 'homogenization', dtype_id, filespace_id, dset_id, ierr)
if (ierr < 0) call IO_error(1,ext_msg='results_mapping_materialpoint: h5dcreate_f') if (ierr < 0) call IO_error(1,ext_msg='results_mapping_homogenization: h5dcreate_f')
call h5dwrite_f(dset_id, name_id, reshape(label(pack(homogenizationAtMaterialpoint,.true.)),myShape), & call h5dwrite_f(dset_id, name_id, reshape(label(pack(homogenizationAtMaterialpoint,.true.)),myShape), &
myShape, ierr, file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id) myShape, ierr, file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id)
if (ierr < 0) call IO_error(1,ext_msg='results_mapping_materialpoint: h5dwrite_f/name_id') if (ierr < 0) call IO_error(1,ext_msg='results_mapping_homogenization: h5dwrite_f/name_id')
call h5dwrite_f(dset_id, position_id, reshape(pack(memberAtGlobal,.true.),myShape), & call h5dwrite_f(dset_id, position_id, reshape(pack(memberAtGlobal,.true.),myShape), &
myShape, ierr, file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id) myShape, ierr, file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id)
if (ierr < 0) call IO_error(1,ext_msg='results_mapping_materialpoint: h5dwrite_f/position_id') if (ierr < 0) call IO_error(1,ext_msg='results_mapping_homogenization: h5dwrite_f/position_id')
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! close all ! close all
@ -776,7 +776,7 @@ subroutine results_mapping_materialpoint(homogenizationAt,memberAtLocal,label)
! for backward compatibility ! for backward compatibility
call results_setLink('/mapping/homogenization','/mapping/cellResults/materialpoint') call results_setLink('/mapping/homogenization','/mapping/cellResults/materialpoint')
end subroutine results_mapping_materialpoint end subroutine results_mapping_homogenization
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------

16
src/source_damage_anisoBrittle.f90
View File

@ -25,7 +25,7 @@ submodule (constitutive:constitutive_damage) source_damage_anisoBrittle
output output
end type tParameters end type tParameters
type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance) type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstances)
contains contains
@ -45,19 +45,19 @@ module function source_damage_anisoBrittle_init(source_length) result(mySources)
phase, & phase, &
sources, & sources, &
src src
integer :: Ninstance,sourceOffset,NipcMyPhase,p integer :: Ninstances,sourceOffset,Nconstituents,p
integer, dimension(:), allocatable :: N_cl integer, dimension(:), allocatable :: N_cl
character(len=pStringLen) :: extmsg = '' character(len=pStringLen) :: extmsg = ''
print'(/,a)', ' <<<+- source_damage_anisoBrittle init -+>>>' print'(/,a)', ' <<<+- source_damage_anisoBrittle init -+>>>'
mySources = source_active('damage_anisoBrittle',source_length) mySources = source_active('damage_anisoBrittle',source_length)
Ninstance = count(mySources) Ninstances = count(mySources)
print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT) print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
if(Ninstance == 0) return if(Ninstances == 0) return
phases => config_material%get('phase') phases => config_material%get('phase')
allocate(param(Ninstance)) allocate(param(Ninstances))
allocate(source_damage_anisoBrittle_offset (phases%length), source=0) allocate(source_damage_anisoBrittle_offset (phases%length), source=0)
allocate(source_damage_anisoBrittle_instance(phases%length), source=0) allocate(source_damage_anisoBrittle_instance(phases%length), source=0)
@ -100,8 +100,8 @@ module function source_damage_anisoBrittle_init(source_length) result(mySources)
if (any(prm%g_crit < 0.0_pReal)) extmsg = trim(extmsg)//' g_crit' if (any(prm%g_crit < 0.0_pReal)) extmsg = trim(extmsg)//' g_crit'
if (any(prm%s_crit < 0.0_pReal)) extmsg = trim(extmsg)//' s_crit' if (any(prm%s_crit < 0.0_pReal)) extmsg = trim(extmsg)//' s_crit'
NipcMyPhase = count(material_phaseAt==p) * discretization_nIP Nconstituents = count(material_phaseAt==p) * discretization_nIPs
call constitutive_allocateState(sourceState(p)%p(sourceOffset),NipcMyPhase,1,1,0) call constitutive_allocateState(sourceState(p)%p(sourceOffset),Nconstituents,1,1,0)
sourceState(p)%p(sourceOffset)%atol = src%get_asFloat('anisobrittle_atol',defaultVal=1.0e-3_pReal) sourceState(p)%p(sourceOffset)%atol = src%get_asFloat('anisobrittle_atol',defaultVal=1.0e-3_pReal)
if(any(sourceState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_atol' if(any(sourceState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_atol'

16
src/source_damage_anisoDuctile.f90
View File

@ -19,7 +19,7 @@ submodule(constitutive:constitutive_damage) source_damage_anisoDuctile
output output
end type tParameters end type tParameters
type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance) type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstances)
contains contains
@ -39,19 +39,19 @@ module function source_damage_anisoDuctile_init(source_length) result(mySources)
pl, & pl, &
sources, & sources, &
src src
integer :: Ninstance,sourceOffset,NipcMyPhase,p integer :: Ninstances,sourceOffset,Nconstituents,p
integer, dimension(:), allocatable :: N_sl integer, dimension(:), allocatable :: N_sl
character(len=pStringLen) :: extmsg = '' character(len=pStringLen) :: extmsg = ''
print'(/,a)', ' <<<+- source_damage_anisoDuctile init -+>>>' print'(/,a)', ' <<<+- source_damage_anisoDuctile init -+>>>'
mySources = source_active('damage_anisoDuctile',source_length) mySources = source_active('damage_anisoDuctile',source_length)
Ninstance = count(mySources) Ninstances = count(mySources)
print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT) print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
if(Ninstance == 0) return if(Ninstances == 0) return
phases => config_material%get('phase') phases => config_material%get('phase')
allocate(param(Ninstance)) allocate(param(Ninstances))
allocate(source_damage_anisoDuctile_offset (phases%length), source=0) allocate(source_damage_anisoDuctile_offset (phases%length), source=0)
allocate(source_damage_anisoDuctile_instance(phases%length), source=0) allocate(source_damage_anisoDuctile_instance(phases%length), source=0)
@ -84,8 +84,8 @@ module function source_damage_anisoDuctile_init(source_length) result(mySources)
if (prm%q <= 0.0_pReal) extmsg = trim(extmsg)//' q' if (prm%q <= 0.0_pReal) extmsg = trim(extmsg)//' q'
if (any(prm%gamma_crit < 0.0_pReal)) extmsg = trim(extmsg)//' gamma_crit' if (any(prm%gamma_crit < 0.0_pReal)) extmsg = trim(extmsg)//' gamma_crit'
NipcMyPhase=count(material_phaseAt==p) * discretization_nIP Nconstituents=count(material_phaseAt==p) * discretization_nIPs
call constitutive_allocateState(sourceState(p)%p(sourceOffset),NipcMyPhase,1,1,0) call constitutive_allocateState(sourceState(p)%p(sourceOffset),Nconstituents,1,1,0)
sourceState(p)%p(sourceOffset)%atol = src%get_asFloat('anisoDuctile_atol',defaultVal=1.0e-3_pReal) sourceState(p)%p(sourceOffset)%atol = src%get_asFloat('anisoDuctile_atol',defaultVal=1.0e-3_pReal)
if(any(sourceState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisoductile_atol' if(any(sourceState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisoductile_atol'

16
src/source_damage_isoBrittle.f90
View File

@ -17,7 +17,7 @@ submodule(constitutive:constitutive_damage) source_damage_isoBrittle
output output
end type tParameters end type tParameters
type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance) type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstances)
contains contains
@ -36,18 +36,18 @@ module function source_damage_isoBrittle_init(source_length) result(mySources)
phase, & phase, &
sources, & sources, &
src src
integer :: Ninstance,sourceOffset,NipcMyPhase,p integer :: Ninstances,sourceOffset,Nconstituents,p
character(len=pStringLen) :: extmsg = '' character(len=pStringLen) :: extmsg = ''
print'(/,a)', ' <<<+- source_damage_isoBrittle init -+>>>' print'(/,a)', ' <<<+- source_damage_isoBrittle init -+>>>'
mySources = source_active('damage_isoBrittle',source_length) mySources = source_active('damage_isoBrittle',source_length)
Ninstance = count(mySources) Ninstances = count(mySources)
print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT) print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
if(Ninstance == 0) return if(Ninstances == 0) return
phases => config_material%get('phase') phases => config_material%get('phase')
allocate(param(Ninstance)) allocate(param(Ninstances))
allocate(source_damage_isoBrittle_offset (phases%length), source=0) allocate(source_damage_isoBrittle_offset (phases%length), source=0)
allocate(source_damage_isoBrittle_instance(phases%length), source=0) allocate(source_damage_isoBrittle_instance(phases%length), source=0)
@ -73,8 +73,8 @@ module function source_damage_isoBrittle_init(source_length) result(mySources)
! sanity checks ! sanity checks
if (prm%W_crit <= 0.0_pReal) extmsg = trim(extmsg)//' W_crit' if (prm%W_crit <= 0.0_pReal) extmsg = trim(extmsg)//' W_crit'
NipcMyPhase = count(material_phaseAt==p) * discretization_nIP Nconstituents = count(material_phaseAt==p) * discretization_nIPs
call constitutive_allocateState(sourceState(p)%p(sourceOffset),NipcMyPhase,1,1,1) call constitutive_allocateState(sourceState(p)%p(sourceOffset),Nconstituents,1,1,1)
sourceState(p)%p(sourceOffset)%atol = src%get_asFloat('isoBrittle_atol',defaultVal=1.0e-3_pReal) sourceState(p)%p(sourceOffset)%atol = src%get_asFloat('isoBrittle_atol',defaultVal=1.0e-3_pReal)
if(any(sourceState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isobrittle_atol' if(any(sourceState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isobrittle_atol'

16
src/source_damage_isoDuctile.f90
View File

@ -18,7 +18,7 @@ submodule (constitutive:constitutive_damage) source_damage_isoDuctile
output output
end type tParameters end type tParameters
type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance) type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstances)
contains contains
@ -38,18 +38,18 @@ module function source_damage_isoDuctile_init(source_length) result(mySources)
phase, & phase, &
sources, & sources, &
src src
integer :: Ninstance,sourceOffset,NipcMyPhase,p integer :: Ninstances,sourceOffset,Nconstituents,p
character(len=pStringLen) :: extmsg = '' character(len=pStringLen) :: extmsg = ''
print'(/,a)', ' <<<+- source_damage_isoDuctile init -+>>>' print'(/,a)', ' <<<+- source_damage_isoDuctile init -+>>>'
mySources = source_active('damage_isoDuctile',source_length) mySources = source_active('damage_isoDuctile',source_length)
Ninstance = count(mySources) Ninstances = count(mySources)
print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT) print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
if(Ninstance == 0) return if(Ninstances == 0) return
phases => config_material%get('phase') phases => config_material%get('phase')
allocate(param(Ninstance)) allocate(param(Ninstances))
allocate(source_damage_isoDuctile_offset (phases%length), source=0) allocate(source_damage_isoDuctile_offset (phases%length), source=0)
allocate(source_damage_isoDuctile_instance(phases%length), source=0) allocate(source_damage_isoDuctile_instance(phases%length), source=0)
@ -77,8 +77,8 @@ module function source_damage_isoDuctile_init(source_length) result(mySources)
if (prm%q <= 0.0_pReal) extmsg = trim(extmsg)//' q' if (prm%q <= 0.0_pReal) extmsg = trim(extmsg)//' q'
if (prm%gamma_crit <= 0.0_pReal) extmsg = trim(extmsg)//' gamma_crit' if (prm%gamma_crit <= 0.0_pReal) extmsg = trim(extmsg)//' gamma_crit'
NipcMyPhase=count(material_phaseAt==p) * discretization_nIP Nconstituents=count(material_phaseAt==p) * discretization_nIPs
call constitutive_allocateState(sourceState(p)%p(sourceOffset),NipcMyPhase,1,1,0) call constitutive_allocateState(sourceState(p)%p(sourceOffset),Nconstituents,1,1,0)
sourceState(p)%p(sourceOffset)%atol = src%get_asFloat('isoDuctile_atol',defaultVal=1.0e-3_pReal) sourceState(p)%p(sourceOffset)%atol = src%get_asFloat('isoDuctile_atol',defaultVal=1.0e-3_pReal)
if(any(sourceState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isoductile_atol' if(any(sourceState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isoductile_atol'

18
src/source_thermal_dissipation.f90
View File

@ -15,7 +15,7 @@ submodule(constitutive:constitutive_thermal) source_thermal_dissipation
kappa !< TAYLOR-QUINNEY factor kappa !< TAYLOR-QUINNEY factor
end type tParameters end type tParameters
type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance) type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstances)
contains contains
@ -35,17 +35,17 @@ module function source_thermal_dissipation_init(source_length) result(mySources)
phase, & phase, &
sources, & sources, &
src src
integer :: Ninstance,sourceOffset,NipcMyPhase,p integer :: Ninstances,sourceOffset,Nconstituents,p
print'(/,a)', ' <<<+- source_thermal_dissipation init -+>>>' print'(/,a)', ' <<<+- source_thermal_dissipation init -+>>>'
mySources = source_active('thermal_dissipation',source_length) mySources = source_active('thermal_dissipation',source_length)
Ninstance = count(mySources) Ninstances = count(mySources)
print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT) print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
if(Ninstance == 0) return if(Ninstances == 0) return
phases => config_material%get('phase') phases => config_material%get('phase')
allocate(param(Ninstance)) allocate(param(Ninstances))
allocate(source_thermal_dissipation_offset (phases%length), source=0) allocate(source_thermal_dissipation_offset (phases%length), source=0)
allocate(source_thermal_dissipation_instance(phases%length), source=0) allocate(source_thermal_dissipation_instance(phases%length), source=0)
@ -61,8 +61,8 @@ module function source_thermal_dissipation_init(source_length) result(mySources)
src => sources%get(sourceOffset) src => sources%get(sourceOffset)
prm%kappa = src%get_asFloat('kappa') prm%kappa = src%get_asFloat('kappa')
NipcMyPhase = count(material_phaseAt==p) * discretization_nIP Nconstituents = count(material_phaseAt==p) * discretization_nIPs
call constitutive_allocateState(sourceState(p)%p(sourceOffset),NipcMyPhase,0,0,0) call constitutive_allocateState(sourceState(p)%p(sourceOffset),Nconstituents,0,0,0)
end associate end associate
endif endif
@ -74,7 +74,7 @@ end function source_thermal_dissipation_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Ninstances dissipation rate !> @brief Ninstancess dissipation rate
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDot_dT, Tstar, Lp, phase) module subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDot_dT, Tstar, Lp, phase)

16
src/source_thermal_externalheat.f90
View File

@ -19,7 +19,7 @@ submodule(constitutive:constitutive_thermal) source_thermal_externalheat
nIntervals nIntervals
end type tParameters end type tParameters
type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance) type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstances)
contains contains
@ -39,17 +39,17 @@ module function source_thermal_externalheat_init(source_length) result(mySources
phase, & phase, &
sources, & sources, &
src src
integer :: Ninstance,sourceOffset,NipcMyPhase,p integer :: Ninstances,sourceOffset,Nconstituents,p
print'(/,a)', ' <<<+- source_thermal_externalHeat init -+>>>' print'(/,a)', ' <<<+- source_thermal_externalHeat init -+>>>'
mySources = source_active('thermal_externalheat',source_length) mySources = source_active('thermal_externalheat',source_length)
Ninstance = count(mySources) Ninstances = count(mySources)
print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT) print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
if(Ninstance == 0) return if(Ninstances == 0) return
phases => config_material%get('phase') phases => config_material%get('phase')
allocate(param(Ninstance)) allocate(param(Ninstances))
allocate(source_thermal_externalheat_offset (phases%length), source=0) allocate(source_thermal_externalheat_offset (phases%length), source=0)
allocate(source_thermal_externalheat_instance(phases%length), source=0) allocate(source_thermal_externalheat_instance(phases%length), source=0)
@ -69,8 +69,8 @@ module function source_thermal_externalheat_init(source_length) result(mySources
prm%f_T = src%get_asFloats('f_T',requiredSize = size(prm%t_n)) prm%f_T = src%get_asFloats('f_T',requiredSize = size(prm%t_n))
NipcMyPhase = count(material_phaseAt==p) * discretization_nIP Nconstituents = count(material_phaseAt==p) * discretization_nIPs
call constitutive_allocateState(sourceState(p)%p(sourceOffset),NipcMyPhase,1,1,0) call constitutive_allocateState(sourceState(p)%p(sourceOffset),Nconstituents,1,1,0)
end associate end associate
endif endif

32
src/thermal_adiabatic.f90
View File

@ -40,7 +40,7 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine thermal_adiabatic_init subroutine thermal_adiabatic_init
integer :: maxNinstance,h,NofMyHomog integer :: maxNinstances,h,Nmaterialpoints
class(tNode), pointer :: & class(tNode), pointer :: &
material_homogenization, & material_homogenization, &
homog, & homog, &
@ -48,13 +48,13 @@ subroutine thermal_adiabatic_init
print'(/,a)', ' <<<+- thermal_adiabatic init -+>>>'; flush(6) print'(/,a)', ' <<<+- thermal_adiabatic init -+>>>'; flush(6)
maxNinstance = count(thermal_type == THERMAL_adiabatic_ID) maxNinstances = count(thermal_type == THERMAL_adiabatic_ID)
if (maxNinstance == 0) return if (maxNinstances == 0) return
allocate(param(maxNinstance)) allocate(param(maxNinstances))
material_homogenization => config_material%get('homogenization') material_homogenization => config_material%get('homogenization')
do h = 1, material_Nhomogenization do h = 1, size(material_name_homogenization)
if (thermal_type(h) /= THERMAL_adiabatic_ID) cycle if (thermal_type(h) /= THERMAL_adiabatic_ID) cycle
homog => material_homogenization%get(h) homog => material_homogenization%get(h)
homogThermal => homog%get('thermal') homogThermal => homog%get('thermal')
@ -67,17 +67,17 @@ subroutine thermal_adiabatic_init
prm%output = homogThermal%get_asStrings('output',defaultVal=emptyStringArray) prm%output = homogThermal%get_asStrings('output',defaultVal=emptyStringArray)
#endif #endif
NofMyHomog=count(material_homogenizationAt==h) Nmaterialpoints=count(material_homogenizationAt==h)
thermalState(h)%sizeState = 1 thermalState(h)%sizeState = 1
allocate(thermalState(h)%state0 (1,NofMyHomog), source=thermal_initialT(h)) allocate(thermalState(h)%state0 (1,Nmaterialpoints), source=thermal_initialT(h))
allocate(thermalState(h)%subState0(1,NofMyHomog), source=thermal_initialT(h)) allocate(thermalState(h)%subState0(1,Nmaterialpoints), source=thermal_initialT(h))
allocate(thermalState(h)%state (1,NofMyHomog), source=thermal_initialT(h)) allocate(thermalState(h)%state (1,Nmaterialpoints), source=thermal_initialT(h))
thermalMapping(h)%p => material_homogenizationMemberAt thermalMapping(h)%p => material_homogenizationMemberAt
deallocate(temperature(h)%p) deallocate(temperature(h)%p)
temperature(h)%p => thermalState(h)%state(1,:) temperature(h)%p => thermalState(h)%state(1,:)
deallocate(temperatureRate(h)%p) deallocate(temperatureRate(h)%p)
allocate (temperatureRate(h)%p(NofMyHomog), source=0.0_pReal) allocate (temperatureRate(h)%p(Nmaterialpoints), source=0.0_pReal)
end associate end associate
enddo enddo
@ -145,8 +145,8 @@ subroutine thermal_adiabatic_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el)
homog = material_homogenizationAt(el) homog = material_homogenizationAt(el)
call constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, crystallite_S, crystallite_Lp, ip, el) call constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, crystallite_S, crystallite_Lp, ip, el)
Tdot = Tdot/real(homogenization_Ngrains(homog),pReal) Tdot = Tdot/real(homogenization_Nconstituents(homog),pReal)
dTdot_dT = dTdot_dT/real(homogenization_Ngrains(homog),pReal) dTdot_dT = dTdot_dT/real(homogenization_Nconstituents(homog),pReal)
end subroutine thermal_adiabatic_getSourceAndItsTangent end subroutine thermal_adiabatic_getSourceAndItsTangent
@ -167,13 +167,13 @@ function thermal_adiabatic_getSpecificHeat(ip,el)
thermal_adiabatic_getSpecificHeat = 0.0_pReal thermal_adiabatic_getSpecificHeat = 0.0_pReal
do grain = 1, homogenization_Ngrains(material_homogenizationAt(el)) do grain = 1, homogenization_Nconstituents(material_homogenizationAt(el))
thermal_adiabatic_getSpecificHeat = thermal_adiabatic_getSpecificHeat & thermal_adiabatic_getSpecificHeat = thermal_adiabatic_getSpecificHeat &
+ lattice_c_p(material_phaseAt(grain,el)) + lattice_c_p(material_phaseAt(grain,el))
enddo enddo
thermal_adiabatic_getSpecificHeat = thermal_adiabatic_getSpecificHeat & thermal_adiabatic_getSpecificHeat = thermal_adiabatic_getSpecificHeat &
/ real(homogenization_Ngrains(material_homogenizationAt(el)),pReal) / real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal)
end function thermal_adiabatic_getSpecificHeat end function thermal_adiabatic_getSpecificHeat
@ -193,13 +193,13 @@ function thermal_adiabatic_getMassDensity(ip,el)
thermal_adiabatic_getMassDensity = 0.0_pReal thermal_adiabatic_getMassDensity = 0.0_pReal
do grain = 1, homogenization_Ngrains(material_homogenizationAt(el)) do grain = 1, homogenization_Nconstituents(material_homogenizationAt(el))
thermal_adiabatic_getMassDensity = thermal_adiabatic_getMassDensity & thermal_adiabatic_getMassDensity = thermal_adiabatic_getMassDensity &
+ lattice_rho(material_phaseAt(grain,el)) + lattice_rho(material_phaseAt(grain,el))
enddo enddo
thermal_adiabatic_getMassDensity = thermal_adiabatic_getMassDensity & thermal_adiabatic_getMassDensity = thermal_adiabatic_getMassDensity &
/ real(homogenization_Ngrains(material_homogenizationAt(el)),pReal) / real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal)
end function thermal_adiabatic_getMassDensity end function thermal_adiabatic_getMassDensity

36
src/thermal_conduction.f90
View File

@ -41,7 +41,7 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine thermal_conduction_init subroutine thermal_conduction_init
integer :: Ninstance,NofMyHomog,h integer :: Ninstances,Nmaterialpoints,h
class(tNode), pointer :: & class(tNode), pointer :: &
material_homogenization, & material_homogenization, &
homog, & homog, &
@ -49,11 +49,11 @@ subroutine thermal_conduction_init
print'(/,a)', ' <<<+- thermal_conduction init -+>>>'; flush(6) print'(/,a)', ' <<<+- thermal_conduction init -+>>>'; flush(6)
Ninstance = count(thermal_type == THERMAL_conduction_ID) Ninstances = count(thermal_type == THERMAL_conduction_ID)
allocate(param(Ninstance)) allocate(param(Ninstances))
material_homogenization => config_material%get('homogenization') material_homogenization => config_material%get('homogenization')
do h = 1, material_Nhomogenization do h = 1, size(material_name_homogenization)
if (thermal_type(h) /= THERMAL_conduction_ID) cycle if (thermal_type(h) /= THERMAL_conduction_ID) cycle
homog => material_homogenization%get(h) homog => material_homogenization%get(h)
homogThermal => homog%get('thermal') homogThermal => homog%get('thermal')
@ -65,17 +65,17 @@ subroutine thermal_conduction_init
prm%output = homogThermal%get_asStrings('output',defaultVal=emptyStringArray) prm%output = homogThermal%get_asStrings('output',defaultVal=emptyStringArray)
#endif #endif
NofMyHomog=count(material_homogenizationAt==h) Nmaterialpoints=count(material_homogenizationAt==h)
thermalState(h)%sizeState = 0 thermalState(h)%sizeState = 0
allocate(thermalState(h)%state0 (0,NofMyHomog)) allocate(thermalState(h)%state0 (0,Nmaterialpoints))
allocate(thermalState(h)%subState0(0,NofMyHomog)) allocate(thermalState(h)%subState0(0,Nmaterialpoints))
allocate(thermalState(h)%state (0,NofMyHomog)) allocate(thermalState(h)%state (0,Nmaterialpoints))
thermalMapping(h)%p => material_homogenizationMemberAt thermalMapping(h)%p => material_homogenizationMemberAt
deallocate(temperature (h)%p) deallocate(temperature (h)%p)
allocate (temperature (h)%p(NofMyHomog), source=thermal_initialT(h)) allocate (temperature (h)%p(Nmaterialpoints), source=thermal_initialT(h))
deallocate(temperatureRate(h)%p) deallocate(temperatureRate(h)%p)
allocate (temperatureRate(h)%p(NofMyHomog), source=0.0_pReal) allocate (temperatureRate(h)%p(Nmaterialpoints), source=0.0_pReal)
end associate end associate
enddo enddo
@ -104,8 +104,8 @@ subroutine thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el)
homog = material_homogenizationAt(el) homog = material_homogenizationAt(el)
call constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, crystallite_S,crystallite_Lp ,ip, el) call constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, crystallite_S,crystallite_Lp ,ip, el)
Tdot = Tdot/real(homogenization_Ngrains(homog),pReal) Tdot = Tdot/real(homogenization_Nconstituents(homog),pReal)
dTdot_dT = dTdot_dT/real(homogenization_Ngrains(homog),pReal) dTdot_dT = dTdot_dT/real(homogenization_Nconstituents(homog),pReal)
end subroutine thermal_conduction_getSourceAndItsTangent end subroutine thermal_conduction_getSourceAndItsTangent
@ -125,13 +125,13 @@ function thermal_conduction_getConductivity(ip,el)
thermal_conduction_getConductivity = 0.0_pReal thermal_conduction_getConductivity = 0.0_pReal
do grain = 1, homogenization_Ngrains(material_homogenizationAt(el)) do grain = 1, homogenization_Nconstituents(material_homogenizationAt(el))
thermal_conduction_getConductivity = thermal_conduction_getConductivity + & thermal_conduction_getConductivity = thermal_conduction_getConductivity + &
crystallite_push33ToRef(grain,ip,el,lattice_K(:,:,material_phaseAt(grain,el))) crystallite_push33ToRef(grain,ip,el,lattice_K(:,:,material_phaseAt(grain,el)))
enddo enddo
thermal_conduction_getConductivity = thermal_conduction_getConductivity & thermal_conduction_getConductivity = thermal_conduction_getConductivity &
/ real(homogenization_Ngrains(material_homogenizationAt(el)),pReal) / real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal)
end function thermal_conduction_getConductivity end function thermal_conduction_getConductivity
@ -151,13 +151,13 @@ function thermal_conduction_getSpecificHeat(ip,el)
thermal_conduction_getSpecificHeat = 0.0_pReal thermal_conduction_getSpecificHeat = 0.0_pReal
do grain = 1, homogenization_Ngrains(material_homogenizationAt(el)) do grain = 1, homogenization_Nconstituents(material_homogenizationAt(el))
thermal_conduction_getSpecificHeat = thermal_conduction_getSpecificHeat & thermal_conduction_getSpecificHeat = thermal_conduction_getSpecificHeat &
+ lattice_c_p(material_phaseAt(grain,el)) + lattice_c_p(material_phaseAt(grain,el))
enddo enddo
thermal_conduction_getSpecificHeat = thermal_conduction_getSpecificHeat & thermal_conduction_getSpecificHeat = thermal_conduction_getSpecificHeat &
/ real(homogenization_Ngrains(material_homogenizationAt(el)),pReal) / real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal)
end function thermal_conduction_getSpecificHeat end function thermal_conduction_getSpecificHeat
@ -178,13 +178,13 @@ function thermal_conduction_getMassDensity(ip,el)
thermal_conduction_getMassDensity = 0.0_pReal thermal_conduction_getMassDensity = 0.0_pReal
do grain = 1, homogenization_Ngrains(material_homogenizationAt(el)) do grain = 1, homogenization_Nconstituents(material_homogenizationAt(el))
thermal_conduction_getMassDensity = thermal_conduction_getMassDensity & thermal_conduction_getMassDensity = thermal_conduction_getMassDensity &
+ lattice_rho(material_phaseAt(grain,el)) + lattice_rho(material_phaseAt(grain,el))
enddo enddo
thermal_conduction_getMassDensity = thermal_conduction_getMassDensity & thermal_conduction_getMassDensity = thermal_conduction_getMassDensity &
/ real(homogenization_Ngrains(material_homogenizationAt(el)),pReal) / real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal)
end function thermal_conduction_getMassDensity end function thermal_conduction_getMassDensity

12
src/thermal_isothermal.f90
View File

@ -16,18 +16,18 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine thermal_isothermal_init subroutine thermal_isothermal_init
integer :: h,NofMyHomog integer :: h,Nmaterialpoints
print'(/,a)', ' <<<+- thermal_isothermal init -+>>>'; flush(6) print'(/,a)', ' <<<+- thermal_isothermal init -+>>>'; flush(6)
do h = 1, material_Nhomogenization do h = 1, size(material_name_homogenization)
if (thermal_type(h) /= THERMAL_isothermal_ID) cycle if (thermal_type(h) /= THERMAL_isothermal_ID) cycle
NofMyHomog = count(material_homogenizationAt == h) Nmaterialpoints = count(material_homogenizationAt == h)
thermalState(h)%sizeState = 0 thermalState(h)%sizeState = 0
allocate(thermalState(h)%state0 (0,NofMyHomog)) allocate(thermalState(h)%state0 (0,Nmaterialpoints))
allocate(thermalState(h)%subState0(0,NofMyHomog)) allocate(thermalState(h)%subState0(0,Nmaterialpoints))
allocate(thermalState(h)%state (0,NofMyHomog)) allocate(thermalState(h)%state (0,Nmaterialpoints))
deallocate(temperature (h)%p) deallocate(temperature (h)%p)
allocate (temperature (h)%p(1), source=thermal_initialT(h)) allocate (temperature (h)%p(1), source=thermal_initialT(h))