diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 7b33574de..d82c7c8ac 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -213,6 +213,14 @@ Post_OrientationConversion: - master - release +Post_OrientationAverageMisorientation: + stage: postprocessing + script: + - OrientationAverageMisorientation/test.py + except: + - master + - release + ################################################################################################### grid_mech_compile_Intel: stage: compilePETSc diff --git a/PRIVATE b/PRIVATE index 76d367e29..183d7a3a3 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 76d367e297c054de78f5568074470b047427395b +Subproject commit 183d7a3a3bafa0a308c3eac858ca03c08fc03d50 diff --git a/VERSION b/VERSION index 00ca40975..4b8c1b29a 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.3-152-gd74599d3 +v2.0.3-252-g0e77be25 diff --git a/processing/post/DADF5_postResults.py b/processing/post/DADF5_postResults.py new file mode 100755 index 000000000..caf09d536 --- /dev/null +++ b/processing/post/DADF5_postResults.py @@ -0,0 +1,83 @@ +#!/usr/bin/env python3 +# -*- coding: UTF-8 no BOM -*- + +import os +import numpy as np +import argparse +import damask + +scriptName = os.path.splitext(os.path.basename(__file__))[0] +scriptID = ' '.join([scriptName,damask.version]) + +# -------------------------------------------------------------------- +# MAIN +# -------------------------------------------------------------------- +parser = argparse.ArgumentParser() + +#ToDo: We need to decide on a way of handling arguments of variable lentght +#https://stackoverflow.com/questions/15459997/passing-integer-lists-to-python + +#parser.add_argument('--version', action='version', version='%(prog)s {}'.format(scriptID)) +parser.add_argument('filenames', nargs='+', + help='DADF5 files') +parser.add_argument('-d','--dir', dest='dir',default='postProc',metavar='string', + help='name of subdirectory to hold output') + +options = parser.parse_args() + +options.labels = ['Fe','Fp','xi_sl'] + +# --- loop over input files ------------------------------------------------------------------------ + +for filename in options.filenames: + results = damask.DADF5(filename) + + if not results.structured: continue + delta = results.size/results.grid*0.5 + x, y, z = np.meshgrid(np.linspace(delta[2],results.size[2]-delta[2],results.grid[2]), + np.linspace(delta[1],results.size[1]-delta[1],results.grid[1]), + np.linspace(delta[0],results.size[0]-delta[0],results.grid[0]), + indexing = 'ij') + + coords = np.concatenate((z[:,:,:,None],y[:,:,:,None],x[:,:,:,None]),axis = 3) + + for i,inc in enumerate(results.increments): + print('Output step {}/{}'.format(i+1,len(results.increments))) + + header = '1 header\n' + + data = np.array([inc['inc'] for j in range(np.product(results.grid))]).reshape([np.product(results.grid),1]) + header+= 'inc' + + data = np.concatenate((data,np.array([j+1 for j in range(np.product(results.grid))]).reshape([np.product(results.grid),1])),1) + header+=' node' + + coords = coords.reshape([np.product(results.grid),3]) + data = np.concatenate((data,coords),1) + header+=' 1_pos 2_pos 3_pos' + + results.active['increments'] = [inc] + for label in options.labels: + for o in results.c_output_types: + results.active['c_output_types'] = [o] + for c in results.constituents: + results.active['constituents'] = [c] + x = results.get_dataset_location(label) + if len(x) == 0: + continue + label = x[0].split('/')[-1] + array = results.read_dataset(x,0) + d = np.product(np.shape(array)[1:]) + array = np.reshape(array,[np.product(results.grid),d]) + data = np.concatenate((data,array),1) + + header+= ''.join([' {}_{}'.format(j+1,label) for j in range(d)]) + + + dirname = os.path.abspath(os.path.join(os.path.dirname(filename),options.dir)) + try: + os.mkdir(dirname) + except FileExistsError: + pass + file_out = '{}_inc{:04d}.txt'.format(filename.split('.')[0],inc['inc']) + np.savetxt(os.path.join(dirname,file_out),data,header=header,comments='') diff --git a/processing/post/DADF5_vtk_cells.py b/processing/post/DADF5_vtk_cells.py index 3bbf9fd45..dc1177488 100755 --- a/processing/post/DADF5_vtk_cells.py +++ b/processing/post/DADF5_vtk_cells.py @@ -21,6 +21,8 @@ parser = argparse.ArgumentParser() #parser.add_argument('--version', action='version', version='%(prog)s {}'.format(scriptID)) parser.add_argument('filenames', nargs='+', help='DADF5 files') +parser.add_argument('-d','--dir', dest='dir',default='postProc',metavar='string', + help='name of subdirectory to hold output') options = parser.parse_args() @@ -80,12 +82,17 @@ for filename in options.filenames: if results.structured: writer = vtk.vtkXMLRectilinearGridWriter() + + dirname = os.path.abspath(os.path.join(os.path.dirname(filename),options.dir)) + try: + os.mkdir(dirname) + except FileExistsError: + pass + file_out = '{}_inc{:04d}.{}'.format(filename.split('.')[0],inc['inc'],writer.GetDefaultFileExtension()) + writer.SetCompressorTypeToZLib() writer.SetDataModeToBinary() - writer.SetFileName(os.path.join(os.path.split(filename)[0], - os.path.splitext(os.path.split(filename)[1])[0] + - '_inc{:04d}'.format(i) + # ToDo: adjust to length of increments - '.' + writer.GetDefaultFileExtension())) + writer.SetFileName(os.path.join(dirname,file_out)) if results.structured: writer.SetInputData(rGrid) diff --git a/python/damask/dadf5.py b/python/damask/dadf5.py index 4214f4922..841029af8 100644 --- a/python/damask/dadf5.py +++ b/python/damask/dadf5.py @@ -20,14 +20,14 @@ class DADF5(): with h5py.File(filename,'r') as f: - if f.attrs['DADF5-major'] != 0 or f.attrs['DADF5-minor'] != 1: + if f.attrs['DADF5-major'] != 0 or f.attrs['DADF5-minor'] != 2: raise TypeError('Unsupported DADF5 version {} '.format(f.attrs['DADF5-version'])) - self.structured = 'grid' in f['mapping'].attrs.keys() + self.structured = 'grid' in f['geometry'].attrs.keys() if self.structured: - self.grid = f['mapping'].attrs['grid'] - self.size = f['mapping'].attrs['size'] + self.grid = f['geometry'].attrs['grid'] + self.size = f['geometry'].attrs['size'] r=re.compile('inc[0-9]+') self.increments = [{'inc': int(u[3:]), diff --git a/python/damask/orientation.py b/python/damask/orientation.py index 252f39420..b19037b0c 100644 --- a/python/damask/orientation.py +++ b/python/damask/orientation.py @@ -3,211 +3,8 @@ import math import numpy as np from . import Lambert - -P = -1 - -#################################################################################################### -class Quaternion: - u""" - Quaternion with basic operations - - q is the real part, p = (x, y, z) are the imaginary parts. - Defintion of multiplication depends on variable P, P ∉ {-1,1}. - """ - - def __init__(self, - q = 0.0, - p = np.zeros(3,dtype=float)): - """Initializes to identity unless specified""" - self.q = q - self.p = np.array(p) - - - def __copy__(self): - """Copy""" - return self.__class__(q=self.q, - p=self.p.copy()) - - copy = __copy__ - - - def __iter__(self): - """Components""" - return iter(self.asList()) - - def asArray(self): - """As numpy array""" - return np.array((self.q,self.p[0],self.p[1],self.p[2])) - - def asList(self): - return [self.q]+list(self.p) - - - def __repr__(self): - """Readable string""" - return 'Quaternion: (real={q:+.6f}, imag=<{p[0]:+.6f}, {p[1]:+.6f}, {p[2]:+.6f}>)'.format(q=self.q,p=self.p) - - - def __add__(self, other): - """Addition""" - if isinstance(other, Quaternion): - return self.__class__(q=self.q + other.q, - p=self.p + other.p) - else: - return NotImplemented - - def __iadd__(self, other): - """In-place addition""" - if isinstance(other, Quaternion): - self.q += other.q - self.p += other.p - return self - else: - return NotImplemented - - def __pos__(self): - """Unary positive operator""" - return self - - - def __sub__(self, other): - """Subtraction""" - if isinstance(other, Quaternion): - return self.__class__(q=self.q - other.q, - p=self.p - other.p) - else: - return NotImplemented - - def __isub__(self, other): - """In-place subtraction""" - if isinstance(other, Quaternion): - self.q -= other.q - self.p -= other.p - return self - else: - return NotImplemented - - def __neg__(self): - """Unary positive operator""" - self.q *= -1.0 - self.p *= -1.0 - return self - - - def __mul__(self, other): - """Multiplication with quaternion or scalar""" - if isinstance(other, Quaternion): - return self.__class__(q=self.q*other.q - np.dot(self.p,other.p), - p=self.q*other.p + other.q*self.p + P * np.cross(self.p,other.p)) - elif isinstance(other, (int, float)): - return self.__class__(q=self.q*other, - p=self.p*other) - else: - return NotImplemented - - def __imul__(self, other): - """In-place multiplication with quaternion or scalar""" - if isinstance(other, Quaternion): - self.q = self.q*other.q - np.dot(self.p,other.p) - self.p = self.q*other.p + other.q*self.p + P * np.cross(self.p,other.p) - return self - elif isinstance(other, (int, float)): - self *= other - return self - else: - return NotImplemented - - - def __truediv__(self, other): - """Divsion with quaternion or scalar""" - if isinstance(other, Quaternion): - s = other.conjugate()/abs(other)**2. - return self.__class__(q=self.q * s, - p=self.p * s) - elif isinstance(other, (int, float)): - self.q /= other - self.p /= other - return self - else: - return NotImplemented - - def __itruediv__(self, other): - """In-place divsion with quaternion or scalar""" - if isinstance(other, Quaternion): - s = other.conjugate()/abs(other)**2. - self *= s - return self - elif isinstance(other, (int, float)): - self.q /= other - return self - else: - return NotImplemented - - - def __pow__(self, exponent): - """Power""" - if isinstance(exponent, (int, float)): - omega = np.acos(self.q) - return self.__class__(q= np.cos(exponent*omega), - p=self.p * np.sin(exponent*omega)/np.sin(omega)) - else: - return NotImplemented - - def __ipow__(self, exponent): - """In-place power""" - if isinstance(exponent, (int, float)): - omega = np.acos(self.q) - self.q = np.cos(exponent*omega) - self.p *= np.sin(exponent*omega)/np.sin(omega) - else: - return NotImplemented - - - def __abs__(self): - """Norm""" - return math.sqrt(self.q ** 2 + np.dot(self.p,self.p)) - - magnitude = __abs__ - - - def __eq__(self,other): - """Equal (sufficiently close) to each other""" - return np.isclose(( self-other).magnitude(),0.0) \ - or np.isclose((-self-other).magnitude(),0.0) - - def __ne__(self,other): - """Not equal (sufficiently close) to each other""" - return not self.__eq__(other) - - - def normalize(self): - d = self.magnitude() - if d > 0.0: - self.q /= d - self.p /= d - return self - - def normalized(self): - return self.copy().normalize() - - - def conjugate(self): - self.p = -self.p - return self - - def conjugated(self): - return self.copy().conjugate() - - - def homomorph(self): - if self.q < 0.0: - self.q = -self.q - self.p = -self.p - return self - - def homomorphed(self): - return self.copy().homomorph() - +from .quaternion import Quaternion +from .quaternion import P #################################################################################################### @@ -247,7 +44,13 @@ class Rotation: self.quaternion = quaternion.copy() else: self.quaternion = Quaternion(q=quaternion[0],p=quaternion[1:4]) - self.quaternion.homomorph() # ToDo: Needed? + + def __copy__(self): + """Copy""" + return self.__class__(self.quaternion) + + copy = __copy__ + def __repr__(self): """Value in selected representation""" @@ -257,6 +60,78 @@ class Rotation: 'Bunge Eulers / deg: {}'.format('\t'.join(list(map(str,self.asEulers(degrees=True)))) ), ]) + def __mul__(self, other): + """ + Multiplication + + Rotation: Details needed (active/passive), rotation of (3,3,3,3)-matrix should be considered + """ + if isinstance(other, Rotation): # rotate a rotation + return self.__class__(self.quaternion * other.quaternion).standardize() + elif isinstance(other, np.ndarray): + if other.shape == (3,): # rotate a single (3)-vector + ( x, y, z) = self.quaternion.p + (Vx,Vy,Vz) = other[0:3] + A = self.quaternion.q*self.quaternion.q - np.dot(self.quaternion.p,self.quaternion.p) + B = 2.0 * (x*Vx + y*Vy + z*Vz) + C = 2.0 * P*self.quaternion.q + + return np.array([ + A*Vx + B*x + C*(y*Vz - z*Vy), + A*Vy + B*y + C*(z*Vx - x*Vz), + A*Vz + B*z + C*(x*Vy - y*Vx), + ]) + elif other.shape == (3,3,): # rotate a single (3x3)-matrix + return np.dot(self.asMatrix(),np.dot(other,self.asMatrix().T)) + elif other.shape == (3,3,3,3,): + raise NotImplementedError + else: + return NotImplemented + elif isinstance(other, tuple): # used to rotate a meshgrid-tuple + ( x, y, z) = self.quaternion.p + (Vx,Vy,Vz) = other[0:3] + A = self.quaternion.q*self.quaternion.q - np.dot(self.quaternion.p,self.quaternion.p) + B = 2.0 * (x*Vx + y*Vy + z*Vz) + C = 2.0 * P*self.quaternion.q + + return np.array([ + A*Vx + B*x + C*(y*Vz - z*Vy), + A*Vy + B*y + C*(z*Vx - x*Vz), + A*Vz + B*z + C*(x*Vy - y*Vx), + ]) + else: + return NotImplemented + + + def inverse(self): + """In-place inverse rotation/backward rotation""" + self.quaternion.conjugate() + return self + + def inversed(self): + """Inverse rotation/backward rotation""" + return self.copy().inverse() + + + def standardize(self): + """In-place quaternion representation with positive q""" + if self.quaternion.q < 0.0: self.quaternion.homomorph() + return self + + def standardized(self): + """Quaternion representation with positive q""" + return self.copy().standardize() + + + def misorientation(self,other): + """Misorientation""" + return other*self.inversed() + + + def average(self,other): + """Calculate the average rotation""" + return Rotation.fromAverage([self,other]) + ################################################################################################ # convert to different orientation representations (numpy arrays) @@ -283,9 +158,10 @@ class Rotation: """Rotation matrix""" return qu2om(self.quaternion.asArray()) - def asRodrigues(self): + def asRodrigues(self,vector=False): """Rodrigues-Frank vector: ([n_1, n_2, n_3], tan(ω/2))""" - return qu2ro(self.quaternion.asArray()) + ro = qu2ro(self.quaternion.asArray()) + return ro[:3]*ro[3] if vector else ro def asHomochoric(self): """Homochoric vector: (h_1, h_2, h_3)""" @@ -294,6 +170,10 @@ class Rotation: def asCubochoric(self): return qu2cu(self.quaternion.asArray()) + def asM(self): + """Intermediate representation supporting quaternion averaging (see F. Landis Markley et al.)""" + return self.quaternion.asM() + ################################################################################################ # static constructors. The input data needs to follow the convention, options allow to @@ -402,63 +282,47 @@ class Rotation: return cls(ho2qu(ho)) - def __mul__(self, other): + @classmethod + def fromAverage(cls, + rotations, + weights = []): """ - Multiplication - - Rotation: Details needed (active/passive), rotation of (3,3,3,3)-matrix should be considered + Average rotation + + ref: F. Landis Markley, Yang Cheng, John Lucas Crassidis, and Yaakov Oshman. + Averaging Quaternions, + Journal of Guidance, Control, and Dynamics, Vol. 30, No. 4 (2007), pp. 1193-1197. + doi: 10.2514/1.28949 + + usage: input a list of rotations and optional weights """ - if isinstance(other, Rotation): # rotate a rotation - return self.__class__((self.quaternion * other.quaternion).asArray()) - elif isinstance(other, np.ndarray): - if other.shape == (3,): # rotate a single (3)-vector - ( x, y, z) = self.quaternion.p - (Vx,Vy,Vz) = other[0:3] - A = self.quaternion.q*self.quaternion.q - np.dot(self.quaternion.p,self.quaternion.p) - B = 2.0 * (x*Vx + y*Vy + z*Vz) - C = 2.0 * P*self.quaternion.q + if not all(isinstance(item, Rotation) for item in rotations): + raise TypeError("Only instances of Rotation can be averaged.") - return np.array([ - A*Vx + B*x + C*(y*Vz - z*Vy), - A*Vy + B*y + C*(z*Vx - x*Vz), - A*Vz + B*z + C*(x*Vy - y*Vx), - ]) - elif other.shape == (3,3,): # rotate a single (3x3)-matrix - return np.dot(self.asMatrix(),np.dot(other,self.asMatrix().T)) - elif other.shape == (3,3,3,3): - raise NotImplementedError - else: - return NotImplemented - elif isinstance(other, tuple): # used to rotate a meshgrid-tuple - ( x, y, z) = self.quaternion.p - (Vx,Vy,Vz) = other[0:3] - A = self.quaternion.q*self.quaternion.q - np.dot(self.quaternion.p,self.quaternion.p) - B = 2.0 * (x*Vx + y*Vy + z*Vz) - C = 2.0 * P*self.quaternion.q + N = len(rotations) + if weights == [] or not weights: + weights = np.ones(N,dtype='i') - return np.array([ - A*Vx + B*x + C*(y*Vz - z*Vy), - A*Vy + B*y + C*(z*Vx - x*Vz), - A*Vz + B*z + C*(x*Vy - y*Vx), - ]) - else: - return NotImplemented - - - def inverse(self): - """Inverse rotation/backward rotation""" - self.quaternion.conjugate() - return self - - def inversed(self): - """In-place inverse rotation/backward rotation""" - return self.__class__(self.quaternion.conjugated()) + for i,(r,n) in enumerate(zip(rotations,weights)): + M = r.asM() * n if i == 0 \ + else M + r.asM() * n # noqa add (multiples) of this rotation to average noqa + eig, vec = np.linalg.eig(M/N) + + return cls.fromQuaternion(np.real(vec.T[eig.argmax()]),acceptHomomorph = True) + + + @classmethod + def fromRandom(cls): + r = np.random.random(3) + A = np.sqrt(r[2]) + B = np.sqrt(1.0-r[2]) + w = np.cos(2.0*np.pi*r[0])*A + x = np.sin(2.0*np.pi*r[1])*B + y = np.cos(2.0*np.pi*r[1])*B + z = np.sin(2.0*np.pi*r[0])*A + return cls.fromQuaternion([w,x,y,z],acceptHomomorph=True) - def misorientation(self,other): - """Misorientation""" - return self.__class__(other.quaternion*self.quaternion.conjugated()) - # ****************************************************************************************** class Symmetry: @@ -503,60 +367,60 @@ class Symmetry: otherOrder = Symmetry.lattices.index(other.lattice) return (myOrder > otherOrder) - (myOrder < otherOrder) - def symmetryOperations(self): - """List of symmetry operations as quaternions.""" + def symmetryOperations(self,members=[]): + """List (or single element) of symmetry operations as rotations.""" if self.lattice == 'cubic': symQuats = [ - [ 1.0, 0.0, 0.0, 0.0 ], - [ 0.0, 1.0, 0.0, 0.0 ], - [ 0.0, 0.0, 1.0, 0.0 ], - [ 0.0, 0.0, 0.0, 1.0 ], - [ 0.0, 0.0, 0.5*math.sqrt(2), 0.5*math.sqrt(2) ], - [ 0.0, 0.0, 0.5*math.sqrt(2),-0.5*math.sqrt(2) ], - [ 0.0, 0.5*math.sqrt(2), 0.0, 0.5*math.sqrt(2) ], - [ 0.0, 0.5*math.sqrt(2), 0.0, -0.5*math.sqrt(2) ], - [ 0.0, 0.5*math.sqrt(2),-0.5*math.sqrt(2), 0.0 ], - [ 0.0, -0.5*math.sqrt(2),-0.5*math.sqrt(2), 0.0 ], - [ 0.5, 0.5, 0.5, 0.5 ], - [-0.5, 0.5, 0.5, 0.5 ], - [-0.5, 0.5, 0.5, -0.5 ], - [-0.5, 0.5, -0.5, 0.5 ], - [-0.5, -0.5, 0.5, 0.5 ], - [-0.5, -0.5, 0.5, -0.5 ], - [-0.5, -0.5, -0.5, 0.5 ], - [-0.5, 0.5, -0.5, -0.5 ], - [-0.5*math.sqrt(2), 0.0, 0.0, 0.5*math.sqrt(2) ], - [ 0.5*math.sqrt(2), 0.0, 0.0, 0.5*math.sqrt(2) ], - [-0.5*math.sqrt(2), 0.0, 0.5*math.sqrt(2), 0.0 ], - [-0.5*math.sqrt(2), 0.0, -0.5*math.sqrt(2), 0.0 ], - [-0.5*math.sqrt(2), 0.5*math.sqrt(2), 0.0, 0.0 ], - [-0.5*math.sqrt(2),-0.5*math.sqrt(2), 0.0, 0.0 ], + [ 1.0, 0.0, 0.0, 0.0 ], + [ 0.0, 1.0, 0.0, 0.0 ], + [ 0.0, 0.0, 1.0, 0.0 ], + [ 0.0, 0.0, 0.0, 1.0 ], + [ 0.0, 0.0, 0.5*np.sqrt(2), 0.5*np.sqrt(2) ], + [ 0.0, 0.0, 0.5*np.sqrt(2),-0.5*np.sqrt(2) ], + [ 0.0, 0.5*np.sqrt(2), 0.0, 0.5*np.sqrt(2) ], + [ 0.0, 0.5*np.sqrt(2), 0.0, -0.5*np.sqrt(2) ], + [ 0.0, 0.5*np.sqrt(2),-0.5*np.sqrt(2), 0.0 ], + [ 0.0, -0.5*np.sqrt(2),-0.5*np.sqrt(2), 0.0 ], + [ 0.5, 0.5, 0.5, 0.5 ], + [-0.5, 0.5, 0.5, 0.5 ], + [-0.5, 0.5, 0.5, -0.5 ], + [-0.5, 0.5, -0.5, 0.5 ], + [-0.5, -0.5, 0.5, 0.5 ], + [-0.5, -0.5, 0.5, -0.5 ], + [-0.5, -0.5, -0.5, 0.5 ], + [-0.5, 0.5, -0.5, -0.5 ], + [-0.5*np.sqrt(2), 0.0, 0.0, 0.5*np.sqrt(2) ], + [ 0.5*np.sqrt(2), 0.0, 0.0, 0.5*np.sqrt(2) ], + [-0.5*np.sqrt(2), 0.0, 0.5*np.sqrt(2), 0.0 ], + [-0.5*np.sqrt(2), 0.0, -0.5*np.sqrt(2), 0.0 ], + [-0.5*np.sqrt(2), 0.5*np.sqrt(2), 0.0, 0.0 ], + [-0.5*np.sqrt(2),-0.5*np.sqrt(2), 0.0, 0.0 ], ] elif self.lattice == 'hexagonal': symQuats = [ - [ 1.0,0.0,0.0,0.0 ], - [-0.5*math.sqrt(3), 0.0, 0.0,-0.5 ], - [ 0.5, 0.0, 0.0, 0.5*math.sqrt(3) ], - [ 0.0,0.0,0.0,1.0 ], - [-0.5, 0.0, 0.0, 0.5*math.sqrt(3) ], - [-0.5*math.sqrt(3), 0.0, 0.0, 0.5 ], - [ 0.0,1.0,0.0,0.0 ], - [ 0.0,-0.5*math.sqrt(3), 0.5, 0.0 ], - [ 0.0, 0.5,-0.5*math.sqrt(3), 0.0 ], - [ 0.0,0.0,1.0,0.0 ], - [ 0.0,-0.5,-0.5*math.sqrt(3), 0.0 ], - [ 0.0, 0.5*math.sqrt(3), 0.5, 0.0 ], + [ 1.0, 0.0, 0.0, 0.0 ], + [-0.5*np.sqrt(3), 0.0, 0.0, -0.5 ], + [ 0.5, 0.0, 0.0, 0.5*np.sqrt(3) ], + [ 0.0, 0.0, 0.0, 1.0 ], + [-0.5, 0.0, 0.0, 0.5*np.sqrt(3) ], + [-0.5*np.sqrt(3), 0.0, 0.0, 0.5 ], + [ 0.0, 1.0, 0.0, 0.0 ], + [ 0.0, -0.5*np.sqrt(3), 0.5, 0.0 ], + [ 0.0, 0.5, -0.5*np.sqrt(3), 0.0 ], + [ 0.0, 0.0, 1.0, 0.0 ], + [ 0.0, -0.5, -0.5*np.sqrt(3), 0.0 ], + [ 0.0, 0.5*np.sqrt(3), 0.5, 0.0 ], ] elif self.lattice == 'tetragonal': symQuats = [ - [ 1.0,0.0,0.0,0.0 ], - [ 0.0,1.0,0.0,0.0 ], - [ 0.0,0.0,1.0,0.0 ], - [ 0.0,0.0,0.0,1.0 ], - [ 0.0, 0.5*math.sqrt(2), 0.5*math.sqrt(2), 0.0 ], - [ 0.0,-0.5*math.sqrt(2), 0.5*math.sqrt(2), 0.0 ], - [ 0.5*math.sqrt(2), 0.0, 0.0, 0.5*math.sqrt(2) ], - [-0.5*math.sqrt(2), 0.0, 0.0, 0.5*math.sqrt(2) ], + [ 1.0, 0.0, 0.0, 0.0 ], + [ 0.0, 1.0, 0.0, 0.0 ], + [ 0.0, 0.0, 1.0, 0.0 ], + [ 0.0, 0.0, 0.0, 1.0 ], + [ 0.0, 0.5*np.sqrt(2), 0.5*np.sqrt(2), 0.0 ], + [ 0.0, -0.5*np.sqrt(2), 0.5*np.sqrt(2), 0.0 ], + [ 0.5*np.sqrt(2), 0.0, 0.0, 0.5*np.sqrt(2) ], + [-0.5*np.sqrt(2), 0.0, 0.0, 0.5*np.sqrt(2) ], ] elif self.lattice == 'orthorhombic': symQuats = [ @@ -570,16 +434,28 @@ class Symmetry: [ 1.0,0.0,0.0,0.0 ], ] - return [Rotation(q) for q in symQuats] + symOps = list(map(Rotation, + np.array(symQuats)[np.atleast_1d(members) if members != [] else range(len(symQuats))])) + try: + iter(members) # asking for (even empty) list of members? + except TypeError: + return symOps[0] # no, return rotation object + else: + return symOps # yes, return list of rotations - def inFZ(self,R): + def inFZ(self,rodrigues): """ Check whether given Rodrigues vector falls into fundamental zone of own symmetry. Fundamental zone in Rodrigues space is point symmetric around origin. """ - Rabs = abs(R[0:3]*R[3]) + if (len(rodrigues) != 3): + raise ValueError('Input is not a Rodriques-Frank vector.\n') + + if np.any(rodrigues == np.inf): return False + + Rabs = abs(rodrigues) if self.lattice == 'cubic': return math.sqrt(2.0)-1.0 >= Rabs[0] \ @@ -609,14 +485,10 @@ class Symmetry: Representation of Orientation and Disorientation Data for Cubic, Hexagonal, Tetragonal and Orthorhombic Crystals Acta Cryst. (1991). A47, 780-789 """ - if isinstance(rodrigues, Quaternion): - R = rodrigues.asRodrigues() # translate accidentially passed quaternion - else: - R = rodrigues + if (len(rodrigues) != 3): + raise ValueError('Input is not a Rodriques-Frank vector.\n') + R = rodrigues - if R.shape[0]==4: # transition old (length not stored separately) to new - R = (R[0:3]*R[3]) - epsilon = 0.0 if self.lattice == 'cubic': return R[0] >= R[1]+epsilon and R[1] >= R[2]+epsilon and R[2] >= epsilon @@ -864,7 +736,7 @@ class Lattice: # Greninger--Troiano' orientation relationship for fcc <-> bcc transformation # from Y. He et al. Journal of Applied Crystallography 39 (2006) 72-81 - GTdash = {'mapping':{'fcc':0,'bcc':1}, + GTprime = {'mapping':{'fcc':0,'bcc':1}, 'planes': np.array([ [[ 7, 17, 17],[ 12, 5, 17]], [[ 17, 7, 17],[ 17, 12, 5]], @@ -977,7 +849,7 @@ class Lattice: [[ 1, 1, 0],[ 1, 1, -1]]],dtype='float')} # Bain orientation relationship for fcc <-> bcc transformation - # from Y. He et al. Journal of Applied Crystallography 39 (2006) 72-81 + # from Y. He et al./Journal of Applied Crystallography (2006). 39, 72-81 Bain = {'mapping':{'fcc':0,'bcc':1}, 'planes': np.array([ [[ 1, 0, 0],[ 1, 0, 0]], @@ -990,7 +862,7 @@ class Lattice: def relationOperations(self,model): - models={'KS':self.KS, 'GT':self.GT, "GT'":self.GTdash, + models={'KS':self.KS, 'GT':self.GT, "GT'":self.GTprime, 'NW':self.NW, 'Pitsch': self.Pitsch, 'Bain':self.Bain} try: relationship = models[model] @@ -1050,7 +922,8 @@ class Orientation: def disorientation(self, other, - SST = True): + SST = True, + symmetries = False): """ Disorientation between myself and given other orientation. @@ -1058,33 +931,42 @@ class Orientation: (Currently requires same symmetry for both orientations. Look into A. Heinz and P. Neumann 1991 for cases with differing sym.) """ - #if self.lattice.symmetry != other.lattice.symmetry: - # raise NotImplementedError('disorientation between different symmetry classes not supported yet.') + if self.lattice.symmetry != other.lattice.symmetry: + raise NotImplementedError('disorientation between different symmetry classes not supported yet.') - mis = other.rotation*self.rotation.inversed() - mySymEqs = self.equivalentOrientations() if SST else self.equivalentOrientations()[:1] # take all or only first sym operation + mySymEqs = self.equivalentOrientations() if SST else self.equivalentOrientations([0]) # take all or only first sym operation otherSymEqs = other.equivalentOrientations() for i,sA in enumerate(mySymEqs): + aInv = sA.rotation.inversed() for j,sB in enumerate(otherSymEqs): - r = sB.rotation*mis*sA.rotation.inversed() + b = sB.rotation + r = b*aInv for k in range(2): - r.inversed() - breaker = self.lattice.symmetry.inFZ(r.asRodrigues()) \ - and (not SST or other.lattice.symmetry.inDisorientationSST(r.asRodrigues())) + r.inverse() + breaker = self.lattice.symmetry.inFZ(r.asRodrigues(vector=True)) \ + and (not SST or other.lattice.symmetry.inDisorientationSST(r.asRodrigues(vector=True))) if breaker: break if breaker: break if breaker: break - return r + return (Orientation(r,self.lattice), i,j, k == 1) if symmetries else r # disorientation ... + # ... own sym, other sym, + # self-->other: True, self<--other: False + def inFZ(self): - return self.lattice.symmetry.inFZ(self.rotation.asRodrigues()) + return self.lattice.symmetry.inFZ(self.rotation.asRodrigues(vector=True)) - def equivalentOrientations(self): + def equivalentOrientations(self,members=[]): """List of orientations which are symmetrically equivalent""" - return [self.__class__(q*self.rotation,self.lattice) \ - for q in self.lattice.symmetry.symmetryOperations()] + try: + iter(members) # asking for (even empty) list of members? + except TypeError: + return self.__class__(self.lattice.symmetry.symmetryOperations(members)*self.rotation,self.lattice) # no, return rotation object + else: + return [self.__class__(q*self.rotation,self.lattice) \ + for q in self.lattice.symmetry.symmetryOperations(members)] # yes, return list of rotations def relatedOrientations(self,model): """List of orientations related by the given orientation relationship""" @@ -1094,7 +976,7 @@ class Orientation: def reduced(self): """Transform orientation to fall into fundamental zone according to symmetry""" for me in self.equivalentOrientations(): - if self.lattice.symmetry.inFZ(me.rotation.asRodrigues()): break + if self.lattice.symmetry.inFZ(me.rotation.asRodrigues(vector=True)): break return self.__class__(me.rotation,self.lattice) @@ -1125,37 +1007,28 @@ class Orientation: return color - # @classmethod - # def average(cls, - # orientations, - # multiplicity = []): - # """ - # Average orientation + @classmethod + def fromAverage(cls, + orientations, + weights = []): + """Create orientation from average of list of orientations""" + if not all(isinstance(item, Orientation) for item in orientations): + raise TypeError("Only instances of Orientation can be averaged.") - # ref: F. Landis Markley, Yang Cheng, John Lucas Crassidis, and Yaakov Oshman. - # Averaging Quaternions, - # Journal of Guidance, Control, and Dynamics, Vol. 30, No. 4 (2007), pp. 1193-1197. - # doi: 10.2514/1.28949 - # usage: - # a = Orientation(Eulers=np.radians([10, 10, 0]), symmetry='hexagonal') - # b = Orientation(Eulers=np.radians([20, 0, 0]), symmetry='hexagonal') - # avg = Orientation.average([a,b]) - # """ - # if not all(isinstance(item, Orientation) for item in orientations): - # raise TypeError("Only instances of Orientation can be averaged.") + closest = [] + ref = orientations[0] + for o in orientations: + closest.append(o.equivalentOrientations( + ref.disorientation(o, + SST = False, # select (o[ther]'s) sym orientation + symmetries = True)[2]).rotation) # with lowest misorientation - # N = len(orientations) - # if multiplicity == [] or not multiplicity: - # multiplicity = np.ones(N,dtype='i') + return Orientation(Rotation.fromAverage(closest,weights),ref.lattice) - # reference = orientations[0] # take first as reference - # for i,(o,n) in enumerate(zip(orientations,multiplicity)): - # closest = o.equivalentOrientations(reference.disorientation(o,SST = False)[2])[0] # select sym orientation with lowest misorientation - # M = closest.quaternion.asM() * n if i == 0 else M + closest.quaternion.asM() * n # noqa add (multiples) of this orientation to average noqa - # eig, vec = np.linalg.eig(M/N) - # return Orientation(quaternion = Quaternion(quat = np.real(vec.T[eig.argmax()])), - # symmetry = reference.symmetry.lattice) + def average(self,other): + """Calculate the average rotation""" + return Orientation.fromAverage([self,other]) #################################################################################################### @@ -1195,7 +1068,167 @@ def isone(a): def iszero(a): return np.isclose(a,0.0,atol=1.0e-12,rtol=0.0) +#---------- quaternion ---------- + +def qu2om(qu): + """Quaternion to orientation matrix""" + qq = qu[0]**2-(qu[1]**2 + qu[2]**2 + qu[3]**2) + om = np.diag(qq + 2.0*np.array([qu[1],qu[2],qu[3]])**2) + + om[1,0] = 2.0*(qu[2]*qu[1]+qu[0]*qu[3]) + om[0,1] = 2.0*(qu[1]*qu[2]-qu[0]*qu[3]) + om[2,1] = 2.0*(qu[3]*qu[2]+qu[0]*qu[1]) + om[1,2] = 2.0*(qu[2]*qu[3]-qu[0]*qu[1]) + om[0,2] = 2.0*(qu[1]*qu[3]+qu[0]*qu[2]) + om[2,0] = 2.0*(qu[3]*qu[1]-qu[0]*qu[2]) + return om if P > 0.0 else om.T + +def qu2eu(qu): + """Quaternion to Euler angles""" + q03 = qu[0]**2+qu[3]**2 + q12 = qu[1]**2+qu[2]**2 + chi = np.sqrt(q03*q12) + + if iszero(chi): + eu = np.array([np.arctan2(-P*2.0*qu[0]*qu[3],qu[0]**2-qu[3]**2), 0.0, 0.0]) if iszero(q12) else \ + np.array([np.arctan2(2.0*qu[1]*qu[2],qu[1]**2-qu[2]**2), np.pi, 0.0]) + else: + eu = np.array([np.arctan2((-P*qu[0]*qu[2]+qu[1]*qu[3])*chi, (-P*qu[0]*qu[1]-qu[2]*qu[3])*chi ), + np.arctan2( 2.0*chi, q03-q12 ), + np.arctan2(( P*qu[0]*qu[2]+qu[1]*qu[3])*chi, (-P*qu[0]*qu[1]+qu[2]*qu[3])*chi )]) + + # reduce Euler angles to definition range, i.e a lower limit of 0.0 + eu = np.where(eu<0, (eu+2.0*np.pi)%np.array([2.0*np.pi,np.pi,2.0*np.pi]),eu) + return eu + + +def qu2ax(qu): + """ + Quaternion to axis angle + + Modified version of the original formulation, should be numerically more stable + """ + if iszero(qu[1]**2+qu[2]**2+qu[3]**2): # set axis to [001] if the angle is 0/360 + ax = [ 0.0, 0.0, 1.0, 0.0 ] + elif not iszero(qu[0]): + s = np.sign(qu[0])/np.sqrt(qu[1]**2+qu[2]**2+qu[3]**2) + omega = 2.0 * np.arccos(np.clip(qu[0],-1.0,1.0)) + ax = [ qu[1]*s, qu[2]*s, qu[3]*s, omega ] + else: + ax = [ qu[1], qu[2], qu[3], np.pi] + + return np.array(ax) + + +def qu2ro(qu): + """Quaternion to Rodrigues vector""" + if iszero(qu[0]): + ro = [qu[1], qu[2], qu[3], np.inf] + else: + s = np.linalg.norm([qu[1],qu[2],qu[3]]) + ro = [0.0,0.0,P,0.0] if iszero(s) else \ + [ qu[1]/s, qu[2]/s, qu[3]/s, np.tan(np.arccos(np.clip(qu[0],-1.0,1.0)))] # avoid numerical difficulties + + return np.array(ro) + + +def qu2ho(qu): + """Quaternion to homochoric""" + omega = 2.0 * np.arccos(np.clip(qu[0],-1.0,1.0)) # avoid numerical difficulties + + if iszero(omega): + ho = np.array([ 0.0, 0.0, 0.0 ]) + else: + ho = np.array([qu[1], qu[2], qu[3]]) + f = 0.75 * ( omega - np.sin(omega) ) + ho = ho/np.linalg.norm(ho) * f**(1./3.) + + return ho + + +def qu2cu(qu): + """Quaternion to cubochoric""" + return ho2cu(qu2ho(qu)) + + +#---------- orientation matrix ---------- + +def om2qu(om): + """ + Orientation matrix to quaternion + + The original formulation (direct conversion) had (numerical?) issues + """ + return eu2qu(om2eu(om)) + + +def om2eu(om): + """Orientation matrix to Euler angles""" + if abs(om[2,2]) < 1.0: + zeta = 1.0/np.sqrt(1.0-om[2,2]**2) + eu = np.array([np.arctan2(om[2,0]*zeta,-om[2,1]*zeta), + np.arccos(om[2,2]), + np.arctan2(om[0,2]*zeta, om[1,2]*zeta)]) + else: + eu = np.array([np.arctan2( om[0,1],om[0,0]), np.pi*0.5*(1-om[2,2]),0.0]) # following the paper, not the reference implementation + + # reduce Euler angles to definition range, i.e a lower limit of 0.0 + eu = np.where(eu<0, (eu+2.0*np.pi)%np.array([2.0*np.pi,np.pi,2.0*np.pi]),eu) + return eu + + +def om2ax(om): + """Orientation matrix to axis angle""" + ax=np.empty(4) + + # first get the rotation angle + t = 0.5*(om.trace() -1.0) + ax[3] = np.arccos(np.clip(t,-1.0,1.0)) + + if iszero(ax[3]): + ax = [ 0.0, 0.0, 1.0, 0.0] + else: + w,vr = np.linalg.eig(om) + # next, find the eigenvalue (1,0j) + i = np.where(np.isclose(w,1.0+0.0j))[0][0] + ax[0:3] = np.real(vr[0:3,i]) + diagDelta = np.array([om[1,2]-om[2,1],om[2,0]-om[0,2],om[0,1]-om[1,0]]) + ax[0:3] = np.where(iszero(diagDelta), ax[0:3],np.abs(ax[0:3])*np.sign(-P*diagDelta)) + + return np.array(ax) + + +def om2ro(om): + """Orientation matrix to Rodriques vector""" + return eu2ro(om2eu(om)) + + +def om2ho(om): + """Orientation matrix to homochoric""" + return ax2ho(om2ax(om)) + + +def om2cu(om): + """Orientation matrix to cubochoric""" + return ho2cu(om2ho(om)) + + +#---------- Euler angles ---------- + +def eu2qu(eu): + """Euler angles to quaternion""" + ee = 0.5*eu + cPhi = np.cos(ee[1]) + sPhi = np.sin(ee[1]) + qu = np.array([ cPhi*np.cos(ee[0]+ee[2]), + -P*sPhi*np.cos(ee[0]-ee[2]), + -P*sPhi*np.sin(ee[0]-ee[2]), + -P*cPhi*np.sin(ee[0]+ee[2]) ]) + if qu[0] < 0.0: qu*=-1 + return qu + + def eu2om(eu): """Euler angles to orientation matrix""" c = np.cos(eu) @@ -1241,32 +1274,28 @@ def eu2ro(eu): return ro -def eu2qu(eu): - """Euler angles to quaternion""" - ee = 0.5*eu - cPhi = np.cos(ee[1]) - sPhi = np.sin(ee[1]) - qu = np.array([ cPhi*np.cos(ee[0]+ee[2]), - -P*sPhi*np.cos(ee[0]-ee[2]), - -P*sPhi*np.sin(ee[0]-ee[2]), - -P*cPhi*np.sin(ee[0]+ee[2]) ]) - #if qu[0] < 0.0: qu.homomorph() !ToDo: Check with original - return qu +def eu2ho(eu): + """Euler angles to homochoric""" + return ax2ho(eu2ax(eu)) -def om2eu(om): - """Euler angles to orientation matrix""" - if isone(om[2,2]**2): - eu = np.array([np.arctan2( om[0,1],om[0,0]), np.pi*0.5*(1-om[2,2]),0.0]) # following the paper, not the reference implementation +def eu2cu(eu): + """Euler angles to cubochoric""" + return ho2cu(eu2ho(eu)) + + +#---------- axis angle ---------- + +def ax2qu(ax): + """Axis angle to quaternion""" + if iszero(ax[3]): + qu = np.array([ 1.0, 0.0, 0.0, 0.0 ]) else: - zeta = 1.0/np.sqrt(1.0-om[2,2]**2) - eu = np.array([np.arctan2(om[2,0]*zeta,-om[2,1]*zeta), - np.arccos(om[2,2]), - np.arctan2(om[0,2]*zeta, om[1,2]*zeta)]) + c = np.cos(ax[3]*0.5) + s = np.sin(ax[3]*0.5) + qu = np.array([ c, ax[0]*s, ax[1]*s, ax[2]*s ]) - # reduce Euler angles to definition range, i.e a lower limit of 0.0 - eu = np.where(eu<0, (eu+2.0*np.pi)%np.array([2.0*np.pi,np.pi,2.0*np.pi]),eu) - return eu + return qu def ax2om(ax): @@ -1283,24 +1312,23 @@ def ax2om(ax): return om if P < 0.0 else om.T + +def ax2eu(ax): + """Orientation matrix to Euler angles""" + return om2eu(ax2om(ax)) -def qu2eu(qu): - """Quaternion to Euler angles""" - q03 = qu[0]**2+qu[3]**2 - q12 = qu[1]**2+qu[2]**2 - chi = np.sqrt(q03*q12) - - if iszero(chi): - eu = np.array([np.arctan2(-P*2.0*qu[0]*qu[3],qu[0]**2-qu[3]**2), 0.0, 0.0]) if iszero(q12) else \ - np.array([np.arctan2(2.0*qu[1]*qu[2],qu[1]**2-qu[2]**2), np.pi, 0.0]) + +def ax2ro(ax): + """Axis angle to Rodrigues vector""" + if iszero(ax[3]): + ro = [ 0.0, 0.0, P, 0.0 ] else: - eu = np.array([np.arctan2((-P*qu[0]*qu[2]+qu[1]*qu[3])*chi, (-P*qu[0]*qu[1]-qu[2]*qu[3])*chi ), - np.arctan2( 2.0*chi, q03-q12 ), - np.arctan2(( P*qu[0]*qu[2]+qu[1]*qu[3])*chi, (-P*qu[0]*qu[1]+qu[2]*qu[3])*chi )]) - - # reduce Euler angles to definition range, i.e a lower limit of 0.0 - eu = np.where(eu<0, (eu+2.0*np.pi)%np.array([2.0*np.pi,np.pi,2.0*np.pi]),eu) - return eu + ro = [ax[0], ax[1], ax[2]] + # 180 degree case + ro += [np.inf] if np.isclose(ax[3],np.pi,atol=1.0e-15,rtol=0.0) else \ + [np.tan(ax[3]*0.5)] + + return np.array(ro) def ax2ho(ax): @@ -1310,6 +1338,77 @@ def ax2ho(ax): return ho +def ax2cu(ax): + """Axis angle to cubochoric""" + return ho2cu(ax2ho(ax)) + + +#---------- Rodrigues--Frank ---------- + +def ro2qu(ro): + """Rodrigues vector to quaternion""" + return ax2qu(ro2ax(ro)) + + +def ro2om(ro): + """Rodgrigues vector to orientation matrix""" + return ax2om(ro2ax(ro)) + + +def ro2eu(ro): + """Rodrigues vector to orientation matrix""" + return om2eu(ro2om(ro)) + + +def ro2ax(ro): + """Rodrigues vector to axis angle""" + ta = ro[3] + + if iszero(ta): + ax = [ 0.0, 0.0, 1.0, 0.0 ] + elif not np.isfinite(ta): + ax = [ ro[0], ro[1], ro[2], np.pi ] + else: + angle = 2.0*np.arctan(ta) + ta = 1.0/np.linalg.norm(ro[0:3]) + ax = [ ro[0]/ta, ro[1]/ta, ro[2]/ta, angle ] + + return np.array(ax) + + +def ro2ho(ro): + """Rodrigues vector to homochoric""" + if iszero(np.sum(ro[0:3]**2.0)): + ho = [ 0.0, 0.0, 0.0 ] + else: + f = 2.0*np.arctan(ro[3]) -np.sin(2.0*np.arctan(ro[3])) if np.isfinite(ro[3]) else np.pi + ho = ro[0:3] * (0.75*f)**(1.0/3.0) + + return np.array(ho) + + +def ro2cu(ro): + """Rodrigues vector to cubochoric""" + return ho2cu(ro2ho(ro)) + + +#---------- homochoric ---------- + +def ho2qu(ho): + """Homochoric to quaternion""" + return ax2qu(ho2ax(ho)) + + +def ho2om(ho): + """Homochoric to orientation matrix""" + return ax2om(ho2ax(ho)) + + +def ho2eu(ho): + """Homochoric to Euler angles""" + return ax2eu(ho2ax(ho)) + + def ho2ax(ho): """Homochoric to axis angle""" tfit = np.array([+1.0000000000018852, -0.5000000002194847, @@ -1336,135 +1435,9 @@ def ho2ax(ho): return ax -def om2ax(om): - """Orientation matrix to axis angle""" - ax=np.empty(4) - - # first get the rotation angle - t = 0.5*(om.trace() -1.0) - ax[3] = np.arccos(np.clip(t,-1.0,1.0)) - - if iszero(ax[3]): - ax = [ 0.0, 0.0, 1.0, 0.0] - else: - w,vr = np.linalg.eig(om) - # next, find the eigenvalue (1,0j) - i = np.where(np.isclose(w,1.0+0.0j))[0][0] - ax[0:3] = np.real(vr[0:3,i]) - diagDelta = np.array([om[1,2]-om[2,1],om[2,0]-om[0,2],om[0,1]-om[1,0]]) - ax[0:3] = np.where(iszero(diagDelta), ax[0:3],np.abs(ax[0:3])*np.sign(-P*diagDelta)) - - return np.array(ax) - - -def ro2ax(ro): - """Rodrigues vector to axis angle""" - ta = ro[3] - - if iszero(ta): - ax = [ 0.0, 0.0, 1.0, 0.0 ] - elif not np.isfinite(ta): - ax = [ ro[0], ro[1], ro[2], np.pi ] - else: - angle = 2.0*np.arctan(ta) - ta = 1.0/np.linalg.norm(ro[0:3]) - ax = [ ro[0]/ta, ro[1]/ta, ro[2]/ta, angle ] - - return np.array(ax) - - -def ax2ro(ax): - """Axis angle to Rodrigues vector""" - if iszero(ax[3]): - ro = [ 0.0, 0.0, P, 0.0 ] - else: - ro = [ax[0], ax[1], ax[2]] - # 180 degree case - ro += [np.inf] if np.isclose(ax[3],np.pi,atol=1.0e-15,rtol=0.0) else \ - [np.tan(ax[3]*0.5)] - - return np.array(ro) - - -def ax2qu(ax): - """Axis angle to quaternion""" - if iszero(ax[3]): - qu = np.array([ 1.0, 0.0, 0.0, 0.0 ]) - else: - c = np.cos(ax[3]*0.5) - s = np.sin(ax[3]*0.5) - qu = np.array([ c, ax[0]*s, ax[1]*s, ax[2]*s ]) - - return qu - - -def ro2ho(ro): - """Rodrigues vector to homochoric""" - if iszero(np.sum(ro[0:3]**2.0)): - ho = [ 0.0, 0.0, 0.0 ] - else: - f = 2.0*np.arctan(ro[3]) -np.sin(2.0*np.arctan(ro[3])) if np.isfinite(ro[3]) else np.pi - ho = ro[0:3] * (0.75*f)**(1.0/3.0) - - return np.array(ho) - - -def qu2om(qu): - """Quaternion to orientation matrix""" - qq = qu[0]**2-(qu[1]**2 + qu[2]**2 + qu[3]**2) - om = np.diag(qq + 2.0*np.array([qu[1],qu[2],qu[3]])**2) - - om[1,0] = 2.0*(qu[2]*qu[1]+qu[0]*qu[3]) - om[0,1] = 2.0*(qu[1]*qu[2]-qu[0]*qu[3]) - om[2,1] = 2.0*(qu[3]*qu[2]+qu[0]*qu[1]) - om[1,2] = 2.0*(qu[2]*qu[3]-qu[0]*qu[1]) - om[0,2] = 2.0*(qu[1]*qu[3]+qu[0]*qu[2]) - om[2,0] = 2.0*(qu[3]*qu[1]-qu[0]*qu[2]) - return om if P > 0.0 else om.T - - -def qu2ax(qu): - """ - Quaternion to axis angle - - Modified version of the original formulation, should be numerically more stable - """ - if isone(abs(qu[0])): # set axis to [001] if the angle is 0/360 - ax = [ 0.0, 0.0, 1.0, 0.0 ] - elif not iszero(qu[0]): - omega = 2.0 * np.arccos(qu[0]) - s = np.sign(qu[0])/np.sqrt(qu[1]**2+qu[2]**2+qu[3]**2) - ax = [ qu[1]*s, qu[2]*s, qu[3]*s, omega ] - else: - ax = [ qu[1], qu[2], qu[3], np.pi] - - return np.array(ax) - - -def qu2ro(qu): - """Quaternion to Rodrigues vector""" - if iszero(qu[0]): - ro = [qu[1], qu[2], qu[3], np.inf] - else: - s = np.linalg.norm([qu[1],qu[2],qu[3]]) - ro = [0.0,0.0,P,0.0] if iszero(s) else \ - [ qu[1]/s, qu[2]/s, qu[3]/s, np.tan(np.arccos(np.clip(qu[0],-1.0,1.0)))] # avoid numerical difficulties - - return np.array(ro) - - -def qu2ho(qu): - """Quaternion to homochoric""" - omega = 2.0 * np.arccos(np.clip(qu[0],-1.0,1.0)) # avoid numerical difficulties - - if iszero(omega): - ho = np.array([ 0.0, 0.0, 0.0 ]) - else: - ho = np.array([qu[1], qu[2], qu[3]]) - f = 0.75 * ( omega - np.sin(omega) ) - ho = ho/np.linalg.norm(ho) * f**(1./3.) - - return ho +def ho2ro(ho): + """Axis angle to Rodriques vector""" + return ax2ro(ho2ax(ho)) def ho2cu(ho): @@ -1472,103 +1445,11 @@ def ho2cu(ho): return Lambert.BallToCube(ho) -def cu2ho(cu): - """Cubochoric to homochoric""" - return Lambert.CubeToBall(cu) +#---------- cubochoric ---------- - -def ro2eu(ro): - """Rodrigues vector to orientation matrix""" - return om2eu(ro2om(ro)) - - -def eu2ho(eu): - """Euler angles to homochoric""" - return ax2ho(eu2ax(eu)) - - -def om2ro(om): - """Orientation matrix to Rodriques vector""" - return eu2ro(om2eu(om)) - - -def om2ho(om): - """Orientation matrix to homochoric""" - return ax2ho(om2ax(om)) - - -def ax2eu(ax): - """Orientation matrix to Euler angles""" - return om2eu(ax2om(ax)) - - -def ro2om(ro): - """Rodgrigues vector to orientation matrix""" - return ax2om(ro2ax(ro)) - - -def ro2qu(ro): - """Rodrigues vector to quaternion""" - return ax2qu(ro2ax(ro)) - - -def ho2eu(ho): - """Homochoric to Euler angles""" - return ax2eu(ho2ax(ho)) - - -def ho2om(ho): - """Homochoric to orientation matrix""" - return ax2om(ho2ax(ho)) - - -def ho2ro(ho): - """Axis angle to Rodriques vector""" - return ax2ro(ho2ax(ho)) - - -def ho2qu(ho): - """Homochoric to quaternion""" - return ax2qu(ho2ax(ho)) - - -def eu2cu(eu): - """Euler angles to cubochoric""" - return ho2cu(eu2ho(eu)) - - -def om2cu(om): - """Orientation matrix to cubochoric""" - return ho2cu(om2ho(om)) - - -def om2qu(om): - """ - Orientation matrix to quaternion - - The original formulation (direct conversion) had numerical issues - """ - return ax2qu(om2ax(om)) - - -def ax2cu(ax): - """Axis angle to cubochoric""" - return ho2cu(ax2ho(ax)) - - -def ro2cu(ro): - """Rodrigues vector to cubochoric""" - return ho2cu(ro2ho(ro)) - - -def qu2cu(qu): - """Quaternion to cubochoric""" - return ho2cu(qu2ho(qu)) - - -def cu2eu(cu): - """Cubochoric to Euler angles""" - return ho2eu(cu2ho(cu)) +def cu2qu(cu): + """Cubochoric to quaternion""" + return ho2qu(cu2ho(cu)) def cu2om(cu): @@ -1576,6 +1457,11 @@ def cu2om(cu): return ho2om(cu2ho(cu)) +def cu2eu(cu): + """Cubochoric to Euler angles""" + return ho2eu(cu2ho(cu)) + + def cu2ax(cu): """Cubochoric to axis angle""" return ho2ax(cu2ho(cu)) @@ -1586,6 +1472,6 @@ def cu2ro(cu): return ho2ro(cu2ho(cu)) -def cu2qu(cu): - """Cubochoric to quaternion""" - return ho2qu(cu2ho(cu)) +def cu2ho(cu): + """Cubochoric to homochoric""" + return Lambert.CubeToBall(cu) diff --git a/python/damask/quaternion.py b/python/damask/quaternion.py new file mode 100644 index 000000000..fd681fe9b --- /dev/null +++ b/python/damask/quaternion.py @@ -0,0 +1,214 @@ +# -*- coding: UTF-8 no BOM -*- + +import numpy as np + +P = -1 # convention (see DOI:10.1088/0965-0393/23/8/083501) + +#################################################################################################### +class Quaternion: + u""" + Quaternion with basic operations + + q is the real part, p = (x, y, z) are the imaginary parts. + Definition of multiplication depends on variable P, P ∈ {-1,1}. + """ + + def __init__(self, + q = 0.0, + p = np.zeros(3,dtype=float)): + """Initializes to identity unless specified""" + self.q = float(q) + self.p = np.array(p,dtype=float) + + + def __copy__(self): + """Copy""" + return self.__class__(q=self.q, + p=self.p.copy()) + + copy = __copy__ + + + def __iter__(self): + """Components""" + return iter(self.asList()) + + def __repr__(self): + """Readable string""" + return 'Quaternion: (real={q:+.6f}, imag=<{p[0]:+.6f}, {p[1]:+.6f}, {p[2]:+.6f}>)'.format(q=self.q,p=self.p) + + + def __add__(self, other): + """Addition""" + if isinstance(other, Quaternion): + return self.__class__(q=self.q + other.q, + p=self.p + other.p) + else: + return NotImplemented + + def __iadd__(self, other): + """In-place addition""" + if isinstance(other, Quaternion): + self.q += other.q + self.p += other.p + return self + else: + return NotImplemented + + def __pos__(self): + """Unary positive operator""" + return self + + + def __sub__(self, other): + """Subtraction""" + if isinstance(other, Quaternion): + return self.__class__(q=self.q - other.q, + p=self.p - other.p) + else: + return NotImplemented + + def __isub__(self, other): + """In-place subtraction""" + if isinstance(other, Quaternion): + self.q -= other.q + self.p -= other.p + return self + else: + return NotImplemented + + def __neg__(self): + """Unary negative operator""" + return self * -1.0 + + + def __mul__(self, other): + """Multiplication with quaternion or scalar""" + if isinstance(other, Quaternion): + return self.__class__(q=self.q*other.q - np.dot(self.p,other.p), + p=self.q*other.p + other.q*self.p + P * np.cross(self.p,other.p)) + elif isinstance(other, (int, float)): + return self.__class__(q=self.q*other, + p=self.p*other) + else: + return NotImplemented + + def __imul__(self, other): + """In-place multiplication with quaternion or scalar""" + if isinstance(other, Quaternion): + self.q = self.q*other.q - np.dot(self.p,other.p) + self.p = self.q*other.p + other.q*self.p + P * np.cross(self.p,other.p) + return self + elif isinstance(other, (int, float)): + self.q *= other + self.p *= other + return self + else: + return NotImplemented + + + def __truediv__(self, other): + """Divsion with quaternion or scalar""" + if isinstance(other, Quaternion): + s = other.conjugate()/abs(other)**2. + return self.__class__(q=self.q * s, + p=self.p * s) + elif isinstance(other, (int, float)): + return self.__class__(q=self.q / other, + p=self.p / other) + else: + return NotImplemented + + def __itruediv__(self, other): + """In-place divsion with quaternion or scalar""" + if isinstance(other, Quaternion): + s = other.conjugate()/abs(other)**2. + self *= s + return self + elif isinstance(other, (int, float)): + self.q /= other + self.p /= other + return self + else: + return NotImplemented + + + def __pow__(self, exponent): + """Power""" + if isinstance(exponent, (int, float)): + omega = np.acos(self.q) + return self.__class__(q= np.cos(exponent*omega), + p=self.p * np.sin(exponent*omega)/np.sin(omega)) + else: + return NotImplemented + + def __ipow__(self, exponent): + """In-place power""" + if isinstance(exponent, (int, float)): + omega = np.acos(self.q) + self.q = np.cos(exponent*omega) + self.p *= np.sin(exponent*omega)/np.sin(omega) + return self + else: + return NotImplemented + + + def __abs__(self): + """Norm""" + return np.sqrt(self.q ** 2 + np.dot(self.p,self.p)) + + magnitude = __abs__ + + + def __eq__(self,other): + """Equal (sufficiently close) to each other""" + return np.isclose(( self-other).magnitude(),0.0) \ + or np.isclose((-self-other).magnitude(),0.0) + + def __ne__(self,other): + """Not equal (sufficiently close) to each other""" + return not self.__eq__(other) + + + def asM(self): + """Intermediate representation useful for quaternion averaging (see F. Landis Markley et al.)""" + return np.outer(self.asArray(),self.asArray()) + + def asArray(self): + """As numpy array""" + return np.array((self.q,self.p[0],self.p[1],self.p[2])) + + def asList(self): + """As list""" + return [self.q]+list(self.p) + + + def normalize(self): + """Normalizes in-place""" + d = self.magnitude() + if d > 0.0: self /= d + return self + + def normalized(self): + """Returns normalized copy""" + return self.copy().normalize() + + + def conjugate(self): + """Conjugates in-place""" + self.p = -self.p + return self + + def conjugated(self): + """Returns conjugated copy""" + return self.copy().conjugate() + + + def homomorph(self): + """Homomorphs in-place""" + self *= -1.0 + return self + + def homomorphed(self): + """Returns homomorphed copy""" + return self.copy().homomorph() diff --git a/src/CPFEM.f90 b/src/CPFEM.f90 index d2eaa7979..a1e562c24 100644 --- a/src/CPFEM.f90 +++ b/src/CPFEM.f90 @@ -47,7 +47,8 @@ module CPFEM public :: & CPFEM_general, & - CPFEM_initAll + CPFEM_initAll, & + CPFEM_results contains @@ -633,4 +634,35 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt end subroutine CPFEM_general + +!-------------------------------------------------------------------------------------------------- +!> @brief triggers writing of the results +!-------------------------------------------------------------------------------------------------- +subroutine CPFEM_results(inc,time) + use prec, only: & + pInt +#ifdef DAMASK_HDF5 + use results + use HDF5_utilities +#endif + use constitutive, only: & + constitutive_results + use crystallite, only: & + crystallite_results + + implicit none + integer(pInt), intent(in) :: inc + real(pReal), intent(in) :: time + +#ifdef DAMASK_HDF5 + call results_openJobFile + call results_addIncrement(inc,time) + call constitutive_results + call crystallite_results + call results_removeLink('current') ! ToDo: put this into closeJobFile + call results_closeJobFile +#endif + +end subroutine CPFEM_results + end module CPFEM diff --git a/src/CPFEM2.f90 b/src/CPFEM2.f90 index 45e423fe3..bf8deeb78 100644 --- a/src/CPFEM2.f90 +++ b/src/CPFEM2.f90 @@ -257,7 +257,7 @@ subroutine CPFEM_age() write(6,'(a)') '<< CPFEM >> writing restart variables of last converged step to hdf5 file' write(rankStr,'(a1,i0)')'_',worldrank - fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5','w') + fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5','a') call HDF5_write(fileHandle,material_phase, 'recordedPhase') call HDF5_write(fileHandle,crystallite_F0, 'convergedF') @@ -290,6 +290,7 @@ subroutine CPFEM_age() end subroutine CPFEM_age + !-------------------------------------------------------------------------------------------------- !> @brief triggers writing of the results !-------------------------------------------------------------------------------------------------- diff --git a/src/DAMASK_interface.f90 b/src/DAMASK_interface.f90 index 87a73af37..ebe625ef6 100644 --- a/src/DAMASK_interface.f90 +++ b/src/DAMASK_interface.f90 @@ -45,7 +45,7 @@ contains !> @brief initializes the solver by interpreting the command line arguments. Also writes !! information on computation to screen !-------------------------------------------------------------------------------------------------- -subroutine DAMASK_interface_init() +subroutine DAMASK_interface_init use, intrinsic :: & iso_fortran_env use, intrinsic :: & @@ -130,11 +130,11 @@ subroutine DAMASK_interface_init() call MPI_Init_Thread(MPI_THREAD_FUNNELED,threadLevel,mpi_err) if (mpi_err /= 0) call quit(1) if (threadLevel @details Marc subroutines used: !> @details - hypela2 !> @details - plotv +!> @details - uedinc !> @details - flux !> @details - quit !> @details Marc common blocks included: @@ -273,26 +274,20 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, & outdatedByNewInc = .false. ! no aging of state calcMode = .false. ! pretend last step was collection lastLovl = lovl ! pretend that this is NOT the first after a lovl change - !$OMP CRITICAL (write2out) - write(6,'(a,i6,1x,i2)') '<< HYPELA2 >> start of analysis..! ',m(1),nn - flush(6) - !$OMP END CRITICAL (write2out) + write(6,'(a,i6,1x,i2)') '<< HYPELA2 >> start of analysis..! ',m(1),nn + flush(6) else if (inc - theInc > 1) then ! >> restart of broken analysis << lastIncConverged = .false. ! no Jacobian backup outdatedByNewInc = .false. ! no aging of state calcMode = .true. ! pretend last step was calculation - !$OMP CRITICAL (write2out) - write(6,'(a,i6,1x,i2)') '<< HYPELA2 >> restart of analysis..! ',m(1),nn - flush(6) - !$OMP END CRITICAL (write2out) + write(6,'(a,i6,1x,i2)') '<< HYPELA2 >> restart of analysis..! ',m(1),nn + flush(6) else ! >> just the next inc << lastIncConverged = .true. ! request Jacobian backup outdatedByNewInc = .true. ! request aging of state calcMode = .true. ! assure last step was calculation - !$OMP CRITICAL (write2out) - write(6,'(a,i6,1x,i2)') '<< HYPELA2 >> new increment..! ',m(1),nn - flush(6) - !$OMP END CRITICAL (write2out) + write(6,'(a,i6,1x,i2)') '<< HYPELA2 >> new increment..! ',m(1),nn + flush(6) endif else if ( timinc < theDelta ) then ! >> cutBack << lastIncConverged = .false. ! no Jacobian backup @@ -300,10 +295,8 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, & terminallyIll = .false. cycleCounter = -1 ! first calc step increments this to cycle = 0 calcMode = .true. ! pretend last step was calculation - !$OMP CRITICAL (write2out) - write(6,'(a,i6,1x,i2)') '<< HYPELA2 >> cutback detected..! ',m(1),nn - flush(6) - !$OMP END CRITICAL (write2out) + write(6,'(a,i6,1x,i2)') '<< HYPELA2 >> cutback detected..! ',m(1),nn + flush(6) endif ! convergence treatment end @@ -365,7 +358,6 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, & lastLovl = lovl ! record lovl call CPFEM_general(computationMode,usePingPong,ffn,ffn1,t(1),timinc,m(1),nn,stress,ddsdde) - ! Mandel: 11, 22, 33, SQRT(2)*12, SQRT(2)*23, SQRT(2)*13 ! Marc: 11, 22, 33, 12, 23, 13 ! Marc: 11, 22, 33, 12 @@ -407,6 +399,26 @@ subroutine flux(f,ts,n,time) end subroutine flux +!-------------------------------------------------------------------------------------------------- +!> @brief sets user defined output variables for Marc +!> @details select a variable contour plotting (user subroutine). +!-------------------------------------------------------------------------------------------------- +subroutine uedinc(inc,incsub) + use prec, only: & + pReal, & + pInt + use CPFEM, only: & + CPFEM_results + + implicit none + integer, intent(in) :: inc, incsub +#include QUOTE(PASTE(./MarcInclude/creeps,Marc4DAMASK)) ! creeps is needed for timinc (time increment) + + call CPFEM_results(inc,cptim) + +end subroutine uedinc + + !-------------------------------------------------------------------------------------------------- !> @brief sets user defined output variables for Marc !> @details select a variable contour plotting (user subroutine). diff --git a/src/grid/DAMASK_grid.f90 b/src/grid/DAMASK_grid.f90 index f2f52bb2f..2adf72f89 100644 --- a/src/grid/DAMASK_grid.f90 +++ b/src/grid/DAMASK_grid.f90 @@ -77,6 +77,7 @@ program DAMASK_spectral use grid_mech_FEM use grid_damage_spectral use grid_thermal_spectral + use HDF5_utilities use results implicit none @@ -155,8 +156,6 @@ program DAMASK_spectral write(6,'(/,a)') ' Shanthraj et al., Handbook of Mechanics of Materials, 2019' write(6,'(a)') ' https://doi.org/10.1007/978-981-10-6855-3_80' - call results_openJobFile() - call results_closeJobFile() !-------------------------------------------------------------------------------------------------- ! initialize field solver information nActiveFields = 1 @@ -301,7 +300,7 @@ program DAMASK_spectral reportAndCheck: if (worldrank == 0) then write (loadcase_string, '(i6)' ) currentLoadCase - write(6,'(1x,a,i6)') 'load case: ', currentLoadCase + write(6,'(/,1x,a,i6)') 'load case: ', currentLoadCase if (.not. newLoadCase%followFormerTrajectory) write(6,'(2x,a)') 'drop guessing along trajectory' if (newLoadCase%deformation%myType == 'l') then do j = 1_pInt, 3_pInt @@ -348,10 +347,10 @@ program DAMASK_spectral if (newLoadCase%time < 0.0_pReal) errorID = 834_pInt ! negative time increment write(6,'(2x,a,f12.6)') 'time: ', newLoadCase%time if (newLoadCase%incs < 1_pInt) errorID = 835_pInt ! non-positive incs count - write(6,'(2x,a,i5)') 'increments: ', newLoadCase%incs + write(6,'(2x,a,i5)') 'increments: ', newLoadCase%incs if (newLoadCase%outputfrequency < 1_pInt) errorID = 836_pInt ! non-positive result frequency - write(6,'(2x,a,i5)') 'output frequency: ', newLoadCase%outputfrequency - write(6,'(2x,a,i5,/)') 'restart frequency: ', newLoadCase%restartfrequency + write(6,'(2x,a,i5)') 'output frequency: ', newLoadCase%outputfrequency + write(6,'(2x,a,i5)') 'restart frequency: ', newLoadCase%restartfrequency if (errorID > 0_pInt) call IO_error(error_ID = errorID, ext_msg = loadcase_string) ! exit with error message endif reportAndCheck loadCases = [loadCases,newLoadCase] ! load case is ok, append it @@ -359,8 +358,9 @@ program DAMASK_spectral close(fileUnit) call results_openJobFile - call results_addAttribute('grid',grid,'mapping') - call results_addAttribute('size',geomSize,'mapping') + call HDF5_closeGroup(results_addGroup('geometry')) + call results_addAttribute('grid',grid,'geometry') + call results_addAttribute('size',geomSize,'geometry') call results_closeJobFile !-------------------------------------------------------------------------------------------------- diff --git a/src/grid/grid_damage_spectral.f90 b/src/grid/grid_damage_spectral.f90 index c5e9a254b..3ce37c5ff 100644 --- a/src/grid/grid_damage_spectral.f90 +++ b/src/grid/grid_damage_spectral.f90 @@ -64,7 +64,6 @@ subroutine grid_damage_spectral_init worldsize, & petsc_options - implicit none PetscInt, dimension(worldsize) :: localK integer :: i, j, k, cell DM :: damage_grid @@ -164,7 +163,6 @@ function grid_damage_spectral_solution(timeinc,timeinc_old,loadCaseTime) result( use damage_nonlocal, only: & damage_nonlocal_putNonLocalDamage - implicit none real(pReal), intent(in) :: & timeinc, & !< increment in time for current solution timeinc_old, & !< increment in time of last increment @@ -236,7 +234,6 @@ subroutine grid_damage_spectral_forward damage_nonlocal_getDiffusion33, & damage_nonlocal_getMobility - implicit none integer :: i, j, k, cell DM :: dm_local PetscScalar, dimension(:,:,:), pointer :: x_scal @@ -301,7 +298,6 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr) damage_nonlocal_getDiffusion33, & damage_nonlocal_getMobility - implicit none DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: & in PetscScalar, dimension( & diff --git a/src/grid/grid_mech_FEM.f90 b/src/grid/grid_mech_FEM.f90 index e31d93637..2feec76df 100644 --- a/src/grid/grid_mech_FEM.f90 +++ b/src/grid/grid_mech_FEM.f90 @@ -7,65 +7,66 @@ module grid_mech_FEM #include #include - use PETScdmda - use PETScsnes - use prec, only: & - pInt, & - pReal - use math, only: & - math_I3 - use spectral_utilities, only: & - tSolutionState, & - tSolutionParams - - implicit none - private + use DAMASK_interface + use HDF5_utilities + use PETScdmda + use PETScsnes + use prec, only: & + pReal + use math, only: & + math_I3 + use spectral_utilities, only: & + tSolutionState, & + tSolutionParams + + implicit none + private !-------------------------------------------------------------------------------------------------- ! derived types - type(tSolutionParams), private :: params + type(tSolutionParams), private :: params !-------------------------------------------------------------------------------------------------- ! PETSc data - DM, private :: mech_grid - SNES, private :: mech_snes - Vec, private :: solution_current, solution_lastInc, solution_rate + DM, private :: mech_grid + SNES, private :: mech_snes + Vec, private :: solution_current, solution_lastInc, solution_rate !-------------------------------------------------------------------------------------------------- ! common pointwise data - real(pReal), private, dimension(:,:,:,:,:), allocatable :: F, P_current, F_lastInc - real(pReal), private :: detJ - real(pReal), private, dimension(3) :: delta - real(pReal), private, dimension(3,8) :: BMat - real(pReal), private, dimension(8,8) :: HGMat - PetscInt, private :: xstart,ystart,zstart,xend,yend,zend + real(pReal), private, dimension(:,:,:,:,:), allocatable :: F, P_current, F_lastInc + real(pReal), private :: detJ + real(pReal), private, dimension(3) :: delta + real(pReal), private, dimension(3,8) :: BMat + real(pReal), private, dimension(8,8) :: HGMat + PetscInt, private :: xstart,ystart,zstart,xend,yend,zend !-------------------------------------------------------------------------------------------------- ! stress, stiffness and compliance average etc. - real(pReal), private, dimension(3,3) :: & - F_aimDot = 0.0_pReal, & !< assumed rate of average deformation gradient - F_aim = math_I3, & !< current prescribed deformation gradient - F_aim_lastIter = math_I3, & - F_aim_lastInc = math_I3, & !< previous average deformation gradient - P_av = 0.0_pReal !< average 1st Piola--Kirchhoff stress - - character(len=1024), private :: incInfo !< time and increment information - - real(pReal), private, dimension(3,3,3,3) :: & - C_volAvg = 0.0_pReal, & !< current volume average stiffness - C_volAvgLastInc = 0.0_pReal, & !< previous volume average stiffness - S = 0.0_pReal !< current compliance (filled up with zeros) - - real(pReal), private :: & - err_BC !< deviation from stress BC - - integer(pInt), private :: & - totalIter = 0_pInt !< total iteration in current increment - - public :: & - grid_mech_FEM_init, & - grid_mech_FEM_solution, & - grid_mech_FEM_forward + real(pReal), private, dimension(3,3) :: & + F_aimDot = 0.0_pReal, & !< assumed rate of average deformation gradient + F_aim = math_I3, & !< current prescribed deformation gradient + F_aim_lastIter = math_I3, & + F_aim_lastInc = math_I3, & !< previous average deformation gradient + P_av = 0.0_pReal !< average 1st Piola--Kirchhoff stress + + character(len=1024), private :: incInfo !< time and increment information + + real(pReal), private, dimension(3,3,3,3) :: & + C_volAvg = 0.0_pReal, & !< current volume average stiffness + C_volAvgLastInc = 0.0_pReal, & !< previous volume average stiffness + S = 0.0_pReal !< current compliance (filled up with zeros) + + real(pReal), private :: & + err_BC !< deviation from stress BC + + integer, private :: & + totalIter = 0 !< total iteration in current increment + + public :: & + grid_mech_FEM_init, & + grid_mech_FEM_solution, & + grid_mech_FEM_forward contains @@ -77,172 +78,163 @@ subroutine grid_mech_FEM_init IO_intOut, & IO_error, & IO_open_jobFile_binary - use FEsolving, only: & - restartInc - use numerics, only: & - worldrank, & - worldsize, & - petsc_options - use homogenization, only: & - materialpoint_F0 - use DAMASK_interface, only: & - getSolverJobName - use spectral_utilities, only: & - utilities_constitutiveResponse, & - utilities_updateIPcoords, & - wgt - use mesh, only: & - geomSize, & - grid, & - grid3 - use math, only: & - math_invSym3333 - - implicit none - real(pReal) :: HGCoeff = 0e-2_pReal - PetscInt, dimension(:), allocatable :: localK - real(pReal), dimension(3,3) :: & - temp33_Real = 0.0_pReal - real(pReal), dimension(4,8) :: & - HGcomp = reshape([ 1.0_pReal, 1.0_pReal, 1.0_pReal,-1.0_pReal, & - 1.0_pReal,-1.0_pReal,-1.0_pReal, 1.0_pReal, & - -1.0_pReal, 1.0_pReal,-1.0_pReal, 1.0_pReal, & - -1.0_pReal,-1.0_pReal, 1.0_pReal,-1.0_pReal, & - -1.0_pReal,-1.0_pReal, 1.0_pReal, 1.0_pReal, & - -1.0_pReal, 1.0_pReal,-1.0_pReal,-1.0_pReal, & - 1.0_pReal,-1.0_pReal,-1.0_pReal,-1.0_pReal, & - 1.0_pReal, 1.0_pReal, 1.0_pReal, 1.0_pReal], [4,8]) - PetscErrorCode :: ierr - integer(pInt) :: rank - integer :: fileUnit - character(len=1024) :: rankStr - real(pReal), dimension(3,3,3,3) :: devNull - PetscScalar, pointer, dimension(:,:,:,:) :: & - u_current,u_lastInc - - write(6,'(/,a)') ' <<<+- grid_mech_FEM init -+>>>' + use FEsolving, only: & + restartInc + use numerics, only: & + worldrank, & + worldsize, & + petsc_options + use homogenization, only: & + materialpoint_F0 + use DAMASK_interface, only: & + getSolverJobName + use spectral_utilities, only: & + utilities_constitutiveResponse, & + utilities_updateIPcoords, & + wgt + use mesh, only: & + geomSize, & + grid, & + grid3 + use math, only: & + math_invSym3333 + + real(pReal) :: HGCoeff = 0e-2_pReal + PetscInt, dimension(:), allocatable :: localK + real(pReal), dimension(3,3) :: & + temp33_Real = 0.0_pReal + real(pReal), dimension(4,8) :: & + HGcomp = reshape([ 1.0_pReal, 1.0_pReal, 1.0_pReal,-1.0_pReal, & + 1.0_pReal,-1.0_pReal,-1.0_pReal, 1.0_pReal, & + -1.0_pReal, 1.0_pReal,-1.0_pReal, 1.0_pReal, & + -1.0_pReal,-1.0_pReal, 1.0_pReal,-1.0_pReal, & + -1.0_pReal,-1.0_pReal, 1.0_pReal, 1.0_pReal, & + -1.0_pReal, 1.0_pReal,-1.0_pReal,-1.0_pReal, & + 1.0_pReal,-1.0_pReal,-1.0_pReal,-1.0_pReal, & + 1.0_pReal, 1.0_pReal, 1.0_pReal, 1.0_pReal], [4,8]) + PetscErrorCode :: ierr + integer :: rank + integer(HID_T) :: fileHandle + character(len=1024) :: rankStr + real(pReal), dimension(3,3,3,3) :: devNull + PetscScalar, pointer, dimension(:,:,:,:) :: & + u_current,u_lastInc + + write(6,'(/,a)') ' <<<+- grid_mech_FEM init -+>>>' !-------------------------------------------------------------------------------------------------- ! set default and user defined options for PETSc - call PETScOptionsInsertString(PETSC_NULL_OPTIONS,'-mech_snes_type newtonls -mech_ksp_type fgmres & - &-mech_ksp_max_it 25 -mech_pc_type ml -mech_mg_levels_ksp_type chebyshev',ierr) - CHKERRQ(ierr) - call PETScOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_options),ierr) - CHKERRQ(ierr) + call PETScOptionsInsertString(PETSC_NULL_OPTIONS,'-mech_snes_type newtonls -mech_ksp_type fgmres & + &-mech_ksp_max_it 25 -mech_pc_type ml -mech_mg_levels_ksp_type chebyshev',ierr) + CHKERRQ(ierr) + call PETScOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_options),ierr) + CHKERRQ(ierr) !-------------------------------------------------------------------------------------------------- ! allocate global fields - allocate (F (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) - allocate (P_current (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) - allocate (F_lastInc (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) + allocate(F (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) + allocate(P_current (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) + allocate(F_lastInc (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) !-------------------------------------------------------------------------------------------------- ! initialize solver specific parts of PETSc - call SNESCreate(PETSC_COMM_WORLD,mech_snes,ierr); CHKERRQ(ierr) - call SNESSetOptionsPrefix(mech_snes,'mech_',ierr);CHKERRQ(ierr) - allocate(localK(worldsize), source = 0); localK(worldrank+1) = grid3 - do rank = 1, worldsize - call MPI_Bcast(localK(rank),1,MPI_INTEGER,rank-1,PETSC_COMM_WORLD,ierr) - enddo - call DMDACreate3d(PETSC_COMM_WORLD, & - DM_BOUNDARY_PERIODIC, DM_BOUNDARY_PERIODIC, DM_BOUNDARY_PERIODIC, & - DMDA_STENCIL_BOX, & - grid(1),grid(2),grid(3), & - 1, 1, worldsize, & - 3, 1, & - [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 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) - call DMSNESSetFunctionLocal(mech_grid,formResidual,PETSC_NULL_SNES,ierr) - CHKERRQ(ierr) - call DMSNESSetJacobianLocal(mech_grid,formJacobian,PETSC_NULL_SNES,ierr) - CHKERRQ(ierr) - call SNESSetConvergenceTest(mech_snes,converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,ierr) - CHKERRQ(ierr) ! specify custom convergence check function "_converged" - call SNESSetMaxLinearSolveFailures(mech_snes, huge(1), ierr); CHKERRQ(ierr) ! ignore linear solve failures - call SNESSetFromOptions(mech_snes,ierr); CHKERRQ(ierr) ! pull it all together with additional cli arguments + call SNESCreate(PETSC_COMM_WORLD,mech_snes,ierr); CHKERRQ(ierr) + call SNESSetOptionsPrefix(mech_snes,'mech_',ierr);CHKERRQ(ierr) + allocate(localK(worldsize), source = 0); localK(worldrank+1) = grid3 + do rank = 1, worldsize + call MPI_Bcast(localK(rank),1,MPI_INTEGER,rank-1,PETSC_COMM_WORLD,ierr) + enddo + call DMDACreate3d(PETSC_COMM_WORLD, & + DM_BOUNDARY_PERIODIC, DM_BOUNDARY_PERIODIC, DM_BOUNDARY_PERIODIC, & + DMDA_STENCIL_BOX, & + grid(1),grid(2),grid(3), & + 1, 1, worldsize, & + 3, 1, & + [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 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) + call DMSNESSetFunctionLocal(mech_grid,formResidual,PETSC_NULL_SNES,ierr) + CHKERRQ(ierr) + call DMSNESSetJacobianLocal(mech_grid,formJacobian,PETSC_NULL_SNES,ierr) + CHKERRQ(ierr) + call SNESSetConvergenceTest(mech_snes,converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,ierr) ! specify custom convergence check function "_converged" + CHKERRQ(ierr) + call SNESSetMaxLinearSolveFailures(mech_snes, huge(1), ierr); CHKERRQ(ierr) ! ignore linear solve failures + call SNESSetFromOptions(mech_snes,ierr); CHKERRQ(ierr) ! pull it all together with additional cli arguments !-------------------------------------------------------------------------------------------------- ! init fields - call VecSet(solution_current,0.0_pReal,ierr);CHKERRQ(ierr) - call VecSet(solution_lastInc,0.0_pReal,ierr);CHKERRQ(ierr) - call VecSet(solution_rate ,0.0_pReal,ierr);CHKERRQ(ierr) - call DMDAVecGetArrayF90(mech_grid,solution_current,u_current,ierr); CHKERRQ(ierr) - call DMDAVecGetArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr); CHKERRQ(ierr) - - call DMDAGetCorners(mech_grid,xstart,ystart,zstart,xend,yend,zend,ierr) ! local grid extent - CHKERRQ(ierr) - xend = xstart+xend-1 - yend = ystart+yend-1 - zend = zstart+zend-1 - delta = geomSize/real(grid,pReal) ! grid spacing - detJ = product(delta) ! cell volume - - BMat = reshape(real([-1.0_pReal/delta(1),-1.0_pReal/delta(2),-1.0_pReal/delta(3), & - 1.0_pReal/delta(1),-1.0_pReal/delta(2),-1.0_pReal/delta(3), & - -1.0_pReal/delta(1), 1.0_pReal/delta(2),-1.0_pReal/delta(3), & - 1.0_pReal/delta(1), 1.0_pReal/delta(2),-1.0_pReal/delta(3), & - -1.0_pReal/delta(1),-1.0_pReal/delta(2), 1.0_pReal/delta(3), & - 1.0_pReal/delta(1),-1.0_pReal/delta(2), 1.0_pReal/delta(3), & - -1.0_pReal/delta(1), 1.0_pReal/delta(2), 1.0_pReal/delta(3), & - 1.0_pReal/delta(1), 1.0_pReal/delta(2), 1.0_pReal/delta(3)],pReal), [3,8])/4.0_pReal ! shape function derivative matrix - - HGMat = matmul(transpose(HGcomp),HGcomp) & - * HGCoeff*(delta(1)*delta(2) + delta(2)*delta(3) + delta(3)*delta(1))/16.0_pReal ! hourglass stabilization matrix + call VecSet(solution_current,0.0_pReal,ierr);CHKERRQ(ierr) + call VecSet(solution_lastInc,0.0_pReal,ierr);CHKERRQ(ierr) + call VecSet(solution_rate ,0.0_pReal,ierr);CHKERRQ(ierr) + call DMDAVecGetArrayF90(mech_grid,solution_current,u_current,ierr); CHKERRQ(ierr) + call DMDAVecGetArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr); CHKERRQ(ierr) + + call DMDAGetCorners(mech_grid,xstart,ystart,zstart,xend,yend,zend,ierr) ! local grid extent + CHKERRQ(ierr) + xend = xstart+xend-1 + yend = ystart+yend-1 + zend = zstart+zend-1 + delta = geomSize/real(grid,pReal) ! grid spacing + detJ = product(delta) ! cell volume + + BMat = reshape(real([-1.0_pReal/delta(1),-1.0_pReal/delta(2),-1.0_pReal/delta(3), & + 1.0_pReal/delta(1),-1.0_pReal/delta(2),-1.0_pReal/delta(3), & + -1.0_pReal/delta(1), 1.0_pReal/delta(2),-1.0_pReal/delta(3), & + 1.0_pReal/delta(1), 1.0_pReal/delta(2),-1.0_pReal/delta(3), & + -1.0_pReal/delta(1),-1.0_pReal/delta(2), 1.0_pReal/delta(3), & + 1.0_pReal/delta(1),-1.0_pReal/delta(2), 1.0_pReal/delta(3), & + -1.0_pReal/delta(1), 1.0_pReal/delta(2), 1.0_pReal/delta(3), & + 1.0_pReal/delta(1), 1.0_pReal/delta(2), 1.0_pReal/delta(3)],pReal), [3,8])/4.0_pReal ! shape function derivative matrix + + HGMat = matmul(transpose(HGcomp),HGcomp) & + * HGCoeff*(delta(1)*delta(2) + delta(2)*delta(3) + delta(3)*delta(1))/16.0_pReal ! hourglass stabilization matrix !-------------------------------------------------------------------------------------------------- ! init fields - restart: if (restartInc > 0) then - write(6,'(/,a,'//IO_intOut(restartInc)//',a)') 'reading values of increment ', restartInc, ' from file' - - fileUnit = IO_open_jobFile_binary('F_aim') - read(fileUnit) F_aim; close(fileUnit) - fileUnit = IO_open_jobFile_binary('F_aim_lastInc') - read(fileUnit) F_aim_lastInc; close(fileUnit) - fileUnit = IO_open_jobFile_binary('F_aimDot') - read(fileUnit) F_aimDot; close(fileUnit) - - write(rankStr,'(a1,i0)')'_',worldrank - - fileUnit = IO_open_jobFile_binary('F'//trim(rankStr)) - read(fileUnit) F; close (fileUnit) - fileUnit = IO_open_jobFile_binary('F_lastInc'//trim(rankStr)) - read(fileUnit) F_lastInc; close (fileUnit) - fileUnit = IO_open_jobFile_binary('u'//trim(rankStr)) - read(fileUnit) u_current; close (fileUnit) - fileUnit = IO_open_jobFile_binary('u_lastInc'//trim(rankStr)) - read(fileUnit) u_lastInc; close (fileUnit) - - elseif (restartInc == 0) then restart - F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity - F = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) - endif restart - materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent - call Utilities_updateIPcoords(F) - call Utilities_constitutiveResponse(P_current,temp33_Real,C_volAvg,devNull, & ! stress field, stress avg, global average of stiffness and (min+max)/2 - F, & ! target F - 0.0_pReal, & ! time increment - math_I3) ! no rotation of boundary condition - call DMDAVecRestoreArrayF90(mech_grid,solution_current,u_current,ierr) - CHKERRQ(ierr) - call DMDAVecRestoreArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr) - CHKERRQ(ierr) - - restartRead: if (restartInc > 0_pInt) then - write(6,'(/,a,'//IO_intOut(restartInc)//',a)') 'reading more values of increment ', restartInc, ' from file' - fileUnit = IO_open_jobFile_binary('C_volAvg') - read(fileUnit) C_volAvg; close(fileUnit) - fileUnit = IO_open_jobFile_binary('C_volAvgLastInv') - read(fileUnit) C_volAvgLastInc; close(fileUnit) - endif restartRead + restart: if (restartInc > 0) then + write(6,'(/,a,'//IO_intOut(restartInc)//',a)') 'reading values of increment ', restartInc, ' from file' + + write(rankStr,'(a1,i0)')'_',worldrank + fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5') + + call HDF5_read(fileHandle,F_aim, 'F_aim') + call HDF5_read(fileHandle,F_aim_lastInc,'F_aim_lastInc') + call HDF5_read(fileHandle,F_aimDot, 'F_aimDot') + call HDF5_read(fileHandle,F, 'F') + call HDF5_read(fileHandle,F_lastInc, 'F_lastInc') + call HDF5_read(fileHandle,u_current, 'u') + call HDF5_read(fileHandle,u_lastInc, 'u_lastInc') + + elseif (restartInc == 0) then restart + F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity + F = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) + endif restart + materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent + call utilities_updateIPcoords(F) + call utilities_constitutiveResponse(P_current,temp33_Real,C_volAvg,devNull, & ! stress field, stress avg, global average of stiffness and (min+max)/2 + F, & ! target F + 0.0_pReal, & ! time increment + math_I3) ! no rotation of boundary condition + call DMDAVecRestoreArrayF90(mech_grid,solution_current,u_current,ierr) + CHKERRQ(ierr) + call DMDAVecRestoreArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr) + CHKERRQ(ierr) + + restartRead: if (restartInc > 0) then + write(6,'(/,a,'//IO_intOut(restartInc)//',a)') 'reading more values of increment ', restartInc, ' from file' + call HDF5_read(fileHandle,C_volAvg, 'C_volAvg') + call HDF5_read(fileHandle,C_volAvgLastInc,'C_volAvgLastInc') + call HDF5_closeFile(fileHandle) + endif restartRead end subroutine grid_mech_FEM_init @@ -251,60 +243,57 @@ end subroutine grid_mech_FEM_init !> @brief solution for the FEM scheme with internal iterations !-------------------------------------------------------------------------------------------------- function grid_mech_FEM_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC) result(solution) - use IO, only: & - IO_error - use spectral_utilities, only: & - tBoundaryCondition, & - utilities_maskedCompliance - use FEsolving, only: & - restartWrite, & - terminallyIll - - implicit none + use IO, only: & + IO_error + use spectral_utilities, only: & + tBoundaryCondition, & + utilities_maskedCompliance + use FEsolving, only: & + restartWrite, & + terminallyIll !-------------------------------------------------------------------------------------------------- ! input data for solution - character(len=*), intent(in) :: & - incInfoIn - real(pReal), intent(in) :: & - timeinc, & !< time increment of current solution - timeinc_old !< time increment of last successful increment - type(tBoundaryCondition), intent(in) :: & - stress_BC - real(pReal), dimension(3,3), intent(in) :: rotation_BC - type(tSolutionState) :: & - solution - + character(len=*), intent(in) :: & + incInfoIn + real(pReal), intent(in) :: & + timeinc, & !< time increment of current solution + timeinc_old !< time increment of last successful increment + type(tBoundaryCondition), intent(in) :: & + stress_BC + real(pReal), dimension(3,3), intent(in) :: rotation_BC + type(tSolutionState) :: & + solution !-------------------------------------------------------------------------------------------------- ! PETSc Data - PetscErrorCode :: ierr - SNESConvergedReason :: reason - - incInfo = incInfoIn + PetscErrorCode :: ierr + SNESConvergedReason :: reason + + incInfo = incInfoIn !-------------------------------------------------------------------------------------------------- ! update stiffness (and gamma operator) - S = Utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg) + S = utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg) !-------------------------------------------------------------------------------------------------- ! set module wide available data - params%stress_mask = stress_BC%maskFloat - params%stress_BC = stress_BC%values - params%rotation_BC = rotation_BC - params%timeinc = timeinc - params%timeincOld = timeinc_old + params%stress_mask = stress_BC%maskFloat + params%stress_BC = stress_BC%values + params%rotation_BC = rotation_BC + params%timeinc = timeinc + params%timeincOld = timeinc_old !-------------------------------------------------------------------------------------------------- ! solve BVP - call SNESsolve(mech_snes,PETSC_NULL_VEC,solution_current,ierr);CHKERRQ(ierr) + call SNESsolve(mech_snes,PETSC_NULL_VEC,solution_current,ierr);CHKERRQ(ierr) !-------------------------------------------------------------------------------------------------- ! check convergence - call SNESGetConvergedReason(mech_snes,reason,ierr);CHKERRQ(ierr) - - solution%converged = reason > 0 - solution%iterationsNeeded = totalIter - solution%termIll = terminallyIll - terminallyIll = .false. + call SNESGetConvergedReason(mech_snes,reason,ierr);CHKERRQ(ierr) + + solution%converged = reason > 0 + solution%iterationsNeeded = totalIter + solution%termIll = terminallyIll + terminallyIll = .false. end function grid_mech_FEM_solution @@ -335,61 +324,52 @@ subroutine grid_mech_FEM_forward(guess,timeinc,timeinc_old,loadCaseTime,deformat use FEsolving, only: & restartWrite - implicit none logical, intent(in) :: & guess real(pReal), intent(in) :: & timeinc_old, & timeinc, & - loadCaseTime !< remaining time of current load case + loadCaseTime !< remaining time of current load case type(tBoundaryCondition), intent(in) :: & stress_BC, & deformation_BC real(pReal), dimension(3,3), intent(in) :: & rotation_BC PetscErrorCode :: ierr - integer :: fileUnit + integer(HID_T) :: fileHandle character(len=32) :: rankStr PetscScalar, pointer, dimension(:,:,:,:) :: & u_current,u_lastInc - call DMDAVecGetArrayF90(mech_grid,solution_current,u_current,ierr); CHKERRQ(ierr) - call DMDAVecGetArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr); CHKERRQ(ierr) + call DMDAVecGetArrayF90(mech_grid,solution_current,u_current,ierr); CHKERRQ(ierr) + call DMDAVecGetArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr); CHKERRQ(ierr) if (cutBack) then C_volAvg = C_volAvgLastInc ! QUESTION: where is this required? else !-------------------------------------------------------------------------------------------------- ! restart information for spectral solver - if (restartWrite) then ! QUESTION: where is this logical properly set? - write(6,'(/,a)') ' writing converged results for restart' - flush(6) - - if (worldrank == 0) then - fileUnit = IO_open_jobFile_binary('C_volAvg','w') - write(fileUnit) C_volAvg; close(fileUnit) - fileUnit = IO_open_jobFile_binary('C_volAvgLastInv','w') - write(fileUnit) C_volAvgLastInc; close(fileUnit) - fileUnit = IO_open_jobFile_binary('F_aim','w') - write(fileUnit) F_aim; close(fileUnit) - fileUnit = IO_open_jobFile_binary('F_aim_lastInc','w') - write(fileUnit) F_aim_lastInc; close(fileUnit) - fileUnit = IO_open_jobFile_binary('F_aimDot','w') - write(fileUnit) F_aimDot; close(fileUnit) - endif - + if (restartWrite) then + write(6,'(/,a)') ' writing converged results for restart';flush(6) + write(rankStr,'(a1,i0)')'_',worldrank - fileUnit = IO_open_jobFile_binary('F'//trim(rankStr),'w') - write(fileUnit) F; close (fileUnit) - fileUnit = IO_open_jobFile_binary('F_lastInc'//trim(rankStr),'w') - write(fileUnit) F_lastInc; close (fileUnit) - fileUnit = IO_open_jobFile_binary('u'//trim(rankStr),'w') - write(fileUnit) u_current; close (fileUnit) - fileUnit = IO_open_jobFile_binary('u_lastInc'//trim(rankStr),'w') - write(fileUnit) u_lastInc; close (fileUnit) + fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5','w') + + call HDF5_write(fileHandle,F_aim, 'F_aim') + call HDF5_write(fileHandle,F_aim_lastInc, 'F_aim_lastInc') + call HDF5_write(fileHandle,F_aimDot, 'F_aimDot') + call HDF5_write(fileHandle,F, 'F') + call HDF5_write(fileHandle,F_lastInc, 'F_lastInc') + call HDF5_write(fileHandle,u_current, 'u') + call HDF5_write(fileHandle,u_lastInc, 'u_lastInc') + + call HDF5_write(fileHandle,C_volAvg, 'C_volAvg') + call HDF5_write(fileHandle,C_volAvgLastInc,'C_volAvgLastInc') + + call HDF5_closeFile(fileHandle) endif - call CPFEM_age() ! age state and kinematics + call CPFEM_age ! age state and kinematics call utilities_updateIPcoords(F) C_volAvgLastInc = C_volAvg @@ -399,13 +379,13 @@ subroutine grid_mech_FEM_forward(guess,timeinc,timeinc_old,loadCaseTime,deformat !-------------------------------------------------------------------------------------------------- ! calculate rate for aim - if (deformation_BC%myType=='l') then ! calculate F_aimDot from given L and current F + if (deformation_BC%myType=='l') then ! calculate F_aimDot from given L and current F F_aimDot = & F_aimDot + deformation_BC%maskFloat * matmul(deformation_BC%values, F_aim_lastInc) - elseif(deformation_BC%myType=='fdot') then ! F_aimDot is prescribed + elseif(deformation_BC%myType=='fdot') then ! F_aimDot is prescribed F_aimDot = & F_aimDot + deformation_BC%maskFloat * deformation_BC%values - elseif (deformation_BC%myType=='f') then ! aim at end of load case is prescribed + elseif (deformation_BC%myType=='f') then ! aim at end of load case is prescribed F_aimDot = & F_aimDot + deformation_BC%maskFloat * (deformation_BC%values - F_aim_lastInc)/loadCaseTime endif @@ -420,20 +400,20 @@ subroutine grid_mech_FEM_forward(guess,timeinc,timeinc_old,loadCaseTime,deformat endif call VecCopy(solution_current,solution_lastInc,ierr); CHKERRQ(ierr) - F_lastInc = F ! winding F forward - materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent + F_lastInc = F ! winding F forward + materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent endif !-------------------------------------------------------------------------------------------------- ! update average and local deformation gradients - F_aim = F_aim_lastInc + F_aimDot * timeinc - call VecAXPY(solution_current,timeinc,solution_rate,ierr); CHKERRQ(ierr) + F_aim = F_aim_lastInc + F_aimDot * timeinc + call VecAXPY(solution_current,timeinc,solution_rate,ierr); CHKERRQ(ierr) - call DMDAVecRestoreArrayF90(mech_grid,solution_current,u_current,ierr) - CHKERRQ(ierr) - call DMDAVecRestoreArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr) - CHKERRQ(ierr) + call DMDAVecRestoreArrayF90(mech_grid,solution_current,u_current,ierr) + CHKERRQ(ierr) + call DMDAVecRestoreArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr) + CHKERRQ(ierr) end subroutine grid_mech_FEM_forward @@ -441,59 +421,57 @@ end subroutine grid_mech_FEM_forward !-------------------------------------------------------------------------------------------------- !> @brief convergence check !-------------------------------------------------------------------------------------------------- -subroutine converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr) -use mesh -use spectral_utilities - use numerics, only: & - itmax, & - itmin, & - err_div_tolRel, & - err_div_tolAbs, & - err_stress_tolRel, & - err_stress_tolAbs - use FEsolving, only: & - terminallyIll +subroutine converged(snes_local,PETScIter,devNull1,devNull2,fnorm,reason,dummy,ierr) + use mesh + use spectral_utilities + use numerics, only: & + itmax, & + itmin, & + err_div_tolRel, & + err_div_tolAbs, & + err_stress_tolRel, & + err_stress_tolAbs + use FEsolving, only: & + terminallyIll - implicit none - SNES :: snes_local - PetscInt :: PETScIter - PetscReal :: & - xnorm, & ! not used - snorm, & ! not used - fnorm - SNESConvergedReason :: reason - PetscObject :: dummy - PetscErrorCode :: ierr - real(pReal) :: & - err_div, & - divTol, & - BCTol + SNES :: snes_local + PetscInt, intent(in) :: PETScIter + PetscReal, intent(in) :: & + devNull1, & + devNull2, & + fnorm + SNESConvergedReason :: reason + PetscObject :: dummy + PetscErrorCode :: ierr + real(pReal) :: & + err_div, & + divTol, & + BCTol - err_div = fnorm*sqrt(wgt)*geomSize(1)/scaledGeomSize(1)/detJ - divTol = max(maxval(abs(P_av))*err_div_tolRel ,err_div_tolAbs) - BCTol = max(maxval(abs(P_av))*err_stress_tolRel,err_stress_tolAbs) + err_div = fnorm*sqrt(wgt)*geomSize(1)/scaledGeomSize(1)/detJ + divTol = max(maxval(abs(P_av))*err_div_tolRel ,err_div_tolAbs) + BCTol = max(maxval(abs(P_av))*err_stress_tolRel,err_stress_tolAbs) - - if ((totalIter >= itmin .and. & - all([ err_div/divTol, & - err_BC /BCTol ] < 1.0_pReal)) & - .or. terminallyIll) then - reason = 1 - elseif (totalIter >= itmax) then - reason = -1 - else - reason = 0 - endif + if ((totalIter >= itmin .and. & + all([ err_div/divTol, & + err_BC /BCTol ] < 1.0_pReal)) & + .or. terminallyIll) then + reason = 1 + elseif (totalIter >= itmax) then + reason = -1 + else + reason = 0 + endif !-------------------------------------------------------------------------------------------------- ! report - write(6,'(1/,a)') ' ... reporting .............................................................' - write(6,'(1/,a,f12.2,a,es8.2,a,es9.2,a)') ' error divergence = ', & - err_div/divTol, ' (',err_div,' / m, tol = ',divTol,')' - write(6,'(a,f12.2,a,es8.2,a,es9.2,a)') ' error stress BC = ', & - err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')' - write(6,'(/,a)') ' ===========================================================================' - flush(6) + write(6,'(1/,a)') ' ... reporting .............................................................' + write(6,'(1/,a,f12.2,a,es8.2,a,es9.2,a)') ' error divergence = ', & + err_div/divTol, ' (',err_div,' / m, tol = ',divTol,')' + write(6,'(a,f12.2,a,es8.2,a,es9.2,a)') ' error stress BC = ', & + err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')' + write(6,'(/,a)') ' ===========================================================================' + flush(6) end subroutine converged @@ -501,136 +479,136 @@ end subroutine converged !-------------------------------------------------------------------------------------------------- !> @brief forms the residual vector !-------------------------------------------------------------------------------------------------- -subroutine formResidual(da_local,x_local,f_local,dummy,ierr) - use numerics, only: & - itmax, & - itmin - use numerics, only: & - worldrank - use mesh, only: & - grid - use math, only: & - math_rotate_backward33, & - math_mul3333xx33 - use debug, only: & - debug_level, & - debug_spectral, & - debug_spectralRotation - use spectral_utilities, only: & - utilities_constitutiveResponse - use IO, only: & - IO_intOut - use FEsolving, only: & - terminallyIll - use homogenization, only: & - materialpoint_dPdF +subroutine formResidual(da_local,x_local, & + f_local,dummy,ierr) + use numerics, only: & + itmax, & + itmin + use numerics, only: & + worldrank + use mesh, only: & + grid + use math, only: & + math_rotate_backward33, & + math_mul3333xx33 + use debug, only: & + debug_level, & + debug_spectral, & + debug_spectralRotation + use spectral_utilities, only: & + utilities_constitutiveResponse + use IO, only: & + IO_intOut + use FEsolving, only: & + terminallyIll + use homogenization, only: & + materialpoint_dPdF - implicit none - DM :: da_local - Vec :: x_local, f_local - PetscScalar, pointer,dimension(:,:,:,:) :: x_scal, f_scal - PetscScalar, dimension(8,3) :: x_elem, f_elem - PetscInt :: i, ii, j, jj, k, kk, ctr, ele - real(pReal), dimension(3,3) :: & - deltaF_aim - PetscInt :: & - PETScIter, & - nfuncs - PetscObject :: dummy - PetscErrorCode :: ierr - real(pReal), dimension(3,3,3,3) :: devNull + DM :: da_local + Vec :: x_local, f_local + PetscScalar, pointer,dimension(:,:,:,:) :: x_scal, f_scal + PetscScalar, dimension(8,3) :: x_elem, f_elem + PetscInt :: i, ii, j, jj, k, kk, ctr, ele + real(pReal), dimension(3,3) :: & + deltaF_aim + PetscInt :: & + PETScIter, & + nfuncs + PetscObject :: dummy + PetscErrorCode :: ierr + real(pReal), dimension(3,3,3,3) :: devNull - call SNESGetNumberFunctionEvals(mech_snes,nfuncs,ierr); CHKERRQ(ierr) - call SNESGetIterationNumber(mech_snes,PETScIter,ierr); CHKERRQ(ierr) + call SNESGetNumberFunctionEvals(mech_snes,nfuncs,ierr); CHKERRQ(ierr) + call SNESGetIterationNumber(mech_snes,PETScIter,ierr); CHKERRQ(ierr) - if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1_pInt ! new increment + if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1 ! new increment !-------------------------------------------------------------------------------------------------- ! begin of new iteration - newIteration: if (totalIter <= PETScIter) then - totalIter = totalIter + 1_pInt - write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') & - trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter+1, '≤', itmax - if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) & - write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & - ' deformation gradient aim (lab) =', transpose(math_rotate_backward33(F_aim,params%rotation_BC)) - write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & - ' deformation gradient aim =', transpose(F_aim) - flush(6) - endif newIteration + newIteration: if (totalIter <= PETScIter) then + totalIter = totalIter + 1 + write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') & + trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter+1, '≤', itmax + if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) & + write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & + ' deformation gradient aim (lab) =', transpose(math_rotate_backward33(F_aim,params%rotation_BC)) + write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & + ' deformation gradient aim =', transpose(F_aim) + flush(6) + endif newIteration !-------------------------------------------------------------------------------------------------- ! get deformation gradient - call DMDAVecGetArrayF90(da_local,x_local,x_scal,ierr);CHKERRQ(ierr) - do k = zstart, zend; do j = ystart, yend; do i = xstart, xend - ctr = 0 - do kk = 0, 1; do jj = 0, 1; do ii = 0, 1 - ctr = ctr + 1 - x_elem(ctr,1:3) = x_scal(0:2,i+ii,j+jj,k+kk) - enddo; enddo; enddo - ii = i-xstart+1; jj = j-ystart+1; kk = k-zstart+1 - F(1:3,1:3,ii,jj,kk) = math_rotate_backward33(F_aim,params%rotation_BC) + transpose(matmul(BMat,x_elem)) - enddo; enddo; enddo - call DMDAVecRestoreArrayF90(da_local,x_local,x_scal,ierr);CHKERRQ(ierr) + call DMDAVecGetArrayF90(da_local,x_local,x_scal,ierr);CHKERRQ(ierr) + do k = zstart, zend; do j = ystart, yend; do i = xstart, xend + ctr = 0 + do kk = 0, 1; do jj = 0, 1; do ii = 0, 1 + ctr = ctr + 1 + x_elem(ctr,1:3) = x_scal(0:2,i+ii,j+jj,k+kk) + enddo; enddo; enddo + ii = i-xstart+1; jj = j-ystart+1; kk = k-zstart+1 + F(1:3,1:3,ii,jj,kk) = math_rotate_backward33(F_aim,params%rotation_BC) + transpose(matmul(BMat,x_elem)) + enddo; enddo; enddo + call DMDAVecRestoreArrayF90(da_local,x_local,x_scal,ierr);CHKERRQ(ierr) !-------------------------------------------------------------------------------------------------- ! evaluate constitutive response - call Utilities_constitutiveResponse(P_current,& - P_av,C_volAvg,devNull, & - F,params%timeinc,params%rotation_BC) - call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr) + call Utilities_constitutiveResponse(P_current,& + P_av,C_volAvg,devNull, & + F,params%timeinc,params%rotation_BC) + call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr) !-------------------------------------------------------------------------------------------------- ! stress BC handling - F_aim_lastIter = F_aim - deltaF_aim = math_mul3333xx33(S, P_av - params%stress_BC) - F_aim = F_aim - deltaF_aim - err_BC = maxval(abs(params%stress_mask * (P_av - params%stress_BC))) ! mask = 0.0 when no stress bc + F_aim_lastIter = F_aim + deltaF_aim = math_mul3333xx33(S, P_av - params%stress_BC) + F_aim = F_aim - deltaF_aim + err_BC = maxval(abs(params%stress_mask * (P_av - params%stress_BC))) ! mask = 0.0 when no stress bc !-------------------------------------------------------------------------------------------------- ! constructing residual - call VecSet(f_local,0.0_pReal,ierr);CHKERRQ(ierr) - call DMDAVecGetArrayF90(da_local,f_local,f_scal,ierr);CHKERRQ(ierr) - call DMDAVecGetArrayF90(da_local,x_local,x_scal,ierr);CHKERRQ(ierr) - ele = 0 - do k = zstart, zend; do j = ystart, yend; do i = xstart, xend - ctr = 0 - do kk = 0, 1; do jj = 0, 1; do ii = 0, 1 - ctr = ctr + 1 - x_elem(ctr,1:3) = x_scal(0:2,i+ii,j+jj,k+kk) - enddo; enddo; enddo - ii = i-xstart+1; jj = j-ystart+1; kk = k-zstart+1 - ele = ele + 1 - f_elem = matmul(transpose(BMat),transpose(P_current(1:3,1:3,ii,jj,kk)))*detJ + & - matmul(HGMat,x_elem)*(materialpoint_dPdF(1,1,1,1,1,ele) + & - materialpoint_dPdF(2,2,2,2,1,ele) + & - materialpoint_dPdF(3,3,3,3,1,ele))/3.0_pReal - ctr = 0 - do kk = 0, 1; do jj = 0, 1; do ii = 0, 1 - ctr = ctr + 1 - f_scal(0:2,i+ii,j+jj,k+kk) = f_scal(0:2,i+ii,j+jj,k+kk) + f_elem(ctr,1:3) - enddo; enddo; enddo - enddo; enddo; enddo - call DMDAVecRestoreArrayF90(da_local,x_local,x_scal,ierr);CHKERRQ(ierr) - call DMDAVecRestoreArrayF90(da_local,f_local,f_scal,ierr);CHKERRQ(ierr) + call VecSet(f_local,0.0_pReal,ierr);CHKERRQ(ierr) + call DMDAVecGetArrayF90(da_local,f_local,f_scal,ierr);CHKERRQ(ierr) + call DMDAVecGetArrayF90(da_local,x_local,x_scal,ierr);CHKERRQ(ierr) + ele = 0 + do k = zstart, zend; do j = ystart, yend; do i = xstart, xend + ctr = 0 + do kk = 0, 1; do jj = 0, 1; do ii = 0, 1 + ctr = ctr + 1 + x_elem(ctr,1:3) = x_scal(0:2,i+ii,j+jj,k+kk) + enddo; enddo; enddo + ii = i-xstart+1; jj = j-ystart+1; kk = k-zstart+1 + ele = ele + 1 + f_elem = matmul(transpose(BMat),transpose(P_current(1:3,1:3,ii,jj,kk)))*detJ + & + matmul(HGMat,x_elem)*(materialpoint_dPdF(1,1,1,1,1,ele) + & + materialpoint_dPdF(2,2,2,2,1,ele) + & + materialpoint_dPdF(3,3,3,3,1,ele))/3.0_pReal + ctr = 0 + do kk = 0, 1; do jj = 0, 1; do ii = 0, 1 + ctr = ctr + 1 + f_scal(0:2,i+ii,j+jj,k+kk) = f_scal(0:2,i+ii,j+jj,k+kk) + f_elem(ctr,1:3) + enddo; enddo; enddo + enddo; enddo; enddo + call DMDAVecRestoreArrayF90(da_local,x_local,x_scal,ierr);CHKERRQ(ierr) + call DMDAVecRestoreArrayF90(da_local,f_local,f_scal,ierr);CHKERRQ(ierr) !-------------------------------------------------------------------------------------------------- ! applying boundary conditions - call DMDAVecGetArrayF90(da_local,f_local,f_scal,ierr);CHKERRQ(ierr) - if (zstart == 0) then - f_scal(0:2,xstart,ystart,zstart) = 0.0 - f_scal(0:2,xend+1,ystart,zstart) = 0.0 - f_scal(0:2,xstart,yend+1,zstart) = 0.0 - f_scal(0:2,xend+1,yend+1,zstart) = 0.0 - endif - if (zend + 1 == grid(3)) then - f_scal(0:2,xstart,ystart,zend+1) = 0.0 - f_scal(0:2,xend+1,ystart,zend+1) = 0.0 - f_scal(0:2,xstart,yend+1,zend+1) = 0.0 - f_scal(0:2,xend+1,yend+1,zend+1) = 0.0 - endif - call DMDAVecRestoreArrayF90(da_local,f_local,f_scal,ierr);CHKERRQ(ierr) + call DMDAVecGetArrayF90(da_local,f_local,f_scal,ierr);CHKERRQ(ierr) + if (zstart == 0) then + f_scal(0:2,xstart,ystart,zstart) = 0.0 + f_scal(0:2,xend+1,ystart,zstart) = 0.0 + f_scal(0:2,xstart,yend+1,zstart) = 0.0 + f_scal(0:2,xend+1,yend+1,zstart) = 0.0 + endif + if (zend + 1 == grid(3)) then + f_scal(0:2,xstart,ystart,zend+1) = 0.0 + f_scal(0:2,xend+1,ystart,zend+1) = 0.0 + f_scal(0:2,xstart,yend+1,zend+1) = 0.0 + f_scal(0:2,xend+1,yend+1,zend+1) = 0.0 + endif + call DMDAVecRestoreArrayF90(da_local,f_local,f_scal,ierr);CHKERRQ(ierr) end subroutine formResidual @@ -639,97 +617,96 @@ end subroutine formResidual !> @brief forms the FEM stiffness matrix !-------------------------------------------------------------------------------------------------- subroutine formJacobian(da_local,x_local,Jac_pre,Jac,dummy,ierr) - use mesh, only: & - mesh_ipCoordinates - use homogenization, only: & - materialpoint_dPdF - - implicit none - - DM :: da_local - Vec :: x_local, coordinates - Mat :: Jac_pre, Jac - MatStencil,dimension(4,24) :: row, col - PetscScalar,pointer,dimension(:,:,:,:) :: x_scal - PetscScalar,dimension(24,24) :: K_ele - PetscScalar,dimension(9,24) :: BMatFull - PetscInt :: i, ii, j, jj, k, kk, ctr, ele - PetscInt,dimension(3) :: rows - PetscScalar :: diag - PetscObject :: dummy - MatNullSpace :: matnull - PetscErrorCode :: ierr + use mesh, only: & + mesh_ipCoordinates + use homogenization, only: & + materialpoint_dPdF - BMatFull = 0.0 - BMatFull(1:3,1 :8 ) = BMat - BMatFull(4:6,9 :16) = BMat - BMatFull(7:9,17:24) = BMat - call MatSetOption(Jac,MAT_KEEP_NONZERO_PATTERN,PETSC_TRUE,ierr); CHKERRQ(ierr) - call MatSetOption(Jac,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE,ierr); CHKERRQ(ierr) - call MatZeroEntries(Jac,ierr); CHKERRQ(ierr) - ele = 0 - do k = zstart, zend; do j = ystart, yend; do i = xstart, xend - ctr = 0 - do kk = 0, 1; do jj = 0, 1; do ii = 0, 1 - ctr = ctr + 1 - col(MatStencil_i,ctr ) = i+ii - col(MatStencil_j,ctr ) = j+jj - col(MatStencil_k,ctr ) = k+kk - col(MatStencil_c,ctr ) = 0 - col(MatStencil_i,ctr+8 ) = i+ii - col(MatStencil_j,ctr+8 ) = j+jj - col(MatStencil_k,ctr+8 ) = k+kk - col(MatStencil_c,ctr+8 ) = 1 - col(MatStencil_i,ctr+16) = i+ii - col(MatStencil_j,ctr+16) = j+jj - col(MatStencil_k,ctr+16) = k+kk - col(MatStencil_c,ctr+16) = 2 - enddo; enddo; enddo - row = col - ele = ele + 1 - K_ele = 0.0 - K_ele(1 :8 ,1 :8 ) = HGMat*(materialpoint_dPdF(1,1,1,1,1,ele) + & - materialpoint_dPdF(2,2,2,2,1,ele) + & - materialpoint_dPdF(3,3,3,3,1,ele))/3.0_pReal - K_ele(9 :16,9 :16) = HGMat*(materialpoint_dPdF(1,1,1,1,1,ele) + & - materialpoint_dPdF(2,2,2,2,1,ele) + & - materialpoint_dPdF(3,3,3,3,1,ele))/3.0_pReal - K_ele(17:24,17:24) = HGMat*(materialpoint_dPdF(1,1,1,1,1,ele) + & - materialpoint_dPdF(2,2,2,2,1,ele) + & - materialpoint_dPdF(3,3,3,3,1,ele))/3.0_pReal - K_ele = K_ele + & - matmul(transpose(BMatFull), & - matmul(reshape(reshape(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,ele), & - shape=[3,3,3,3], order=[2,1,4,3]),shape=[9,9]),BMatFull))*detJ - call MatSetValuesStencil(Jac,24,row,24,col,K_ele,ADD_VALUES,ierr) - CHKERRQ(ierr) - enddo; enddo; enddo - call MatAssemblyBegin(Jac,MAT_FINAL_ASSEMBLY,ierr); CHKERRQ(ierr) - call MatAssemblyEnd(Jac,MAT_FINAL_ASSEMBLY,ierr); CHKERRQ(ierr) - call MatAssemblyBegin(Jac_pre,MAT_FINAL_ASSEMBLY,ierr); CHKERRQ(ierr) - call MatAssemblyEnd(Jac_pre,MAT_FINAL_ASSEMBLY,ierr); CHKERRQ(ierr) + + DM :: da_local + Vec :: x_local, coordinates + Mat :: Jac_pre, Jac + MatStencil,dimension(4,24) :: row, col + PetscScalar,pointer,dimension(:,:,:,:) :: x_scal + PetscScalar,dimension(24,24) :: K_ele + PetscScalar,dimension(9,24) :: BMatFull + PetscInt :: i, ii, j, jj, k, kk, ctr, ele + PetscInt,dimension(3) :: rows + PetscScalar :: diag + PetscObject :: dummy + MatNullSpace :: matnull + PetscErrorCode :: ierr + + BMatFull = 0.0 + BMatFull(1:3,1 :8 ) = BMat + BMatFull(4:6,9 :16) = BMat + BMatFull(7:9,17:24) = BMat + call MatSetOption(Jac,MAT_KEEP_NONZERO_PATTERN,PETSC_TRUE,ierr); CHKERRQ(ierr) + call MatSetOption(Jac,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE,ierr); CHKERRQ(ierr) + call MatZeroEntries(Jac,ierr); CHKERRQ(ierr) + ele = 0 + do k = zstart, zend; do j = ystart, yend; do i = xstart, xend + ctr = 0 + do kk = 0, 1; do jj = 0, 1; do ii = 0, 1 + ctr = ctr + 1 + col(MatStencil_i,ctr ) = i+ii + col(MatStencil_j,ctr ) = j+jj + col(MatStencil_k,ctr ) = k+kk + col(MatStencil_c,ctr ) = 0 + col(MatStencil_i,ctr+8 ) = i+ii + col(MatStencil_j,ctr+8 ) = j+jj + col(MatStencil_k,ctr+8 ) = k+kk + col(MatStencil_c,ctr+8 ) = 1 + col(MatStencil_i,ctr+16) = i+ii + col(MatStencil_j,ctr+16) = j+jj + col(MatStencil_k,ctr+16) = k+kk + col(MatStencil_c,ctr+16) = 2 + enddo; enddo; enddo + row = col + ele = ele + 1 + K_ele = 0.0 + K_ele(1 :8 ,1 :8 ) = HGMat*(materialpoint_dPdF(1,1,1,1,1,ele) + & + materialpoint_dPdF(2,2,2,2,1,ele) + & + materialpoint_dPdF(3,3,3,3,1,ele))/3.0_pReal + K_ele(9 :16,9 :16) = HGMat*(materialpoint_dPdF(1,1,1,1,1,ele) + & + materialpoint_dPdF(2,2,2,2,1,ele) + & + materialpoint_dPdF(3,3,3,3,1,ele))/3.0_pReal + K_ele(17:24,17:24) = HGMat*(materialpoint_dPdF(1,1,1,1,1,ele) + & + materialpoint_dPdF(2,2,2,2,1,ele) + & + materialpoint_dPdF(3,3,3,3,1,ele))/3.0_pReal + K_ele = K_ele + & + matmul(transpose(BMatFull), & + matmul(reshape(reshape(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,ele), & + shape=[3,3,3,3], order=[2,1,4,3]),shape=[9,9]),BMatFull))*detJ + call MatSetValuesStencil(Jac,24,row,24,col,K_ele,ADD_VALUES,ierr) + CHKERRQ(ierr) + enddo; enddo; enddo + call MatAssemblyBegin(Jac,MAT_FINAL_ASSEMBLY,ierr); CHKERRQ(ierr) + call MatAssemblyEnd(Jac,MAT_FINAL_ASSEMBLY,ierr); CHKERRQ(ierr) + call MatAssemblyBegin(Jac_pre,MAT_FINAL_ASSEMBLY,ierr); CHKERRQ(ierr) + call MatAssemblyEnd(Jac_pre,MAT_FINAL_ASSEMBLY,ierr); CHKERRQ(ierr) !-------------------------------------------------------------------------------------------------- ! applying boundary conditions - rows = [0, 1, 2] - diag = (C_volAvg(1,1,1,1)/delta(1)**2.0_pReal + & - C_volAvg(2,2,2,2)/delta(2)**2.0_pReal + & - C_volAvg(3,3,3,3)/delta(3)**2.0_pReal)*detJ - call MatZeroRowsColumns(Jac,size(rows),rows,diag,PETSC_NULL_VEC,PETSC_NULL_VEC,ierr) - CHKERRQ(ierr) - call DMGetGlobalVector(da_local,coordinates,ierr);CHKERRQ(ierr) - call DMDAVecGetArrayF90(da_local,coordinates,x_scal,ierr);CHKERRQ(ierr) - ele = 0 - do k = zstart, zend; do j = ystart, yend; do i = xstart, xend - ele = ele + 1 - x_scal(0:2,i,j,k) = mesh_ipCoordinates(1:3,1,ele) - enddo; enddo; enddo - call DMDAVecRestoreArrayF90(da_local,coordinates,x_scal,ierr);CHKERRQ(ierr) ! initialize to undeformed coordinates (ToDo: use ip coordinates) - call MatNullSpaceCreateRigidBody(coordinates,matnull,ierr);CHKERRQ(ierr) ! get rigid body deformation modes - call DMRestoreGlobalVector(da_local,coordinates,ierr);CHKERRQ(ierr) - call MatSetNullSpace(Jac,matnull,ierr); CHKERRQ(ierr) - call MatSetNearNullSpace(Jac,matnull,ierr); CHKERRQ(ierr) - call MatNullSpaceDestroy(matnull,ierr); CHKERRQ(ierr) + rows = [0, 1, 2] + diag = (C_volAvg(1,1,1,1)/delta(1)**2.0_pReal + & + C_volAvg(2,2,2,2)/delta(2)**2.0_pReal + & + C_volAvg(3,3,3,3)/delta(3)**2.0_pReal)*detJ + call MatZeroRowsColumns(Jac,size(rows),rows,diag,PETSC_NULL_VEC,PETSC_NULL_VEC,ierr) + CHKERRQ(ierr) + call DMGetGlobalVector(da_local,coordinates,ierr);CHKERRQ(ierr) + call DMDAVecGetArrayF90(da_local,coordinates,x_scal,ierr);CHKERRQ(ierr) + ele = 0 + do k = zstart, zend; do j = ystart, yend; do i = xstart, xend + ele = ele + 1 + x_scal(0:2,i,j,k) = mesh_ipCoordinates(1:3,1,ele) + enddo; enddo; enddo + call DMDAVecRestoreArrayF90(da_local,coordinates,x_scal,ierr);CHKERRQ(ierr) ! initialize to undeformed coordinates (ToDo: use ip coordinates) + call MatNullSpaceCreateRigidBody(coordinates,matnull,ierr);CHKERRQ(ierr) ! get rigid body deformation modes + call DMRestoreGlobalVector(da_local,coordinates,ierr);CHKERRQ(ierr) + call MatSetNullSpace(Jac,matnull,ierr); CHKERRQ(ierr) + call MatSetNearNullSpace(Jac,matnull,ierr); CHKERRQ(ierr) + call MatNullSpaceDestroy(matnull,ierr); CHKERRQ(ierr) end subroutine formJacobian diff --git a/src/grid/grid_mech_spectral_basic.f90 b/src/grid/grid_mech_spectral_basic.f90 index 99839e50f..2daebefbd 100644 --- a/src/grid/grid_mech_spectral_basic.f90 +++ b/src/grid/grid_mech_spectral_basic.f90 @@ -7,6 +7,8 @@ module grid_mech_spectral_basic #include #include + use DAMASK_interface + use HDF5_utilities use PETScdmda use PETScsnes use prec, only: & @@ -25,8 +27,7 @@ module grid_mech_spectral_basic type(tSolutionParams), private :: params type, private :: tNumerics - logical :: & - update_gamma !< update gamma operator with current stiffness + logical :: update_gamma !< update gamma operator with current stiffness end type tNumerics type(tNumerics) :: num ! numerics parameters. Better name? @@ -40,8 +41,8 @@ module grid_mech_spectral_basic !-------------------------------------------------------------------------------------------------- ! common pointwise data real(pReal), private, dimension(:,:,:,:,:), allocatable :: & - F_lastInc, & - Fdot + F_lastInc, & !< field of previous compatible deformation gradients + Fdot !< field of assumed rate of compatible deformation gradient !-------------------------------------------------------------------------------------------------- ! stress, stiffness and compliance average etc. @@ -99,23 +100,22 @@ subroutine grid_mech_spectral_basic_init use spectral_utilities, only: & utilities_constitutiveResponse, & utilities_updateGamma, & - utilities_updateIPcoords, & - wgt + utilities_updateIPcoords use mesh, only: & grid, & grid3 use math, only: & math_invSym3333 - implicit none real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: P real(pReal), dimension(3,3) :: & temp33_Real = 0.0_pReal PetscErrorCode :: ierr PetscScalar, pointer, dimension(:,:,:,:) :: & - F ! pointer to solution data + F ! pointer to solution data PetscInt, dimension(worldsize) :: localK + integer(HID_T) :: fileHandle integer :: fileUnit character(len=1024) :: rankStr @@ -163,7 +163,7 @@ subroutine grid_mech_spectral_basic_init call DMcreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr) ! global solution vector (grid x 9, i.e. every def grad tensor) call DMDASNESsetFunctionLocal(da,INSERT_VALUES,formResidual,PETSC_NULL_SNES,ierr) ! residual vector of same shape as solution vector CHKERRQ(ierr) - call SNESsetConvergenceTest(snes,converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,ierr)! specify custom convergence check function "converged" + call SNESsetConvergenceTest(snes,converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,ierr) ! specify custom convergence check function "converged" CHKERRQ(ierr) call SNESsetFromOptions(snes,ierr); CHKERRQ(ierr) ! pull it all together with additional CLI arguments @@ -174,19 +174,14 @@ subroutine grid_mech_spectral_basic_init restart: if (restartInc > 0) then write(6,'(/,a,'//IO_intOut(restartInc)//',a)') ' reading values of increment ', restartInc, ' from file' - fileUnit = IO_open_jobFile_binary('F_aim') - read(fileUnit) F_aim; close(fileUnit) - fileUnit = IO_open_jobFile_binary('F_aim_lastInc') - read(fileUnit) F_aim_lastInc; close(fileUnit) - fileUnit = IO_open_jobFile_binary('F_aimDot') - read(fileUnit) F_aimDot; close(fileUnit) - write(rankStr,'(a1,i0)')'_',worldrank + fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5') - fileUnit = IO_open_jobFile_binary('F'//trim(rankStr)) - read(fileUnit) F; close (fileUnit) - fileUnit = IO_open_jobFile_binary('F_lastInc'//trim(rankStr)) - read(fileUnit) F_lastInc; close (fileUnit) + call HDF5_read(fileHandle,F_aim, 'F_aim') + call HDF5_read(fileHandle,F_aim_lastInc,'F_aim_lastInc') + call HDF5_read(fileHandle,F_aimDot, 'F_aimDot') + call HDF5_read(fileHandle,F, 'F') + call HDF5_read(fileHandle,F_lastInc, 'F_lastInc') elseif (restartInc == 0) then restart F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity @@ -203,15 +198,15 @@ subroutine grid_mech_spectral_basic_init restartRead: if (restartInc > 0) then write(6,'(/,a,'//IO_intOut(restartInc)//',a)') 'reading more values of increment ', restartInc, ' from file' - fileUnit = IO_open_jobFile_binary('C_volAvg') - read(fileUnit) C_volAvg; close(fileUnit) - fileUnit = IO_open_jobFile_binary('C_volAvgLastInv') - read(fileUnit) C_volAvgLastInc; close(fileUnit) + call HDF5_read(fileHandle,C_volAvg, 'C_volAvg') + call HDF5_read(fileHandle,C_volAvgLastInc,'C_volAvgLastInc') + call HDF5_closeFile(fileHandle) + fileUnit = IO_open_jobFile_binary('C_ref') read(fileUnit) C_minMaxAvg; close(fileUnit) endif restartRead - call Utilities_updateGamma(C_minMaxAvg,.true.) + call utilities_updateGamma(C_minMaxAvg,.true.) end subroutine grid_mech_spectral_basic_init @@ -228,8 +223,6 @@ function grid_mech_spectral_basic_solution(incInfoIn,timeinc,timeinc_old,stress_ restartWrite, & terminallyIll - implicit none - !-------------------------------------------------------------------------------------------------- ! input data for solution character(len=*), intent(in) :: & @@ -251,8 +244,8 @@ function grid_mech_spectral_basic_solution(incInfoIn,timeinc,timeinc_old,stress_ !-------------------------------------------------------------------------------------------------- ! update stiffness (and gamma operator) - S = Utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg) - if (num%update_gamma) call Utilities_updateGamma(C_minMaxAvg,restartWrite) + S = utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg) + if (num%update_gamma) call utilities_updateGamma(C_minMaxAvg,restartWrite) !-------------------------------------------------------------------------------------------------- ! set module wide available data @@ -306,7 +299,6 @@ subroutine grid_mech_spectral_basic_forward(guess,timeinc,timeinc_old,loadCaseTi use FEsolving, only: & restartWrite - implicit none logical, intent(in) :: & guess real(pReal), intent(in) :: & @@ -321,7 +313,7 @@ subroutine grid_mech_spectral_basic_forward(guess,timeinc,timeinc_old,loadCaseTi PetscErrorCode :: ierr PetscScalar, dimension(:,:,:,:), pointer :: F - integer :: fileUnit + integer(HID_T) :: fileHandle character(len=32) :: rankStr call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) @@ -331,29 +323,24 @@ subroutine grid_mech_spectral_basic_forward(guess,timeinc,timeinc_old,loadCaseTi C_minMaxAvg = C_minMaxAvgLastInc ! QUESTION: where is this required? else !-------------------------------------------------------------------------------------------------- - ! restart information for spectral solver - if (restartWrite) then ! QUESTION: where is this logical properly set? - write(6,'(/,a)') ' writing converged results for restart' - flush(6) - - if (worldrank == 0) then - fileUnit = IO_open_jobFile_binary('C_volAvg','w') - write(fileUnit) C_volAvg; close(fileUnit) - fileUnit = IO_open_jobFile_binary('C_volAvgLastInv','w') - write(fileUnit) C_volAvgLastInc; close(fileUnit) - fileUnit = IO_open_jobFile_binary('F_aim','w') - write(fileUnit) F_aim; close(fileUnit) - fileUnit = IO_open_jobFile_binary('F_aim_lastInc','w') - write(fileUnit) F_aim_lastInc; close(fileUnit) - fileUnit = IO_open_jobFile_binary('F_aimDot','w') - write(fileUnit) F_aimDot; close(fileUnit) - endif - + ! restart information for spectral solver + if (restartWrite) then + write(6,'(/,a)') ' writing converged results for restart';flush(6) + write(rankStr,'(a1,i0)')'_',worldrank - fileUnit = IO_open_jobFile_binary('F'//trim(rankStr),'w') - write(fileUnit) F; close (fileUnit) - fileUnit = IO_open_jobFile_binary('F_lastInc'//trim(rankStr),'w') - write(fileUnit) F_lastInc; close (fileUnit) + fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5','w') + + call HDF5_write(fileHandle,F_aim, 'F_aim') + call HDF5_write(fileHandle,F_aim_lastInc,'F_aim_lastInc') + call HDF5_write(fileHandle,F_aimDot, 'F_aimDot') + call HDF5_write(fileHandle,F, 'F') + call HDF5_write(fileHandle,F_lastInc, 'F_lastInc') + + call HDF5_write(fileHandle,C_volAvg, 'C_volAvg') + call HDF5_write(fileHandle,C_volAvgLastInc,'C_volAvgLastInc') + call HDF5_write(fileHandle,C_minMaxAvg, 'C_minMaxAvg') + + call HDF5_closeFile(fileHandle) endif call CPFEM_age ! age state and kinematics @@ -399,7 +386,7 @@ end subroutine grid_mech_spectral_basic_forward !-------------------------------------------------------------------------------------------------- !> @brief convergence check !-------------------------------------------------------------------------------------------------- -subroutine converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr) +subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dummy,ierr) use numerics, only: & itmax, & itmin, & @@ -410,13 +397,12 @@ subroutine converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr) use FEsolving, only: & terminallyIll - implicit none SNES :: snes_local - PetscInt :: PETScIter - PetscReal :: & - xnorm, & ! not used - snorm, & ! not used - fnorm ! not used + PetscInt, intent(in) :: PETScIter + PetscReal, intent(in) :: & + devNull1, & + devNull2, & + devNull3 SNESConvergedReason :: reason PetscObject :: dummy PetscErrorCode :: ierr @@ -452,7 +438,7 @@ end subroutine converged !-------------------------------------------------------------------------------------------------- -!> @brief forms the basic residual vector +!> @brief forms the residual vector !-------------------------------------------------------------------------------------------------- subroutine formResidual(in, F, & residuum, dummy, ierr) @@ -481,7 +467,6 @@ subroutine formResidual(in, F, & use FEsolving, only: & terminallyIll - implicit none DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: in !< DMDA info (needs to be named "in" for macros like XRANGE to work) PetscScalar, dimension(3,3,XG_RANGE,YG_RANGE,ZG_RANGE), & intent(in) :: F !< deformation gradient field @@ -515,7 +500,7 @@ subroutine formResidual(in, F, & !-------------------------------------------------------------------------------------------------- ! evaluate constitutive response - call Utilities_constitutiveResponse(residuum, & ! "residuum" gets field of first PK stress (to save memory) + call utilities_constitutiveResponse(residuum, & ! "residuum" gets field of first PK stress (to save memory) P_av,C_volAvg,C_minMaxAvg, & F,params%timeinc,params%rotation_BC) call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr) diff --git a/src/grid/grid_mech_spectral_polarisation.f90 b/src/grid/grid_mech_spectral_polarisation.f90 index aff4913b1..3bd30a360 100644 --- a/src/grid/grid_mech_spectral_polarisation.f90 +++ b/src/grid/grid_mech_spectral_polarisation.f90 @@ -7,78 +7,79 @@ module grid_mech_spectral_polarisation #include #include - use PETScdmda - use PETScsnes - use prec, only: & - pReal - use math, only: & - math_I3 - use spectral_utilities, only: & - tSolutionState, & - tSolutionParams - - implicit none - private + use DAMASK_interface + use HDF5_utilities + use PETScdmda + use PETScsnes + use prec, only: & + pReal + use math, only: & + math_I3 + use spectral_utilities, only: & + tSolutionState, & + tSolutionParams + + implicit none + private !-------------------------------------------------------------------------------------------------- ! derived types type(tSolutionParams), private :: params type, private :: tNumerics - logical :: & - update_gamma !< update gamma operator with current stiffness + logical :: update_gamma !< update gamma operator with current stiffness end type tNumerics type(tNumerics) :: num ! numerics parameters. Better name? !-------------------------------------------------------------------------------------------------- ! PETSc data - DM, private :: da - SNES, private :: snes - Vec, private :: solution_vec + DM, private :: da + SNES, private :: snes + Vec, private :: solution_vec !-------------------------------------------------------------------------------------------------- ! common pointwise data - real(pReal), private, dimension(:,:,:,:,:), allocatable :: & - F_lastInc, & !< field of previous compatible deformation gradients - F_tau_lastInc, & !< field of previous incompatible deformation gradient - Fdot, & !< field of assumed rate of compatible deformation gradient - F_tauDot !< field of assumed rate of incopatible deformation gradient + real(pReal), private, dimension(:,:,:,:,:), allocatable :: & + F_lastInc, & !< field of previous compatible deformation gradients + F_tau_lastInc, & !< field of previous incompatible deformation gradient + Fdot, & !< field of assumed rate of compatible deformation gradient + F_tauDot !< field of assumed rate of incopatible deformation gradient !-------------------------------------------------------------------------------------------------- ! stress, stiffness and compliance average etc. - real(pReal), private, dimension(3,3) :: & - F_aimDot = 0.0_pReal, & !< assumed rate of average deformation gradient - F_aim = math_I3, & !< current prescribed deformation gradient - F_aim_lastInc = math_I3, & !< previous average deformation gradient - F_av = 0.0_pReal, & !< average incompatible def grad field - P_av = 0.0_pReal !< average 1st Piola--Kirchhoff stress - - character(len=1024), private :: incInfo !< time and increment information - real(pReal), private, dimension(3,3,3,3) :: & - C_volAvg = 0.0_pReal, & !< current volume average stiffness - C_volAvgLastInc = 0.0_pReal, & !< previous volume average stiffness - C_minMaxAvg = 0.0_pReal, & !< current (min+max)/2 stiffness - C_minMaxAvgLastInc = 0.0_pReal, & !< previous (min+max)/2 stiffness - S = 0.0_pReal, & !< current compliance (filled up with zeros) - C_scale = 0.0_pReal, & - S_scale = 0.0_pReal - - real(pReal), private :: & - err_BC, & !< deviation from stress BC - err_curl, & !< RMS of curl of F - err_div !< RMS of div of P + real(pReal), private, dimension(3,3) :: & + F_aimDot = 0.0_pReal, & !< assumed rate of average deformation gradient + F_aim = math_I3, & !< current prescribed deformation gradient + F_aim_lastInc = math_I3, & !< previous average deformation gradient + F_av = 0.0_pReal, & !< average incompatible def grad field + P_av = 0.0_pReal !< average 1st Piola--Kirchhoff stress - integer, private :: & - totalIter = 0 !< total iteration in current increment + character(len=1024), private :: incInfo !< time and increment information + real(pReal), private, dimension(3,3,3,3) :: & + C_volAvg = 0.0_pReal, & !< current volume average stiffness + C_volAvgLastInc = 0.0_pReal, & !< previous volume average stiffness + C_minMaxAvg = 0.0_pReal, & !< current (min+max)/2 stiffness + C_minMaxAvgLastInc = 0.0_pReal, & !< previous (min+max)/2 stiffness + S = 0.0_pReal, & !< current compliance (filled up with zeros) + C_scale = 0.0_pReal, & + S_scale = 0.0_pReal - public :: & - grid_mech_spectral_polarisation_init, & - grid_mech_spectral_polarisation_solution, & - grid_mech_spectral_polarisation_forward - private :: & - converged, & - formResidual + real(pReal), private :: & + err_BC, & !< deviation from stress BC + err_curl, & !< RMS of curl of F + err_div !< RMS of div of P + + integer, private :: & + totalIter = 0 !< total iteration in current increment + + public :: & + grid_mech_spectral_polarisation_init, & + grid_mech_spectral_polarisation_solution, & + grid_mech_spectral_polarisation_forward + private :: & + converged, & + formResidual contains @@ -86,149 +87,141 @@ contains !> @brief allocates all necessary fields and fills them with data, potentially from restart info !-------------------------------------------------------------------------------------------------- subroutine grid_mech_spectral_polarisation_init - use IO, only: & - IO_intOut, & - IO_error, & - IO_open_jobFile_binary - use FEsolving, only: & - restartInc - use config, only :& - config_numerics - use numerics, only: & - worldrank, & - worldsize, & - petsc_options - use homogenization, only: & - materialpoint_F0 - use DAMASK_interface, only: & - getSolverJobName - use spectral_utilities, only: & - utilities_constitutiveResponse, & - utilities_updateGamma, & - utilities_updateIPcoords, & - wgt - use mesh, only: & - grid, & - grid3 - use math, only: & - math_invSym3333 - - implicit none - real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: P - real(pReal), dimension(3,3) :: & - temp33_Real = 0.0_pReal - - PetscErrorCode :: ierr - PetscScalar, pointer, dimension(:,:,:,:) :: & - FandF_tau, & ! overall pointer to solution data - F, & ! specific (sub)pointer - F_tau ! specific (sub)pointer - PetscInt, dimension(worldsize) :: localK - integer :: fileUnit - character(len=1024) :: rankStr + use IO, only: & + IO_intOut, & + IO_error, & + IO_open_jobFile_binary + use FEsolving, only: & + restartInc + use config, only :& + config_numerics + use numerics, only: & + worldrank, & + worldsize, & + petsc_options + use homogenization, only: & + materialpoint_F0 + use DAMASK_interface, only: & + getSolverJobName + use spectral_utilities, only: & + utilities_constitutiveResponse, & + utilities_updateGamma, & + utilities_updateIPcoords + use mesh, only: & + grid, & + grid3 + use math, only: & + math_invSym3333 + + real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: P + real(pReal), dimension(3,3) :: & + temp33_Real = 0.0_pReal - write(6,'(/,a)') ' <<<+- grid_mech_spectral_polarisation init -+>>>' - - write(6,'(/,a)') ' Shanthraj et al., International Journal of Plasticity 66:31–45, 2015' - write(6,'(a)') ' https://doi.org/10.1016/j.ijplas.2014.02.006' + PetscErrorCode :: ierr + PetscScalar, pointer, dimension(:,:,:,:) :: & + FandF_tau, & ! overall pointer to solution data + F, & ! specific (sub)pointer + F_tau ! specific (sub)pointer + PetscInt, dimension(worldsize) :: localK + integer(HID_T) :: fileHandle + integer :: fileUnit + character(len=1024) :: rankStr + + write(6,'(/,a)') ' <<<+- grid_mech_spectral_polarisation init -+>>>' - num%update_gamma = config_numerics%getInt('update_gamma',defaultVal=0) > 0 + write(6,'(/,a)') ' Shanthraj et al., International Journal of Plasticity 66:31–45, 2015' + write(6,'(a)') ' https://doi.org/10.1016/j.ijplas.2014.02.006' + + num%update_gamma = config_numerics%getInt('update_gamma',defaultVal=0) > 0 !-------------------------------------------------------------------------------------------------- ! set default and user defined options for PETSc - call PETScOptionsInsertString(PETSC_NULL_OPTIONS,'-mech_snes_type ngmres',ierr) - CHKERRQ(ierr) - call PETScOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_options),ierr) - CHKERRQ(ierr) + call PETScOptionsInsertString(PETSC_NULL_OPTIONS,'-mech_snes_type ngmres',ierr) + CHKERRQ(ierr) + call PETScOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_options),ierr) + CHKERRQ(ierr) !-------------------------------------------------------------------------------------------------- ! allocate global fields - allocate (F_lastInc (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) - allocate (Fdot (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) - allocate (F_tau_lastInc(3,3,grid(1),grid(2),grid3),source = 0.0_pReal) - allocate (F_tauDot (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) + allocate(F_lastInc (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) + allocate(Fdot (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) + allocate(F_tau_lastInc(3,3,grid(1),grid(2),grid3),source = 0.0_pReal) + allocate(F_tauDot (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) !-------------------------------------------------------------------------------------------------- ! initialize solver specific parts of PETSc - call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr) - call SNESSetOptionsPrefix(snes,'mech_',ierr);CHKERRQ(ierr) - localK = 0 - localK(worldrank+1) = grid3 - call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,ierr) - call DMDACreate3d(PETSC_COMM_WORLD, & - DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & ! cut off stencil at boundary - DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point - grid(1),grid(2),grid(3), & ! global grid - 1 , 1, worldsize, & - 18, 0, & ! #dof (F tensor), ghost boundary width (domain overlap) - [grid(1)],[grid(2)],localK, & ! local grid - da,ierr) ! handle, error - CHKERRQ(ierr) - call SNESSetDM(snes,da,ierr); CHKERRQ(ierr) ! connect snes to da - call DMsetFromOptions(da,ierr); CHKERRQ(ierr) - call DMsetUp(da,ierr); CHKERRQ(ierr) - call DMcreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr) ! global solution vector (grid x 18, i.e. every def grad tensor) - call DMDASNESsetFunctionLocal(da,INSERT_VALUES,formResidual,PETSC_NULL_SNES,ierr) ! residual vector of same shape as solution vector - CHKERRQ(ierr) - call SNESsetConvergenceTest(snes,converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,ierr) ! specify custom convergence check function "converged" - CHKERRQ(ierr) - call SNESsetFromOptions(snes,ierr); CHKERRQ(ierr) ! pull it all together with additional CLI arguments + call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr) + call SNESSetOptionsPrefix(snes,'mech_',ierr);CHKERRQ(ierr) + localK = 0 + localK(worldrank+1) = grid3 + call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,ierr) + call DMDACreate3d(PETSC_COMM_WORLD, & + DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & ! cut off stencil at boundary + DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point + grid(1),grid(2),grid(3), & ! global grid + 1 , 1, worldsize, & + 18, 0, & ! #dof (F tensor), ghost boundary width (domain overlap) + [grid(1)],[grid(2)],localK, & ! local grid + da,ierr) ! handle, error + CHKERRQ(ierr) + call SNESSetDM(snes,da,ierr); CHKERRQ(ierr) ! connect snes to da + call DMsetFromOptions(da,ierr); CHKERRQ(ierr) + call DMsetUp(da,ierr); CHKERRQ(ierr) + call DMcreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr) ! global solution vector (grid x 18, i.e. every def grad tensor) + call DMDASNESsetFunctionLocal(da,INSERT_VALUES,formResidual,PETSC_NULL_SNES,ierr) ! residual vector of same shape as solution vector + CHKERRQ(ierr) + call SNESsetConvergenceTest(snes,converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,ierr) ! specify custom convergence check function "converged" + CHKERRQ(ierr) + call SNESsetFromOptions(snes,ierr); CHKERRQ(ierr) ! pull it all together with additional CLI arguments !-------------------------------------------------------------------------------------------------- ! init fields - call DMDAVecGetArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) ! places pointer on PETSc data - F => FandF_tau( 0: 8,:,:,:) - F_tau => FandF_tau( 9:17,:,:,:) + call DMDAVecGetArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) ! places pointer on PETSc data + F => FandF_tau( 0: 8,:,:,:) + F_tau => FandF_tau( 9:17,:,:,:) + + restart: if (restartInc > 0) then + write(6,'(/,a,'//IO_intOut(restartInc)//',a)') ' reading values of increment ', restartInc, ' from file' + + write(rankStr,'(a1,i0)')'_',worldrank + fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5') + + call HDF5_read(fileHandle,F_aim, 'F_aim') + call HDF5_read(fileHandle,F_aim_lastInc,'F_aim_lastInc') + call HDF5_read(fileHandle,F_aimDot, 'F_aimDot') + call HDF5_read(fileHandle,F, 'F') + call HDF5_read(fileHandle,F_lastInc, 'F_lastInc') + call HDF5_read(fileHandle,F_tau, 'F_tau') + call HDF5_read(fileHandle,F_tau_lastInc,'F_tau_lastInc') + + elseif (restartInc == 0) then restart + F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity + F = reshape(F_lastInc,[9,grid(1),grid(2),grid3]) + F_tau = 2.0_pReal*F + F_tau_lastInc = 2.0_pReal*F_lastInc + endif restart + + materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent + call Utilities_updateIPcoords(reshape(F,shape(F_lastInc))) + call Utilities_constitutiveResponse(P,temp33_Real,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2 + reshape(F,shape(F_lastInc)), & ! target F + 0.0_pReal, & ! time increment + math_I3) ! no rotation of boundary condition + call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) ! deassociate pointer + + restartRead: if (restartInc > 0) then + write(6,'(/,a,'//IO_intOut(restartInc)//',a)') ' reading more values of increment ', restartInc, ' from file' + call HDF5_read(fileHandle,C_volAvg, 'C_volAvg') + call HDF5_read(fileHandle,C_volAvgLastInc,'C_volAvgLastInc') + call HDF5_closeFile(fileHandle) - restart: if (restartInc > 0) then - write(6,'(/,a,'//IO_intOut(restartInc)//',a)') ' reading values of increment ', restartInc, ' from file' - - fileUnit = IO_open_jobFile_binary('F_aim') - read(fileUnit) F_aim; close(fileUnit) - fileUnit = IO_open_jobFile_binary('F_aim_lastInc') - read(fileUnit) F_aim_lastInc; close(fileUnit) - fileUnit = IO_open_jobFile_binary('F_aimDot') - read(fileUnit) F_aimDot; close(fileUnit) - - write(rankStr,'(a1,i0)')'_',worldrank - - fileUnit = IO_open_jobFile_binary('F'//trim(rankStr)) - read(fileUnit) F; close (fileUnit) - fileUnit = IO_open_jobFile_binary('F_lastInc'//trim(rankStr)) - read(fileUnit) F_lastInc; close (fileUnit) - fileUnit = IO_open_jobFile_binary('F_tau'//trim(rankStr)) - read(fileUnit) F_tau; close (fileUnit) - fileUnit = IO_open_jobFile_binary('F_tau_lastInc'//trim(rankStr)) - read(fileUnit) F_tau_lastInc; close (fileUnit) - - elseif (restartInc == 0) then restart - F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity - F = reshape(F_lastInc,[9,grid(1),grid(2),grid3]) - F_tau = 2.0_pReal*F - F_tau_lastInc = 2.0_pReal*F_lastInc - endif restart - - materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent - call Utilities_updateIPcoords(reshape(F,shape(F_lastInc))) - call Utilities_constitutiveResponse(P,temp33_Real,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2 - reshape(F,shape(F_lastInc)), & ! target F - 0.0_pReal, & ! time increment - math_I3) ! no rotation of boundary condition - call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) ! deassociate pointer - - restartRead: if (restartInc > 0) then - write(6,'(/,a,'//IO_intOut(restartInc)//',a)') ' reading more values of increment ', restartInc, ' from file' - fileUnit = IO_open_jobFile_binary('C_volAvg') - read(fileUnit) C_volAvg; close(fileUnit) - fileUnit = IO_open_jobFile_binary('C_volAvgLastInv') - read(fileUnit) C_volAvgLastInc; close(fileUnit) - fileUnit = IO_open_jobFile_binary('C_ref') - read(fileUnit) C_minMaxAvg; close(fileUnit) - endif restartRead - - call Utilities_updateGamma(C_minMaxAvg,.true.) - C_scale = C_minMaxAvg - S_scale = math_invSym3333(C_minMaxAvg) + fileUnit = IO_open_jobFile_binary('C_ref') + read(fileUnit) C_minMaxAvg; close(fileUnit) + endif restartRead + + call utilities_updateGamma(C_minMaxAvg,.true.) + C_scale = C_minMaxAvg + S_scale = math_invSym3333(C_minMaxAvg) end subroutine grid_mech_spectral_polarisation_init @@ -237,66 +230,64 @@ end subroutine grid_mech_spectral_polarisation_init !> @brief solution for the Polarisation scheme with internal iterations !-------------------------------------------------------------------------------------------------- function grid_mech_spectral_polarisation_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC) result(solution) - use math, only: & - math_invSym3333 - use spectral_utilities, only: & - tBoundaryCondition, & - utilities_maskedCompliance, & - utilities_updateGamma - use FEsolving, only: & - restartWrite, & - terminallyIll - - implicit none + use math, only: & + math_invSym3333 + use spectral_utilities, only: & + tBoundaryCondition, & + utilities_maskedCompliance, & + utilities_updateGamma + use FEsolving, only: & + restartWrite, & + terminallyIll !-------------------------------------------------------------------------------------------------- ! input data for solution - character(len=*), intent(in) :: & - incInfoIn - real(pReal), intent(in) :: & - timeinc, & !< time increment of current solution - timeinc_old !< time increment of last successful increment - type(tBoundaryCondition), intent(in) :: & - stress_BC - real(pReal), dimension(3,3), intent(in) :: rotation_BC - type(tSolutionState) :: & - solution + character(len=*), intent(in) :: & + incInfoIn + real(pReal), intent(in) :: & + timeinc, & !< time increment of current solution + timeinc_old !< time increment of last successful increment + type(tBoundaryCondition), intent(in) :: & + stress_BC + real(pReal), dimension(3,3), intent(in) :: rotation_BC + type(tSolutionState) :: & + solution !-------------------------------------------------------------------------------------------------- ! PETSc Data - PetscErrorCode :: ierr - SNESConvergedReason :: reason - - incInfo = incInfoIn + PetscErrorCode :: ierr + SNESConvergedReason :: reason + + incInfo = incInfoIn !-------------------------------------------------------------------------------------------------- ! update stiffness (and gamma operator) - S = Utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg) - if (num%update_gamma) then - call utilities_updateGamma(C_minMaxAvg,restartWrite) - C_scale = C_minMaxAvg - S_scale = math_invSym3333(C_minMaxAvg) - endif + S = utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg) + if (num%update_gamma) then + call utilities_updateGamma(C_minMaxAvg,restartWrite) + C_scale = C_minMaxAvg + S_scale = math_invSym3333(C_minMaxAvg) + endif !-------------------------------------------------------------------------------------------------- ! set module wide available data - params%stress_mask = stress_BC%maskFloat - params%stress_BC = stress_BC%values - params%rotation_BC = rotation_BC - params%timeinc = timeinc - params%timeincOld = timeinc_old + params%stress_mask = stress_BC%maskFloat + params%stress_BC = stress_BC%values + params%rotation_BC = rotation_BC + params%timeinc = timeinc + params%timeincOld = timeinc_old !-------------------------------------------------------------------------------------------------- ! solve BVP - call SNESsolve(snes,PETSC_NULL_VEC,solution_vec,ierr); CHKERRQ(ierr) + call SNESsolve(snes,PETSC_NULL_VEC,solution_vec,ierr); CHKERRQ(ierr) !-------------------------------------------------------------------------------------------------- ! check convergence - call SNESGetConvergedReason(snes,reason,ierr); CHKERRQ(ierr) - - solution%converged = reason > 0 - solution%iterationsNeeded = totalIter - solution%termIll = terminallyIll - terminallyIll = .false. + call SNESGetConvergedReason(snes,reason,ierr); CHKERRQ(ierr) + + solution%converged = reason > 0 + solution%iterationsNeeded = totalIter + solution%termIll = terminallyIll + terminallyIll = .false. end function grid_mech_spectral_polarisation_solution @@ -330,13 +321,12 @@ subroutine grid_mech_spectral_polarisation_forward(guess,timeinc,timeinc_old,loa use FEsolving, only: & restartWrite - implicit none logical, intent(in) :: & guess real(pReal), intent(in) :: & timeinc_old, & timeinc, & - loadCaseTime !< remaining time of current load case + loadCaseTime !< remaining time of current load case type(tBoundaryCondition), intent(in) :: & stress_BC, & deformation_BC @@ -347,7 +337,7 @@ subroutine grid_mech_spectral_polarisation_forward(guess,timeinc,timeinc_old,loa integer :: i, j, k real(pReal), dimension(3,3) :: F_lambda33 - integer :: fileUnit + integer(HID_T) :: fileHandle character(len=32) :: rankStr call DMDAVecGetArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) @@ -355,37 +345,29 @@ subroutine grid_mech_spectral_polarisation_forward(guess,timeinc,timeinc_old,loa F_tau => FandF_tau( 9:17,:,:,:) if (cutBack) then - C_volAvg = C_volAvgLastInc ! QUESTION: where is this required? - C_minMaxAvg = C_minMaxAvgLastInc ! QUESTION: where is this required? + C_volAvg = C_volAvgLastInc ! QUESTION: where is this required? + C_minMaxAvg = C_minMaxAvgLastInc ! QUESTION: where is this required? else !-------------------------------------------------------------------------------------------------- ! restart information for spectral solver - if (restartWrite) then ! QUESTION: where is this logical properly set? - write(6,'(/,a)') ' writing converged results for restart' - flush(6) - - if (worldrank == 0) then - fileUnit = IO_open_jobFile_binary('C_volAvg','w') - write(fileUnit) C_volAvg; close(fileUnit) - fileUnit = IO_open_jobFile_binary('C_volAvgLastInv','w') - write(fileUnit) C_volAvgLastInc; close(fileUnit) - fileUnit = IO_open_jobFile_binary('F_aim','w') - write(fileUnit) F_aim; close(fileUnit) - fileUnit = IO_open_jobFile_binary('F_aim_lastInc','w') - write(fileUnit) F_aim_lastInc; close(fileUnit) - fileUnit = IO_open_jobFile_binary('F_aimDot','w') - write(fileUnit) F_aimDot; close(fileUnit) - endif + if (restartWrite) then + write(6,'(/,a)') ' writing converged results for restart';flush(6) write(rankStr,'(a1,i0)')'_',worldrank - fileUnit = IO_open_jobFile_binary('F'//trim(rankStr),'w') - write(fileUnit) F; close (fileUnit) - fileUnit = IO_open_jobFile_binary('F_lastInc'//trim(rankStr),'w') - write(fileUnit) F_lastInc; close (fileUnit) - fileUnit = IO_open_jobFile_binary('F_tau'//trim(rankStr),'w') - write(fileUnit) F_tau; close (fileUnit) - fileUnit = IO_open_jobFile_binary('F_tau_lastInc'//trim(rankStr),'w') - write(fileUnit) F_tau_lastInc; close (fileUnit) + fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5','w') + + call HDF5_write(fileHandle,F_aim, 'F_aim') + call HDF5_write(fileHandle,F_aim_lastInc, 'F_aim_lastInc') + call HDF5_write(fileHandle,F_aimDot, 'F_aimDot') + call HDF5_write(fileHandle,F, 'F') + call HDF5_write(fileHandle,F_lastInc, 'F_lastInc') + call HDF5_write(fileHandle,F_tau, 'F_tau') + call HDF5_write(fileHandle,F_tau_lastInc, 'F_tau_lastInc') + + call HDF5_write(fileHandle,C_volAvg, 'C_volAvg') + call HDF5_write(fileHandle,C_volAvgLastInc,'C_volAvgLastInc') + + call HDF5_closeFile(fileHandle) endif call CPFEM_age ! age state and kinematics @@ -399,13 +381,13 @@ subroutine grid_mech_spectral_polarisation_forward(guess,timeinc,timeinc_old,loa !-------------------------------------------------------------------------------------------------- ! calculate rate for aim - if (deformation_BC%myType=='l') then ! calculate F_aimDot from given L and current F + if (deformation_BC%myType=='l') then ! calculate F_aimDot from given L and current F F_aimDot = & F_aimDot + deformation_BC%maskFloat * matmul(deformation_BC%values, F_aim_lastInc) - elseif(deformation_BC%myType=='fdot') then ! F_aimDot is prescribed + elseif(deformation_BC%myType=='fdot') then ! F_aimDot is prescribed F_aimDot = & F_aimDot + deformation_BC%maskFloat * deformation_BC%values - elseif (deformation_BC%myType=='f') then ! aim at end of load case is prescribed + elseif (deformation_BC%myType=='f') then ! aim at end of load case is prescribed F_aimDot = & F_aimDot + deformation_BC%maskFloat * (deformation_BC%values - F_aim_lastInc)/loadCaseTime endif @@ -417,33 +399,33 @@ subroutine grid_mech_spectral_polarisation_forward(guess,timeinc,timeinc_old,loa F_tauDot = utilities_calculateRate(guess, & F_tau_lastInc,reshape(F_tau,[3,3,grid(1),grid(2),grid3]), timeinc_old, & math_rotate_backward33(F_aimDot,rotation_BC)) - F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid3]) ! winding F forward - F_tau_lastInc = reshape(F_tau, [3,3,grid(1),grid(2),grid3]) ! winding F_tau forward - materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent + F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid3]) ! winding F forward + F_tau_lastInc = reshape(F_tau, [3,3,grid(1),grid(2),grid3]) ! winding F_tau forward + materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent endif !-------------------------------------------------------------------------------------------------- ! update average and local deformation gradients F_aim = F_aim_lastInc + F_aimDot * timeinc - F = reshape(utilities_forwardField(timeinc,F_lastInc,Fdot, & ! estimate of F at end of time+timeinc that matches rotated F_aim on average + F = reshape(utilities_forwardField(timeinc,F_lastInc,Fdot, & ! estimate of F at end of time+timeinc that matches rotated F_aim on average math_rotate_backward33(F_aim,rotation_BC)),& [9,grid(1),grid(2),grid3]) - if (guess) then - F_tau = reshape(Utilities_forwardField(timeinc,F_tau_lastInc,F_taudot), & - [9,grid(1),grid(2),grid3]) ! does not have any average value as boundary condition - else - do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1) - F_lambda33 = reshape(F_tau(1:9,i,j,k)-F(1:9,i,j,k),[3,3]) - F_lambda33 = math_mul3333xx33(S_scale,matmul(F_lambda33, & - math_mul3333xx33(C_scale,& - matmul(transpose(F_lambda33),& - F_lambda33)-math_I3))*0.5_pReal)& - + math_I3 - F_tau(1:9,i,j,k) = reshape(F_lambda33,[9])+F(1:9,i,j,k) - enddo; enddo; enddo - endif - - call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) + if (guess) then + F_tau = reshape(Utilities_forwardField(timeinc,F_tau_lastInc,F_taudot), & + [9,grid(1),grid(2),grid3]) ! does not have any average value as boundary condition + else + do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1) + F_lambda33 = reshape(F_tau(1:9,i,j,k)-F(1:9,i,j,k),[3,3]) + F_lambda33 = math_mul3333xx33(S_scale,matmul(F_lambda33, & + math_mul3333xx33(C_scale,& + matmul(transpose(F_lambda33),& + F_lambda33)-math_I3))*0.5_pReal)& + + math_I3 + F_tau(1:9,i,j,k) = reshape(F_lambda33,[9])+F(1:9,i,j,k) + enddo; enddo; enddo + endif + + call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) end subroutine grid_mech_spectral_polarisation_forward @@ -451,208 +433,207 @@ end subroutine grid_mech_spectral_polarisation_forward !-------------------------------------------------------------------------------------------------- !> @brief convergence check !-------------------------------------------------------------------------------------------------- -subroutine converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr) - use numerics, only: & - itmax, & - itmin, & - err_div_tolRel, & - err_div_tolAbs, & - err_curl_tolRel, & - err_curl_tolAbs, & - err_stress_tolRel, & - err_stress_tolAbs - use FEsolving, only: & - terminallyIll - - implicit none - SNES :: snes_local - PetscInt :: PETScIter - PetscReal :: & - xnorm, & ! not used - snorm, & ! not used - fnorm ! not used - SNESConvergedReason :: reason - PetscObject :: dummy - PetscErrorCode :: ierr - real(pReal) :: & - curlTol, & - divTol, & - BCTol - - curlTol = max(maxval(abs(F_aim-math_I3))*err_curl_tolRel ,err_curl_tolAbs) - divTol = max(maxval(abs(P_av)) *err_div_tolRel ,err_div_tolAbs) - BCTol = max(maxval(abs(P_av)) *err_stress_tolRel,err_stress_tolAbs) - - if ((totalIter >= itmin .and. & - all([ err_div /divTol, & - err_curl/curlTol, & - err_BC /BCTol ] < 1.0_pReal)) & - .or. terminallyIll) then - reason = 1 - elseif (totalIter >= itmax) then - reason = -1 - else - reason = 0 - endif +subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dummy,ierr) + use numerics, only: & + itmax, & + itmin, & + err_div_tolRel, & + err_div_tolAbs, & + err_curl_tolRel, & + err_curl_tolAbs, & + err_stress_tolRel, & + err_stress_tolAbs + use FEsolving, only: & + terminallyIll + + SNES :: snes_local + PetscInt, intent(in) :: PETScIter + PetscReal, intent(in) :: & + devNull1, & + devNull2, & + devNull3 + SNESConvergedReason :: reason + PetscObject :: dummy + PetscErrorCode :: ierr + real(pReal) :: & + curlTol, & + divTol, & + BCTol + + curlTol = max(maxval(abs(F_aim-math_I3))*err_curl_tolRel ,err_curl_tolAbs) + divTol = max(maxval(abs(P_av)) *err_div_tolRel ,err_div_tolAbs) + BCTol = max(maxval(abs(P_av)) *err_stress_tolRel,err_stress_tolAbs) + + if ((totalIter >= itmin .and. & + all([ err_div /divTol, & + err_curl/curlTol, & + err_BC /BCTol ] < 1.0_pReal)) & + .or. terminallyIll) then + reason = 1 + elseif (totalIter >= itmax) then + reason = -1 + else + reason = 0 + endif !-------------------------------------------------------------------------------------------------- ! report - write(6,'(1/,a)') ' ... reporting .............................................................' - write(6,'(1/,a,f12.2,a,es8.2,a,es9.2,a)') ' error divergence = ', & - err_div/divTol, ' (',err_div, ' / m, tol = ',divTol,')' - write(6, '(a,f12.2,a,es8.2,a,es9.2,a)') ' error curl = ', & - err_curl/curlTol,' (',err_curl,' -, tol = ',curlTol,')' - write(6, '(a,f12.2,a,es8.2,a,es9.2,a)') ' error BC = ', & - err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')' - write(6,'(/,a)') ' ===========================================================================' - flush(6) + write(6,'(1/,a)') ' ... reporting .............................................................' + write(6,'(1/,a,f12.2,a,es8.2,a,es9.2,a)') ' error divergence = ', & + err_div/divTol, ' (',err_div, ' / m, tol = ',divTol,')' + write(6, '(a,f12.2,a,es8.2,a,es9.2,a)') ' error curl = ', & + err_curl/curlTol,' (',err_curl,' -, tol = ',curlTol,')' + write(6, '(a,f12.2,a,es8.2,a,es9.2,a)') ' error BC = ', & + err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')' + write(6,'(/,a)') ' ===========================================================================' + flush(6) end subroutine converged !-------------------------------------------------------------------------------------------------- -!> @brief forms the polarisation residual vector +!> @brief forms the residual vector !-------------------------------------------------------------------------------------------------- subroutine formResidual(in, FandF_tau, & residuum, dummy,ierr) - use numerics, only: & - itmax, & - itmin, & - polarAlpha, & - polarBeta - use mesh, only: & - grid, & - grid3 - use math, only: & - math_rotate_forward33, & - math_rotate_backward33, & - math_mul3333xx33, & - math_invSym3333 - use debug, only: & - debug_level, & - debug_spectral, & - debug_spectralRotation - use spectral_utilities, only: & - wgt, & - tensorField_real, & - utilities_FFTtensorForward, & - utilities_fourierGammaConvolution, & - utilities_FFTtensorBackward, & - utilities_constitutiveResponse, & - utilities_divergenceRMS, & - utilities_curlRMS - use IO, only: & - IO_intOut - use homogenization, only: & - materialpoint_dPdF - use FEsolving, only: & - terminallyIll + use numerics, only: & + itmax, & + itmin, & + polarAlpha, & + polarBeta + use mesh, only: & + grid, & + grid3 + use math, only: & + math_rotate_forward33, & + math_rotate_backward33, & + math_mul3333xx33, & + math_invSym3333 + use debug, only: & + debug_level, & + debug_spectral, & + debug_spectralRotation + use spectral_utilities, only: & + wgt, & + tensorField_real, & + utilities_FFTtensorForward, & + utilities_fourierGammaConvolution, & + utilities_FFTtensorBackward, & + utilities_constitutiveResponse, & + utilities_divergenceRMS, & + utilities_curlRMS + use IO, only: & + IO_intOut + use homogenization, only: & + materialpoint_dPdF + use FEsolving, only: & + terminallyIll - implicit none - DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: in !< DMDA info (needs to be named "in" for macros like XRANGE to work) - PetscScalar, dimension(3,3,2,XG_RANGE,YG_RANGE,ZG_RANGE), & - target, intent(in) :: FandF_tau - PetscScalar, dimension(3,3,2,X_RANGE,Y_RANGE,Z_RANGE),& - target, intent(out) :: residuum !< residuum field - PetscScalar, pointer, dimension(:,:,:,:,:) :: & - F, & - F_tau, & - residual_F, & - residual_F_tau - PetscInt :: & - PETScIter, & - nfuncs - PetscObject :: dummy - PetscErrorCode :: ierr - integer :: & - i, j, k, e + DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: in !< DMDA info (needs to be named "in" for macros like XRANGE to work) + PetscScalar, dimension(3,3,2,XG_RANGE,YG_RANGE,ZG_RANGE), & + target, intent(in) :: FandF_tau + PetscScalar, dimension(3,3,2,X_RANGE,Y_RANGE,Z_RANGE),& + target, intent(out) :: residuum !< residuum field + PetscScalar, pointer, dimension(:,:,:,:,:) :: & + F, & + F_tau, & + residual_F, & + residual_F_tau + PetscInt :: & + PETScIter, & + nfuncs + PetscObject :: dummy + PetscErrorCode :: ierr + integer :: & + i, j, k, e - F => FandF_tau(1:3,1:3,1,& - XG_RANGE,YG_RANGE,ZG_RANGE) - F_tau => FandF_tau(1:3,1:3,2,& - XG_RANGE,YG_RANGE,ZG_RANGE) - residual_F => residuum(1:3,1:3,1,& - X_RANGE, Y_RANGE, Z_RANGE) - residual_F_tau => residuum(1:3,1:3,2,& - X_RANGE, Y_RANGE, Z_RANGE) + F => FandF_tau(1:3,1:3,1,& + XG_RANGE,YG_RANGE,ZG_RANGE) + F_tau => FandF_tau(1:3,1:3,2,& + XG_RANGE,YG_RANGE,ZG_RANGE) + residual_F => residuum(1:3,1:3,1,& + X_RANGE, Y_RANGE, Z_RANGE) + residual_F_tau => residuum(1:3,1:3,2,& + X_RANGE, Y_RANGE, Z_RANGE) - F_av = sum(sum(sum(F,dim=5),dim=4),dim=3) * wgt - call MPI_Allreduce(MPI_IN_PLACE,F_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) - - call SNESGetNumberFunctionEvals(snes,nfuncs,ierr); CHKERRQ(ierr) - call SNESGetIterationNumber(snes,PETScIter,ierr); CHKERRQ(ierr) + F_av = sum(sum(sum(F,dim=5),dim=4),dim=3) * wgt + call MPI_Allreduce(MPI_IN_PLACE,F_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) + + call SNESGetNumberFunctionEvals(snes,nfuncs,ierr); CHKERRQ(ierr) + call SNESGetIterationNumber(snes,PETScIter,ierr); CHKERRQ(ierr) - if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1 ! new increment + if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1 ! new increment !-------------------------------------------------------------------------------------------------- ! begin of new iteration - newIteration: if (totalIter <= PETScIter) then - totalIter = totalIter + 1 - write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') & - trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax - if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) & - write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & - ' deformation gradient aim (lab) =', transpose(math_rotate_backward33(F_aim,params%rotation_BC)) - write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & - ' deformation gradient aim =', transpose(F_aim) - flush(6) - endif newIteration + newIteration: if (totalIter <= PETScIter) then + totalIter = totalIter + 1 + write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') & + trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax + if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) & + write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & + ' deformation gradient aim (lab) =', transpose(math_rotate_backward33(F_aim,params%rotation_BC)) + write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & + ' deformation gradient aim =', transpose(F_aim) + flush(6) + endif newIteration !-------------------------------------------------------------------------------------------------- ! - tensorField_real = 0.0_pReal - do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1) - tensorField_real(1:3,1:3,i,j,k) = & - polarBeta*math_mul3333xx33(C_scale,F(1:3,1:3,i,j,k) - math_I3) -& - polarAlpha*matmul(F(1:3,1:3,i,j,k), & - math_mul3333xx33(C_scale,F_tau(1:3,1:3,i,j,k) - F(1:3,1:3,i,j,k) - math_I3)) - enddo; enddo; enddo + tensorField_real = 0.0_pReal + do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1) + tensorField_real(1:3,1:3,i,j,k) = & + polarBeta*math_mul3333xx33(C_scale,F(1:3,1:3,i,j,k) - math_I3) -& + polarAlpha*matmul(F(1:3,1:3,i,j,k), & + math_mul3333xx33(C_scale,F_tau(1:3,1:3,i,j,k) - F(1:3,1:3,i,j,k) - math_I3)) + enddo; enddo; enddo !-------------------------------------------------------------------------------------------------- ! doing convolution in Fourier space - call utilities_FFTtensorForward - call utilities_fourierGammaConvolution(math_rotate_backward33(polarBeta*F_aim,params%rotation_BC)) - call utilities_FFTtensorBackward + call utilities_FFTtensorForward + call utilities_fourierGammaConvolution(math_rotate_backward33(polarBeta*F_aim,params%rotation_BC)) + call utilities_FFTtensorBackward !-------------------------------------------------------------------------------------------------- ! constructing residual - residual_F_tau = polarBeta*F - tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) + residual_F_tau = polarBeta*F - tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) !-------------------------------------------------------------------------------------------------- ! evaluate constitutive response - call utilities_constitutiveResponse(residual_F, & ! "residuum" gets field of first PK stress (to save memory) - P_av,C_volAvg,C_minMaxAvg, & - F - residual_F_tau/polarBeta,params%timeinc,params%rotation_BC) - call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr) + call utilities_constitutiveResponse(residual_F, & ! "residuum" gets field of first PK stress (to save memory) + P_av,C_volAvg,C_minMaxAvg, & + F - residual_F_tau/polarBeta,params%timeinc,params%rotation_BC) + call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr) !-------------------------------------------------------------------------------------------------- ! stress BC handling - F_aim = F_aim - math_mul3333xx33(S, ((P_av - params%stress_BC))) ! S = 0.0 for no bc - err_BC = maxval(abs((1.0_pReal-params%stress_mask) * math_mul3333xx33(C_scale,F_aim & - -math_rotate_forward33(F_av,params%rotation_BC)) + & - params%stress_mask * (P_av-params%stress_BC))) ! mask = 0.0 for no bc + F_aim = F_aim - math_mul3333xx33(S, ((P_av - params%stress_BC))) ! S = 0.0 for no bc + err_BC = maxval(abs((1.0_pReal-params%stress_mask) * math_mul3333xx33(C_scale,F_aim & + -math_rotate_forward33(F_av,params%rotation_BC)) + & + params%stress_mask * (P_av-params%stress_BC))) ! mask = 0.0 for no bc ! calculate divergence - tensorField_real = 0.0_pReal - tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = residual_F !< stress field in disguise - call utilities_FFTtensorForward - err_div = Utilities_divergenceRMS() !< root mean squared error in divergence of stress + tensorField_real = 0.0_pReal + tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = residual_F !< stress field in disguise + call utilities_FFTtensorForward + err_div = Utilities_divergenceRMS() !< root mean squared error in divergence of stress + !-------------------------------------------------------------------------------------------------- ! constructing residual - e = 0 - do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1) - e = e + 1 - residual_F(1:3,1:3,i,j,k) = & - math_mul3333xx33(math_invSym3333(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,e) + C_scale), & - residual_F(1:3,1:3,i,j,k) - matmul(F(1:3,1:3,i,j,k), & - math_mul3333xx33(C_scale,F_tau(1:3,1:3,i,j,k) - F(1:3,1:3,i,j,k) - math_I3))) & - + residual_F_tau(1:3,1:3,i,j,k) - enddo; enddo; enddo + e = 0 + do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1) + e = e + 1 + residual_F(1:3,1:3,i,j,k) = & + math_mul3333xx33(math_invSym3333(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,e) + C_scale), & + residual_F(1:3,1:3,i,j,k) - matmul(F(1:3,1:3,i,j,k), & + math_mul3333xx33(C_scale,F_tau(1:3,1:3,i,j,k) - F(1:3,1:3,i,j,k) - math_I3))) & + + residual_F_tau(1:3,1:3,i,j,k) + enddo; enddo; enddo !-------------------------------------------------------------------------------------------------- ! calculating curl - tensorField_real = 0.0_pReal - tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = F - call utilities_FFTtensorForward - err_curl = Utilities_curlRMS() + tensorField_real = 0.0_pReal + tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = F + call utilities_FFTtensorForward + err_curl = Utilities_curlRMS() end subroutine formResidual diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90 index 31740feb9..e899fd89a 100644 --- a/src/grid/grid_thermal_spectral.f90 +++ b/src/grid/grid_thermal_spectral.f90 @@ -69,7 +69,6 @@ subroutine grid_thermal_spectral_init worldsize, & petsc_options - implicit none PetscInt, dimension(worldsize) :: localK integer :: i, j, k, cell DM :: thermal_grid @@ -167,7 +166,6 @@ function grid_thermal_spectral_solution(timeinc,timeinc_old,loadCaseTime) result use thermal_conduction, only: & thermal_conduction_putTemperatureAndItsRate - implicit none real(pReal), intent(in) :: & timeinc, & !< increment in time for current solution timeinc_old, & !< increment in time of last increment @@ -242,7 +240,6 @@ subroutine grid_thermal_spectral_forward thermal_conduction_getMassDensity, & thermal_conduction_getSpecificHeat - implicit none integer :: i, j, k, cell DM :: dm_local PetscScalar, dimension(:,:,:), pointer :: x_scal @@ -311,7 +308,6 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr) thermal_conduction_getMassDensity, & thermal_conduction_getSpecificHeat - implicit none DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: & in PetscScalar, dimension( & diff --git a/src/material.f90 b/src/material.f90 index 383462ae1..64c3c2529 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -183,14 +183,11 @@ module material integer(pInt), private :: & microstructure_maxNconstituents, & !< max number of constituents in any phase - texture_maxNgauss, & !< max number of Gauss components in any texture - texture_maxNfiber !< max number of Fiber components in any texture + texture_maxNgauss !< max number of Gauss components in any texture integer(pInt), dimension(:), allocatable, private :: & microstructure_Nconstituents, & !< number of constituents in each microstructure - texture_symmetry, & !< number of symmetric orientations per texture - texture_Ngauss, & !< number of Gauss components per texture - texture_Nfiber !< number of Fiber components per texture + texture_Ngauss !< number of Gauss components per texture integer(pInt), dimension(:,:), allocatable, private :: & microstructure_phase, & !< phase IDs of each microstructure @@ -200,9 +197,7 @@ module material microstructure_fraction !< vol fraction of each constituent in microstructure real(pReal), dimension(:,:,:), allocatable, private :: & - material_volume, & !< volume of each grain,IP,element texture_Gauss, & !< data of each Gauss component - texture_Fiber, & !< data of each Fiber component texture_transformation !< transformation for each texture logical, dimension(:), allocatable, private :: & @@ -807,31 +802,27 @@ subroutine material_parseTexture math_det33 implicit none - integer(pInt) :: section, gauss, fiber, j, t, i + integer(pInt) :: section, gauss, j, t, i character(len=65536), dimension(:), allocatable :: strings ! Values for given key in material config integer(pInt), dimension(:), allocatable :: chunkPos - allocate(texture_symmetry(size(config_texture)), source=1_pInt) allocate(texture_Ngauss(size(config_texture)), source=0_pInt) - allocate(texture_Nfiber(size(config_texture)), source=0_pInt) do t=1_pInt, size(config_texture) - texture_Ngauss(t) = config_texture(t)%countKeys('(gauss)') & - + config_texture(t)%countKeys('(random)') - texture_Nfiber(t) = config_texture(t)%countKeys('(fiber)') + texture_Ngauss(t) = config_texture(t)%countKeys('(gauss)') + if (config_texture(t)%keyExists('symmetry')) call IO_error(147,ext_msg='symmetry') + if (config_texture(t)%keyExists('(random)')) call IO_error(147,ext_msg='(random)') + if (config_texture(t)%keyExists('(fiber)')) call IO_error(147,ext_msg='(fiber)') enddo texture_maxNgauss = maxval(texture_Ngauss) - texture_maxNfiber = maxval(texture_Nfiber) allocate(texture_Gauss (5,texture_maxNgauss,size(config_texture)), source=0.0_pReal) - allocate(texture_Fiber (6,texture_maxNfiber,size(config_texture)), source=0.0_pReal) allocate(texture_transformation(3,3,size(config_texture)), source=0.0_pReal) texture_transformation = spread(math_I3,3,size(config_texture)) do t=1_pInt, size(config_texture) section = t gauss = 0_pInt - fiber = 0_pInt if (config_texture(t)%keyExists('axes')) then strings = config_texture(t)%getStrings('axes') @@ -856,10 +847,6 @@ subroutine material_parseTexture if(dNeq(math_det33(texture_transformation(1:3,1:3,t)),1.0_pReal)) call IO_error(157_pInt,t) endif - if (config_texture(t)%keyExists('symmetry')) call IO_error(147,ext_msg='symmetry') - if (config_texture(t)%keyExists('(random)')) call IO_error(147,ext_msg='(random)') - if (config_texture(t)%keyExists('(fiber)')) call IO_error(147,ext_msg='(fiber)') - if (config_texture(t)%keyExists('(gauss)')) then gauss = gauss + 1_pInt strings = config_texture(t)%getStrings('(gauss)',raw= .true.) @@ -873,10 +860,6 @@ subroutine material_parseTexture texture_Gauss(2,gauss,t) = IO_floatValue(strings(i),chunkPos,j+1_pInt)*inRad case('phi2') texture_Gauss(3,gauss,t) = IO_floatValue(strings(i),chunkPos,j+1_pInt)*inRad - case('scatter') - texture_Gauss(4,gauss,t) = IO_floatValue(strings(i),chunkPos,j+1_pInt)*inRad - case('fraction') - texture_Gauss(5,gauss,t) = IO_floatValue(strings(i),chunkPos,j+1_pInt) end select enddo enddo @@ -984,28 +967,19 @@ end subroutine material_allocateSourceState !! calculates the volume of the grains and deals with texture components !-------------------------------------------------------------------------------------------------- subroutine material_populateGrains - use prec, only: & - dEq use math, only: & + math_EulertoR, & math_RtoEuler, & - math_EulerToR, & math_mul33x33, & math_range use mesh, only: & - theMesh, & - mesh_ipVolume + theMesh use config, only: & config_homogenization, & config_microstructure, & - config_deallocate, & - homogenization_name, & - microstructure_name + config_deallocate use IO, only: & IO_error - use debug, only: & - debug_level, & - debug_material, & - debug_levelBasic implicit none integer(pInt), dimension (:,:), allocatable :: Ngrains @@ -1015,21 +989,17 @@ subroutine material_populateGrains randomOrder real(pReal), dimension (microstructure_maxNconstituents) :: & rndArray - real(pReal), dimension (:), allocatable :: volumeOfGrain real(pReal), dimension (:,:), allocatable :: orientationOfGrain real(pReal), dimension (3) :: orientation - real(pReal), dimension (3,3) :: symOrientation integer(pInt), dimension (:), allocatable :: phaseOfGrain, textureOfGrain - integer(pInt) :: t,e,i,g,j,m,c,r,homog,micro,sgn,hme, myDebug, & + integer(pInt) :: t,e,i,g,j,m,c,r,homog,micro,sgn,hme, & phaseID,textureID,dGrains,myNgrains,myNorientations,myNconstituents, & - grain,constituentGrain,ipGrain,symExtension, ip + grain,constituentGrain,ipGrain,ip real(pReal) :: deviation,extreme,rnd integer(pInt), dimension (:,:), allocatable :: Nelems ! counts number of elements in homog, micro array type(group_int), dimension (:,:), allocatable :: elemsOfHomogMicro ! lists element number in homog, micro array - myDebug = debug_level(debug_material) - allocate(material_volume(homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%Nelems), source=0.0_pReal) allocate(material_phase(homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%Nelems), source=0_pInt) allocate(material_texture(homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%Nelems), source=0_pInt) allocate(material_EulerAngles(3,homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%Nelems),source=0.0_pReal) @@ -1038,6 +1008,25 @@ subroutine material_populateGrains allocate(Nelems (size(config_homogenization),size(config_microstructure)), source=0_pInt) + do e = 1, theMesh%Nelems + do i = 1, theMesh%elem%nIPs + homog = theMesh%homogenizationAt(e) + micro = theMesh%microstructureAt(e) + do c = 1, homogenization_Ngrains(homog) + material_phase(c,i,e) = microstructure_phase(c,micro) + material_texture(c,i,e) = microstructure_texture(c,micro) + material_EulerAngles(1:3,c,i,e) = texture_Gauss(1:3,1,material_texture(c,i,e)) + material_EulerAngles(1:3,c,i,e) = math_RtoEuler( & ! translate back to Euler angles + math_mul33x33( & ! pre-multiply + math_EulertoR(material_EulerAngles(1:3,c,i,e)), & ! face-value orientation + texture_transformation(1:3,1:3,material_texture(c,i,e)) & ! and transformation matrix + ) & + ) + enddo + enddo + enddo + + return !-------------------------------------------------------------------------------------------------- ! precounting of elements for each homog/micro pair do e = 1_pInt, theMesh%Nelems @@ -1075,47 +1064,17 @@ subroutine material_populateGrains elemsOfHomogMicro(homog,micro)%p(Nelems(homog,micro)) = e ! remember elements active in this homog/micro pair enddo elementLooping - allocate(volumeOfGrain(maxval(Ngrains)), source=0.0_pReal) ! reserve memory for maximum case allocate(phaseOfGrain(maxval(Ngrains)), source=0_pInt) ! reserve memory for maximum case allocate(textureOfGrain(maxval(Ngrains)), source=0_pInt) ! reserve memory for maximum case allocate(orientationOfGrain(3,maxval(Ngrains)),source=0.0_pReal) ! reserve memory for maximum case - if (iand(myDebug,debug_levelBasic) /= 0_pInt) then - write(6,'(/,a/)') ' MATERIAL grain population' - write(6,'(a32,1x,a32,1x,a6)') 'homogenization_name','microstructure_name','grain#' - endif + homogenizationLoop: do homog = 1_pInt,size(config_homogenization) dGrains = homogenization_Ngrains(homog) ! grain number per material point - microstructureLoop: do micro = 1_pInt,size(config_microstructure) ! all pairs of homog and micro + microstructureLoop: do micro = 1_pInt,size(config_microstructure) ! all pairs of homog and micro activePair: if (Ngrains(homog,micro) > 0_pInt) then myNgrains = Ngrains(homog,micro) ! assign short name for total number of grains to populate myNconstituents = microstructure_Nconstituents(micro) ! assign short name for number of constituents - if (iand(myDebug,debug_levelBasic) /= 0_pInt) & - write(6,'(/,a32,1x,a32,1x,i6)') homogenization_name(homog),microstructure_name(micro),myNgrains - - -!-------------------------------------------------------------------------------------------------- -! calculate volume of each grain - - volumeOfGrain = 0.0_pReal - grain = 0_pInt - - do hme = 1_pInt, Nelems(homog,micro) - e = elemsOfHomogMicro(homog,micro)%p(hme) ! my combination of homog and micro, only perform calculations for elements with homog, micro combinations which is indexed in cpElemsindex - if (microstructure_elemhomo(micro)) then ! homogeneous distribution of grains over each element's IPs - volumeOfGrain(grain+1_pInt:grain+dGrains) = sum(mesh_ipVolume(1:theMesh%elem%nIPs,e))/& - real(dGrains,pReal) ! each grain combines size of all IPs in that element - grain = grain + dGrains ! wind forward by Ngrains@IP - else - forall (i = 1_pInt:theMesh%elem%nIPs) & ! loop over IPs - volumeOfGrain(grain+(i-1)*dGrains+1_pInt:grain+i*dGrains) = & - mesh_ipVolume(i,e)/real(dGrains,pReal) ! assign IPvolume/Ngrains@IP to all grains of IP - grain = grain + theMesh%elem%nIPs * dGrains ! wind forward by Nips*Ngrains@IP - endif - enddo - - if (grain /= myNgrains) & - call IO_error(0,el = homog,ip = micro,ext_msg = 'inconsistent grain count after volume calc') !-------------------------------------------------------------------------------------------------- ! divide myNgrains as best over constituents @@ -1165,8 +1124,7 @@ subroutine material_populateGrains phaseOfGrain (grain+1_pInt:grain+NgrainsOfConstituent(i)) = phaseID ! assign resp. phase textureOfGrain(grain+1_pInt:grain+NgrainsOfConstituent(i)) = textureID ! assign resp. texture - myNorientations = ceiling(real(NgrainsOfConstituent(i),pReal)/& - real(texture_symmetry(textureID),pReal),pInt) ! max number of unique orientations (excl. symmetry) + myNorientations = ceiling(real(NgrainsOfConstituent(i),pReal)/1.0,pInt) ! max number of unique orientations (excl. symmetry) !-------------------------------------------------------------------------------------------------- ! has texture components @@ -1196,7 +1154,7 @@ subroutine material_populateGrains do j = 1_pInt,NgrainsOfConstituent(i)-1_pInt ! walk thru grains of current constituent call random_number(rnd) - t = nint(rnd*real(NgrainsOfConstituent(i)-j,pReal)+real(j,pReal)+0.5_pReal,pInt) ! select a grain in remaining list + t = nint(rnd*real(NgrainsOfConstituent(i)-j,pReal)+real(j,pReal)+0.5_pReal,pInt) ! select a grain in remaining list m = phaseOfGrain(grain+t) ! exchange current with random phaseOfGrain(grain+t) = phaseOfGrain(grain+j) phaseOfGrain(grain+j) = m @@ -1244,7 +1202,6 @@ subroutine material_populateGrains currentGrainOfConstituent(c))) ipGrain = ipGrain + 1_pInt ! advance IP grain counter currentGrainOfConstituent(c) = currentGrainOfConstituent(c) + 1_pInt ! advance index of grain population for constituent c - material_volume(ipGrain,i,e) = volumeOfGrain(grain+currentGrainOfConstituent(c)) ! assign properties material_phase(ipGrain,i,e) = phaseOfGrain(grain+currentGrainOfConstituent(c)) material_texture(ipGrain,i,e) = textureOfGrain(grain+currentGrainOfConstituent(c)) material_EulerAngles(1:3,ipGrain,i,e) = orientationOfGrain(1:3,grain+currentGrainOfConstituent(c)) @@ -1254,7 +1211,6 @@ subroutine material_populateGrains grain = sum(NgrainsOfConstituent(1:c-1_pInt)) ! figure out actual starting index in overall/consecutive grain population do ipGrain = ipGrain + 1_pInt, dGrains ! ensure last constituent fills up to dGrains currentGrainOfConstituent(c) = currentGrainOfConstituent(c) + 1_pInt - material_volume(ipGrain,i,e) = volumeOfGrain(grain+currentGrainOfConstituent(c)) material_phase(ipGrain,i,e) = phaseOfGrain(grain+currentGrainOfConstituent(c)) material_texture(ipGrain,i,e) = textureOfGrain(grain+currentGrainOfConstituent(c)) material_EulerAngles(1:3,ipGrain,i,e) = orientationOfGrain(1:3,grain+currentGrainOfConstituent(c)) @@ -1263,7 +1219,6 @@ subroutine material_populateGrains enddo do i = i, theMesh%elem%nIPs ! loop over IPs to (possibly) distribute copies from first IP - material_volume (1_pInt:dGrains,i,e) = material_volume (1_pInt:dGrains,1,e) material_phase (1_pInt:dGrains,i,e) = material_phase (1_pInt:dGrains,1,e) material_texture(1_pInt:dGrains,i,e) = material_texture(1_pInt:dGrains,1,e) material_EulerAngles(1:3,1_pInt:dGrains,i,e) = material_EulerAngles(1:3,1_pInt:dGrains,1,e) @@ -1275,7 +1230,7 @@ subroutine material_populateGrains enddo homogenizationLoop deallocate(texture_transformation) - deallocate(elemsOfHomogMicro) + call config_deallocate('material.config/microstructure') end subroutine material_populateGrains diff --git a/src/plastic_disloUCLA.f90 b/src/plastic_disloUCLA.f90 index efab7bc0b..56d910011 100644 --- a/src/plastic_disloUCLA.f90 +++ b/src/plastic_disloUCLA.f90 @@ -121,15 +121,7 @@ subroutine plastic_disloUCLA_init() math_expand use IO, only: & IO_error - use material, only: & - phase_plasticity, & - phase_plasticityInstance, & - phase_Noutput, & - material_allocatePlasticState, & - PLASTICITY_DISLOUCLA_label, & - PLASTICITY_DISLOUCLA_ID, & - material_phase, & - plasticState + use material use config, only: & config_phase use lattice diff --git a/src/plastic_dislotwin.f90 b/src/plastic_dislotwin.f90 index e8aca7f89..7db1f5f7f 100644 --- a/src/plastic_dislotwin.f90 +++ b/src/plastic_dislotwin.f90 @@ -185,15 +185,7 @@ subroutine plastic_dislotwin_init PI use IO, only: & IO_error - use material, only: & - phase_plasticity, & - phase_plasticityInstance, & - phase_Noutput, & - material_allocatePlasticState, & - PLASTICITY_DISLOTWIN_label, & - PLASTICITY_DISLOTWIN_ID, & - material_phase, & - plasticState + use material use config, only: & config_phase use lattice diff --git a/src/plastic_isotropic.f90 b/src/plastic_isotropic.f90 index 4f2892f57..c572f0ded 100644 --- a/src/plastic_isotropic.f90 +++ b/src/plastic_isotropic.f90 @@ -91,19 +91,7 @@ subroutine plastic_isotropic_init debug_levelBasic use IO, only: & IO_error - use material, only: & - phase_plasticity, & - phase_plasticityInstance, & - phase_Noutput, & - material_allocatePlasticState, & - PLASTICITY_ISOTROPIC_label, & - PLASTICITY_ISOTROPIC_ID, & - material_phase, & - plasticState -#ifdef DEBUG - use material, only: & - phasememberAt -#endif + use material use config, only: & config_phase use lattice @@ -464,7 +452,7 @@ function plastic_isotropic_postResults(Mp,instance,of) result(postResults) c = c + 1 case (dot_gamma_ID) postResults(c+1) = prm%dot_gamma_0 & - * (sqrt(1.5_pReal) * norm_Mp /(prm%M * stt%xi(of)))**prm%n + * (sqrt(1.5_pReal) * norm_Mp /(prm%M * stt%xi(of)))**prm%n c = c + 1 end select @@ -479,11 +467,10 @@ end function plastic_isotropic_postResults !> @brief writes results to HDF5 output file !-------------------------------------------------------------------------------------------------- subroutine plastic_isotropic_results(instance,group) -#if defined(PETSc) || defined(DAMASK_HDF5) - use results, only: & - results_writeDataset +#if defined(PETSc) || defined(DAMASKHDF5) + use results - integer, intent(in) :: instance + integer, intent(in) :: instance character(len=*), intent(in) :: group integer :: o @@ -496,10 +483,9 @@ subroutine plastic_isotropic_results(instance,group) end select enddo outputsLoop end associate - #else - integer, intent(in) :: instance - character(len=*), intent(in) :: group + integer, intent(in) :: instance + character(len=*) :: group #endif end subroutine plastic_isotropic_results diff --git a/src/plastic_kinematichardening.f90 b/src/plastic_kinematichardening.f90 index a255572e1..861b98da3 100644 --- a/src/plastic_kinematichardening.f90 +++ b/src/plastic_kinematichardening.f90 @@ -112,19 +112,7 @@ subroutine plastic_kinehardening_init math_expand use IO, only: & IO_error - use material, only: & - phase_plasticity, & - phase_plasticityInstance, & - phase_Noutput, & - material_allocatePlasticState, & - PLASTICITY_kinehardening_label, & - PLASTICITY_kinehardening_ID, & - material_phase, & - plasticState -#ifdef DEBUG - use material, only: & - phasememberAt -#endif + use material use config, only: & config_phase use lattice diff --git a/src/plastic_none.f90 b/src/plastic_none.f90 index 6f2ef2230..4b14266f1 100644 --- a/src/plastic_none.f90 +++ b/src/plastic_none.f90 @@ -23,13 +23,7 @@ subroutine plastic_none_init debug_level, & debug_constitutive, & debug_levelBasic - use material, only: & - phase_plasticity, & - material_allocatePlasticState, & - PLASTICITY_NONE_label, & - PLASTICITY_NONE_ID, & - material_phase, & - plasticState + use material integer :: & Ninstance, & diff --git a/src/plastic_nonlocal.f90 b/src/plastic_nonlocal.f90 index caaa0e4a2..6097bbbc8 100644 --- a/src/plastic_nonlocal.f90 +++ b/src/plastic_nonlocal.f90 @@ -248,15 +248,7 @@ subroutine plastic_nonlocal_init debug_levelBasic use mesh, only: & theMesh - use material, only: & - phase_plasticity, & - phase_plasticityInstance, & - phase_Noutput, & - PLASTICITY_NONLOCAL_label, & - PLASTICITY_NONLOCAL_ID, & - plasticState, & - material_phase, & - material_allocatePlasticState + use material use config use lattice diff --git a/src/plastic_phenopowerlaw.f90 b/src/plastic_phenopowerlaw.f90 index f8ebae68d..196129f64 100644 --- a/src/plastic_phenopowerlaw.f90 +++ b/src/plastic_phenopowerlaw.f90 @@ -116,15 +116,7 @@ subroutine plastic_phenopowerlaw_init math_expand use IO, only: & IO_error - use material, only: & - phase_plasticity, & - phase_plasticityInstance, & - phase_Noutput, & - material_allocatePlasticState, & - PLASTICITY_PHENOPOWERLAW_LABEL, & - PLASTICITY_PHENOPOWERLAW_ID, & - material_phase, & - plasticState + use material use config, only: & config_phase use lattice diff --git a/src/quit.f90 b/src/quit.f90 index 0215737e5..63184c113 100644 --- a/src/quit.f90 +++ b/src/quit.f90 @@ -34,13 +34,13 @@ subroutine quit(stop_id) #endif call date_and_time(values = dateAndTime) - write(6,'(/,a)') 'DAMASK terminated on:' - write(6,'(a,2(i2.2,a),i4.4)') 'Date: ',dateAndTime(3),'/',& - dateAndTime(2),'/',& - dateAndTime(1) - write(6,'(a,2(i2.2,a),i2.2)') 'Time: ',dateAndTime(5),':',& - dateAndTime(6),':',& - dateAndTime(7) + write(6,'(/,a)') ' DAMASK terminated on:' + write(6,'(a,2(i2.2,a),i4.4)') ' Date: ',dateAndTime(3),'/',& + dateAndTime(2),'/',& + dateAndTime(1) + write(6,'(a,2(i2.2,a),i2.2)') ' Time: ',dateAndTime(5),':',& + dateAndTime(6),':',& + dateAndTime(7) if (stop_id == 0 .and. ierr == 0 .and. error == 0) stop 0 ! normal termination if (stop_id == 2 .and. ierr == 0 .and. error == 0) stop 2 ! not all incs converged diff --git a/src/results.f90 b/src/results.f90 index 516c64552..c5582b927 100644 --- a/src/results.f90 +++ b/src/results.f90 @@ -66,9 +66,9 @@ subroutine results_init write(6,'(a)') ' https://doi.org/10.1007/s40192-018-0118-7' resultsFile = HDF5_openFile(trim(getSolverJobName())//'.hdf5','w',.true.) - call HDF5_addAttribute(resultsFile,'DADF5-version',0.1) + call HDF5_addAttribute(resultsFile,'DADF5-version',0.2) call HDF5_addAttribute(resultsFile,'DADF5-major',0) - call HDF5_addAttribute(resultsFile,'DADF5-minor',1) + call HDF5_addAttribute(resultsFile,'DADF5-minor',2) call HDF5_addAttribute(resultsFile,'DAMASK',DAMASKVERSION) call get_command(commandLine) call HDF5_addAttribute(resultsFile,'call',trim(commandLine)) @@ -103,7 +103,7 @@ end subroutine results_closeJobFile !-------------------------------------------------------------------------------------------------- -!> @brief closes the results file +!> @brief creates the group of increment and adds time as attribute to the file !-------------------------------------------------------------------------------------------------- subroutine results_addIncrement(inc,time) @@ -250,6 +250,8 @@ subroutine results_writeScalarDataset_real(group,dataset,label,description,SIuni #ifdef PETSc call HDF5_write(groupHandle,dataset,label,.true.) +#else + call HDF5_write(groupHandle,dataset,label,.false.) #endif if (HDF5_objectExists(groupHandle,label)) & @@ -277,6 +279,8 @@ subroutine results_writeVectorDataset_real(group,dataset,label,description,SIuni #ifdef PETSc call HDF5_write(groupHandle,dataset,label,.true.) +#else + call HDF5_write(groupHandle,dataset,label,.false.) #endif if (HDF5_objectExists(groupHandle,label)) & @@ -305,6 +309,8 @@ subroutine results_writeTensorDataset_real(group,dataset,label,description,SIuni #ifdef PETSc call HDF5_write(groupHandle,dataset,label,.true.) +#else + call HDF5_write(groupHandle,dataset,label,.false.) #endif if (HDF5_objectExists(groupHandle,label)) & @@ -333,6 +339,8 @@ subroutine results_writeVectorDataset_int(group,dataset,label,description,SIunit #ifdef PETSc call HDF5_write(groupHandle,dataset,label,.true.) +#else + call HDF5_write(groupHandle,dataset,label,.false.) #endif if (HDF5_objectExists(groupHandle,label)) & @@ -347,7 +355,7 @@ end subroutine results_writeVectorDataset_int !-------------------------------------------------------------------------------------------------- -!> @brief stores a vector dataset in a group +!> @brief stores a tensor dataset in a group !-------------------------------------------------------------------------------------------------- subroutine results_writeTensorDataset_int(group,dataset,label,description,SIunit) @@ -361,6 +369,8 @@ subroutine results_writeTensorDataset_int(group,dataset,label,description,SIunit #ifdef PETSc call HDF5_write(groupHandle,dataset,label,.true.) +#else + call HDF5_write(groupHandle,dataset,label,.false.) #endif if (HDF5_objectExists(groupHandle,label)) & @@ -375,7 +385,7 @@ end subroutine results_writeTensorDataset_int !-------------------------------------------------------------------------------------------------- -!> @brief stores a vector dataset in a group +!> @brief stores a scalar dataset in a group !-------------------------------------------------------------------------------------------------- subroutine results_writeScalarDataset_rotation(group,dataset,label,description,lattice_structure) use rotations, only: & @@ -391,6 +401,8 @@ subroutine results_writeScalarDataset_rotation(group,dataset,label,description,l #ifdef PETSc call HDF5_write(groupHandle,dataset,label,.true.) +#else + call HDF5_write(groupHandle,dataset,label,.false.) #endif if (HDF5_objectExists(groupHandle,label)) & diff --git a/src/rotations.f90 b/src/rotations.f90 index 7c29b03bc..69529ed24 100644 --- a/src/rotations.f90 +++ b/src/rotations.f90 @@ -77,7 +77,6 @@ contains !--------------------------------------------------------------------------------------------------- function asQuaternion(self) - implicit none class(rotation), intent(in) :: self real(pReal), dimension(4) :: asQuaternion @@ -87,7 +86,6 @@ end function asQuaternion !--------------------------------------------------------------------------------------------------- function asEulerAngles(self) - implicit none class(rotation), intent(in) :: self real(pReal), dimension(3) :: asEulerAngles @@ -97,7 +95,6 @@ end function asEulerAngles !--------------------------------------------------------------------------------------------------- function asAxisAnglePair(self) - implicit none class(rotation), intent(in) :: self real(pReal), dimension(4) :: asAxisAnglePair @@ -107,7 +104,6 @@ end function asAxisAnglePair !--------------------------------------------------------------------------------------------------- function asRotationMatrix(self) - implicit none class(rotation), intent(in) :: self real(pReal), dimension(3,3) :: asRotationMatrix @@ -117,7 +113,6 @@ end function asRotationMatrix !--------------------------------------------------------------------------------------------------- function asRodriguesFrankVector(self) - implicit none class(rotation), intent(in) :: self real(pReal), dimension(4) :: asRodriguesFrankVector @@ -127,7 +122,6 @@ end function asRodriguesFrankVector !--------------------------------------------------------------------------------------------------- function asHomochoric(self) - implicit none class(rotation), intent(in) :: self real(pReal), dimension(3) :: asHomochoric @@ -140,7 +134,6 @@ end function asHomochoric !--------------------------------------------------------------------------------------------------- subroutine fromRotationMatrix(self,om) - implicit none class(rotation), intent(out) :: self real(pReal), dimension(3,3), intent(in) :: om @@ -158,7 +151,6 @@ function rotVector(self,v,active) use prec, only: & dEq - implicit none real(pReal), dimension(3) :: rotVector class(rotation), intent(in) :: self real(pReal), intent(in), dimension(3) :: v @@ -198,7 +190,6 @@ end function rotVector !--------------------------------------------------------------------------------------------------- function rotTensor(self,m,active) - implicit none real(pReal), dimension(3,3) :: rotTensor class(rotation), intent(in) :: self real(pReal), intent(in), dimension(3,3) :: m @@ -226,7 +217,6 @@ end function rotTensor !--------------------------------------------------------------------------------------------------- function misorientation(self,other) - implicit none type(rotation) :: misorientation class(rotation), intent(in) :: self, other @@ -235,469 +225,12 @@ function misorientation(self,other) end function misorientation -!--------------------------------------------------------------------------------------------------- -! The following routines convert between different representations -!--------------------------------------------------------------------------------------------------- - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief Euler angles to orientation matrix -!--------------------------------------------------------------------------------------------------- -pure function eu2om(eu) result(om) - use prec, only: & - dEq0 - - implicit none - real(pReal), intent(in), dimension(3) :: eu - real(pReal), dimension(3,3) :: om - - 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) - om(2,1) = -c(1)*s(3)-s(1)*c(3)*c(2) - om(2,2) = -s(1)*s(3)+c(1)*c(3)*c(2) - om(2,3) = c(3)*s(2) - 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 - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief convert euler to axis angle -!--------------------------------------------------------------------------------------------------- -pure function eu2ax(eu) result(ax) - use prec, only: & - dEq0, & - dEq - use math, only: & - PI - - implicit none - 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) - sigma = 0.5*(eu(1)+eu(3)) - delta = 0.5*(eu(1)-eu(3)) - tau = sqrt(t**2+sin(sigma)**2) - - alpha = merge(PI, 2.0*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 - ax(1:3) = -P/tau * [ t*cos(delta), t*sin(delta), sin(sigma) ] ! passive axis-angle pair so a minus sign in front - ax(4) = alpha - if (alpha < 0.0) ax = -ax ! ensure alpha is positive - end if - -end function eu2ax - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief Euler angles to Rodrigues vector -!--------------------------------------------------------------------------------------------------- -pure function eu2ro(eu) result(ro) - use prec, only: & - dEq0 - use, intrinsic :: IEEE_ARITHMETIC, only: & - IEEE_value, & - IEEE_positive_inf - use math, only: & - PI - - implicit none - 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) - elseif(dEq0(ro(4))) then - ro = [ 0.0_pReal, 0.0_pReal, P, 0.0_pReal ] - else - ro(4) = tan(ro(4)*0.5) - end if - -end function eu2ro - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief Euler angles to unit quaternion -!--------------------------------------------------------------------------------------------------- -pure function eu2qu(eu) result(qu) - - implicit none - real(pReal), intent(in), dimension(3) :: eu - type(quaternion) :: qu - real(pReal), dimension(3) :: ee - real(pReal) :: cPhi, sPhi - - ee = 0.5_pReal*eu - - cPhi = cos(ee(2)) - sPhi = sin(ee(2)) - - qu = quaternion([ cPhi*cos(ee(1)+ee(3)), & - -P*sPhi*cos(ee(1)-ee(3)), & - -P*sPhi*sin(ee(1)-ee(3)), & - -P*cPhi*sin(ee(1)+ee(3))]) - if(qu%w < 0.0_pReal) qu = qu%homomorphed() - -end function eu2qu - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief orientation matrix to Euler angles -!--------------------------------------------------------------------------------------------------- -pure function om2eu(om) result(eu) - use math, only: & - PI - - implicit none - 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 - eu = [ atan2( om(1,2),om(1,1)), 0.5*PI*(1-om(3,3)),0.0_pReal ] - else - 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)] - 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 - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief convert axis angle pair to orientation matrix -!--------------------------------------------------------------------------------------------------- -pure function ax2om(ax) result(om) - - implicit none - real(pReal), intent(in), dimension(4) :: ax - real(pReal), dimension(3,3) :: om - - real(pReal) :: q, c, s, omc - integer :: i - - c = cos(ax(4)) - s = sin(ax(4)) - omc = 1.0-c - - forall(i=1:3) om(i,i) = ax(i)**2*omc + c - - 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) - - if (P > 0.0) om = transpose(om) - -end function ax2om - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief convert unit quaternion to Euler angles -!--------------------------------------------------------------------------------------------------- -pure function qu2eu(qu) result(eu) - use prec, only: & - dEq0 - use math, only: & - PI - - implicit none - type(quaternion), intent(in) :: qu - real(pReal), dimension(3) :: eu - - real(pReal) :: q12, q03, chi, chiInv - - q03 = qu%w**2+qu%z**2 - q12 = qu%x**2+qu%y**2 - chi = sqrt(q03*q12) - - degenerated: if (dEq0(chi)) then - eu = merge([atan2(-P*2.0*qu%w*qu%z,qu%w**2-qu%z**2), 0.0_pReal, 0.0_pReal], & - [atan2(2.0*qu%x*qu%y,qu%x**2-qu%y**2), PI, 0.0_pReal], & - dEq0(q12)) - else degenerated - chiInv = 1.0/chi - eu = [atan2((-P*qu%w*qu%y+qu%x*qu%z)*chi, (-P*qu%w*qu%x-qu%y*qu%z)*chi ), & - atan2( 2.0*chi, q03-q12 ), & - atan2(( P*qu%w*qu%y+qu%x*qu%z)*chi, (-P*qu%w*qu%x+qu%y*qu%z)*chi )] - endif degenerated - where(eu<0.0_pReal) eu = mod(eu+2.0_pReal*PI,[2.0_pReal*PI,PI,2.0_pReal*PI]) - -end function qu2eu - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief convert axis angle pair to homochoric -!--------------------------------------------------------------------------------------------------- -pure function ax2ho(ax) result(ho) - - implicit none - real(pReal), intent(in), dimension(4) :: ax - real(pReal), dimension(3) :: ho - - real(pReal) :: f - - f = 0.75 * ( ax(4) - sin(ax(4)) ) - f = f**(1.0/3.0) - ho = ax(1:3) * f - -end function ax2ho - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief convert homochoric to axis angle pair -!--------------------------------------------------------------------------------------------------- -pure function ho2ax(ho) result(ax) - use prec, only: & - dEq0 - - implicit none - 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, & - +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 - hm = hm*hmag_squared - s = s + tfit(i) * hm - end do - ax = [ho/sqrt(hmag_squared), 2.0_pReal*acos(s)] - end if - -end function ho2ax - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief convert orientation matrix to axis angle pair -!--------------------------------------------------------------------------------------------------- -function om2ax(om) result(ax) - use prec, only: & - dEq0, & - cEq, & - dNeq0 - use IO, only: & - IO_error - use math, only: & - math_clip, & - math_trace33 - - implicit none - real(pReal), intent(in) :: om(3,3) - real(pReal) :: ax(4) - - real(pReal) :: t - real(pReal), dimension(3) :: Wr, Wi - real(pReal), dimension(10) :: WORK - real(pReal), dimension(3,3) :: VR, devNull, o - integer :: INFO, LWORK, i - - external :: dgeev,sgeev - - o = om - - ! first get the rotation angle - t = 0.5_pReal * (math_trace33(om) - 1.0) - ax(4) = acos(math_clip(t,-1.0_pReal,1.0_pReal)) - - if (dEq0(ax(4))) then - ax(1:3) = [ 0.0, 0.0, 1.0 ] - else - ! set some initial LAPACK variables - INFO = 0 - ! first initialize the parameters for the LAPACK DGEEV routines - LWORK = 20 - - ! call the eigenvalue solver - call dgeev('N','V',3,o,3,Wr,Wi,devNull,3,VR,3,WORK,LWORK,INFO) - if (INFO /= 0) call IO_error(0,ext_msg='Error in om2ax DGEEV return not zero') - i = maxloc(merge(1.0_pReal,0.0_pReal,cEq(cmplx(Wr,Wi,pReal),cmplx(1.0_pReal,0.0_pReal,pReal),tol=1.0e-14_pReal)),dim=1) ! poor substitute for findloc - ax(1:3) = VR(1:3,i) - where ( dNeq0([om(2,3)-om(3,2), om(3,1)-om(1,3), om(1,2)-om(2,1)])) & - ax(1:3) = sign(ax(1:3),-P *[om(2,3)-om(3,2), om(3,1)-om(1,3), om(1,2)-om(2,1)]) - endif - -end function om2ax - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief convert Rodrigues vector to axis angle pair -!--------------------------------------------------------------------------------------------------- -pure function ro2ax(ro) result(ax) - use, intrinsic :: IEEE_ARITHMETIC, only: & - IEEE_is_finite - use prec, only: & - dEq0 - use math, only: & - PI - - implicit none - real(pReal), intent(in), dimension(4) :: ro - real(pReal), dimension(4) :: ax - - real(pReal) :: ta, angle - - ta = ro(4) - - if (dEq0(ta)) then - ax = [ 0.0, 0.0, 1.0, 0.0 ] - elseif (.not. IEEE_is_finite(ta)) then - ax = [ ro(1), ro(2), ro(3), PI ] - else - angle = 2.0*atan(ta) - ta = 1.0/norm2(ro(1:3)) - ax = [ ro(1)/ta, ro(2)/ta, ro(3)/ta, angle ] - end if - -end function ro2ax - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief convert axis angle pair to Rodrigues vector -!--------------------------------------------------------------------------------------------------- -pure function ax2ro(ax) result(ro) - use, intrinsic :: IEEE_ARITHMETIC, only: & - IEEE_value, & - IEEE_positive_inf - use prec, only: & - dEq0 - use math, only: & - PI - - implicit none - real(pReal), intent(in), dimension(4) :: ax - real(pReal), dimension(4) :: ro - - real(pReal), parameter :: thr = 1.0E-7 - - if (dEq0(ax(4))) then - ro = [ 0.0_pReal, 0.0_pReal, P, 0.0_pReal ] - 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 ),abs(ax(4)-PI) < thr) - end if - -end function ax2ro - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief convert axis angle pair to quaternion -!--------------------------------------------------------------------------------------------------- -pure function ax2qu(ax) result(qu) - use prec, only: & - dEq0 - - implicit none - real(pReal), intent(in), dimension(4) :: ax - type(quaternion) :: qu - - real(pReal) :: c, s - - - if (dEq0(ax(4))) then - qu = quaternion([ 1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal ]) - else - c = cos(ax(4)*0.5) - s = sin(ax(4)*0.5) - qu = quaternion([ c, ax(1)*s, ax(2)*s, ax(3)*s ]) - end if - -end function ax2qu - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief convert Rodrigues vector to homochoric -!--------------------------------------------------------------------------------------------------- -pure function ro2ho(ro) result(ho) - use, intrinsic :: IEEE_ARITHMETIC, only: & - IEEE_is_finite - use prec, only: & - dEq0 - use math, only: & - PI - - implicit none - 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, 0.0, 0.0 ] - else - f = merge(2.0*atan(ro(4)) - sin(2.0*atan(ro(4))),PI, IEEE_is_finite(ro(4))) - ho = ro(1:3) * (0.75_pReal*f)**(1.0/3.0) - end if - -end function ro2ho - - !--------------------------------------------------------------------------------------------------- !> @author Marc De Graef, Carnegie Mellon University !> @brief convert unit quaternion to rotation matrix !--------------------------------------------------------------------------------------------------- pure function qu2om(qu) result(om) - implicit none type(quaternion), intent(in) :: qu real(pReal), dimension(3,3) :: om @@ -724,39 +257,36 @@ end function qu2om !--------------------------------------------------------------------------------------------------- !> @author Marc De Graef, Carnegie Mellon University -!> @brief convert rotation matrix to a unit quaternion +!> @brief convert unit quaternion to Euler angles !--------------------------------------------------------------------------------------------------- -function om2qu(om) result(qu) +pure function qu2eu(qu) result(eu) use prec, only: & - dEq + dEq0 + use math, only: & + PI - implicit none - real(pReal), intent(in), dimension(3,3) :: om - type(quaternion) :: qu + type(quaternion), intent(in) :: qu + real(pReal), dimension(3) :: eu - real(pReal), dimension(4) :: qu_A - real(pReal), dimension(4) :: s - - s = [+om(1,1) +om(2,2) +om(3,3) +1.0_pReal, & - +om(1,1) -om(2,2) -om(3,3) +1.0_pReal, & - -om(1,1) +om(2,2) -om(3,3) +1.0_pReal, & - -om(1,1) -om(2,2) +om(3,3) +1.0_pReal] - - qu_A = sqrt(max(s,0.0_pReal))*0.5_pReal*[1.0_pReal,P,P,P] - qu_A = qu_A/norm2(qu_A) - - if(any(dEq(abs(qu_A),1.0_pReal,1.0e-15_pReal))) & - where (.not.(dEq(abs(qu_A),1.0_pReal,1.0e-15_pReal))) qu_A = 0.0_pReal - - if (om(3,2) < om(2,3)) qu_A(2) = -qu_A(2) - if (om(1,3) < om(3,1)) qu_A(3) = -qu_A(3) - if (om(2,1) < om(1,2)) qu_A(4) = -qu_A(4) + real(pReal) :: q12, q03, chi, chiInv - qu = quaternion(qu_A) - !qu_A = om2ax(om) - !if(any(qu_A(1:3) * [qu%x,qu%y,qu%z] < 0.0)) print*, 'sign error' + q03 = qu%w**2+qu%z**2 + q12 = qu%x**2+qu%y**2 + chi = sqrt(q03*q12) + + degenerated: if (dEq0(chi)) then + eu = merge([atan2(-P*2.0*qu%w*qu%z,qu%w**2-qu%z**2), 0.0_pReal, 0.0_pReal], & + [atan2(2.0*qu%x*qu%y,qu%x**2-qu%y**2), PI, 0.0_pReal], & + dEq0(q12)) + else degenerated + chiInv = 1.0/chi + eu = [atan2((-P*qu%w*qu%y+qu%x*qu%z)*chi, (-P*qu%w*qu%x-qu%y*qu%z)*chi ), & + atan2( 2.0*chi, q03-q12 ), & + atan2(( P*qu%w*qu%y+qu%x*qu%z)*chi, (-P*qu%w*qu%x+qu%y*qu%z)*chi )] + endif degenerated + where(eu<0.0_pReal) eu = mod(eu+2.0_pReal*PI,[2.0_pReal*PI,PI,2.0_pReal*PI]) -end function om2qu +end function qu2eu !--------------------------------------------------------------------------------------------------- @@ -771,15 +301,13 @@ pure function qu2ax(qu) result(ax) PI, & math_clip - implicit none type(quaternion), intent(in) :: qu real(pReal), dimension(4) :: ax real(pReal) :: omega, s - ! if the angle equals zero, then we return the rotation axis as [001] - if (dEq0(sqrt(qu%x**2+qu%y**2+qu%z**2))) then - ax = [ 0.0_pReal, 0.0_pReal, 1.0_pReal, 0.0_pReal ] + if (dEq0(qu%x**2+qu%y**2+qu%z**2)) then + ax = [ 0.0_pReal, 0.0_pReal, 1.0_pReal, 0.0_pReal ] ! axis = [001] elseif (dNeq0(qu%w)) then s = sign(1.0_pReal,qu%w)/sqrt(qu%x**2+qu%y**2+qu%z**2) omega = 2.0_pReal * acos(math_clip(qu%w,-1.0_pReal,1.0_pReal)) @@ -835,7 +363,6 @@ pure function qu2ho(qu) result(ho) use math, only: & math_clip - implicit none type(quaternion), intent(in) :: qu real(pReal), dimension(3) :: ho @@ -856,66 +383,109 @@ end function qu2ho !--------------------------------------------------------------------------------------------------- !> @author Marc De Graef, Carnegie Mellon University -!> @brief convert homochoric to cubochoric +!> @brief convert unit quaternion to cubochoric !--------------------------------------------------------------------------------------------------- -function ho2cu(ho) result(cu) - use Lambert, only: & - LambertBallToCube +function qu2cu(qu) result(cu) + + type(quaternion), intent(in) :: qu + real(pReal), dimension(3) :: cu - implicit none - real(pReal), intent(in), dimension(3) :: ho - real(pReal), dimension(3) :: cu + cu = ho2cu(qu2ho(qu)) - cu = LambertBallToCube(ho) +end function qu2cu -end function ho2cu + +!--------------------------------------------------------------------------------------------------- +!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH +!> @brief convert rotation matrix to cubochoric +!> @details the original formulation (direct conversion) had (numerical?) issues +!--------------------------------------------------------------------------------------------------- +function om2qu(om) result(qu) + + real(pReal), intent(in), dimension(3,3) :: om + type(quaternion) :: qu + + qu = eu2qu(om2eu(om)) + +end function om2qu !--------------------------------------------------------------------------------------------------- !> @author Marc De Graef, Carnegie Mellon University -!> @brief convert cubochoric to homochoric +!> @brief orientation matrix to Euler angles !--------------------------------------------------------------------------------------------------- -function cu2ho(cu) result(ho) - use Lambert, only: & - LambertCubeToBall +pure function om2eu(om) result(eu) + use math, only: & + PI - implicit none - real(pReal), intent(in), dimension(3) :: cu - real(pReal), dimension(3) :: ho + 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 + eu = [ atan2( om(1,2),om(1,1)), 0.5*PI*(1-om(3,3)),0.0_pReal ] + end if - ho = LambertCubeToBall(cu) - -end function cu2ho + where(eu<0.0_pReal) eu = mod(eu+2.0_pReal*PI,[2.0_pReal*PI,PI,2.0_pReal*PI]) + +end function om2eu !--------------------------------------------------------------------------------------------------- !> @author Marc De Graef, Carnegie Mellon University -!> @brief convert Rodrigues vector to Euler angles +!> @brief convert orientation matrix to axis angle pair !--------------------------------------------------------------------------------------------------- -pure function ro2eu(ro) result(eu) +function om2ax(om) result(ax) + use prec, only: & + dEq0, & + cEq, & + dNeq0 + use IO, only: & + IO_error + use math, only: & + math_clip, & + math_trace33 - implicit none - real(pReal), intent(in), dimension(4) :: ro - real(pReal), dimension(3) :: eu - - eu = om2eu(ro2om(ro)) + real(pReal), intent(in) :: om(3,3) + real(pReal) :: ax(4) + + real(pReal) :: t + real(pReal), dimension(3) :: Wr, Wi + real(pReal), dimension(10) :: WORK + real(pReal), dimension(3,3) :: VR, devNull, o + integer :: INFO, LWORK, i + + external :: dgeev,sgeev + + o = om + + ! first get the rotation angle + t = 0.5_pReal * (math_trace33(om) - 1.0) + ax(4) = acos(math_clip(t,-1.0_pReal,1.0_pReal)) + + if (dEq0(ax(4))) then + ax(1:3) = [ 0.0, 0.0, 1.0 ] + else + ! set some initial LAPACK variables + INFO = 0 + ! first initialize the parameters for the LAPACK DGEEV routines + LWORK = 20 + + ! call the eigenvalue solver + call dgeev('N','V',3,o,3,Wr,Wi,devNull,3,VR,3,WORK,LWORK,INFO) + if (INFO /= 0) call IO_error(0,ext_msg='Error in om2ax DGEEV return not zero') + i = maxloc(merge(1.0_pReal,0.0_pReal,cEq(cmplx(Wr,Wi,pReal),cmplx(1.0_pReal,0.0_pReal,pReal),tol=1.0e-14_pReal)),dim=1) ! poor substitute for findloc + ax(1:3) = VR(1:3,i) + where ( dNeq0([om(2,3)-om(3,2), om(3,1)-om(1,3), om(1,2)-om(2,1)])) & + ax(1:3) = sign(ax(1:3),-P *[om(2,3)-om(3,2), om(3,1)-om(1,3), om(1,2)-om(2,1)]) + endif -end function ro2eu - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief convert Euler angles to homochoric -!--------------------------------------------------------------------------------------------------- -pure function eu2ho(eu) result(ho) - - implicit none - real(pReal), intent(in), dimension(3) :: eu - real(pReal), dimension(3) :: ho - - ho = ax2ho(eu2ax(eu)) - -end function eu2ho +end function om2ax !--------------------------------------------------------------------------------------------------- @@ -924,7 +494,6 @@ end function eu2ho !--------------------------------------------------------------------------------------------------- pure function om2ro(om) result(ro) - implicit none real(pReal), intent(in), dimension(3,3) :: om real(pReal), dimension(4) :: ro @@ -939,133 +508,12 @@ end function om2ro !--------------------------------------------------------------------------------------------------- function om2ho(om) result(ho) - implicit none real(pReal), intent(in), dimension(3,3) :: om real(pReal), dimension(3) :: ho ho = ax2ho(om2ax(om)) -end function om2ho - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief convert axis angle pair to Euler angles -!--------------------------------------------------------------------------------------------------- -pure function ax2eu(ax) result(eu) - - implicit none - real(pReal), intent(in), dimension(4) :: ax - real(pReal), dimension(3) :: eu - - eu = om2eu(ax2om(ax)) - -end function ax2eu - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief convert Rodrigues vector to rotation matrix -!--------------------------------------------------------------------------------------------------- -pure function ro2om(ro) result(om) - - implicit none - real(pReal), intent(in), dimension(4) :: ro - real(pReal), dimension(3,3) :: om - - om = ax2om(ro2ax(ro)) - -end function ro2om - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief convert Rodrigues vector to unit quaternion -!--------------------------------------------------------------------------------------------------- -pure function ro2qu(ro) result(qu) - - implicit none - real(pReal), intent(in), dimension(4) :: ro - type(quaternion) :: qu - - qu = ax2qu(ro2ax(ro)) - -end function ro2qu - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief convert homochoric to Euler angles -!--------------------------------------------------------------------------------------------------- -pure function ho2eu(ho) result(eu) - - implicit none - real(pReal), intent(in), dimension(3) :: ho - real(pReal), dimension(3) :: eu - - eu = ax2eu(ho2ax(ho)) - -end function ho2eu - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief convert homochoric to rotation matrix -!--------------------------------------------------------------------------------------------------- -pure function ho2om(ho) result(om) - - implicit none - real(pReal), intent(in), dimension(3) :: ho - real(pReal), dimension(3,3) :: om - - om = ax2om(ho2ax(ho)) - -end function ho2om - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief convert homochoric to Rodrigues vector -!--------------------------------------------------------------------------------------------------- -pure function ho2ro(ho) result(ro) - - implicit none - real(pReal), intent(in), dimension(3) :: ho - real(pReal), dimension(4) :: ro - - ro = ax2ro(ho2ax(ho)) - -end function ho2ro - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief convert homochoric to unit quaternion -!--------------------------------------------------------------------------------------------------- -pure function ho2qu(ho) result(qu) - - implicit none - real(pReal), intent(in), dimension(3) :: ho - type(quaternion) :: qu - - qu = ax2qu(ho2ax(ho)) - -end function ho2qu - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief convert Euler angles to cubochoric -!--------------------------------------------------------------------------------------------------- -function eu2cu(eu) result(cu) - - implicit none - real(pReal), intent(in), dimension(3) :: eu - real(pReal), dimension(3) :: cu - - cu = ho2cu(eu2ho(eu)) - -end function eu2cu +end function om2ho !--------------------------------------------------------------------------------------------------- @@ -1074,7 +522,6 @@ end function eu2cu !--------------------------------------------------------------------------------------------------- function om2cu(om) result(cu) - implicit none real(pReal), intent(in), dimension(3,3) :: om real(pReal), dimension(3) :: cu @@ -1083,13 +530,279 @@ function om2cu(om) result(cu) end function om2cu +!--------------------------------------------------------------------------------------------------- +!> @author Marc De Graef, Carnegie Mellon University +!> @brief Euler angles to unit quaternion +!--------------------------------------------------------------------------------------------------- +pure function eu2qu(eu) result(qu) + + real(pReal), intent(in), dimension(3) :: eu + type(quaternion) :: qu + real(pReal), dimension(3) :: ee + real(pReal) :: cPhi, sPhi + + ee = 0.5_pReal*eu + + cPhi = cos(ee(2)) + sPhi = sin(ee(2)) + + qu = quaternion([ cPhi*cos(ee(1)+ee(3)), & + -P*sPhi*cos(ee(1)-ee(3)), & + -P*sPhi*sin(ee(1)-ee(3)), & + -P*cPhi*sin(ee(1)+ee(3))]) + if(qu%w < 0.0_pReal) qu = qu%homomorphed() + +end function eu2qu + + +!--------------------------------------------------------------------------------------------------- +!> @author Marc De Graef, Carnegie Mellon University +!> @brief Euler angles to orientation matrix +!--------------------------------------------------------------------------------------------------- +pure function eu2om(eu) result(om) + use prec, only: & + dEq0 + + real(pReal), intent(in), dimension(3) :: eu + real(pReal), dimension(3,3) :: om + + 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) + om(2,1) = -c(1)*s(3)-s(1)*c(3)*c(2) + om(2,2) = -s(1)*s(3)+c(1)*c(3)*c(2) + om(2,3) = c(3)*s(2) + 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 + + +!--------------------------------------------------------------------------------------------------- +!> @author Marc De Graef, Carnegie Mellon University +!> @brief convert euler to axis angle +!--------------------------------------------------------------------------------------------------- +pure function eu2ax(eu) result(ax) + use prec, only: & + dEq0, & + dEq + use math, only: & + PI + + 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) + sigma = 0.5*(eu(1)+eu(3)) + delta = 0.5*(eu(1)-eu(3)) + tau = sqrt(t**2+sin(sigma)**2) + + alpha = merge(PI, 2.0*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 + ax(1:3) = -P/tau * [ t*cos(delta), t*sin(delta), sin(sigma) ] ! passive axis-angle pair so a minus sign in front + ax(4) = alpha + if (alpha < 0.0) ax = -ax ! ensure alpha is positive + end if + +end function eu2ax + + +!--------------------------------------------------------------------------------------------------- +!> @author Marc De Graef, Carnegie Mellon University +!> @brief Euler angles to Rodrigues vector +!--------------------------------------------------------------------------------------------------- +pure function eu2ro(eu) result(ro) + use prec, only: & + dEq0 + use, intrinsic :: IEEE_ARITHMETIC, only: & + IEEE_value, & + IEEE_positive_inf + use math, only: & + PI + + 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) + elseif(dEq0(ro(4))) then + ro = [ 0.0_pReal, 0.0_pReal, P, 0.0_pReal ] + else + ro(4) = tan(ro(4)*0.5) + end if + +end function eu2ro + + +!--------------------------------------------------------------------------------------------------- +!> @author Marc De Graef, Carnegie Mellon University +!> @brief convert Euler angles to homochoric +!--------------------------------------------------------------------------------------------------- +pure function eu2ho(eu) result(ho) + + real(pReal), intent(in), dimension(3) :: eu + real(pReal), dimension(3) :: ho + + ho = ax2ho(eu2ax(eu)) + +end function eu2ho + + +!--------------------------------------------------------------------------------------------------- +!> @author Marc De Graef, Carnegie Mellon University +!> @brief convert Euler angles to cubochoric +!--------------------------------------------------------------------------------------------------- +function eu2cu(eu) result(cu) + + real(pReal), intent(in), dimension(3) :: eu + real(pReal), dimension(3) :: cu + + cu = ho2cu(eu2ho(eu)) + +end function eu2cu + + +!--------------------------------------------------------------------------------------------------- +!> @author Marc De Graef, Carnegie Mellon University +!> @brief convert axis angle pair to quaternion +!--------------------------------------------------------------------------------------------------- +pure function ax2qu(ax) result(qu) + use prec, only: & + dEq0 + + real(pReal), intent(in), dimension(4) :: ax + type(quaternion) :: qu + + real(pReal) :: c, s + + + if (dEq0(ax(4))) then + qu = quaternion([ 1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal ]) + else + c = cos(ax(4)*0.5) + s = sin(ax(4)*0.5) + qu = quaternion([ c, ax(1)*s, ax(2)*s, ax(3)*s ]) + end if + +end function ax2qu + + +!--------------------------------------------------------------------------------------------------- +!> @author Marc De Graef, Carnegie Mellon University +!> @brief convert axis angle pair to orientation matrix +!--------------------------------------------------------------------------------------------------- +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 + integer :: i + + c = cos(ax(4)) + s = sin(ax(4)) + omc = 1.0-c + + forall(i=1:3) om(i,i) = ax(i)**2*omc + c + + 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) + + if (P > 0.0) om = transpose(om) + +end function ax2om + + +!--------------------------------------------------------------------------------------------------- +!> @author Marc De Graef, Carnegie Mellon University +!> @brief convert axis angle pair to Euler angles +!--------------------------------------------------------------------------------------------------- +pure function ax2eu(ax) result(eu) + + real(pReal), intent(in), dimension(4) :: ax + real(pReal), dimension(3) :: eu + + eu = om2eu(ax2om(ax)) + +end function ax2eu + + +!--------------------------------------------------------------------------------------------------- +!> @author Marc De Graef, Carnegie Mellon University +!> @brief convert axis angle pair to Rodrigues vector +!--------------------------------------------------------------------------------------------------- +pure function ax2ro(ax) result(ro) + use, intrinsic :: IEEE_ARITHMETIC, only: & + IEEE_value, & + IEEE_positive_inf + use prec, only: & + dEq0 + use math, only: & + PI + + real(pReal), intent(in), dimension(4) :: ax + real(pReal), dimension(4) :: ro + + real(pReal), parameter :: thr = 1.0E-7 + + if (dEq0(ax(4))) then + ro = [ 0.0_pReal, 0.0_pReal, P, 0.0_pReal ] + 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 ),abs(ax(4)-PI) < thr) + end if + +end function ax2ro + + +!--------------------------------------------------------------------------------------------------- +!> @author Marc De Graef, Carnegie Mellon University +!> @brief convert axis angle pair to homochoric +!--------------------------------------------------------------------------------------------------- +pure function ax2ho(ax) result(ho) + + real(pReal), intent(in), dimension(4) :: ax + real(pReal), dimension(3) :: ho + + real(pReal) :: f + + f = 0.75 * ( ax(4) - sin(ax(4)) ) + f = f**(1.0/3.0) + ho = ax(1:3) * f + +end function ax2ho + + !--------------------------------------------------------------------------------------------------- !> @author Marc De Graef, Carnegie Mellon University !> @brief convert axis angle pair to cubochoric !--------------------------------------------------------------------------------------------------- function ax2cu(ax) result(cu) - implicit none real(pReal), intent(in), dimension(4) :: ax real(pReal), dimension(3) :: cu @@ -1098,13 +811,113 @@ function ax2cu(ax) result(cu) end function ax2cu +!--------------------------------------------------------------------------------------------------- +!> @author Marc De Graef, Carnegie Mellon University +!> @brief convert Rodrigues vector to unit quaternion +!--------------------------------------------------------------------------------------------------- +pure function ro2qu(ro) result(qu) + + real(pReal), intent(in), dimension(4) :: ro + type(quaternion) :: qu + + qu = ax2qu(ro2ax(ro)) + +end function ro2qu + + +!--------------------------------------------------------------------------------------------------- +!> @author Marc De Graef, Carnegie Mellon University +!> @brief convert Rodrigues vector to rotation matrix +!--------------------------------------------------------------------------------------------------- +pure function ro2om(ro) result(om) + + real(pReal), intent(in), dimension(4) :: ro + real(pReal), dimension(3,3) :: om + + om = ax2om(ro2ax(ro)) + +end function ro2om + + +!--------------------------------------------------------------------------------------------------- +!> @author Marc De Graef, Carnegie Mellon University +!> @brief convert Rodrigues vector to Euler angles +!--------------------------------------------------------------------------------------------------- +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 + + +!--------------------------------------------------------------------------------------------------- +!> @author Marc De Graef, Carnegie Mellon University +!> @brief convert Rodrigues vector to axis angle pair +!--------------------------------------------------------------------------------------------------- +pure function ro2ax(ro) result(ax) + use, intrinsic :: IEEE_ARITHMETIC, only: & + IEEE_is_finite + use prec, only: & + dEq0 + use math, only: & + PI + + real(pReal), intent(in), dimension(4) :: ro + real(pReal), dimension(4) :: ax + + real(pReal) :: ta, angle + + ta = ro(4) + + if (dEq0(ta)) then + ax = [ 0.0, 0.0, 1.0, 0.0 ] + elseif (.not. IEEE_is_finite(ta)) then + ax = [ ro(1), ro(2), ro(3), PI ] + else + angle = 2.0*atan(ta) + ta = 1.0/norm2(ro(1:3)) + ax = [ ro(1)/ta, ro(2)/ta, ro(3)/ta, angle ] + end if + +end function ro2ax + + +!--------------------------------------------------------------------------------------------------- +!> @author Marc De Graef, Carnegie Mellon University +!> @brief convert Rodrigues vector to homochoric +!--------------------------------------------------------------------------------------------------- +pure function ro2ho(ro) result(ho) + use, intrinsic :: IEEE_ARITHMETIC, only: & + IEEE_is_finite + use prec, only: & + dEq0 + use math, only: & + PI + + 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, 0.0, 0.0 ] + else + f = merge(2.0*atan(ro(4)) - sin(2.0*atan(ro(4))),PI, IEEE_is_finite(ro(4))) + ho = ro(1:3) * (0.75_pReal*f)**(1.0/3.0) + end if + +end function ro2ho + + !--------------------------------------------------------------------------------------------------- !> @author Marc De Graef, Carnegie Mellon University !> @brief convert Rodrigues vector to cubochoric !--------------------------------------------------------------------------------------------------- function ro2cu(ro) result(cu) - implicit none real(pReal), intent(in), dimension(4) :: ro real(pReal), dimension(3) :: cu @@ -1115,32 +928,130 @@ end function ro2cu !--------------------------------------------------------------------------------------------------- !> @author Marc De Graef, Carnegie Mellon University -!> @brief convert unit quaternion to cubochoric +!> @brief convert homochoric to unit quaternion !--------------------------------------------------------------------------------------------------- -function qu2cu(qu) result(cu) - - implicit none - type(quaternion), intent(in) :: qu - real(pReal), dimension(3) :: cu +pure function ho2qu(ho) result(qu) - cu = ho2cu(qu2ho(qu)) + real(pReal), intent(in), dimension(3) :: ho + type(quaternion) :: qu -end function qu2cu + qu = ax2qu(ho2ax(ho)) + +end function ho2qu !--------------------------------------------------------------------------------------------------- !> @author Marc De Graef, Carnegie Mellon University -!> @brief convert cubochoric to Euler angles +!> @brief convert homochoric to rotation matrix !--------------------------------------------------------------------------------------------------- -function cu2eu(cu) result(eu) +pure function ho2om(ho) result(om) - implicit none - real(pReal), intent(in), dimension(3) :: cu + real(pReal), intent(in), dimension(3) :: ho + real(pReal), dimension(3,3) :: om + + om = ax2om(ho2ax(ho)) + +end function ho2om + + +!--------------------------------------------------------------------------------------------------- +!> @author Marc De Graef, Carnegie Mellon University +!> @brief convert homochoric to Euler angles +!--------------------------------------------------------------------------------------------------- +pure function ho2eu(ho) result(eu) + + real(pReal), intent(in), dimension(3) :: ho real(pReal), dimension(3) :: eu - eu = ho2eu(cu2ho(cu)) + eu = ax2eu(ho2ax(ho)) -end function cu2eu +end function ho2eu + + +!--------------------------------------------------------------------------------------------------- +!> @author Marc De Graef, Carnegie Mellon University +!> @brief convert homochoric to axis angle pair +!--------------------------------------------------------------------------------------------------- +pure function ho2ax(ho) result(ax) + use prec, only: & + dEq0 + + 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, & + +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 + hm = hm*hmag_squared + s = s + tfit(i) * hm + end do + ax = [ho/sqrt(hmag_squared), 2.0_pReal*acos(s)] + end if + +end function ho2ax + + +!--------------------------------------------------------------------------------------------------- +!> @author Marc De Graef, Carnegie Mellon University +!> @brief convert homochoric to Rodrigues vector +!--------------------------------------------------------------------------------------------------- +pure function ho2ro(ho) result(ro) + + real(pReal), intent(in), dimension(3) :: ho + real(pReal), dimension(4) :: ro + + ro = ax2ro(ho2ax(ho)) + +end function ho2ro + + +!--------------------------------------------------------------------------------------------------- +!> @author Marc De Graef, Carnegie Mellon University +!> @brief convert homochoric to cubochoric +!--------------------------------------------------------------------------------------------------- +function ho2cu(ho) result(cu) + use Lambert, only: & + LambertBallToCube + + real(pReal), intent(in), dimension(3) :: ho + real(pReal), dimension(3) :: cu + + cu = LambertBallToCube(ho) + +end function ho2cu + + +!--------------------------------------------------------------------------------------------------- +!> @author Marc De Graef, Carnegie Mellon University +!> @brief convert cubochoric to unit quaternion +!--------------------------------------------------------------------------------------------------- +function cu2qu(cu) result(qu) + + real(pReal), intent(in), dimension(3) :: cu + type(quaternion) :: qu + + qu = ho2qu(cu2ho(cu)) + +end function cu2qu !--------------------------------------------------------------------------------------------------- @@ -1149,7 +1060,6 @@ end function cu2eu !--------------------------------------------------------------------------------------------------- function cu2om(cu) result(om) - implicit none real(pReal), intent(in), dimension(3) :: cu real(pReal), dimension(3,3) :: om @@ -1158,13 +1068,26 @@ function cu2om(cu) result(om) end function cu2om +!--------------------------------------------------------------------------------------------------- +!> @author Marc De Graef, Carnegie Mellon University +!> @brief convert cubochoric to Euler angles +!--------------------------------------------------------------------------------------------------- +function cu2eu(cu) result(eu) + + real(pReal), intent(in), dimension(3) :: cu + real(pReal), dimension(3) :: eu + + eu = ho2eu(cu2ho(cu)) + +end function cu2eu + + !--------------------------------------------------------------------------------------------------- !> @author Marc De Graef, Carnegie Mellon University !> @brief convert cubochoric to axis angle pair !--------------------------------------------------------------------------------------------------- function cu2ax(cu) result(ax) - implicit none real(pReal), intent(in), dimension(3) :: cu real(pReal), dimension(4) :: ax @@ -1179,7 +1102,6 @@ end function cu2ax !--------------------------------------------------------------------------------------------------- function cu2ro(cu) result(ro) - implicit none real(pReal), intent(in), dimension(3) :: cu real(pReal), dimension(4) :: ro @@ -1190,16 +1112,18 @@ end function cu2ro !--------------------------------------------------------------------------------------------------- !> @author Marc De Graef, Carnegie Mellon University -!> @brief convert cubochoric to unit quaternion +!> @brief convert cubochoric to homochoric !--------------------------------------------------------------------------------------------------- -function cu2qu(cu) result(qu) +function cu2ho(cu) result(ho) + use Lambert, only: & + LambertCubeToBall - implicit none real(pReal), intent(in), dimension(3) :: cu - type(quaternion) :: qu + real(pReal), dimension(3) :: ho - qu = ho2qu(cu2ho(cu)) + ho = LambertCubeToBall(cu) + +end function cu2ho -end function cu2qu end module rotations