Merge branch 'development' into python-module
This commit is contained in:
commit
54d1c35c45
|
@ -433,6 +433,15 @@ J2_plasticBehavior:
|
||||||
- master
|
- master
|
||||||
- release
|
- release
|
||||||
|
|
||||||
|
Marc_elementLib:
|
||||||
|
stage: marc
|
||||||
|
script:
|
||||||
|
- module load $IntelMarc $HDF5Marc $MSC
|
||||||
|
- Marc_elementLib/test.py
|
||||||
|
except:
|
||||||
|
- master
|
||||||
|
- release
|
||||||
|
|
||||||
###################################################################################################
|
###################################################################################################
|
||||||
grid_all_example:
|
grid_all_example:
|
||||||
stage: example
|
stage: example
|
||||||
|
|
2
PRIVATE
2
PRIVATE
|
@ -1 +1 @@
|
||||||
Subproject commit 65905b5cd467546873be8be2dd7d6217db0012be
|
Subproject commit 83cb0b6ceea2e26c554ad4c6fbe6aabaee6fa27f
|
|
@ -40,8 +40,9 @@ class DADF5():
|
||||||
self.version_major = f.attrs['DADF5-major']
|
self.version_major = f.attrs['DADF5-major']
|
||||||
self.version_minor = f.attrs['DADF5-minor']
|
self.version_minor = f.attrs['DADF5-minor']
|
||||||
|
|
||||||
if self.version_major != 0 or not 2 <= self.version_minor <= 5:
|
if self.version_major != 0 or not 2 <= self.version_minor <= 6:
|
||||||
raise TypeError('Unsupported DADF5 version {} '.format(f.attrs['DADF5-version']))
|
raise TypeError('Unsupported DADF5 version {}.{} '.format(f.attrs['DADF5_version_major'],
|
||||||
|
f.attrs['DADF5_version_minor']))
|
||||||
|
|
||||||
self.structured = 'grid' in f['geometry'].attrs.keys()
|
self.structured = 'grid' in f['geometry'].attrs.keys()
|
||||||
|
|
||||||
|
@ -886,8 +887,27 @@ class DADF5():
|
||||||
vtk_geom = vtk.vtkUnstructuredGrid()
|
vtk_geom = vtk.vtkUnstructuredGrid()
|
||||||
vtk_geom.SetPoints(nodes)
|
vtk_geom.SetPoints(nodes)
|
||||||
vtk_geom.Allocate(f['/geometry/T_c'].shape[0])
|
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']:
|
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':
|
elif mode == 'Point':
|
||||||
Points = vtk.vtkPoints()
|
Points = vtk.vtkPoints()
|
||||||
Vertices = vtk.vtkCellArray()
|
Vertices = vtk.vtkCellArray()
|
||||||
|
|
|
@ -3,9 +3,8 @@ import math
|
||||||
import numpy as np
|
import numpy as np
|
||||||
|
|
||||||
from . import Lambert
|
from . import Lambert
|
||||||
from .quaternion import Quaternion
|
|
||||||
from .quaternion import P
|
|
||||||
|
|
||||||
|
P = -1
|
||||||
|
|
||||||
####################################################################################################
|
####################################################################################################
|
||||||
class Rotation:
|
class Rotation:
|
||||||
|
@ -26,7 +25,8 @@ class Rotation:
|
||||||
Convention 4: Euler angle triplets are implemented using the Bunge convention,
|
Convention 4: Euler angle triplets are implemented using the Bunge convention,
|
||||||
with the angular ranges as [0, 2π],[0, π],[0, 2π].
|
with the angular ranges as [0, 2π],[0, π],[0, 2π].
|
||||||
Convention 5: The rotation angle ω is limited to the interval [0, π].
|
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
|
Usage
|
||||||
-----
|
-----
|
||||||
|
@ -49,10 +49,7 @@ class Rotation:
|
||||||
Unit quaternion that follows the conventions. Use .fromQuaternion to perform a sanity check.
|
Unit quaternion that follows the conventions. Use .fromQuaternion to perform a sanity check.
|
||||||
|
|
||||||
"""
|
"""
|
||||||
if isinstance(quaternion,Quaternion):
|
|
||||||
self.quaternion = quaternion.copy()
|
self.quaternion = quaternion.copy()
|
||||||
else:
|
|
||||||
self.quaternion = Quaternion(q=quaternion[0],p=quaternion[1:4])
|
|
||||||
|
|
||||||
def __copy__(self):
|
def __copy__(self):
|
||||||
"""Copy."""
|
"""Copy."""
|
||||||
|
@ -64,9 +61,9 @@ class Rotation:
|
||||||
def __repr__(self):
|
def __repr__(self):
|
||||||
"""Orientation displayed as unit quaternion, rotation matrix, and Bunge-Euler angles."""
|
"""Orientation displayed as unit quaternion, rotation matrix, and Bunge-Euler angles."""
|
||||||
return '\n'.join([
|
return '\n'.join([
|
||||||
'{}'.format(self.quaternion),
|
'Quaternion: (real={:.3f}, imag=<{:+.3f}, {:+.3f}, {:+.3f}>)'.format(*(self.quaternion)),
|
||||||
'Matrix:\n{}'.format( '\n'.join(['\t'.join(list(map(str,self.asMatrix()[i,:]))) for i in range(3)]) ),
|
'Matrix:\n{}'.format(self.asMatrix()),
|
||||||
'Bunge Eulers / deg: {}'.format('\t'.join(list(map(str,self.asEulers(degrees=True)))) ),
|
'Bunge Eulers / deg: ({:3.2f}, {:3.2f}, {:3.2f})'.format(*self.asEulers(degrees=True)),
|
||||||
])
|
])
|
||||||
|
|
||||||
def __mul__(self, other):
|
def __mul__(self, other):
|
||||||
|
@ -85,19 +82,25 @@ class Rotation:
|
||||||
|
|
||||||
"""
|
"""
|
||||||
if isinstance(other, Rotation): # rotate a rotation
|
if isinstance(other, Rotation): # rotate a rotation
|
||||||
return self.__class__(self.quaternion * other.quaternion).standardize()
|
self_q = self.quaternion[0]
|
||||||
elif isinstance(other, np.ndarray):
|
self_p = self.quaternion[1:]
|
||||||
if other.shape == (3,): # rotate a single (3)-vector
|
other_q = other.quaternion[0]
|
||||||
( x, y, z) = self.quaternion.p
|
other_p = other.quaternion[1:]
|
||||||
(Vx,Vy,Vz) = other[0:3]
|
R = self.__class__(np.append(self_q*other_q - np.dot(self_p,other_p),
|
||||||
A = self.quaternion.q*self.quaternion.q - np.dot(self.quaternion.p,self.quaternion.p)
|
self_q*other_p + other_q*self_p + P * np.cross(self_p,other_p)))
|
||||||
B = 2.0 * (x*Vx + y*Vy + z*Vz)
|
return R.standardize()
|
||||||
C = 2.0 * P*self.quaternion.q
|
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([
|
return np.array([
|
||||||
A*Vx + B*x + C*(y*Vz - z*Vy),
|
A*other[0] + B*self.quaternion[1] + C*(self.quaternion[2]*other[2] - self.quaternion[3]*other[1]),
|
||||||
A*Vy + B*y + C*(z*Vx - x*Vz),
|
A*other[1] + B*self.quaternion[2] + C*(self.quaternion[3]*other[0] - self.quaternion[1]*other[2]),
|
||||||
A*Vz + B*z + C*(x*Vy - y*Vx),
|
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
|
elif other.shape == (3,3,): # rotate a single (3x3)-matrix
|
||||||
return np.dot(self.asMatrix(),np.dot(other,self.asMatrix().T))
|
return np.dot(self.asMatrix(),np.dot(other,self.asMatrix().T))
|
||||||
|
@ -105,25 +108,13 @@ class Rotation:
|
||||||
raise NotImplementedError
|
raise NotImplementedError
|
||||||
else:
|
else:
|
||||||
return NotImplemented
|
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:
|
else:
|
||||||
return NotImplemented
|
return NotImplemented
|
||||||
|
|
||||||
|
|
||||||
def inverse(self):
|
def inverse(self):
|
||||||
"""In-place inverse rotation/backward rotation."""
|
"""In-place inverse rotation/backward rotation."""
|
||||||
self.quaternion.conjugate()
|
self.quaternion[1:] *= -1
|
||||||
return self
|
return self
|
||||||
|
|
||||||
def inversed(self):
|
def inversed(self):
|
||||||
|
@ -133,7 +124,7 @@ class Rotation:
|
||||||
|
|
||||||
def standardize(self):
|
def standardize(self):
|
||||||
"""In-place quaternion representation with positive q."""
|
"""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
|
return self
|
||||||
|
|
||||||
def standardized(self):
|
def standardized(self):
|
||||||
|
@ -170,8 +161,7 @@ class Rotation:
|
||||||
################################################################################################
|
################################################################################################
|
||||||
# convert to different orientation representations (numpy arrays)
|
# convert to different orientation representations (numpy arrays)
|
||||||
|
|
||||||
def asQuaternion(self,
|
def asQuaternion(self):
|
||||||
quaternion = False):
|
|
||||||
"""
|
"""
|
||||||
Unit quaternion [q, p_1, p_2, p_3] unless quaternion == True: damask.quaternion object.
|
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 quaternion as DAMASK object.
|
||||||
|
|
||||||
"""
|
"""
|
||||||
return self.quaternion if quaternion else self.quaternion.asArray()
|
return self.quaternion
|
||||||
|
|
||||||
def asEulers(self,
|
def asEulers(self,
|
||||||
degrees = False):
|
degrees = False):
|
||||||
|
@ -194,7 +184,7 @@ class Rotation:
|
||||||
return angles in degrees.
|
return angles in degrees.
|
||||||
|
|
||||||
"""
|
"""
|
||||||
eu = qu2eu(self.quaternion.asArray())
|
eu = qu2eu(self.quaternion)
|
||||||
if degrees: eu = np.degrees(eu)
|
if degrees: eu = np.degrees(eu)
|
||||||
return eu
|
return eu
|
||||||
|
|
||||||
|
@ -212,13 +202,13 @@ class Rotation:
|
||||||
return tuple of axis and angle.
|
return tuple of axis and angle.
|
||||||
|
|
||||||
"""
|
"""
|
||||||
ax = qu2ax(self.quaternion.asArray())
|
ax = qu2ax(self.quaternion)
|
||||||
if degrees: ax[3] = np.degrees(ax[3])
|
if degrees: ax[3] = np.degrees(ax[3])
|
||||||
return (ax[:3],np.degrees(ax[3])) if pair else ax
|
return (ax[:3],np.degrees(ax[3])) if pair else ax
|
||||||
|
|
||||||
def asMatrix(self):
|
def asMatrix(self):
|
||||||
"""Rotation matrix."""
|
"""Rotation matrix."""
|
||||||
return qu2om(self.quaternion.asArray())
|
return qu2om(self.quaternion)
|
||||||
|
|
||||||
def asRodrigues(self,
|
def asRodrigues(self,
|
||||||
vector = False):
|
vector = False):
|
||||||
|
@ -231,16 +221,16 @@ class Rotation:
|
||||||
return as actual Rodrigues--Frank vector, i.e. rotation axis scaled by tan(ω/2).
|
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
|
return ro[:3]*ro[3] if vector else ro
|
||||||
|
|
||||||
def asHomochoric(self):
|
def asHomochoric(self):
|
||||||
"""Homochoric vector: (h_1, h_2, h_3)."""
|
"""Homochoric vector: (h_1, h_2, h_3)."""
|
||||||
return qu2ho(self.quaternion.asArray())
|
return qu2ho(self.quaternion)
|
||||||
|
|
||||||
def asCubochoric(self):
|
def asCubochoric(self):
|
||||||
"""Cubochoric vector: (c_1, c_2, c_3)."""
|
"""Cubochoric vector: (c_1, c_2, c_3)."""
|
||||||
return qu2cu(self.quaternion.asArray())
|
return qu2cu(self.quaternion)
|
||||||
|
|
||||||
def asM(self):
|
def asM(self):
|
||||||
"""
|
"""
|
||||||
|
@ -252,15 +242,14 @@ class Rotation:
|
||||||
https://doi.org/10.2514/1.28949
|
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
|
# static constructors. The input data needs to follow the convention, options allow to
|
||||||
# relax these convections
|
# relax these convections
|
||||||
@classmethod
|
@staticmethod
|
||||||
def fromQuaternion(cls,
|
def fromQuaternion(quaternion,
|
||||||
quaternion,
|
|
||||||
acceptHomomorph = False,
|
acceptHomomorph = False,
|
||||||
P = -1):
|
P = -1):
|
||||||
|
|
||||||
|
@ -275,11 +264,10 @@ class Rotation:
|
||||||
if not np.isclose(np.linalg.norm(qu), 1.0):
|
if not np.isclose(np.linalg.norm(qu), 1.0):
|
||||||
raise ValueError('Quaternion is not of unit length.\n{} {} {} {}'.format(*qu))
|
raise ValueError('Quaternion is not of unit length.\n{} {} {} {}'.format(*qu))
|
||||||
|
|
||||||
return cls(qu)
|
return Rotation(qu)
|
||||||
|
|
||||||
@classmethod
|
@staticmethod
|
||||||
def fromEulers(cls,
|
def fromEulers(eulers,
|
||||||
eulers,
|
|
||||||
degrees = False):
|
degrees = False):
|
||||||
|
|
||||||
eu = eulers if isinstance(eulers, np.ndarray) and eulers.dtype == np.dtype(float) \
|
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:
|
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))
|
raise ValueError('Euler angles outside of [0..2π],[0..π],[0..2π].\n{} {} {}.'.format(*eu))
|
||||||
|
|
||||||
return cls(eu2qu(eu))
|
return Rotation(eu2qu(eu))
|
||||||
|
|
||||||
@classmethod
|
@staticmethod
|
||||||
def fromAxisAngle(cls,
|
def fromAxisAngle(angleAxis,
|
||||||
angleAxis,
|
|
||||||
degrees = False,
|
degrees = False,
|
||||||
normalise = False,
|
normalise = False,
|
||||||
P = -1):
|
P = -1):
|
||||||
|
@ -307,11 +294,10 @@ class Rotation:
|
||||||
if not np.isclose(np.linalg.norm(ax[0:3]), 1.0):
|
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]))
|
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
|
@staticmethod
|
||||||
def fromBasis(cls,
|
def fromBasis(basis,
|
||||||
basis,
|
|
||||||
orthonormal = True,
|
orthonormal = True,
|
||||||
reciprocal = False,
|
reciprocal = False,
|
||||||
):
|
):
|
||||||
|
@ -330,18 +316,16 @@ class Rotation:
|
||||||
or not np.isclose(np.dot(om[2],om[0]), 0.0):
|
or not np.isclose(np.dot(om[2],om[0]), 0.0):
|
||||||
raise ValueError('matrix is not orthogonal.\n{}'.format(om))
|
raise ValueError('matrix is not orthogonal.\n{}'.format(om))
|
||||||
|
|
||||||
return cls(om2qu(om))
|
return Rotation(om2qu(om))
|
||||||
|
|
||||||
@classmethod
|
@staticmethod
|
||||||
def fromMatrix(cls,
|
def fromMatrix(om,
|
||||||
om,
|
|
||||||
):
|
):
|
||||||
|
|
||||||
return cls.fromBasis(om)
|
return Rotation.fromBasis(om)
|
||||||
|
|
||||||
@classmethod
|
@staticmethod
|
||||||
def fromRodrigues(cls,
|
def fromRodrigues(rodrigues,
|
||||||
rodrigues,
|
|
||||||
normalise = False,
|
normalise = False,
|
||||||
P = -1):
|
P = -1):
|
||||||
|
|
||||||
|
@ -354,22 +338,20 @@ class Rotation:
|
||||||
if ro[3] < 0.0:
|
if ro[3] < 0.0:
|
||||||
raise ValueError('Rodriques rotation angle not positive.\n'.format(ro[3]))
|
raise ValueError('Rodriques rotation angle not positive.\n'.format(ro[3]))
|
||||||
|
|
||||||
return cls(ro2qu(ro))
|
return Rotation(ro2qu(ro))
|
||||||
|
|
||||||
@classmethod
|
@staticmethod
|
||||||
def fromHomochoric(cls,
|
def fromHomochoric(homochoric,
|
||||||
homochoric,
|
|
||||||
P = -1):
|
P = -1):
|
||||||
|
|
||||||
ho = homochoric if isinstance(homochoric, np.ndarray) and homochoric.dtype == np.dtype(float) \
|
ho = homochoric if isinstance(homochoric, np.ndarray) and homochoric.dtype == np.dtype(float) \
|
||||||
else np.array(homochoric,dtype=float)
|
else np.array(homochoric,dtype=float)
|
||||||
if P > 0: ho *= -1 # convert from P=1 to P=-1
|
if P > 0: ho *= -1 # convert from P=1 to P=-1
|
||||||
|
|
||||||
return cls(ho2qu(ho))
|
return Rotation(ho2qu(ho))
|
||||||
|
|
||||||
@classmethod
|
@staticmethod
|
||||||
def fromCubochoric(cls,
|
def fromCubochoric(cubochoric,
|
||||||
cubochoric,
|
|
||||||
P = -1):
|
P = -1):
|
||||||
|
|
||||||
cu = cubochoric if isinstance(cubochoric, np.ndarray) and cubochoric.dtype == np.dtype(float) \
|
cu = cubochoric if isinstance(cubochoric, np.ndarray) and cubochoric.dtype == np.dtype(float) \
|
||||||
|
@ -377,12 +359,11 @@ class Rotation:
|
||||||
ho = cu2ho(cu)
|
ho = cu2ho(cu)
|
||||||
if P > 0: ho *= -1 # convert from P=1 to P=-1
|
if P > 0: ho *= -1 # convert from P=1 to P=-1
|
||||||
|
|
||||||
return cls(ho2qu(ho))
|
return Rotation(ho2qu(ho))
|
||||||
|
|
||||||
|
|
||||||
@classmethod
|
@staticmethod
|
||||||
def fromAverage(cls,
|
def fromAverage(rotations,
|
||||||
rotations,
|
|
||||||
weights = []):
|
weights = []):
|
||||||
"""
|
"""
|
||||||
Average rotation.
|
Average rotation.
|
||||||
|
@ -412,19 +393,18 @@ class Rotation:
|
||||||
else M + r.asM() * n # noqa add (multiples) of this rotation to average noqa
|
else M + r.asM() * n # noqa add (multiples) of this rotation to average noqa
|
||||||
eig, vec = np.linalg.eig(M/N)
|
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
|
@staticmethod
|
||||||
def fromRandom(cls):
|
def fromRandom():
|
||||||
r = np.random.random(3)
|
r = np.random.random(3)
|
||||||
A = np.sqrt(r[2])
|
A = np.sqrt(r[2])
|
||||||
B = np.sqrt(1.0-r[2])
|
B = np.sqrt(1.0-r[2])
|
||||||
w = np.cos(2.0*np.pi*r[0])*A
|
return Rotation(np.array([np.cos(2.0*np.pi*r[0])*A,
|
||||||
x = np.sin(2.0*np.pi*r[1])*B
|
np.sin(2.0*np.pi*r[1])*B,
|
||||||
y = np.cos(2.0*np.pi*r[1])*B
|
np.cos(2.0*np.pi*r[1])*B,
|
||||||
z = np.sin(2.0*np.pi*r[0])*A
|
np.sin(2.0*np.pi*r[0])*A])).standardize()
|
||||||
return cls.fromQuaternion([w,x,y,z],acceptHomomorph=True)
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
@ -1193,9 +1173,8 @@ class Orientation:
|
||||||
return color
|
return color
|
||||||
|
|
||||||
|
|
||||||
@classmethod
|
@staticmethod
|
||||||
def fromAverage(cls,
|
def fromAverage(orientations,
|
||||||
orientations,
|
|
||||||
weights = []):
|
weights = []):
|
||||||
"""Create orientation from average of list of orientations."""
|
"""Create orientation from average of list of orientations."""
|
||||||
if not all(isinstance(item, Orientation) for item in orientations):
|
if not all(isinstance(item, Orientation) for item in orientations):
|
||||||
|
|
|
@ -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()
|
|
|
@ -1,4 +1,5 @@
|
||||||
import os
|
import os
|
||||||
|
from itertools import permutations
|
||||||
|
|
||||||
import pytest
|
import pytest
|
||||||
import numpy as np
|
import numpy as np
|
||||||
|
@ -6,6 +7,7 @@ import numpy as np
|
||||||
import damask
|
import damask
|
||||||
from damask import Rotation
|
from damask import Rotation
|
||||||
from damask import Orientation
|
from damask import Orientation
|
||||||
|
from damask import Lattice
|
||||||
|
|
||||||
n = 1000
|
n = 1000
|
||||||
|
|
||||||
|
@ -45,19 +47,37 @@ class TestRotation:
|
||||||
def test_Homochoric(self,default):
|
def test_Homochoric(self,default):
|
||||||
for rot in default:
|
for rot in default:
|
||||||
assert np.allclose(rot.asRodrigues(),
|
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):
|
def test_Cubochoric(self,default):
|
||||||
for rot in default:
|
for rot in default:
|
||||||
assert np.allclose(rot.asHomochoric(),
|
assert np.allclose(rot.asHomochoric(),
|
||||||
Rotation.fromCubochoric(rot.asCubochoric()).asHomochoric(),rtol=5.e-5)
|
Rotation.fromCubochoric(rot.asCubochoric()).asHomochoric())
|
||||||
|
|
||||||
def test_Quaternion(self,default):
|
def test_Quaternion(self,default):
|
||||||
for rot in default:
|
for rot in default:
|
||||||
assert np.allclose(rot.asCubochoric(),
|
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('model',['Bain','KS','GT','GT_prime','NW','Pitsch'])
|
||||||
@pytest.mark.parametrize('lattice',['fcc','bcc'])
|
@pytest.mark.parametrize('lattice',['fcc','bcc'])
|
||||||
def test_relationship_forward_backward(self,model,lattice):
|
def test_relationship_forward_backward(self,model,lattice):
|
||||||
|
|
|
@ -377,6 +377,7 @@ subroutine CPFEM_results(inc,time)
|
||||||
call constitutive_results
|
call constitutive_results
|
||||||
call crystallite_results
|
call crystallite_results
|
||||||
call homogenization_results
|
call homogenization_results
|
||||||
|
call discretization_results
|
||||||
call results_finalizeIncrement
|
call results_finalizeIncrement
|
||||||
call results_closeJobFile
|
call results_closeJobFile
|
||||||
|
|
||||||
|
|
|
@ -350,6 +350,9 @@ subroutine flux(f,ts,n,time)
|
||||||
|
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
!> @brief trigger writing of results
|
!> @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)
|
subroutine uedinc(inc,incsub)
|
||||||
use prec
|
use prec
|
||||||
|
@ -357,8 +360,12 @@ subroutine uedinc(inc,incsub)
|
||||||
|
|
||||||
implicit none
|
implicit none
|
||||||
integer, intent(in) :: inc, incsub
|
integer, intent(in) :: inc, incsub
|
||||||
|
integer, save :: inc_written
|
||||||
#include QUOTE(PASTE(./MarcInclude/creeps,Marc4DAMASK)) ! creeps is needed for timinc (time increment)
|
#include QUOTE(PASTE(./MarcInclude/creeps,Marc4DAMASK)) ! creeps is needed for timinc (time increment)
|
||||||
|
|
||||||
|
if (inc > inc_written) then
|
||||||
call CPFEM_results(inc,cptim)
|
call CPFEM_results(inc,cptim)
|
||||||
|
inc_written = inc
|
||||||
|
endif
|
||||||
|
|
||||||
end subroutine uedinc
|
end subroutine uedinc
|
||||||
|
|
118
src/IO.f90
118
src/IO.f90
|
@ -37,15 +37,12 @@ module IO
|
||||||
#if defined(Marc4DAMASK) || defined(Abaqus)
|
#if defined(Marc4DAMASK) || defined(Abaqus)
|
||||||
public :: &
|
public :: &
|
||||||
IO_open_inputFile, &
|
IO_open_inputFile, &
|
||||||
IO_countContinuousIntValues, &
|
IO_continuousIntValues
|
||||||
IO_continuousIntValues, &
|
#endif
|
||||||
#if defined(Abaqus)
|
#if defined(Abaqus)
|
||||||
|
public :: &
|
||||||
IO_extractValue, &
|
IO_extractValue, &
|
||||||
IO_countDataLines
|
IO_countDataLines
|
||||||
#elif defined(Marc4DAMASK)
|
|
||||||
IO_fixedNoEFloatValue, &
|
|
||||||
IO_fixedIntValue
|
|
||||||
#endif
|
|
||||||
#endif
|
#endif
|
||||||
|
|
||||||
contains
|
contains
|
||||||
|
@ -57,7 +54,7 @@ contains
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
subroutine IO_init
|
subroutine IO_init
|
||||||
|
|
||||||
write(6,'(/,a)') ' <<<+- IO init -+>>>'
|
write(6,'(/,a)') ' <<<+- IO init -+>>>'; flush(6)
|
||||||
|
|
||||||
end subroutine IO_init
|
end subroutine IO_init
|
||||||
|
|
||||||
|
@ -444,57 +441,6 @@ integer function IO_intValue(string,chunkPos,myChunk)
|
||||||
end function IO_intValue
|
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
|
!> @brief changes characters in string to lower case
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
|
@ -916,62 +862,6 @@ end function IO_countDataLines
|
||||||
#endif
|
#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.
|
!> @brief return integer list corresponding to items in consecutive lines.
|
||||||
!! First integer in array is counter
|
!! First integer in array is counter
|
||||||
|
|
|
@ -51,7 +51,7 @@ subroutine discretization_init(homogenizationAt,microstructureAt,&
|
||||||
integer, optional, intent(in) :: &
|
integer, optional, intent(in) :: &
|
||||||
sharedNodesBeginn
|
sharedNodesBeginn
|
||||||
|
|
||||||
write(6,'(/,a)') ' <<<+- discretization init -+>>>'
|
write(6,'(/,a)') ' <<<+- discretization init -+>>>'; flush(6)
|
||||||
|
|
||||||
discretization_nElem = size(microstructureAt,1)
|
discretization_nElem = size(microstructureAt,1)
|
||||||
discretization_nIP = size(IPcoords0,2)/discretization_nElem
|
discretization_nIP = size(IPcoords0,2)/discretization_nElem
|
||||||
|
@ -81,15 +81,15 @@ subroutine discretization_results
|
||||||
|
|
||||||
real(pReal), dimension(:,:), allocatable :: u
|
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) &
|
u = discretization_NodeCoords (1:3,:discretization_sharedNodesBeginn) &
|
||||||
- discretization_NodeCoords0(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 &
|
u = discretization_IPcoords &
|
||||||
- discretization_IPcoords0
|
- 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
|
end subroutine discretization_results
|
||||||
|
|
||||||
|
|
336
src/element.f90
336
src/element.f90
|
@ -22,6 +22,8 @@ module element
|
||||||
NcellnodesPerCell, &
|
NcellnodesPerCell, &
|
||||||
nIPs, &
|
nIPs, &
|
||||||
nIPneighbors
|
nIPneighbors
|
||||||
|
character(len=:), allocatable :: &
|
||||||
|
vtkType
|
||||||
integer, dimension(:,:), allocatable :: &
|
integer, dimension(:,:), allocatable :: &
|
||||||
Cell, & !< intra-element (cell) nodes that constitute a cell
|
Cell, & !< intra-element (cell) nodes that constitute a cell
|
||||||
IPneighbor, &
|
IPneighbor, &
|
||||||
|
@ -36,10 +38,10 @@ module element
|
||||||
end type tElement
|
end type tElement
|
||||||
|
|
||||||
|
|
||||||
integer, parameter, private :: &
|
integer, parameter :: &
|
||||||
NELEMTYPE = 13
|
NELEMTYPE = 13
|
||||||
|
|
||||||
integer, dimension(NELEMTYPE), parameter, private :: NNODE = &
|
integer, dimension(NELEMTYPE), parameter :: NNODE = &
|
||||||
[ &
|
[ &
|
||||||
3, & ! 2D 3node 1ip
|
3, & ! 2D 3node 1ip
|
||||||
6, & ! 2D 6node 3ip
|
6, & ! 2D 6node 3ip
|
||||||
|
@ -57,7 +59,7 @@ module element
|
||||||
20 & ! 3D 20node 27ip
|
20 & ! 3D 20node 27ip
|
||||||
] !< number of nodes that constitute a specific type of element
|
] !< number of nodes that constitute a specific type of element
|
||||||
|
|
||||||
integer, dimension(NELEMTYPE), parameter, public :: GEOMTYPE = &
|
integer, dimension(NELEMTYPE), parameter :: GEOMTYPE = &
|
||||||
[ &
|
[ &
|
||||||
1, &
|
1, &
|
||||||
2, &
|
2, &
|
||||||
|
@ -74,7 +76,7 @@ module element
|
||||||
10 &
|
10 &
|
||||||
] !< geometry type of particular element type
|
] !< geometry type of particular element type
|
||||||
|
|
||||||
integer, dimension(maxval(GEOMTYPE)), parameter, private :: NCELLNODE = &
|
integer, dimension(maxval(GEOMTYPE)), parameter :: NCELLNODE = &
|
||||||
[ &
|
[ &
|
||||||
3, &
|
3, &
|
||||||
7, &
|
7, &
|
||||||
|
@ -88,7 +90,7 @@ module element
|
||||||
64 &
|
64 &
|
||||||
] !< number of cell nodes in a specific geometry type
|
] !< number of cell nodes in a specific geometry type
|
||||||
|
|
||||||
integer, dimension(maxval(GEOMTYPE)), parameter, private :: NIP = &
|
integer, dimension(maxval(GEOMTYPE)), parameter :: NIP = &
|
||||||
[ &
|
[ &
|
||||||
1, &
|
1, &
|
||||||
3, &
|
3, &
|
||||||
|
@ -102,7 +104,7 @@ module element
|
||||||
27 &
|
27 &
|
||||||
] !< number of IPs in a specific geometry type
|
] !< number of IPs in a specific geometry type
|
||||||
|
|
||||||
integer, dimension(maxval(GEOMTYPE)), parameter, private :: CELLTYPE = &
|
integer, dimension(maxval(GEOMTYPE)), parameter :: CELLTYPE = &
|
||||||
[ &
|
[ &
|
||||||
1, & ! 2D 3node
|
1, & ! 2D 3node
|
||||||
2, & ! 2D 4node
|
2, & ! 2D 4node
|
||||||
|
@ -116,7 +118,7 @@ module element
|
||||||
4 & ! 3D 8node
|
4 & ! 3D 8node
|
||||||
] !< cell type that is used by each geometry type
|
] !< 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
|
3, & ! 2D 3node
|
||||||
4, & ! 2D 4node
|
4, & ! 2D 4node
|
||||||
|
@ -124,7 +126,7 @@ module element
|
||||||
6 & ! 3D 8node
|
6 & ! 3D 8node
|
||||||
] !< number of ip neighbors / cell faces in a specific cell type
|
] !< 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 3node
|
||||||
2, & ! 2D 4node
|
2, & ! 2D 4node
|
||||||
|
@ -132,7 +134,7 @@ module element
|
||||||
4 & ! 3D 8node
|
4 & ! 3D 8node
|
||||||
] !< number of cell nodes in a specific cell type
|
] !< 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
|
3, & ! 2D 3node
|
||||||
4, & ! 2D 4node
|
4, & ! 2D 4node
|
||||||
|
@ -146,27 +148,39 @@ module element
|
||||||
! Positive integers denote an intra-element IP identifier.
|
! Positive integers denote an intra-element IP identifier.
|
||||||
! Negative integers denote the interface behind which the neighboring (extra-element) IP will be located.
|
! 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([&
|
reshape([&
|
||||||
-2,-3,-1 &
|
-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([&
|
reshape([&
|
||||||
2,-3, 3,-1, &
|
2,-3, 3,-1, &
|
||||||
-2, 1, 3,-1, &
|
-2, 1, 3,-1, &
|
||||||
2,-3,-2, 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([&
|
reshape([&
|
||||||
2,-4, 3,-1, &
|
2,-4, 3,-1, &
|
||||||
-2, 1, 4,-1, &
|
-2, 1, 4,-1, &
|
||||||
4,-4,-3, 1, &
|
4,-4,-3, 1, &
|
||||||
-2, 3,-3, 2 &
|
-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([&
|
reshape([&
|
||||||
2,-4, 4,-1, &
|
2,-4, 4,-1, &
|
||||||
3, 1, 5,-1, &
|
3, 1, 5,-1, &
|
||||||
|
@ -177,22 +191,34 @@ module element
|
||||||
8,-4,-3, 4, &
|
8,-4,-3, 4, &
|
||||||
9, 7,-3, 5, &
|
9, 7,-3, 5, &
|
||||||
-2, 8,-3, 6 &
|
-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([&
|
reshape([&
|
||||||
-1,-2,-3,-4 &
|
-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([&
|
reshape([&
|
||||||
2,-4, 3,-2, 4,-1, &
|
2,-4, 3,-2, 4,-1, &
|
||||||
-2, 1, 3,-2, 4,-1, &
|
-2, 1, 3,-2, 4,-1, &
|
||||||
2,-4,-3, 1, 4,-1, &
|
2,-4,-3, 1, 4,-1, &
|
||||||
2,-4, 3,-2,-3, 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([&
|
reshape([&
|
||||||
2,-4, 3,-2, 4,-1, &
|
2,-4, 3,-2, 4,-1, &
|
||||||
-3, 1, 3,-2, 5,-1, &
|
-3, 1, 3,-2, 5,-1, &
|
||||||
|
@ -200,14 +226,22 @@ module element
|
||||||
5,-4, 6,-2,-5, 1, &
|
5,-4, 6,-2,-5, 1, &
|
||||||
-3, 4, 6,-2,-5, 2, &
|
-3, 4, 6,-2,-5, 2, &
|
||||||
5,-4,-3, 4,-5, 3 &
|
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([&
|
reshape([&
|
||||||
-3,-5,-4,-2,-6,-1 &
|
-3,-5,-4,-2,-6,-1 &
|
||||||
],[nIPneighbor(cellType(8)),nIP(8)])
|
#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([&
|
reshape([&
|
||||||
2,-5, 3,-2, 5,-1, &
|
2,-5, 3,-2, 5,-1, &
|
||||||
-3, 1, 4,-2, 6,-1, &
|
-3, 1, 4,-2, 6,-1, &
|
||||||
|
@ -217,9 +251,13 @@ module element
|
||||||
-3, 5, 8,-2,-6, 2, &
|
-3, 5, 8,-2,-6, 2, &
|
||||||
8,-5,-4, 5,-6, 3, &
|
8,-5,-4, 5,-6, 3, &
|
||||||
-3, 7,-4, 6,-6, 4 &
|
-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([&
|
reshape([&
|
||||||
2,-5, 4,-2,10,-1, &
|
2,-5, 4,-2,10,-1, &
|
||||||
3, 1, 5,-2,11,-1, &
|
3, 1, 5,-2,11,-1, &
|
||||||
|
@ -248,17 +286,25 @@ module element
|
||||||
26,-5,-4,22,-6,16, &
|
26,-5,-4,22,-6,16, &
|
||||||
27,25,-4,23,-6,17, &
|
27,25,-4,23,-6,17, &
|
||||||
-3,26,-4,24,-6,18 &
|
-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([&
|
reshape([&
|
||||||
1, 0, 0, &
|
1, 0, 0, &
|
||||||
0, 1, 0, &
|
0, 1, 0, &
|
||||||
0, 0, 1 &
|
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([&
|
reshape([&
|
||||||
1, 0, 0, 0, 0, 0, &
|
1, 0, 0, 0, 0, 0, &
|
||||||
0, 1, 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, 1, 0, &
|
||||||
0, 0, 0, 0, 0, 1, &
|
0, 0, 0, 0, 0, 1, &
|
||||||
1, 1, 1, 2, 2, 2 &
|
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([&
|
reshape([&
|
||||||
1, 0, 0, 0, &
|
1, 0, 0, 0, &
|
||||||
0, 1, 0, 0, &
|
0, 1, 0, 0, &
|
||||||
|
@ -280,9 +330,13 @@ module element
|
||||||
0, 0, 1, 1, &
|
0, 0, 1, 1, &
|
||||||
1, 0, 0, 1, &
|
1, 0, 0, 1, &
|
||||||
1, 1, 1, 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([&
|
reshape([&
|
||||||
1, 0, 0, 0, 0, 0, 0, 0, &
|
1, 0, 0, 0, 0, 0, 0, 0, &
|
||||||
0, 1, 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, 4, 1, 1, 8, 8, 2, 2, &
|
||||||
1, 1, 4, 1, 2, 8, 8, 2, &
|
1, 1, 4, 1, 2, 8, 8, 2, &
|
||||||
1, 1, 1, 4, 2, 2, 8, 8 &
|
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([&
|
reshape([&
|
||||||
1, 0, 0, 0, 0, 0, 0, 0, &
|
1, 0, 0, 0, 0, 0, 0, 0, &
|
||||||
0, 1, 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, 1, 0, &
|
||||||
0, 0, 0, 0, 0, 0, 0, 1, &
|
0, 0, 0, 0, 0, 0, 0, 1, &
|
||||||
1, 1, 1, 1, 2, 2, 2, 2 &
|
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([&
|
reshape([&
|
||||||
1, 0, 0, 0, &
|
1, 0, 0, 0, &
|
||||||
0, 1, 0, 0, &
|
0, 1, 0, 0, &
|
||||||
0, 0, 1, 0, &
|
0, 0, 1, 0, &
|
||||||
0, 0, 0, 1 &
|
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([&
|
reshape([&
|
||||||
1, 0, 0, 0, 0, &
|
1, 0, 0, 0, 0, &
|
||||||
0, 1, 0, 0, 0, &
|
0, 1, 0, 0, 0, &
|
||||||
|
@ -340,9 +406,13 @@ module element
|
||||||
0, 1, 1, 1, 0, &
|
0, 1, 1, 1, 0, &
|
||||||
1, 0, 1, 1, 0, &
|
1, 0, 1, 1, 0, &
|
||||||
0, 0, 0, 0, 1 &
|
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([&
|
reshape([&
|
||||||
1, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
|
1, 0, 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, &
|
||||||
|
@ -359,9 +429,13 @@ module element
|
||||||
0, 1, 1, 1, 0, 2, 0, 0, 2, 2, &
|
0, 1, 1, 1, 0, 2, 0, 0, 2, 2, &
|
||||||
1, 0, 1, 1, 0, 0, 2, 2, 0, 2, &
|
1, 0, 1, 1, 0, 0, 2, 2, 0, 2, &
|
||||||
3, 3, 3, 3, 4, 4, 4, 4, 4, 4 &
|
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([&
|
reshape([&
|
||||||
1, 0, 0, 0, 0, 0, &
|
1, 0, 0, 0, 0, 0, &
|
||||||
0, 1, 0, 0, 0, 0, &
|
0, 1, 0, 0, 0, 0, &
|
||||||
|
@ -384,9 +458,13 @@ module element
|
||||||
1, 0, 1, 1, 0, 1, &
|
1, 0, 1, 1, 0, 1, &
|
||||||
0, 0, 0, 1, 1, 1, &
|
0, 0, 0, 1, 1, 1, &
|
||||||
1, 1, 1, 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([&
|
reshape([&
|
||||||
1, 0, 0, 0, 0, 0, 0, 0, &
|
1, 0, 0, 0, 0, 0, 0, 0, &
|
||||||
0, 1, 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, 1, 0, 0, &
|
||||||
0, 0, 0, 0, 0, 0, 1, 0, &
|
0, 0, 0, 0, 0, 0, 1, 0, &
|
||||||
0, 0, 0, 0, 0, 0, 0, 1 &
|
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([&
|
reshape([&
|
||||||
1, 0, 0, 0, 0, 0, 0, 0, & !
|
1, 0, 0, 0, 0, 0, 0, 0, & !
|
||||||
0, 1, 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
|
1, 0, 0, 1, 1, 0, 0, 1, & ! 25
|
||||||
0, 0, 0, 0, 1, 1, 1, 1, & !
|
0, 0, 0, 0, 1, 1, 1, 1, & !
|
||||||
1, 1, 1, 1, 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([&
|
reshape([&
|
||||||
1, 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, 0, & !
|
||||||
0, 1, 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
|
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, & !
|
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 & !
|
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([&
|
reshape([&
|
||||||
1, 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, 0, & !
|
||||||
0, 1, 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,20 +616,32 @@ module element
|
||||||
4, 8, 4, 3, 8,24, 8, 4, 12,12, 4, 4, 32,32,12,12, 12,32,12, 4, & !
|
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, & !
|
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 & !
|
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([&
|
reshape([&
|
||||||
1,2,3 &
|
1,2,3 &
|
||||||
|
#if !defined(__GFORTRAN__)
|
||||||
|
],shape(CELL1))
|
||||||
|
#else
|
||||||
],[NCELLNODEPERCELL(CELLTYPE(1)),NIP(1)])
|
],[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([&
|
reshape([&
|
||||||
1, 4, 7, 6, &
|
1, 4, 7, 6, &
|
||||||
2, 5, 7, 4, &
|
2, 5, 7, 4, &
|
||||||
3, 6, 7, 5 &
|
3, 6, 7, 5 &
|
||||||
|
#if !defined(__GFORTRAN__)
|
||||||
|
],shape(CELL2))
|
||||||
|
#else
|
||||||
],[NCELLNODEPERCELL(CELLTYPE(2)),NIP(2)])
|
],[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([&
|
reshape([&
|
||||||
|
@ -547,7 +649,11 @@ module element
|
||||||
5, 2, 6, 9, &
|
5, 2, 6, 9, &
|
||||||
8, 9, 7, 4, &
|
8, 9, 7, 4, &
|
||||||
9, 6, 3, 7 &
|
9, 6, 3, 7 &
|
||||||
|
#if !defined(__GFORTRAN__)
|
||||||
|
],shape(CELL3))
|
||||||
|
#else
|
||||||
],[NCELLNODEPERCELL(CELLTYPE(3)),NIP(3)])
|
],[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([&
|
reshape([&
|
||||||
|
@ -560,12 +666,20 @@ module element
|
||||||
11,16,10, 4, &
|
11,16,10, 4, &
|
||||||
16,15, 9,10, &
|
16,15, 9,10, &
|
||||||
15, 8, 3, 9 &
|
15, 8, 3, 9 &
|
||||||
|
#if !defined(__GFORTRAN__)
|
||||||
|
],shape(CELL4))
|
||||||
|
#else
|
||||||
],[NCELLNODEPERCELL(CELLTYPE(4)),NIP(4)])
|
],[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([&
|
reshape([&
|
||||||
1, 2, 3, 4 &
|
1, 2, 3, 4 &
|
||||||
|
#if !defined(__GFORTRAN__)
|
||||||
|
],shape(CELL5))
|
||||||
|
#else
|
||||||
],[NCELLNODEPERCELL(CELLTYPE(5)),NIP(5)])
|
],[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([&
|
reshape([&
|
||||||
|
@ -573,7 +687,11 @@ module element
|
||||||
5, 2, 6,11,12, 9,13,15, &
|
5, 2, 6,11,12, 9,13,15, &
|
||||||
7,11, 6, 3,14,15,13,10, &
|
7,11, 6, 3,14,15,13,10, &
|
||||||
8,12,15, 4, 4, 9,13,10 &
|
8,12,15, 4, 4, 9,13,10 &
|
||||||
|
#if !defined(__GFORTRAN__)
|
||||||
|
],shape(CELL6))
|
||||||
|
#else
|
||||||
],[NCELLNODEPERCELL(CELLTYPE(6)),NIP(6)])
|
],[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([&
|
reshape([&
|
||||||
|
@ -583,12 +701,20 @@ module element
|
||||||
10,17,21,19, 4,13,20,15, &
|
10,17,21,19, 4,13,20,15, &
|
||||||
17,11,18,21,13, 5,14,20, &
|
17,11,18,21,13, 5,14,20, &
|
||||||
19,21,18,12,15,20,14, 6 &
|
19,21,18,12,15,20,14, 6 &
|
||||||
|
#if !defined(__GFORTRAN__)
|
||||||
|
],shape(CELL7))
|
||||||
|
#else
|
||||||
],[NCELLNODEPERCELL(CELLTYPE(7)),NIP(7)])
|
],[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([&
|
reshape([&
|
||||||
1, 2, 3, 4, 5, 6, 7, 8 &
|
1, 2, 3, 4, 5, 6, 7, 8 &
|
||||||
|
#if !defined(__GFORTRAN__)
|
||||||
|
],shape(CELL8))
|
||||||
|
#else
|
||||||
],[NCELLNODEPERCELL(CELLTYPE(8)),NIP(8)])
|
],[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([&
|
reshape([&
|
||||||
|
@ -600,7 +726,11 @@ module element
|
||||||
22,14,23,27,17, 6,18,26, &
|
22,14,23,27,17, 6,18,26, &
|
||||||
25,27,24,16,20,26,19, 8, &
|
25,27,24,16,20,26,19, 8, &
|
||||||
27,23,15,24,26,18, 7,19 &
|
27,23,15,24,26,18, 7,19 &
|
||||||
|
#if !defined(__GFORTRAN__)
|
||||||
|
],shape(CELL9))
|
||||||
|
#else
|
||||||
],[NCELLNODEPERCELL(CELLTYPE(9)),NIP(9)])
|
],[NCELLNODEPERCELL(CELLTYPE(9)),NIP(9)])
|
||||||
|
#endif
|
||||||
|
|
||||||
integer, dimension(NCELLNODEPERCELL(CELLTYPE(10)),NIP(10)), parameter :: CELL10 = &
|
integer, dimension(NCELLNODEPERCELL(CELLTYPE(10)),NIP(10)), parameter :: CELL10 = &
|
||||||
reshape([&
|
reshape([&
|
||||||
|
@ -631,7 +761,11 @@ module element
|
||||||
51,64,50,24,31,56,30, 8, &
|
51,64,50,24,31,56,30, 8, &
|
||||||
64,63,49,50,56,55,29,30, &
|
64,63,49,50,56,55,29,30, &
|
||||||
63,48,23,49,55,28, 7,29 &
|
63,48,23,49,55,28, 7,29 &
|
||||||
|
#if !defined(__GFORTRAN__)
|
||||||
|
],shape(CELL10))
|
||||||
|
#else
|
||||||
],[NCELLNODEPERCELL(CELLTYPE(10)),NIP(10)])
|
],[NCELLNODEPERCELL(CELLTYPE(10)),NIP(10)])
|
||||||
|
#endif
|
||||||
|
|
||||||
|
|
||||||
integer, dimension(NCELLNODEPERCELLFACE(1),NIPNEIGHBOR(1)), parameter :: CELLFACE1 = &
|
integer, dimension(NCELLNODEPERCELLFACE(1),NIPNEIGHBOR(1)), parameter :: CELLFACE1 = &
|
||||||
|
@ -639,7 +773,11 @@ module element
|
||||||
2,3, &
|
2,3, &
|
||||||
3,1, &
|
3,1, &
|
||||||
1,2 &
|
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 = &
|
integer, dimension(NCELLNODEPERCELLFACE(2),NIPNEIGHBOR(2)), parameter :: CELLFACE2 = &
|
||||||
reshape([&
|
reshape([&
|
||||||
|
@ -647,7 +785,11 @@ module element
|
||||||
4,1, &
|
4,1, &
|
||||||
3,4, &
|
3,4, &
|
||||||
1,2 &
|
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 = &
|
integer, dimension(NCELLNODEPERCELLFACE(3),NIPNEIGHBOR(3)), parameter :: CELLFACE3 = &
|
||||||
reshape([&
|
reshape([&
|
||||||
|
@ -655,7 +797,11 @@ module element
|
||||||
1,2,4, &
|
1,2,4, &
|
||||||
2,3,4, &
|
2,3,4, &
|
||||||
1,4,3 &
|
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 = &
|
integer, dimension(NCELLNODEPERCELLFACE(4),NIPNEIGHBOR(4)), parameter :: CELLFACE4 = &
|
||||||
reshape([&
|
reshape([&
|
||||||
|
@ -665,85 +811,95 @@ module element
|
||||||
1,2,6,5, &
|
1,2,6,5, &
|
||||||
5,6,7,8, &
|
5,6,7,8, &
|
||||||
1,4,3,2 &
|
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
|
contains
|
||||||
|
|
||||||
|
|
||||||
|
!---------------------------------------------------------------------------------------------------
|
||||||
|
!> define properties of an element
|
||||||
|
!---------------------------------------------------------------------------------------------------
|
||||||
subroutine tElement_init(self,elemType)
|
subroutine tElement_init(self,elemType)
|
||||||
|
|
||||||
class(tElement) :: self
|
class(tElement) :: self
|
||||||
integer, intent(in) :: elemType
|
integer, intent(in) :: elemType
|
||||||
|
|
||||||
self%elemType = elemType
|
self%elemType = elemType
|
||||||
|
|
||||||
self%Nnodes = Nnode (self%elemType)
|
self%Nnodes = NNODE (self%elemType)
|
||||||
self%geomType = geomType (self%elemType)
|
self%geomType = GEOMTYPE(self%elemType)
|
||||||
|
|
||||||
select case (self%elemType)
|
select case (self%elemType)
|
||||||
case(1)
|
case(1)
|
||||||
self%cellNodeParentNodeWeights = cellNodeParentNodeWeights1
|
self%cellNodeParentNodeWeights = CELLNODEPARENTNODEWEIGHTS1
|
||||||
case(2)
|
case(2)
|
||||||
self%cellNodeParentNodeWeights = cellNodeParentNodeWeights2
|
self%cellNodeParentNodeWeights = CELLNODEPARENTNODEWEIGHTS2
|
||||||
case(3)
|
case(3)
|
||||||
self%cellNodeParentNodeWeights = cellNodeParentNodeWeights3
|
self%cellNodeParentNodeWeights = CELLNODEPARENTNODEWEIGHTS3
|
||||||
case(4)
|
case(4)
|
||||||
self%cellNodeParentNodeWeights = cellNodeParentNodeWeights4
|
self%cellNodeParentNodeWeights = CELLNODEPARENTNODEWEIGHTS4
|
||||||
case(5)
|
case(5)
|
||||||
self%cellNodeParentNodeWeights = cellNodeParentNodeWeights5
|
self%cellNodeParentNodeWeights = CELLNODEPARENTNODEWEIGHTS5
|
||||||
case(6)
|
case(6)
|
||||||
self%cellNodeParentNodeWeights = cellNodeParentNodeWeights6
|
self%cellNodeParentNodeWeights = CELLNODEPARENTNODEWEIGHTS6
|
||||||
case(7)
|
case(7)
|
||||||
self%cellNodeParentNodeWeights = cellNodeParentNodeWeights7
|
self%cellNodeParentNodeWeights = CELLNODEPARENTNODEWEIGHTS7
|
||||||
case(8)
|
case(8)
|
||||||
self%cellNodeParentNodeWeights = cellNodeParentNodeWeights8
|
self%cellNodeParentNodeWeights = CELLNODEPARENTNODEWEIGHTS8
|
||||||
case(9)
|
case(9)
|
||||||
self%cellNodeParentNodeWeights = cellNodeParentNodeWeights9
|
self%cellNodeParentNodeWeights = CELLNODEPARENTNODEWEIGHTS9
|
||||||
case(10)
|
case(10)
|
||||||
self%cellNodeParentNodeWeights = cellNodeParentNodeWeights10
|
self%cellNodeParentNodeWeights = CELLNODEPARENTNODEWEIGHTS10
|
||||||
case(11)
|
case(11)
|
||||||
self%cellNodeParentNodeWeights = cellNodeParentNodeWeights11
|
self%cellNodeParentNodeWeights = CELLNODEPARENTNODEWEIGHTS11
|
||||||
case(12)
|
case(12)
|
||||||
self%cellNodeParentNodeWeights = cellNodeParentNodeWeights12
|
self%cellNodeParentNodeWeights = CELLNODEPARENTNODEWEIGHTS12
|
||||||
case(13)
|
case(13)
|
||||||
self%cellNodeParentNodeWeights = cellNodeParentNodeWeights13
|
self%cellNodeParentNodeWeights = CELLNODEPARENTNODEWEIGHTS13
|
||||||
case default
|
case default
|
||||||
call IO_error(0,ext_msg='invalid element type')
|
call IO_error(0,ext_msg='invalid element type')
|
||||||
end select
|
end select
|
||||||
|
|
||||||
|
|
||||||
self%NcellNodes = NcellNode (self%geomType)
|
self%NcellNodes = NCELLNODE(self%geomType)
|
||||||
self%nIPs = nIP (self%geomType)
|
self%nIPs = NIP (self%geomType)
|
||||||
self%cellType = cellType (self%geomType)
|
self%cellType = CELLTYPE (self%geomType)
|
||||||
|
|
||||||
select case (self%geomType)
|
select case (self%geomType)
|
||||||
case(1)
|
case(1)
|
||||||
self%IPneighbor = IPneighbor1
|
self%IPneighbor = IPNEIGHBOR1
|
||||||
self%cell = CELL1
|
self%cell = CELL1
|
||||||
case(2)
|
case(2)
|
||||||
self%IPneighbor = IPneighbor2
|
self%IPneighbor = IPNEIGHBOR2
|
||||||
self%cell = CELL2
|
self%cell = CELL2
|
||||||
case(3)
|
case(3)
|
||||||
self%IPneighbor = IPneighbor3
|
self%IPneighbor = IPNEIGHBOR3
|
||||||
self%cell = CELL3
|
self%cell = CELL3
|
||||||
case(4)
|
case(4)
|
||||||
self%IPneighbor = IPneighbor4
|
self%IPneighbor = IPNEIGHBOR4
|
||||||
self%cell = CELL4
|
self%cell = CELL4
|
||||||
case(5)
|
case(5)
|
||||||
self%IPneighbor = IPneighbor5
|
self%IPneighbor = IPNEIGHBOR5
|
||||||
self%cell = CELL5
|
self%cell = CELL5
|
||||||
case(6)
|
case(6)
|
||||||
self%IPneighbor = IPneighbor6
|
self%IPneighbor = IPNEIGHBOR6
|
||||||
self%cell = CELL6
|
self%cell = CELL6
|
||||||
case(7)
|
case(7)
|
||||||
self%IPneighbor = IPneighbor7
|
self%IPneighbor = IPNEIGHBOR7
|
||||||
self%cell = CELL7
|
self%cell = CELL7
|
||||||
case(8)
|
case(8)
|
||||||
self%IPneighbor = IPneighbor8
|
self%IPneighbor = IPNEIGHBOR8
|
||||||
self%cell = CELL8
|
self%cell = CELL8
|
||||||
case(9)
|
case(9)
|
||||||
self%IPneighbor = IPneighbor9
|
self%IPneighbor = IPNEIGHBOR9
|
||||||
self%cell = CELL9
|
self%cell = CELL9
|
||||||
case(10)
|
case(10)
|
||||||
self%IPneighbor = IPneighbor10
|
self%IPneighbor = IPNEIGHBOR10
|
||||||
self%cell = CELL10
|
self%cell = CELL10
|
||||||
end select
|
end select
|
||||||
|
|
||||||
|
@ -752,17 +908,21 @@ subroutine tElement_init(self,elemType)
|
||||||
select case(self%cellType)
|
select case(self%cellType)
|
||||||
case(1)
|
case(1)
|
||||||
self%cellFace = CELLFACE1
|
self%cellFace = CELLFACE1
|
||||||
|
self%vtkType = 'TRIANGLE'
|
||||||
case(2)
|
case(2)
|
||||||
self%cellFace = CELLFACE2
|
self%cellFace = CELLFACE2
|
||||||
|
self%vtkType = 'QUAD'
|
||||||
case(3)
|
case(3)
|
||||||
self%cellFace = CELLFACE3
|
self%cellFace = CELLFACE3
|
||||||
|
self%vtkType = 'TETRA'
|
||||||
case(4)
|
case(4)
|
||||||
self%cellFace = CELLFACE4
|
self%cellFace = CELLFACE4
|
||||||
|
self%vtkType = 'HEXAHEDRON'
|
||||||
end select
|
end select
|
||||||
|
|
||||||
self%nIPneighbors = size(self%IPneighbor,1)
|
self%nIPneighbors = size(self%IPneighbor,1)
|
||||||
|
|
||||||
write(6,'(/,a)') ' <<<+- element_init -+>>>'
|
write(6,'(/,a)') ' <<<+- element_init -+>>>'; flush(6)
|
||||||
|
|
||||||
write(6,*) ' element type: ',self%elemType
|
write(6,*) ' element type: ',self%elemType
|
||||||
write(6,*) ' geom type: ',self%geomType
|
write(6,*) ' geom type: ',self%geomType
|
||||||
|
|
|
@ -72,7 +72,7 @@ subroutine mesh_init(ip,el)
|
||||||
real(pReal), dimension(:,:,:,:),allocatable :: &
|
real(pReal), dimension(:,:,:,:),allocatable :: &
|
||||||
unscaledNormals
|
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
|
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_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')
|
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
|
FEsolving_execElem = [1,nElems]
|
||||||
allocate(FEsolving_execIP(2,nElems), source=1) ! parallel loop bounds set to comprise from first IP...
|
FEsolving_execIP = spread([1,elem%nIPs],2,nElems)
|
||||||
FEsolving_execIP(2,:) = elem%nIPs
|
|
||||||
|
|
||||||
allocate(calcMode(elem%nIPs,nElems),source=.false.) ! pretend to have collected what first call is asking (F = I)
|
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"
|
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,&
|
call discretization_init(microstructureAt,homogenizationAt,&
|
||||||
ip_reshaped,&
|
ip_reshaped,&
|
||||||
node0_elem)
|
node0_cell)
|
||||||
|
|
||||||
call writeGeometry(0,connectivity_elem,&
|
call writeGeometry(elem,connectivity_elem,&
|
||||||
reshape(connectivity_cell,[elem%NcellNodesPerCell,elem%nIPs*nElems]),&
|
reshape(connectivity_cell,[elem%NcellNodesPerCell,elem%nIPs*nElems]),&
|
||||||
node0_cell,ip_reshaped)
|
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, &
|
connectivity_elem,connectivity_cell, &
|
||||||
coordinates_nodes,coordinates_points)
|
coordinates_nodes,coordinates_points)
|
||||||
|
|
||||||
integer, intent(in) :: elemType
|
type(tElement), intent(in) :: &
|
||||||
|
elem
|
||||||
integer, dimension(:,:), intent(in) :: &
|
integer, dimension(:,:), intent(in) :: &
|
||||||
connectivity_elem, &
|
connectivity_elem, &
|
||||||
connectivity_cell
|
connectivity_cell
|
||||||
|
@ -150,20 +150,24 @@ subroutine writeGeometry(elemType, &
|
||||||
connectivity_temp = connectivity_cell
|
connectivity_temp = connectivity_cell
|
||||||
call results_writeDataset('geometry',connectivity_temp,'T_c', &
|
call results_writeDataset('geometry',connectivity_temp,'T_c', &
|
||||||
'connectivity of the cells','-')
|
'connectivity of the cells','-')
|
||||||
|
call results_addAttribute('VTK_TYPE',elem%vtkType,'geometry/T_c')
|
||||||
|
|
||||||
coordinates_temp = coordinates_nodes
|
coordinates_temp = coordinates_nodes
|
||||||
call results_writeDataset('geometry',coordinates_temp,'x_n', &
|
call results_writeDataset('geometry',coordinates_temp,'x_n', &
|
||||||
'coordinates of the nodes','m')
|
'initial coordinates of the nodes','m')
|
||||||
|
|
||||||
coordinates_temp = coordinates_points
|
coordinates_temp = coordinates_points
|
||||||
call results_writeDataset('geometry',coordinates_temp,'x_p', &
|
call results_writeDataset('geometry',coordinates_temp,'x_p', &
|
||||||
'coordinates of the material points','m')
|
'initial coordinates of the materialpoints','m')
|
||||||
|
|
||||||
call results_closeJobFile
|
call results_closeJobFile
|
||||||
|
|
||||||
end subroutine writeGeometry
|
end subroutine writeGeometry
|
||||||
|
|
||||||
|
|
||||||
|
!--------------------------------------------------------------------------------------------------
|
||||||
|
!> @brief Read mesh from marc input file
|
||||||
|
!--------------------------------------------------------------------------------------------------
|
||||||
subroutine inputRead(elem,node0_elem,connectivity_elem,microstructureAt,homogenizationAt)
|
subroutine inputRead(elem,node0_elem,connectivity_elem,microstructureAt,homogenizationAt)
|
||||||
|
|
||||||
type(tElement), intent(out) :: elem
|
type(tElement), intent(out) :: elem
|
||||||
|
@ -171,7 +175,6 @@ subroutine inputRead(elem,node0_elem,connectivity_elem,microstructureAt,homogeni
|
||||||
node0_elem !< node x,y,z coordinates (initially!)
|
node0_elem !< node x,y,z coordinates (initially!)
|
||||||
integer, dimension(:,:), allocatable, intent(out) :: &
|
integer, dimension(:,:), allocatable, intent(out) :: &
|
||||||
connectivity_elem
|
connectivity_elem
|
||||||
|
|
||||||
integer, dimension(:), allocatable, intent(out) :: &
|
integer, dimension(:), allocatable, intent(out) :: &
|
||||||
microstructureAt, &
|
microstructureAt, &
|
||||||
homogenizationAt
|
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 IO_open_inputFile(FILEUNIT) ! ToDo: It would be better to use fileContent
|
||||||
|
|
||||||
call inputRead_mapElemSets(nameElemSet,mapElemSet,&
|
call inputRead_mapElemSets(nameElemSet,mapElemSet,&
|
||||||
FILEUNIT)
|
FILEUNIT,inputFile)
|
||||||
|
|
||||||
|
call inputRead_elemType(elem, &
|
||||||
|
nElems,inputFile)
|
||||||
|
|
||||||
allocate (mesh_mapFEtoCPelem(2,nElems), source=0)
|
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)
|
allocate (mesh_mapFEtoCPnode(2,Nnodes), source=0)
|
||||||
call inputRead_mapNodes(inputFile)
|
call inputRead_mapNodes(inputFile)
|
||||||
|
|
||||||
call inputRead_elemType(elem, &
|
|
||||||
nElems,inputFile)
|
|
||||||
call inputRead_elemNodes(node0_elem, &
|
call inputRead_elemNodes(node0_elem, &
|
||||||
Nnodes,inputFile)
|
Nnodes,inputFile)
|
||||||
|
|
||||||
|
@ -347,47 +352,63 @@ end subroutine inputRead_NnodesAndElements
|
||||||
!> @brief Count overall number of element sets in mesh.
|
!> @brief Count overall number of element sets in mesh.
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
subroutine inputRead_NelemSets(nElemSets,maxNelemInSet,&
|
subroutine inputRead_NelemSets(nElemSets,maxNelemInSet,&
|
||||||
fileUnit)
|
fileContent)
|
||||||
|
|
||||||
integer, intent(out) :: nElemSets, maxNelemInSet
|
integer, intent(out) :: nElemSets, maxNelemInSet
|
||||||
integer, intent(in) :: fileUnit
|
character(len=pStringLen), dimension(:), intent(in) :: fileContent !< file content, separated per lines
|
||||||
|
|
||||||
integer, allocatable, dimension(:) :: chunkPos
|
integer, allocatable, dimension(:) :: chunkPos
|
||||||
character(len=300) :: line
|
integer :: i,l,elemInCurrentSet
|
||||||
|
|
||||||
nElemSets = 0
|
nElemSets = 0
|
||||||
maxNelemInSet = 0
|
maxNelemInSet = 0
|
||||||
|
|
||||||
rewind(fileUnit)
|
do l = 1, size(fileContent)
|
||||||
do
|
chunkPos = IO_stringPos(fileContent(l))
|
||||||
read (fileUnit,'(A300)',END=620) line
|
if ( IO_lc(IO_StringValue(fileContent(l),chunkPos,1)) == 'define' .and. &
|
||||||
chunkPos = IO_stringPos(line)
|
IO_lc(IO_StringValue(fileContent(l),chunkPos,2)) == 'element' ) then
|
||||||
|
|
||||||
if ( IO_lc(IO_StringValue(line,chunkPos,1)) == 'define' .and. &
|
|
||||||
IO_lc(IO_StringValue(line,chunkPos,2)) == 'element' ) then
|
|
||||||
nElemSets = nElemSets + 1
|
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
|
endif
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
620 end subroutine inputRead_NelemSets
|
end subroutine inputRead_NelemSets
|
||||||
|
|
||||||
|
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
!> @brief map element sets
|
!> @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
|
character(len=64), dimension(:), allocatable, intent(out) :: nameElemSet
|
||||||
integer, dimension(:,:), allocatable, intent(out) :: mapElemSet
|
integer, dimension(:,:), allocatable, intent(out) :: mapElemSet
|
||||||
integer, intent(in) :: fileUnit
|
integer, intent(in) :: fileUnit
|
||||||
|
character(len=pStringLen), dimension(:), intent(in) :: fileContent !< file content, separated per lines
|
||||||
|
|
||||||
integer, allocatable, dimension(:) :: chunkPos
|
integer, allocatable, dimension(:) :: chunkPos
|
||||||
character(len=300) :: line
|
character(len=300) :: line
|
||||||
integer :: elemSet, NelemSets, maxNelemInSet
|
integer :: elemSet, NelemSets, maxNelemInSet
|
||||||
|
|
||||||
|
|
||||||
call inputRead_NelemSets(NelemSets,maxNelemInSet,fileUnit)
|
call inputRead_NelemSets(NelemSets,maxNelemInSet,fileContent)
|
||||||
allocate(nameElemSet(NelemSets)); nameElemSet = 'n/a'
|
allocate(nameElemSet(NelemSets)); nameElemSet = 'n/a'
|
||||||
allocate(mapElemSet(1+maxNelemInSet,NelemSets),source=0)
|
allocate(mapElemSet(1+maxNelemInSet,NelemSets),source=0)
|
||||||
elemSet = 0
|
elemSet = 0
|
||||||
|
@ -411,62 +432,32 @@ subroutine inputRead_mapElemSets(nameElemSet,mapElemSet,fileUnit)
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
!> @brief Maps elements from FE ID to internal (consecutive) representation.
|
!> @brief Maps elements from FE ID to internal (consecutive) representation.
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
subroutine inputRead_mapElems(tableStyle,nameElemSet,mapElemSet,fileFormatVersion,matNumber,fileUnit)
|
subroutine inputRead_mapElems(nNodes,nElem,fileContent)
|
||||||
|
|
||||||
integer, intent(in) :: fileUnit,tableStyle,fileFormatVersion
|
integer, intent(in) :: &
|
||||||
integer, dimension(:), intent(in) :: matNumber
|
nElem, &
|
||||||
character(len=64), intent(in), dimension(:) :: nameElemSet
|
nNodes !< number of nodes per element
|
||||||
integer, dimension(:,:), intent(in) :: &
|
character(len=pStringLen), dimension(:), intent(in) :: fileContent !< file content, separated per lines
|
||||||
mapElemSet
|
|
||||||
|
|
||||||
integer, allocatable, dimension(:) :: chunkPos
|
integer, allocatable, dimension(:) :: chunkPos
|
||||||
character(len=300) :: line, &
|
integer :: i,j,l,nNodesAlreadyRead
|
||||||
tmp
|
|
||||||
|
|
||||||
integer, dimension(:), allocatable :: contInts
|
do l = 1, size(fileContent)
|
||||||
integer :: i,cpElem
|
chunkPos = IO_stringPos(fileContent(l))
|
||||||
|
if( IO_lc(IO_stringValue(fileContent(l),chunkPos,1)) == 'connectivity' ) then
|
||||||
allocate(contInts(size(mesh_mapFEtoCPelem,2)+1))
|
j = 0
|
||||||
|
do i = 1,nElem
|
||||||
cpElem = 0
|
chunkPos = IO_stringPos(fileContent(l+1+i+j))
|
||||||
contInts = 0
|
mesh_mapFEtoCPelem(:,i) = [IO_intValue(fileContent(l+1+i+j),chunkPos,1),i]
|
||||||
rewind(fileUnit)
|
nNodesAlreadyRead = chunkPos(1) - 2
|
||||||
do
|
do while(nNodesAlreadyRead < nNodes) ! read on if not all nodes in one line
|
||||||
read (fileUnit,'(A300)',END=620) line
|
j = j + 1
|
||||||
chunkPos = IO_stringPos(line)
|
chunkPos = IO_stringPos(fileContent(l+1+i+j))
|
||||||
Marc2016andEarlier: if (fileFormatVersion < 13) then
|
nNodesAlreadyRead = nNodesAlreadyRead + chunkPos(1)
|
||||||
if( IO_lc(IO_stringValue(line,chunkPos,1)) == 'hypoelastic' ) then
|
enddo
|
||||||
skipLines: do i=1,3+TableStyle
|
enddo
|
||||||
read (fileUnit,'(A300)') line
|
|
||||||
enddo skipLines
|
|
||||||
contInts = IO_continuousIntValues(fileUnit,size(mesh_mapFEtoCPelem,2),nameElemSet,&
|
|
||||||
mapElemSet,size(nameElemSet))
|
|
||||||
exit
|
exit
|
||||||
endif
|
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
|
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
call math_sort(mesh_mapFEtoCPelem)
|
call math_sort(mesh_mapFEtoCPelem)
|
||||||
|
@ -488,7 +479,8 @@ subroutine inputRead_mapNodes(fileContent)
|
||||||
chunkPos = IO_stringPos(fileContent(l))
|
chunkPos = IO_stringPos(fileContent(l))
|
||||||
if( IO_lc(IO_stringValue(fileContent(l),chunkPos,1)) == 'coordinates' ) then
|
if( IO_lc(IO_stringValue(fileContent(l),chunkPos,1)) == 'coordinates' ) then
|
||||||
do i = 1,size(mesh_mapFEtoCPnode,2)
|
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
|
enddo
|
||||||
exit
|
exit
|
||||||
endif
|
endif
|
||||||
|
@ -509,7 +501,6 @@ subroutine inputRead_elemNodes(nodes, &
|
||||||
integer, intent(in) :: nNode
|
integer, intent(in) :: nNode
|
||||||
character(len=pStringLen), dimension(:), intent(in) :: fileContent !< file content, separated per lines
|
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, allocatable, dimension(:) :: chunkPos
|
||||||
integer :: i,j,m,l
|
integer :: i,j,m,l
|
||||||
|
|
||||||
|
@ -519,9 +510,10 @@ subroutine inputRead_elemNodes(nodes, &
|
||||||
chunkPos = IO_stringPos(fileContent(l))
|
chunkPos = IO_stringPos(fileContent(l))
|
||||||
if( IO_lc(IO_stringValue(fileContent(l),chunkPos,1)) == 'coordinates' ) then
|
if( IO_lc(IO_stringValue(fileContent(l),chunkPos,1)) == 'coordinates' ) then
|
||||||
do i=1,nNode
|
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
|
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
|
||||||
enddo
|
enddo
|
||||||
exit
|
exit
|
||||||
|
@ -709,7 +701,7 @@ subroutine inputRead_microstructureAndHomogenization(microstructureAt,homogeniza
|
||||||
read (fileUnit,'(A300)',END=630) line ! read line with value of state var
|
read (fileUnit,'(A300)',END=630) line ! read line with value of state var
|
||||||
chunkPos = IO_stringPos(line)
|
chunkPos = IO_stringPos(line)
|
||||||
do while (scan(IO_stringValue(line,chunkPos,1),'+-',back=.true.)>1) ! is noEfloat value?
|
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
|
if (initialcondTableStyle == 2) then
|
||||||
read (fileUnit,'(A300)',END=630) line ! read extra line
|
read (fileUnit,'(A300)',END=630) line ! read extra line
|
||||||
read (fileUnit,'(A300)',END=630) line ! read extra line
|
read (fileUnit,'(A300)',END=630) line ! read extra line
|
||||||
|
|
|
@ -48,10 +48,7 @@ module quaternions
|
||||||
procedure, private :: pow_scal__
|
procedure, private :: pow_scal__
|
||||||
generic, public :: operator(**) => pow_quat__, pow_scal__
|
generic, public :: operator(**) => pow_quat__, pow_scal__
|
||||||
|
|
||||||
procedure, public :: abs__
|
procedure, public :: abs => abs__
|
||||||
procedure, public :: dot_product__
|
|
||||||
procedure, public :: exp__
|
|
||||||
procedure, public :: log__
|
|
||||||
procedure, public :: conjg => conjg__
|
procedure, public :: conjg => conjg__
|
||||||
procedure, public :: real => real__
|
procedure, public :: real => real__
|
||||||
procedure, public :: aimag => aimag__
|
procedure, public :: aimag => aimag__
|
||||||
|
@ -317,17 +314,17 @@ end function pow_scal__
|
||||||
!---------------------------------------------------------------------------------------------------
|
!---------------------------------------------------------------------------------------------------
|
||||||
!> take exponential
|
!> 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
|
real(pReal) :: absImag
|
||||||
|
|
||||||
absImag = norm2(aimag(self))
|
absImag = norm2(aimag(a))
|
||||||
|
|
||||||
exp__ = merge(exp(self%w) * [ cos(absImag), &
|
exp__ = merge(exp(a%w) * [ cos(absImag), &
|
||||||
self%x/absImag * sin(absImag), &
|
a%x/absImag * sin(absImag), &
|
||||||
self%y/absImag * sin(absImag), &
|
a%y/absImag * sin(absImag), &
|
||||||
self%z/absImag * sin(absImag)], &
|
a%z/absImag * sin(absImag)], &
|
||||||
IEEE_value(1.0_pReal,IEEE_SIGNALING_NAN), &
|
IEEE_value(1.0_pReal,IEEE_SIGNALING_NAN), &
|
||||||
dNeq0(absImag))
|
dNeq0(absImag))
|
||||||
|
|
||||||
|
@ -337,17 +334,17 @@ end function exp__
|
||||||
!---------------------------------------------------------------------------------------------------
|
!---------------------------------------------------------------------------------------------------
|
||||||
!> take logarithm
|
!> 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
|
real(pReal) :: absImag
|
||||||
|
|
||||||
absImag = norm2(aimag(self))
|
absImag = norm2(aimag(a))
|
||||||
|
|
||||||
log__ = merge([log(abs(self)), &
|
log__ = merge([log(abs(a)), &
|
||||||
self%x/absImag * acos(self%w/abs(self)), &
|
a%x/absImag * acos(a%w/abs(a)), &
|
||||||
self%y/absImag * acos(self%w/abs(self)), &
|
a%y/absImag * acos(a%w/abs(a)), &
|
||||||
self%z/absImag * acos(self%w/abs(self))], &
|
a%z/absImag * acos(a%w/abs(a))], &
|
||||||
IEEE_value(1.0_pReal,IEEE_SIGNALING_NAN), &
|
IEEE_value(1.0_pReal,IEEE_SIGNALING_NAN), &
|
||||||
dNeq0(absImag))
|
dNeq0(absImag))
|
||||||
|
|
||||||
|
@ -357,11 +354,11 @@ end function log__
|
||||||
!---------------------------------------------------------------------------------------------------
|
!---------------------------------------------------------------------------------------------------
|
||||||
!> return norm
|
!> 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__
|
end function abs__
|
||||||
|
|
||||||
|
@ -381,11 +378,11 @@ end function dot_product__
|
||||||
!---------------------------------------------------------------------------------------------------
|
!---------------------------------------------------------------------------------------------------
|
||||||
!> take conjugate complex
|
!> 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__
|
end function conjg__
|
||||||
|
|
||||||
|
|
|
@ -70,7 +70,7 @@ subroutine results_init
|
||||||
|
|
||||||
resultsFile = HDF5_openFile(trim(getSolverJobName())//'.hdf5','w',.true.)
|
resultsFile = HDF5_openFile(trim(getSolverJobName())//'.hdf5','w',.true.)
|
||||||
call results_addAttribute('DADF5_version_major',0)
|
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 results_addAttribute('DAMASK_version',DAMASKVERSION)
|
||||||
call get_command(commandLine)
|
call get_command(commandLine)
|
||||||
call results_addAttribute('call',trim(commandLine))
|
call results_addAttribute('call',trim(commandLine))
|
||||||
|
|
|
@ -42,7 +42,8 @@
|
||||||
! Convention 4: Euler angle triplets are implemented using the Bunge convention,
|
! Convention 4: Euler angle triplets are implemented using the Bunge convention,
|
||||||
! with the angular ranges as [0, 2π],[0, π],[0, 2π]
|
! with the angular ranges as [0, 2π],[0, π],[0, 2π]
|
||||||
! Convention 5: the rotation angle ω is limited to the interval [0, π]
|
! 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
|
module rotations
|
||||||
|
@ -77,6 +78,7 @@ module rotations
|
||||||
procedure, public :: rotTensor4
|
procedure, public :: rotTensor4
|
||||||
procedure, public :: rotTensor4sym
|
procedure, public :: rotTensor4sym
|
||||||
procedure, public :: misorientation
|
procedure, public :: misorientation
|
||||||
|
procedure, public :: standardize
|
||||||
end type rotation
|
end type rotation
|
||||||
|
|
||||||
public :: &
|
public :: &
|
||||||
|
@ -92,7 +94,7 @@ contains
|
||||||
subroutine rotations_init
|
subroutine rotations_init
|
||||||
|
|
||||||
call quaternions_init
|
call quaternions_init
|
||||||
write(6,'(/,a)') ' <<<+- rotations init -+>>>'
|
write(6,'(/,a)') ' <<<+- rotations init -+>>>'; flush(6)
|
||||||
call unitTest
|
call unitTest
|
||||||
|
|
||||||
end subroutine rotations_init
|
end subroutine rotations_init
|
||||||
|
@ -237,7 +239,6 @@ end subroutine fromMatrix
|
||||||
|
|
||||||
!---------------------------------------------------------------------------------------------------
|
!---------------------------------------------------------------------------------------------------
|
||||||
!> @brief: Rotate a rotation
|
!> @brief: Rotate a rotation
|
||||||
!> ToDo: completly untested
|
|
||||||
!---------------------------------------------------------------------------------------------------
|
!---------------------------------------------------------------------------------------------------
|
||||||
pure elemental function rotRot__(self,R) result(rRot)
|
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
|
class(rotation), intent(in) :: self,R
|
||||||
|
|
||||||
rRot = rotation(self%q*R%q)
|
rRot = rotation(self%q*R%q)
|
||||||
|
call rRot%standardize()
|
||||||
|
|
||||||
end function rotRot__
|
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
|
!> @author Marc De Graef, Carnegie Mellon University
|
||||||
!> @brief rotate a vector passively (default) or actively
|
!> @brief rotate a vector passively (default) or actively
|
||||||
|
@ -375,7 +389,7 @@ pure elemental function misorientation(self,other)
|
||||||
type(rotation) :: misorientation
|
type(rotation) :: misorientation
|
||||||
class(rotation), intent(in) :: self, other
|
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
|
end function misorientation
|
||||||
|
|
||||||
|
|
Loading…
Reference in New Issue