From 29cbf1304b2ed5420fd3aa9f3c1fe1f040191cd1 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 12 Jan 2022 07:48:09 +0100 Subject: [PATCH 1/6] grain growth not maintained and has issues the grain growth model is based on the Voronoi Implicit Interface Method (https://doi.org/10.1016/j.jcp.2012.04.004). The last step in this algorithm is the assignment of the new phase/material ID to the voxels in the 'thick boundary' which is done with distance_transform_edt from ndimage. This problem can have multiple solution and can lead to the translation of grains. In the original publication, the position of the boundary is calculated with subvoxel resolution by solving the eikonal equation. The following python packages might help: https://pypi.org/project/eikonalfm https://pypi.org/project/scikit-fmm https://github.com/malcolmw/pykonal --- processing/legacy/geom_grainGrowth.py | 178 -------------------------- 1 file changed, 178 deletions(-) delete mode 100755 processing/legacy/geom_grainGrowth.py diff --git a/processing/legacy/geom_grainGrowth.py b/processing/legacy/geom_grainGrowth.py deleted file mode 100755 index d969b415c..000000000 --- a/processing/legacy/geom_grainGrowth.py +++ /dev/null @@ -1,178 +0,0 @@ -#!/usr/bin/env python3 - -import os -import sys -from io import StringIO -from optparse import OptionParser - -import numpy as np -from scipy import ndimage - -import damask - - -scriptName = os.path.splitext(os.path.basename(__file__))[0] -scriptID = ' '.join([scriptName,damask.version]) - - -getInterfaceEnergy = lambda A,B: np.float32((A != B)*1.0) # 1.0 if A & B are distinct, 0.0 otherwise -struc = ndimage.generate_binary_structure(3,1) # 3D von Neumann neighborhood - - -#-------------------------------------------------------------------------------------------------- -# MAIN -#-------------------------------------------------------------------------------------------------- - -parser = OptionParser(option_class=damask.extendableOption, usage='%prog option(s) [geomfile(s)]', description = """ -Smoothen interface roughness by simulated curvature flow. -This is achieved by the diffusion of each initially sharply bounded grain volume within the periodic domain -up to a given distance 'd' voxels. -The final geometry is assembled by selecting at each voxel that grain index for which the concentration remains largest. -References 10.1073/pnas.1111557108 (10.1006/jcph.1994.1105) - -""", version = scriptID) - -parser.add_option('-d', '--distance', - dest = 'd', - type = 'float', metavar = 'float', - help = 'diffusion distance in voxels [%default]') -parser.add_option('-N', '--iterations', - dest = 'N', - type = 'int', metavar = 'int', - help = 'curvature flow iterations [%default]') -parser.add_option('-i', '--immutable', - action = 'extend', dest = 'immutable', metavar = '', - help = 'list of immutable material indices') -parser.add_option('--ndimage', - dest = 'ndimage', action='store_true', - help = 'use ndimage.gaussian_filter in lieu of explicit FFT') - -parser.set_defaults(d = 1, - N = 1, - immutable = [], - ndimage = False, - ) - -(options, filenames) = parser.parse_args() - -options.immutable = list(map(int,options.immutable)) - - -if filenames == []: filenames = [None] - -for name in filenames: - damask.util.report(scriptName,name) - - geom = damask.Grid.load(StringIO(''.join(sys.stdin.read())) if name is None else name) - - grid_original = geom.cells - damask.util.croak(geom) - material = np.tile(geom.material,np.where(grid_original == 1, 2,1)) # make one copy along dimensions with grid == 1 - grid = np.array(material.shape) - -# --- initialize support data --------------------------------------------------------------------- - -# store a copy of the initial material indices to find locations of immutable indices - material_original = np.copy(material) - - if not options.ndimage: - X,Y,Z = np.mgrid[0:grid[0],0:grid[1],0:grid[2]] - - # Calculates gaussian weights for simulating 3d diffusion - gauss = np.exp(-(X*X + Y*Y + Z*Z)/(2.0*options.d*options.d),dtype=np.float32) \ - /np.power(2.0*np.pi*options.d*options.d,(3.0 - np.count_nonzero(grid_original == 1))/2.,dtype=np.float32) - - gauss[:,:,:grid[2]//2:-1] = gauss[:,:,1:(grid[2]+1)//2] # trying to cope with uneven (odd) grid size - gauss[:,:grid[1]//2:-1,:] = gauss[:,1:(grid[1]+1)//2,:] - gauss[:grid[0]//2:-1,:,:] = gauss[1:(grid[0]+1)//2,:,:] - gauss = np.fft.rfftn(gauss).astype(np.complex64) - - for smoothIter in range(options.N): - - interfaceEnergy = np.zeros(material.shape,dtype=np.float32) - for i in (-1,0,1): - for j in (-1,0,1): - for k in (-1,0,1): - # assign interfacial energy to all voxels that have a differing neighbor (in Moore neighborhood) - interfaceEnergy = np.maximum(interfaceEnergy, - getInterfaceEnergy(material,np.roll(np.roll(np.roll( - material,i,axis=0), j,axis=1), k,axis=2))) - - # periodically extend interfacial energy array by half a grid size in positive and negative directions - periodic_interfaceEnergy = np.tile(interfaceEnergy,(3,3,3))[grid[0]//2:-grid[0]//2, - grid[1]//2:-grid[1]//2, - grid[2]//2:-grid[2]//2] - - # transform bulk volume (i.e. where interfacial energy remained zero), store index of closest boundary voxel - index = ndimage.morphology.distance_transform_edt(periodic_interfaceEnergy == 0., - return_distances = False, - return_indices = True) - - # want array index of nearest voxel on periodically extended boundary - periodic_bulkEnergy = periodic_interfaceEnergy[index[0], - index[1], - index[2]].reshape(2*grid) # fill bulk with energy of nearest interface - - if options.ndimage: - periodic_diffusedEnergy = ndimage.gaussian_filter( - np.where(ndimage.morphology.binary_dilation(periodic_interfaceEnergy > 0., - structure = struc, - iterations = int(round(options.d*2.))-1, # fat boundary - ), - periodic_bulkEnergy, # ...and zero everywhere else - 0.), - sigma = options.d) - else: - diffusedEnergy = np.fft.irfftn(np.fft.rfftn( - np.where( - ndimage.morphology.binary_dilation(interfaceEnergy > 0., - structure = struc, - iterations = int(round(options.d*2.))-1),# fat boundary - periodic_bulkEnergy[grid[0]//2:-grid[0]//2, # retain filled energy on fat boundary... - grid[1]//2:-grid[1]//2, - grid[2]//2:-grid[2]//2], # ...and zero everywhere else - 0.)).astype(np.complex64) * - gauss).astype(np.float32) - - periodic_diffusedEnergy = np.tile(diffusedEnergy,(3,3,3))[grid[0]//2:-grid[0]//2, - grid[1]//2:-grid[1]//2, - grid[2]//2:-grid[2]//2] # periodically extend the smoothed bulk energy - - - # transform voxels close to interface region - index = ndimage.morphology.distance_transform_edt(periodic_diffusedEnergy >= 0.95*np.amax(periodic_diffusedEnergy), - return_distances = False, - return_indices = True) # want index of closest bulk grain - - periodic_material = np.tile(material,(3,3,3))[grid[0]//2:-grid[0]//2, - grid[1]//2:-grid[1]//2, - grid[2]//2:-grid[2]//2] # periodically extend the geometry - - material = periodic_material[index[0], - index[1], - index[2]].reshape(2*grid)[grid[0]//2:-grid[0]//2, - grid[1]//2:-grid[1]//2, - grid[2]//2:-grid[2]//2] # extent grains into interface region - - # replace immutable materials with closest mutable ones - index = ndimage.morphology.distance_transform_edt(np.in1d(material,options.immutable).reshape(grid), - return_distances = False, - return_indices = True) - material = material[index[0], - index[1], - index[2]] - - immutable = np.zeros(material.shape, dtype=np.bool) - # find locations where immutable materials have been in original structure - for micro in options.immutable: - immutable += material_original == micro - - # undo any changes involving immutable materials - material = np.where(immutable, material_original,material) - - damask.Grid(material = material[0:grid_original[0],0:grid_original[1],0:grid_original[2]], - size = geom.size, - origin = geom.origin, - comments = geom.comments + [scriptID + ' ' + ' '.join(sys.argv[1:])], - )\ - .save(sys.stdout if name is None else name) From 6dcf6b972ccb06bcd0cb754e5d1a6705348b015e Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 2 Feb 2022 09:10:24 +0100 Subject: [PATCH 2/6] complaints from prospector (PEP8) --- python/tests/test_Orientation.py | 2 +- python/tests/test_Result.py | 16 ++++++++-------- python/tests/test_Rotation.py | 20 ++++++++++---------- python/tests/test_VTK.py | 16 ++++++++-------- 4 files changed, 27 insertions(+), 27 deletions(-) diff --git a/python/tests/test_Orientation.py b/python/tests/test_Orientation.py index 2a32adea4..bb4a14ed7 100644 --- a/python/tests/test_Orientation.py +++ b/python/tests/test_Orientation.py @@ -287,7 +287,7 @@ class TestOrientation: @pytest.mark.parametrize('family',crystal_families) @pytest.mark.parametrize('proper',[True,False]) def test_in_SST(self,family,proper): - assert Orientation(family=family).in_SST(np.zeros(3),proper) + assert Orientation(family=family).in_SST(np.zeros(3),proper) @pytest.mark.parametrize('function',['in_SST','IPF_color']) def test_invalid_argument(self,function): diff --git a/python/tests/test_Result.py b/python/tests/test_Result.py index 1da56dedf..71a4ce802 100644 --- a/python/tests/test_Result.py +++ b/python/tests/test_Result.py @@ -367,13 +367,13 @@ class TestResult: @pytest.mark.parametrize('mode',['cell','node']) def test_coordinates(self,default,mode): - if mode == 'cell': - a = grid_filters.coordinates0_point(default.cells,default.size,default.origin) - b = default.coordinates0_point.reshape(tuple(default.cells)+(3,),order='F') - elif mode == 'node': - a = grid_filters.coordinates0_node(default.cells,default.size,default.origin) - b = default.coordinates0_node.reshape(tuple(default.cells+1)+(3,),order='F') - assert np.allclose(a,b) + if mode == 'cell': + a = grid_filters.coordinates0_point(default.cells,default.size,default.origin) + b = default.coordinates0_point.reshape(tuple(default.cells)+(3,),order='F') + elif mode == 'node': + a = grid_filters.coordinates0_node(default.cells,default.size,default.origin) + b = default.coordinates0_node.reshape(tuple(default.cells+1)+(3,),order='F') + assert np.allclose(a,b) @pytest.mark.parametrize('output',['F','*',['P'],['P','F']],ids=range(4)) @pytest.mark.parametrize('fname',['12grains6x7x8_tensionY.hdf5'],ids=range(1)) @@ -421,7 +421,7 @@ class TestResult: def test_XDMF_datatypes(self,tmp_path,single_phase,update,ref_path): for shape in [('scalar',()),('vector',(3,)),('tensor',(3,3)),('matrix',(12,))]: for dtype in ['f4','f8','i1','i2','i4','i8','u1','u2','u4','u8']: - single_phase.add_calculation(f"np.ones(np.shape(#F#)[0:1]+{shape[1]},'{dtype}')",f'{shape[0]}_{dtype}') + single_phase.add_calculation(f"np.ones(np.shape(#F#)[0:1]+{shape[1]},'{dtype}')",f'{shape[0]}_{dtype}') fname = os.path.splitext(os.path.basename(single_phase.fname))[0]+'.xdmf' os.chdir(tmp_path) single_phase.export_XDMF() diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py index a431bc64b..fda986c2f 100644 --- a/python/tests/test_Rotation.py +++ b/python/tests/test_Rotation.py @@ -1076,19 +1076,19 @@ class TestRotation: def test_from_fiber_component(self,N,sigma): p = [] for run in range(5): - alpha = np.random.random()*2*np.pi,np.arccos(np.random.random()) - beta = np.random.random()*2*np.pi,np.arccos(np.random.random()) + alpha = np.random.random()*2*np.pi,np.arccos(np.random.random()) + beta = np.random.random()*2*np.pi,np.arccos(np.random.random()) - f_in_C = np.array([np.sin(alpha[0])*np.cos(alpha[1]), np.sin(alpha[0])*np.sin(alpha[1]), np.cos(alpha[0])]) - f_in_S = np.array([np.sin(beta[0] )*np.cos(beta[1] ), np.sin(beta[0] )*np.sin(beta[1] ), np.cos(beta[0] )]) - ax = np.append(np.cross(f_in_C,f_in_S), - np.arccos(np.dot(f_in_C,f_in_S))) - n = Rotation.from_axis_angle(ax if ax[3] > 0.0 else ax*-1.0 ,normalize=True) # rotation to align fiber axis in crystal and sample system + f_in_C = np.array([np.sin(alpha[0])*np.cos(alpha[1]), np.sin(alpha[0])*np.sin(alpha[1]), np.cos(alpha[0])]) + f_in_S = np.array([np.sin(beta[0] )*np.cos(beta[1] ), np.sin(beta[0] )*np.sin(beta[1] ), np.cos(beta[0] )]) + ax = np.append(np.cross(f_in_C,f_in_S), - np.arccos(np.dot(f_in_C,f_in_S))) + n = Rotation.from_axis_angle(ax if ax[3] > 0.0 else ax*-1.0 ,normalize=True) # rotation to align fiber axis in crystal and sample system - o = Rotation.from_fiber_component(alpha,beta,np.radians(sigma),N,False) - angles = np.arccos(np.clip(np.dot(o@np.broadcast_to(f_in_S,(N,3)),n@f_in_S),-1,1)) - dist = np.array(angles) * (np.random.randint(0,2,N)*2-1) + o = Rotation.from_fiber_component(alpha,beta,np.radians(sigma),N,False) + angles = np.arccos(np.clip(np.dot(o@np.broadcast_to(f_in_S,(N,3)),n@f_in_S),-1,1)) + dist = np.array(angles) * (np.random.randint(0,2,N)*2-1) - p.append(stats.normaltest(dist)[1]) + p.append(stats.normaltest(dist)[1]) sigma_out = np.degrees(np.std(dist)) p = np.average(p) diff --git a/python/tests/test_VTK.py b/python/tests/test_VTK.py index e59409a20..9f0dcc7cf 100644 --- a/python/tests/test_VTK.py +++ b/python/tests/test_VTK.py @@ -173,11 +173,11 @@ class TestVTK: polyData = VTK.from_poly_data(points) polyData.add(points,'coordinates') if update: - polyData.save(ref_path/'polyData') + polyData.save(ref_path/'polyData') else: - reference = VTK.load(ref_path/'polyData.vtp') - assert polyData.__repr__() == reference.__repr__() and \ - np.allclose(polyData.get('coordinates'),points) + reference = VTK.load(ref_path/'polyData.vtp') + assert polyData.__repr__() == reference.__repr__() and \ + np.allclose(polyData.get('coordinates'),points) @pytest.mark.xfail(int(vtk.vtkVersion.GetVTKVersion().split('.')[0])<8, reason='missing METADATA') def test_compare_reference_rectilinearGrid(self,update,ref_path,tmp_path): @@ -189,8 +189,8 @@ class TestVTK: rectilinearGrid.add(np.ascontiguousarray(c),'cell') rectilinearGrid.add(np.ascontiguousarray(n),'node') if update: - rectilinearGrid.save(ref_path/'rectilinearGrid') + rectilinearGrid.save(ref_path/'rectilinearGrid') else: - reference = VTK.load(ref_path/'rectilinearGrid.vtr') - assert rectilinearGrid.__repr__() == reference.__repr__() and \ - np.allclose(rectilinearGrid.get('cell'),c) + reference = VTK.load(ref_path/'rectilinearGrid.vtr') + assert rectilinearGrid.__repr__() == reference.__repr__() and \ + np.allclose(rectilinearGrid.get('cell'),c) From 4b2e104f0303735ccf2a89bd5c5811d3eafcfb54 Mon Sep 17 00:00:00 2001 From: Daniel Otto de Mentock Date: Wed, 2 Feb 2022 12:52:59 +0100 Subject: [PATCH 3/6] added new type for rng_seed in _typehints module --- python/damask/_typehints.py | 2 ++ python/damask/seeds.py | 7 ++++--- python/damask/util.py | 4 ++-- 3 files changed, 8 insertions(+), 5 deletions(-) diff --git a/python/damask/_typehints.py b/python/damask/_typehints.py index 67e920957..b8232a752 100644 --- a/python/damask/_typehints.py +++ b/python/damask/_typehints.py @@ -9,3 +9,5 @@ import numpy as np FloatSequence = Union[np.ndarray,Sequence[float]] IntSequence = Union[np.ndarray,Sequence[int]] FileHandle = Union[TextIO, str, Path] + +NumpyRngSeed = Union[int, IntSequence, np.random.SeedSequence, np.random.BitGenerator, np.random.Generator] diff --git a/python/damask/seeds.py b/python/damask/seeds.py index 47e34d871..943d1bfd9 100644 --- a/python/damask/seeds.py +++ b/python/damask/seeds.py @@ -1,3 +1,4 @@ + """Functionality for generation of seed points for Voronoi or Laguerre tessellation.""" from typing import Tuple as _Tuple @@ -5,7 +6,7 @@ from typing import Tuple as _Tuple from scipy import spatial as _spatial import numpy as _np -from ._typehints import FloatSequence as _FloatSequence, IntSequence as _IntSequence +from ._typehints import FloatSequence as _FloatSequence, IntSequence as _IntSequence, NumpyRngSeed from . import util as _util from . import grid_filters as _grid_filters @@ -13,7 +14,7 @@ from . import grid_filters as _grid_filters def from_random(size: _FloatSequence, N_seeds: int, cells: _IntSequence = None, - rng_seed=None) -> _np.ndarray: + rng_seed: NumpyRngSeed = None) -> _np.ndarray: """ Place seeds randomly in space. @@ -53,7 +54,7 @@ def from_Poisson_disc(size: _FloatSequence, N_candidates: int, distance: float, periodic: bool = True, - rng_seed=None) -> _np.ndarray: + rng_seed: NumpyRngSeed = None) -> _np.ndarray: """ Place seeds according to a Poisson disc distribution. diff --git a/python/damask/util.py b/python/damask/util.py index b4c287884..fc3f8ea8f 100644 --- a/python/damask/util.py +++ b/python/damask/util.py @@ -16,7 +16,7 @@ import numpy as np import h5py from . import version -from ._typehints import IntSequence, FloatSequence +from ._typehints import FloatSequence, NumpyRngSeed # limit visibility __all__=[ @@ -396,7 +396,7 @@ def execution_stamp(class_name: str, def hybrid_IA(dist: np.ndarray, N: int, - rng_seed: Union[int, IntSequence] = None) -> np.ndarray: + rng_seed: NumpyRngSeed = None) -> np.ndarray: """ Hybrid integer approximation. From 8fead8e306158ff8466fbd784e1d85ac4b8831d5 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 2 Feb 2022 13:28:17 +0100 Subject: [PATCH 4/6] relaxed tolerances for optimized code, the last digit might differ --- src/polynomials.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/polynomials.f90 b/src/polynomials.f90 index 49a46e2e4..155f54bf8 100644 --- a/src/polynomials.f90 +++ b/src/polynomials.f90 @@ -163,7 +163,7 @@ subroutine selfTest 'T_ref: '//trim(adjustl(x_ref_s))//IO_EOL Dict => YAML_parse_str(trim(YAML_s)) p2 = polynomial(dict%asDict(),'C','T') - if (dNeq(p1%at(x),p2%at(x),1.0e-12_pReal)) error stop 'polynomials: init' + if (dNeq(p1%at(x),p2%at(x),1.0e-10_pReal)) error stop 'polynomials: init' p1 = polynomial(coef*[0.0_pReal,1.0_pReal,0.0_pReal],x_ref) if (dNeq(p1%at(x_ref+x),-p1%at(x_ref-x),1.0e-10_pReal)) error stop 'polynomials: eval(odd)' From d868a240b242948a88f414b2b6a56d6c694ef49f Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 2 Feb 2022 16:45:13 +0000 Subject: [PATCH 5/6] bugix: change of behavior --- PRIVATE | 2 +- ...hase_mechanical_eigen_thermalexpansion.f90 | 45 +++++++------------ src/phase_mechanical_elastic.f90 | 4 +- 3 files changed, 20 insertions(+), 31 deletions(-) diff --git a/PRIVATE b/PRIVATE index 5598ec489..419ba7b6f 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 5598ec4892b0921fccf63e8570f9fcd8e0365716 +Subproject commit 419ba7b6f6989a14e0c91c1f1f6621daf8215b71 diff --git a/src/phase_mechanical_eigen_thermalexpansion.f90 b/src/phase_mechanical_eigen_thermalexpansion.f90 index 70bf84cd8..a5d9868a8 100644 --- a/src/phase_mechanical_eigen_thermalexpansion.f90 +++ b/src/phase_mechanical_eigen_thermalexpansion.f90 @@ -8,10 +8,9 @@ submodule(phase:eigen) thermalexpansion integer, dimension(:), allocatable :: kinematics_thermal_expansion_instance type :: tParameters - real(pReal) :: & - T_ref - real(pReal), dimension(3,3,3) :: & - A = 0.0_pReal + type(tPolynomial) :: & + A_11, & + A_33 end type tParameters type(tParameters), dimension(:), allocatable :: param @@ -34,7 +33,7 @@ module function thermalexpansion_init(kinematics_length) result(myKinematics) phase, & mech, & kinematics, & - kinematic_type + myConfig print'(/,1x,a)', '<<<+- phase:mechanical:eigen:thermalexpansion init -+>>>' @@ -56,21 +55,13 @@ module function thermalexpansion_init(kinematics_length) result(myKinematics) do k = 1, kinematics%length if (myKinematics(k,p)) then associate(prm => param(kinematics_thermal_expansion_instance(p))) - kinematic_type => kinematics%get(k) - prm%T_ref = kinematic_type%get_asFloat('T_ref', defaultVal=T_ROOM) + myConfig => kinematics%get(k) + + prm%A_11 = polynomial(myConfig%asDict(),'A_11','T') + if (any(phase_lattice(p) == ['hP','tI'])) & + prm%A_33 = polynomial(myConfig%asDict(),'A_33','T') - 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,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,3) = kinematic_type%get_asFloat('A_33,T^2',defaultVal=0.0_pReal) - end if - do i=1, size(prm%A,3) - prm%A(1:3,1:3,i) = lattice_symmetrize_33(prm%A(1:3,1:3,i),phase_lattice(p)) - end do end associate end if end do @@ -91,22 +82,20 @@ module subroutine thermalexpansion_LiAndItsTangent(Li, dLi_dTstar, ph,me) dLi_dTstar !< derivative of Li with respect to Tstar (4th-order tensor defined to be zero) real(pReal) :: T, dot_T + real(pReal), dimension(3,3) :: A T = thermal_T(ph,me) dot_T = thermal_dot_T(ph,me) associate(prm => param(kinematics_thermal_expansion_instance(ph))) - Li = dot_T * ( & - prm%A(1:3,1:3,1) & ! constant coefficient - + prm%A(1:3,1:3,2)*(T - prm%T_ref) & ! 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.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 & - ) + + A = 0.0_pReal + A(1,1) = prm%A_11%at(T) + if (any(phase_lattice(ph) == ['hP','tI'])) A(3,3) = prm%A_33%at(T) + A = lattice_symmetrize_33(A,phase_lattice(ph)) + Li = dot_T * A + end associate dLi_dTstar = 0.0_pReal diff --git a/src/phase_mechanical_elastic.f90 b/src/phase_mechanical_elastic.f90 index c57849b1f..084cc75bc 100644 --- a/src/phase_mechanical_elastic.f90 +++ b/src/phase_mechanical_elastic.f90 @@ -48,7 +48,7 @@ module subroutine elastic_init(phases) prm%C_11 = polynomial(elastic%asDict(),'C_11','T') prm%C_12 = polynomial(elastic%asDict(),'C_12','T') prm%C_44 = polynomial(elastic%asDict(),'C_44','T') - + if (any(phase_lattice(ph) == ['hP','tI'])) then prm%C_13 = polynomial(elastic%asDict(),'C_13','T') prm%C_33 = polynomial(elastic%asDict(),'C_33','T') @@ -77,7 +77,7 @@ pure module function elastic_C66(ph,en) result(C66) associate(prm => param(ph)) - + C66 = 0.0_pReal T = thermal_T(ph,en) From 7078b5ec8709ad78dff64813d0e290c583c13979 Mon Sep 17 00:00:00 2001 From: Test User Date: Wed, 2 Feb 2022 22:07:12 +0100 Subject: [PATCH 6/6] [skip ci] updated version information after successful test of v3.0.0-alpha5-552-ga6e78c5b6 --- python/damask/VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/damask/VERSION b/python/damask/VERSION index 8d2c8c56d..0f9864e69 100644 --- a/python/damask/VERSION +++ b/python/damask/VERSION @@ -1 +1 @@ -v3.0.0-alpha5-545-gad74f5dbe +v3.0.0-alpha5-552-ga6e78c5b6