diff --git a/VERSION b/VERSION index 6f8955f5c..29f305cbe 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.3-2402-gb88f5ec0 +v2.0.3-2464-g90f93d23 diff --git a/python/damask/_Lambert.py b/python/damask/_Lambert.py deleted file mode 100644 index a823764e9..000000000 --- a/python/damask/_Lambert.py +++ /dev/null @@ -1,164 +0,0 @@ -#################################################################################################### -# Code below available according to the following conditions on -# https://github.com/MarDiehl/3Drotations -#################################################################################################### -# Copyright (c) 2017-2019, Martin Diehl/Max-Planck-Institut für Eisenforschung GmbH -# Copyright (c) 2013-2014, Marc De Graef/Carnegie Mellon University -# All rights reserved. -# -# Redistribution and use in source and binary forms, with or without modification, are -# permitted provided that the following conditions are met: -# -# - Redistributions of source code must retain the above copyright notice, this list -# of conditions and the following disclaimer. -# - Redistributions in binary form must reproduce the above copyright notice, this -# list of conditions and the following disclaimer in the documentation and/or -# other materials provided with the distribution. -# - Neither the names of Marc De Graef, Carnegie Mellon University nor the names -# of its contributors may be used to endorse or promote products derived from -# this software without specific prior written permission. -# -# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" -# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE -# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE -# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE -# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL -# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR -# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER -# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, -# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE -# USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -#################################################################################################### - -import numpy as np - -sc = np.pi**(1./6.)/6.**(1./6.) -beta = np.pi**(5./6.)/6.**(1./6.)/2. -R1 = (3.*np.pi/4.)**(1./3.) - -def cube_to_ball(cube): - """ - Map a point in a uniform refinable cubical grid to a point on a uniform refinable grid on a ball. - - Parameters - ---------- - cube : numpy.ndarray - coordinates of a point in a uniform refinable cubical grid. - - References - ---------- - D. Roşca et al., Modelling and Simulation in Materials Science and Engineering 22:075013, 2014 - https://doi.org/10.1088/0965-0393/22/7/075013 - - """ - cube_ = np.clip(cube,None,np.pi**(2./3.) * 0.5) if np.isclose(np.abs(np.max(cube)),np.pi**(2./3.) * 0.5,atol=1e-6) else cube - - # transform to the sphere grid via the curved square, and intercept the zero point - if np.allclose(cube_,0.0,rtol=0.0,atol=1.0e-16): - ball = np.zeros(3) - else: - # get pyramide and scale by grid parameter ratio - p = _get_order(cube_) - XYZ = cube_[p[0]] * sc - - # intercept all the points along the z-axis - if np.allclose(XYZ[0:2],0.0,rtol=0.0,atol=1.0e-16): - ball = np.array([0.0, 0.0, np.sqrt(6.0/np.pi) * XYZ[2]]) - else: - order = [1,0] if np.abs(XYZ[1]) <= np.abs(XYZ[0]) else [0,1] - q = np.pi/12.0 * XYZ[order[0]]/XYZ[order[1]] - c = np.cos(q) - s = np.sin(q) - q = R1*2.0**0.25/beta * XYZ[order[1]] / np.sqrt(np.sqrt(2.0)-c) - T = np.array([ (np.sqrt(2.0)*c - 1.0), np.sqrt(2.0) * s]) * q - - # transform to sphere grid (inverse Lambert) - # note that there is no need to worry about dividing by zero, since XYZ[2] can not become zero - c = np.sum(T**2) - s = c * np.pi/24.0 /XYZ[2]**2 - c = c * np.sqrt(np.pi/24.0)/XYZ[2] - q = np.sqrt( 1.0 - s ) - ball = np.array([ T[order[1]] * q, T[order[0]] * q, np.sqrt(6.0/np.pi) * XYZ[2] - c ]) - - # reverse the coordinates back to the regular order according to the original pyramid number - ball = ball[p[1]] - - return ball - - -def ball_to_cube(ball): - """ - Map a point on a uniform refinable grid on a ball to a point in a uniform refinable cubical grid. - - Parameters - ---------- - ball : numpy.ndarray - coordinates of a point on a uniform refinable grid on a ball. - - References - ---------- - D. Roşca et al., Modelling and Simulation in Materials Science and Engineering 22:075013, 2014 - https://doi.org/10.1088/0965-0393/22/7/075013 - - """ - ball_ = ball/np.linalg.norm(ball)*R1 if np.isclose(np.linalg.norm(ball),R1,atol=1e-6) else ball - rs = np.linalg.norm(ball_) - - if np.allclose(ball_,0.0,rtol=0.0,atol=1.0e-16): - cube = np.zeros(3) - else: - p = _get_order(ball_) - xyz3 = ball_[p[0]] - - # inverse M_3 - xyz2 = xyz3[0:2] * np.sqrt( 2.0*rs/(rs+np.abs(xyz3[2])) ) - - # inverse M_2 - qxy = np.sum(xyz2**2) - - if np.isclose(qxy,0.0,rtol=0.0,atol=1.0e-16): - Tinv = np.zeros(2) - else: - q2 = qxy + np.max(np.abs(xyz2))**2 - sq2 = np.sqrt(q2) - q = (beta/np.sqrt(2.0)/R1) * np.sqrt(q2*qxy/(q2-np.max(np.abs(xyz2))*sq2)) - tt = np.clip((np.min(np.abs(xyz2))**2+np.max(np.abs(xyz2))*sq2)/np.sqrt(2.0)/qxy,-1.0,1.0) - Tinv = np.array([1.0,np.arccos(tt)/np.pi*12.0]) if np.abs(xyz2[1]) <= np.abs(xyz2[0]) else \ - np.array([np.arccos(tt)/np.pi*12.0,1.0]) - Tinv = q * np.where(xyz2<0.0,-Tinv,Tinv) - - # inverse M_1 - cube = np.array([ Tinv[0], Tinv[1], (-1.0 if xyz3[2] < 0.0 else 1.0) * rs / np.sqrt(6.0/np.pi) ]) /sc - # reverse the coordinates back to the regular order according to the original pyramid number - cube = cube[p[1]] - - return cube - - -def _get_order(xyz): - """ - Get order of the coordinates. - - Depending on the pyramid in which the point is located, the order need to be adjusted. - - Parameters - ---------- - xyz : numpy.ndarray - coordinates of a point on a uniform refinable grid on a ball or - in a uniform refinable cubical grid. - - References - ---------- - D. Roşca et al., Modelling and Simulation in Materials Science and Engineering 22:075013, 2014 - https://doi.org/10.1088/0965-0393/22/7/075013 - - """ - if (abs(xyz[0])<= xyz[2]) and (abs(xyz[1])<= xyz[2]) or \ - (abs(xyz[0])<=-xyz[2]) and (abs(xyz[1])<=-xyz[2]): - return [[0,1,2],[0,1,2]] - elif (abs(xyz[2])<= xyz[0]) and (abs(xyz[1])<= xyz[0]) or \ - (abs(xyz[2])<=-xyz[0]) and (abs(xyz[1])<=-xyz[0]): - return [[1,2,0],[2,0,1]] - elif (abs(xyz[0])<= xyz[1]) and (abs(xyz[2])<= xyz[1]) or \ - (abs(xyz[0])<=-xyz[1]) and (abs(xyz[2])<=-xyz[1]): - return [[2,0,1],[1,2,0]] diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index cb5b2fbd1..059f148a4 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -1,10 +1,14 @@ import numpy as np -from ._Lambert import ball_to_cube, cube_to_ball from . import mechanics _P = -1 +# parameters for conversion from/to cubochoric +_sc = np.pi**(1./6.)/6.**(1./6.) +_beta = np.pi**(5./6.)/6.**(1./6.)/2. +_R1 = (3.*np.pi/4.)**(1./3.) + def iszero(a): return np.isclose(a,0.0,atol=1.0e-12,rtol=0.0) @@ -288,7 +292,7 @@ class Rotation: """Homochoric vector: (h_1, h_2, h_3).""" return Rotation.qu2ho(self.quaternion) - def asCubochoric(self): + def as_cubochoric(self): """Cubochoric vector: (c_1, c_2, c_3).""" return Rotation.qu2cu(self.quaternion) @@ -312,6 +316,7 @@ class Rotation: asMatrix = as_matrix asRodrigues = as_Rodrigues asHomochoric = as_homochoric + asCubochoric = as_cubochoric ################################################################################################ # Static constructors. The input data needs to follow the conventions, options allow to @@ -433,7 +438,7 @@ class Rotation: return Rotation(Rotation.ho2qu(ho)) @staticmethod - def fromCubochoric(cubochoric, + def from_cubochoric(cubochoric, P = -1): cu = np.array(cubochoric,dtype=float) @@ -508,6 +513,7 @@ class Rotation: fromMatrix = from_matrix fromRodrigues = from_Rodrigues fromHomochoric = from_homochoric + fromCubochoric = from_cubochoric fromRandom = from_random #################################################################################################### @@ -1098,12 +1104,71 @@ class Rotation: @staticmethod def ho2cu(ho): - """Homochoric vector to cubochoric vector.""" - if len(ho.shape) == 1: - return ball_to_cube(ho) - else: - raise NotImplementedError('Support for multiple rotations missing') + """ + Homochoric vector to cubochoric vector. + References + ---------- + D. Roşca et al., Modelling and Simulation in Materials Science and Engineering 22:075013, 2014 + https://doi.org/10.1088/0965-0393/22/7/075013 + + """ + if len(ho.shape) == 1: + rs = np.linalg.norm(ho) + + if np.allclose(ho,0.0,rtol=0.0,atol=1.0e-16): + cu = np.zeros(3) + else: + xyz3 = ho[Rotation._get_pyramid_order(ho,'forward')] + + # inverse M_3 + xyz2 = xyz3[0:2] * np.sqrt( 2.0*rs/(rs+np.abs(xyz3[2])) ) + + # inverse M_2 + qxy = np.sum(xyz2**2) + + if np.isclose(qxy,0.0,rtol=0.0,atol=1.0e-16): + Tinv = np.zeros(2) + else: + q2 = qxy + np.max(np.abs(xyz2))**2 + sq2 = np.sqrt(q2) + q = (_beta/np.sqrt(2.0)/_R1) * np.sqrt(q2*qxy/(q2-np.max(np.abs(xyz2))*sq2)) + tt = np.clip((np.min(np.abs(xyz2))**2+np.max(np.abs(xyz2))*sq2)/np.sqrt(2.0)/qxy,-1.0,1.0) + Tinv = np.array([1.0,np.arccos(tt)/np.pi*12.0]) if np.abs(xyz2[1]) <= np.abs(xyz2[0]) else \ + np.array([np.arccos(tt)/np.pi*12.0,1.0]) + Tinv = q * np.where(xyz2<0.0,-Tinv,Tinv) + + # inverse M_1 + cu = np.array([ Tinv[0], Tinv[1], (-1.0 if xyz3[2] < 0.0 else 1.0) * rs / np.sqrt(6.0/np.pi) ]) /_sc + cu = cu[Rotation._get_pyramid_order(ho,'backward')] + else: + rs = np.linalg.norm(ho,axis=-1,keepdims=True) + + xyz3 = np.take_along_axis(ho,Rotation._get_pyramid_order(ho,'forward'),-1) + + with np.errstate(invalid='ignore',divide='ignore'): + # inverse M_3 + xyz2 = xyz3[...,0:2] * np.sqrt( 2.0*rs/(rs+np.abs(xyz3[...,2:3])) ) + qxy = np.sum(xyz2**2,axis=-1,keepdims=True) + + q2 = qxy + np.max(np.abs(xyz2),axis=-1,keepdims=True)**2 + sq2 = np.sqrt(q2) + q = (_beta/np.sqrt(2.0)/_R1) * np.sqrt(q2*qxy/(q2-np.max(np.abs(xyz2),axis=-1,keepdims=True)*sq2)) + tt = np.clip((np.min(np.abs(xyz2),axis=-1,keepdims=True)**2\ + +np.max(np.abs(xyz2),axis=-1,keepdims=True)*sq2)/np.sqrt(2.0)/qxy,-1.0,1.0) + T_inv = np.where(np.abs(xyz2[...,1:2]) <= np.abs(xyz2[...,0:1]), + np.block([np.ones_like(tt),np.arccos(tt)/np.pi*12.0]), + np.block([np.arccos(tt)/np.pi*12.0,np.ones_like(tt)]))*q + T_inv[xyz2<0.0] *= -1.0 + T_inv[np.broadcast_to(np.isclose(qxy,0.0,rtol=0.0,atol=1.0e-12),T_inv.shape)] = 0.0 + cu = np.block([T_inv, np.where(xyz3[...,2:3]<0.0,-np.ones_like(xyz3[...,2:3]),np.ones_like(xyz3[...,2:3])) \ + * rs/np.sqrt(6.0/np.pi), + ])/ _sc + + cu[np.isclose(np.sum(np.abs(ho),axis=-1),0.0,rtol=0.0,atol=1.0e-16)] = 0.0 + cu = np.take_along_axis(cu,Rotation._get_pyramid_order(ho,'backward'),-1) + + return cu #---------- Cubochoric ---------- @staticmethod @@ -1133,8 +1198,110 @@ class Rotation: @staticmethod def cu2ho(cu): - """Cubochoric vector to homochoric vector.""" + """ + Cubochoric vector to homochoric vector. + + References + ---------- + D. Roşca et al., Modelling and Simulation in Materials Science and Engineering 22:075013, 2014 + https://doi.org/10.1088/0965-0393/22/7/075013 + + """ if len(cu.shape) == 1: - return cube_to_ball(cu) + # transform to the sphere grid via the curved square, and intercept the zero point + if np.allclose(cu,0.0,rtol=0.0,atol=1.0e-16): + ho = np.zeros(3) + else: + # get pyramide and scale by grid parameter ratio + XYZ = cu[Rotation._get_pyramid_order(cu,'forward')] * _sc + + # intercept all the points along the z-axis + if np.allclose(XYZ[0:2],0.0,rtol=0.0,atol=1.0e-16): + ho = np.array([0.0, 0.0, np.sqrt(6.0/np.pi) * XYZ[2]]) + else: + order = [1,0] if np.abs(XYZ[1]) <= np.abs(XYZ[0]) else [0,1] + q = np.pi/12.0 * XYZ[order[0]]/XYZ[order[1]] + c = np.cos(q) + s = np.sin(q) + q = _R1*2.0**0.25/_beta * XYZ[order[1]] / np.sqrt(np.sqrt(2.0)-c) + T = np.array([ (np.sqrt(2.0)*c - 1.0), np.sqrt(2.0) * s]) * q + + # transform to sphere grid (inverse Lambert) + # note that there is no need to worry about dividing by zero, since XYZ[2] can not become zero + c = np.sum(T**2) + s = c * np.pi/24.0 /XYZ[2]**2 + c = c * np.sqrt(np.pi/24.0)/XYZ[2] + + q = np.sqrt( 1.0 - s ) + ho = np.array([ T[order[1]] * q, T[order[0]] * q, np.sqrt(6.0/np.pi) * XYZ[2] - c ]) + + ho = ho[Rotation._get_pyramid_order(cu,'backward')] else: - raise NotImplementedError('Support for multiple rotations missing') + with np.errstate(invalid='ignore',divide='ignore'): + # get pyramide and scale by grid parameter ratio + XYZ = np.take_along_axis(cu,Rotation._get_pyramid_order(cu,'forward'),-1) * _sc + order = np.abs(XYZ[...,1:2]) <= np.abs(XYZ[...,0:1]) + q = np.pi/12.0 * np.where(order,XYZ[...,1:2],XYZ[...,0:1]) \ + / np.where(order,XYZ[...,0:1],XYZ[...,1:2]) + c = np.cos(q) + s = np.sin(q) + q = _R1*2.0**0.25/_beta/ np.sqrt(np.sqrt(2.0)-c) \ + * np.where(order,XYZ[...,0:1],XYZ[...,1:2]) + + T = np.block([ (np.sqrt(2.0)*c - 1.0), np.sqrt(2.0) * s]) * q + + # transform to sphere grid (inverse Lambert) + c = np.sum(T**2,axis=-1,keepdims=True) + s = c * np.pi/24.0 /XYZ[...,2:3]**2 + c = c * np.sqrt(np.pi/24.0)/XYZ[...,2:3] + q = np.sqrt( 1.0 - s) + + ho = np.where(np.isclose(np.sum(np.abs(XYZ[...,0:2]),axis=-1,keepdims=True),0.0,rtol=0.0,atol=1.0e-16), + np.block([np.zeros_like(XYZ[...,0:2]),np.sqrt(6.0/np.pi) *XYZ[...,2:3]]), + np.block([np.where(order,T[...,0:1],T[...,1:2])*q, + np.where(order,T[...,1:2],T[...,0:1])*q, + np.sqrt(6.0/np.pi) * XYZ[...,2:3] - c]) + ) + + ho[np.isclose(np.sum(np.abs(cu),axis=-1),0.0,rtol=0.0,atol=1.0e-16)] = 0.0 + ho = np.take_along_axis(ho,Rotation._get_pyramid_order(cu,'backward'),-1) + + return ho + + + @staticmethod + def _get_pyramid_order(xyz,direction=None): + """ + Get order of the coordinates. + + Depending on the pyramid in which the point is located, the order need to be adjusted. + + Parameters + ---------- + xyz : numpy.ndarray + coordinates of a point on a uniform refinable grid on a ball or + in a uniform refinable cubical grid. + + References + ---------- + D. Roşca et al., Modelling and Simulation in Materials Science and Engineering 22:075013, 2014 + https://doi.org/10.1088/0965-0393/22/7/075013 + + """ + order = {'forward':np.array([[0,1,2],[1,2,0],[2,0,1]]), + 'backward':np.array([[0,1,2],[2,0,1],[1,2,0]])} + if len(xyz.shape) == 1: + if np.maximum(abs(xyz[0]),abs(xyz[1])) <= xyz[2] or \ + np.maximum(abs(xyz[0]),abs(xyz[1])) <=-xyz[2]: + p = 0 + elif np.maximum(abs(xyz[1]),abs(xyz[2])) <= xyz[0] or \ + np.maximum(abs(xyz[1]),abs(xyz[2])) <=-xyz[0]: + p = 1 + elif np.maximum(abs(xyz[2]),abs(xyz[0])) <= xyz[1] or \ + np.maximum(abs(xyz[2]),abs(xyz[0])) <=-xyz[1]: + p = 2 + else: + p = np.where(np.maximum(np.abs(xyz[...,0]),np.abs(xyz[...,1])) <= np.abs(xyz[...,2]),0, + np.where(np.maximum(np.abs(xyz[...,1]),np.abs(xyz[...,2])) <= np.abs(xyz[...,0]),1,2)) + + return order[direction][p] diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py index 302e895ee..b7442035f 100644 --- a/python/tests/test_Rotation.py +++ b/python/tests/test_Rotation.py @@ -78,9 +78,9 @@ def default(): specials_scatter /= np.linalg.norm(specials_scatter,axis=1).reshape(-1,1) specials_scatter[specials_scatter[:,0]<0]*=-1 - return [Rotation.fromQuaternion(s) for s in specials] + \ - [Rotation.fromQuaternion(s) for s in specials_scatter] + \ - [Rotation.fromRandom() for _ in range(n-len(specials)-len(specials_scatter))] + return [Rotation.from_quaternion(s) for s in specials] + \ + [Rotation.from_quaternion(s) for s in specials_scatter] + \ + [Rotation.from_random() for _ in range(n-len(specials)-len(specials_scatter))] @pytest.fixture def reference_dir(reference_dir_base): @@ -92,41 +92,41 @@ class TestRotation: def test_Eulers(self,default): for rot in default: - m = rot.asQuaternion() - o = Rotation.fromEulers(rot.asEulers()).asQuaternion() + m = rot.as_quaternion() + o = Rotation.from_Eulers(rot.as_Eulers()).as_quaternion() ok = np.allclose(m,o,atol=atol) - if np.isclose(rot.asQuaternion()[0],0.0,atol=atol): + if np.isclose(rot.as_quaternion()[0],0.0,atol=atol): ok = ok or np.allclose(m*-1.,o,atol=atol) - print(m,o,rot.asQuaternion()) + print(m,o,rot.as_quaternion()) assert ok and np.isclose(np.linalg.norm(o),1.0) def test_AxisAngle(self,default): for rot in default: - m = rot.asEulers() - o = Rotation.fromAxisAngle(rot.asAxisAngle()).asEulers() + m = rot.as_Eulers() + o = Rotation.from_axis_angle(rot.as_axis_angle()).as_Eulers() u = np.array([np.pi*2,np.pi,np.pi*2]) ok = np.allclose(m,o,atol=atol) ok = ok or np.allclose(np.where(np.isclose(m,u),m-u,m),np.where(np.isclose(o,u),o-u,o),atol=atol) if np.isclose(m[1],0.0,atol=atol) or np.isclose(m[1],np.pi,atol=atol): sum_phi = np.unwrap([m[0]+m[2],o[0]+o[2]]) ok = ok or np.isclose(sum_phi[0],sum_phi[1],atol=atol) - print(m,o,rot.asQuaternion()) + print(m,o,rot.as_quaternion()) assert ok and (np.zeros(3)-1.e-9 <= o).all() and (o <= np.array([np.pi*2.,np.pi,np.pi*2.])+1.e-9).all() def test_Matrix(self,default): for rot in default: - m = rot.asAxisAngle() - o = Rotation.fromAxisAngle(rot.asAxisAngle()).asAxisAngle() + m = rot.as_axis_angle() + o = Rotation.from_axis_angle(rot.as_axis_angle()).as_axis_angle() ok = np.allclose(m,o,atol=atol) if np.isclose(m[3],np.pi,atol=atol): ok = ok or np.allclose(m*np.array([-1.,-1.,-1.,1.]),o,atol=atol) - print(m,o,rot.asQuaternion()) + print(m,o,rot.as_quaternion()) assert ok and np.isclose(np.linalg.norm(o[:3]),1.0) and o[3]<=np.pi++1.e-9 def test_Rodrigues(self,default): for rot in default: - m = rot.asMatrix() - o = Rotation.fromRodrigues(rot.asRodrigues()).asMatrix() + m = rot.as_matrix() + o = Rotation.from_Rodrigues(rot.as_Rodrigues()).as_matrix() ok = np.allclose(m,o,atol=atol) print(m,o) assert ok and np.isclose(np.linalg.det(o),1.0) @@ -134,27 +134,27 @@ class TestRotation: def test_Homochoric(self,default): cutoff = np.tan(np.pi*.5*(1.-1e-4)) for rot in default: - m = rot.asRodrigues() - o = Rotation.fromHomochoric(rot.asHomochoric()).asRodrigues() + m = rot.as_Rodrigues() + o = Rotation.from_homochoric(rot.as_homochoric()).as_Rodrigues() ok = np.allclose(np.clip(m,None,cutoff),np.clip(o,None,cutoff),atol=atol) ok = ok or np.isclose(m[3],0.0,atol=atol) - print(m,o,rot.asQuaternion()) + print(m,o,rot.as_quaternion()) assert ok and np.isclose(np.linalg.norm(o[:3]),1.0) def test_Cubochoric(self,default): for rot in default: - m = rot.asHomochoric() - o = Rotation.fromCubochoric(rot.asCubochoric()).asHomochoric() + m = rot.as_homochoric() + o = Rotation.from_cubochoric(rot.as_cubochoric()).as_homochoric() ok = np.allclose(m,o,atol=atol) - print(m,o,rot.asQuaternion()) + print(m,o,rot.as_quaternion()) assert ok and np.linalg.norm(o) < (3.*np.pi/4.)**(1./3.) + 1.e-9 def test_Quaternion(self,default): for rot in default: - m = rot.asCubochoric() - o = Rotation.fromQuaternion(rot.asQuaternion()).asCubochoric() + m = rot.as_cubochoric() + o = Rotation.from_quaternion(rot.as_quaternion()).as_cubochoric() ok = np.allclose(m,o,atol=atol) - print(m,o,rot.asQuaternion()) + print(m,o,rot.as_quaternion()) assert ok and o.max() < np.pi**(2./3.)*0.5+1.e-9 @pytest.mark.parametrize('function',[Rotation.from_quaternion, @@ -185,9 +185,11 @@ class TestRotation: Rotation.qu2eu, Rotation.qu2ax, Rotation.qu2ro, - Rotation.qu2ho]) + Rotation.qu2ho, + Rotation.qu2cu + ]) def test_quaternion_vectorization(self,default,conversion): - qu = np.array([rot.asQuaternion() for rot in default]) + qu = np.array([rot.as_quaternion() for rot in default]) conversion(qu.reshape(qu.shape[0]//2,-1,4)) co = conversion(qu) for q,c in zip(qu,co): @@ -199,9 +201,10 @@ class TestRotation: Rotation.om2ax, Rotation.om2ro, Rotation.om2ho, + Rotation.om2cu ]) def test_matrix_vectorization(self,default,conversion): - om = np.array([rot.asMatrix() for rot in default]) + om = np.array([rot.as_matrix() for rot in default]) conversion(om.reshape(om.shape[0]//2,-1,3,3)) co = conversion(om) for o,c in zip(om,co): @@ -213,9 +216,10 @@ class TestRotation: Rotation.eu2ax, Rotation.eu2ro, Rotation.eu2ho, + Rotation.eu2cu ]) def test_Euler_vectorization(self,default,conversion): - eu = np.array([rot.asEulers() for rot in default]) + eu = np.array([rot.as_Eulers() for rot in default]) conversion(eu.reshape(eu.shape[0]//2,-1,3)) co = conversion(eu) for e,c in zip(eu,co): @@ -227,9 +231,10 @@ class TestRotation: Rotation.ax2eu, Rotation.ax2ro, Rotation.ax2ho, + Rotation.ax2cu ]) def test_axisAngle_vectorization(self,default,conversion): - ax = np.array([rot.asAxisAngle() for rot in default]) + ax = np.array([rot.as_axis_angle() for rot in default]) conversion(ax.reshape(ax.shape[0]//2,-1,4)) co = conversion(ax) for a,c in zip(ax,co): @@ -242,9 +247,10 @@ class TestRotation: Rotation.ro2eu, Rotation.ro2ax, Rotation.ro2ho, + Rotation.ro2cu ]) def test_Rodrigues_vectorization(self,default,conversion): - ro = np.array([rot.asRodrigues() for rot in default]) + ro = np.array([rot.as_Rodrigues() for rot in default]) conversion(ro.reshape(ro.shape[0]//2,-1,4)) co = conversion(ro) for r,c in zip(ro,co): @@ -256,11 +262,41 @@ class TestRotation: Rotation.ho2eu, Rotation.ho2ax, Rotation.ho2ro, + Rotation.ho2cu ]) def test_homochoric_vectorization(self,default,conversion): - ho = np.array([rot.asHomochoric() for rot in default]) + ho = np.array([rot.as_homochoric() for rot in default]) conversion(ho.reshape(ho.shape[0]//2,-1,3)) co = conversion(ho) for h,c in zip(ho,co): print(h,c) assert np.allclose(conversion(h),c) + + @pytest.mark.parametrize('conversion',[Rotation.cu2qu, + Rotation.cu2om, + Rotation.cu2eu, + Rotation.cu2ax, + Rotation.cu2ro, + Rotation.cu2ho + ]) + def test_cubochoric_vectorization(self,default,conversion): + cu = np.array([rot.as_cubochoric() for rot in default]) + conversion(cu.reshape(cu.shape[0]//2,-1,3)) + co = conversion(cu) + for u,c in zip(cu,co): + print(u,c) + assert np.allclose(conversion(u),c) + + @pytest.mark.parametrize('direction',['forward', + 'backward']) + def test_pyramid_vectorization(self,direction): + p = np.random.rand(n,3) + o = Rotation._get_pyramid_order(p,direction) + for i,o_i in enumerate(o): + assert np.all(o_i==Rotation._get_pyramid_order(p[i],direction)) + + def test_pyramid_invariant(self): + a = np.random.rand(n,3) + f = Rotation._get_pyramid_order(a,'forward') + b = Rotation._get_pyramid_order(a,'backward') + assert np.all(np.take_along_axis(np.take_along_axis(a,f,-1),b,-1) == a) diff --git a/python/tests/test_grid_filters.py b/python/tests/test_grid_filters.py index eb359006a..ab60c0446 100644 --- a/python/tests/test_grid_filters.py +++ b/python/tests/test_grid_filters.py @@ -101,6 +101,7 @@ class TestGridFilters: with pytest.raises(ValueError): function(uneven) + @pytest.mark.parametrize('mode',[True,False]) @pytest.mark.parametrize('function',[grid_filters.node_coord0_gridSizeOrigin, grid_filters.cell_coord0_gridSizeOrigin]) @@ -120,3 +121,186 @@ class TestGridFilters: grid = np.random.randint(8,32,(3)) F = np.broadcast_to(np.eye(3), tuple(grid)+(3,3)) assert all(grid_filters.regrid(size,F,grid) == np.arange(grid.prod())) + + + @pytest.mark.parametrize('differential_operator',[grid_filters.curl, + grid_filters.divergence, + grid_filters.gradient]) + def test_differential_operator_constant(self,differential_operator): + size = np.random.random(3)+1.0 + grid = np.random.randint(8,32,(3)) + shapes = { + grid_filters.curl: [(3,),(3,3)], + grid_filters.divergence:[(3,),(3,3)], + grid_filters.gradient: [(1,),(3,)] + } + for shape in shapes[differential_operator]: + field = np.ones(tuple(grid)+shape)*np.random.random()*1.0e5 + assert np.allclose(differential_operator(size,field),0.0) + + + grad_test_data = [ + (['np.sin(np.pi*2*nodes[...,0]/size[0])', '0.0', '0.0'], + ['np.cos(np.pi*2*nodes[...,0]/size[0])*np.pi*2/size[0]', '0.0', '0.0', + '0.0', '0.0', '0.0', + '0.0', '0.0', '0.0']), + + (['0.0', 'np.cos(np.pi*2*nodes[...,1]/size[1])', '0.0' ], + ['0.0', '0.0', '0.0', + '0.0', '-np.pi*2/size[1]*np.sin(np.pi*2*nodes[...,1]/size[1])', '0.0', + '0.0', '0.0', '0.0' ]), + + (['1.0', '0.0', '2.0*np.cos(np.pi*2*nodes[...,2]/size[2])'], + ['0.0', '0.0', '0.0', + '0.0', '0.0', '0.0', + '0.0', '0.0', '-2.0*np.pi*2/size[2]*np.sin(np.pi*2*nodes[...,2]/size[2])']), + + (['np.cos(np.pi*2*nodes[...,2]/size[2])', '3.0', 'np.sin(np.pi*2*nodes[...,2]/size[2])'], + ['0.0', '0.0', '-np.sin(np.pi*2*nodes[...,2]/size[2])*np.pi*2/size[2]', + '0.0', '0.0', '0.0', + '0.0', '0.0', ' np.cos(np.pi*2*nodes[...,2]/size[2])*np.pi*2/size[2]']), + + (['np.sin(np.pi*2*nodes[...,0]/size[0])', + 'np.sin(np.pi*2*nodes[...,1]/size[1])', + 'np.sin(np.pi*2*nodes[...,2]/size[2])'], + ['np.cos(np.pi*2*nodes[...,0]/size[0])*np.pi*2/size[0]', '0.0', '0.0', + '0.0', 'np.cos(np.pi*2*nodes[...,1]/size[1])*np.pi*2/size[1]', '0.0', + '0.0', '0.0', 'np.cos(np.pi*2*nodes[...,2]/size[2])*np.pi*2/size[2]']), + + (['np.sin(np.pi*2*nodes[...,0]/size[0])'], + ['np.cos(np.pi*2*nodes[...,0]/size[0])*np.pi*2/size[0]', '0.0', '0.0']), + + (['8.0'], + ['0.0', '0.0', '0.0' ]) + ] + + @pytest.mark.parametrize('field_def,grad_def',grad_test_data) + def test_grad(self,field_def,grad_def): + size = np.random.random(3)+1.0 + grid = np.random.randint(8,32,(3)) + + nodes = grid_filters.cell_coord0(grid,size) + my_locals = locals() # needed for list comprehension + + field = np.stack([np.broadcast_to(eval(f,globals(),my_locals),grid) for f in field_def],axis=-1) + field = field.reshape(tuple(grid) + ((3,) if len(field_def)==3 else (1,))) + grad = np.stack([np.broadcast_to(eval(c,globals(),my_locals),grid) for c in grad_def], axis=-1) + grad = grad.reshape(tuple(grid) + ((3,3) if len(grad_def)==9 else (3,))) + + assert np.allclose(grad,grid_filters.gradient(size,field)) + + + curl_test_data = [ + (['np.sin(np.pi*2*nodes[...,2]/size[2])', '0.0', '0.0', + '0.0', '0.0', '0.0', + '0.0', '0.0', '0.0'], + ['0.0' , '0.0', '0.0', + 'np.cos(np.pi*2*nodes[...,2]/size[2])*np.pi*2/size[2]', '0.0', '0.0', + '0.0', '0.0', '0.0']), + + (['np.cos(np.pi*2*nodes[...,1]/size[1])', '0.0', '0.0', + '0.0', '0.0', '0.0', + 'np.cos(np.pi*2*nodes[...,0]/size[0])', '0.0', '0.0'], + ['0.0', '0.0', '0.0', + '0.0', '0.0', '0.0', + 'np.sin(np.pi*2*nodes[...,1]/size[1])*np.pi*2/size[1]', '0.0', '0.0']), + + (['np.sin(np.pi*2*nodes[...,0]/size[0])','np.cos(np.pi*2*nodes[...,1]/size[1])','np.sin(np.pi*2*nodes[...,2]/size[2])', + 'np.sin(np.pi*2*nodes[...,0]/size[0])','np.cos(np.pi*2*nodes[...,1]/size[1])','np.sin(np.pi*2*nodes[...,2]/size[2])', + 'np.sin(np.pi*2*nodes[...,0]/size[0])','np.cos(np.pi*2*nodes[...,1]/size[1])','np.sin(np.pi*2*nodes[...,2]/size[2])'], + ['0.0', '0.0', '0.0', + '0.0', '0.0', '0.0', + '0.0', '0.0', '0.0']), + + (['5.0', '0.0', '0.0', + '0.0', '0.0', '0.0', + '0.0', '0.0', '2*np.cos(np.pi*2*nodes[...,1]/size[1])'], + ['0.0', '0.0', '-2*np.pi*2/size[1]*np.sin(np.pi*2*nodes[...,1]/size[1])', + '0.0', '0.0', '0.0', + '0.0', '0.0', '0.0']), + + ([ '4*np.sin(np.pi*2*nodes[...,2]/size[2])', + '8*np.sin(np.pi*2*nodes[...,0]/size[0])', + '16*np.sin(np.pi*2*nodes[...,1]/size[1])'], + ['16*np.pi*2/size[1]*np.cos(np.pi*2*nodes[...,1]/size[1])', + '4*np.pi*2/size[2]*np.cos(np.pi*2*nodes[...,2]/size[2])', + '8*np.pi*2/size[0]*np.cos(np.pi*2*nodes[...,0]/size[0])']), + + (['0.0', + 'np.cos(np.pi*2*nodes[...,0]/size[0])+5*np.cos(np.pi*2*nodes[...,2]/size[2])', + '0.0'], + ['5*np.sin(np.pi*2*nodes[...,2]/size[2])*np.pi*2/size[2]', + '0.0', + '-np.sin(np.pi*2*nodes[...,0]/size[0])*np.pi*2/size[0]']) + ] + + @pytest.mark.parametrize('field_def,curl_def',curl_test_data) + def test_curl(self,field_def,curl_def): + size = np.random.random(3)+1.0 + grid = np.random.randint(8,32,(3)) + + nodes = grid_filters.cell_coord0(grid,size) + my_locals = locals() # needed for list comprehension + + field = np.stack([np.broadcast_to(eval(f,globals(),my_locals),grid) for f in field_def],axis=-1) + field = field.reshape(tuple(grid) + ((3,3) if len(field_def)==9 else (3,))) + curl = np.stack([np.broadcast_to(eval(c,globals(),my_locals),grid) for c in curl_def], axis=-1) + curl = curl.reshape(tuple(grid) + ((3,3) if len(curl_def)==9 else (3,))) + + assert np.allclose(curl,grid_filters.curl(size,field)) + + + div_test_data =[ + (['np.sin(np.pi*2*nodes[...,0]/size[0])', '0.0', '0.0', + '0.0' , '0.0', '0.0', + '0.0' , '0.0', '0.0'], + ['np.cos(np.pi*2*nodes[...,0]/size[0])*np.pi*2/size[0]','0.0', '0.0']), + + (['0.0', '0.0', '0.0', + '0.0', 'np.cos(np.pi*2*nodes[...,1]/size[1])', '0.0', + '0.0', '0.0', '0.0'], + ['0.0', '-np.sin(np.pi*2*nodes[...,1]/size[1])*np.pi*2/size[1]', '0.0']), + + (['1.0', '0.0', '0.0', + '0.0', '0.0', '0.0', + '0.0', '0.0', '2*np.cos(np.pi*2*nodes[...,2]/size[2])' ], + ['0.0', '0.0', '-2.0*np.pi*2/size[2]*np.sin(np.pi*2*nodes[...,2]/size[2])'] + ), + + ([ '23.0', '0.0', 'np.sin(np.pi*2*nodes[...,2]/size[2])', + '0.0', '100.0', 'np.sin(np.pi*2*nodes[...,2]/size[2])', + '0.0', '0.0', 'np.sin(np.pi*2*nodes[...,2]/size[2])'], + ['np.cos(np.pi*2*nodes[...,2]/size[2])*np.pi*2/size[2]',\ + 'np.cos(np.pi*2*nodes[...,2]/size[2])*np.pi*2/size[2]', \ + 'np.cos(np.pi*2*nodes[...,2]/size[2])*np.pi*2/size[2]']), + + (['400.0', '0.0', '0.0', + 'np.sin(np.pi*2*nodes[...,0]/size[0])', 'np.sin(np.pi*2*nodes[...,1]/size[1])', 'np.sin(np.pi*2*nodes[...,2]/size[2])', + '0.0', '10.0', '6.0'], + ['0.0','np.sum(np.cos(np.pi*2*nodes/size)*np.pi*2/size,axis=-1)', '0.0' ]), + + (['np.sin(np.pi*2*nodes[...,0]/size[0])', '0.0', '0.0'], + ['np.cos(np.pi*2*nodes[...,0]/size[0])*np.pi*2/size[0]',]), + + (['0.0', 'np.cos(np.pi*2*nodes[...,1]/size[1])', '0.0' ], + ['-np.sin(np.pi*2*nodes[...,1]/size[1])*np.pi*2/size[1]']) + ] + + @pytest.mark.parametrize('field_def,div_def',div_test_data) + + def test_div(self,field_def,div_def): + size = np.random.random(3)+1.0 + grid = np.random.randint(8,32,(3)) + + nodes = grid_filters.cell_coord0(grid,size) + my_locals = locals() # needed for list comprehension + + field = np.stack([np.broadcast_to(eval(f,globals(),my_locals),grid) for f in field_def],axis=-1) + field = field.reshape(tuple(grid) + ((3,3) if len(field_def)==9 else (3,))) + div = np.stack([np.broadcast_to(eval(c,globals(),my_locals),grid) for c in div_def], axis=-1) + if len(div_def)==3: + div = div.reshape(tuple(grid) + ((3,))) + else: + div=div.reshape(tuple(grid)) + + assert np.allclose(div,grid_filters.divergence(size,field)) diff --git a/src/CPFEM.f90 b/src/CPFEM.f90 index 5394fc90c..2987d55f9 100644 --- a/src/CPFEM.f90 +++ b/src/CPFEM.f90 @@ -10,6 +10,7 @@ module CPFEM use FEsolving use math use rotations + use YAML_types use discretization_marc use material use config @@ -82,6 +83,7 @@ subroutine CPFEM_initAll(el,ip) call config_init call math_init call rotations_init + call YAML_types_init call HDF5_utilities_init call results_init call discretization_marc_init(ip, el) diff --git a/src/CPFEM2.f90 b/src/CPFEM2.f90 index 357bcce9f..9f61e4ebd 100644 --- a/src/CPFEM2.f90 +++ b/src/CPFEM2.f90 @@ -11,6 +11,7 @@ module CPFEM2 use FEsolving use math use rotations + use YAML_types use material use lattice use IO @@ -50,6 +51,7 @@ subroutine CPFEM_initAll call config_init call math_init call rotations_init + call YAML_types_init call lattice_init call HDF5_utilities_init call results_init diff --git a/src/IO.f90 b/src/IO.f90 index 6e467cb89..1bd2df833 100644 --- a/src/IO.f90 +++ b/src/IO.f90 @@ -29,9 +29,13 @@ module IO IO_getTag, & IO_stringPos, & IO_stringValue, & - IO_floatValue, & IO_intValue, & + IO_floatValue, & IO_lc, & + IO_rmComment, & + IO_stringAsInt, & + IO_stringAsFloat, & + IO_stringAsBool, & IO_error, & IO_warning @@ -250,7 +254,7 @@ integer function IO_intValue(string,chunkPos,myChunk) integer, dimension(:), intent(in) :: chunkPos !< positions of start and end of each tag/chunk in given string integer, intent(in) :: myChunk !< position number of desired chunk - IO_intValue = verifyIntValue(IO_stringValue(string,chunkPos,myChunk)) + IO_intValue = IO_stringAsInt(IO_stringValue(string,chunkPos,myChunk)) end function IO_intValue @@ -264,7 +268,7 @@ real(pReal) function IO_floatValue(string,chunkPos,myChunk) integer, dimension(:), intent(in) :: chunkPos !< positions of start and end of each tag/chunk in given string integer, intent(in) :: myChunk !< position number of desired chunk - IO_floatValue = verifyFloatValue(IO_stringValue(string,chunkPos,myChunk)) + IO_floatValue = IO_stringAsFloat(IO_stringValue(string,chunkPos,myChunk)) end function IO_floatValue @@ -294,6 +298,88 @@ pure function IO_lc(string) end function IO_lc +!-------------------------------------------------------------------------------------------------- +! @brief Remove comments (characters beyond '#') and trailing space +! ToDo: Discuss name (the trim aspect is not clear) +!-------------------------------------------------------------------------------------------------- +function IO_rmComment(line) + + character(len=*), intent(in) :: line + character(len=:), allocatable :: IO_rmComment + integer :: split + + split = index(line,IO_COMMENT) + + if (split == 0) then + IO_rmComment = trim(line) + else + IO_rmComment = trim(line(:split-1)) + endif + +end function IO_rmComment + + +!-------------------------------------------------------------------------------------------------- +!> @brief return verified integer value in given string +!-------------------------------------------------------------------------------------------------- +integer function IO_stringAsInt(string) + + character(len=*), intent(in) :: string !< string for conversion to int value + + integer :: readStatus + character(len=*), parameter :: VALIDCHARS = '0123456789+- ' + + valid: if (verify(string,VALIDCHARS) == 0) then + read(string,*,iostat=readStatus) IO_stringAsInt + if (readStatus /= 0) call IO_error(111,ext_msg=string) + else valid + IO_stringAsInt = 0 + call IO_error(111,ext_msg=string) + endif valid + +end function IO_stringAsInt + + +!-------------------------------------------------------------------------------------------------- +!> @brief return verified float value in given string +!-------------------------------------------------------------------------------------------------- +real(pReal) function IO_stringAsFloat(string) + + character(len=*), intent(in) :: string !< string for conversion to float value + + integer :: readStatus + character(len=*), parameter :: VALIDCHARS = '0123456789eE.+- ' + + valid: if (verify(string,VALIDCHARS) == 0) then + read(string,*,iostat=readStatus) IO_stringAsFloat + if (readStatus /= 0) call IO_error(112,ext_msg=string) + else valid + IO_stringAsFloat = 0.0_pReal + call IO_error(112,ext_msg=string) + endif valid + +end function IO_stringAsFloat + + +!-------------------------------------------------------------------------------------------------- +!> @brief return verified logical value in given string +!-------------------------------------------------------------------------------------------------- +logical function IO_stringAsBool(string) + + character(len=*), intent(in) :: string !< string for conversion to int value + + if (trim(adjustl(string)) == 'True') then + IO_stringAsBool = .true. + elseif (trim(adjustl(string)) == 'False') then + IO_stringAsBool = .false. + else + IO_stringAsBool = .false. + call IO_error(113,ext_msg=string) + endif + +end function IO_stringAsBool + + !-------------------------------------------------------------------------------------------------- !> @brief write error statements to standard out and terminate the Marc/spectral run with exit #9xxx !-------------------------------------------------------------------------------------------------- @@ -335,7 +421,8 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg) msg = 'invalid character for int:' case (112) msg = 'invalid character for float:' - + case (113) + msg = 'invalid character for logical:' !-------------------------------------------------------------------------------------------------- ! lattice error messages case (130) @@ -606,51 +693,6 @@ subroutine IO_warning(warning_ID,el,ip,g,ext_msg) end subroutine IO_warning -!-------------------------------------------------------------------------------------------------- -! internal helper functions - -!-------------------------------------------------------------------------------------------------- -!> @brief returns verified integer value in given string -!-------------------------------------------------------------------------------------------------- -integer function verifyIntValue(string) - - character(len=*), intent(in) :: string !< string for conversion to int value - - integer :: readStatus - character(len=*), parameter :: VALIDCHARS = '0123456789+- ' - - valid: if (verify(string,VALIDCHARS) == 0) then - read(string,*,iostat=readStatus) verifyIntValue - if (readStatus /= 0) call IO_error(111,ext_msg=string) - else valid - verifyIntValue = 0 - call IO_error(111,ext_msg=string) - endif valid - -end function verifyIntValue - - -!-------------------------------------------------------------------------------------------------- -!> @brief returns verified float value in given string -!-------------------------------------------------------------------------------------------------- -real(pReal) function verifyFloatValue(string) - - character(len=*), intent(in) :: string !< string for conversion to float value - - integer :: readStatus - character(len=*), parameter :: VALIDCHARS = '0123456789eE.+- ' - - valid: if (verify(string,VALIDCHARS) == 0) then - read(string,*,iostat=readStatus) verifyFloatValue - if (readStatus /= 0) call IO_error(112,ext_msg=string) - else valid - verifyFloatValue = 0.0_pReal - call IO_error(112,ext_msg=string) - endif valid - -end function verifyFloatValue - - !-------------------------------------------------------------------------------------------------- !> @brief check correctness of some IO functions !-------------------------------------------------------------------------------------------------- @@ -659,14 +701,19 @@ subroutine unitTest integer, dimension(:), allocatable :: chunkPos character(len=:), allocatable :: str - if(dNeq(1.0_pReal, verifyFloatValue('1.0'))) call IO_error(0,ext_msg='verifyFloatValue') - if(dNeq(1.0_pReal, verifyFloatValue('1e0'))) call IO_error(0,ext_msg='verifyFloatValue') - if(dNeq(0.1_pReal, verifyFloatValue('1e-1'))) call IO_error(0,ext_msg='verifyFloatValue') + if(dNeq(1.0_pReal, IO_stringAsFloat('1.0'))) call IO_error(0,ext_msg='IO_stringAsFloat') + if(dNeq(1.0_pReal, IO_stringAsFloat('1e0'))) call IO_error(0,ext_msg='IO_stringAsFloat') + if(dNeq(0.1_pReal, IO_stringAsFloat('1e-1'))) call IO_error(0,ext_msg='IO_stringAsFloat') - if(3112019 /= verifyIntValue( '3112019')) call IO_error(0,ext_msg='verifyIntValue') - if(3112019 /= verifyIntValue(' 3112019')) call IO_error(0,ext_msg='verifyIntValue') - if(-3112019 /= verifyIntValue('-3112019')) call IO_error(0,ext_msg='verifyIntValue') - if(3112019 /= verifyIntValue('+3112019 ')) call IO_error(0,ext_msg='verifyIntValue') + if(3112019 /= IO_stringAsInt( '3112019')) call IO_error(0,ext_msg='IO_stringAsInt') + if(3112019 /= IO_stringAsInt(' 3112019')) call IO_error(0,ext_msg='IO_stringAsInt') + if(-3112019 /= IO_stringAsInt('-3112019')) call IO_error(0,ext_msg='IO_stringAsInt') + if(3112019 /= IO_stringAsInt('+3112019 ')) call IO_error(0,ext_msg='IO_stringAsInt') + + if(.not. IO_stringAsBool(' True')) call IO_error(0,ext_msg='IO_stringAsBool') + if(.not. IO_stringAsBool(' True ')) call IO_error(0,ext_msg='IO_stringAsBool') + if( IO_stringAsBool(' False')) call IO_error(0,ext_msg='IO_stringAsBool') + if( IO_stringAsBool('False')) call IO_error(0,ext_msg='IO_stringAsBool') if(any([1,1,1] /= IO_stringPos('a'))) call IO_error(0,ext_msg='IO_stringPos') if(any([2,2,3,5,5] /= IO_stringPos(' aa b'))) call IO_error(0,ext_msg='IO_stringPos') @@ -683,6 +730,21 @@ subroutine unitTest if(.not. IO_isBlank(' #isBlank')) call IO_error(0,ext_msg='IO_isBlank/2') if( IO_isBlank(' i#s')) call IO_error(0,ext_msg='IO_isBlank/3') + str = IO_rmComment('#') + if (str /= '' .or. len(str) /= 0) call IO_error(0,ext_msg='IO_rmComment/1') + str = IO_rmComment(' #') + if (str /= '' .or. len(str) /= 0) call IO_error(0,ext_msg='IO_rmComment/2') + str = IO_rmComment(' # ') + if (str /= '' .or. len(str) /= 0) call IO_error(0,ext_msg='IO_rmComment/3') + str = IO_rmComment(' # a') + if (str /= '' .or. len(str) /= 0) call IO_error(0,ext_msg='IO_rmComment/4') + str = IO_rmComment(' # a') + if (str /= '' .or. len(str) /= 0) call IO_error(0,ext_msg='IO_rmComment/5') + str = IO_rmComment(' a#') + if (str /= ' a' .or. len(str) /= 2) call IO_error(0,ext_msg='IO_rmComment/6') + str = IO_rmComment(' ab #') + if (str /= ' ab'.or. len(str) /= 3) call IO_error(0,ext_msg='IO_rmComment/7') + end subroutine unitTest end module IO diff --git a/src/Lambert.f90 b/src/Lambert.f90 deleted file mode 100644 index 932fe221b..000000000 --- a/src/Lambert.f90 +++ /dev/null @@ -1,201 +0,0 @@ -! ################################################################### -! Copyright (c) 2013-2015, Marc De Graef/Carnegie Mellon University -! Modified 2017-2019, Martin Diehl/Max-Planck-Institut für Eisenforschung GmbH -! All rights reserved. -! -! Redistribution and use in source and binary forms, with or without modification, are -! permitted provided that the following conditions are met: -! -! - Redistributions of source code must retain the above copyright notice, this list -! of conditions and the following disclaimer. -! - Redistributions in binary form must reproduce the above copyright notice, this -! list of conditions and the following disclaimer in the documentation and/or -! other materials provided with the distribution. -! - Neither the names of Marc De Graef, Carnegie Mellon University nor the names -! of its contributors may be used to endorse or promote products derived from -! this software without specific prior written permission. -! -! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" -! AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE -! IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE -! ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE -! LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL -! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR -! SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER -! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, -! OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE -! USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -! ################################################################### - -!-------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH -!> @brief Mapping homochoric <-> cubochoric -! -!> @details -!> D. Rosca, A. Morawiec, and M. De Graef. “A new method of constructing a grid -!> in the space of 3D rotations and its applications to texture analysis”. -!> Modeling and Simulations in Materials Science and Engineering 22, 075013 (2014). -!-------------------------------------------------------------------------- -module Lambert - use prec - use math - - implicit none - private - - real(pReal), parameter :: & - SPI = sqrt(PI), & - PREF = sqrt(6.0_pReal/PI), & - A = PI**(5.0_pReal/6.0_pReal)/6.0_pReal**(1.0_pReal/6.0_pReal), & - AP = PI**(2.0_pReal/3.0_pReal), & - SC = A/AP, & - BETA = A/2.0_pReal, & - R1 = (3.0_pReal*PI/4.0_pReal)**(1.0_pReal/3.0_pReal), & - R2 = sqrt(2.0_pReal), & - PI12 = PI/12.0_pReal, & - PREK = R1 * 2.0_pReal**(1.0_pReal/4.0_pReal)/BETA - - public :: & - Lambert_CubeToBall, & - Lambert_BallToCube - -contains - - -!-------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH -!> @brief map from 3D cubic grid to 3D ball -!-------------------------------------------------------------------------- -pure function Lambert_CubeToBall(cube) result(ball) - - real(pReal), intent(in), dimension(3) :: cube - real(pReal), dimension(3) :: ball, LamXYZ, XYZ - real(pReal), dimension(2) :: T - real(pReal) :: c, s, q - real(pReal), parameter :: eps = 1.0e-8_pReal - integer, dimension(3,2) :: p - integer, dimension(2) :: order - - if (maxval(abs(cube)) > AP/2.0+eps) then - ball = IEEE_value(cube,IEEE_positive_inf) - return - end if - - ! transform to the sphere grid via the curved square, and intercept the zero point - center: if (all(dEq0(cube))) then - ball = 0.0_pReal - else center - ! get pyramide and scale by grid parameter ratio - p = GetPyramidOrder(cube) - XYZ = cube(p(:,1)) * sc - - ! intercept all the points along the z-axis - special: if (all(dEq0(XYZ(1:2)))) then - LamXYZ = [ 0.0_pReal, 0.0_pReal, pref * XYZ(3) ] - else special - order = merge( [2,1], [1,2], abs(XYZ(2)) <= abs(XYZ(1))) ! order of absolute values of XYZ - q = PI12 * XYZ(order(1))/XYZ(order(2)) ! smaller by larger - c = cos(q) - s = sin(q) - q = prek * XYZ(order(2))/ sqrt(R2-c) - T = [ (R2*c - 1.0), R2 * s] * q - - ! transform to sphere grid (inverse Lambert) - ! [note that there is no need to worry about dividing by zero, since XYZ(3) can not become zero] - c = sum(T**2) - s = Pi * c/(24.0*XYZ(3)**2) - c = sPi * c / sqrt(24.0_pReal) / XYZ(3) - q = sqrt( 1.0 - s ) - LamXYZ = [ T(order(2)) * q, T(order(1)) * q, pref * XYZ(3) - c ] - endif special - - ! reverse the coordinates back to order according to the original pyramid number - ball = LamXYZ(p(:,2)) - - endif center - -end function Lambert_CubeToBall - - -!-------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH -!> @brief map from 3D ball to 3D cubic grid -!-------------------------------------------------------------------------- -pure function Lambert_BallToCube(xyz) result(cube) - - real(pReal), intent(in), dimension(3) :: xyz - real(pReal), dimension(3) :: cube, xyz1, xyz3 - real(pReal), dimension(2) :: Tinv, xyz2 - real(pReal) :: rs, qxy, q2, sq2, q, tt - integer, dimension(3,2) :: p - - rs = norm2(xyz) - if (rs > R1) then - cube = IEEE_value(cube,IEEE_positive_inf) - return - endif - - center: if (all(dEq0(xyz))) then - cube = 0.0_pReal - else center - p = GetPyramidOrder(xyz) - xyz3 = xyz(p(:,1)) - - ! inverse M_3 - xyz2 = xyz3(1:2) * sqrt( 2.0*rs/(rs+abs(xyz3(3))) ) - - ! inverse M_2 - qxy = sum(xyz2**2) - - special: if (dEq0(qxy)) then - Tinv = 0.0_pReal - else special - q2 = qxy + maxval(abs(xyz2))**2 - sq2 = sqrt(q2) - q = (beta/R2/R1) * sqrt(q2*qxy/(q2-maxval(abs(xyz2))*sq2)) - tt = (minval(abs(xyz2))**2+maxval(abs(xyz2))*sq2)/R2/qxy - Tinv = q * sign(1.0_pReal,xyz2) * merge([ 1.0_pReal, acos(math_clip(tt,-1.0_pReal,1.0_pReal))/PI12], & - [ acos(math_clip(tt,-1.0_pReal,1.0_pReal))/PI12, 1.0_pReal], & - abs(xyz2(2)) <= abs(xyz2(1))) - endif special - - ! inverse M_1 - xyz1 = [ Tinv(1), Tinv(2), sign(1.0_pReal,xyz3(3)) * rs / pref ] /sc - - ! reverse the coordinates back to order according to the original pyramid number - cube = xyz1(p(:,2)) - - endif center - -end function Lambert_BallToCube - - -!-------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH -!> @brief determine to which pyramid a point in a cubic grid belongs -!-------------------------------------------------------------------------- -pure function GetPyramidOrder(xyz) - - real(pReal),intent(in),dimension(3) :: xyz - integer, dimension(3,2) :: GetPyramidOrder - - if (((abs(xyz(1)) <= xyz(3)).and.(abs(xyz(2)) <= xyz(3))) .or. & - ((abs(xyz(1)) <= -xyz(3)).and.(abs(xyz(2)) <= -xyz(3)))) then - GetPyramidOrder = reshape([[1,2,3],[1,2,3]],[3,2]) - else if (((abs(xyz(3)) <= xyz(1)).and.(abs(xyz(2)) <= xyz(1))) .or. & - ((abs(xyz(3)) <= -xyz(1)).and.(abs(xyz(2)) <= -xyz(1)))) then - GetPyramidOrder = reshape([[2,3,1],[3,1,2]],[3,2]) - else if (((abs(xyz(1)) <= xyz(2)).and.(abs(xyz(3)) <= xyz(2))) .or. & - ((abs(xyz(1)) <= -xyz(2)).and.(abs(xyz(3)) <= -xyz(2)))) then - GetPyramidOrder = reshape([[3,1,2],[2,3,1]],[3,2]) - else - GetPyramidOrder = -1 ! should be impossible, but might simplify debugging - end if - -end function GetPyramidOrder - -end module Lambert diff --git a/src/YAML_types.f90 b/src/YAML_types.f90 new file mode 100644 index 000000000..07541bc1a --- /dev/null +++ b/src/YAML_types.f90 @@ -0,0 +1,1044 @@ +!-------------------------------------------------------------------------------------------------- +!> @author Sharan Roongta, Max-Planck-Institut für Eisenforschung GmbH +!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH +!> @brief Data types to create a scalar, a list, and a dictionary/hash +!> @details module describes the various functions to store and get the yaml data. +!! A node is the base class for scalar, list and dictionary, list items and dictionary entries point +!! to a node. +!-------------------------------------------------------------------------------------------------- + +module YAML_types + + use IO + use prec + + implicit none + + private + + public :: & + tNode, & + tScalar, & + tDict, & + tList, & + YAML_types_init + + type, abstract :: tNode + integer :: length = 0 + contains + procedure(asFormattedString), deferred :: asFormattedString + procedure :: & + asScalar => tNode_asScalar + procedure :: & + asList => tNode_asList + procedure :: & + asDict => tNode_asDict + procedure :: & + tNode_get_byIndex => tNode_get_byIndex + procedure :: & + tNode_get_byIndex_asFloat => tNode_get_byIndex_asFloat + procedure :: & + tNode_get_byIndex_asFloats => tNode_get_byIndex_asFloats + procedure :: & + tNode_get_byIndex_asInt => tNode_get_byIndex_asInt + procedure :: & + tNode_get_byIndex_asInts => tNode_get_byIndex_asInts + procedure :: & + tNode_get_byIndex_asBool => tNode_get_byIndex_asBool + procedure :: & + tNode_get_byIndex_asBools => tNode_get_byIndex_asBools + procedure :: & + tNode_get_byIndex_asString => tNode_get_byIndex_asString + procedure :: & + tNode_get_byIndex_asStrings => tNode_get_byIndex_asStrings + procedure :: & + tNode_get_byKey => tNode_get_byKey + procedure :: & + tNode_get_byKey_asFloat => tNode_get_byKey_asFloat + procedure :: & + tNode_get_byKey_asFloats => tNode_get_byKey_asFloats + procedure :: & + tNode_get_byKey_asInt => tNode_get_byKey_asInt + procedure :: & + tNode_get_byKey_asInts => tNode_get_byKey_asInts + procedure :: & + tNode_get_byKey_asBool => tNode_get_byKey_asBool + procedure :: & + tNode_get_byKey_asBools => tNode_get_byKey_asBools + procedure :: & + tNode_get_byKey_asString => tNode_get_byKey_asString + procedure :: & + tNode_get_byKey_asStrings => tNode_get_byKey_asStrings + procedure :: & + getIndex => tNode_get_byKey_asIndex + + generic :: & + get => tNode_get_byIndex, & + tNode_get_byKey + generic :: & + get_asFloat => tNode_get_byIndex_asFloat, & + tNode_get_byKey_asFloat + generic :: & + get_asFloats => tNode_get_byIndex_asFloats, & + tNode_get_byKey_asFloats + generic :: & + get_asInt => tNode_get_byIndex_asInt, & + tNode_get_byKey_asInt + generic :: & + get_asInts => tNode_get_byIndex_asInts, & + tNode_get_byKey_asInts + generic :: & + get_asBool => tNode_get_byIndex_asBool, & + tNode_get_byKey_asBool + generic :: & + get_asBools => tNode_get_byIndex_asBools, & + tNode_get_byKey_asBools + generic :: & + get_asString => tNode_get_byIndex_asString, & + tNode_get_byKey_asString + generic :: & + get_asStrings => tNode_get_byIndex_asStrings, & + tNode_get_byKey_asStrings + end type tNode + + + type, extends(tNode) :: tScalar + + character(len=:), allocatable, private :: value + + contains + procedure :: asFormattedString => tScalar_asFormattedString + procedure :: & + asFloat => tScalar_asFloat + procedure :: & + asInt => tScalar_asInt + procedure :: & + asBool => tScalar_asBool + procedure :: & + asString => tScalar_asString + end type tScalar + + type, extends(tNode) :: tList + + class(tItem), pointer :: first => null() + + contains + procedure :: asFormattedString => tList_asFormattedString + procedure :: append => tList_append + procedure :: & + asFloats => tList_asFloats + procedure :: & + asInts => tList_asInts + procedure :: & + asBools => tList_asBools + procedure :: & + asStrings => tList_asStrings + final :: tList_finalize + end type tList + + type, extends(tList) :: tDict + contains + procedure :: asFormattedString => tDict_asFormattedString + procedure :: set => tDict_set + end type tDict + + + type :: tItem + character(len=:), allocatable :: key + class(tNode), pointer :: node => null() + class(tItem), pointer :: next => null() + + contains + final :: tItem_finalize + end type tItem + + abstract interface + + recursive function asFormattedString(self,indent) + import tNode + character(len=:), allocatable :: asFormattedString + class(tNode), intent(in), target :: self + integer, intent(in), optional :: indent + end function asFormattedString + + end interface + + interface tScalar + module procedure tScalar_init__ + end interface tScalar + + interface assignment (=) + module procedure tScalar_assign__ + end interface assignment (=) + +contains + +!-------------------------------------------------------------------------------------------------- +!> @brief do sanity checks +!-------------------------------------------------------------------------------------------------- +subroutine YAML_types_init + + write(6,'(/,a)') ' <<<+- YAML_types init -+>>>' + + call unitTest + +end subroutine YAML_types_init + + +!-------------------------------------------------------------------------------------------------- +!> @brief check correctness of some type bound procedures +!-------------------------------------------------------------------------------------------------- +subroutine unitTest + + class(tNode), pointer :: s1,s2 + allocate(tScalar::s1) + allocate(tScalar::s2) + select type(s1) + class is(tScalar) + s1 = '1' + if(s1%asInt() /= 1) call IO_error(0,ext_msg='tScalar_asInt') + if(dNeq(s1%asFloat(),1.0_pReal)) call IO_error(0,ext_msg='tScalar_asFloat') + s1 = 'True' + if(.not. s1%asBool()) call IO_error(0,ext_msg='tScalar_asBool') + if(s1%asString() /= 'True') call IO_error(0,ext_msg='tScalar_asString') + end select + + block + class(tNode), pointer :: l1, l2, n + select type(s1) + class is(tScalar) + s1 = '2' + endselect + + select type(s2) + class is(tScalar) + s2 = '3' + endselect + + allocate(tList::l1) + select type(l1) + class is(tList) + call l1%append(s1) + call l1%append(s2) + n => l1 + if(any(l1%asInts() /= [2,3])) call IO_error(0,ext_msg='tList_asInts') + if(any(dNeq(l1%asFloats(),[2.0_pReal,3.0_pReal]))) call IO_error(0,ext_msg='tList_asFloats') + if(n%get_asInt(1) /= 2) call IO_error(0,ext_msg='byIndex_asInt') + if(dNeq(n%get_asFloat(2),3.0_pReal)) call IO_error(0,ext_msg='byIndex_asFloat') + endselect + + allocate(tList::l2) + select type(l2) + class is(tList) + call l2%append(l1) + if(any(l2%get_asInts(1) /= [2,3])) call IO_error(0,ext_msg='byIndex_asInts') + if(any(dNeq(l2%get_asFloats(1),[2.0_pReal,3.0_pReal]))) call IO_error(0,ext_msg='byIndex_asFloats') + n => l2 + end select + deallocate(n) + end block + + block + type(tList), target :: l1 + type(tScalar),pointer :: s3,s4 + class(tNode), pointer :: n + + allocate(tScalar::s1) + allocate(tScalar::s2) + s3 => s1%asScalar() + s4 => s2%asScalar() + s3 = 'True' + s4 = 'False' + + call l1%append(s1) + call l1%append(s2) + n => l1 + + if(any(l1%asBools() .neqv. [.true., .false.])) call IO_error(0,ext_msg='tList_asBools') + if(any(l1%asStrings() /= ['True ','False'])) call IO_error(0,ext_msg='tList_asStrings') + if(n%get_asBool(2)) call IO_error(0,ext_msg='byIndex_asBool') + if(n%get_asString(1) /= 'True') call IO_error(0,ext_msg='byIndex_asString') + end block + +end subroutine unitTest + + +!--------------------------------------------------------------------------------------------------- +!> @brief init from string +!--------------------------------------------------------------------------------------------------- +type(tScalar) pure function tScalar_init__(value) + + character(len=*), intent(in) :: value + + tScalar_init__%value =value + +end function tScalar_init__ + + +!--------------------------------------------------------------------------------------------------- +!> @brief set value from string +!--------------------------------------------------------------------------------------------------- +elemental pure subroutine tScalar_assign__(self,value) + + type(tScalar), intent(out) :: self + character(len=*), intent(in) :: value + + self%value = value + +end subroutine tScalar_assign__ + + +!-------------------------------------------------------------------------------------------------- +!> @brief Type guard, guarantee scalar +!-------------------------------------------------------------------------------------------------- +function tNode_asScalar(self) result(scalar) + + class(tNode), intent(in), target :: self + class(tScalar), pointer :: scalar + + select type(self) + class is(tScalar) + scalar => self + class default + call IO_error(0) + end select + +end function tNode_asScalar + + +!-------------------------------------------------------------------------------------------------- +!> @brief Type guard, guarantee list +!-------------------------------------------------------------------------------------------------- +function tNode_asList(self) result(list) + + class(tNode), intent(in), target :: self + class(tList), pointer :: list + + select type(self) + class is(tList) + list => self + class default + call IO_error(0) + end select + +end function tNode_asList + + +!-------------------------------------------------------------------------------------------------- +!> @brief Type guard, guarantee dict +!-------------------------------------------------------------------------------------------------- +function tNode_asDict(self) result(dict) + + class(tNode), intent(in), target :: self + class(tDict), pointer :: dict + + select type(self) + class is(tDict) + dict => self + class default + call IO_error(0) + end select + +end function tNode_asDict + + +!-------------------------------------------------------------------------------------------------- +!> @brief Access by index +!-------------------------------------------------------------------------------------------------- +function tNode_get_byIndex(self,i) result(node) + + class(tNode), intent(in), target :: self + integer, intent(in) :: i + class(tNode), pointer :: node + + class(tList), pointer :: self_ + class(tItem), pointer :: item + integer :: j + + self_ => self%asList() + if(i < 1 .or. i > self_%length) call IO_error(0) + + j = 1 + item => self_%first + do while(j item%next + j = j + 1 + enddo + node => item%node + +end function tNode_get_byIndex + + +!-------------------------------------------------------------------------------------------------- +!> @brief Access by index and convert to float +!-------------------------------------------------------------------------------------------------- +function tNode_get_byIndex_asFloat(self,i) result(nodeAsFloat) + + class(tNode), intent(in), target :: self + integer, intent(in) :: i + real(pReal) :: nodeAsFloat + + class(tNode), pointer :: node + type(tScalar), pointer :: scalar + + node => self%get(i) + scalar => node%asScalar() + nodeAsFloat = scalar%asFloat() + +end function tNode_get_byIndex_asFloat + + +!-------------------------------------------------------------------------------------------------- +!> @brief Access by index and convert to int +!-------------------------------------------------------------------------------------------------- +function tNode_get_byIndex_asInt(self,i) result(nodeAsInt) + + class(tNode), intent(in), target :: self + integer, intent(in) :: i + integer :: nodeAsInt + + class(tNode), pointer :: node + type(tScalar), pointer :: scalar + + node => self%get(i) + scalar => node%asScalar() + nodeAsInt = scalar%asInt() + +end function tNode_get_byIndex_asInt + + +!-------------------------------------------------------------------------------------------------- +!> @brief Access by index and convert to bool +!-------------------------------------------------------------------------------------------------- +function tNode_get_byIndex_asBool(self,i) result(nodeAsBool) + + class(tNode), intent(in), target :: self + integer, intent(in) :: i + logical :: nodeAsBool + + class(tNode), pointer :: node + type(tScalar), pointer :: scalar + + node => self%get(i) + scalar => node%asScalar() + nodeAsBool = scalar%asBool() + +end function tNode_get_byIndex_asBool + + +!-------------------------------------------------------------------------------------------------- +!> @brief Access by index and convert to string +!-------------------------------------------------------------------------------------------------- +function tNode_get_byIndex_asString(self,i) result(nodeAsString) + + class(tNode), intent(in), target :: self + integer, intent(in) :: i + character(len=:), allocatable :: nodeAsString + + class(tNode), pointer :: node + type(tScalar), pointer :: scalar + + node => self%get(i) + scalar => node%asScalar() + nodeAsString = scalar%asString() + +end function tNode_get_byIndex_asString + + +!-------------------------------------------------------------------------------------------------- +!> @brief Access by index and convert to float array +!-------------------------------------------------------------------------------------------------- +function tNode_get_byIndex_asFloats(self,i) result(nodeAsFloats) + + class(tNode), intent(in), target :: self + integer, intent(in) :: i + real(pReal), dimension(:), allocatable :: nodeAsFloats + + class(tNode), pointer :: node + class(tList), pointer :: list + + node => self%get(i) + list => node%asList() + nodeAsFloats = list%asFloats() + +end function tNode_get_byIndex_asFloats + + +!-------------------------------------------------------------------------------------------------- +!> @brief Access by index and convert to int array +!-------------------------------------------------------------------------------------------------- +function tNode_get_byIndex_asInts(self,i) result(nodeAsInts) + + class(tNode), intent(in), target :: self + integer, intent(in) :: i + integer, dimension(:), allocatable :: nodeAsInts + + class(tNode), pointer :: node + class(tList), pointer :: list + + node => self%get(i) + list => node%asList() + nodeAsInts = list%asInts() + +end function tNode_get_byIndex_asInts + + +!-------------------------------------------------------------------------------------------------- +!> @brief Access by index and convert to bool array +!-------------------------------------------------------------------------------------------------- +function tNode_get_byIndex_asBools(self,i) result(nodeAsBools) + + class(tNode), intent(in), target :: self + integer, intent(in) :: i + logical, dimension(:), allocatable :: nodeAsBools + + class(tNode), pointer :: node + class(tList), pointer :: list + + node => self%get(i) + list => node%asList() + nodeAsBools = list%asBools() + +end function tNode_get_byIndex_asBools + + +!-------------------------------------------------------------------------------------------------- +!> @brief Access by index and convert to string array +!-------------------------------------------------------------------------------------------------- +function tNode_get_byIndex_asStrings(self,i) result(nodeAsStrings) + + class(tNode), intent(in), target :: self + integer, intent(in) :: i + character(len=:), allocatable, dimension(:) :: nodeAsStrings + + class(tNode), pointer :: node + type(tList), pointer :: list + + node => self%get(i) + list => node%asList() + nodeAsStrings = list%asStrings() + +end function tNode_get_byIndex_asStrings + + +!-------------------------------------------------------------------------------------------------- +!> @brief Access by index +!-------------------------------------------------------------------------------------------------- +function tNode_get_byKey(self,k) result(node) + + class(tNode), intent(in), target :: self + character(len=*), intent(in) :: k + class(tNode), pointer :: node + + type(tDict), pointer :: self_ + type(tItem), pointer :: item + integer :: j + + self_ => self%asDict() + + j = 1 + item => self_%first + do while(j <= self_%length) + if (item%key == k) exit + item => item%next + j = j + 1 + enddo + if (.not. item%key == k) call IO_error(0) + node => item%node + +end function tNode_get_byKey + + +!-------------------------------------------------------------------------------------------------- +!> @brief Access by key and convert to float +!-------------------------------------------------------------------------------------------------- +function tNode_get_byKey_asFloat(self,k) result(nodeAsFloat) + + class(tNode), intent(in), target :: self + character(len=*), intent(in) :: k + real(pReal) :: nodeAsFloat + + class(tNode), pointer :: node + type(tScalar), pointer :: scalar + + node => self%get(k) + scalar => node%asScalar() + nodeAsFloat = scalar%asFloat() + +end function tNode_get_byKey_asFloat + + +!-------------------------------------------------------------------------------------------------- +!> @brief Access by key and convert to int +!-------------------------------------------------------------------------------------------------- +function tNode_get_byKey_asInt(self,k) result(nodeAsInt) + + class(tNode), intent(in), target :: self + character(len=*), intent(in) :: k + integer :: nodeAsInt + + class(tNode), pointer :: node + type(tScalar), pointer :: scalar + + node => self%get(k) + scalar => node%asScalar() + nodeAsInt = scalar%asInt() + +end function tNode_get_byKey_asInt + + +!-------------------------------------------------------------------------------------------------- +!> @brief Access by key and convert to bool +!-------------------------------------------------------------------------------------------------- +function tNode_get_byKey_asBool(self,k) result(nodeAsBool) + + class(tNode), intent(in), target :: self + character(len=*), intent(in) :: k + logical :: nodeAsBool + + class(tNode), pointer :: node + type(tScalar), pointer :: scalar + + node => self%get(k) + scalar => node%asScalar() + nodeAsBool = scalar%asBool() + +end function tNode_get_byKey_asBool + + +!-------------------------------------------------------------------------------------------------- +!> @brief Access by key and convert to string +!-------------------------------------------------------------------------------------------------- +function tNode_get_byKey_asString(self,k) result(nodeAsString) + + class(tNode), intent(in), target :: self + character(len=*), intent(in) :: k + character(len=:), allocatable :: nodeAsString + + class(tNode), pointer :: node + type(tScalar), pointer :: scalar + + node => self%get(k) + scalar => node%asScalar() + nodeAsString = scalar%asString() + +end function tNode_get_byKey_asString + + +!-------------------------------------------------------------------------------------------------- +!> @brief Access by key and convert to float array +!-------------------------------------------------------------------------------------------------- +function tNode_get_byKey_asFloats(self,k) result(nodeAsFloats) + + class(tNode), intent(in), target :: self + character(len=*), intent(in) :: k + real(pReal), dimension(:), allocatable :: nodeAsFloats + + class(tNode), pointer :: node + type(tList), pointer :: list + + node => self%get(k) + list => node%asList() + nodeAsFloats = list%asFloats() + +end function tNode_get_byKey_asFloats + + +!-------------------------------------------------------------------------------------------------- +!> @brief Access by key and convert to int array +!-------------------------------------------------------------------------------------------------- +function tNode_get_byKey_asInts(self,k) result(nodeAsInts) + + class(tNode), intent(in), target :: self + character(len=*), intent(in) :: k + integer, dimension(:), allocatable :: nodeAsInts + + class(tNode), pointer :: node + type(tList), pointer :: list + + node => self%get(k) + list => node%asList() + nodeAsInts = list%asInts() + +end function tNode_get_byKey_asInts + + +!-------------------------------------------------------------------------------------------------- +!> @brief Access by key and convert to bool array +!-------------------------------------------------------------------------------------------------- +function tNode_get_byKey_asBools(self,k) result(nodeAsBools) + + class(tNode), intent(in), target :: self + character(len=*), intent(in) :: k + logical, dimension(:), allocatable :: nodeAsBools + + class(tNode), pointer :: node + type(tList), pointer :: list + + node => self%get(k) + list => node%asList() + nodeAsBools = list%asBools() + +end function tNode_get_byKey_asBools + + +!-------------------------------------------------------------------------------------------------- +!> @brief Access by key and convert to string array +!-------------------------------------------------------------------------------------------------- +function tNode_get_byKey_asStrings(self,k) result(nodeAsStrings) + + class(tNode), intent(in), target :: self + character(len=*), intent(in) :: k + character(len=:), allocatable, dimension(:) :: nodeAsStrings + + class(tNode), pointer :: node + type(tList), pointer :: list + + node => self%get(k) + list => node%asList() + nodeAsStrings = list%asStrings() + +end function tNode_get_byKey_asStrings + + +!-------------------------------------------------------------------------------------------------- +!> @brief Returns the index of a key in a dictionary +!-------------------------------------------------------------------------------------------------- +function tNode_get_byKey_asIndex(self,key) result(keyIndex) + + class(tNode), intent(in), target :: self + character(len=*), intent(in) :: key + + integer :: keyIndex + integer :: i + type(tDict), pointer :: dict + type(tItem), pointer :: item + + dict => self%asDict() + item => dict%first + keyIndex = -1 + do i = 1, dict%length + if(key == item%key) then + keyIndex = i + exit + else + item => item%next + endif + enddo + +end function tNode_get_byKey_asIndex + + +!-------------------------------------------------------------------------------------------------- +!> @brief Scalar as string (YAML block style) +!-------------------------------------------------------------------------------------------------- +recursive function tScalar_asFormattedString(self,indent) + + character(len=:), allocatable :: tScalar_asFormattedString + class (tScalar), intent(in), target :: self + integer, intent(in), optional :: indent + + tScalar_asFormattedString = trim(self%value)//IO_EOL + +end function tScalar_asFormattedString + + +!-------------------------------------------------------------------------------------------------- +!> @brief List as string (YAML block style) +!-------------------------------------------------------------------------------------------------- +recursive function tList_asFormattedString(self,indent) result(str) + + class (tList),intent(in),target :: self + integer, intent(in),optional :: indent + + type (tItem), pointer :: item + character(len=:), allocatable :: str + integer :: i, indent_ + + str = '' + if(present(indent)) then + indent_ = indent + else + indent_ = 0 + endif + + item => self%first + do i = 1, self%length + if(i /= 1) str = str//repeat(' ',indent_) + str = str//'- '//item%node%asFormattedString(indent_+2) + item => item%next + end do + +end function tList_asFormattedString + + +!-------------------------------------------------------------------------------------------------- +!> @brief Dictionary as string (YAML block style) +!-------------------------------------------------------------------------------------------------- +recursive function tDict_asFormattedString(self,indent) result(str) + + class (tDict),intent(in),target :: self + integer, intent(in),optional :: indent + + type (tItem),pointer :: item + character(len=:), allocatable :: str + integer :: i, indent_ + + str = '' + if(present(indent)) then + indent_ = indent + else + indent_ = 0 + endif + + item => self%first + do i = 1, self%length + if(i /= 1) str = str//repeat(' ',indent_) + select type(node_1 =>item%node) + class is(tScalar) + str = str//trim(item%key)//': '//item%node%asFormattedString(indent_+len_trim(item%key)+2) + class default + str = str//trim(item%key)//':'//IO_EOL//repeat(' ',indent_+2)//item%node%asFormattedString(indent_+2) + endselect + item => item%next + end do + +end function tDict_asFormattedString + + +!-------------------------------------------------------------------------------------------------- +!> @brief Convert to float +!-------------------------------------------------------------------------------------------------- +function tScalar_asFloat(self) + + class(tScalar), intent(in), target :: self + real(pReal) :: tScalar_asFloat + + tScalar_asFloat = IO_stringAsFloat(self%value) + +end function tScalar_asFloat + + +!-------------------------------------------------------------------------------------------------- +!> @brief Convert to int +!-------------------------------------------------------------------------------------------------- +function tScalar_asInt(self) + + class(tScalar), intent(in), target :: self + integer :: tScalar_asInt + + tScalar_asInt = IO_stringAsInt(self%value) + +end function tScalar_asInt + + +!-------------------------------------------------------------------------------------------------- +!> @brief Convert to bool +!-------------------------------------------------------------------------------------------------- +function tScalar_asBool(self) + + class(tScalar), intent(in), target :: self + logical :: tScalar_asBool + + tScalar_asBool = IO_stringAsBool(self%value) + +end function tScalar_asBool + + +!-------------------------------------------------------------------------------------------------- +!> @brief Convert to string +!-------------------------------------------------------------------------------------------------- +function tScalar_asString(self) + + class(tScalar), intent(in), target :: self + character(len=:), allocatable :: tScalar_asString + + tScalar_asString = self%value + +end function tScalar_asString + + +!-------------------------------------------------------------------------------------------------- +!> @brief Convert to float array +!-------------------------------------------------------------------------------------------------- +function tList_asFloats(self) + + class(tList), intent(in), target :: self + real(pReal), dimension(:), allocatable :: tList_asFloats + + integer :: i + type(tItem), pointer :: item + type(tScalar), pointer :: scalar + + allocate(tList_asFloats(self%length)) + item => self%first + do i = 1, self%length + scalar => item%node%asScalar() + tList_asFloats(i) = scalar%asFloat() + item => item%next + enddo + +end function tList_asFloats + + +!-------------------------------------------------------------------------------------------------- +!> @brief Convert to int array +!-------------------------------------------------------------------------------------------------- +function tList_asInts(self) + + class(tList), intent(in), target :: self + integer, dimension(:), allocatable :: tList_asInts + + integer :: i + type(tItem), pointer :: item + type(tScalar), pointer :: scalar + + allocate(tList_asInts(self%length)) + item => self%first + do i = 1, self%length + scalar => item%node%asScalar() + tList_asInts(i) = scalar%asInt() + item => item%next + enddo + +end function tList_asInts + + +!-------------------------------------------------------------------------------------------------- +!> @brief Convert to bool array +!-------------------------------------------------------------------------------------------------- +function tList_asBools(self) + + class(tList), intent(in), target :: self + logical, dimension(:), allocatable :: tList_asBools + + integer :: i + type(tItem), pointer :: item + type(tScalar), pointer :: scalar + + allocate(tList_asBools(self%length)) + item => self%first + do i = 1, self%length + scalar => item%node%asScalar() + tList_asBools(i) = scalar%asBool() + item => item%next + enddo + +end function tList_asBools + + +!-------------------------------------------------------------------------------------------------- +!> @brief Convert to string array +!-------------------------------------------------------------------------------------------------- +function tList_asStrings(self) + + class(tList), intent(in), target :: self + character(len=:), allocatable, dimension(:) :: tList_asStrings + + integer :: i,len_max + type(tItem), pointer :: item + type(tScalar), pointer :: scalar + + len_max = 0 + allocate(character(len=pStringLen) :: tList_asStrings(self%length)) + item => self%first + do i = 1, self%length + scalar => item%node%asScalar() + tList_asStrings(i) = scalar%asString() + len_max = max(len_max, len_trim(tList_asStrings(i))) + item => item%next + enddo + + !ToDo: trim to len_max + +end function tList_asStrings + + +!-------------------------------------------------------------------------------------------------- +!> @brief Append element +!-------------------------------------------------------------------------------------------------- +subroutine tList_append(self,node) + + class(tList), intent(inout) :: self + class(tNode), intent(in), target :: node + + type(tItem), pointer :: item + + if (.not. associated(self%first)) then + allocate(self%first) + item => self%first + else + item => self%first + do while (associated(item%next)) + item => item%next + end do + allocate(item%next) + item => item%next + end if + + item%node => node + self%length = self%length + 1 + +end subroutine tList_append + + +!-------------------------------------------------------------------------------------------------- +!> @brief Set the value of a key (either replace or add new) +!-------------------------------------------------------------------------------------------------- +subroutine tDict_set(self,key,node) + + class (tDict), intent(inout) :: self + character(len=*), intent(in) :: key + class(tNode), intent(in), target :: node + + type(tItem), pointer :: item + + if (.not. associated(self%first)) then + allocate(self%first) + item => self%first + self%length = 1 + else + item => self%first + searchExisting: do while (associated(item%next)) + if (item%key == key) exit + item => item%next + end do searchExisting + if (.not. item%key == key) then + allocate(item%next) + item => item%next + self%length = self%length + 1 + end if + end if + + item%key = key + item%node => node + +end subroutine tDict_set + + +!-------------------------------------------------------------------------------------------------- +!> @brief empties lists and dicts and free associated memory +!> @details called when variable goes out of scope. +!-------------------------------------------------------------------------------------------------- +recursive subroutine tList_finalize(self) + + type (tList),intent(inout) :: self + + deallocate(self%first) + +end subroutine tList_finalize + + +!-------------------------------------------------------------------------------------------------- +!> @brief empties nodes and frees associated memory +!-------------------------------------------------------------------------------------------------- +recursive subroutine tItem_finalize(self) + + type(tItem),intent(inout) :: self + + deallocate(self%node) + if(associated(self%next)) deallocate(self%next) + +end subroutine tItem_finalize + +end module YAML_types diff --git a/src/commercialFEM_fileList.f90 b/src/commercialFEM_fileList.f90 index 64ad3e1d7..f8a5ad96a 100644 --- a/src/commercialFEM_fileList.f90 +++ b/src/commercialFEM_fileList.f90 @@ -7,12 +7,12 @@ #include "numerics.f90" #include "debug.f90" #include "list.f90" +#include "YAML_types.f90" #include "future.f90" #include "config.f90" #include "LAPACK_interface.f90" #include "math.f90" #include "quaternions.f90" -#include "Lambert.f90" #include "rotations.f90" #include "FEsolving.f90" #include "element.f90" diff --git a/src/constitutive_plastic_nonlocal.f90 b/src/constitutive_plastic_nonlocal.f90 index 77f1556f6..65fea24b2 100644 --- a/src/constitutive_plastic_nonlocal.f90 +++ b/src/constitutive_plastic_nonlocal.f90 @@ -962,39 +962,24 @@ module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature,timestep, & integer :: & ph, & - neighbor_instance, & !< instance of my neighbor's plasticity ns, & !< short notation for the total number of active slip systems c, & !< character of dislocation - n, & !< index of my current neighbor - neighbor_el, & !< element number of my neighbor - neighbor_ip, & !< integration point of my neighbor - neighbor_n, & !< neighbor index pointing to me when looking from my neighbor - opposite_neighbor, & !< index of my opposite neighbor - opposite_ip, & !< ip of my opposite neighbor - opposite_el, & !< element index of my opposite neighbor - opposite_n, & !< neighbor index pointing to me when looking from my opposite neighbor t, & !< type of dislocation - no,& !< neighbor offset shortcut - np,& !< neighbor phase shortcut - topp, & !< type of dislocation with opposite sign to t s !< index of my current slip system real(pReal), dimension(param(instance)%sum_N_sl,10) :: & rho, & rho0, & !< dislocation density at beginning of time step rhoDot, & !< density evolution rhoDotMultiplication, & !< density evolution by multiplication - rhoDotFlux, & !< density evolution by flux rhoDotSingle2DipoleGlide, & !< density evolution by dipole formation (by glide) rhoDotAthermalAnnihilation, & !< density evolution by athermal annihilation rhoDotThermalAnnihilation !< density evolution by thermal annihilation real(pReal), dimension(param(instance)%sum_N_sl,8) :: & rhoSgl, & !< current single dislocation densities (positive/negative screw and edge without dipoles) - neighbor_rhoSgl0, & !< current single dislocation densities of neighboring ip (positive/negative screw and edge without dipoles) my_rhoSgl0 !< single dislocation densities of central ip (positive/negative screw and edge without dipoles) real(pReal), dimension(param(instance)%sum_N_sl,4) :: & v, & !< current dislocation glide velocity v0, & - neighbor_v0, & !< dislocation glide velocity of enighboring ip gdot !< shear rates real(pReal), dimension(param(instance)%sum_N_sl) :: & tau, & !< current resolved shear stress @@ -1003,23 +988,7 @@ module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature,timestep, & rhoDip, & !< current dipole dislocation densities (screw and edge dipoles) dLower, & !< minimum stable dipole distance for edges and screws dUpper !< current maximum stable dipole distance for edges and screws - real(pReal), dimension(3,param(instance)%sum_N_sl,4) :: & - m !< direction of dislocation motion - real(pReal), dimension(3,3) :: & - my_F, & !< my total deformation gradient - neighbor_F, & !< total deformation gradient of my neighbor - my_Fe, & !< my elastic deformation gradient - neighbor_Fe, & !< elastic deformation gradient of my neighbor - Favg !< average total deformation gradient of me and my neighbor - real(pReal), dimension(3) :: & - normal_neighbor2me, & !< interface normal pointing from my neighbor to me in neighbor's lattice configuration - normal_neighbor2me_defConf, & !< interface normal pointing from my neighbor to me in shared deformed configuration - normal_me2neighbor, & !< interface normal pointing from me to my neighbor in my lattice configuration - normal_me2neighbor_defConf !< interface normal pointing from me to my neighbor in shared deformed configuration real(pReal) :: & - area, & !< area of the current interface - transmissivity, & !< overall transmissivity of dislocation flux to neighboring material point - lineLength, & !< dislocation line length leaving the current interface selfDiffusion !< self diffusion ph = material_phaseAt(1,el) @@ -1092,6 +1061,172 @@ module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature,timestep, & * sqrt(stt%rho_forest(:,of)) / prm%lambda0 / prm%burgers, 2, 4) endif isBCC + forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,instance),of) + + + !**************************************************************************** + !*** calculate dipole formation and annihilation + + !*** formation by glide + do c = 1,2 + rhoDotSingle2DipoleGlide(:,2*c-1) = -2.0_pReal * dUpper(:,c) / prm%burgers & + * ( rhoSgl(:,2*c-1) * abs(gdot(:,2*c)) & ! negative mobile --> positive mobile + + rhoSgl(:,2*c) * abs(gdot(:,2*c-1)) & ! positive mobile --> negative mobile + + abs(rhoSgl(:,2*c+4)) * abs(gdot(:,2*c-1))) ! positive mobile --> negative immobile + + rhoDotSingle2DipoleGlide(:,2*c) = -2.0_pReal * dUpper(:,c) / prm%burgers & + * ( rhoSgl(:,2*c-1) * abs(gdot(:,2*c)) & ! negative mobile --> positive mobile + + rhoSgl(:,2*c) * abs(gdot(:,2*c-1)) & ! positive mobile --> negative mobile + + abs(rhoSgl(:,2*c+3)) * abs(gdot(:,2*c))) ! negative mobile --> positive immobile + + rhoDotSingle2DipoleGlide(:,2*c+3) = -2.0_pReal * dUpper(:,c) / prm%burgers & + * rhoSgl(:,2*c+3) * abs(gdot(:,2*c)) ! negative mobile --> positive immobile + + rhoDotSingle2DipoleGlide(:,2*c+4) = -2.0_pReal * dUpper(:,c) / prm%burgers & + * rhoSgl(:,2*c+4) * abs(gdot(:,2*c-1)) ! positive mobile --> negative immobile + + rhoDotSingle2DipoleGlide(:,c+8) = abs(rhoDotSingle2DipoleGlide(:,2*c+3)) & + + abs(rhoDotSingle2DipoleGlide(:,2*c+4)) & + - rhoDotSingle2DipoleGlide(:,2*c-1) & + - rhoDotSingle2DipoleGlide(:,2*c) + enddo + + + !*** athermal annihilation + rhoDotAthermalAnnihilation = 0.0_pReal + forall (c=1:2) & + rhoDotAthermalAnnihilation(:,c+8) = -2.0_pReal * dLower(:,c) / prm%burgers & + * ( 2.0_pReal * (rhoSgl(:,2*c-1) * abs(gdot(:,2*c)) + rhoSgl(:,2*c) * abs(gdot(:,2*c-1))) & ! was single hitting single + + 2.0_pReal * (abs(rhoSgl(:,2*c+3)) * abs(gdot(:,2*c)) + abs(rhoSgl(:,2*c+4)) * abs(gdot(:,2*c-1))) & ! was single hitting immobile single or was immobile single hit by single + + rhoDip(:,c) * (abs(gdot(:,2*c-1)) + abs(gdot(:,2*c)))) ! single knocks dipole constituent + + ! annihilated screw dipoles leave edge jogs behind on the colinear system + if (lattice_structure(ph) == LATTICE_fcc_ID) & + forall (s = 1:ns, prm%colinearSystem(s) > 0) & + rhoDotAthermalAnnihilation(prm%colinearSystem(s),1:2) = - rhoDotAthermalAnnihilation(s,10) & + * 0.25_pReal * sqrt(stt%rho_forest(s,of)) * (dUpper(s,2) + dLower(s,2)) * prm%edgeJogFactor + + + !*** thermally activated annihilation of edge dipoles by climb + rhoDotThermalAnnihilation = 0.0_pReal + selfDiffusion = prm%Dsd0 * exp(-prm%selfDiffusionEnergy / (kB * Temperature)) + vClimb = prm%atomicVolume * selfDiffusion * prm%mu & + / ( kB * Temperature * PI * (1.0_pReal-prm%nu) * (dUpper(:,1) + dLower(:,1))) + forall (s = 1:ns, dUpper(s,1) > dLower(s,1)) & + rhoDotThermalAnnihilation(s,9) = max(- 4.0_pReal * rhoDip(s,1) * vClimb(s) / (dUpper(s,1) - dLower(s,1)), & + - rhoDip(s,1) / timestep - rhoDotAthermalAnnihilation(s,9) & + - rhoDotSingle2DipoleGlide(s,9)) ! make sure that we do not annihilate more dipoles than we have + + rhoDot = rhoDotFlux(F,Fp,timestep, instance,of,ip,el) & + + rhoDotMultiplication & + + rhoDotSingle2DipoleGlide & + + rhoDotAthermalAnnihilation & + + rhoDotThermalAnnihilation + + + if ( any(rho(:,mob) + rhoDot(:,1:4) * timestep < -prm%atol_rho) & + .or. any(rho(:,dip) + rhoDot(:,9:10) * timestep < -prm%atol_rho)) then +#ifdef DEBUG + if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0) then + write(6,'(a,i5,a,i2)') '<< CONST >> evolution rate leads to negative density at el ',el,' ip ',ip + write(6,'(a)') '<< CONST >> enforcing cutback !!!' + endif +#endif + plasticState(ph)%dotState = IEEE_value(1.0_pReal,IEEE_quiet_NaN) + else + dot%rho(:,of) = pack(rhoDot,.true.) + dot%gamma(:,of) = sum(gdot,2) + endif + + end associate + +end subroutine plastic_nonlocal_dotState + + +!--------------------------------------------------------------------------------------------------- +!> @brief calculates the rate of change of microstructure +!--------------------------------------------------------------------------------------------------- +function rhoDotFlux(F,Fp,timestep, instance,of,ip,el) + + real(pReal), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem), intent(in) :: & + F, & !< elastic deformation gradient + Fp !< plastic deformation gradient + real(pReal), intent(in) :: & + timestep !< substepped crystallite time increment + integer, intent(in) :: & + instance, & + of, & + ip, & !< current integration point + el !< current element number + + integer :: & + ph, & + neighbor_instance, & !< instance of my neighbor's plasticity + ns, & !< short notation for the total number of active slip systems + c, & !< character of dislocation + n, & !< index of my current neighbor + neighbor_el, & !< element number of my neighbor + neighbor_ip, & !< integration point of my neighbor + neighbor_n, & !< neighbor index pointing to me when looking from my neighbor + opposite_neighbor, & !< index of my opposite neighbor + opposite_ip, & !< ip of my opposite neighbor + opposite_el, & !< element index of my opposite neighbor + opposite_n, & !< neighbor index pointing to me when looking from my opposite neighbor + t, & !< type of dislocation + no,& !< neighbor offset shortcut + np,& !< neighbor phase shortcut + topp, & !< type of dislocation with opposite sign to t + s !< index of my current slip system + real(pReal), dimension(param(instance)%sum_N_sl,10) :: & + rho, & + rho0, & !< dislocation density at beginning of time step + rhoDotFlux !< density evolution by flux + real(pReal), dimension(param(instance)%sum_N_sl,8) :: & + rhoSgl, & !< current single dislocation densities (positive/negative screw and edge without dipoles) + neighbor_rhoSgl0, & !< current single dislocation densities of neighboring ip (positive/negative screw and edge without dipoles) + my_rhoSgl0 !< single dislocation densities of central ip (positive/negative screw and edge without dipoles) + real(pReal), dimension(param(instance)%sum_N_sl,4) :: & + v, & !< current dislocation glide velocity + v0, & + neighbor_v0, & !< dislocation glide velocity of enighboring ip + gdot !< shear rates + real(pReal), dimension(3,param(instance)%sum_N_sl,4) :: & + m !< direction of dislocation motion + real(pReal), dimension(3,3) :: & + my_F, & !< my total deformation gradient + neighbor_F, & !< total deformation gradient of my neighbor + my_Fe, & !< my elastic deformation gradient + neighbor_Fe, & !< elastic deformation gradient of my neighbor + Favg !< average total deformation gradient of me and my neighbor + real(pReal), dimension(3) :: & + normal_neighbor2me, & !< interface normal pointing from my neighbor to me in neighbor's lattice configuration + normal_neighbor2me_defConf, & !< interface normal pointing from my neighbor to me in shared deformed configuration + normal_me2neighbor, & !< interface normal pointing from me to my neighbor in my lattice configuration + normal_me2neighbor_defConf !< interface normal pointing from me to my neighbor in shared deformed configuration + real(pReal) :: & + area, & !< area of the current interface + transmissivity, & !< overall transmissivity of dislocation flux to neighboring material point + lineLength !< dislocation line length leaving the current interface + + ph = material_phaseAt(1,el) + + associate(prm => param(instance), & + dst => microstructure(instance), & + dot => dotState(instance), & + stt => state(instance)) + ns = prm%sum_N_sl + + gdot = 0.0_pReal + + rho = getRho(instance,of,ip,el) + rhoSgl = rho(:,sgl) + rho0 = getRho0(instance,of,ip,el) + my_rhoSgl0 = rho0(:,sgl) + + forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,instance),of) !ToDo: MD: I think we should use state0 here + gdot = rhoSgl(:,1:4) * v * spread(prm%burgers,2,4) + + forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,instance),of) !**************************************************************************** @@ -1114,7 +1249,7 @@ module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature,timestep, & write(6,'(a)') '<< CONST >> enforcing cutback !!!' endif #endif - plasticState(ph)%dotState = IEEE_value(1.0_pReal,IEEE_quiet_NaN) ! -> return NaN and, hence, enforce cutback + rhoDotFlux = IEEE_value(1.0_pReal,IEEE_quiet_NaN) ! enforce cutback return endif @@ -1240,108 +1375,9 @@ module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature,timestep, & enddo neighbors endif - - - !**************************************************************************** - !*** calculate dipole formation and annihilation - - !*** formation by glide - do c = 1,2 - rhoDotSingle2DipoleGlide(:,2*c-1) = -2.0_pReal * dUpper(:,c) / prm%burgers & - * ( rhoSgl(:,2*c-1) * abs(gdot(:,2*c)) & ! negative mobile --> positive mobile - + rhoSgl(:,2*c) * abs(gdot(:,2*c-1)) & ! positive mobile --> negative mobile - + abs(rhoSgl(:,2*c+4)) * abs(gdot(:,2*c-1))) ! positive mobile --> negative immobile - - rhoDotSingle2DipoleGlide(:,2*c) = -2.0_pReal * dUpper(:,c) / prm%burgers & - * ( rhoSgl(:,2*c-1) * abs(gdot(:,2*c)) & ! negative mobile --> positive mobile - + rhoSgl(:,2*c) * abs(gdot(:,2*c-1)) & ! positive mobile --> negative mobile - + abs(rhoSgl(:,2*c+3)) * abs(gdot(:,2*c))) ! negative mobile --> positive immobile - - rhoDotSingle2DipoleGlide(:,2*c+3) = -2.0_pReal * dUpper(:,c) / prm%burgers & - * rhoSgl(:,2*c+3) * abs(gdot(:,2*c)) ! negative mobile --> positive immobile - - rhoDotSingle2DipoleGlide(:,2*c+4) = -2.0_pReal * dUpper(:,c) / prm%burgers & - * rhoSgl(:,2*c+4) * abs(gdot(:,2*c-1)) ! positive mobile --> negative immobile - - rhoDotSingle2DipoleGlide(:,c+8) = abs(rhoDotSingle2DipoleGlide(:,2*c+3)) & - + abs(rhoDotSingle2DipoleGlide(:,2*c+4)) & - - rhoDotSingle2DipoleGlide(:,2*c-1) & - - rhoDotSingle2DipoleGlide(:,2*c) - enddo - - - !*** athermal annihilation - rhoDotAthermalAnnihilation = 0.0_pReal - forall (c=1:2) & - rhoDotAthermalAnnihilation(:,c+8) = -2.0_pReal * dLower(:,c) / prm%burgers & - * ( 2.0_pReal * (rhoSgl(:,2*c-1) * abs(gdot(:,2*c)) + rhoSgl(:,2*c) * abs(gdot(:,2*c-1))) & ! was single hitting single - + 2.0_pReal * (abs(rhoSgl(:,2*c+3)) * abs(gdot(:,2*c)) + abs(rhoSgl(:,2*c+4)) * abs(gdot(:,2*c-1))) & ! was single hitting immobile single or was immobile single hit by single - + rhoDip(:,c) * (abs(gdot(:,2*c-1)) + abs(gdot(:,2*c)))) ! single knocks dipole constituent - - ! annihilated screw dipoles leave edge jogs behind on the colinear system - if (lattice_structure(ph) == LATTICE_fcc_ID) & - forall (s = 1:ns, prm%colinearSystem(s) > 0) & - rhoDotAthermalAnnihilation(prm%colinearSystem(s),1:2) = - rhoDotAthermalAnnihilation(s,10) & - * 0.25_pReal * sqrt(stt%rho_forest(s,of)) * (dUpper(s,2) + dLower(s,2)) * prm%edgeJogFactor - - - !*** thermally activated annihilation of edge dipoles by climb - rhoDotThermalAnnihilation = 0.0_pReal - selfDiffusion = prm%Dsd0 * exp(-prm%selfDiffusionEnergy / (kB * Temperature)) - vClimb = prm%atomicVolume * selfDiffusion * prm%mu & - / ( kB * Temperature * PI * (1.0_pReal-prm%nu) * (dUpper(:,1) + dLower(:,1))) - forall (s = 1:ns, dUpper(s,1) > dLower(s,1)) & - rhoDotThermalAnnihilation(s,9) = max(- 4.0_pReal * rhoDip(s,1) * vClimb(s) / (dUpper(s,1) - dLower(s,1)), & - - rhoDip(s,1) / timestep - rhoDotAthermalAnnihilation(s,9) & - - rhoDotSingle2DipoleGlide(s,9)) ! make sure that we do not annihilate more dipoles than we have - - rhoDot = rhoDotFlux & - + rhoDotMultiplication & - + rhoDotSingle2DipoleGlide & - + rhoDotAthermalAnnihilation & - + rhoDotThermalAnnihilation - -#ifdef DEBUG - if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0 & - .and. ((debug_e == el .and. debug_i == ip)& - .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0 )) then - write(6,'(a,/,4(12x,12(e12.5,1x),/))') '<< CONST >> dislocation multiplication', & - rhoDotMultiplication(:,1:4) * timestep - write(6,'(a,/,8(12x,12(e12.5,1x),/))') '<< CONST >> dislocation flux', & - rhoDotFlux(:,1:8) * timestep - write(6,'(a,/,10(12x,12(e12.5,1x),/))') '<< CONST >> dipole formation by glide', & - rhoDotSingle2DipoleGlide * timestep - write(6,'(a,/,10(12x,12(e12.5,1x),/))') '<< CONST >> athermal dipole annihilation', & - rhoDotAthermalAnnihilation * timestep - write(6,'(a,/,2(12x,12(e12.5,1x),/))') '<< CONST >> thermally activated dipole annihilation', & - rhoDotThermalAnnihilation(:,9:10) * timestep - write(6,'(a,/,10(12x,12(e12.5,1x),/))') '<< CONST >> total density change', & - rhoDot * timestep - write(6,'(a,/,10(12x,12(f12.5,1x),/))') '<< CONST >> relative density change', & - rhoDot(:,1:8) * timestep / (abs(stt%rho(:,sgl))+1.0e-10), & - rhoDot(:,9:10) * timestep / (stt%rho(:,dip)+1.0e-10) - write(6,*) - endif -#endif - - - if ( any(rho(:,mob) + rhoDot(:,1:4) * timestep < -prm%atol_rho) & - .or. any(rho(:,dip) + rhoDot(:,9:10) * timestep < -prm%atol_rho)) then -#ifdef DEBUG - if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0) then - write(6,'(a,i5,a,i2)') '<< CONST >> evolution rate leads to negative density at el ',el,' ip ',ip - write(6,'(a)') '<< CONST >> enforcing cutback !!!' - endif -#endif - plasticState(ph)%dotState = IEEE_value(1.0_pReal,IEEE_quiet_NaN) - else - dot%rho(:,of) = pack(rhoDot,.true.) - dot%gamma(:,of) = sum(gdot,2) - endif - end associate -end subroutine plastic_nonlocal_dotState +end function rhoDotFlux !-------------------------------------------------------------------------------------------------- diff --git a/src/grid/grid_mech_FEM.f90 b/src/grid/grid_mech_FEM.f90 index 9cdf2f5b5..c8bef7ff2 100644 --- a/src/grid/grid_mech_FEM.f90 +++ b/src/grid/grid_mech_FEM.f90 @@ -132,11 +132,11 @@ subroutine grid_mech_FEM_init [grid(1)],[grid(2)],localK, & mech_grid,ierr) CHKERRQ(ierr) - call DMDASetUniformCoordinates(mech_grid,0.0_pReal,geomSize(1),0.0_pReal,geomSize(2),0.0_pReal,geomSize(3),ierr) - CHKERRQ(ierr) call SNESSetDM(mech_snes,mech_grid,ierr); CHKERRQ(ierr) call DMsetFromOptions(mech_grid,ierr); CHKERRQ(ierr) call DMsetUp(mech_grid,ierr); CHKERRQ(ierr) + call DMDASetUniformCoordinates(mech_grid,0.0_pReal,geomSize(1),0.0_pReal,geomSize(2),0.0_pReal,geomSize(3),ierr) + CHKERRQ(ierr) call DMCreateGlobalVector(mech_grid,solution_current,ierr); CHKERRQ(ierr) call DMCreateGlobalVector(mech_grid,solution_lastInc,ierr); CHKERRQ(ierr) call DMCreateGlobalVector(mech_grid,solution_rate ,ierr); CHKERRQ(ierr) diff --git a/src/rotations.f90 b/src/rotations.f90 index c20b9e905..fea642170 100644 --- a/src/rotations.f90 +++ b/src/rotations.f90 @@ -3,27 +3,27 @@ ! Modified 2017-2020, Martin Diehl/Max-Planck-Institut für Eisenforschung GmbH ! All rights reserved. ! -! Redistribution and use in source and binary forms, with or without modification, are +! Redistribution and use in source and binary forms, with or without modification, are ! permitted provided that the following conditions are met: ! -! - Redistributions of source code must retain the above copyright notice, this list +! - Redistributions of source code must retain the above copyright notice, this list ! of conditions and the following disclaimer. -! - Redistributions in binary form must reproduce the above copyright notice, this -! list of conditions and the following disclaimer in the documentation and/or +! - Redistributions in binary form must reproduce the above copyright notice, this +! list of conditions and the following disclaimer in the documentation and/or ! other materials provided with the distribution. -! - Neither the names of Marc De Graef, Carnegie Mellon University nor the names -! of its contributors may be used to endorse or promote products derived from +! - Neither the names of Marc De Graef, Carnegie Mellon University nor the names +! of its contributors may be used to endorse or promote products derived from ! this software without specific prior written permission. ! -! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" -! AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE -! IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE -! ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE -! LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL -! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR -! SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER -! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, -! OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +! AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +! IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +! ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +! LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +! SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, +! OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE ! USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ! ################################################################### @@ -31,7 +31,7 @@ !> @author Marc De Graef, Carnegie Mellon University !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @brief rotation storage and conversion -!> @details: rotation is internally stored as quaternion. It can be inialized from different +!> @details: rotation is internally stored as quaternion. It can be inialized from different !> representations and also returns itself in different representations. ! ! All methods and naming conventions based on Rowenhorst_etal2015 @@ -50,9 +50,8 @@ module rotations use prec use IO use math - use Lambert use quaternions - + implicit none private @@ -80,7 +79,19 @@ module rotations procedure, public :: misorientation procedure, public :: standardize end type rotation - + + real(pReal), parameter :: & + SPI = sqrt(PI), & + PREF = sqrt(6.0_pReal/PI), & + A = PI**(5.0_pReal/6.0_pReal)/6.0_pReal**(1.0_pReal/6.0_pReal), & + AP = PI**(2.0_pReal/3.0_pReal), & + SC = A/AP, & + BETA = A/2.0_pReal, & + R1 = (3.0_pReal*PI/4.0_pReal)**(1.0_pReal/3.0_pReal), & + R2 = sqrt(2.0_pReal), & + PI12 = PI/12.0_pReal, & + PREK = R1 * 2.0_pReal**(1.0_pReal/4.0_pReal)/BETA + public :: & rotations_init, & eu2om @@ -106,16 +117,16 @@ pure function asQuaternion(self) class(rotation), intent(in) :: self real(pReal), dimension(4) :: asQuaternion - + asQuaternion = self%q%asArray() end function asQuaternion !--------------------------------------------------------------------------------------------------- pure function asEulers(self) - + class(rotation), intent(in) :: self real(pReal), dimension(3) :: asEulers - + asEulers = qu2eu(self%q%asArray()) end function asEulers @@ -124,16 +135,16 @@ pure function asAxisAngle(self) class(rotation), intent(in) :: self real(pReal), dimension(4) :: asAxisAngle - + asAxisAngle = qu2ax(self%q%asArray()) end function asAxisAngle !--------------------------------------------------------------------------------------------------- pure function asMatrix(self) - + class(rotation), intent(in) :: self real(pReal), dimension(3,3) :: asMatrix - + asMatrix = qu2om(self%q%asArray()) end function asMatrix @@ -142,20 +153,20 @@ pure function asRodrigues(self) class(rotation), intent(in) :: self real(pReal), dimension(4) :: asRodrigues - + asRodrigues = qu2ro(self%q%asArray()) - + end function asRodrigues !--------------------------------------------------------------------------------------------------- pure function asHomochoric(self) class(rotation), intent(in) :: self real(pReal), dimension(3) :: asHomochoric - + asHomochoric = qu2ho(self%q%asArray()) end function asHomochoric - + !--------------------------------------------------------------------------------------------------- ! Initialize rotation from different representations !--------------------------------------------------------------------------------------------------- @@ -207,7 +218,7 @@ subroutine fromAxisAngle(self,ax,degrees,P) else angle = merge(ax(4)*INRAD,ax(4),degrees) endif - + if (.not. present(P)) then axis = ax(1:3) else @@ -217,7 +228,7 @@ subroutine fromAxisAngle(self,ax,degrees,P) if(dNeq(norm2(axis),1.0_pReal) .or. angle < 0.0_pReal .or. angle > PI) & call IO_error(402,ext_msg='fromAxisAngle') - + self%q = ax2qu([axis,angle]) end subroutine fromAxisAngle @@ -240,10 +251,10 @@ end subroutine fromMatrix !> @brief: Rotate a rotation !--------------------------------------------------------------------------------------------------- pure elemental function rotRot__(self,R) result(rRot) - + type(rotation) :: rRot class(rotation), intent(in) :: self,R - + rRot = rotation(self%q*R%q) call rRot%standardize() @@ -251,12 +262,12 @@ end function rotRot__ !--------------------------------------------------------------------------------------------------- -!> @brief quaternion representation with positive q +!> @brief quaternion representation with positive q !--------------------------------------------------------------------------------------------------- pure elemental subroutine standardize(self) class(rotation), intent(inout) :: self - + if (real(self%q) < 0.0_pReal) self%q = self%q%homomorphed() end subroutine standardize @@ -267,22 +278,22 @@ end subroutine standardize !> @brief rotate a vector passively (default) or actively !--------------------------------------------------------------------------------------------------- pure function rotVector(self,v,active) result(vRot) - + real(pReal), dimension(3) :: vRot class(rotation), intent(in) :: self real(pReal), intent(in), dimension(3) :: v logical, intent(in), optional :: active - + real(pReal), dimension(3) :: v_normed type(quaternion) :: q logical :: passive - + if (present(active)) then passive = .not. active else passive = .true. endif - + if (dEq0(norm2(v))) then vRot = v else @@ -304,12 +315,12 @@ end function rotVector !> @details: rotation is based on rotation matrix !--------------------------------------------------------------------------------------------------- pure function rotTensor2(self,T,active) result(tRot) - + real(pReal), dimension(3,3) :: tRot class(rotation), intent(in) :: self real(pReal), intent(in), dimension(3,3) :: T logical, intent(in), optional :: active - + logical :: passive if (present(active)) then @@ -317,7 +328,7 @@ pure function rotTensor2(self,T,active) result(tRot) else passive = .true. endif - + if (passive) then tRot = matmul(matmul(self%asMatrix(),T),transpose(self%asMatrix())) else @@ -339,7 +350,7 @@ pure function rotTensor4(self,T,active) result(tRot) class(rotation), intent(in) :: self real(pReal), intent(in), dimension(3,3,3,3) :: T logical, intent(in), optional :: active - + real(pReal), dimension(3,3) :: R integer :: i,j,k,l,m,n,o,p @@ -370,7 +381,7 @@ pure function rotTensor4sym(self,T,active) result(tRot) class(rotation), intent(in) :: self real(pReal), intent(in), dimension(6,6) :: T logical, intent(in), optional :: active - + if (present(active)) then tRot = math_sym3333to66(rotTensor4(self,math_66toSym3333(T),active)) else @@ -384,10 +395,10 @@ end function rotTensor4sym !> @brief misorientation !--------------------------------------------------------------------------------------------------- pure elemental function misorientation(self,other) - + type(rotation) :: misorientation class(rotation), intent(in) :: self, other - + misorientation%q = other%q * conjg(self%q) end function misorientation @@ -401,7 +412,7 @@ pure function qu2om(qu) result(om) real(pReal), intent(in), dimension(4) :: qu real(pReal), dimension(3,3) :: om - + real(pReal) :: qq qq = qu(1)**2-sum(qu(2:4)**2) @@ -431,13 +442,13 @@ pure function qu2eu(qu) result(eu) real(pReal), intent(in), dimension(4) :: qu real(pReal), dimension(3) :: eu - + real(pReal) :: q12, q03, chi - + q03 = qu(1)**2+qu(4)**2 q12 = qu(2)**2+qu(3)**2 chi = sqrt(q03*q12) - + degenerated: if (dEq0(q12)) then eu = [atan2(-P*2.0_pReal*qu(1)*qu(4),qu(1)**2-qu(4)**2), 0.0_pReal, 0.0_pReal] elseif (dEq0(q03)) then @@ -460,7 +471,7 @@ pure function qu2ax(qu) result(ax) real(pReal), intent(in), dimension(4) :: qu real(pReal), dimension(4) :: ax - + real(pReal) :: omega, s if (dEq0(sum(qu(2:4)**2))) then @@ -481,13 +492,13 @@ end function qu2ax !> @brief convert unit quaternion to Rodrigues vector !--------------------------------------------------------------------------------------------------- pure function qu2ro(qu) result(ro) - + real(pReal), intent(in), dimension(4) :: qu real(pReal), dimension(4) :: ro - + real(pReal) :: s real(pReal), parameter :: thr = 1.0e-8_pReal - + if (abs(qu(1)) < thr) then ro = [qu(2), qu(3), qu(4), IEEE_value(1.0_pReal,IEEE_positive_inf)] else @@ -497,7 +508,7 @@ pure function qu2ro(qu) result(ro) else ro = [qu(2)/s,qu(3)/s,qu(4)/s, tan(acos(math_clip(qu(1),-1.0_pReal,1.0_pReal)))] endif - + end if end function qu2ro @@ -511,12 +522,12 @@ pure function qu2ho(qu) result(ho) real(pReal), intent(in), dimension(4) :: qu real(pReal), dimension(3) :: ho - + real(pReal) :: omega, f omega = 2.0 * acos(math_clip(qu(1),-1.0_pReal,1.0_pReal)) - - if (dEq0(omega)) then + + if (dEq0(omega,tol=1.e-5_pReal)) then ho = [ 0.0_pReal, 0.0_pReal, 0.0_pReal ] else ho = qu(2:4) @@ -532,7 +543,7 @@ end function qu2ho !> @brief convert unit quaternion to cubochoric !--------------------------------------------------------------------------------------------------- pure function qu2cu(qu) result(cu) - + real(pReal), intent(in), dimension(4) :: qu real(pReal), dimension(3) :: cu @@ -565,18 +576,18 @@ pure function om2eu(om) result(eu) real(pReal), intent(in), dimension(3,3) :: om real(pReal), dimension(3) :: eu real(pReal) :: zeta - + if (abs(om(3,3)) < 1.0_pReal) then zeta = 1.0_pReal/sqrt(1.0_pReal-om(3,3)**2.0_pReal) eu = [atan2(om(3,1)*zeta,-om(3,2)*zeta), & acos(om(3,3)), & atan2(om(1,3)*zeta, om(2,3)*zeta)] - else + else eu = [atan2(om(1,2),om(1,1)), 0.5_pReal*PI*(1.0_pReal-om(3,3)),0.0_pReal ] end if where(eu<0.0_pReal) eu = mod(eu+2.0_pReal*PI,[2.0_pReal*PI,PI,2.0_pReal*PI]) - + end function om2eu @@ -588,19 +599,19 @@ function om2ax(om) result(ax) real(pReal), intent(in), dimension(3,3) :: om real(pReal), dimension(4) :: ax - + real(pReal) :: t real(pReal), dimension(3) :: Wr, Wi real(pReal), dimension((64+2)*3) :: work real(pReal), dimension(3,3) :: VR, devNull, om_ integer :: ierr, i - + om_ = om - + ! first get the rotation angle t = 0.5_pReal * (math_trace33(om) - 1.0_pReal) ax(4) = acos(math_clip(t,-1.0_pReal,1.0_pReal)) - + if (dEq0(ax(4))) then ax(1:3) = [ 0.0_pReal, 0.0_pReal, 1.0_pReal ] else @@ -674,7 +685,7 @@ pure function eu2qu(eu) result(qu) real(pReal) :: cPhi, sPhi ee = 0.5_pReal*eu - + cPhi = cos(ee(2)) sPhi = sin(ee(2)) @@ -692,15 +703,15 @@ end function eu2qu !> @brief Euler angles to orientation matrix !--------------------------------------------------------------------------------------------------- pure function eu2om(eu) result(om) - + real(pReal), intent(in), dimension(3) :: eu real(pReal), dimension(3,3) :: om - - real(pReal), dimension(3) :: c, s - + + real(pReal), dimension(3) :: c, s + c = cos(eu) s = sin(eu) - + om(1,1) = c(1)*c(3)-s(1)*s(3)*c(2) om(1,2) = s(1)*c(3)+c(1)*s(3)*c(2) om(1,3) = s(3)*s(2) @@ -710,7 +721,7 @@ pure function eu2om(eu) result(om) om(3,1) = s(1)*s(2) om(3,2) = -c(1)*s(2) om(3,3) = c(2) - + where(dEq0(om)) om = 0.0_pReal end function eu2om @@ -721,19 +732,19 @@ end function eu2om !> @brief convert euler to axis angle !--------------------------------------------------------------------------------------------------- pure function eu2ax(eu) result(ax) - + real(pReal), intent(in), dimension(3) :: eu real(pReal), dimension(4) :: ax - + real(pReal) :: t, delta, tau, alpha, sigma - + t = tan(eu(2)*0.5_pReal) sigma = 0.5_pReal*(eu(1)+eu(3)) delta = 0.5_pReal*(eu(1)-eu(3)) tau = sqrt(t**2+sin(sigma)**2) - + alpha = merge(PI, 2.0_pReal*atan(tau/cos(sigma)), dEq(sigma,PI*0.5_pReal,tol=1.0e-15_pReal)) - + if (dEq0(alpha)) then ! return a default identity axis-angle pair ax = [ 0.0_pReal, 0.0_pReal, 1.0_pReal, 0.0_pReal ] else @@ -741,7 +752,7 @@ pure function eu2ax(eu) result(ax) ax(4) = alpha if (alpha < 0.0_pReal) ax = -ax ! ensure alpha is positive end if - + end function eu2ax @@ -753,7 +764,7 @@ pure function eu2ro(eu) result(ro) real(pReal), intent(in), dimension(3) :: eu real(pReal), dimension(4) :: ro - + ro = eu2ax(eu) if (ro(4) >= PI) then ro(4) = IEEE_value(ro(4),IEEE_positive_inf) @@ -762,7 +773,7 @@ pure function eu2ro(eu) result(ro) else ro(4) = tan(ro(4)*0.5_pReal) end if - + end function eu2ro @@ -799,7 +810,7 @@ end function eu2cu !> @brief convert axis angle pair to quaternion !--------------------------------------------------------------------------------------------------- pure function ax2qu(ax) result(qu) - + real(pReal), intent(in), dimension(4) :: ax real(pReal), dimension(4) :: qu @@ -825,7 +836,7 @@ pure function ax2om(ax) result(om) real(pReal), intent(in), dimension(4) :: ax real(pReal), dimension(3,3) :: om - + real(pReal) :: q, c, s, omc c = cos(ax(4)) @@ -839,11 +850,11 @@ pure function ax2om(ax) result(om) q = omc*ax(1)*ax(2) om(1,2) = q + s*ax(3) om(2,1) = q - s*ax(3) - + q = omc*ax(2)*ax(3) om(2,3) = q + s*ax(1) om(3,2) = q - s*ax(1) - + q = omc*ax(3)*ax(1) om(3,1) = q + s*ax(2) om(1,3) = q - s*ax(2) @@ -875,12 +886,12 @@ pure function ax2ro(ax) result(ro) real(pReal), intent(in), dimension(4) :: ax real(pReal), dimension(4) :: ro - + real(pReal), parameter :: thr = 1.0e-7_pReal - + if (dEq0(ax(4))) then ro = [ 0.0_pReal, 0.0_pReal, P, 0.0_pReal ] - else + else ro(1:3) = ax(1:3) ! we need to deal with the 180 degree case ro(4) = merge(IEEE_value(ro(4),IEEE_positive_inf),tan(ax(4)*0.5_pReal),abs(ax(4)-PI) < thr) @@ -897,9 +908,9 @@ pure function ax2ho(ax) result(ho) real(pReal), intent(in), dimension(4) :: ax real(pReal), dimension(3) :: ho - + real(pReal) :: f - + f = 0.75_pReal * ( ax(4) - sin(ax(4)) ) f = f**(1.0_pReal/3.0_pReal) ho = ax(1:3) * f @@ -916,7 +927,7 @@ function ax2cu(ax) result(cu) real(pReal), intent(in), dimension(4) :: ax real(pReal), dimension(3) :: cu - cu = ho2cu(ax2ho(ax)) + cu = ho2cu(ax2ho(ax)) end function ax2cu @@ -929,7 +940,7 @@ pure function ro2qu(ro) result(qu) real(pReal), intent(in), dimension(4) :: ro real(pReal), dimension(4) :: qu - + qu = ax2qu(ro2ax(ro)) end function ro2qu @@ -957,7 +968,7 @@ pure function ro2eu(ro) result(eu) real(pReal), intent(in), dimension(4) :: ro real(pReal), dimension(3) :: eu - + eu = om2eu(ro2om(ro)) end function ro2eu @@ -971,14 +982,14 @@ pure function ro2ax(ro) result(ax) real(pReal), intent(in), dimension(4) :: ro real(pReal), dimension(4) :: ax - + real(pReal) :: ta, angle - + ta = ro(4) - + if (.not. IEEE_is_finite(ta)) then ax = [ ro(1), ro(2), ro(3), PI ] - elseif (dEq0(ta)) then + elseif (dEq0(ta)) then ax = [ 0.0_pReal, 0.0_pReal, 1.0_pReal, 0.0_pReal ] else angle = 2.0_pReal*atan(ta) @@ -997,9 +1008,9 @@ pure function ro2ho(ro) result(ho) real(pReal), intent(in), dimension(4) :: ro real(pReal), dimension(3) :: ho - + real(pReal) :: f - + if (dEq0(norm2(ro(1:3)))) then ho = [ 0.0_pReal, 0.0_pReal, 0.0_pReal ] else @@ -1074,26 +1085,26 @@ pure function ho2ax(ho) result(ax) real(pReal), intent(in), dimension(3) :: ho real(pReal), dimension(4) :: ax - + integer :: i real(pReal) :: hmag_squared, s, hm real(pReal), parameter, dimension(16) :: & - tfit = [ 1.0000000000018852_pReal, -0.5000000002194847_pReal, & - -0.024999992127593126_pReal, -0.003928701544781374_pReal, & - -0.0008152701535450438_pReal, -0.0002009500426119712_pReal, & - -0.00002397986776071756_pReal, -0.00008202868926605841_pReal, & - +0.00012448715042090092_pReal, -0.0001749114214822577_pReal, & - +0.0001703481934140054_pReal, -0.00012062065004116828_pReal, & - +0.000059719705868660826_pReal, -0.00001980756723965647_pReal, & + tfit = [ 1.0000000000018852_pReal, -0.5000000002194847_pReal, & + -0.024999992127593126_pReal, -0.003928701544781374_pReal, & + -0.0008152701535450438_pReal, -0.0002009500426119712_pReal, & + -0.00002397986776071756_pReal, -0.00008202868926605841_pReal, & + +0.00012448715042090092_pReal, -0.0001749114214822577_pReal, & + +0.0001703481934140054_pReal, -0.00012062065004116828_pReal, & + +0.000059719705868660826_pReal, -0.00001980756723965647_pReal, & +0.000003953714684212874_pReal, -0.00000036555001439719544_pReal ] - + ! normalize h and store the magnitude hmag_squared = sum(ho**2.0_pReal) if (dEq0(hmag_squared)) then ax = [ 0.0_pReal, 0.0_pReal, 1.0_pReal, 0.0_pReal ] else hm = hmag_squared - + ! convert the magnitude to the rotation angle s = tfit(1) + tfit(2) * hmag_squared do i=3,16 @@ -1120,16 +1131,56 @@ pure function ho2ro(ho) result(ro) end function ho2ro -!--------------------------------------------------------------------------------------------------- +!-------------------------------------------------------------------------- !> @author Marc De Graef, Carnegie Mellon University +!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @brief convert homochoric to cubochoric -!--------------------------------------------------------------------------------------------------- +!-------------------------------------------------------------------------- pure function ho2cu(ho) result(cu) - real(pReal), intent(in), dimension(3) :: ho - real(pReal), dimension(3) :: cu + real(pReal), intent(in), dimension(3) :: ho + real(pReal), dimension(3) :: cu, xyz1, xyz3 + real(pReal), dimension(2) :: Tinv, xyz2 + real(pReal) :: rs, qxy, q2, sq2, q, tt + integer, dimension(3,2) :: p - cu = Lambert_BallToCube(ho) + rs = norm2(ho) + if (rs > R1+1.e-6_pReal) then + cu = IEEE_value(cu,IEEE_positive_inf) + return + endif + + center: if (all(dEq0(ho))) then + cu = 0.0_pReal + else center + p = GetPyramidOrder(ho) + xyz3 = ho(p(:,1)) + + ! inverse M_3 + xyz2 = xyz3(1:2) * sqrt( 2.0*rs/(rs+abs(xyz3(3))) ) + + ! inverse M_2 + qxy = sum(xyz2**2) + + special: if (dEq0(qxy)) then + Tinv = 0.0_pReal + else special + q2 = qxy + maxval(abs(xyz2))**2 + sq2 = sqrt(q2) + q = (beta/R2/R1) * sqrt(q2*qxy/(q2-maxval(abs(xyz2))*sq2)) + tt = (minval(abs(xyz2))**2+maxval(abs(xyz2))*sq2)/R2/qxy + Tinv = q * sign(1.0_pReal,xyz2) * merge([ 1.0_pReal, acos(math_clip(tt,-1.0_pReal,1.0_pReal))/PI12], & + [ acos(math_clip(tt,-1.0_pReal,1.0_pReal))/PI12, 1.0_pReal], & + abs(xyz2(2)) <= abs(xyz2(1))) + endif special + + ! inverse M_1 + xyz1 = [ Tinv(1), Tinv(2), sign(1.0_pReal,xyz3(3)) * rs / pref ] /sc + + ! reverse the coordinates back to order according to the original pyramid number + cu = xyz1(p(:,2)) + + endif center end function ho2cu @@ -1204,25 +1255,93 @@ pure function cu2ro(cu) result(ro) end function cu2ro -!--------------------------------------------------------------------------------------------------- +!-------------------------------------------------------------------------- !> @author Marc De Graef, Carnegie Mellon University -!> @brief convert cubochoric to homochoric -!--------------------------------------------------------------------------------------------------- +!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH +!> @brief map from 3D cubic grid to 3D ball +!-------------------------------------------------------------------------- pure function cu2ho(cu) result(ho) - real(pReal), intent(in), dimension(3) :: cu - real(pReal), dimension(3) :: ho + real(pReal), intent(in), dimension(3) :: cu + real(pReal), dimension(3) :: ho, LamXYZ, XYZ + real(pReal), dimension(2) :: T + real(pReal) :: c, s, q + real(pReal), parameter :: eps = 1.0e-8_pReal + integer, dimension(3,2) :: p + integer, dimension(2) :: order - ho = Lambert_CubeToBall(cu) + if (maxval(abs(cu)) > AP/2.0+eps) then + ho = IEEE_value(cu,IEEE_positive_inf) + return + end if + + ! transform to the sphere grid via the curved square, and intercept the zero point + center: if (all(dEq0(cu))) then + ho = 0.0_pReal + else center + ! get pyramide and scale by grid parameter ratio + p = GetPyramidOrder(cu) + XYZ = cu(p(:,1)) * sc + + ! intercept all the points along the z-axis + special: if (all(dEq0(XYZ(1:2)))) then + LamXYZ = [ 0.0_pReal, 0.0_pReal, pref * XYZ(3) ] + else special + order = merge( [2,1], [1,2], abs(XYZ(2)) <= abs(XYZ(1))) ! order of absolute values of XYZ + q = PI12 * XYZ(order(1))/XYZ(order(2)) ! smaller by larger + c = cos(q) + s = sin(q) + q = prek * XYZ(order(2))/ sqrt(R2-c) + T = [ (R2*c - 1.0), R2 * s] * q + + ! transform to sphere grid (inverse Lambert) + ! [note that there is no need to worry about dividing by zero, since XYZ(3) can not become zero] + c = sum(T**2) + s = Pi * c/(24.0*XYZ(3)**2) + c = sPi * c / sqrt(24.0_pReal) / XYZ(3) + q = sqrt( 1.0 - s ) + LamXYZ = [ T(order(2)) * q, T(order(1)) * q, pref * XYZ(3) - c ] + endif special + + ! reverse the coordinates back to order according to the original pyramid number + ho = LamXYZ(p(:,2)) + + endif center end function cu2ho +!-------------------------------------------------------------------------- +!> @author Marc De Graef, Carnegie Mellon University +!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH +!> @brief determine to which pyramid a point in a cubic grid belongs +!-------------------------------------------------------------------------- +pure function GetPyramidOrder(xyz) + + real(pReal),intent(in),dimension(3) :: xyz + integer, dimension(3,2) :: GetPyramidOrder + + if (((abs(xyz(1)) <= xyz(3)).and.(abs(xyz(2)) <= xyz(3))) .or. & + ((abs(xyz(1)) <= -xyz(3)).and.(abs(xyz(2)) <= -xyz(3)))) then + GetPyramidOrder = reshape([[1,2,3],[1,2,3]],[3,2]) + else if (((abs(xyz(3)) <= xyz(1)).and.(abs(xyz(2)) <= xyz(1))) .or. & + ((abs(xyz(3)) <= -xyz(1)).and.(abs(xyz(2)) <= -xyz(1)))) then + GetPyramidOrder = reshape([[2,3,1],[3,1,2]],[3,2]) + else if (((abs(xyz(1)) <= xyz(2)).and.(abs(xyz(3)) <= xyz(2))) .or. & + ((abs(xyz(1)) <= -xyz(2)).and.(abs(xyz(3)) <= -xyz(2)))) then + GetPyramidOrder = reshape([[3,1,2],[2,3,1]],[3,2]) + else + GetPyramidOrder = -1 ! should be impossible, but might simplify debugging + end if + +end function GetPyramidOrder + + !-------------------------------------------------------------------------------------------------- !> @brief check correctness of some rotations functions !-------------------------------------------------------------------------------------------------- subroutine unitTest - + type(rotation) :: R real(pReal), dimension(4) :: qu, ax, ro real(pReal), dimension(3) :: x, eu, ho, v3 @@ -1233,7 +1352,7 @@ subroutine unitTest integer :: i do i=1,10 - + msg = '' #if defined(__GFORTRAN__) && __GNUC__<9 @@ -1307,15 +1426,15 @@ subroutine unitTest #endif call R%fromMatrix(om) - + call random_number(v3) if(all(dNeq(R%rotVector(R%rotVector(v3),active=.true.),v3,1.0e-12_pReal))) & msg = trim(msg)//'rotVector,' - + call random_number(t33) if(all(dNeq(R%rotTensor2(R%rotTensor2(t33),active=.true.),t33,1.0e-12_pReal))) & msg = trim(msg)//'rotTensor2,' - + call random_number(t3333) if(all(dNeq(R%rotTensor4(R%rotTensor4(t3333),active=.true.),t3333,1.0e-12_pReal))) & msg = trim(msg)//'rotTensor4,'