diff --git a/PRIVATE b/PRIVATE index 76bb51348..bc6de828c 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 76bb51348de75207d483d369628670e5ae51dca9 +Subproject commit bc6de828cc4ee9c941b37113ca49fcf51abd3512 diff --git a/examples/config/phase/mechanical/plastic/phenopowerlaw_Mg.yaml b/examples/config/phase/mechanical/plastic/phenopowerlaw_Mg.yaml index b29075173..9b06f92d9 100644 --- a/examples/config/phase/mechanical/plastic/phenopowerlaw_Mg.yaml +++ b/examples/config/phase/mechanical/plastic/phenopowerlaw_Mg.yaml @@ -6,12 +6,12 @@ references: output: [xi_sl, xi_tw] -N_sl: [3, 3, 0, 6, 0, 6] # basal, 1. prism, -, 1. pyr, -, 2. pyr -N_tw: [6, 0, 6] # tension, -, compression +N_sl: [3, 3, 0, 6, 0, 6] # basal, prism, -, 1. pyr, -, 2. pyr +N_tw: [6, 0, 6] # tension, -, compression -xi_0_sl: [10.e+6, 55.e+6, 0., 60.e+6, 0., 60.e+6] +xi_0_sl: [10.e+6, 55.e+6, 0., 60.e+6, 0., 60.e+6] xi_inf_sl: [40.e+6, 135.e+6, 0., 150.e+6, 0., 150.e+6] -xi_0_tw: [40.e+6, 0., 60.e+6] +xi_0_tw: [40.e+6, 0., 60.e+6] a_sl: 2.25 dot_gamma_0_sl: 0.001 @@ -21,9 +21,18 @@ n_tw: 20 f_sat_sl-tw: 10.0 h_0_sl-sl: 500.0e+6 -h_0_tw-tw: 50.0e+6 +h_0_tw-tw: 50.0e+6 h_0_tw-sl: 150.0e+6 -h_sl-sl: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1] -h_tw-tw: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1] -h_tw-sl: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1] -h_sl-tw: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1] +h_sl-sl: [+1.0, 1.0, 1.0, 1.0, 1.0, 1.0, -1.0, -1.0, -1.0, -1.0, + -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, + -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, + +1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, + +1.0, 1.0] # unused entries are indicated by -1.0 +h_tw-tw: [+1.0, 1.0, -1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, + -1.0, 1.0] # unused entries are indicated by -1.0 +h_tw-sl: [+1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, -1.0, -1.0, + -1.0, -1.0, 1.0, -1.0, 1.0, -1.0, -1.0, -1.0, -1.0, -1.0, + +1.0, -1.0, 1.0, -1.0] # unused entries are indicated by -1.0 +h_sl-tw: [+1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, -1.0, -1.0, + -1.0, -1.0, 1.0, -1.0, 1.0, -1.0, -1.0, -1.0, -1.0, -1.0, + +1.0, -1.0, 1.0] # unused entries are indicated by -1.0 diff --git a/examples/config/phase/mechanical/plastic/phenopowerlaw_Ti.yaml b/examples/config/phase/mechanical/plastic/phenopowerlaw_Ti.yaml index 804cf4541..85fd2ee8e 100644 --- a/examples/config/phase/mechanical/plastic/phenopowerlaw_Ti.yaml +++ b/examples/config/phase/mechanical/plastic/phenopowerlaw_Ti.yaml @@ -8,7 +8,7 @@ references: https://doi.org/10.1016/j.actamat.2017.05.015 output: [gamma_sl] -N_sl: [3, 3, 0, 0, 12] # basal, 1. prism, -, -, 2. pyr +N_sl: [3, 3, 0, 0, 12] # basal, 1. prism, -, -, 1. pyr n_sl: 20 a_sl: 2.0 dot_gamma_0_sl: 0.001 @@ -20,4 +20,6 @@ xi_inf_sl: [568.e+6, 150.e+7, 0.0, 0.0, 3420.e+6] # L. Wang et al. : # xi_0_sl: [127.e+6, 96.e+6, 0.0, 0.0, 240.e+6] -h_sl-sl: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1] +h_sl-sl: [+1.0, 1.0, 1.0, 1.0, 1.0, 1.0, -1.0, -1.0, -1.0, -1.0, + -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, + +1.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0, 1.0] # unused entries are indicated by -1.0 diff --git a/python/damask/VERSION b/python/damask/VERSION index 3b4b0236b..cbf66449a 100644 --- a/python/damask/VERSION +++ b/python/damask/VERSION @@ -1 +1 @@ -v3.0.0-alpha5-135-ga3f74994b +v3.0.0-alpha5-191-gf32a78861 diff --git a/python/damask/_colormap.py b/python/damask/_colormap.py index 8c897e6ab..3060e99e0 100644 --- a/python/damask/_colormap.py +++ b/python/damask/_colormap.py @@ -2,6 +2,8 @@ import os import json import functools import colorsys +from pathlib import Path +from typing import Sequence, Union, TextIO import numpy as np import matplotlib as mpl @@ -14,9 +16,9 @@ from PIL import Image from . import util from . import Table -_eps = 216./24389. -_kappa = 24389./27. -_ref_white = np.array([.95047, 1.00000, 1.08883]) # Observer = 2, Illuminant = D65 +_EPS = 216./24389. +_KAPPA = 24389./27. +_REF_WHITE = np.array([.95047, 1.00000, 1.08883]) # Observer = 2, Illuminant = D65 # ToDo (if needed) # - support alpha channel (paraview/ASCII/input) @@ -39,20 +41,20 @@ class Colormap(mpl.colors.ListedColormap): """ - def __add__(self,other): + def __add__(self, other: "Colormap") -> "Colormap": """Concatenate.""" return Colormap(np.vstack((self.colors,other.colors)), f'{self.name}+{other.name}') - def __iadd__(self,other): + def __iadd__(self, other: "Colormap") -> "Colormap": """Concatenate (in-place).""" return self.__add__(other) - def __invert__(self): + def __invert__(self) -> "Colormap": """Reverse.""" return self.reversed() - def __repr__(self): + def __repr__(self) -> str: """Show as matplotlib figure.""" fig = plt.figure(self.name,figsize=(5,.5)) ax1 = fig.add_axes([0, 0, 1, 1]) @@ -64,7 +66,11 @@ class Colormap(mpl.colors.ListedColormap): @staticmethod - def from_range(low,high,name='DAMASK colormap',N=256,model='rgb'): + def from_range(low: Sequence[float], + high: Sequence[float], + name: str = 'DAMASK colormap', + N: int = 256, + model: str = 'rgb') -> "Colormap": """ Create a perceptually uniform colormap between given (inclusive) bounds. @@ -145,7 +151,7 @@ class Colormap(mpl.colors.ListedColormap): @staticmethod - def from_predefined(name,N=256): + def from_predefined(name: str, N: int = 256) -> "Colormap": """ Select from a set of predefined colormaps. @@ -185,7 +191,10 @@ class Colormap(mpl.colors.ListedColormap): return Colormap.from_range(definition['low'],definition['high'],name,N) - def shade(self,field,bounds=None,gap=None): + def shade(self, + field: np.ndarray, + bounds: Sequence[float] = None, + gap: float = None) -> Image: """ Generate PIL image of 2D field using colormap. @@ -226,7 +235,7 @@ class Colormap(mpl.colors.ListedColormap): mode='RGBA') - def reversed(self,name=None): + def reversed(self, name: str = None) -> "Colormap": """ Reverse. @@ -251,7 +260,7 @@ class Colormap(mpl.colors.ListedColormap): return Colormap(np.array(rev.colors),rev.name[:-4] if rev.name.endswith('_r_r') else rev.name) - def _get_file_handle(self,fname,extension): + def _get_file_handle(self, fname: Union[TextIO, str, Path, None], suffix: str) -> TextIO: """ Provide file handle. @@ -259,8 +268,7 @@ class Colormap(mpl.colors.ListedColormap): ---------- fname : file, str, pathlib.Path, or None Filename or filehandle, will be name of the colormap+extension if None. - - extension: str + suffix: str Extension of the filename. Returns @@ -270,17 +278,14 @@ class Colormap(mpl.colors.ListedColormap): """ if fname is None: - fhandle = open(self.name.replace(' ','_')+'.'+extension,'w',newline='\n') + return open(self.name.replace(' ','_')+suffix, 'w', newline='\n') + elif isinstance(fname, (str, Path)): + return open(fname, 'w', newline='\n') else: - try: - fhandle = open(fname,'w',newline='\n') - except TypeError: - fhandle = fname - - return fhandle + return fname - def save_paraview(self,fname=None): + def save_paraview(self, fname: Union[TextIO, str, Path] = None): """ Save as JSON file for use in Paraview. @@ -303,10 +308,10 @@ class Colormap(mpl.colors.ListedColormap): 'RGBPoints':colors }] - json.dump(out,self._get_file_handle(fname,'json'),indent=4) + json.dump(out,self._get_file_handle(fname,'.json'),indent=4) - def save_ASCII(self,fname=None): + def save_ASCII(self, fname: Union[TextIO, str, Path] = None): """ Save as ASCII file. @@ -319,10 +324,10 @@ class Colormap(mpl.colors.ListedColormap): """ labels = {'RGBA':4} if self.colors.shape[1] == 4 else {'RGB': 3} t = Table(self.colors,labels,f'Creator: {util.execution_stamp("Colormap")}') - t.save(self._get_file_handle(fname,'txt')) + t.save(self._get_file_handle(fname,'.txt')) - def save_GOM(self,fname=None): + def save_GOM(self, fname: Union[TextIO, str, Path] = None): """ Save as ASCII file for use in GOM Aramis. @@ -340,10 +345,10 @@ class Colormap(mpl.colors.ListedColormap): + ' '.join([f' 0 {c[0]} {c[1]} {c[2]} 255 1' for c in reversed((self.colors*255).astype(int))]) \ + '\n' - self._get_file_handle(fname,'legend').write(GOM_str) + self._get_file_handle(fname,'.legend').write(GOM_str) - def save_gmsh(self,fname=None): + def save_gmsh(self, fname: Union[TextIO, str, Path] = None): """ Save as ASCII file for use in gmsh. @@ -358,11 +363,13 @@ class Colormap(mpl.colors.ListedColormap): gmsh_str = 'View.ColorTable = {\n' \ +'\n'.join([f'{c[0]},{c[1]},{c[2]},' for c in self.colors[:,:3]*255]) \ +'\n}\n' - self._get_file_handle(fname,'msh').write(gmsh_str) + self._get_file_handle(fname,'.msh').write(gmsh_str) @staticmethod - def _interpolate_msh(frac,low,high): + def _interpolate_msh(frac, + low: np.ndarray, + high: np.ndarray) -> np.ndarray: """ Interpolate in Msh color space. @@ -439,31 +446,31 @@ class Colormap(mpl.colors.ListedColormap): @staticmethod - def _hsv2rgb(hsv): + def _hsv2rgb(hsv: np.ndarray) -> np.ndarray: """H(ue) S(aturation) V(alue) to R(red) G(reen) B(lue).""" return np.array(colorsys.hsv_to_rgb(hsv[0]/360.,hsv[1],hsv[2])) @staticmethod - def _rgb2hsv(rgb): + def _rgb2hsv(rgb: np.ndarray) -> np.ndarray: """R(ed) G(reen) B(lue) to H(ue) S(aturation) V(alue).""" h,s,v = colorsys.rgb_to_hsv(rgb[0],rgb[1],rgb[2]) return np.array([h*360,s,v]) @staticmethod - def _hsl2rgb(hsl): + def _hsl2rgb(hsl: np.ndarray) -> np.ndarray: """H(ue) S(aturation) L(uminance) to R(red) G(reen) B(lue).""" return np.array(colorsys.hls_to_rgb(hsl[0]/360.,hsl[2],hsl[1])) @staticmethod - def _rgb2hsl(rgb): + def _rgb2hsl(rgb: np.ndarray) -> np.ndarray: """R(ed) G(reen) B(lue) to H(ue) S(aturation) L(uminance).""" h,l,s = colorsys.rgb_to_hls(rgb[0],rgb[1],rgb[2]) return np.array([h*360,s,l]) @staticmethod - def _xyz2rgb(xyz): + def _xyz2rgb(xyz: np.ndarray) -> np.ndarray: """ CIE Xyz to R(ed) G(reen) B(lue). @@ -483,7 +490,7 @@ class Colormap(mpl.colors.ListedColormap): return np.clip(rgb,0.,1.) @staticmethod - def _rgb2xyz(rgb): + def _rgb2xyz(rgb: np.ndarray) -> np.ndarray: """ R(ed) G(reen) B(lue) to CIE Xyz. @@ -501,7 +508,7 @@ class Colormap(mpl.colors.ListedColormap): @staticmethod - def _lab2xyz(lab,ref_white=None): + def _lab2xyz(lab: np.ndarray, ref_white: np.ndarray = None) -> np.ndarray: """ CIE Lab to CIE Xyz. @@ -514,13 +521,13 @@ class Colormap(mpl.colors.ListedColormap): f_z = (lab[0]+16.)/116. - lab[2]/200. return np.array([ - f_x**3. if f_x**3. > _eps else (116.*f_x-16.)/_kappa, - ((lab[0]+16.)/116.)**3 if lab[0]>_kappa*_eps else lab[0]/_kappa, - f_z**3. if f_z**3. > _eps else (116.*f_z-16.)/_kappa - ])*(ref_white if ref_white is not None else _ref_white) + f_x**3. if f_x**3. > _EPS else (116.*f_x-16.)/_KAPPA, + ((lab[0]+16.)/116.)**3 if lab[0]>_KAPPA*_EPS else lab[0]/_KAPPA, + f_z**3. if f_z**3. > _EPS else (116.*f_z-16.)/_KAPPA + ])*(ref_white if ref_white is not None else _REF_WHITE) @staticmethod - def _xyz2lab(xyz,ref_white=None): + def _xyz2lab(xyz: np.ndarray, ref_white: np.ndarray = None) -> np.ndarray: """ CIE Xyz to CIE Lab. @@ -529,8 +536,8 @@ class Colormap(mpl.colors.ListedColormap): http://www.brucelindbloom.com/index.html?Eqn_Lab_to_XYZ.html """ - ref_white = ref_white if ref_white is not None else _ref_white - f = np.where(xyz/ref_white > _eps,(xyz/ref_white)**(1./3.),(_kappa*xyz/ref_white+16.)/116.) + ref_white = ref_white if ref_white is not None else _REF_WHITE + f = np.where(xyz/ref_white > _EPS,(xyz/ref_white)**(1./3.),(_KAPPA*xyz/ref_white+16.)/116.) return np.array([ 116.0 * f[1] - 16.0, @@ -540,7 +547,7 @@ class Colormap(mpl.colors.ListedColormap): @staticmethod - def _lab2msh(lab): + def _lab2msh(lab: np.ndarray) -> np.ndarray: """ CIE Lab to Msh. @@ -558,7 +565,7 @@ class Colormap(mpl.colors.ListedColormap): ]) @staticmethod - def _msh2lab(msh): + def _msh2lab(msh: np.ndarray) -> np.ndarray: """ Msh to CIE Lab. @@ -575,29 +582,29 @@ class Colormap(mpl.colors.ListedColormap): ]) @staticmethod - def _lab2rgb(lab): + def _lab2rgb(lab: np.ndarray) -> np.ndarray: return Colormap._xyz2rgb(Colormap._lab2xyz(lab)) @staticmethod - def _rgb2lab(rgb): + def _rgb2lab(rgb: np.ndarray) -> np.ndarray: return Colormap._xyz2lab(Colormap._rgb2xyz(rgb)) @staticmethod - def _msh2rgb(msh): + def _msh2rgb(msh: np.ndarray) -> np.ndarray: return Colormap._lab2rgb(Colormap._msh2lab(msh)) @staticmethod - def _rgb2msh(rgb): + def _rgb2msh(rgb: np.ndarray) -> np.ndarray: return Colormap._lab2msh(Colormap._rgb2lab(rgb)) @staticmethod - def _hsv2msh(hsv): + def _hsv2msh(hsv: np.ndarray) -> np.ndarray: return Colormap._rgb2msh(Colormap._hsv2rgb(hsv)) @staticmethod - def _hsl2msh(hsl): + def _hsl2msh(hsl: np.ndarray) -> np.ndarray: return Colormap._rgb2msh(Colormap._hsl2rgb(hsl)) @staticmethod - def _xyz2msh(xyz): + def _xyz2msh(xyz: np.ndarray) -> np.ndarray: return Colormap._lab2msh(Colormap._xyz2lab(xyz)) diff --git a/python/tests/test_Colormap.py b/python/tests/test_Colormap.py index 81bf8d2f6..342165134 100644 --- a/python/tests/test_Colormap.py +++ b/python/tests/test_Colormap.py @@ -77,6 +77,12 @@ class TestColormap: # xyz2msh assert np.allclose(Colormap._xyz2msh(xyz),msh,atol=1.e-6,rtol=0) + @pytest.mark.parametrize('low,high',[((0,0,0),(1,1,1)), + ([0,0,0],[1,1,1])]) + def test_from_range_types(self,low,high): + a = Colormap.from_range(low,high) + b = Colormap.from_range(np.array(low),np.array(high)) + assert np.all(a.colors == b.colors) @pytest.mark.parametrize('format',['ASCII','paraview','GOM','gmsh']) @pytest.mark.parametrize('model',['rgb','hsv','hsl','xyz','lab','msh']) diff --git a/src/commercialFEM_fileList.f90 b/src/commercialFEM_fileList.f90 index e1d53ca83..e67149dea 100644 --- a/src/commercialFEM_fileList.f90 +++ b/src/commercialFEM_fileList.f90 @@ -4,6 +4,7 @@ !> @details List of files needed by MSC.Marc !-------------------------------------------------------------------------------------------------- #include "parallelization.f90" +#include "constants.f90" #include "IO.f90" #include "YAML_types.f90" #include "YAML_parse.f90" diff --git a/src/constants.f90 b/src/constants.f90 new file mode 100644 index 000000000..dd26ce78c --- /dev/null +++ b/src/constants.f90 @@ -0,0 +1,15 @@ +!-------------------------------------------------------------------------------------------------- +!> @author Martin Diehl, KU Leuven +!> @brief physical constants +!-------------------------------------------------------------------------------------------------- +module constants + use prec + + implicit none + public + + real(pReal), parameter :: & + T_ROOM = 300.0_pReal, & !< Room temperature in K + K_B = 1.38e-23_pReal !< Boltzmann constant in J/Kelvin + +end module constants diff --git a/src/lattice.f90 b/src/lattice.f90 index d90215c7e..77a1beff8 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -219,18 +219,18 @@ module lattice 2, -1, -1, 0, 0, 1, -1, 0, & -1, 2, -1, 0, -1, 0, 1, 0, & -1, -1, 2, 0, 1, -1, 0, 0, & - ! <-11.0>{11.0}/2nd order prismatic compound systems (plane normal independent of c/a-ratio) + ! <-11.0>{11.0}/2. order prismatic compound systems (plane normal independent of c/a-ratio) -1, 1, 0, 0, 1, 1, -2, 0, & 0, -1, 1, 0, -2, 1, 1, 0, & 1, 0, -1, 0, 1, -2, 1, 0, & - ! <-1-1.0>{-11.1}/1st order pyramidal systems (direction independent of c/a-ratio) + ! <-1-1.0>{-11.1}/1. order pyramidal systems (direction independent of c/a-ratio) -1, 2, -1, 0, 1, 0, -1, 1, & -2, 1, 1, 0, 0, 1, -1, 1, & -1, -1, 2, 0, -1, 1, 0, 1, & 1, -2, 1, 0, -1, 0, 1, 1, & 2, -1, -1, 0, 0, -1, 1, 1, & 1, 1, -2, 0, 1, -1, 0, 1, & - ! <11.3>{-10.1}/1st order pyramidal systems (direction independent of c/a-ratio) + ! <11.3>{-10.1}/1. order pyramidal systems (direction independent of c/a-ratio) -2, 1, 1, 3, 1, 0, -1, 1, & -1, -1, 2, 3, 1, 0, -1, 1, & -1, -1, 2, 3, 0, 1, -1, 1, & @@ -243,7 +243,7 @@ module lattice -1, 2, -1, 3, 0, -1, 1, 1, & -1, 2, -1, 3, 1, -1, 0, 1, & -2, 1, 1, 3, 1, -1, 0, 1, & - ! <11.3>{-1-1.2}/2nd order pyramidal systems + ! <11.3>{-1-1.2}/2. order pyramidal systems -1, -1, 2, 3, 1, 1, -2, 2, & 1, -2, 1, 3, -1, 2, -1, 2, & 2, -1, -1, 3, -2, 1, 1, 2, & @@ -491,7 +491,7 @@ end function lattice_characteristicShear_Twin !-------------------------------------------------------------------------------------------------- -!> @brief Rotated elasticity matrices for twinning in 66-vector notation +!> @brief Rotated elasticity matrices for twinning in 6x6-matrix notation !-------------------------------------------------------------------------------------------------- function lattice_C66_twin(Ntwin,C66,lattice,CoverA) @@ -529,7 +529,7 @@ end function lattice_C66_twin !-------------------------------------------------------------------------------------------------- -!> @brief Rotated elasticity matrices for transformation in 66-vector notation +!> @brief Rotated elasticity matrices for transformation in 6x6-matrix notation !-------------------------------------------------------------------------------------------------- function lattice_C66_trans(Ntrans,C_parent66,lattice_target, & cOverA_trans,a_bcc,a_fcc) @@ -753,22 +753,23 @@ function lattice_interaction_SlipBySlip(Nslip,interactionValues,lattice) result( integer, dimension(HEX_NSLIP,HEX_NSLIP), parameter :: & HEX_INTERACTIONSLIPSLIP = reshape( [& + ! basal prism 2. prism 1. pyr 1. pyr 2. pyr 1, 2, 2, 3, 3, 3, 7, 7, 7, 13,13,13,13,13,13, 21,21,21,21,21,21,21,21,21,21,21,21, 31,31,31,31,31,31, & ! -----> acting (forest) - 2, 1, 2, 3, 3, 3, 7, 7, 7, 13,13,13,13,13,13, 21,21,21,21,21,21,21,21,21,21,21,21, 31,31,31,31,31,31, & ! | + 2, 1, 2, 3, 3, 3, 7, 7, 7, 13,13,13,13,13,13, 21,21,21,21,21,21,21,21,21,21,21,21, 31,31,31,31,31,31, & ! | basal 2, 2, 1, 3, 3, 3, 7, 7, 7, 13,13,13,13,13,13, 21,21,21,21,21,21,21,21,21,21,21,21, 31,31,31,31,31,31, & ! | ! v 6, 6, 6, 4, 5, 5, 8, 8, 8, 14,14,14,14,14,14, 22,22,22,22,22,22,22,22,22,22,22,22, 32,32,32,32,32,32, & ! reacting (primary) - 6, 6, 6, 5, 4, 5, 8, 8, 8, 14,14,14,14,14,14, 22,22,22,22,22,22,22,22,22,22,22,22, 32,32,32,32,32,32, & + 6, 6, 6, 5, 4, 5, 8, 8, 8, 14,14,14,14,14,14, 22,22,22,22,22,22,22,22,22,22,22,22, 32,32,32,32,32,32, & ! prism 6, 6, 6, 5, 5, 4, 8, 8, 8, 14,14,14,14,14,14, 22,22,22,22,22,22,22,22,22,22,22,22, 32,32,32,32,32,32, & 12,12,12, 11,11,11, 9,10,10, 15,15,15,15,15,15, 23,23,23,23,23,23,23,23,23,23,23,23, 33,33,33,33,33,33, & - 12,12,12, 11,11,11, 10, 9,10, 15,15,15,15,15,15, 23,23,23,23,23,23,23,23,23,23,23,23, 33,33,33,33,33,33, & + 12,12,12, 11,11,11, 10, 9,10, 15,15,15,15,15,15, 23,23,23,23,23,23,23,23,23,23,23,23, 33,33,33,33,33,33, & ! 2. prism 12,12,12, 11,11,11, 10,10, 9, 15,15,15,15,15,15, 23,23,23,23,23,23,23,23,23,23,23,23, 33,33,33,33,33,33, & 20,20,20, 19,19,19, 18,18,18, 16,17,17,17,17,17, 24,24,24,24,24,24,24,24,24,24,24,24, 34,34,34,34,34,34, & 20,20,20, 19,19,19, 18,18,18, 17,16,17,17,17,17, 24,24,24,24,24,24,24,24,24,24,24,24, 34,34,34,34,34,34, & 20,20,20, 19,19,19, 18,18,18, 17,17,16,17,17,17, 24,24,24,24,24,24,24,24,24,24,24,24, 34,34,34,34,34,34, & - 20,20,20, 19,19,19, 18,18,18, 17,17,17,16,17,17, 24,24,24,24,24,24,24,24,24,24,24,24, 34,34,34,34,34,34, & + 20,20,20, 19,19,19, 18,18,18, 17,17,17,16,17,17, 24,24,24,24,24,24,24,24,24,24,24,24, 34,34,34,34,34,34, & ! 1. pyr 20,20,20, 19,19,19, 18,18,18, 17,17,17,17,16,17, 24,24,24,24,24,24,24,24,24,24,24,24, 34,34,34,34,34,34, & 20,20,20, 19,19,19, 18,18,18, 17,17,17,17,17,16, 24,24,24,24,24,24,24,24,24,24,24,24, 34,34,34,34,34,34, & @@ -778,7 +779,7 @@ function lattice_interaction_SlipBySlip(Nslip,interactionValues,lattice) result( 30,30,30, 29,29,29, 28,28,28, 27,27,27,27,27,27, 26,26,26,25,26,26,26,26,26,26,26,26, 35,35,35,35,35,35, & 30,30,30, 29,29,29, 28,28,28, 27,27,27,27,27,27, 26,26,26,26,25,26,26,26,26,26,26,26, 35,35,35,35,35,35, & 30,30,30, 29,29,29, 28,28,28, 27,27,27,27,27,27, 26,26,26,26,26,25,26,26,26,26,26,26, 35,35,35,35,35,35, & - 30,30,30, 29,29,29, 28,28,28, 27,27,27,27,27,27, 26,26,26,26,26,26,25,26,26,26,26,26, 35,35,35,35,35,35, & + 30,30,30, 29,29,29, 28,28,28, 27,27,27,27,27,27, 26,26,26,26,26,26,25,26,26,26,26,26, 35,35,35,35,35,35, & ! 1. pyr 30,30,30, 29,29,29, 28,28,28, 27,27,27,27,27,27, 26,26,26,26,26,26,26,25,26,26,26,26, 35,35,35,35,35,35, & 30,30,30, 29,29,29, 28,28,28, 27,27,27,27,27,27, 26,26,26,26,26,26,26,26,25,26,26,26, 35,35,35,35,35,35, & 30,30,30, 29,29,29, 28,28,28, 27,27,27,27,27,27, 26,26,26,26,26,26,26,26,26,25,26,26, 35,35,35,35,35,35, & @@ -788,7 +789,7 @@ function lattice_interaction_SlipBySlip(Nslip,interactionValues,lattice) result( 42,42,42, 41,41,41, 40,40,40, 39,39,39,39,39,39, 38,38,38,38,38,38,38,38,38,38,38,38, 36,37,37,37,37,37, & 42,42,42, 41,41,41, 40,40,40, 39,39,39,39,39,39, 38,38,38,38,38,38,38,38,38,38,38,38, 37,36,37,37,37,37, & 42,42,42, 41,41,41, 40,40,40, 39,39,39,39,39,39, 38,38,38,38,38,38,38,38,38,38,38,38, 37,37,36,37,37,37, & - 42,42,42, 41,41,41, 40,40,40, 39,39,39,39,39,39, 38,38,38,38,38,38,38,38,38,38,38,38, 37,37,37,36,37,37, & + 42,42,42, 41,41,41, 40,40,40, 39,39,39,39,39,39, 38,38,38,38,38,38,38,38,38,38,38,38, 37,37,37,36,37,37, & ! 2. pyr 42,42,42, 41,41,41, 40,40,40, 39,39,39,39,39,39, 38,38,38,38,38,38,38,38,38,38,38,38, 37,37,37,37,36,37, & 42,42,42, 41,41,41, 40,40,40, 39,39,39,39,39,39, 38,38,38,38,38,38,38,38,38,38,38,38, 37,37,37,37,37,36 & ],shape(HEX_INTERACTIONSLIPSLIP)) !< Slip-slip interaction types for hex (onion peel naming scheme) @@ -934,31 +935,32 @@ function lattice_interaction_TwinByTwin(Ntwin,interactionValues,lattice) result( !< 3: other interaction integer, dimension(HEX_NTWIN,HEX_NTWIN), parameter :: & HEX_INTERACTIONTWINTWIN = reshape( [& + ! <-10.1>{10.2} <11.6>{-1-1.1} <10.-2>{10.1} <11.-3>{11.2} 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 7, 7, 7, 7, 7, 7, 13,13,13,13,13,13, & ! -----> acting 2, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 7, 7, 7, 7, 7, 7, 13,13,13,13,13,13, & ! | 2, 2, 1, 2, 2, 2, 3, 3, 3, 3, 3, 3, 7, 7, 7, 7, 7, 7, 13,13,13,13,13,13, & ! | - 2, 2, 2, 1, 2, 2, 3, 3, 3, 3, 3, 3, 7, 7, 7, 7, 7, 7, 13,13,13,13,13,13, & ! v + 2, 2, 2, 1, 2, 2, 3, 3, 3, 3, 3, 3, 7, 7, 7, 7, 7, 7, 13,13,13,13,13,13, & ! v <-10.1>{10.2} 2, 2, 2, 2, 1, 2, 3, 3, 3, 3, 3, 3, 7, 7, 7, 7, 7, 7, 13,13,13,13,13,13, & ! reacting 2, 2, 2, 2, 2, 1, 3, 3, 3, 3, 3, 3, 7, 7, 7, 7, 7, 7, 13,13,13,13,13,13, & 6, 6, 6, 6, 6, 6, 4, 5, 5, 5, 5, 5, 8, 8, 8, 8, 8, 8, 14,14,14,14,14,14, & 6, 6, 6, 6, 6, 6, 5, 4, 5, 5, 5, 5, 8, 8, 8, 8, 8, 8, 14,14,14,14,14,14, & 6, 6, 6, 6, 6, 6, 5, 5, 4, 5, 5, 5, 8, 8, 8, 8, 8, 8, 14,14,14,14,14,14, & - 6, 6, 6, 6, 6, 6, 5, 5, 5, 4, 5, 5, 8, 8, 8, 8, 8, 8, 14,14,14,14,14,14, & + 6, 6, 6, 6, 6, 6, 5, 5, 5, 4, 5, 5, 8, 8, 8, 8, 8, 8, 14,14,14,14,14,14, & ! <11.6>{-1-1.1} 6, 6, 6, 6, 6, 6, 5, 5, 5, 5, 4, 5, 8, 8, 8, 8, 8, 8, 14,14,14,14,14,14, & 6, 6, 6, 6, 6, 6, 5, 5, 5, 5, 5, 4, 8, 8, 8, 8, 8, 8, 14,14,14,14,14,14, & 12,12,12,12,12,12, 11,11,11,11,11,11, 9,10,10,10,10,10, 15,15,15,15,15,15, & 12,12,12,12,12,12, 11,11,11,11,11,11, 10, 9,10,10,10,10, 15,15,15,15,15,15, & 12,12,12,12,12,12, 11,11,11,11,11,11, 10,10, 9,10,10,10, 15,15,15,15,15,15, & - 12,12,12,12,12,12, 11,11,11,11,11,11, 10,10,10, 9,10,10, 15,15,15,15,15,15, & + 12,12,12,12,12,12, 11,11,11,11,11,11, 10,10,10, 9,10,10, 15,15,15,15,15,15, & ! <10.-2>{10.1} 12,12,12,12,12,12, 11,11,11,11,11,11, 10,10,10,10, 9,10, 15,15,15,15,15,15, & 12,12,12,12,12,12, 11,11,11,11,11,11, 10,10,10,10,10, 9, 15,15,15,15,15,15, & 20,20,20,20,20,20, 19,19,19,19,19,19, 18,18,18,18,18,18, 16,17,17,17,17,17, & 20,20,20,20,20,20, 19,19,19,19,19,19, 18,18,18,18,18,18, 17,16,17,17,17,17, & 20,20,20,20,20,20, 19,19,19,19,19,19, 18,18,18,18,18,18, 17,17,16,17,17,17, & - 20,20,20,20,20,20, 19,19,19,19,19,19, 18,18,18,18,18,18, 17,17,17,16,17,17, & + 20,20,20,20,20,20, 19,19,19,19,19,19, 18,18,18,18,18,18, 17,17,17,16,17,17, & ! <11.-3>{11.2} 20,20,20,20,20,20, 19,19,19,19,19,19, 18,18,18,18,18,18, 17,17,17,17,16,17, & 20,20,20,20,20,20, 19,19,19,19,19,19, 18,18,18,18,18,18, 17,17,17,17,17,16 & ],shape(HEX_INTERACTIONTWINTWIN)) !< Twin-twin interaction types for hex @@ -1012,7 +1014,7 @@ function lattice_interaction_TransByTrans(Ntrans,interactionValues,lattice) resu 2,2,2,2,2,2,2,2,2,1,1,1 & ],shape(FCC_INTERACTIONTRANSTRANS)) !< Trans-trans interaction types for fcc - if(lattice == 'cF') then + if (lattice == 'cF') then interactionTypes = FCC_INTERACTIONTRANSTRANS NtransMax = FCC_NTRANSSYSTEM else @@ -1125,22 +1127,23 @@ function lattice_interaction_SlipByTwin(Nslip,Ntwin,interactionValues,lattice) r integer, dimension(HEX_NTWIN,HEX_NSLIP), parameter :: & HEX_INTERACTIONSLIPTWIN = reshape( [& + ! <-10.1>{10.2} <11.6>{-1-1.1} <10.-2>{10.1} <11.-3>{11.2} 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, & ! ----> twin (acting) - 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, & ! | + 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, & ! | basal 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, & ! | ! v 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, & ! slip (reacting) - 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, & + 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, & ! prism 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, & 9, 9, 9, 9, 9, 9, 10,10,10,10,10,10, 11,11,11,11,11,11, 12,12,12,12,12,12, & - 9, 9, 9, 9, 9, 9, 10,10,10,10,10,10, 11,11,11,11,11,11, 12,12,12,12,12,12, & + 9, 9, 9, 9, 9, 9, 10,10,10,10,10,10, 11,11,11,11,11,11, 12,12,12,12,12,12, & ! 2.prism 9, 9, 9, 9, 9, 9, 10,10,10,10,10,10, 11,11,11,11,11,11, 12,12,12,12,12,12, & 13,13,13,13,13,13, 14,14,14,14,14,14, 15,15,15,15,15,15, 16,16,16,16,16,16, & 13,13,13,13,13,13, 14,14,14,14,14,14, 15,15,15,15,15,15, 16,16,16,16,16,16, & 13,13,13,13,13,13, 14,14,14,14,14,14, 15,15,15,15,15,15, 16,16,16,16,16,16, & - 13,13,13,13,13,13, 14,14,14,14,14,14, 15,15,15,15,15,15, 16,16,16,16,16,16, & + 13,13,13,13,13,13, 14,14,14,14,14,14, 15,15,15,15,15,15, 16,16,16,16,16,16, & ! 1. pyr 13,13,13,13,13,13, 14,14,14,14,14,14, 15,15,15,15,15,15, 16,16,16,16,16,16, & 13,13,13,13,13,13, 14,14,14,14,14,14, 15,15,15,15,15,15, 16,16,16,16,16,16, & @@ -1150,7 +1153,7 @@ function lattice_interaction_SlipByTwin(Nslip,Ntwin,interactionValues,lattice) r 17,17,17,17,17,17, 18,18,18,18,18,18, 19,19,19,19,19,19, 20,20,20,20,20,20, & 17,17,17,17,17,17, 18,18,18,18,18,18, 19,19,19,19,19,19, 20,20,20,20,20,20, & 17,17,17,17,17,17, 18,18,18,18,18,18, 19,19,19,19,19,19, 20,20,20,20,20,20, & - 17,17,17,17,17,17, 18,18,18,18,18,18, 19,19,19,19,19,19, 20,20,20,20,20,20, & + 17,17,17,17,17,17, 18,18,18,18,18,18, 19,19,19,19,19,19, 20,20,20,20,20,20, & ! 1. pyr 17,17,17,17,17,17, 18,18,18,18,18,18, 19,19,19,19,19,19, 20,20,20,20,20,20, & 17,17,17,17,17,17, 18,18,18,18,18,18, 19,19,19,19,19,19, 20,20,20,20,20,20, & 17,17,17,17,17,17, 18,18,18,18,18,18, 19,19,19,19,19,19, 20,20,20,20,20,20, & @@ -1160,7 +1163,7 @@ function lattice_interaction_SlipByTwin(Nslip,Ntwin,interactionValues,lattice) r 21,21,21,21,21,21, 22,22,22,22,22,22, 23,23,23,23,23,23, 24,24,24,24,24,24, & 21,21,21,21,21,21, 22,22,22,22,22,22, 23,23,23,23,23,23, 24,24,24,24,24,24, & 21,21,21,21,21,21, 22,22,22,22,22,22, 23,23,23,23,23,23, 24,24,24,24,24,24, & - 21,21,21,21,21,21, 22,22,22,22,22,22, 23,23,23,23,23,23, 24,24,24,24,24,24, & + 21,21,21,21,21,21, 22,22,22,22,22,22, 23,23,23,23,23,23, 24,24,24,24,24,24, & ! 2. pyr 21,21,21,21,21,21, 22,22,22,22,22,22, 23,23,23,23,23,23, 24,24,24,24,24,24, & 21,21,21,21,21,21, 22,22,22,22,22,22, 23,23,23,23,23,23, 24,24,24,24,24,24 & ],shape(HEX_INTERACTIONSLIPTWIN)) !< Slip-twin interaction types for hex @@ -1264,31 +1267,32 @@ function lattice_interaction_TwinBySlip(Ntwin,Nslip,interactionValues,lattice) r integer, dimension(HEX_NSLIP,HEX_NTWIN), parameter :: & HEX_INTERACTIONTWINSLIP = reshape( [& + ! basal prism 2. prism 1. pyr 1. pyr 2. pyr 1, 1, 1, 5, 5, 5, 9, 9, 9, 13,13,13,13,13,13, 17,17,17,17,17,17,17,17,17,17,17,17, 21,21,21,21,21,21, & ! ----> slip (acting) 1, 1, 1, 5, 5, 5, 9, 9, 9, 13,13,13,13,13,13, 17,17,17,17,17,17,17,17,17,17,17,17, 21,21,21,21,21,21, & ! | 1, 1, 1, 5, 5, 5, 9, 9, 9, 13,13,13,13,13,13, 17,17,17,17,17,17,17,17,17,17,17,17, 21,21,21,21,21,21, & ! | - 1, 1, 1, 5, 5, 5, 9, 9, 9, 13,13,13,13,13,13, 17,17,17,17,17,17,17,17,17,17,17,17, 21,21,21,21,21,21, & ! v + 1, 1, 1, 5, 5, 5, 9, 9, 9, 13,13,13,13,13,13, 17,17,17,17,17,17,17,17,17,17,17,17, 21,21,21,21,21,21, & ! v <-10.1>{10.2} 1, 1, 1, 5, 5, 5, 9, 9, 9, 13,13,13,13,13,13, 17,17,17,17,17,17,17,17,17,17,17,17, 21,21,21,21,21,21, & ! twin (reacting) 1, 1, 1, 5, 5, 5, 9, 9, 9, 13,13,13,13,13,13, 17,17,17,17,17,17,17,17,17,17,17,17, 21,21,21,21,21,21, & 2, 2, 2, 6, 6, 6, 10,10,10, 14,14,14,14,14,14, 18,18,18,18,18,18,18,18,18,18,18,18, 22,22,22,22,22,22, & 2, 2, 2, 6, 6, 6, 10,10,10, 14,14,14,14,14,14, 18,18,18,18,18,18,18,18,18,18,18,18, 22,22,22,22,22,22, & 2, 2, 2, 6, 6, 6, 10,10,10, 14,14,14,14,14,14, 18,18,18,18,18,18,18,18,18,18,18,18, 22,22,22,22,22,22, & - 2, 2, 2, 6, 6, 6, 10,10,10, 14,14,14,14,14,14, 18,18,18,18,18,18,18,18,18,18,18,18, 22,22,22,22,22,22, & + 2, 2, 2, 6, 6, 6, 10,10,10, 14,14,14,14,14,14, 18,18,18,18,18,18,18,18,18,18,18,18, 22,22,22,22,22,22, & ! <11.6>{-1-1.1} 2, 2, 2, 6, 6, 6, 10,10,10, 14,14,14,14,14,14, 18,18,18,18,18,18,18,18,18,18,18,18, 22,22,22,22,22,22, & 2, 2, 2, 6, 6, 6, 10,10,10, 14,14,14,14,14,14, 18,18,18,18,18,18,18,18,18,18,18,18, 22,22,22,22,22,22, & 3, 3, 3, 7, 7, 7, 11,11,11, 15,15,15,15,15,15, 19,19,19,19,19,19,19,19,19,19,19,19, 23,23,23,23,23,23, & 3, 3, 3, 7, 7, 7, 11,11,11, 15,15,15,15,15,15, 19,19,19,19,19,19,19,19,19,19,19,19, 23,23,23,23,23,23, & 3, 3, 3, 7, 7, 7, 11,11,11, 15,15,15,15,15,15, 19,19,19,19,19,19,19,19,19,19,19,19, 23,23,23,23,23,23, & - 3, 3, 3, 7, 7, 7, 11,11,11, 15,15,15,15,15,15, 19,19,19,19,19,19,19,19,19,19,19,19, 23,23,23,23,23,23, & + 3, 3, 3, 7, 7, 7, 11,11,11, 15,15,15,15,15,15, 19,19,19,19,19,19,19,19,19,19,19,19, 23,23,23,23,23,23, & ! <10.-2>{10.1} 3, 3, 3, 7, 7, 7, 11,11,11, 15,15,15,15,15,15, 19,19,19,19,19,19,19,19,19,19,19,19, 23,23,23,23,23,23, & 3, 3, 3, 7, 7, 7, 11,11,11, 15,15,15,15,15,15, 19,19,19,19,19,19,19,19,19,19,19,19, 23,23,23,23,23,23, & 4, 4, 4, 8, 8, 8, 12,12,12, 16,16,16,16,16,16, 20,20,20,20,20,20,20,20,20,20,20,20, 24,24,24,24,24,24, & 4, 4, 4, 8, 8, 8, 12,12,12, 16,16,16,16,16,16, 20,20,20,20,20,20,20,20,20,20,20,20, 24,24,24,24,24,24, & 4, 4, 4, 8, 8, 8, 12,12,12, 16,16,16,16,16,16, 20,20,20,20,20,20,20,20,20,20,20,20, 24,24,24,24,24,24, & - 4, 4, 4, 8, 8, 8, 12,12,12, 16,16,16,16,16,16, 20,20,20,20,20,20,20,20,20,20,20,20, 24,24,24,24,24,24, & + 4, 4, 4, 8, 8, 8, 12,12,12, 16,16,16,16,16,16, 20,20,20,20,20,20,20,20,20,20,20,20, 24,24,24,24,24,24, & ! <11.-3>{11.2} 4, 4, 4, 8, 8, 8, 12,12,12, 16,16,16,16,16,16, 20,20,20,20,20,20,20,20,20,20,20,20, 24,24,24,24,24,24, & 4, 4, 4, 8, 8, 8, 12,12,12, 16,16,16,16,16,16, 20,20,20,20,20,20,20,20,20,20,20,20, 24,24,24,24,24,24 & ],shape(HEX_INTERACTIONTWINSLIP)) !< Twin-slip interaction types for hex @@ -1413,7 +1417,7 @@ end function lattice_SchmidMatrix_twin !-------------------------------------------------------------------------------------------------- -!> @brief Schmid matrix for twinning +!> @brief Schmid matrix for transformation !> details only active twin systems are considered !-------------------------------------------------------------------------------------------------- function lattice_SchmidMatrix_trans(Ntrans,lattice_target,cOverA,a_bcc,a_fcc) result(SchmidMatrix) @@ -1539,7 +1543,7 @@ end function lattice_slip_transverse !-------------------------------------------------------------------------------------------------- -!> @brief Labels for slip systems +!> @brief Labels of slip systems !> details only active slip systems are considered !-------------------------------------------------------------------------------------------------- function lattice_labels_slip(Nslip,lattice) result(labels) @@ -2089,12 +2093,12 @@ function lattice_equivalent_nu(C,assumption) result(nu) real(pReal), dimension(6,6) :: S - if (IO_lc(assumption) == 'voigt') then + if (IO_lc(assumption) == 'voigt') then K = (C(1,1)+C(2,2)+C(3,3) +2.0_pReal*(C(1,2)+C(2,3)+C(1,3))) & / 9.0_pReal - elseif(IO_lc(assumption) == 'reuss') then - call math_invert(S,error,C) ! ToDo: correct for Voigt? - if(error) error stop 'matrix inversion failed' + elseif (IO_lc(assumption) == 'reuss') then + call math_invert(S,error,C) + if (error) error stop 'matrix inversion failed' K = 1.0_pReal & / (S(1,1)+S(2,2)+S(3,3) +2.0_pReal*(S(1,2)+S(2,3)+S(1,3))) else @@ -2102,7 +2106,7 @@ function lattice_equivalent_nu(C,assumption) result(nu) endif mu = lattice_equivalent_mu(C,assumption) - nu = (1.5_pReal*K -mu)/(3.0_pReal*K+mu) + nu = (1.5_pReal*K-mu)/(3.0_pReal*K+mu) end function lattice_equivalent_nu @@ -2121,12 +2125,12 @@ function lattice_equivalent_mu(C,assumption) result(mu) real(pReal), dimension(6,6) :: S - if (IO_lc(assumption) == 'voigt') then + if (IO_lc(assumption) == 'voigt') then mu = (1.0_pReal*(C(1,1)+C(2,2)+C(3,3)) -1.0_pReal*(C(1,2)+C(2,3)+C(1,3)) +3.0_pReal*(C(4,4)+C(5,5)+C(6,6))) & / 15.0_pReal - elseif(IO_lc(assumption) == 'reuss') then - call math_invert(S,error,C) ! ToDo: correct for Voigt? - if(error) error stop 'matrix inversion failed' + elseif (IO_lc(assumption) == 'reuss') then + call math_invert(S,error,C) + if (error) error stop 'matrix inversion failed' mu = 15.0_pReal & / (4.0_pReal*(S(1,1)+S(2,2)+S(3,3)) -4.0_pReal*(S(1,2)+S(2,3)+S(1,3)) +3.0_pReal*(S(4,4)+S(5,5)+S(6,6))) else @@ -2155,7 +2159,7 @@ subroutine selfTest system = reshape([1.0_pReal+r(1),0.0_pReal,0.0_pReal, 0.0_pReal,1.0_pReal+r(2),0.0_pReal],[6,1]) CoSy = buildCoordinateSystem([1],[1],system,'cF',0.0_pReal) - if(any(dNeq(CoSy(1:3,1:3,1),math_I3))) error stop 'buildCoordinateSystem' + if (any(dNeq(CoSy(1:3,1:3,1),math_I3))) error stop 'buildCoordinateSystem' do i = 1, 10 call random_number(C) @@ -2201,13 +2205,13 @@ subroutine selfTest C(1,1) = C(1,1) + C(1,2) + 0.1_pReal C(4,4) = 0.5_pReal * (C(1,1) - C(1,2)) C = lattice_symmetrize_C66(C,'cI') - if(dNeq(C(4,4),lattice_equivalent_mu(C,'voigt'),1.0e-12_pReal)) error stop 'equivalent_mu/voigt' - if(dNeq(C(4,4),lattice_equivalent_mu(C,'reuss'),1.0e-12_pReal)) error stop 'equivalent_mu/reuss' + if (dNeq(C(4,4),lattice_equivalent_mu(C,'voigt'),1.0e-12_pReal)) error stop 'equivalent_mu/voigt' + if (dNeq(C(4,4),lattice_equivalent_mu(C,'reuss'),1.0e-12_pReal)) error stop 'equivalent_mu/reuss' lambda = C(1,2) - if(dNeq(lambda*0.5_pReal/(lambda+lattice_equivalent_mu(C,'voigt')), & + if (dNeq(lambda*0.5_pReal/(lambda+lattice_equivalent_mu(C,'voigt')), & lattice_equivalent_nu(C,'voigt'),1.0e-12_pReal)) error stop 'equivalent_nu/voigt' - if(dNeq(lambda*0.5_pReal/(lambda+lattice_equivalent_mu(C,'reuss')), & + if (dNeq(lambda*0.5_pReal/(lambda+lattice_equivalent_mu(C,'reuss')), & lattice_equivalent_nu(C,'reuss'),1.0e-12_pReal)) error stop 'equivalent_nu/reuss' end subroutine selfTest diff --git a/src/math.f90 b/src/math.f90 index aa81a0ed3..85192d077 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -15,15 +15,15 @@ module math implicit none public #if __INTEL_COMPILER >= 1900 - ! do not make use associated entities available to other modules + ! do not make use of associated entities available to other modules private :: & IO, & config #endif real(pReal), parameter :: PI = acos(-1.0_pReal) !< ratio of a circle's circumference to its diameter - real(pReal), parameter :: INDEG = 180.0_pReal/PI !< conversion from radian into degree - real(pReal), parameter :: INRAD = PI/180.0_pReal !< conversion from degree into radian + real(pReal), parameter :: INDEG = 180.0_pReal/PI !< conversion from radian to degree + real(pReal), parameter :: INRAD = PI/180.0_pReal !< conversion from degree to radian complex(pReal), parameter :: TWOPIIMG = cmplx(0.0_pReal,2.0_pReal*PI) !< Re(0.0), Im(2xPi) real(pReal), dimension(3,3), parameter :: & @@ -132,19 +132,19 @@ pure recursive subroutine math_sort(a, istart, iend, sortDim) integer, intent(in),optional :: istart,iend, sortDim integer :: ipivot,s,e,d - if(present(istart)) then + if (present(istart)) then s = istart else s = lbound(a,2) endif - if(present(iend)) then + if (present(iend)) then e = iend else e = ubound(a,2) endif - if(present(sortDim)) then + if (present(sortDim)) then d = sortDim else d = 1 @@ -467,7 +467,7 @@ pure function math_inv33(A) logical :: error call math_invert33(math_inv33,DetA,error,A) - if(error) math_inv33 = 0.0_pReal + if (error) math_inv33 = 0.0_pReal end function math_inv33 @@ -698,7 +698,7 @@ pure function math_sym33to6(m33,weighted) integer :: i - if(present(weighted)) then + if (present(weighted)) then w = merge(NRMMANDEL,1.0_pReal,weighted) else w = NRMMANDEL @@ -725,7 +725,7 @@ pure function math_6toSym33(v6,weighted) integer :: i - if(present(weighted)) then + if (present(weighted)) then w = merge(INVNRMMANDEL,1.0_pReal,weighted) else w = INVNRMMANDEL @@ -802,11 +802,11 @@ pure function math_sym3333to66(m3333,weighted) integer :: i,j - if(present(weighted)) then + if (present(weighted)) then w = merge(NRMMANDEL,1.0_pReal,weighted) else w = NRMMANDEL - endif + end if #ifndef __INTEL_COMPILER do concurrent(i=1:6, j=1:6) @@ -822,7 +822,7 @@ end function math_sym3333to66 !-------------------------------------------------------------------------------------------------- -!> @brief convert 66 matrix into symmetric 3x3x3x3 matrix +!> @brief convert 6x6 matrix into symmetric 3x3x3x3 matrix !> @details Weighted conversion (default) rearranges according to Nye and weights shear ! components according to Mandel. Advisable for matrix operations. ! Unweighted conversion only rearranges order according to Nye @@ -837,7 +837,7 @@ pure function math_66toSym3333(m66,weighted) integer :: i,j - if(present(weighted)) then + if (present(weighted)) then w = merge(INVNRMMANDEL,1.0_pReal,weighted) else w = INVNRMMANDEL @@ -1047,7 +1047,7 @@ subroutine math_eigh33(w,v,m) v(2,2) + m(2, 3) * w(1), & (m(1,1) - w(1)) * (m(2,2) - w(1)) - v(3,2)] norm = norm2(v(1:3, 1)) - fallback1: if(norm < threshold) then + fallback1: if (norm < threshold) then call math_eigh(w,v,error,m) else fallback1 v(1:3,1) = v(1:3, 1) / norm @@ -1055,7 +1055,7 @@ subroutine math_eigh33(w,v,m) v(2,2) + m(2, 3) * w(2), & (m(1,1) - w(2)) * (m(2,2) - w(2)) - v(3,2)] norm = norm2(v(1:3, 2)) - fallback2: if(norm < threshold) then + fallback2: if (norm < threshold) then call math_eigh(w,v,error,m) else fallback2 v(1:3,2) = v(1:3, 2) / norm @@ -1092,7 +1092,7 @@ pure function math_rotationalPart(F) result(R) I_F = [math_trace33(F), 0.5*(math_trace33(F)**2 - math_trace33(matmul(F,F)))] x = math_clip(I_C(1)**2 -3.0_pReal*I_C(2),0.0_pReal)**(3.0_pReal/2.0_pReal) - if(dNeq0(x)) then + if (dNeq0(x)) then Phi = acos(math_clip((I_C(1)**3 -4.5_pReal*I_C(1)*I_C(2) +13.5_pReal*I_C(3))/x,-1.0_pReal,1.0_pReal)) lambda = I_C(1) +(2.0_pReal * sqrt(math_clip(I_C(1)**2-3.0_pReal*I_C(2),0.0_pReal))) & *cos((Phi-2.0_pReal * PI*[1.0_pReal,2.0_pReal,3.0_pReal])/3.0_pReal) @@ -1154,7 +1154,7 @@ function math_eigvalsh33(m) - 2.0_pReal/27.0_pReal*I(1)**3.0_pReal & - I(3) ! different from http://arxiv.org/abs/physics/0610206 (this formulation was in DAMASK) - if(all(abs([P,Q]) < TOL)) then + if (all(abs([P,Q]) < TOL)) then math_eigvalsh33 = math_eigvalsh(m) else rho=sqrt(-3.0_pReal*P**3.0_pReal)/9.0_pReal @@ -1277,7 +1277,7 @@ real(pReal) pure elemental function math_clip(a, left, right) if (present(left)) math_clip = max(left,math_clip) if (present(right)) math_clip = min(right,math_clip) if (present(left) .and. present(right)) then - if(left>right) error stop 'left > right' + if (left>right) error stop 'left > right' endif end function math_clip @@ -1323,38 +1323,38 @@ subroutine selfTest error stop 'math_expand [1,2] by [1,2,3] => [1,2,2,1,1,1]' call math_sort(sort_in_,1,3,2) - if(any(sort_in_ /= sort_out_)) & + if (any(sort_in_ /= sort_out_)) & error stop 'math_sort' - if(any(math_range(5) /= range_out_)) & + if (any(math_range(5) /= range_out_)) & error stop 'math_range' - if(any(dNeq(math_exp33(math_I3,0),math_I3))) & + if (any(dNeq(math_exp33(math_I3,0),math_I3))) & error stop 'math_exp33(math_I3,1)' - if(any(dNeq(math_exp33(math_I3,128),exp(1.0_pReal)*math_I3))) & + if (any(dNeq(math_exp33(math_I3,128),exp(1.0_pReal)*math_I3))) & error stop 'math_exp33(math_I3,128)' call random_number(v9) - if(any(dNeq(math_33to9(math_9to33(v9)),v9))) & + if (any(dNeq(math_33to9(math_9to33(v9)),v9))) & error stop 'math_33to9/math_9to33' call random_number(t99) - if(any(dNeq(math_3333to99(math_99to3333(t99)),t99))) & + if (any(dNeq(math_3333to99(math_99to3333(t99)),t99))) & error stop 'math_3333to99/math_99to3333' call random_number(v6) - if(any(dNeq(math_sym33to6(math_6toSym33(v6)),v6))) & + if (any(dNeq(math_sym33to6(math_6toSym33(v6)),v6))) & error stop 'math_sym33to6/math_6toSym33' call random_number(t66) - if(any(dNeq(math_sym3333to66(math_66toSym3333(t66)),t66,1.0e-15_pReal))) & + if (any(dNeq(math_sym3333to66(math_66toSym3333(t66)),t66,1.0e-15_pReal))) & error stop 'math_sym3333to66/math_66toSym3333' - if(any(dNeq(math_3333toVoigt66(math_Voigt66to3333(t66)),t66,1.0e-15_pReal))) & + if (any(dNeq(math_3333toVoigt66(math_Voigt66to3333(t66)),t66,1.0e-15_pReal))) & error stop 'math_3333toVoigt66/math_Voigt66to3333' call random_number(v6) - if(any(dNeq0(math_6toSym33(v6) - math_symmetric33(math_6toSym33(v6))))) & + if (any(dNeq0(math_6toSym33(v6) - math_symmetric33(math_6toSym33(v6))))) & error stop 'math_symmetric33' call random_number(v3_1) @@ -1362,34 +1362,34 @@ subroutine selfTest call random_number(v3_3) call random_number(v3_4) - if(dNeq(abs(dot_product(math_cross(v3_1-v3_4,v3_2-v3_4),v3_3-v3_4))/6.0, & + if (dNeq(abs(dot_product(math_cross(v3_1-v3_4,v3_2-v3_4),v3_3-v3_4))/6.0, & math_volTetrahedron(v3_1,v3_2,v3_3,v3_4),tol=1.0e-12_pReal)) & error stop 'math_volTetrahedron' call random_number(t33) - if(dNeq(math_det33(math_symmetric33(t33)),math_detSym33(math_symmetric33(t33)),tol=1.0e-12_pReal)) & + if (dNeq(math_det33(math_symmetric33(t33)),math_detSym33(math_symmetric33(t33)),tol=1.0e-12_pReal)) & error stop 'math_det33/math_detSym33' - if(any(dNeq(t33+transpose(t33),math_mul3333xx33(math_identity4th(),t33+transpose(t33))))) & + if (any(dNeq(t33+transpose(t33),math_mul3333xx33(math_identity4th(),t33+transpose(t33))))) & error stop 'math_mul3333xx33/math_identity4th' - if(any(dNeq0(math_eye(3),math_inv33(math_I3)))) & + if (any(dNeq0(math_eye(3),math_inv33(math_I3)))) & error stop 'math_inv33(math_I3)' do while(abs(math_det33(t33))<1.0e-9_pReal) call random_number(t33) enddo - if(any(dNeq0(matmul(t33,math_inv33(t33)) - math_eye(3),tol=1.0e-9_pReal))) & + if (any(dNeq0(matmul(t33,math_inv33(t33)) - math_eye(3),tol=1.0e-9_pReal))) & error stop 'math_inv33' call math_invert33(t33_2,det,e,t33) - if(any(dNeq0(matmul(t33,t33_2) - math_eye(3),tol=1.0e-9_pReal)) .or. e) & + if (any(dNeq0(matmul(t33,t33_2) - math_eye(3),tol=1.0e-9_pReal)) .or. e) & error stop 'math_invert33: T:T^-1 != I' - if(dNeq(det,math_det33(t33),tol=1.0e-12_pReal)) & + if (dNeq(det,math_det33(t33),tol=1.0e-12_pReal)) & error stop 'math_invert33 (determinant)' call math_invert(t33_2,e,t33) - if(any(dNeq0(matmul(t33,t33_2) - math_eye(3),tol=1.0e-9_pReal)) .or. e) & + if (any(dNeq0(matmul(t33,t33_2) - math_eye(3),tol=1.0e-9_pReal)) .or. e) & error stop 'math_invert t33' do while(math_det33(t33)<1.0e-2_pReal) ! O(det(F)) = 1 @@ -1397,7 +1397,7 @@ subroutine selfTest enddo t33_2 = math_rotationalPart(transpose(t33)) t33 = math_rotationalPart(t33) - if(any(dNeq0(matmul(t33_2,t33) - math_I3,tol=1.0e-10_pReal))) & + if (any(dNeq0(matmul(t33_2,t33) - math_I3,tol=1.0e-10_pReal))) & error stop 'math_rotationalPart' call random_number(r) @@ -1405,33 +1405,33 @@ subroutine selfTest txx = math_eye(d) allocate(txx_2(d,d)) call math_invert(txx_2,e,txx) - if(any(dNeq0(txx_2,txx) .or. e)) & + if (any(dNeq0(txx_2,txx) .or. e)) & error stop 'math_invert(txx)/math_eye' call math_invert(t99_2,e,t99) ! not sure how likely it is that we get a singular matrix - if(any(dNeq0(matmul(t99_2,t99)-math_eye(9),tol=1.0e-9_pReal)) .or. e) & + if (any(dNeq0(matmul(t99_2,t99)-math_eye(9),tol=1.0e-9_pReal)) .or. e) & error stop 'math_invert(t99)' - if(any(dNeq(math_clip([4.0_pReal,9.0_pReal],5.0_pReal,6.5_pReal),[5.0_pReal,6.5_pReal]))) & + if (any(dNeq(math_clip([4.0_pReal,9.0_pReal],5.0_pReal,6.5_pReal),[5.0_pReal,6.5_pReal]))) & error stop 'math_clip' - if(math_factorial(10) /= 3628800) & + if (math_factorial(10) /= 3628800) & error stop 'math_factorial' - if(math_binomial(49,6) /= 13983816) & + if (math_binomial(49,6) /= 13983816) & error stop 'math_binomial' - if(math_multinomial([1,2,3,4]) /= 12600) & + if (math_multinomial([1,2,3,4]) /= 12600) & error stop 'math_multinomial' ijk = cshift([1,2,3],int(r*1.0e2_pReal)) - if(dNeq(math_LeviCivita(ijk(1),ijk(2),ijk(3)),+1.0_pReal)) & + if (dNeq(math_LeviCivita(ijk(1),ijk(2),ijk(3)),+1.0_pReal)) & error stop 'math_LeviCivita(even)' ijk = cshift([3,2,1],int(r*2.0e2_pReal)) - if(dNeq(math_LeviCivita(ijk(1),ijk(2),ijk(3)),-1.0_pReal)) & + if (dNeq(math_LeviCivita(ijk(1),ijk(2),ijk(3)),-1.0_pReal)) & error stop 'math_LeviCivita(odd)' ijk = cshift([2,2,1],int(r*2.0e2_pReal)) - if(dNeq0(math_LeviCivita(ijk(1),ijk(2),ijk(3)))) & + if (dNeq0(math_LeviCivita(ijk(1),ijk(2),ijk(3)))) & error stop 'math_LeviCivita' end subroutine selfTest diff --git a/src/phase.f90 b/src/phase.f90 index 22c35416b..214b0f5fa 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -5,6 +5,7 @@ !-------------------------------------------------------------------------------------------------- module phase use prec + use constants use math use rotations use IO diff --git a/src/phase_mechanical_eigen_thermalexpansion.f90 b/src/phase_mechanical_eigen_thermalexpansion.f90 index 1c08140ae..a84dd7797 100644 --- a/src/phase_mechanical_eigen_thermalexpansion.f90 +++ b/src/phase_mechanical_eigen_thermalexpansion.f90 @@ -58,14 +58,14 @@ module function thermalexpansion_init(kinematics_length) result(myKinematics) associate(prm => param(kinematics_thermal_expansion_instance(p))) kinematic_type => kinematics%get(k) - prm%T_ref = kinematic_type%get_asFloat('T_ref', defaultVal=0.0_pReal) + prm%T_ref = kinematic_type%get_asFloat('T_ref', defaultVal=T_ROOM) prm%A(1,1,1) = kinematic_type%get_asFloat('A_11') - prm%A(1,1,2) = kinematic_type%get_asFloat('A_11,T',defaultVal=0.0_pReal) + prm%A(1,1,2) = kinematic_type%get_asFloat('A_11,T', defaultVal=0.0_pReal) prm%A(1,1,3) = kinematic_type%get_asFloat('A_11,T^2',defaultVal=0.0_pReal) if (any(phase_lattice(p) == ['hP','tI'])) then prm%A(3,3,1) = kinematic_type%get_asFloat('A_33') - prm%A(3,3,2) = kinematic_type%get_asFloat('A_33,T',defaultVal=0.0_pReal) + prm%A(3,3,2) = kinematic_type%get_asFloat('A_33,T', defaultVal=0.0_pReal) prm%A(3,3,3) = kinematic_type%get_asFloat('A_33,T^2',defaultVal=0.0_pReal) end if do i=1, size(prm%A,3) @@ -98,14 +98,14 @@ module subroutine thermalexpansion_LiAndItsTangent(Li, dLi_dTstar, ph,me) associate(prm => param(kinematics_thermal_expansion_instance(ph))) Li = dot_T * ( & - prm%A(1:3,1:3,1)*(T - prm%T_ref)**0 & ! constant coefficient - + prm%A(1:3,1:3,2)*(T - prm%T_ref)**1 & ! linear coefficient - + prm%A(1:3,1:3,3)*(T - prm%T_ref)**2 & ! quadratic coefficient + prm%A(1:3,1:3,1) & ! constant coefficient + + prm%A(1:3,1:3,2)*(T - prm%T_ref)**1 & ! linear coefficient + + prm%A(1:3,1:3,3)*(T - prm%T_ref)**2 & ! quadratic coefficient ) / & (1.0_pReal & - + prm%A(1:3,1:3,1)*(T - prm%T_ref)**1 / 1. & - + prm%A(1:3,1:3,2)*(T - prm%T_ref)**2 / 2. & - + prm%A(1:3,1:3,3)*(T - prm%T_ref)**3 / 3. & + + prm%A(1:3,1:3,1)*(T - prm%T_ref)**1 / 1.0_pReal & + + prm%A(1:3,1:3,2)*(T - prm%T_ref)**2 / 2.0_pReal & + + prm%A(1:3,1:3,3)*(T - prm%T_ref)**3 / 3.0_pReal & ) end associate dLi_dTstar = 0.0_pReal diff --git a/src/phase_mechanical_elastic.f90 b/src/phase_mechanical_elastic.f90 index bb0422bd1..ea113ddfb 100644 --- a/src/phase_mechanical_elastic.f90 +++ b/src/phase_mechanical_elastic.f90 @@ -1,13 +1,15 @@ submodule(phase:mechanical) elastic type :: tParameters - real(pReal) :: & + real(pReal),dimension(3) :: & C_11 = 0.0_pReal, & C_12 = 0.0_pReal, & C_13 = 0.0_pReal, & C_33 = 0.0_pReal, & C_44 = 0.0_pReal, & C_66 = 0.0_pReal + real(pReal) :: & + T_ref end type tParameters type(tParameters), allocatable, dimension(:) :: param @@ -15,7 +17,7 @@ submodule(phase:mechanical) elastic contains !-------------------------------------------------------------------------------------------------- -!> @brief Initialize elasticity. +!> @brief initialize elasticity !-------------------------------------------------------------------------------------------------- module subroutine elastic_init(phases) @@ -28,7 +30,7 @@ module subroutine elastic_init(phases) phase, & mech, & elastic - + logical :: thermal_active print'(/,1x,a)', '<<<+- phase:mechanical:elastic init -+>>>' print'(/,1x,a)', '<<<+- phase:mechanical:elastic:Hooke init -+>>>' @@ -45,15 +47,35 @@ module subroutine elastic_init(phases) associate(prm => param(ph)) - prm%C_11 = elastic%get_asFloat('C_11') - prm%C_12 = elastic%get_asFloat('C_12') - prm%C_44 = elastic%get_asFloat('C_44') + prm%T_ref = elastic%get_asFloat('T_ref', defaultVal=T_ROOM) + prm%C_11(1) = elastic%get_asFloat('C_11') + prm%C_11(2) = elastic%get_asFloat('C_11,T', defaultVal=0.0_pReal) + prm%C_11(3) = elastic%get_asFloat('C_11,T^2',defaultVal=0.0_pReal) + + prm%C_12(1) = elastic%get_asFloat('C_12') + prm%C_12(2) = elastic%get_asFloat('C_12,T', defaultVal=0.0_pReal) + prm%C_12(3) = elastic%get_asFloat('C_12,T^2',defaultVal=0.0_pReal) + + prm%C_44(1) = elastic%get_asFloat('C_44') + prm%C_44(2) = elastic%get_asFloat('C_44,T', defaultVal=0.0_pReal) + prm%C_44(3) = elastic%get_asFloat('C_44,T^2',defaultVal=0.0_pReal) + if (any(phase_lattice(ph) == ['hP','tI'])) then - prm%C_13 = elastic%get_asFloat('C_13') - prm%C_33 = elastic%get_asFloat('C_33') + prm%C_13(1) = elastic%get_asFloat('C_13') + prm%C_13(2) = elastic%get_asFloat('C_13,T', defaultVal=0.0_pReal) + prm%C_13(3) = elastic%get_asFloat('C_13,T^2',defaultVal=0.0_pReal) + + prm%C_33(1) = elastic%get_asFloat('C_33') + prm%C_33(2) = elastic%get_asFloat('C_33,T', defaultVal=0.0_pReal) + prm%C_33(3) = elastic%get_asFloat('C_33,T^2',defaultVal=0.0_pReal) + end if + + if (phase_lattice(ph) == 'tI') then + prm%C_66(1) = elastic%get_asFloat('C_66') + prm%C_66(2) = elastic%get_asFloat('C_66,T', defaultVal=0.0_pReal) + prm%C_66(3) = elastic%get_asFloat('C_66,T^2',defaultVal=0.0_pReal) end if - if (phase_lattice(ph) == 'tI') prm%C_66 = elastic%get_asFloat('C_66') end associate end do @@ -62,28 +84,51 @@ end subroutine elastic_init !-------------------------------------------------------------------------------------------------- -!> @brief Return 6x6 elasticity tensor. +!> @brief return 6x6 elasticity tensor !-------------------------------------------------------------------------------------------------- module function elastic_C66(ph,en) result(C66) integer, intent(in) :: & ph, & en + real(pReal), dimension(6,6) :: C66 + real(pReal) :: T associate(prm => param(ph)) C66 = 0.0_pReal - C66(1,1) = prm%C_11 - C66(1,2) = prm%C_12 - C66(4,4) = prm%C_44 + T = thermal_T(ph,en) + + C66(1,1) = prm%C_11(1) & + + prm%C_11(2)*(T - prm%T_ref)**1 & + + prm%C_11(3)*(T - prm%T_ref)**2 + + C66(1,2) = prm%C_12(1) & + + prm%C_12(2)*(T - prm%T_ref)**1 & + + prm%C_12(3)*(T - prm%T_ref)**2 + + C66(4,4) = prm%C_44(1) & + + prm%C_44(2)*(T - prm%T_ref)**1 & + + prm%C_44(3)*(T - prm%T_ref)**2 + if (any(phase_lattice(ph) == ['hP','tI'])) then - C66(1,3) = prm%C_13 - C66(3,3) = prm%C_33 + C66(1,3) = prm%C_13(1) & + + prm%C_13(2)*(T - prm%T_ref)**1 & + + prm%C_13(3)*(T - prm%T_ref)**2 + + C66(3,3) = prm%C_33(1) & + + prm%C_33(2)*(T - prm%T_ref)**1 & + + prm%C_33(3)*(T - prm%T_ref)**2 + end if - if (phase_lattice(ph) == 'tI') C66(6,6) = prm%C_66 + if (phase_lattice(ph) == 'tI') then + C66(6,6) = prm%C_66(1) & + + prm%C_66(2)*(T - prm%T_ref)**1 & + + prm%C_66(3)*(T - prm%T_ref)**2 + end if C66 = lattice_symmetrize_C66(C66,phase_lattice(ph)) @@ -93,7 +138,7 @@ end function elastic_C66 !-------------------------------------------------------------------------------------------------- -!> @brief Return shear modulus. +!> @brief return shear modulus !-------------------------------------------------------------------------------------------------- module function elastic_mu(ph,en) result(mu) @@ -110,7 +155,7 @@ end function elastic_mu !-------------------------------------------------------------------------------------------------- -!> @brief Return Poisson ratio. +!> @brief return Poisson ratio !-------------------------------------------------------------------------------------------------- module function elastic_nu(ph,en) result(nu) @@ -128,7 +173,7 @@ end function elastic_nu !-------------------------------------------------------------------------------------------------- -!> @brief returns the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to +!> @brief return the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to !> the elastic and intermediate deformation gradients using Hooke's law ! ToDo: Use Voigt matrix directly !-------------------------------------------------------------------------------------------------- diff --git a/src/phase_mechanical_plastic_dislotungsten.f90 b/src/phase_mechanical_plastic_dislotungsten.f90 index 102e009fe..1dab5dfd4 100644 --- a/src/phase_mechanical_plastic_dislotungsten.f90 +++ b/src/phase_mechanical_plastic_dislotungsten.f90 @@ -7,9 +7,6 @@ !-------------------------------------------------------------------------------------------------- submodule(phase:plastic) dislotungsten - real(pReal), parameter :: & - kB = 1.38e-23_pReal !< Boltzmann constant in J/Kelvin - type :: tParameters real(pReal) :: & D = 1.0_pReal, & !< grain size @@ -344,7 +341,7 @@ module subroutine dislotungsten_dotState(Mp,T,ph,en) dot_rho_dip_formation = merge(2.0_pReal*(d_hat-prm%d_caron)*stt%rho_mob(:,en)*dot%gamma_sl(:,en)/prm%b_sl, & 0.0_pReal, & prm%dipoleformation) - v_cl = (3.0_pReal*mu*prm%D_0*exp(-prm%Q_cl/(kB*T))*prm%f_at/(2.0_pReal*PI*kB*T)) & + v_cl = (3.0_pReal*mu*prm%D_0*exp(-prm%Q_cl/(K_B*T))*prm%f_at/(2.0_pReal*PI*K_B*T)) & * (1.0_pReal/(d_hat+prm%d_caron)) dot_rho_dip_climb = (4.0_pReal*v_cl*stt%rho_dip(:,en))/(d_hat-prm%d_caron) ! ToDo: Discuss with Franz: Stress dependency? end where @@ -475,7 +472,7 @@ pure subroutine kinetics(Mp,T,ph,en, & if (present(tau_pos_out)) tau_pos_out = tau_pos if (present(tau_neg_out)) tau_neg_out = tau_neg - associate(BoltzmannRatio => prm%Q_s/(kB*T), & + associate(BoltzmannRatio => prm%Q_s/(K_B*T), & b_rho_half => stt%rho_mob(:,en) * prm%b_sl * 0.5_pReal, & effectiveLength => dst%Lambda_sl(:,en) - prm%w) diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90 index de027ae44..034b3f503 100644 --- a/src/phase_mechanical_plastic_dislotwin.f90 +++ b/src/phase_mechanical_plastic_dislotwin.f90 @@ -9,9 +9,6 @@ !-------------------------------------------------------------------------------------------------- submodule(phase:plastic) dislotwin - real(pReal), parameter :: & - kB = 1.38e-23_pReal !< Boltzmann constant in J/Kelvin - type :: tParameters real(pReal) :: & Q_cl = 1.0_pReal, & !< activation energy for dislocation climb @@ -31,7 +28,7 @@ submodule(phase:plastic) dislotwin delta_G = 1.0_pReal, & !< Free energy difference between austensite and martensite i_tr = 1.0_pReal, & !< adjustment parameter to calculate MFP for transformation h = 1.0_pReal, & !< Stack height of hex nucleus - T_ref = 0.0_pReal, & + T_ref = T_ROOM, & a_cI = 1.0_pReal, & a_cF = 1.0_pReal real(pReal), dimension(2) :: & @@ -150,7 +147,7 @@ module function plastic_dislotwin_init() result(myPlasticity) myPlasticity = plastic_active('dislotwin') - if(count(myPlasticity) == 0) return + if (count(myPlasticity) == 0) return print'(/,1x,a)', '<<<+- phase:mechanical:plastic:dislotwin init -+>>>' print'(/,a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT) @@ -173,7 +170,7 @@ module function plastic_dislotwin_init() result(myPlasticity) do ph = 1, phases%length - if(.not. myPlasticity(ph)) cycle + if (.not. myPlasticity(ph)) cycle associate(prm => param(ph), dot => dotState(ph), stt => state(ph), dst => dependentState(ph)) @@ -368,7 +365,7 @@ module function plastic_dislotwin_init() result(myPlasticity) !-------------------------------------------------------------------------------------------------- ! parameters required for several mechanisms and their interactions - if(prm%sum_N_sl + prm%sum_N_tw + prm%sum_N_tw > 0) & + if (prm%sum_N_sl + prm%sum_N_tw + prm%sum_N_tw > 0) & prm%D = pl%get_asFloat('D') if (prm%sum_N_tw + prm%sum_N_tr > 0) & @@ -425,7 +422,7 @@ module function plastic_dislotwin_init() result(myPlasticity) stt%gamma_sl=>plasticState(ph)%state(startIndex:endIndex,:) dot%gamma_sl=>plasticState(ph)%dotState(startIndex:endIndex,:) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal) - if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma' + if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma' startIndex = endIndex + 1 endIndex = endIndex + prm%sum_N_tw @@ -473,9 +470,7 @@ module function plastic_dislotwin_homogenizedC(ph,en) result(homogenizedC) integer, intent(in) :: & ph, en real(pReal), dimension(6,6) :: & - homogenizedC - - real(pReal), dimension(6,6) :: & + homogenizedC, & C real(pReal), dimension(:,:,:), allocatable :: & C66_tw, & @@ -500,8 +495,8 @@ module function plastic_dislotwin_homogenizedC(ph,en) result(homogenizedC) homogenizedC = homogenizedC & + stt%f_tw(i,en)*C66_tw(1:6,1:6,i) end do - end if twinActive - + end if twinActive + transActive: if (prm%sum_N_tr > 0) then C66_tr = lattice_C66_trans(prm%N_tr,C,prm%lattice_tr,0.0_pReal,prm%a_cI,prm%a_cF) do i=1,prm%sum_N_tr @@ -597,9 +592,9 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en) Lp = Lp * f_unrotated dLp_dMp = dLp_dMp * f_unrotated - shearBandingContribution: if(dNeq0(prm%v_sb)) then + shearBandingContribution: if (dNeq0(prm%v_sb)) then - E_kB_T = prm%E_sb/(kB*T) + E_kB_T = prm%E_sb/(K_B*T) call math_eigh33(eigValues,eigVectors,Mp) ! is Mp symmetric by design? do i = 1,6 @@ -696,8 +691,8 @@ module subroutine dislotwin_dotState(Mp,T,ph,en) * (prm%Gamma_sf(1) + prm%Gamma_sf(2) * T) / (mu*prm%b_sl(i)), & 1.0_pReal, & prm%ExtendedDislocations) - v_cl = 2.0_pReal*prm%omega*b_d**2.0_pReal*exp(-prm%Q_cl/(kB*T)) & - * (exp(abs(sigma_cl)*prm%b_sl(i)**3.0_pReal/(kB*T)) - 1.0_pReal) + v_cl = 2.0_pReal*prm%omega*b_d**2.0_pReal*exp(-prm%Q_cl/(K_B*T)) & + * (exp(abs(sigma_cl)*prm%b_sl(i)**3.0_pReal/(K_B*T)) - 1.0_pReal) dot_rho_dip_climb(i) = 4.0_pReal*v_cl*stt%rho_dip(i,en) & / (d_hat-prm%d_caron(i)) @@ -848,7 +843,7 @@ module subroutine plastic_dislotwin_results(ph,group) 'threshold stress for twinning','Pa',prm%systems_tw) case('f_tr') - if(prm%sum_N_tr>0) call results_writeDataset(stt%f_tr,group,trim(prm%output(ou)), & + if (prm%sum_N_tr>0) call results_writeDataset(stt%f_tr,group,trim(prm%output(ou)), & 'martensite volume fraction','m³/m³') end select @@ -909,7 +904,7 @@ pure subroutine kinetics_sl(Mp,T,ph,en, & significantStress: where(tau_eff > tol_math_check) stressRatio = tau_eff/prm%tau_0 StressRatio_p = stressRatio** prm%p - Q_kB_T = prm%Q_sl/(kB*T) + Q_kB_T = prm%Q_sl/(K_B*T) v_wait_inverse = exp(Q_kB_T*(1.0_pReal-StressRatio_p)** prm%q) & / prm%v_0 v_run_inverse = prm%B/(tau_eff*prm%b_sl) @@ -931,8 +926,8 @@ pure subroutine kinetics_sl(Mp,T,ph,en, & end associate - if(present(ddot_gamma_dtau_sl)) ddot_gamma_dtau_sl = ddot_gamma_dtau - if(present(tau_sl)) tau_sl = tau + if (present(ddot_gamma_dtau_sl)) ddot_gamma_dtau_sl = ddot_gamma_dtau + if (present(tau_sl)) tau_sl = tau end subroutine kinetics_sl @@ -982,7 +977,7 @@ pure subroutine kinetics_tw(Mp,T,dot_gamma_sl,ph,en,& Ndot0=(abs(dot_gamma_sl(s1))*(stt%rho_mob(s2,en)+stt%rho_dip(s2,en))+& abs(dot_gamma_sl(s2))*(stt%rho_mob(s1,en)+stt%rho_dip(s1,en)))/& (prm%L_tw*prm%b_sl(i))*& - (1.0_pReal-exp(-prm%V_cs/(kB*T)*(dst%tau_r_tw(i,en)-tau(i)))) + (1.0_pReal-exp(-prm%V_cs/(K_B*T)*(dst%tau_r_tw(i,en)-tau(i)))) else Ndot0=0.0_pReal end if @@ -1002,7 +997,7 @@ pure subroutine kinetics_tw(Mp,T,dot_gamma_sl,ph,en,& end associate - if(present(ddot_gamma_dtau_tw)) ddot_gamma_dtau_tw = ddot_gamma_dtau + if (present(ddot_gamma_dtau_tw)) ddot_gamma_dtau_tw = ddot_gamma_dtau end subroutine kinetics_tw @@ -1051,7 +1046,7 @@ pure subroutine kinetics_tr(Mp,T,dot_gamma_sl,ph,en,& Ndot0=(abs(dot_gamma_sl(s1))*(stt%rho_mob(s2,en)+stt%rho_dip(s2,en))+& abs(dot_gamma_sl(s2))*(stt%rho_mob(s1,en)+stt%rho_dip(s1,en)))/& (prm%L_tr*prm%b_sl(i))*& - (1.0_pReal-exp(-prm%V_cs/(kB*T)*(dst%tau_r_tr(i,en)-tau(i)))) + (1.0_pReal-exp(-prm%V_cs/(K_B*T)*(dst%tau_r_tr(i,en)-tau(i)))) else Ndot0=0.0_pReal end if @@ -1071,7 +1066,7 @@ pure subroutine kinetics_tr(Mp,T,dot_gamma_sl,ph,en,& end associate - if(present(ddot_gamma_dtau_tr)) ddot_gamma_dtau_tr = ddot_gamma_dtau + if (present(ddot_gamma_dtau_tr)) ddot_gamma_dtau_tr = ddot_gamma_dtau end subroutine kinetics_tr diff --git a/src/phase_mechanical_plastic_nonlocal.f90 b/src/phase_mechanical_plastic_nonlocal.f90 index f53a16042..09c5f4d0f 100644 --- a/src/phase_mechanical_plastic_nonlocal.f90 +++ b/src/phase_mechanical_plastic_nonlocal.f90 @@ -19,9 +19,6 @@ submodule(phase:plastic) nonlocal type(tGeometry), dimension(:), allocatable :: geom - real(pReal), parameter :: & - kB = 1.38e-23_pReal !< Boltzmann constant in J/Kelvin - ! storage order of dislocation types integer, dimension(*), parameter :: & sgl = [1,2,3,4,5,6,7,8] !< signed (single) @@ -1094,9 +1091,9 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, & ! thermally activated annihilation of edge dipoles by climb rhoDotThermalAnnihilation = 0.0_pReal - D_SD = prm%D_0 * exp(-prm%Q_cl / (kB * Temperature)) ! eq. 3.53 + D_SD = prm%D_0 * exp(-prm%Q_cl / (K_B * Temperature)) ! eq. 3.53 v_climb = D_SD * mu * prm%V_at & - / (PI * (1.0_pReal-nu) * (dUpper(:,1) + dLower(:,1)) * kB * Temperature) ! eq. 3.54 + / (PI * (1.0_pReal-nu) * (dUpper(:,1) + dLower(:,1)) * K_B * Temperature) ! eq. 3.54 forall (s = 1:prm%sum_N_sl, dUpper(s,1) > dLower(s,1)) & rhoDotThermalAnnihilation(s,9) = max(- 4.0_pReal * rhoDip(s,1) * v_climb(s) / (dUpper(s,1) - dLower(s,1)), & - rhoDip(s,1) / timestep - rhoDotAthermalAnnihilation(s,9) & @@ -1671,9 +1668,9 @@ pure subroutine kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, tauThreshold, c, T, activationEnergy_P = criticalStress_P * activationVolume_P tauRel_P = min(1.0_pReal, tauEff / criticalStress_P) tPeierls = 1.0_pReal / prm%nu_a & - * exp(activationEnergy_P / (kB * T) & + * exp(activationEnergy_P / (K_B * T) & * (1.0_pReal - tauRel_P**prm%p)**prm%q) - dtPeierls_dtau = merge(tPeierls * prm%p * prm%q * activationVolume_P / (kB * T) & + dtPeierls_dtau = merge(tPeierls * prm%p * prm%q * activationVolume_P / (K_B * T) & * (1.0_pReal - tauRel_P**prm%p)**(prm%q-1.0_pReal) * tauRel_P**(prm%p-1.0_pReal), & 0.0_pReal, & tauEff < criticalStress_P) @@ -1685,8 +1682,8 @@ pure subroutine kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, tauThreshold, c, T, criticalStress_S = prm%Q_sol / activationVolume_S tauRel_S = min(1.0_pReal, tauEff / criticalStress_S) tSolidSolution = 1.0_pReal / prm%nu_a & - * exp(prm%Q_sol / (kB * T)* (1.0_pReal - tauRel_S**prm%p)**prm%q) - dtSolidSolution_dtau = merge(tSolidSolution * prm%p * prm%q * activationVolume_S / (kB * T) & + * exp(prm%Q_sol / (K_B * T)* (1.0_pReal - tauRel_S**prm%p)**prm%q) + dtSolidSolution_dtau = merge(tSolidSolution * prm%p * prm%q * activationVolume_S / (K_B * T) & * (1.0_pReal - tauRel_S**prm%p)**(prm%q-1.0_pReal)* tauRel_S**(prm%p-1.0_pReal), & 0.0_pReal, & tauEff < criticalStress_S)