diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 7958db9b8..dda2ed1b6 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -433,6 +433,15 @@ J2_plasticBehavior: - master - release +Marc_elementLib: + stage: marc + script: + - module load $IntelMarc $HDF5Marc $MSC + - Marc_elementLib/test.py + except: + - master + - release + ################################################################################################### grid_all_example: stage: example diff --git a/PRIVATE b/PRIVATE index 65905b5cd..83cb0b6ce 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 65905b5cd467546873be8be2dd7d6217db0012be +Subproject commit 83cb0b6ceea2e26c554ad4c6fbe6aabaee6fa27f diff --git a/VERSION b/VERSION index c0a399148..526119953 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.3-1484-g93fca511 +v2.0.3-1520-g701a9b18 diff --git a/python/damask/dadf5.py b/python/damask/dadf5.py index 5c9ebe08e..0c3a6248a 100644 --- a/python/damask/dadf5.py +++ b/python/damask/dadf5.py @@ -40,8 +40,9 @@ class DADF5(): self.version_major = f.attrs['DADF5-major'] self.version_minor = f.attrs['DADF5-minor'] - if self.version_major != 0 or not 2 <= self.version_minor <= 5: - raise TypeError('Unsupported DADF5 version {} '.format(f.attrs['DADF5-version'])) + if self.version_major != 0 or not 2 <= self.version_minor <= 6: + raise TypeError('Unsupported DADF5 version {}.{} '.format(f.attrs['DADF5_version_major'], + f.attrs['DADF5_version_minor'])) self.structured = 'grid' in f['geometry'].attrs.keys() @@ -886,8 +887,27 @@ class DADF5(): vtk_geom = vtk.vtkUnstructuredGrid() vtk_geom.SetPoints(nodes) vtk_geom.Allocate(f['/geometry/T_c'].shape[0]) + + if self.version_major == 0 and self.version_minor <= 5: + vtk_type = vtk.VTK_HEXAHEDRON + n_nodes = 8 + else: + if f['/geometry/T_c'].attrs['VTK_TYPE'] == b'TRIANGLE': + vtk_type = vtk.VTK_TRIANGLE + n_nodes = 3 + elif f['/geometry/T_c'].attrs['VTK_TYPE'] == b'QUAD': + vtk_type = vtk.VTK_QUAD + n_nodes = 4 + elif f['/geometry/T_c'].attrs['VTK_TYPE'] == b'TETRA': # not tested + vtk_type = vtk.VTK_TETRA + n_nodes = 4 + elif f['/geometry/T_c'].attrs['VTK_TYPE'] == b'HEXAHEDRON': + vtk_type = vtk.VTK_HEXAHEDRON + n_nodes = 8 + for i in f['/geometry/T_c']: - vtk_geom.InsertNextCell(vtk.VTK_HEXAHEDRON,8,i-1) # not for all elements! + vtk_geom.InsertNextCell(n_nodes,vtk_type,i-1) + elif mode == 'Point': Points = vtk.vtkPoints() Vertices = vtk.vtkCellArray() @@ -936,7 +956,7 @@ class DADF5(): vtk_geom.GetCellData().AddArray(vtk_data[-1]) self.set_visible('materialpoints',materialpoints_backup) - + constituents_backup = self.visible['constituents'].copy() self.set_visible('constituents',False) for label in (labels if isinstance(labels,list) else [labels]): diff --git a/python/damask/orientation.py b/python/damask/orientation.py index a86ba9331..dc6a74504 100644 --- a/python/damask/orientation.py +++ b/python/damask/orientation.py @@ -3,9 +3,8 @@ import math import numpy as np from . import Lambert -from .quaternion import Quaternion -from .quaternion import P +P = -1 #################################################################################################### class Rotation: @@ -26,7 +25,8 @@ class Rotation: Convention 4: Euler angle triplets are implemented using the Bunge convention, with the angular ranges as [0, 2π],[0, π],[0, 2π]. Convention 5: The rotation angle ω is limited to the interval [0, π]. - Convention 6: P = -1 (as default). + Convention 6: the real part of a quaternion is positive, Re(q) > 0 + Convention 7: P = -1 (as default). Usage ----- @@ -49,10 +49,7 @@ class Rotation: Unit quaternion that follows the conventions. Use .fromQuaternion to perform a sanity check. """ - if isinstance(quaternion,Quaternion): - self.quaternion = quaternion.copy() - else: - self.quaternion = Quaternion(q=quaternion[0],p=quaternion[1:4]) + self.quaternion = quaternion.copy() def __copy__(self): """Copy.""" @@ -64,9 +61,9 @@ class Rotation: def __repr__(self): """Orientation displayed as unit quaternion, rotation matrix, and Bunge-Euler angles.""" return '\n'.join([ - '{}'.format(self.quaternion), - 'Matrix:\n{}'.format( '\n'.join(['\t'.join(list(map(str,self.asMatrix()[i,:]))) for i in range(3)]) ), - 'Bunge Eulers / deg: {}'.format('\t'.join(list(map(str,self.asEulers(degrees=True)))) ), + 'Quaternion: (real={:.3f}, imag=<{:+.3f}, {:+.3f}, {:+.3f}>)'.format(*(self.quaternion)), + 'Matrix:\n{}'.format(self.asMatrix()), + 'Bunge Eulers / deg: ({:3.2f}, {:3.2f}, {:3.2f})'.format(*self.asEulers(degrees=True)), ]) def __mul__(self, other): @@ -85,19 +82,25 @@ class Rotation: """ 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 + self_q = self.quaternion[0] + self_p = self.quaternion[1:] + other_q = other.quaternion[0] + other_p = other.quaternion[1:] + R = self.__class__(np.append(self_q*other_q - np.dot(self_p,other_p), + self_q*other_p + other_q*self_p + P * np.cross(self_p,other_p))) + return R.standardize() + elif isinstance(other, (tuple,np.ndarray)): + if isinstance(other,tuple) or other.shape == (3,): # rotate a single (3)-vector or meshgrid + A = self.quaternion[0]**2.0 - np.dot(self.quaternion[1:],self.quaternion[1:]) + B = 2.0 * ( self.quaternion[1]*other[0] + + self.quaternion[2]*other[1] + + self.quaternion[3]*other[2]) + C = 2.0 * P*self.quaternion[0] 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), + A*other[0] + B*self.quaternion[1] + C*(self.quaternion[2]*other[2] - self.quaternion[3]*other[1]), + A*other[1] + B*self.quaternion[2] + C*(self.quaternion[3]*other[0] - self.quaternion[1]*other[2]), + A*other[2] + B*self.quaternion[3] + C*(self.quaternion[1]*other[1] - self.quaternion[2]*other[0]), ]) elif other.shape == (3,3,): # rotate a single (3x3)-matrix return np.dot(self.asMatrix(),np.dot(other,self.asMatrix().T)) @@ -105,25 +108,13 @@ class Rotation: 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() + self.quaternion[1:] *= -1 return self def inversed(self): @@ -133,7 +124,7 @@ class Rotation: def standardize(self): """In-place quaternion representation with positive q.""" - if self.quaternion.q < 0.0: self.quaternion.homomorph() + if self.quaternion[0] < 0.0: self.quaternion*=-1 return self def standardized(self): @@ -170,8 +161,7 @@ class Rotation: ################################################################################################ # convert to different orientation representations (numpy arrays) - def asQuaternion(self, - quaternion = False): + def asQuaternion(self): """ Unit quaternion [q, p_1, p_2, p_3] unless quaternion == True: damask.quaternion object. @@ -181,7 +171,7 @@ class Rotation: return quaternion as DAMASK object. """ - return self.quaternion if quaternion else self.quaternion.asArray() + return self.quaternion def asEulers(self, degrees = False): @@ -194,7 +184,7 @@ class Rotation: return angles in degrees. """ - eu = qu2eu(self.quaternion.asArray()) + eu = qu2eu(self.quaternion) if degrees: eu = np.degrees(eu) return eu @@ -212,13 +202,13 @@ class Rotation: return tuple of axis and angle. """ - ax = qu2ax(self.quaternion.asArray()) + ax = qu2ax(self.quaternion) if degrees: ax[3] = np.degrees(ax[3]) return (ax[:3],np.degrees(ax[3])) if pair else ax def asMatrix(self): """Rotation matrix.""" - return qu2om(self.quaternion.asArray()) + return qu2om(self.quaternion) def asRodrigues(self, vector = False): @@ -231,16 +221,16 @@ class Rotation: return as actual Rodrigues--Frank vector, i.e. rotation axis scaled by tan(ω/2). """ - ro = qu2ro(self.quaternion.asArray()) + ro = qu2ro(self.quaternion) return ro[:3]*ro[3] if vector else ro def asHomochoric(self): """Homochoric vector: (h_1, h_2, h_3).""" - return qu2ho(self.quaternion.asArray()) + return qu2ho(self.quaternion) def asCubochoric(self): """Cubochoric vector: (c_1, c_2, c_3).""" - return qu2cu(self.quaternion.asArray()) + return qu2cu(self.quaternion) def asM(self): """ @@ -252,19 +242,18 @@ class Rotation: https://doi.org/10.2514/1.28949 """ - return self.quaternion.asM() + return np.outer(self.quaternion,self.quaternion) ################################################################################################ # static constructors. The input data needs to follow the convention, options allow to # relax these convections - @classmethod - def fromQuaternion(cls, - quaternion, + @staticmethod + def fromQuaternion(quaternion, acceptHomomorph = False, P = -1): - qu = quaternion if isinstance(quaternion, np.ndarray) and quaternion.dtype == np.dtype(float) \ + qu = quaternion if isinstance(quaternion,np.ndarray) and quaternion.dtype == np.dtype(float) \ else np.array(quaternion,dtype=float) if P > 0: qu[1:4] *= -1 # convert from P=1 to P=-1 if qu[0] < 0.0: @@ -275,11 +264,10 @@ class Rotation: if not np.isclose(np.linalg.norm(qu), 1.0): raise ValueError('Quaternion is not of unit length.\n{} {} {} {}'.format(*qu)) - return cls(qu) + return Rotation(qu) - @classmethod - def fromEulers(cls, - eulers, + @staticmethod + def fromEulers(eulers, degrees = False): eu = eulers if isinstance(eulers, np.ndarray) and eulers.dtype == np.dtype(float) \ @@ -288,11 +276,10 @@ class Rotation: if np.any(eu < 0.0) or np.any(eu > 2.0*np.pi) or eu[1] > np.pi: raise ValueError('Euler angles outside of [0..2π],[0..π],[0..2π].\n{} {} {}.'.format(*eu)) - return cls(eu2qu(eu)) + return Rotation(eu2qu(eu)) - @classmethod - def fromAxisAngle(cls, - angleAxis, + @staticmethod + def fromAxisAngle(angleAxis, degrees = False, normalise = False, P = -1): @@ -307,11 +294,10 @@ class Rotation: if not np.isclose(np.linalg.norm(ax[0:3]), 1.0): raise ValueError('Axis angle rotation axis is not of unit length.\n{} {} {}'.format(*ax[0:3])) - return cls(ax2qu(ax)) + return Rotation(ax2qu(ax)) - @classmethod - def fromBasis(cls, - basis, + @staticmethod + def fromBasis(basis, orthonormal = True, reciprocal = False, ): @@ -330,18 +316,16 @@ class Rotation: or not np.isclose(np.dot(om[2],om[0]), 0.0): raise ValueError('matrix is not orthogonal.\n{}'.format(om)) - return cls(om2qu(om)) + return Rotation(om2qu(om)) - @classmethod - def fromMatrix(cls, - om, + @staticmethod + def fromMatrix(om, ): - return cls.fromBasis(om) + return Rotation.fromBasis(om) - @classmethod - def fromRodrigues(cls, - rodrigues, + @staticmethod + def fromRodrigues(rodrigues, normalise = False, P = -1): @@ -354,35 +338,32 @@ class Rotation: if ro[3] < 0.0: raise ValueError('Rodriques rotation angle not positive.\n'.format(ro[3])) - return cls(ro2qu(ro)) + return Rotation(ro2qu(ro)) - @classmethod - def fromHomochoric(cls, - homochoric, - P = -1): + @staticmethod + def fromHomochoric(homochoric, + P = -1): ho = homochoric if isinstance(homochoric, np.ndarray) and homochoric.dtype == np.dtype(float) \ else np.array(homochoric,dtype=float) if P > 0: ho *= -1 # convert from P=1 to P=-1 - return cls(ho2qu(ho)) + return Rotation(ho2qu(ho)) - @classmethod - def fromCubochoric(cls, - cubochoric, - P = -1): + @staticmethod + def fromCubochoric(cubochoric, + P = -1): cu = cubochoric if isinstance(cubochoric, np.ndarray) and cubochoric.dtype == np.dtype(float) \ else np.array(cubochoric,dtype=float) ho = cu2ho(cu) if P > 0: ho *= -1 # convert from P=1 to P=-1 - return cls(ho2qu(ho)) + return Rotation(ho2qu(ho)) - @classmethod - def fromAverage(cls, - rotations, + @staticmethod + def fromAverage(rotations, weights = []): """ Average rotation. @@ -412,19 +393,18 @@ class Rotation: 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) + return Rotation.fromQuaternion(np.real(vec.T[eig.argmax()]),acceptHomomorph = True) - @classmethod - def fromRandom(cls): + @staticmethod + def fromRandom(): 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) + return Rotation(np.array([np.cos(2.0*np.pi*r[0])*A, + np.sin(2.0*np.pi*r[1])*B, + np.cos(2.0*np.pi*r[1])*B, + np.sin(2.0*np.pi*r[0])*A])).standardize() @@ -1193,9 +1173,8 @@ class Orientation: return color - @classmethod - def fromAverage(cls, - orientations, + @staticmethod + def fromAverage(orientations, weights = []): """Create orientation from average of list of orientations.""" if not all(isinstance(item, Orientation) for item in orientations): diff --git a/python/damask/quaternion.py b/python/damask/quaternion.py deleted file mode 100644 index 203a398ef..000000000 --- a/python/damask/quaternion.py +++ /dev/null @@ -1,212 +0,0 @@ -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/python/tests/test_Rotation.py b/python/tests/test_Rotation.py index e352e1c26..012c95469 100644 --- a/python/tests/test_Rotation.py +++ b/python/tests/test_Rotation.py @@ -1,4 +1,5 @@ import os +from itertools import permutations import pytest import numpy as np @@ -6,6 +7,7 @@ import numpy as np import damask from damask import Rotation from damask import Orientation +from damask import Lattice n = 1000 @@ -45,19 +47,37 @@ class TestRotation: def test_Homochoric(self,default): for rot in default: assert np.allclose(rot.asRodrigues(), - Rotation.fromHomochoric(rot.asHomochoric()).asRodrigues(),rtol=5.e-5) + Rotation.fromHomochoric(rot.asHomochoric()).asRodrigues(),rtol=1.e-4) def test_Cubochoric(self,default): for rot in default: assert np.allclose(rot.asHomochoric(), - Rotation.fromCubochoric(rot.asCubochoric()).asHomochoric(),rtol=5.e-5) + Rotation.fromCubochoric(rot.asCubochoric()).asHomochoric()) def test_Quaternion(self,default): for rot in default: assert np.allclose(rot.asCubochoric(), - Rotation.fromQuaternion(rot.asQuaternion()).asCubochoric(),rtol=5.e-5) + Rotation.fromQuaternion(rot.asQuaternion()).asCubochoric()) + @pytest.mark.parametrize('color',[{'label':'red', 'RGB':[1,0,0],'direction':[0,0,1]}, + {'label':'green','RGB':[0,1,0],'direction':[0,1,1]}, + {'label':'blue', 'RGB':[0,0,1],'direction':[1,1,1]}]) + @pytest.mark.parametrize('lattice',['fcc','bcc']) + def test_IPF_cubic(self,default,color,lattice): + cube = damask.Orientation(damask.Rotation(),lattice) + for direction in set(permutations(np.array(color['direction']))): + assert np.allclose(cube.IPFcolor(direction),np.array(color['RGB'])) + + @pytest.mark.parametrize('lattice',Lattice.lattices) + def test_IPF(self,lattice): + direction = np.random.random(3)*2.0-1 + for rot in [Rotation.fromRandom() for r in range(n//100)]: + R = damask.Orientation(rot,lattice) + color = R.IPFcolor(direction) + for equivalent in R.equivalentOrientations(): + assert np.allclose(color,R.IPFcolor(direction)) + @pytest.mark.parametrize('model',['Bain','KS','GT','GT_prime','NW','Pitsch']) @pytest.mark.parametrize('lattice',['fcc','bcc']) def test_relationship_forward_backward(self,model,lattice): diff --git a/src/CPFEM.f90 b/src/CPFEM.f90 index d74155703..044837b74 100644 --- a/src/CPFEM.f90 +++ b/src/CPFEM.f90 @@ -377,6 +377,7 @@ subroutine CPFEM_results(inc,time) call constitutive_results call crystallite_results call homogenization_results + call discretization_results call results_finalizeIncrement call results_closeJobFile diff --git a/src/DAMASK_marc.f90 b/src/DAMASK_marc.f90 index a0a82500f..942c9afeb 100644 --- a/src/DAMASK_marc.f90 +++ b/src/DAMASK_marc.f90 @@ -350,6 +350,9 @@ subroutine flux(f,ts,n,time) !-------------------------------------------------------------------------------------------------- !> @brief trigger writing of results +!> @details uedinc is called before each new increment, not at the end of a converged one. +!> Therefore, storing the last written inc with an 'save' variable is required to avoid writing the +! same increment multiple times. !-------------------------------------------------------------------------------------------------- subroutine uedinc(inc,incsub) use prec @@ -357,8 +360,12 @@ subroutine uedinc(inc,incsub) implicit none integer, intent(in) :: inc, incsub + integer, save :: inc_written #include QUOTE(PASTE(./MarcInclude/creeps,Marc4DAMASK)) ! creeps is needed for timinc (time increment) - call CPFEM_results(inc,cptim) + if (inc > inc_written) then + call CPFEM_results(inc,cptim) + inc_written = inc + endif end subroutine uedinc diff --git a/src/IO.f90 b/src/IO.f90 index 0436e9b19..33924bb65 100644 --- a/src/IO.f90 +++ b/src/IO.f90 @@ -37,15 +37,12 @@ module IO #if defined(Marc4DAMASK) || defined(Abaqus) public :: & IO_open_inputFile, & - IO_countContinuousIntValues, & - IO_continuousIntValues, & + IO_continuousIntValues +#endif #if defined(Abaqus) + public :: & IO_extractValue, & IO_countDataLines -#elif defined(Marc4DAMASK) - IO_fixedNoEFloatValue, & - IO_fixedIntValue -#endif #endif contains @@ -57,7 +54,7 @@ contains !-------------------------------------------------------------------------------------------------- subroutine IO_init - write(6,'(/,a)') ' <<<+- IO init -+>>>' + write(6,'(/,a)') ' <<<+- IO init -+>>>'; flush(6) end subroutine IO_init @@ -444,57 +441,6 @@ integer function IO_intValue(string,chunkPos,myChunk) end function IO_intValue -#ifdef Marc4DAMASK -!-------------------------------------------------------------------------------------------------- -!> @brief reads float x.y+z value at myChunk from format string -!-------------------------------------------------------------------------------------------------- -real(pReal) function IO_fixedNoEFloatValue (string,ends,myChunk) - - character(len=*), intent(in) :: string !< raw input with known ends of each chunk - integer, intent(in) :: myChunk !< position number of desired chunk - integer, dimension(:), intent(in) :: ends !< positions of end of each tag/chunk in given string - character(len=*), parameter :: MYNAME = 'IO_fixedNoEFloatValue ' - character(len=*), parameter :: VALIDBASE = '0123456789.+-' - character(len=*), parameter :: VALIDEXP = '0123456789+-' - - real(pReal) :: base - integer :: expon - integer :: pos_exp - - pos_exp = scan(string(ends(myChunk)+1:ends(myChunk+1)),'+-',back=.true.) - hasExponent: if (pos_exp > 1) then - base = verifyFloatValue(trim(adjustl(string(ends(myChunk)+1:ends(myChunk)+pos_exp-1))),& - VALIDBASE,MYNAME//'(base): ') - expon = verifyIntValue(trim(adjustl(string(ends(myChunk)+pos_exp:ends(myChunk+1)))),& - VALIDEXP,MYNAME//'(exp): ') - else hasExponent - base = verifyFloatValue(trim(adjustl(string(ends(myChunk)+1:ends(myChunk+1)))),& - VALIDBASE,MYNAME//'(base): ') - expon = 0 - endif hasExponent - IO_fixedNoEFloatValue = base*10.0_pReal**real(expon,pReal) - -end function IO_fixedNoEFloatValue - - -!-------------------------------------------------------------------------------------------------- -!> @brief reads integer value at myChunk from fixed format string -!-------------------------------------------------------------------------------------------------- -integer function IO_fixedIntValue(string,ends,myChunk) - - character(len=*), intent(in) :: string !< raw input with known ends of each chunk - integer, intent(in) :: myChunk !< position number of desired chunk - integer, dimension(:), intent(in) :: ends !< positions of end of each tag/chunk in given string - character(len=*), parameter :: MYNAME = 'IO_fixedIntValue: ' - character(len=*), parameter :: VALIDCHARACTERS = '0123456789+-' - - IO_fixedIntValue = verifyIntValue(trim(adjustl(string(ends(myChunk)+1:ends(myChunk+1)))),& - VALIDCHARACTERS,MYNAME) - -end function IO_fixedIntValue -#endif - - !-------------------------------------------------------------------------------------------------- !> @brief changes characters in string to lower case !-------------------------------------------------------------------------------------------------- @@ -916,62 +862,6 @@ end function IO_countDataLines #endif -!-------------------------------------------------------------------------------------------------- -!> @brief count items in consecutive lines depending on lines -!> @details Marc: ints concatenated by "c" as last char or range of values a "to" b -!> Abaqus: triplet of start,stop,inc -!-------------------------------------------------------------------------------------------------- -integer function IO_countContinuousIntValues(fileUnit) - - integer, intent(in) :: fileUnit - -#ifdef Abaqus - integer :: l,c -#endif - integer, allocatable, dimension(:) :: chunkPos - character(len=pStringLen) :: line - - IO_countContinuousIntValues = 0 - line = '' - -#if defined(Marc4DAMASK) - do while (trim(line) /= IO_EOF) - line = IO_read(fileUnit) - chunkPos = IO_stringPos(line) - if (chunkPos(1) < 1) then ! empty line - exit - elseif (IO_lc(IO_stringValue(line,chunkPos,2)) == 'to' ) then ! found range indicator - IO_countContinuousIntValues = 1 + abs( IO_intValue(line,chunkPos,3) & - -IO_intValue(line,chunkPos,1)) - exit ! only one single range indicator allowed - else - IO_countContinuousIntValues = IO_countContinuousIntValues+chunkPos(1)-1 ! add line's count when assuming 'c' - if ( IO_lc(IO_stringValue(line,chunkPos,chunkPos(1))) /= 'c' ) then ! line finished, read last value - IO_countContinuousIntValues = IO_countContinuousIntValues+1 - exit ! data ended - endif - endif - enddo -#elif defined(Abaqus) - c = IO_countDataLines(fileUnit) - do l = 1,c - backspace(fileUnit) - enddo - - l = 1 - do while (trim(line) /= IO_EOF .and. l <= c) ! ToDo: is this correct? - l = l + 1 - line = IO_read(fileUnit) - chunkPos = IO_stringPos(line) - IO_countContinuousIntValues = IO_countContinuousIntValues + 1 + & ! assuming range generation - (IO_intValue(line,chunkPos,2)-IO_intValue(line,chunkPos,1))/& - max(1,IO_intValue(line,chunkPos,3)) - enddo -#endif - -end function IO_countContinuousIntValues - - !-------------------------------------------------------------------------------------------------- !> @brief return integer list corresponding to items in consecutive lines. !! First integer in array is counter diff --git a/src/discretization.f90 b/src/discretization.f90 index 5f9d3f521..ff0f59ef9 100644 --- a/src/discretization.f90 +++ b/src/discretization.f90 @@ -48,10 +48,10 @@ subroutine discretization_init(homogenizationAt,microstructureAt,& real(pReal), dimension(:,:), intent(in) :: & IPcoords0, & NodeCoords0 - integer, optional, intent(in) :: & + integer, optional, intent(in) :: & sharedNodesBeginn - write(6,'(/,a)') ' <<<+- discretization init -+>>>' + write(6,'(/,a)') ' <<<+- discretization init -+>>>'; flush(6) discretization_nElem = size(microstructureAt,1) discretization_nIP = size(IPcoords0,2)/discretization_nElem @@ -81,15 +81,15 @@ subroutine discretization_results real(pReal), dimension(:,:), allocatable :: u - call results_closeGroup(results_addGroup(trim('current/geometry'))) + call results_closeGroup(results_addGroup('current/geometry')) u = discretization_NodeCoords (1:3,:discretization_sharedNodesBeginn) & - discretization_NodeCoords0(1:3,:discretization_sharedNodesBeginn) - call results_writeDataset('current/geometry',u,'u_n','nodal displacements','m') + call results_writeDataset('current/geometry',u,'u_n','displacements of the nodes','m') u = discretization_IPcoords & - discretization_IPcoords0 - call results_writeDataset('current/geometry',u,'u_c','cell center displacements','m') + call results_writeDataset('current/geometry',u,'u_p','displacements of the materialpoints','m') end subroutine discretization_results diff --git a/src/element.f90 b/src/element.f90 index 3a1e3f5a3..a38036f93 100644 --- a/src/element.f90 +++ b/src/element.f90 @@ -22,6 +22,8 @@ module element NcellnodesPerCell, & nIPs, & nIPneighbors + character(len=:), allocatable :: & + vtkType integer, dimension(:,:), allocatable :: & Cell, & !< intra-element (cell) nodes that constitute a cell IPneighbor, & @@ -36,10 +38,10 @@ module element end type tElement - integer, parameter, private :: & + integer, parameter :: & NELEMTYPE = 13 - integer, dimension(NELEMTYPE), parameter, private :: NNODE = & + integer, dimension(NELEMTYPE), parameter :: NNODE = & [ & 3, & ! 2D 3node 1ip 6, & ! 2D 6node 3ip @@ -57,7 +59,7 @@ module element 20 & ! 3D 20node 27ip ] !< number of nodes that constitute a specific type of element - integer, dimension(NELEMTYPE), parameter, public :: GEOMTYPE = & + integer, dimension(NELEMTYPE), parameter :: GEOMTYPE = & [ & 1, & 2, & @@ -74,7 +76,7 @@ module element 10 & ] !< geometry type of particular element type - integer, dimension(maxval(GEOMTYPE)), parameter, private :: NCELLNODE = & + integer, dimension(maxval(GEOMTYPE)), parameter :: NCELLNODE = & [ & 3, & 7, & @@ -88,7 +90,7 @@ module element 64 & ] !< number of cell nodes in a specific geometry type - integer, dimension(maxval(GEOMTYPE)), parameter, private :: NIP = & + integer, dimension(maxval(GEOMTYPE)), parameter :: NIP = & [ & 1, & 3, & @@ -102,7 +104,7 @@ module element 27 & ] !< number of IPs in a specific geometry type - integer, dimension(maxval(GEOMTYPE)), parameter, private :: CELLTYPE = & + integer, dimension(maxval(GEOMTYPE)), parameter :: CELLTYPE = & [ & 1, & ! 2D 3node 2, & ! 2D 4node @@ -116,7 +118,7 @@ module element 4 & ! 3D 8node ] !< cell type that is used by each geometry type - integer, dimension(maxval(CELLTYPE)), parameter, private :: NIPNEIGHBOR = & + integer, dimension(maxval(CELLTYPE)), parameter :: NIPNEIGHBOR = & [ & 3, & ! 2D 3node 4, & ! 2D 4node @@ -124,7 +126,7 @@ module element 6 & ! 3D 8node ] !< number of ip neighbors / cell faces in a specific cell type - integer, dimension(maxval(CELLTYPE)), parameter, private :: NCELLNODEPERCELLFACE = & + integer, dimension(maxval(CELLTYPE)), parameter :: NCELLNODEPERCELLFACE = & [ & 2, & ! 2D 3node 2, & ! 2D 4node @@ -132,7 +134,7 @@ module element 4 & ! 3D 8node ] !< number of cell nodes in a specific cell type - integer, dimension(maxval(CELLTYPE)), parameter, private :: NCELLNODEPERCELL = & + integer, dimension(maxval(CELLTYPE)), parameter :: NCELLNODEPERCELL = & [ & 3, & ! 2D 3node 4, & ! 2D 4node @@ -146,27 +148,39 @@ module element ! Positive integers denote an intra-element IP identifier. ! Negative integers denote the interface behind which the neighboring (extra-element) IP will be located. - integer, dimension(nIPneighbor(cellType(1)),nIP(1)), parameter, private :: IPneighbor1 = & + integer, dimension(NIPNEIGHBOR(CELLTYPE(1)),NIP(1)), parameter :: IPNEIGHBOR1 = & reshape([& -2,-3,-1 & - ],[nIPneighbor(cellType(1)),nIP(1)]) +#if !defined(__GFORTRAN__) + ],shape(IPNEIGHBOR1)) +#else + ],[NIPNEIGHBOR(CELLTYPE(1)),NIP(1)]) +#endif - integer, dimension(nIPneighbor(cellType(2)),nIP(2)), parameter, private :: IPneighbor2 = & + integer, dimension(NIPNEIGHBOR(CELLTYPE(2)),NIP(2)), parameter :: IPNEIGHBOR2 = & reshape([& 2,-3, 3,-1, & -2, 1, 3,-1, & 2,-3,-2, 1 & - ],[nIPneighbor(cellType(2)),nIP(2)]) +#if !defined(__GFORTRAN__) + ],shape(IPNEIGHBOR2)) +#else + ],[NIPNEIGHBOR(CELLTYPE(2)),NIP(2)]) +#endif - integer, dimension(nIPneighbor(cellType(3)),nIP(3)), parameter, private :: IPneighbor3 = & + integer, dimension(NIPNEIGHBOR(CELLTYPE(3)),NIP(3)), parameter :: IPNEIGHBOR3 = & reshape([& 2,-4, 3,-1, & -2, 1, 4,-1, & 4,-4,-3, 1, & -2, 3,-3, 2 & - ],[nIPneighbor(cellType(3)),nIP(3)]) +#if !defined(__GFORTRAN__) + ],shape(IPNEIGHBOR3)) +#else + ],[NIPNEIGHBOR(CELLTYPE(3)),NIP(3)]) +#endif - integer, dimension(nIPneighbor(cellType(4)),nIP(4)), parameter, private :: IPneighbor4 = & + integer, dimension(NIPNEIGHBOR(CELLTYPE(4)),NIP(4)), parameter :: IPNEIGHBOR4 = & reshape([& 2,-4, 4,-1, & 3, 1, 5,-1, & @@ -177,22 +191,34 @@ module element 8,-4,-3, 4, & 9, 7,-3, 5, & -2, 8,-3, 6 & - ],[nIPneighbor(cellType(4)),nIP(4)]) +#if !defined(__GFORTRAN__) + ],shape(IPNEIGHBOR4)) +#else + ],[NIPNEIGHBOR(CELLTYPE(4)),NIP(4)]) +#endif - integer, dimension(nIPneighbor(cellType(5)),nIP(5)), parameter, private :: IPneighbor5 = & + integer, dimension(NIPNEIGHBOR(CELLTYPE(5)),NIP(5)), parameter :: IPNEIGHBOR5 = & reshape([& -1,-2,-3,-4 & - ],[nIPneighbor(cellType(5)),nIP(5)]) +#if !defined(__GFORTRAN__) + ],shape(IPNEIGHBOR5)) +#else + ],[NIPNEIGHBOR(CELLTYPE(5)),NIP(5)]) +#endif - integer, dimension(nIPneighbor(cellType(6)),nIP(6)), parameter, private :: IPneighbor6 = & + integer, dimension(NIPNEIGHBOR(CELLTYPE(6)),NIP(6)), parameter :: IPNEIGHBOR6 = & reshape([& 2,-4, 3,-2, 4,-1, & -2, 1, 3,-2, 4,-1, & 2,-4,-3, 1, 4,-1, & 2,-4, 3,-2,-3, 1 & - ],[nIPneighbor(cellType(6)),nIP(6)]) +#if !defined(__GFORTRAN__) + ],shape(IPNEIGHBOR6)) +#else + ],[NIPNEIGHBOR(CELLTYPE(6)),NIP(6)]) +#endif - integer, dimension(nIPneighbor(cellType(7)),nIP(7)), parameter, private :: IPneighbor7 = & + integer, dimension(NIPNEIGHBOR(CELLTYPE(7)),NIP(7)), parameter :: IPNEIGHBOR7 = & reshape([& 2,-4, 3,-2, 4,-1, & -3, 1, 3,-2, 5,-1, & @@ -200,14 +226,22 @@ module element 5,-4, 6,-2,-5, 1, & -3, 4, 6,-2,-5, 2, & 5,-4,-3, 4,-5, 3 & - ],[nIPneighbor(cellType(7)),nIP(7)]) +#if !defined(__GFORTRAN__) + ],shape(IPNEIGHBOR7)) +#else + ],[NIPNEIGHBOR(CELLTYPE(7)),NIP(7)]) +#endif - integer, dimension(nIPneighbor(cellType(8)),nIP(8)), parameter, private :: IPneighbor8 = & + integer, dimension(NIPNEIGHBOR(CELLTYPE(8)),NIP(8)), parameter :: IPNEIGHBOR8 = & reshape([& - -3,-5,-4,-2,-6,-1 & - ],[nIPneighbor(cellType(8)),nIP(8)]) + -3,-5,-4,-2,-6,-1 & +#if !defined(__GFORTRAN__) + ],shape(IPNEIGHBOR8)) +#else + ],[NIPNEIGHBOR(CELLTYPE(8)),NIP(8)]) +#endif - integer, dimension(nIPneighbor(cellType(9)),nIP(9)), parameter, private :: IPneighbor9 = & + integer, dimension(NIPNEIGHBOR(CELLTYPE(9)),NIP(9)), parameter :: IPNEIGHBOR9 = & reshape([& 2,-5, 3,-2, 5,-1, & -3, 1, 4,-2, 6,-1, & @@ -217,9 +251,13 @@ module element -3, 5, 8,-2,-6, 2, & 8,-5,-4, 5,-6, 3, & -3, 7,-4, 6,-6, 4 & - ],[nIPneighbor(cellType(9)),nIP(9)]) +#if !defined(__GFORTRAN__) + ],shape(IPNEIGHBOR9)) +#else + ],[NIPNEIGHBOR(CELLTYPE(9)),NIP(9)]) +#endif - integer, dimension(nIPneighbor(cellType(10)),nIP(10)), parameter, private :: IPneighbor10 = & + integer, dimension(NIPNEIGHBOR(CELLTYPE(10)),NIP(10)), parameter :: IPNEIGHBOR10 = & reshape([& 2,-5, 4,-2,10,-1, & 3, 1, 5,-2,11,-1, & @@ -248,17 +286,25 @@ module element 26,-5,-4,22,-6,16, & 27,25,-4,23,-6,17, & -3,26,-4,24,-6,18 & - ],[nIPneighbor(cellType(10)),nIP(10)]) +#if !defined(__GFORTRAN__) + ],shape(IPNEIGHBOR10)) +#else + ],[NIPNEIGHBOR(CELLTYPE(10)),NIP(10)]) +#endif - integer, dimension(nNode(1),NcellNode(geomType(1))), parameter :: cellNodeParentNodeWeights1 = & + integer, dimension(NNODE(1),NCELLNODE(GEOMTYPE(1))), parameter :: CELLNODEPARENTNODEWEIGHTS1 = & reshape([& 1, 0, 0, & 0, 1, 0, & 0, 0, 1 & - ],[nNode(1),NcellNode(geomType(1))]) !< 2D 3node 1ip +#if !defined(__GFORTRAN__) + ],shape(CELLNODEPARENTNODEWEIGHTS1)) !< 2D 3node 1ip +#else + ],[NNODE(1),NCELLNODE(GEOMTYPE(1))]) +#endif - integer, dimension(nNode(2),NcellNode(geomType(2))), parameter :: cellNodeParentNodeWeights2 = & + integer, dimension(NNODE(2),NCELLNODE(GEOMTYPE(2))), parameter :: CELLNODEPARENTNODEWEIGHTS2 = & reshape([& 1, 0, 0, 0, 0, 0, & 0, 1, 0, 0, 0, 0, & @@ -267,9 +313,13 @@ module element 0, 0, 0, 0, 1, 0, & 0, 0, 0, 0, 0, 1, & 1, 1, 1, 2, 2, 2 & - ],[nNode(2),NcellNode(geomType(2))]) !< 2D 6node 3ip +#if !defined(__GFORTRAN__) + ],shape(CELLNODEPARENTNODEWEIGHTS2)) !< 2D 6node 3ip +#else + ],[NNODE(2),NCELLNODE(GEOMTYPE(2))]) +#endif - integer, dimension(nNode(3),NcellNode(geomType(3))), parameter :: cellNodeParentNodeWeights3 = & + integer, dimension(NNODE(3),NCELLNODE(GEOMTYPE(3))), parameter :: CELLNODEPARENTNODEWEIGHTS3 = & reshape([& 1, 0, 0, 0, & 0, 1, 0, 0, & @@ -280,9 +330,13 @@ module element 0, 0, 1, 1, & 1, 0, 0, 1, & 1, 1, 1, 1 & - ],[nNode(3),NcellNode(geomType(3))]) !< 2D 6node 3ip +#if !defined(__GFORTRAN__) + ],shape(CELLNODEPARENTNODEWEIGHTS3)) !< 2D 6node 3ip +#else + ],[NNODE(3),NCELLNODE(GEOMTYPE(3))]) +#endif - integer, dimension(nNode(4),NcellNode(geomType(4))), parameter :: cellNodeParentNodeWeights4 = & + integer, dimension(NNODE(4),NCELLNODE(GEOMTYPE(4))), parameter :: CELLNODEPARENTNODEWEIGHTS4 = & reshape([& 1, 0, 0, 0, 0, 0, 0, 0, & 0, 1, 0, 0, 0, 0, 0, 0, & @@ -300,9 +354,13 @@ module element 1, 4, 1, 1, 8, 8, 2, 2, & 1, 1, 4, 1, 2, 8, 8, 2, & 1, 1, 1, 4, 2, 2, 8, 8 & - ],[nNode(4),NcellNode(geomType(4))]) !< 2D 8node 9ip +#if !defined(__GFORTRAN__) + ],shape(CELLNODEPARENTNODEWEIGHTS4)) !< 2D 8node 9ip +#else + ],[NNODE(4),NCELLNODE(GEOMTYPE(4))]) +#endif - integer, dimension(nNode(5),NcellNode(geomType(5))), parameter :: cellNodeParentNodeWeights5 = & + integer, dimension(NNODE(5),NCELLNODE(GEOMTYPE(5))), parameter :: CELLNODEPARENTNODEWEIGHTS5 = & reshape([& 1, 0, 0, 0, 0, 0, 0, 0, & 0, 1, 0, 0, 0, 0, 0, 0, & @@ -313,17 +371,25 @@ module element 0, 0, 0, 0, 0, 0, 1, 0, & 0, 0, 0, 0, 0, 0, 0, 1, & 1, 1, 1, 1, 2, 2, 2, 2 & - ],[nNode(5),NcellNode(geomType(5))]) !< 2D 8node 4ip +#if !defined(__GFORTRAN__) + ],shape(CELLNODEPARENTNODEWEIGHTS5)) !< 2D 8node 4ip +#else + ],[NNODE(5),NCELLNODE(GEOMTYPE(5))]) +#endif - integer, dimension(nNode(6),NcellNode(geomType(6))), parameter :: cellNodeParentNodeWeights6 = & + integer, dimension(NNODE(6),NcellNode(GEOMTYPE(6))), parameter :: CELLNODEPARENTNODEWEIGHTS6 = & reshape([& 1, 0, 0, 0, & 0, 1, 0, 0, & 0, 0, 1, 0, & 0, 0, 0, 1 & - ],[nNode(6),NcellNode(geomType(6))]) !< 3D 4node 1ip +#if !defined(__GFORTRAN__) + ],shape(CELLNODEPARENTNODEWEIGHTS6)) !< 3D 4node 1ip +#else + ],[NNODE(6),NcellNode(GEOMTYPE(6))]) +#endif - integer, dimension(nNode(7),NcellNode(geomType(7))), parameter :: cellNodeParentNodeWeights7 = & + integer, dimension(NNODE(7),NCELLNODE(GEOMTYPE(7))), parameter :: CELLNODEPARENTNODEWEIGHTS7 = & reshape([& 1, 0, 0, 0, 0, & 0, 1, 0, 0, 0, & @@ -340,9 +406,13 @@ module element 0, 1, 1, 1, 0, & 1, 0, 1, 1, 0, & 0, 0, 0, 0, 1 & - ],[nNode(7),NcellNode(geomType(7))]) !< 3D 5node 4ip +#if !defined(__GFORTRAN__) + ],shape(CELLNODEPARENTNODEWEIGHTS7)) !< 3D 5node 4ip +#else + ],[NNODE(7),NCELLNODE(GEOMTYPE(7))]) +#endif - integer, dimension(nNode(8),NcellNode(geomType(8))), parameter :: cellNodeParentNodeWeights8 = & + integer, dimension(NNODE(8),NCELLNODE(GEOMTYPE(8))), parameter :: CELLNODEPARENTNODEWEIGHTS8 = & reshape([& 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, & 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, & @@ -359,9 +429,13 @@ module element 0, 1, 1, 1, 0, 2, 0, 0, 2, 2, & 1, 0, 1, 1, 0, 0, 2, 2, 0, 2, & 3, 3, 3, 3, 4, 4, 4, 4, 4, 4 & - ],[nNode(8),NcellNode(geomType(8))]) !< 3D 10node 4ip +#if !defined(__GFORTRAN__) + ],shape(CELLNODEPARENTNODEWEIGHTS8)) !< 3D 10node 4ip +#else + ],[NNODE(8),NCELLNODE(GEOMTYPE(8))]) +#endif - integer, dimension(nNode(9),NcellNode(geomType(9))), parameter :: cellNodeParentNodeWeights9 = & + integer, dimension(NNODE(9),NCELLNODE(GEOMTYPE(9))), parameter :: CELLNODEPARENTNODEWEIGHTS9 = & reshape([& 1, 0, 0, 0, 0, 0, & 0, 1, 0, 0, 0, 0, & @@ -384,9 +458,13 @@ module element 1, 0, 1, 1, 0, 1, & 0, 0, 0, 1, 1, 1, & 1, 1, 1, 1, 1, 1 & - ],[nNode(9),NcellNode(geomType(9))]) !< 3D 6node 6ip +#if !defined(__GFORTRAN__) + ],shape(CELLNODEPARENTNODEWEIGHTS9)) !< 3D 6node 6ip +#else + ],[NNODE(9),NCELLNODE(GEOMTYPE(9))]) +#endif - integer, dimension(nNode(10),NcellNode(geomType(10))), parameter :: cellNodeParentNodeWeights10 = & + integer, dimension(NNODE(10),NCELLNODE(GEOMTYPE(10))), parameter :: CELLNODEPARENTNODEWEIGHTS10 = & reshape([& 1, 0, 0, 0, 0, 0, 0, 0, & 0, 1, 0, 0, 0, 0, 0, 0, & @@ -396,9 +474,13 @@ module element 0, 0, 0, 0, 0, 1, 0, 0, & 0, 0, 0, 0, 0, 0, 1, 0, & 0, 0, 0, 0, 0, 0, 0, 1 & - ],[nNode(10),NcellNode(geomType(10))]) !< 3D 8node 1ip +#if !defined(__GFORTRAN__) + ],shape(CELLNODEPARENTNODEWEIGHTS10)) !< 3D 8node 1ip +#else + ],[NNODE(10),NCELLNODE(GEOMTYPE(10))]) +#endif - integer, dimension(nNode(11),NcellNode(geomType(11))), parameter :: cellNodeParentNodeWeights11 = & + integer, dimension(NNODE(11),NCELLNODE(GEOMTYPE(11))), parameter :: CELLNODEPARENTNODEWEIGHTS11 = & reshape([& 1, 0, 0, 0, 0, 0, 0, 0, & ! 0, 1, 0, 0, 0, 0, 0, 0, & ! @@ -427,9 +509,13 @@ module element 1, 0, 0, 1, 1, 0, 0, 1, & ! 25 0, 0, 0, 0, 1, 1, 1, 1, & ! 1, 1, 1, 1, 1, 1, 1, 1 & ! - ],[nNode(11),NcellNode(geomType(11))]) !< 3D 8node 8ip +#if !defined(__GFORTRAN__) + ],shape(CELLNODEPARENTNODEWEIGHTS11)) !< 3D 8node 8ip +#else + ],[NNODE(11),NCELLNODE(GEOMTYPE(11))]) +#endif - integer, dimension(nNode(12),NcellNode(geomType(12))), parameter :: cellNodeParentNodeWeights12 = & + integer, dimension(NNODE(12),NCELLNODE(GEOMTYPE(12))), parameter :: CELLNODEPARENTNODEWEIGHTS12 = & reshape([& 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! @@ -458,9 +544,13 @@ module element 1, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 2, 0, 0, 0, 2, 2, 0, 0, 2, & ! 25 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 2, 2, 2, 2, 0, 0, 0, 0, & ! 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4 & ! - ],[nNode(12),NcellNode(geomType(12))]) !< 3D 20node 8ip +#if !defined(__GFORTRAN__) + ],shape(CELLNODEPARENTNODEWEIGHTS12)) !< 3D 20node 8ip +#else + ],[NNODE(12),NCELLNODE(GEOMTYPE(12))]) +#endif - integer, dimension(nNode(13),NcellNode(geomType(13))), parameter :: cellNodeParentNodeWeights13 = & + integer, dimension(NNODE(13),NCELLNODE(GEOMTYPE(13))), parameter :: CELLNODEPARENTNODEWEIGHTS13 = & reshape([& 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! @@ -526,30 +616,46 @@ module element 4, 8, 4, 3, 8,24, 8, 4, 12,12, 4, 4, 32,32,12,12, 12,32,12, 4, & ! 3, 4, 8, 4, 4, 8,24, 8, 4,12,12, 4, 12,32,32,12, 4,12,32,12, & ! 4, 3, 4, 8, 8, 4, 8,24, 4, 4,12,12, 12,12,32,32, 12, 4,12,32 & ! - ],[nNode(13),NcellNode(geomType(13))]) !< 3D 20node 27ip +#if !defined(__GFORTRAN__) + ],shape(CELLNODEPARENTNODEWEIGHTS13)) !< 3D 20node 27ip +#else + ],[NNODE(13),NCELLNODE(GEOMTYPE(13))]) +#endif - integer, dimension(NCELLNODEPERCELL(CELLTYPE(1)),NIP(1)), parameter :: CELL1 = & + integer, dimension(NCELLNODEPERCELL(CELLTYPE(1)),NIP(1)), parameter :: CELL1 = & reshape([& 1,2,3 & +#if !defined(__GFORTRAN__) + ],shape(CELL1)) +#else ],[NCELLNODEPERCELL(CELLTYPE(1)),NIP(1)]) +#endif - integer, dimension(NCELLNODEPERCELL(CELLTYPE(2)),NIP(2)), parameter :: CELL2 = & + integer, dimension(NCELLNODEPERCELL(CELLTYPE(2)),NIP(2)), parameter :: CELL2 = & reshape([& 1, 4, 7, 6, & 2, 5, 7, 4, & 3, 6, 7, 5 & +#if !defined(__GFORTRAN__) + ],shape(CELL2)) +#else ],[NCELLNODEPERCELL(CELLTYPE(2)),NIP(2)]) +#endif - integer, dimension(NCELLNODEPERCELL(CELLTYPE(3)),NIP(3)), parameter :: CELL3 = & + integer, dimension(NCELLNODEPERCELL(CELLTYPE(3)),NIP(3)), parameter :: CELL3 = & reshape([& 1, 5, 9, 8, & 5, 2, 6, 9, & 8, 9, 7, 4, & 9, 6, 3, 7 & +#if !defined(__GFORTRAN__) + ],shape(CELL3)) +#else ],[NCELLNODEPERCELL(CELLTYPE(3)),NIP(3)]) +#endif - integer, dimension(NCELLNODEPERCELL(CELLTYPE(4)),NIP(4)), parameter :: CELL4 = & + integer, dimension(NCELLNODEPERCELL(CELLTYPE(4)),NIP(4)), parameter :: CELL4 = & reshape([& 1, 5,13,12, & 5, 6,14,13, & @@ -560,22 +666,34 @@ module element 11,16,10, 4, & 16,15, 9,10, & 15, 8, 3, 9 & +#if !defined(__GFORTRAN__) + ],shape(CELL4)) +#else ],[NCELLNODEPERCELL(CELLTYPE(4)),NIP(4)]) +#endif - integer, dimension(NCELLNODEPERCELL(CELLTYPE(5)),NIP(5)), parameter :: CELL5 = & + integer, dimension(NCELLNODEPERCELL(CELLTYPE(5)),NIP(5)), parameter :: CELL5 = & reshape([& 1, 2, 3, 4 & +#if !defined(__GFORTRAN__) + ],shape(CELL5)) +#else ],[NCELLNODEPERCELL(CELLTYPE(5)),NIP(5)]) +#endif - integer, dimension(NCELLNODEPERCELL(CELLTYPE(6)),NIP(6)), parameter :: CELL6 = & + integer, dimension(NCELLNODEPERCELL(CELLTYPE(6)),NIP(6)), parameter :: CELL6 = & reshape([& 1, 5,11, 7, 8,12,15,14, & 5, 2, 6,11,12, 9,13,15, & 7,11, 6, 3,14,15,13,10, & 8,12,15, 4, 4, 9,13,10 & +#if !defined(__GFORTRAN__) + ],shape(CELL6)) +#else ],[NCELLNODEPERCELL(CELLTYPE(6)),NIP(6)]) +#endif - integer, dimension(NCELLNODEPERCELL(CELLTYPE(7)),NIP(7)), parameter :: CELL7 = & + integer, dimension(NCELLNODEPERCELL(CELLTYPE(7)),NIP(7)), parameter :: CELL7 = & reshape([& 1, 7,16, 9,10,17,21,19, & 7, 2, 8,16,17,11,18,21, & @@ -583,14 +701,22 @@ module element 10,17,21,19, 4,13,20,15, & 17,11,18,21,13, 5,14,20, & 19,21,18,12,15,20,14, 6 & +#if !defined(__GFORTRAN__) + ],shape(CELL7)) +#else ],[NCELLNODEPERCELL(CELLTYPE(7)),NIP(7)]) +#endif - integer, dimension(NCELLNODEPERCELL(CELLTYPE(8)),NIP(8)), parameter :: CELL8 = & + integer, dimension(NCELLNODEPERCELL(CELLTYPE(8)),NIP(8)), parameter :: CELL8 = & reshape([& 1, 2, 3, 4, 5, 6, 7, 8 & +#if !defined(__GFORTRAN__) + ],shape(CELL8)) +#else ],[NCELLNODEPERCELL(CELLTYPE(8)),NIP(8)]) +#endif - integer, dimension(NCELLNODEPERCELL(CELLTYPE(9)),NIP(9)), parameter :: CELL9 = & + integer, dimension(NCELLNODEPERCELL(CELLTYPE(9)),NIP(9)), parameter :: CELL9 = & reshape([& 1, 9,21,12,13,22,27,25, & 9, 2,10,21,22,14,23,27, & @@ -600,7 +726,11 @@ module element 22,14,23,27,17, 6,18,26, & 25,27,24,16,20,26,19, 8, & 27,23,15,24,26,18, 7,19 & +#if !defined(__GFORTRAN__) + ],shape(CELL9)) +#else ],[NCELLNODEPERCELL(CELLTYPE(9)),NIP(9)]) +#endif integer, dimension(NCELLNODEPERCELL(CELLTYPE(10)),NIP(10)), parameter :: CELL10 = & reshape([& @@ -631,7 +761,11 @@ module element 51,64,50,24,31,56,30, 8, & 64,63,49,50,56,55,29,30, & 63,48,23,49,55,28, 7,29 & +#if !defined(__GFORTRAN__) + ],shape(CELL10)) +#else ],[NCELLNODEPERCELL(CELLTYPE(10)),NIP(10)]) +#endif integer, dimension(NCELLNODEPERCELLFACE(1),NIPNEIGHBOR(1)), parameter :: CELLFACE1 = & @@ -639,7 +773,11 @@ module element 2,3, & 3,1, & 1,2 & - ],[NCELLNODEPERCELLFACE(1),NIPNEIGHBOR(1)]) !< 2D 3node, VTK_TRIANGLE (5) +#if !defined(__GFORTRAN__) + ],shape(CELLFACE1)) !< 2D 3node, VTK_TRIANGLE (5) +#else + ],[NCELLNODEPERCELLFACE(1),NIPNEIGHBOR(1)]) +#endif integer, dimension(NCELLNODEPERCELLFACE(2),NIPNEIGHBOR(2)), parameter :: CELLFACE2 = & reshape([& @@ -647,7 +785,11 @@ module element 4,1, & 3,4, & 1,2 & - ],[NCELLNODEPERCELLFACE(2),NIPNEIGHBOR(2)]) !< 2D 4node, VTK_QUAD (9) +#if !defined(__GFORTRAN__) + ],shape(CELLFACE2)) !< 2D 4node, VTK_QUAD (9) +#else + ],[NCELLNODEPERCELLFACE(2),NIPNEIGHBOR(2)]) +#endif integer, dimension(NCELLNODEPERCELLFACE(3),NIPNEIGHBOR(3)), parameter :: CELLFACE3 = & reshape([& @@ -655,7 +797,11 @@ module element 1,2,4, & 2,3,4, & 1,4,3 & - ],[NCELLNODEPERCELLFACE(3),NIPNEIGHBOR(3)]) !< 3D 4node, VTK_TETRA (10) +#if !defined(__GFORTRAN__) + ],shape(CELLFACE3)) !< 3D 4node, VTK_TETRA (10) +#else + ],[NCELLNODEPERCELLFACE(3),NIPNEIGHBOR(3)]) +#endif integer, dimension(NCELLNODEPERCELLFACE(4),NIPNEIGHBOR(4)), parameter :: CELLFACE4 = & reshape([& @@ -665,114 +811,128 @@ module element 1,2,6,5, & 5,6,7,8, & 1,4,3,2 & - ],[NCELLNODEPERCELLFACE(4),NIPNEIGHBOR(4)]) !< 3D 8node, VTK_HEXAHEDRON (12) +#if !defined(__GFORTRAN__) + ],shape(CELLFACE4)) !< 3D 8node, VTK_HEXAHEDRON (12) +#else + ],[NCELLNODEPERCELLFACE(4),NIPNEIGHBOR(4)]) +#endif contains + +!--------------------------------------------------------------------------------------------------- +!> define properties of an element +!--------------------------------------------------------------------------------------------------- subroutine tElement_init(self,elemType) - class(tElement) :: self - integer, intent(in) :: elemType - self%elemType = elemType - - self%Nnodes = Nnode (self%elemType) - self%geomType = geomType (self%elemType) - select case (self%elemType) - case(1) - self%cellNodeParentNodeWeights = cellNodeParentNodeWeights1 - case(2) - self%cellNodeParentNodeWeights = cellNodeParentNodeWeights2 - case(3) - self%cellNodeParentNodeWeights = cellNodeParentNodeWeights3 - case(4) - self%cellNodeParentNodeWeights = cellNodeParentNodeWeights4 - case(5) - self%cellNodeParentNodeWeights = cellNodeParentNodeWeights5 - case(6) - self%cellNodeParentNodeWeights = cellNodeParentNodeWeights6 - case(7) - self%cellNodeParentNodeWeights = cellNodeParentNodeWeights7 - case(8) - self%cellNodeParentNodeWeights = cellNodeParentNodeWeights8 - case(9) - self%cellNodeParentNodeWeights = cellNodeParentNodeWeights9 - case(10) - self%cellNodeParentNodeWeights = cellNodeParentNodeWeights10 - case(11) - self%cellNodeParentNodeWeights = cellNodeParentNodeWeights11 - case(12) - self%cellNodeParentNodeWeights = cellNodeParentNodeWeights12 - case(13) - self%cellNodeParentNodeWeights = cellNodeParentNodeWeights13 - case default - call IO_error(0,ext_msg='invalid element type') - end select - - - self%NcellNodes = NcellNode (self%geomType) - self%nIPs = nIP (self%geomType) - self%cellType = cellType (self%geomType) - - select case (self%geomType) - case(1) - self%IPneighbor = IPneighbor1 - self%cell = CELL1 - case(2) - self%IPneighbor = IPneighbor2 - self%cell = CELL2 - case(3) - self%IPneighbor = IPneighbor3 - self%cell = CELL3 - case(4) - self%IPneighbor = IPneighbor4 - self%cell = CELL4 - case(5) - self%IPneighbor = IPneighbor5 - self%cell = CELL5 - case(6) - self%IPneighbor = IPneighbor6 - self%cell = CELL6 - case(7) - self%IPneighbor = IPneighbor7 - self%cell = CELL7 - case(8) - self%IPneighbor = IPneighbor8 - self%cell = CELL8 - case(9) - self%IPneighbor = IPneighbor9 - self%cell = CELL9 - case(10) - self%IPneighbor = IPneighbor10 - self%cell = CELL10 - end select - - self%NcellnodesPerCell = NCELLNODEPERCELL(self%cellType) + class(tElement) :: self + integer, intent(in) :: elemType - select case(self%cellType) - case(1) - self%cellFace = CELLFACE1 - case(2) - self%cellFace = CELLFACE2 - case(3) - self%cellFace = CELLFACE3 - case(4) - self%cellFace = CELLFACE4 - end select - - self%nIPneighbors = size(self%IPneighbor,1) - - write(6,'(/,a)') ' <<<+- element_init -+>>>' - - write(6,*) ' element type: ',self%elemType - write(6,*) ' geom type: ',self%geomType - write(6,*) ' cell type: ',self%cellType - write(6,*) ' # node: ',self%Nnodes - write(6,*) ' # IP: ',self%nIPs - write(6,*) ' # cellnode: ',self%Ncellnodes - write(6,*) ' # cellnode/cell: ',self%NcellnodesPerCell - write(6,*) ' # IP neighbor: ',self%nIPneighbors - + self%elemType = elemType + + self%Nnodes = NNODE (self%elemType) + self%geomType = GEOMTYPE(self%elemType) + + select case (self%elemType) + case(1) + self%cellNodeParentNodeWeights = CELLNODEPARENTNODEWEIGHTS1 + case(2) + self%cellNodeParentNodeWeights = CELLNODEPARENTNODEWEIGHTS2 + case(3) + self%cellNodeParentNodeWeights = CELLNODEPARENTNODEWEIGHTS3 + case(4) + self%cellNodeParentNodeWeights = CELLNODEPARENTNODEWEIGHTS4 + case(5) + self%cellNodeParentNodeWeights = CELLNODEPARENTNODEWEIGHTS5 + case(6) + self%cellNodeParentNodeWeights = CELLNODEPARENTNODEWEIGHTS6 + case(7) + self%cellNodeParentNodeWeights = CELLNODEPARENTNODEWEIGHTS7 + case(8) + self%cellNodeParentNodeWeights = CELLNODEPARENTNODEWEIGHTS8 + case(9) + self%cellNodeParentNodeWeights = CELLNODEPARENTNODEWEIGHTS9 + case(10) + self%cellNodeParentNodeWeights = CELLNODEPARENTNODEWEIGHTS10 + case(11) + self%cellNodeParentNodeWeights = CELLNODEPARENTNODEWEIGHTS11 + case(12) + self%cellNodeParentNodeWeights = CELLNODEPARENTNODEWEIGHTS12 + case(13) + self%cellNodeParentNodeWeights = CELLNODEPARENTNODEWEIGHTS13 + case default + call IO_error(0,ext_msg='invalid element type') + end select + + + self%NcellNodes = NCELLNODE(self%geomType) + self%nIPs = NIP (self%geomType) + self%cellType = CELLTYPE (self%geomType) + + select case (self%geomType) + case(1) + self%IPneighbor = IPNEIGHBOR1 + self%cell = CELL1 + case(2) + self%IPneighbor = IPNEIGHBOR2 + self%cell = CELL2 + case(3) + self%IPneighbor = IPNEIGHBOR3 + self%cell = CELL3 + case(4) + self%IPneighbor = IPNEIGHBOR4 + self%cell = CELL4 + case(5) + self%IPneighbor = IPNEIGHBOR5 + self%cell = CELL5 + case(6) + self%IPneighbor = IPNEIGHBOR6 + self%cell = CELL6 + case(7) + self%IPneighbor = IPNEIGHBOR7 + self%cell = CELL7 + case(8) + self%IPneighbor = IPNEIGHBOR8 + self%cell = CELL8 + case(9) + self%IPneighbor = IPNEIGHBOR9 + self%cell = CELL9 + case(10) + self%IPneighbor = IPNEIGHBOR10 + self%cell = CELL10 + end select + + self%NcellnodesPerCell = NCELLNODEPERCELL(self%cellType) + + select case(self%cellType) + case(1) + self%cellFace = CELLFACE1 + self%vtkType = 'TRIANGLE' + case(2) + self%cellFace = CELLFACE2 + self%vtkType = 'QUAD' + case(3) + self%cellFace = CELLFACE3 + self%vtkType = 'TETRA' + case(4) + self%cellFace = CELLFACE4 + self%vtkType = 'HEXAHEDRON' + end select + + self%nIPneighbors = size(self%IPneighbor,1) + + write(6,'(/,a)') ' <<<+- element_init -+>>>'; flush(6) + + write(6,*) ' element type: ',self%elemType + write(6,*) ' geom type: ',self%geomType + write(6,*) ' cell type: ',self%cellType + write(6,*) ' # node: ',self%Nnodes + write(6,*) ' # IP: ',self%nIPs + write(6,*) ' # cellnode: ',self%Ncellnodes + write(6,*) ' # cellnode/cell: ',self%NcellnodesPerCell + write(6,*) ' # IP neighbor: ',self%nIPneighbors + end subroutine tElement_init end module element diff --git a/src/mesh_grid.f90 b/src/mesh_grid.f90 index a24fee9bc..2fd4ee85c 100644 --- a/src/mesh_grid.f90 +++ b/src/mesh_grid.f90 @@ -61,7 +61,7 @@ subroutine mesh_init(ip,el) integer(C_INTPTR_T) :: & devNull, z, z_offset - write(6,'(/,a)') ' <<<+- mesh_grid init -+>>>'; flush(6) + write(6,'(/,a)') ' <<<+- mesh_grid init -+>>>'; flush(6) call readGeom(grid,geomSize,origin,microstructureAt,homogenizationAt) diff --git a/src/mesh_marc.f90 b/src/mesh_marc.f90 index 00fd76326..504f33f23 100644 --- a/src/mesh_marc.f90 +++ b/src/mesh_marc.f90 @@ -72,7 +72,7 @@ subroutine mesh_init(ip,el) real(pReal), dimension(:,:,:,:),allocatable :: & unscaledNormals - write(6,'(/,a)') ' <<<+- mesh init -+>>>' + write(6,'(/,a)') ' <<<+- mesh init -+>>>'; flush(6) mesh_unitlength = numerics_unitlength ! set physical extent of a length unit in mesh @@ -82,9 +82,8 @@ subroutine mesh_init(ip,el) if (debug_e < 1 .or. debug_e > nElems) call IO_error(602,ext_msg='element') if (debug_i < 1 .or. debug_i > elem%nIPs) call IO_error(602,ext_msg='IP') - FEsolving_execElem = [ 1,nElems ] ! parallel loop bounds set to comprise all DAMASK elements - allocate(FEsolving_execIP(2,nElems), source=1) ! parallel loop bounds set to comprise from first IP... - FEsolving_execIP(2,:) = elem%nIPs + FEsolving_execElem = [1,nElems] + FEsolving_execIP = spread([1,elem%nIPs],2,nElems) allocate(calcMode(elem%nIPs,nElems),source=.false.) ! pretend to have collected what first call is asking (F = I) calcMode(ip,mesh_FEasCP('elem',el)) = .true. ! first ip,el needs to be already pingponged to "calc" @@ -103,9 +102,9 @@ subroutine mesh_init(ip,el) call discretization_init(microstructureAt,homogenizationAt,& ip_reshaped,& - node0_elem) + node0_cell) - call writeGeometry(0,connectivity_elem,& + call writeGeometry(elem,connectivity_elem,& reshape(connectivity_cell,[elem%NcellNodesPerCell,elem%nIPs*nElems]),& node0_cell,ip_reshaped) @@ -121,13 +120,14 @@ end subroutine mesh_init !-------------------------------------------------------------------------------------------------- -!> @brief Writes all information needed for the DADF5 geometry +!> @brief Write all information needed for the DADF5 geometry !-------------------------------------------------------------------------------------------------- -subroutine writeGeometry(elemType, & +subroutine writeGeometry(elem, & connectivity_elem,connectivity_cell, & coordinates_nodes,coordinates_points) - integer, intent(in) :: elemType + type(tElement), intent(in) :: & + elem integer, dimension(:,:), intent(in) :: & connectivity_elem, & connectivity_cell @@ -150,28 +150,31 @@ subroutine writeGeometry(elemType, & connectivity_temp = connectivity_cell call results_writeDataset('geometry',connectivity_temp,'T_c', & 'connectivity of the cells','-') + call results_addAttribute('VTK_TYPE',elem%vtkType,'geometry/T_c') coordinates_temp = coordinates_nodes call results_writeDataset('geometry',coordinates_temp,'x_n', & - 'coordinates of the nodes','m') + 'initial coordinates of the nodes','m') coordinates_temp = coordinates_points call results_writeDataset('geometry',coordinates_temp,'x_p', & - 'coordinates of the material points','m') + 'initial coordinates of the materialpoints','m') call results_closeJobFile end subroutine writeGeometry +!-------------------------------------------------------------------------------------------------- +!> @brief Read mesh from marc input file +!-------------------------------------------------------------------------------------------------- subroutine inputRead(elem,node0_elem,connectivity_elem,microstructureAt,homogenizationAt) type(tElement), intent(out) :: elem real(pReal), dimension(:,:), allocatable, intent(out) :: & node0_elem !< node x,y,z coordinates (initially!) - integer, dimension(:,:), allocatable, intent(out) :: & + integer, dimension(:,:), allocatable, intent(out) :: & connectivity_elem - integer, dimension(:), allocatable, intent(out) :: & microstructureAt, & homogenizationAt @@ -207,16 +210,18 @@ subroutine inputRead(elem,node0_elem,connectivity_elem,microstructureAt,homogeni call IO_open_inputFile(FILEUNIT) ! ToDo: It would be better to use fileContent call inputRead_mapElemSets(nameElemSet,mapElemSet,& - FILEUNIT) + FILEUNIT,inputFile) + + call inputRead_elemType(elem, & + nElems,inputFile) allocate (mesh_mapFEtoCPelem(2,nElems), source=0) - call inputRead_mapElems(hypoelasticTableStyle,nameElemSet,mapElemSet,fileFormatVersion,matNumber,FILEUNIT) + call inputRead_mapElems(elem%nNodes,nElems,& + inputFile) allocate (mesh_mapFEtoCPnode(2,Nnodes), source=0) call inputRead_mapNodes(inputFile) - call inputRead_elemType(elem, & - nElems,inputFile) call inputRead_elemNodes(node0_elem, & Nnodes,inputFile) @@ -347,47 +352,63 @@ end subroutine inputRead_NnodesAndElements !> @brief Count overall number of element sets in mesh. !-------------------------------------------------------------------------------------------------- subroutine inputRead_NelemSets(nElemSets,maxNelemInSet,& - fileUnit) + fileContent) - integer, intent(out) :: nElemSets, maxNelemInSet - integer, intent(in) :: fileUnit + integer, intent(out) :: nElemSets, maxNelemInSet + character(len=pStringLen), dimension(:), intent(in) :: fileContent !< file content, separated per lines integer, allocatable, dimension(:) :: chunkPos - character(len=300) :: line + integer :: i,l,elemInCurrentSet nElemSets = 0 maxNelemInSet = 0 - rewind(fileUnit) - do - read (fileUnit,'(A300)',END=620) line - chunkPos = IO_stringPos(line) - - if ( IO_lc(IO_StringValue(line,chunkPos,1)) == 'define' .and. & - IO_lc(IO_StringValue(line,chunkPos,2)) == 'element' ) then + do l = 1, size(fileContent) + chunkPos = IO_stringPos(fileContent(l)) + if ( IO_lc(IO_StringValue(fileContent(l),chunkPos,1)) == 'define' .and. & + IO_lc(IO_StringValue(fileContent(l),chunkPos,2)) == 'element' ) then nElemSets = nElemSets + 1 - maxNelemInSet = max(maxNelemInSet, IO_countContinuousIntValues(fileUnit)) + + chunkPos = IO_stringPos(fileContent(l+1)) + if(IO_lc(IO_StringValue(fileContent(l+1),chunkPos,2)) == 'to' ) then + elemInCurrentSet = 1 + abs( IO_intValue(fileContent(l+1),chunkPos,3) & + -IO_intValue(fileContent(l+1),chunkPos,1)) + else + elemInCurrentSet = 0 + i = 0 + do while (.true.) + i = i + 1 + chunkPos = IO_stringPos(fileContent(l+i)) + elemInCurrentSet = elemInCurrentSet + chunkPos(1) - 1 ! add line's count when assuming 'c' + if(IO_lc(IO_stringValue(fileContent(l+i),chunkPos,chunkPos(1))) /= 'c') then ! line finished, read last value + elemInCurrentSet = elemInCurrentSet + 1 ! data ended + exit + endif + enddo + endif + maxNelemInSet = max(maxNelemInSet, elemInCurrentSet) endif enddo -620 end subroutine inputRead_NelemSets +end subroutine inputRead_NelemSets !-------------------------------------------------------------------------------------------------- !> @brief map element sets !-------------------------------------------------------------------------------------------------- -subroutine inputRead_mapElemSets(nameElemSet,mapElemSet,fileUnit) +subroutine inputRead_mapElemSets(nameElemSet,mapElemSet,fileUnit,fileContent) character(len=64), dimension(:), allocatable, intent(out) :: nameElemSet integer, dimension(:,:), allocatable, intent(out) :: mapElemSet integer, intent(in) :: fileUnit + character(len=pStringLen), dimension(:), intent(in) :: fileContent !< file content, separated per lines integer, allocatable, dimension(:) :: chunkPos character(len=300) :: line integer :: elemSet, NelemSets, maxNelemInSet - call inputRead_NelemSets(NelemSets,maxNelemInSet,fileUnit) + call inputRead_NelemSets(NelemSets,maxNelemInSet,fileContent) allocate(nameElemSet(NelemSets)); nameElemSet = 'n/a' allocate(mapElemSet(1+maxNelemInSet,NelemSets),source=0) elemSet = 0 @@ -411,65 +432,35 @@ subroutine inputRead_mapElemSets(nameElemSet,mapElemSet,fileUnit) !-------------------------------------------------------------------------------------------------- !> @brief Maps elements from FE ID to internal (consecutive) representation. !-------------------------------------------------------------------------------------------------- -subroutine inputRead_mapElems(tableStyle,nameElemSet,mapElemSet,fileFormatVersion,matNumber,fileUnit) - - integer, intent(in) :: fileUnit,tableStyle,fileFormatVersion - integer, dimension(:), intent(in) :: matNumber - character(len=64), intent(in), dimension(:) :: nameElemSet - integer, dimension(:,:), intent(in) :: & - mapElemSet +subroutine inputRead_mapElems(nNodes,nElem,fileContent) + integer, intent(in) :: & + nElem, & + nNodes !< number of nodes per element + character(len=pStringLen), dimension(:), intent(in) :: fileContent !< file content, separated per lines + integer, allocatable, dimension(:) :: chunkPos - character(len=300) :: line, & - tmp - - integer, dimension(:), allocatable :: contInts - integer :: i,cpElem - - allocate(contInts(size(mesh_mapFEtoCPelem,2)+1)) - - cpElem = 0 - contInts = 0 - rewind(fileUnit) - do - read (fileUnit,'(A300)',END=620) line - chunkPos = IO_stringPos(line) - Marc2016andEarlier: if (fileFormatVersion < 13) then - if( IO_lc(IO_stringValue(line,chunkPos,1)) == 'hypoelastic' ) then - skipLines: do i=1,3+TableStyle - read (fileUnit,'(A300)') line - enddo skipLines - contInts = IO_continuousIntValues(fileUnit,size(mesh_mapFEtoCPelem,2),nameElemSet,& - mapElemSet,size(nameElemSet)) - exit - endif - else Marc2016andEarlier - if ( IO_lc(IO_stringValue(line,chunkPos,1)) == 'connectivity') then - read (fileUnit,'(A300)',END=620) line - chunkPos = IO_stringPos(line) - if(any(matNumber==IO_intValue(line,chunkPos,6))) then - do - read (fileUnit,'(A300)',END=620) line - chunkPos = IO_stringPos(line) - tmp = IO_lc(IO_stringValue(line,chunkPos,1)) - if (verify(trim(tmp),"0123456789")/=0) then ! found keyword - exit - else - contInts(1) = contInts(1) + 1 - read (tmp,*) contInts(contInts(1)+1) - endif - enddo - endif - endif - endif Marc2016andEarlier - enddo -620 do i = 1,contInts(1) - cpElem = cpElem+1 - mesh_mapFEtoCPelem(1,cpElem) = contInts(1+i) - mesh_mapFEtoCPelem(2,cpElem) = cpElem + integer :: i,j,l,nNodesAlreadyRead + + do l = 1, size(fileContent) + chunkPos = IO_stringPos(fileContent(l)) + if( IO_lc(IO_stringValue(fileContent(l),chunkPos,1)) == 'connectivity' ) then + j = 0 + do i = 1,nElem + chunkPos = IO_stringPos(fileContent(l+1+i+j)) + mesh_mapFEtoCPelem(:,i) = [IO_intValue(fileContent(l+1+i+j),chunkPos,1),i] + nNodesAlreadyRead = chunkPos(1) - 2 + do while(nNodesAlreadyRead < nNodes) ! read on if not all nodes in one line + j = j + 1 + chunkPos = IO_stringPos(fileContent(l+1+i+j)) + nNodesAlreadyRead = nNodesAlreadyRead + chunkPos(1) + enddo + enddo + exit + endif enddo - -call math_sort(mesh_mapFEtoCPelem) + + call math_sort(mesh_mapFEtoCPelem) end subroutine inputRead_mapElems @@ -488,7 +479,8 @@ subroutine inputRead_mapNodes(fileContent) chunkPos = IO_stringPos(fileContent(l)) if( IO_lc(IO_stringValue(fileContent(l),chunkPos,1)) == 'coordinates' ) then do i = 1,size(mesh_mapFEtoCPnode,2) - mesh_mapFEtoCPnode(1:2,i) = [IO_fixedIntValue (fileContent(l+1+i),[0,10],1),i] ! ToDo: use IO_intValue + chunkPos = IO_stringPos(fileContent(l+1+i)) + mesh_mapFEtoCPnode(:,i) = [IO_intValue(fileContent(l+1+i),chunkPos,1),i] enddo exit endif @@ -509,7 +501,6 @@ subroutine inputRead_elemNodes(nodes, & integer, intent(in) :: nNode character(len=pStringLen), dimension(:), intent(in) :: fileContent !< file content, separated per lines - integer, dimension(5), parameter :: node_ends = [0,10,30,50,70] integer, allocatable, dimension(:) :: chunkPos integer :: i,j,m,l @@ -519,9 +510,10 @@ subroutine inputRead_elemNodes(nodes, & chunkPos = IO_stringPos(fileContent(l)) if( IO_lc(IO_stringValue(fileContent(l),chunkPos,1)) == 'coordinates' ) then do i=1,nNode - m = mesh_FEasCP('node',IO_fixedIntValue(fileContent(l+1+i),node_ends,1)) !ToDo: use IO_intValue + chunkPos = IO_stringPos(fileContent(l+1+i)) + m = mesh_FEasCP('node',IO_intValue(fileContent(l+1+i),chunkPos,1)) do j = 1,3 - nodes(j,m) = mesh_unitlength * IO_fixedNoEFloatValue(fileContent(l+1+i),node_ends,j+1) !ToDo: use IO_floatValue + nodes(j,m) = mesh_unitlength * IO_floatValue(fileContent(l+1+i),chunkPos,j+1) enddo enddo exit @@ -709,7 +701,7 @@ subroutine inputRead_microstructureAndHomogenization(microstructureAt,homogeniza read (fileUnit,'(A300)',END=630) line ! read line with value of state var chunkPos = IO_stringPos(line) do while (scan(IO_stringValue(line,chunkPos,1),'+-',back=.true.)>1) ! is noEfloat value? - myVal = nint(IO_fixedNoEFloatValue(line,[0,20],1),pInt) ! state var's value + myVal = nint(IO_floatValue(line,chunkPos,1)) if (initialcondTableStyle == 2) then read (fileUnit,'(A300)',END=630) line ! read extra line read (fileUnit,'(A300)',END=630) line ! read extra line diff --git a/src/quaternions.f90 b/src/quaternions.f90 index 0ca404880..0215aca6e 100644 --- a/src/quaternions.f90 +++ b/src/quaternions.f90 @@ -48,10 +48,7 @@ module quaternions procedure, private :: pow_scal__ generic, public :: operator(**) => pow_quat__, pow_scal__ - procedure, public :: abs__ - procedure, public :: dot_product__ - procedure, public :: exp__ - procedure, public :: log__ + procedure, public :: abs => abs__ procedure, public :: conjg => conjg__ procedure, public :: real => real__ procedure, public :: aimag => aimag__ @@ -317,17 +314,17 @@ end function pow_scal__ !--------------------------------------------------------------------------------------------------- !> take exponential !--------------------------------------------------------------------------------------------------- -type(quaternion) elemental pure function exp__(self) +type(quaternion) elemental pure function exp__(a) - class(quaternion), intent(in) :: self + class(quaternion), intent(in) :: a real(pReal) :: absImag - absImag = norm2(aimag(self)) + absImag = norm2(aimag(a)) - exp__ = merge(exp(self%w) * [ cos(absImag), & - self%x/absImag * sin(absImag), & - self%y/absImag * sin(absImag), & - self%z/absImag * sin(absImag)], & + exp__ = merge(exp(a%w) * [ cos(absImag), & + a%x/absImag * sin(absImag), & + a%y/absImag * sin(absImag), & + a%z/absImag * sin(absImag)], & IEEE_value(1.0_pReal,IEEE_SIGNALING_NAN), & dNeq0(absImag)) @@ -337,17 +334,17 @@ end function exp__ !--------------------------------------------------------------------------------------------------- !> take logarithm !--------------------------------------------------------------------------------------------------- -type(quaternion) elemental pure function log__(self) +type(quaternion) elemental pure function log__(a) - class(quaternion), intent(in) :: self + class(quaternion), intent(in) :: a real(pReal) :: absImag - absImag = norm2(aimag(self)) + absImag = norm2(aimag(a)) - log__ = merge([log(abs(self)), & - self%x/absImag * acos(self%w/abs(self)), & - self%y/absImag * acos(self%w/abs(self)), & - self%z/absImag * acos(self%w/abs(self))], & + log__ = merge([log(abs(a)), & + a%x/absImag * acos(a%w/abs(a)), & + a%y/absImag * acos(a%w/abs(a)), & + a%z/absImag * acos(a%w/abs(a))], & IEEE_value(1.0_pReal,IEEE_SIGNALING_NAN), & dNeq0(absImag)) @@ -357,11 +354,11 @@ end function log__ !--------------------------------------------------------------------------------------------------- !> return norm !--------------------------------------------------------------------------------------------------- -real(pReal) elemental pure function abs__(a) +real(pReal) elemental pure function abs__(self) - class(quaternion), intent(in) :: a + class(quaternion), intent(in) :: self - abs__ = norm2([a%w,a%x,a%y,a%z]) + abs__ = norm2([self%w,self%x,self%y,self%z]) end function abs__ @@ -381,11 +378,11 @@ end function dot_product__ !--------------------------------------------------------------------------------------------------- !> take conjugate complex !--------------------------------------------------------------------------------------------------- -type(quaternion) elemental pure function conjg__(a) +type(quaternion) elemental pure function conjg__(self) - class(quaternion), intent(in) :: a + class(quaternion), intent(in) :: self - conjg__ = [a%w, -a%x, -a%y, -a%z] + conjg__ = [self%w,-self%x,-self%y,-self%z] end function conjg__ @@ -464,7 +461,7 @@ subroutine unitTest call random_number(qu) qu = (qu-0.5_pReal) * 2.0_pReal q = quaternion(qu) - + q_2= qu if(any(dNeq(q%asArray(),q_2%asArray()))) call IO_error(401,ext_msg='assign_vec__') @@ -479,10 +476,10 @@ subroutine unitTest q_2 = q / 0.5_pReal if(any(dNeq(q_2%asArray(),2.0_pReal*qu))) call IO_error(401,ext_msg='div__') - + q_2 = q * 0.3_pReal if(dNeq0(abs(q)) .and. q_2 == q) call IO_error(401,ext_msg='eq__') - + q_2 = q if(q_2 /= q) call IO_error(401,ext_msg='neq__') @@ -504,12 +501,12 @@ subroutine unitTest if(q /= conjg(q_2)) call IO_error(401,ext_msg='conjg/involution') if(dNeq(q_2%real(), q%real())) call IO_error(401,ext_msg='conjg/real') if(any(dNeq(q_2%aimag(),q%aimag()*(-1.0_pReal)))) call IO_error(401,ext_msg='conjg/aimag') - + if(abs(q) > 0.0_pReal) then q_2 = q * q%inverse() if( dNeq(real(q_2), 1.0_pReal,1.0e-15_pReal)) call IO_error(401,ext_msg='inverse/real') if(any(dNeq0(aimag(q_2), 1.0e-15_pReal))) call IO_error(401,ext_msg='inverse/aimag') - + q_2 = q/abs(q) q_2 = conjg(q_2) - inverse(q_2) if(any(dNeq0(q_2%asArray(),1.0e-15_pReal))) call IO_error(401,ext_msg='inverse/conjg') diff --git a/src/results.f90 b/src/results.f90 index ca2a868fb..e0b5f9898 100644 --- a/src/results.f90 +++ b/src/results.f90 @@ -70,7 +70,7 @@ subroutine results_init resultsFile = HDF5_openFile(trim(getSolverJobName())//'.hdf5','w',.true.) call results_addAttribute('DADF5_version_major',0) - call results_addAttribute('DADF5_version_minor',5) + call results_addAttribute('DADF5_version_minor',6) call results_addAttribute('DAMASK_version',DAMASKVERSION) call get_command(commandLine) call results_addAttribute('call',trim(commandLine)) diff --git a/src/rotations.f90 b/src/rotations.f90 index 5deb02a20..caec1cc74 100644 --- a/src/rotations.f90 +++ b/src/rotations.f90 @@ -42,7 +42,8 @@ ! Convention 4: Euler angle triplets are implemented using the Bunge convention, ! with the angular ranges as [0, 2π],[0, π],[0, 2π] ! Convention 5: the rotation angle ω is limited to the interval [0, π] -! Convention 6: P = -1 +! Convention 6: the real part of a quaternion is positive, Re(q) > 0 +! Convention 7: P = -1 !--------------------------------------------------------------------------------------------------- module rotations @@ -77,6 +78,7 @@ module rotations procedure, public :: rotTensor4 procedure, public :: rotTensor4sym procedure, public :: misorientation + procedure, public :: standardize end type rotation public :: & @@ -92,7 +94,7 @@ contains subroutine rotations_init call quaternions_init - write(6,'(/,a)') ' <<<+- rotations init -+>>>' + write(6,'(/,a)') ' <<<+- rotations init -+>>>'; flush(6) call unitTest end subroutine rotations_init @@ -237,7 +239,6 @@ end subroutine fromMatrix !--------------------------------------------------------------------------------------------------- !> @brief: Rotate a rotation -!> ToDo: completly untested !--------------------------------------------------------------------------------------------------- pure elemental function rotRot__(self,R) result(rRot) @@ -245,10 +246,23 @@ pure elemental function rotRot__(self,R) result(rRot) class(rotation), intent(in) :: self,R rRot = rotation(self%q*R%q) + call rRot%standardize() end function rotRot__ +!--------------------------------------------------------------------------------------------------- +!> @brief quaternion representation with positive q +!--------------------------------------------------------------------------------------------------- +pure elemental subroutine standardize(self) + + class(rotation), intent(inout) :: self + + if (real(self%q) < 0.0_pReal) self%q = self%q%homomorphed() + +end subroutine standardize + + !--------------------------------------------------------------------------------------------------- !> @author Marc De Graef, Carnegie Mellon University !> @brief rotate a vector passively (default) or actively @@ -375,7 +389,7 @@ pure elemental function misorientation(self,other) type(rotation) :: misorientation class(rotation), intent(in) :: self, other - misorientation%q = conjg(self%q) * other%q !ToDo: this is the convention used in math + misorientation%q = other%q * conjg(self%q) end function misorientation