Merge remote-tracking branch 'origin/development' into initial-eigenstrain

.gitlab-ci.yml

@@ -81,7 +81,7 @@ checkout:
     - release
 
 ###################################################################################################
-processing:
+pytest:
   stage: python
   script:
     - cd $DAMASKROOT/python
@@ -91,6 +91,15 @@ processing:
     - master
     - release
 
+mypy:
+  stage: python
+  script:
+    - cd $DAMASKROOT/python
+    - mypy -m damask
+  except:
+    - master
+    - release
+
 
 ###################################################################################################
 compile_grid_Intel:

PRIVATE

@@ -1 +1 @@
-Subproject commit fabe69749425e8a7aceb3b7c2758b40d97d8b809
+Subproject commit 5a769ec759d9dacc1866c35c6663cd0001e198c5

examples/config/Phase_Dislotungsten_W.yaml

@@ -1,26 +0,0 @@
-type: dislotungsten
-
-N_sl: [12]
-
-rho_mob_0: [1.0e+9]
-rho_dip_0: [1.0]
-
-nu_a: [9.1e+11]
-b_sl: [2.72e-10]
-Delta_H_kp,0: [2.61154e-19]  # 1.63 eV, Delta_H0
-
-tau_Peierls: [2.03e+9]
-p_sl: [0.86]
-q_sl: [1.69]
-h: [2.566e-10]
-w: [2.992e-09]
-B: [8.3e-5]
-D_a: 1.0 # d_edge
-
-# climb (disabled)
-D_0: 0.0
-Q_cl: 0.0
-V_cl: [0.0]
-
-h_sl-sl: [0.009, 0.72, 0.009, 0.05, 0.05, 0.06, 0.09]
-a_nonSchmid: [0.938, 0.71, 4.43]

examples/config/phase/mechanical/plastic/dislotungsten_W.yaml

@@ -0,0 +1,35 @@
+type: dislotungsten
+references:
+  - D. Cereceda et al.,
+    International Journal of Plasticity 78:242-265, 2016,
+    http://dx.doi.org/10.1016/j.ijplas.2015.09.002
+  - R. Gröger et al.,
+    Acta Materialia 56(19):5412-5425, 2008,
+    https://doi.org/10.1016/j.actamat.2008.07.037
+output: [Lambda_sl]
+N_sl: [12]
+b_sl: [2.72e-10]
+rho_mob_0: [1.0e+9]  # estimated from section 3.2
+rho_dip_0: [1.0]     # not given
+Q_s: [2.61154e-19]   # 1.63 eV, Delta_H0
+B: [8.3e-5]
+omega: [9.1e+11]     # nu_0
+p_sl: [0.86]
+q_sl: [1.69]
+tau_Peierls: [2.03e+9]
+h: [2.566e-10]
+h_sl-sl: [0.009, 0.72, 0.009, 0.05, 0.05, 0.06, 0.09]
+w: [2.992e-09] # 11b
+
+# values in Cereceda et al. are high, using parameters from Gröger et al.
+a_nonSchmid: [0.0, 0.56, 0.75] # Table 2
+
+# (almost) no annhilation, adjustment needed for simulations beyond the yield point
+i_sl: [1]  # c, eq. (25)
+D: 1.0e+20 # d_g, eq. (25)
+D_a: 1.0 # d_edge = D_a*b
+
+# disable climb (not discussed in Cereceda et al.)
+D_0: 0.0
+f_at: 1
+Q_cl: 1.0

python/damask/VERSION

@@ -1 +1 @@
-v3.0.0-alpha5-31-gddb25ad0e
+v3.0.0-alpha5-64-g8e08af31e

python/damask/_orientation.py

@@ -125,7 +125,7 @@ class Orientation(Rotation,Crystal):
         """Create deep copy."""
         dup = copy.deepcopy(self)
         if rotation is not None:
-            dup.quaternion = Orientation(rotation,family='cubic').quaternion
+            dup.quaternion = Rotation(rotation).quaternion
         return dup
 
     copy = __copy__

python/damask/_rotation.py

@@ -1,3 +1,5 @@
+import copy
+
 import numpy as np
 
 from . import tensor
@@ -85,9 +87,12 @@ class Rotation:
                + str(self.quaternion)
 
 
-    def __copy__(self,**kwargs):
+    def __copy__(self,rotation=None):
         """Create deep copy."""
-        return self.__class__(rotation=kwargs['rotation'] if 'rotation' in kwargs else self.quaternion)
+        dup = copy.deepcopy(self)
+        if rotation is not None:
+            dup.quaternion = Rotation(rotation).quaternion
+        return dup
 
     copy = __copy__
 

python/damask/grid_filters.py

@@ -11,10 +11,14 @@ the following operations are required for tensorial data:
     - D1 = D3.reshape(cells+(-1,)).reshape(-1,9,order='F')
 
 """
+
+from typing import Sequence, Tuple, Union
+
 from scipy import spatial as _spatial
 import numpy as _np
 
-def _ks(size,cells,first_order=False):
+
+def _ks(size: _np.ndarray, cells: Union[_np.ndarray,Sequence[int]], first_order: bool = False) -> _np.ndarray:
     """
     Get wave numbers operator.
 
@@ -41,7 +45,7 @@ def _ks(size,cells,first_order=False):
     return _np.stack(_np.meshgrid(k_sk,k_sj,k_si,indexing = 'ij'), axis=-1)
 
 
-def curl(size,f):
+def curl(size: _np.ndarray, f: _np.ndarray) -> _np.ndarray:
     u"""
     Calculate curl of a vector or tensor field in Fourier space.
 
@@ -72,7 +76,7 @@ def curl(size,f):
     return _np.fft.irfftn(curl_,axes=(0,1,2),s=f.shape[:3])
 
 
-def divergence(size,f):
+def divergence(size: _np.ndarray, f: _np.ndarray) -> _np.ndarray:
     u"""
     Calculate divergence of a vector or tensor field in Fourier space.
 
@@ -99,7 +103,7 @@ def divergence(size,f):
     return _np.fft.irfftn(div_,axes=(0,1,2),s=f.shape[:3])
 
 
-def gradient(size,f):
+def gradient(size: _np.ndarray, f: _np.ndarray) -> _np.ndarray:
     u"""
     Calculate gradient of a scalar or vector fieldin Fourier space.
 
@@ -126,7 +130,9 @@ def gradient(size,f):
     return _np.fft.irfftn(grad_,axes=(0,1,2),s=f.shape[:3])
 
 
-def coordinates0_point(cells,size,origin=_np.zeros(3)):
+def coordinates0_point(cells: Union[ _np.ndarray,Sequence[int]],
+                       size: _np.ndarray,
+                       origin: _np.ndarray = _np.zeros(3)) -> _np.ndarray:
     """
     Cell center positions (undeformed).
 
@@ -145,8 +151,8 @@ def coordinates0_point(cells,size,origin=_np.zeros(3)):
         Undeformed cell center coordinates.
 
     """
-    start = origin        + size/cells*.5
+    start = origin        + size/_np.array(cells)*.5
-    end   = origin + size - size/cells*.5
+    end   = origin + size - size/_np.array(cells)*.5
 
     return _np.stack(_np.meshgrid(_np.linspace(start[0],end[0],cells[0]),
                                   _np.linspace(start[1],end[1],cells[1]),
@@ -154,7 +160,7 @@ def coordinates0_point(cells,size,origin=_np.zeros(3)):
                      axis = -1)
 
 
-def displacement_fluct_point(size,F):
+def displacement_fluct_point(size: _np.ndarray, F: _np.ndarray) -> _np.ndarray:
     """
     Cell center displacement field from fluctuation part of the deformation gradient field.
 
@@ -186,7 +192,7 @@ def displacement_fluct_point(size,F):
     return _np.fft.irfftn(displacement,axes=(0,1,2),s=F.shape[:3])
 
 
-def displacement_avg_point(size,F):
+def displacement_avg_point(size: _np.ndarray, F: _np.ndarray) -> _np.ndarray:
     """
     Cell center displacement field from average part of the deformation gradient field.
 
@@ -207,7 +213,7 @@ def displacement_avg_point(size,F):
     return _np.einsum('ml,ijkl->ijkm',F_avg - _np.eye(3),coordinates0_point(F.shape[:3],size))
 
 
-def displacement_point(size,F):
+def displacement_point(size: _np.ndarray, F: _np.ndarray) -> _np.ndarray:
     """
     Cell center displacement field from deformation gradient field.
 
@@ -227,7 +233,7 @@ def displacement_point(size,F):
     return displacement_avg_point(size,F) + displacement_fluct_point(size,F)
 
 
-def coordinates_point(size,F,origin=_np.zeros(3)):
+def coordinates_point(size: _np.ndarray, F: _np.ndarray, origin: _np.ndarray = _np.zeros(3)) -> _np.ndarray:
     """
     Cell center positions.
 
@@ -249,7 +255,8 @@ def coordinates_point(size,F,origin=_np.zeros(3)):
     return coordinates0_point(F.shape[:3],size,origin) + displacement_point(size,F)
 
 
-def cellsSizeOrigin_coordinates0_point(coordinates0,ordered=True):
+def cellsSizeOrigin_coordinates0_point(coordinates0: _np.ndarray,
+                                       ordered: bool = True) -> Tuple[_np.ndarray,_np.ndarray,_np.ndarray]:
     """
     Return grid 'DNA', i.e. cells, size, and origin from 1D array of point positions.
 
@@ -292,13 +299,15 @@ def cellsSizeOrigin_coordinates0_point(coordinates0,ordered=True):
         raise ValueError('Regular cell spacing violated.')
 
     if ordered and not _np.allclose(coordinates0.reshape(tuple(cells)+(3,),order='F'),
-                                    coordinates0_point(cells,size,origin),atol=atol):
+                                    coordinates0_point(list(cells),size,origin),atol=atol):
         raise ValueError('Input data is not ordered (x fast, z slow).')
 
     return (cells,size,origin)
 
 
-def coordinates0_node(cells,size,origin=_np.zeros(3)):
+def coordinates0_node(cells: Union[_np.ndarray,Sequence[int]],
+                      size: _np.ndarray,
+                      origin: _np.ndarray = _np.zeros(3)) -> _np.ndarray:
     """
     Nodal positions (undeformed).
 
@@ -323,7 +332,7 @@ def coordinates0_node(cells,size,origin=_np.zeros(3)):
                      axis = -1)
 
 
-def displacement_fluct_node(size,F):
+def displacement_fluct_node(size: _np.ndarray, F: _np.ndarray) -> _np.ndarray:
     """
     Nodal displacement field from fluctuation part of the deformation gradient field.
 
@@ -343,7 +352,7 @@ def displacement_fluct_node(size,F):
     return point_to_node(displacement_fluct_point(size,F))
 
 
-def displacement_avg_node(size,F):
+def displacement_avg_node(size: _np.ndarray, F: _np.ndarray) -> _np.ndarray:
     """
     Nodal displacement field from average part of the deformation gradient field.
 
@@ -364,7 +373,7 @@ def displacement_avg_node(size,F):
     return _np.einsum('ml,ijkl->ijkm',F_avg - _np.eye(3),coordinates0_node(F.shape[:3],size))
 
 
-def displacement_node(size,F):
+def displacement_node(size: _np.ndarray, F: _np.ndarray) -> _np.ndarray:
     """
     Nodal displacement field from deformation gradient field.
 
@@ -384,7 +393,7 @@ def displacement_node(size,F):
     return displacement_avg_node(size,F) + displacement_fluct_node(size,F)
 
 
-def coordinates_node(size,F,origin=_np.zeros(3)):
+def coordinates_node(size: _np.ndarray, F: _np.ndarray, origin: _np.ndarray = _np.zeros(3)) -> _np.ndarray:
     """
     Nodal positions.
 
@@ -406,7 +415,8 @@ def coordinates_node(size,F,origin=_np.zeros(3)):
     return coordinates0_node(F.shape[:3],size,origin) + displacement_node(size,F)
 
 
-def cellsSizeOrigin_coordinates0_node(coordinates0,ordered=True):
+def cellsSizeOrigin_coordinates0_node(coordinates0: _np.ndarray,
+                                      ordered: bool = True) -> Tuple[_np.ndarray,_np.ndarray,_np.ndarray]:
     """
     Return grid 'DNA', i.e. cells, size, and origin from 1D array of nodal positions.
 
@@ -441,13 +451,13 @@ def cellsSizeOrigin_coordinates0_node(coordinates0,ordered=True):
         raise ValueError('Regular cell spacing violated.')
 
     if ordered and not _np.allclose(coordinates0.reshape(tuple(cells+1)+(3,),order='F'),
-                                    coordinates0_node(cells,size,origin),atol=atol):
+                                    coordinates0_node(list(cells),size,origin),atol=atol):
         raise ValueError('Input data is not ordered (x fast, z slow).')
 
     return (cells,size,origin)
 
 
-def point_to_node(cell_data):
+def point_to_node(cell_data: _np.ndarray) -> _np.ndarray:
     """
     Interpolate periodic point data to nodal data.
 
@@ -469,7 +479,7 @@ def point_to_node(cell_data):
     return _np.pad(n,((0,1),(0,1),(0,1))+((0,0),)*len(cell_data.shape[3:]),mode='wrap')
 
 
-def node_to_point(node_data):
+def node_to_point(node_data: _np.ndarray) -> _np.ndarray:
     """
     Interpolate periodic nodal data to point data.
 
@@ -491,7 +501,7 @@ def node_to_point(node_data):
     return c[1:,1:,1:]
 
 
-def coordinates0_valid(coordinates0):
+def coordinates0_valid(coordinates0: _np.ndarray) -> bool:
     """
     Check whether coordinates form a regular grid.
 
@@ -513,7 +523,7 @@ def coordinates0_valid(coordinates0):
         return False
 
 
-def regrid(size,F,cells):
+def regrid(size: _np.ndarray, F: _np.ndarray, cells: Union[_np.ndarray,Sequence[int]]) -> _np.ndarray:
     """
     Return mapping from coordinates in deformed configuration to a regular grid.
 

python/damask/mechanics.py

@@ -5,13 +5,15 @@ All routines operate on numpy.ndarrays of shape (...,3,3).
 
 """
 
-from . import tensor as _tensor
+from typing import Sequence
-from . import _rotation
 
 import numpy as _np
 
+from . import tensor as _tensor
+from . import _rotation
 
-def deformation_Cauchy_Green_left(F):
+
+def deformation_Cauchy_Green_left(F: _np.ndarray) -> _np.ndarray:
     """
     Calculate left Cauchy-Green deformation tensor (Finger deformation tensor).
 
@@ -29,7 +31,7 @@ def deformation_Cauchy_Green_left(F):
     return _np.matmul(F,_tensor.transpose(F))
 
 
-def deformation_Cauchy_Green_right(F):
+def deformation_Cauchy_Green_right(F: _np.ndarray) -> _np.ndarray:
     """
     Calculate right Cauchy-Green deformation tensor.
 
@@ -47,7 +49,7 @@ def deformation_Cauchy_Green_right(F):
     return _np.matmul(_tensor.transpose(F),F)
 
 
-def equivalent_strain_Mises(epsilon):
+def equivalent_strain_Mises(epsilon: _np.ndarray) -> _np.ndarray:
     """
     Calculate the Mises equivalent of a strain tensor.
 
@@ -65,7 +67,7 @@ def equivalent_strain_Mises(epsilon):
     return _equivalent_Mises(epsilon,2.0/3.0)
 
 
-def equivalent_stress_Mises(sigma):
+def equivalent_stress_Mises(sigma: _np.ndarray) -> _np.ndarray:
     """
     Calculate the Mises equivalent of a stress tensor.
 
@@ -83,7 +85,7 @@ def equivalent_stress_Mises(sigma):
     return _equivalent_Mises(sigma,3.0/2.0)
 
 
-def maximum_shear(T_sym):
+def maximum_shear(T_sym: _np.ndarray) -> _np.ndarray:
     """
     Calculate the maximum shear component of a symmetric tensor.
 
@@ -102,7 +104,7 @@ def maximum_shear(T_sym):
     return (w[...,0] - w[...,2])*0.5
 
 
-def rotation(T):
+def rotation(T: _np.ndarray) -> _rotation.Rotation:
     """
     Calculate the rotational part of a tensor.
 
@@ -120,7 +122,7 @@ def rotation(T):
     return _rotation.Rotation.from_matrix(_polar_decomposition(T,'R')[0])
 
 
-def strain(F,t,m):
+def strain(F: _np.ndarray, t: str, m: float) -> _np.ndarray:
     """
     Calculate strain tensor (Seth–Hill family).
 
@@ -160,7 +162,7 @@ def strain(F,t,m):
     return eps
 
 
-def stress_Cauchy(P,F):
+def stress_Cauchy(P: _np.ndarray, F: _np.ndarray) -> _np.ndarray:
     """
     Calculate the Cauchy stress (true stress).
 
@@ -182,7 +184,7 @@ def stress_Cauchy(P,F):
     return _tensor.symmetric(_np.einsum('...,...ij,...kj',1.0/_np.linalg.det(F),P,F))
 
 
-def stress_second_Piola_Kirchhoff(P,F):
+def stress_second_Piola_Kirchhoff(P: _np.ndarray, F: _np.ndarray) -> _np.ndarray:
     """
     Calculate the second Piola-Kirchhoff stress.
 
@@ -205,7 +207,7 @@ def stress_second_Piola_Kirchhoff(P,F):
     return _tensor.symmetric(_np.einsum('...ij,...jk',_np.linalg.inv(F),P))
 
 
-def stretch_left(T):
+def stretch_left(T: _np.ndarray) -> _np.ndarray:
     """
     Calculate left stretch of a tensor.
 
@@ -223,7 +225,7 @@ def stretch_left(T):
     return _polar_decomposition(T,'V')[0]
 
 
-def stretch_right(T):
+def stretch_right(T: _np.ndarray) -> _np.ndarray:
     """
     Calculate right stretch of a tensor.
 
@@ -241,7 +243,7 @@ def stretch_right(T):
     return _polar_decomposition(T,'U')[0]
 
 
-def _polar_decomposition(T,requested):
+def _polar_decomposition(T: _np.ndarray, requested: Sequence[str]) -> tuple:
     """
     Perform singular value decomposition.
 
@@ -257,21 +259,21 @@ def _polar_decomposition(T,requested):
     u, _, vh = _np.linalg.svd(T)
     R = _np.einsum('...ij,...jk',u,vh)
 
-    output = ()
+    output = []
     if 'R' in requested:
-        output+=(R,)
+        output+=[R]
     if 'V' in requested:
-        output+=(_np.einsum('...ij,...kj',T,R),)
+        output+=[_np.einsum('...ij,...kj',T,R)]
     if 'U' in requested:
-        output+=(_np.einsum('...ji,...jk',R,T),)
+        output+=[_np.einsum('...ji,...jk',R,T)]
 
     if len(output) == 0:
         raise ValueError('output needs to be out of V, R, U')
 
-    return output
+    return tuple(output)
 
 
-def _equivalent_Mises(T_sym,s):
+def _equivalent_Mises(T_sym: _np.ndarray, s: float) -> _np.ndarray:
     """
     Base equation for Mises equivalent of a stress or strain tensor.
 

python/damask/seeds.py

@@ -1,5 +1,7 @@
 """Functionality for generation of seed points for Voronoi or Laguerre tessellation."""
 
+from typing import Sequence,Tuple
+
 from scipy import spatial as _spatial
 import numpy as _np
 
@@ -7,7 +9,7 @@ from . import util as _util
 from . import grid_filters as _grid_filters
 
 
-def from_random(size,N_seeds,cells=None,rng_seed=None):
+def from_random(size: _np.ndarray, N_seeds: int, cells: _np.ndarray = None, rng_seed=None) -> _np.ndarray:
     """
     Place seeds randomly in space.
 
@@ -41,7 +43,8 @@ def from_random(size,N_seeds,cells=None,rng_seed=None):
     return coords
 
 
-def from_Poisson_disc(size,N_seeds,N_candidates,distance,periodic=True,rng_seed=None):
+def from_Poisson_disc(size: _np.ndarray, N_seeds: int, N_candidates: int, distance: float,
+                      periodic: bool = True, rng_seed=None) -> _np.ndarray:
     """
     Place seeds according to a Poisson disc distribution.
 
@@ -75,18 +78,17 @@ def from_Poisson_disc(size,N_seeds,N_candidates,distance,periodic=True,rng_seed=
     i = 0
     progress = _util._ProgressBar(N_seeds+1,'',50)
     while s < N_seeds:
+        i += 1
         candidates = rng.random((N_candidates,3))*_np.broadcast_to(size,(N_candidates,3))
         tree = _spatial.cKDTree(coords[:s],boxsize=size) if periodic else \
                _spatial.cKDTree(coords[:s])
         distances = tree.query(candidates)[0]
         best = distances.argmax()
         if distances[best] > distance:                                                              # require minimum separation
+            i = 0
             coords[s] = candidates[best]                                                            # maximum separation to existing point cloud
             s += 1
             progress.update(s)
-            i = 0
-        else:
-            i += 1
 
         if i == 100:
             raise ValueError('Seeding not possible')
@@ -94,22 +96,23 @@ def from_Poisson_disc(size,N_seeds,N_candidates,distance,periodic=True,rng_seed=
     return coords
 
 
-def from_grid(grid,selection=None,invert=False,average=False,periodic=True):
+def from_grid(grid, selection: Sequence[int] = None,
+              invert: bool = False, average: bool = False, periodic: bool = True) -> Tuple[_np.ndarray, _np.ndarray]:
     """
     Create seeds from grid description.
 
     Parameters
     ----------
     grid : damask.Grid
-        Grid, from which the material IDs are used as seeds.
+        Grid from which the material IDs are used as seeds.
     selection : iterable of integers, optional
         Material IDs to consider.
     invert : boolean, false
-        Do not consider the material IDs given in selection. Defaults to False.
+        Consider all material IDs except those in selection. Defaults to False.
     average : boolean, optional
         Seed corresponds to center of gravity of material ID cloud.
     periodic : boolean, optional
-        Center of gravity with periodic boundaries.
+        Center of gravity accounts for periodic boundaries.
 
     Returns
     -------

python/damask/tensor.py

@@ -8,7 +8,7 @@ All routines operate on numpy.ndarrays of shape (...,3,3).
 import numpy as _np
 
 
-def deviatoric(T):
+def deviatoric(T: _np.ndarray) -> _np.ndarray:
     """
     Calculate deviatoric part of a tensor.
 
@@ -26,7 +26,7 @@ def deviatoric(T):
     return T - spherical(T,tensor=True)
 
 
-def eigenvalues(T_sym):
+def eigenvalues(T_sym: _np.ndarray) -> _np.ndarray:
     """
     Eigenvalues, i.e. principal components, of a symmetric tensor.
 
@@ -45,7 +45,7 @@ def eigenvalues(T_sym):
     return _np.linalg.eigvalsh(symmetric(T_sym))
 
 
-def eigenvectors(T_sym,RHS=False):
+def eigenvectors(T_sym: _np.ndarray, RHS: bool = False) -> _np.ndarray:
     """
     Eigenvectors of a symmetric tensor.
 
@@ -70,7 +70,7 @@ def eigenvectors(T_sym,RHS=False):
     return v
 
 
-def spherical(T,tensor=True):
+def spherical(T: _np.ndarray, tensor: bool = True) -> _np.ndarray:
     """
     Calculate spherical part of a tensor.
 
@@ -92,7 +92,7 @@ def spherical(T,tensor=True):
     return _np.einsum('...jk,...',_np.eye(3),sph) if tensor else sph
 
 
-def symmetric(T):
+def symmetric(T: _np.ndarray) -> _np.ndarray:
     """
     Symmetrize tensor.
 
@@ -110,7 +110,7 @@ def symmetric(T):
     return (T+transpose(T))*0.5
 
 
-def transpose(T):
+def transpose(T: _np.ndarray) -> _np.ndarray:
     """
     Transpose tensor.
 

python/mypy.ini

@@ -0,0 +1,14 @@
+[mypy-scipy.*]
+ignore_missing_imports = True
+[mypy-h5py.*]
+ignore_missing_imports = True
+[mypy-vtk.*]
+ignore_missing_imports = True
+[mypy-PIL.*]
+ignore_missing_imports = True
+[mypy-matplotlib.*]
+ignore_missing_imports = True
+[mypy-pandas.*]
+ignore_missing_imports = True
+[mypy-wx.*]
+ignore_missing_imports = True

src/IO.f90

@@ -432,7 +432,7 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg)
       msg = 'Nconstituents mismatch between homogenization and material'
 
 !--------------------------------------------------------------------------------------------------
-! material error messages and related messages in mesh
+! material error messages and related messages in geometry
     case (150)
       msg = 'index out of bounds'
     case (153)
@@ -499,6 +499,11 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg)
     case (710)
       msg = 'Closing quotation mark missing in string'
 
+!-------------------------------------------------------------------------------------------------
+! errors related to the mesh solver
+    case (821)
+      msg = 'order not supported'
+
 !-------------------------------------------------------------------------------------------------
 ! errors related to the grid solver
     case (831)

src/mesh/FEM_utilities.f90

@@ -19,6 +19,7 @@ module FEM_utilities
   use IO
   use discretization_mesh
   use homogenization
+  use FEM_quadrature
 
   implicit none
   private
@@ -86,7 +87,9 @@ subroutine FEM_utilities_init
   class(tNode), pointer :: &
     num_mesh, &
     debug_mesh                                                                                      ! pointer to mesh debug options
-  integer :: structOrder                                                                            !< order of displacement shape functions
+  integer :: &
+    p_s, &                                                                                          !< order of shape functions
+    p_i                                                                                             !< integration order (quadrature rule)
   character(len=*), parameter :: &
     PETSCDEBUG = ' -snes_view -snes_monitor '
   PetscErrorCode            :: ierr
@@ -96,7 +99,14 @@ subroutine FEM_utilities_init
   print'(/,a)',   ' <<<+-  FEM_utilities init  -+>>>'
 
   num_mesh    => config_numerics%get('mesh',defaultVal=emptyDict)
-  structOrder =  num_mesh%get_asInt('structOrder', defaultVal = 2)
+
+  p_s = num_mesh%get_asInt('p_s',defaultVal = 2)
+  p_i = num_mesh%get_asInt('p_i',defaultVal = p_s)
+
+  if (p_s < 1_pInt .or. p_s > size(FEM_nQuadrature,2)) &
+    call IO_error(821,ext_msg='shape function order (p_s) out of bounds')
+  if (p_i < max(1_pInt,p_s-1_pInt) .or. p_i > p_s) &
+    call IO_error(821,ext_msg='integration order (p_i) out of bounds')
 
   debug_mesh  => config_debug%get('mesh',defaultVal=emptyList)
   debugPETSc  =  debug_mesh%contains('PETSc')
@@ -119,7 +129,7 @@ subroutine FEM_utilities_init
   CHKERRQ(ierr)
   call PetscOptionsInsertString(PETSC_NULL_OPTIONS,num_mesh%get_asString('PETSc_options',defaultVal=''),ierr)
   CHKERRQ(ierr)
-  write(petsc_optionsOrder,'(a,i0)') '-mechFE_petscspace_degree ', structOrder
+  write(petsc_optionsOrder,'(a,i0)') '-mechFE_petscspace_degree ', p_s
   call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_optionsOrder),ierr)
   CHKERRQ(ierr)
 

src/mesh/discretization_mesh.f90

@@ -85,7 +85,7 @@ subroutine discretization_mesh_init(restart)
     materialAt
   class(tNode), pointer :: &
     num_mesh
-  integer :: integrationOrder                                                                       !< order of quadrature rule required
+  integer :: p_i                                                                                    !< integration order (quadrature rule)
   type(tvec) :: coords_node0
 
   print'(/,a)',   ' <<<+-  discretization_mesh init  -+>>>'
@@ -93,7 +93,7 @@ subroutine discretization_mesh_init(restart)
 !--------------------------------------------------------------------------------
 ! read numerics parameter
   num_mesh => config_numerics%get('mesh',defaultVal=emptyDict)
-  integrationOrder = num_mesh%get_asInt('integrationorder',defaultVal = 2)
+  p_i = num_mesh%get_asInt('p_i',defaultVal = 2)
 
 !---------------------------------------------------------------------------------
 ! read debug parameters
@@ -150,9 +150,9 @@ subroutine discretization_mesh_init(restart)
   call VecGetArrayF90(coords_node0, mesh_node0_temp,ierr)
   CHKERRQ(ierr)
 
-  mesh_maxNips = FEM_nQuadrature(dimPlex,integrationOrder)
+  mesh_maxNips = FEM_nQuadrature(dimPlex,p_i)
 
-  call mesh_FEM_build_ipCoordinates(dimPlex,FEM_quadrature_points(dimPlex,integrationOrder)%p)
+  call mesh_FEM_build_ipCoordinates(dimPlex,FEM_quadrature_points(dimPlex,p_i)%p)
   call mesh_FEM_build_ipVolumes(dimPlex)
 
   allocate(materialAt(mesh_NcpElems))

src/mesh/mesh_mech_FEM.f90

@@ -41,7 +41,7 @@ module mesh_mechanical_FEM
 
   type, private :: tNumerics
     integer :: &
-      integrationOrder, &                                                                           !< order of quadrature rule required
+      p_i, &                                                                                        !< integration order (quadrature rule)
       itmax
     logical :: &
       BBarStabilisation
@@ -118,7 +118,7 @@ subroutine FEM_mechanical_init(fieldBC)
 !-----------------------------------------------------------------------------
 ! read numerical parametes and do sanity checks
   num_mesh => config_numerics%get('mesh',defaultVal=emptyDict)
-  num%integrationOrder  = num_mesh%get_asInt('integrationorder',defaultVal = 2)
+  num%p_i               = num_mesh%get_asInt('p_i',defaultVal = 2)
   num%itmax             = num_mesh%get_asInt('itmax',defaultVal=250)
   num%BBarStabilisation = num_mesh%get_asBool('bbarstabilisation',defaultVal = .false.)
   num%eps_struct_atol   = num_mesh%get_asFloat('eps_struct_atol', defaultVal = 1.0e-10_pReal)
@@ -135,9 +135,9 @@ subroutine FEM_mechanical_init(fieldBC)
 
 !--------------------------------------------------------------------------------------------------
 ! Setup FEM mech discretization
-  qPoints  = FEM_quadrature_points( dimPlex,num%integrationOrder)%p
+  qPoints  = FEM_quadrature_points( dimPlex,num%p_i)%p
-  qWeights = FEM_quadrature_weights(dimPlex,num%integrationOrder)%p
+  qWeights = FEM_quadrature_weights(dimPlex,num%p_i)%p
-  nQuadrature = FEM_nQuadrature(    dimPlex,num%integrationOrder)
+  nQuadrature = FEM_nQuadrature(    dimPlex,num%p_i)
   qPointsP  => qPoints
   qWeightsP => qWeights
   call PetscQuadratureCreate(PETSC_COMM_SELF,mechQuad,ierr); CHKERRQ(ierr)
@@ -146,7 +146,7 @@ subroutine FEM_mechanical_init(fieldBC)
   call PetscQuadratureSetData(mechQuad,dimPlex,nc,nQuadrature,qPointsP,qWeightsP,ierr)
   CHKERRQ(ierr)
   call PetscFECreateDefault(PETSC_COMM_SELF,dimPlex,nc,PETSC_TRUE,prefix, &
-                            num%integrationOrder,mechFE,ierr); CHKERRQ(ierr)
+                            num%p_i,mechFE,ierr); CHKERRQ(ierr)
   call PetscFESetQuadrature(mechFE,mechQuad,ierr); CHKERRQ(ierr)
   call PetscFEGetDimension(mechFE,nBasis,ierr); CHKERRQ(ierr)
   nBasis = nBasis/nc

src/phase_mechanical_plastic_dislotungsten.f90

@@ -24,7 +24,6 @@ submodule(phase:plastic) dislotungsten
       tau_Peierls, &                                                                                !< Peierls stress
       !* mobility law parameters
       Q_s, &                                                                                        !< activation energy for glide [J]
-      v_0, &                                                                                        !< dislocation velocity prefactor [m/s]
       p, &                                                                                          !< p-exponent in glide velocity
       q, &                                                                                          !< q-exponent in glide velocity
       B, &                                                                                          !< friction coefficient
@@ -158,7 +157,6 @@ module function plastic_dislotungsten_init() result(myPlasticity)
 
       rho_mob_0       = pl%get_as1dFloat('rho_mob_0',     requiredSize=size(N_sl))
       rho_dip_0       = pl%get_as1dFloat('rho_dip_0',     requiredSize=size(N_sl))
-      prm%v_0         = pl%get_as1dFloat('v_0',           requiredSize=size(N_sl))
       prm%b_sl        = pl%get_as1dFloat('b_sl',          requiredSize=size(N_sl))
       prm%Q_s         = pl%get_as1dFloat('Q_s',           requiredSize=size(N_sl))
 
@@ -189,18 +187,16 @@ module function plastic_dislotungsten_init() result(myPlasticity)
       prm%w              = math_expand(prm%w,           N_sl)
       prm%omega          = math_expand(prm%omega,       N_sl)
       prm%tau_Peierls    = math_expand(prm%tau_Peierls, N_sl)
-      prm%v_0            = math_expand(prm%v_0,         N_sl)
       prm%B              = math_expand(prm%B,           N_sl)
       prm%i_sl           = math_expand(prm%i_sl,        N_sl)
       prm%f_at           = math_expand(prm%f_at,        N_sl)
       prm%d_caron        = pl%get_asFloat('D_a') * prm%b_sl
 
       ! sanity checks
-      if (    prm%D_0          <= 0.0_pReal)  extmsg = trim(extmsg)//' D_0'
+      if (    prm%D_0          <  0.0_pReal)  extmsg = trim(extmsg)//' D_0'
       if (    prm%Q_cl         <= 0.0_pReal)  extmsg = trim(extmsg)//' Q_cl'
       if (any(rho_mob_0        <  0.0_pReal)) extmsg = trim(extmsg)//' rho_mob_0'
       if (any(rho_dip_0        <  0.0_pReal)) extmsg = trim(extmsg)//' rho_dip_0'
-      if (any(prm%v_0          <  0.0_pReal)) extmsg = trim(extmsg)//' v_0'
       if (any(prm%b_sl         <= 0.0_pReal)) extmsg = trim(extmsg)//' b_sl'
       if (any(prm%Q_s          <= 0.0_pReal)) extmsg = trim(extmsg)//' Q_s'
       if (any(prm%tau_Peierls  <  0.0_pReal)) extmsg = trim(extmsg)//' tau_Peierls'
@@ -211,7 +207,7 @@ module function plastic_dislotungsten_init() result(myPlasticity)
     else slipActive
       rho_mob_0 = emptyRealArray; rho_dip_0 = emptyRealArray
       allocate(prm%b_sl,prm%d_caron,prm%i_sl,prm%f_at,prm%tau_Peierls, &
-               prm%Q_s,prm%v_0,prm%p,prm%q,prm%B,prm%h,prm%w,prm%omega, &
+               prm%Q_s,prm%p,prm%q,prm%B,prm%h,prm%w,prm%omega, &
                source = emptyRealArray)
       allocate(prm%forestProjection(0,0))
       allocate(prm%h_sl_sl         (0,0))
@@ -287,6 +283,7 @@ pure module subroutine dislotungsten_LpAndItsTangent(Lp,dLp_dMp, &
     dot_gamma_pos,dot_gamma_neg, &
     ddot_gamma_dtau_pos,ddot_gamma_dtau_neg
 
+
   Lp = 0.0_pReal
   dLp_dMp = 0.0_pReal
 
@@ -328,6 +325,7 @@ module subroutine dislotungsten_dotState(Mp,T,ph,en)
     dot_rho_dip_climb, &
     d_hat
 
+
   associate(prm => param(ph), stt => state(ph), dot => dotState(ph), dst => dependentState(ph))
 
     call kinetics(Mp,T,ph,en,&
@@ -336,11 +334,11 @@ module subroutine dislotungsten_dotState(Mp,T,ph,en)
 
     dot%gamma_sl(:,en) = abs(dot_gamma_pos+dot_gamma_neg)
 
-  where(dEq0(tau_pos))                                                                              ! ToDo: use avg of +/-
+    where(dEq0((tau_pos+tau_neg)*0.5_pReal))
       dot_rho_dip_formation = 0.0_pReal
       dot_rho_dip_climb     = 0.0_pReal
     else where
-    d_hat = math_clip(3.0_pReal*prm%mu*prm%b_sl/(16.0_pReal*PI*abs(tau_pos)), &                     ! ToDo: use avg of +/-
+      d_hat = math_clip(3.0_pReal*prm%mu*prm%b_sl/(16.0_pReal*PI*abs(tau_pos+tau_neg)*0.5_pReal), &
                         prm%d_caron, &                                                              ! lower limit
                         dst%Lambda_sl(:,en))                                                        ! upper limit
       dot_rho_dip_formation = merge(2.0_pReal*(d_hat-prm%d_caron)*stt%rho_mob(:,en)*dot%gamma_sl(:,en)/prm%b_sl, &
@@ -373,16 +371,17 @@ module subroutine dislotungsten_dependentState(ph,en)
     en
 
   real(pReal), dimension(param(ph)%sum_N_sl) :: &
-    dislocationSpacing
+    Lambda_sl_inv
 
 
   associate(prm => param(ph), stt => state(ph), dst => dependentState(ph))
 
-    dislocationSpacing = sqrt(matmul(prm%forestProjection,stt%rho_mob(:,en)+stt%rho_dip(:,en)))
     dst%tau_pass(:,en) = prm%mu*prm%b_sl &
                        * sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,en)+stt%rho_dip(:,en)))
 
-    dst%Lambda_sl(:,en) = prm%D/(1.0_pReal+prm%D*dislocationSpacing/prm%i_sl)
+    Lambda_sl_inv = 1.0_pReal/prm%D &
+                  + sqrt(matmul(prm%forestProjection,stt%rho_mob(:,en)+stt%rho_dip(:,en)))/prm%i_sl
+    dst%Lambda_sl(:,en) = Lambda_sl_inv**(-1.0_pReal)
 
   end associate
 
@@ -456,14 +455,16 @@ pure subroutine kinetics(Mp,T,ph,en, &
     ddot_gamma_dtau_neg, &
     tau_pos_out, &
     tau_neg_out
+
   real(pReal), dimension(param(ph)%sum_N_sl) :: &
     StressRatio, &
     StressRatio_p,StressRatio_pminus1, &
-    dvel, vel, &
+    dvel, &
-    tau_pos,tau_neg, &
+    tau_pos, tau_neg, tau_eff, &
     t_n, t_k, dtk,dtn
   integer :: j
 
+
   associate(prm => param(ph), stt => state(ph), dst => dependentState(ph))
 
     do j = 1, prm%sum_N_sl
@@ -471,25 +472,25 @@ pure subroutine kinetics(Mp,T,ph,en, &
       tau_neg(j) = math_tensordot(Mp,prm%P_nS_neg(1:3,1:3,j))
     end do
 
-
     if (present(tau_pos_out)) tau_pos_out = tau_pos
     if (present(tau_neg_out)) tau_neg_out = tau_neg
 
     associate(BoltzmannRatio  => prm%Q_s/(kB*T), &
-            dot_gamma_0     => stt%rho_mob(:,en)*prm%b_sl*prm%v_0, &
+              b_rho_half      => stt%rho_mob(:,en) * prm%b_sl * 0.5_pReal, &
               effectiveLength => dst%Lambda_sl(:,en) - prm%w)
 
-  significantPositiveTau: where(abs(tau_pos)-dst%tau_pass(:,en) > tol_math_check)
+      tau_eff = abs(tau_pos)-dst%tau_pass(:,en)
-    StressRatio = (abs(tau_pos)-dst%tau_pass(:,en))/prm%tau_Peierls
+
+      significantPositiveTau: where(tau_eff > tol_math_check)
+        StressRatio = tau_eff/prm%tau_Peierls
         StressRatio_p       = StressRatio** prm%p
         StressRatio_pminus1 = StressRatio**(prm%p-1.0_pReal)
 
-    t_n = prm%b_sl/(exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q)*prm%omega*effectiveLength)
+        t_n = prm%b_sl*exp(BoltzmannRatio*(1.0_pReal-StressRatio_p) ** prm%q) &
-    t_k = effectiveLength * prm%B /(2.0_pReal*prm%b_sl*tau_pos)
+            / (prm%omega*effectiveLength)
+        t_k = effectiveLength * prm%B /(2.0_pReal*prm%b_sl*tau_eff)                                 ! corrected eq. (14)
 
-    vel = prm%h/(t_n + t_k)
+        dot_gamma_pos = b_rho_half * sign(prm%h/(t_n + t_k),tau_pos)
-
-    dot_gamma_pos = dot_gamma_0 * sign(vel,tau_pos) * 0.5_pReal
       else where significantPositiveTau
         dot_gamma_pos = 0.0_pReal
       end where significantPositiveTau
@@ -497,28 +498,29 @@ pure subroutine kinetics(Mp,T,ph,en, &
       if (present(ddot_gamma_dtau_pos)) then
         significantPositiveTau2: where(abs(tau_pos)-dst%tau_pass(:,en) > tol_math_check)
           dtn = -1.0_pReal * t_n * BoltzmannRatio * prm%p * prm%q * (1.0_pReal-StressRatio_p)**(prm%q - 1.0_pReal) &
-          * (StressRatio)**(prm%p - 1.0_pReal) / prm%tau_Peierls
+              * StressRatio_pminus1 / prm%tau_Peierls
           dtk = -1.0_pReal * t_k / tau_pos
 
           dvel = -1.0_pReal * prm%h * (dtk + dtn) / (t_n + t_k)**2.0_pReal
 
-      ddot_gamma_dtau_pos = dot_gamma_0 * dvel* 0.5_pReal
+          ddot_gamma_dtau_pos = b_rho_half * dvel
         else where significantPositiveTau2
           ddot_gamma_dtau_pos = 0.0_pReal
         end where significantPositiveTau2
       end if
 
-  significantNegativeTau: where(abs(tau_neg)-dst%tau_pass(:,en) > tol_math_check)
+      tau_eff = abs(tau_neg)-dst%tau_pass(:,en)
-    StressRatio = (abs(tau_neg)-dst%tau_pass(:,en))/prm%tau_Peierls
+
+      significantNegativeTau: where(tau_eff > tol_math_check)
+        StressRatio = tau_eff/prm%tau_Peierls
         StressRatio_p       = StressRatio** prm%p
         StressRatio_pminus1 = StressRatio**(prm%p-1.0_pReal)
 
-    t_n = prm%b_sl/(exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q)*prm%omega*effectiveLength)
+        t_n = prm%b_sl*exp(BoltzmannRatio*(1.0_pReal-StressRatio_p) ** prm%q) &
-    t_k = effectiveLength * prm%B /(2.0_pReal*prm%b_sl*tau_pos)
+            / (prm%omega*effectiveLength)
+        t_k = effectiveLength * prm%B /(2.0_pReal*prm%b_sl*tau_eff)                                 ! corrected eq. (14)
 
-    vel = prm%h/(t_n + t_k)
+        dot_gamma_neg = b_rho_half * sign(prm%h/(t_n + t_k),tau_neg)
-
-    dot_gamma_neg = dot_gamma_0 * sign(vel,tau_neg) * 0.5_pReal
       else where significantNegativeTau
         dot_gamma_neg = 0.0_pReal
       end where significantNegativeTau
@@ -526,12 +528,12 @@ pure subroutine kinetics(Mp,T,ph,en, &
       if (present(ddot_gamma_dtau_neg)) then
         significantNegativeTau2: where(abs(tau_neg)-dst%tau_pass(:,en) > tol_math_check)
           dtn = -1.0_pReal * t_n * BoltzmannRatio * prm%p * prm%q * (1.0_pReal-StressRatio_p)**(prm%q - 1.0_pReal) &
-          * (StressRatio)**(prm%p - 1.0_pReal) / prm%tau_Peierls
+              * StressRatio_pminus1 / prm%tau_Peierls
           dtk = -1.0_pReal * t_k / tau_neg
 
           dvel = -1.0_pReal * prm%h * (dtk + dtn) / (t_n + t_k)**2.0_pReal
 
-      ddot_gamma_dtau_neg = dot_gamma_0 * dvel * 0.5_pReal
+          ddot_gamma_dtau_neg = b_rho_half * dvel
         else where significantNegativeTau2
           ddot_gamma_dtau_neg = 0.0_pReal
         end where significantNegativeTau2

src/phase_mechanical_plastic_dislotwin.f90

@@ -890,7 +890,8 @@ pure subroutine kinetics_sl(Mp,T,ph,en, &
       stressRatio    = tau_eff/prm%tau_0
       StressRatio_p  = stressRatio** prm%p
       Q_kB_T = prm%Q_sl/(kB*T)
-      v_wait_inverse = prm%v_0**(-1.0_pReal) * exp(Q_kB_T*(1.0_pReal-StressRatio_p)** prm%q)
+      v_wait_inverse = exp(Q_kB_T*(1.0_pReal-StressRatio_p)** prm%q) &
+                     / prm%v_0
       v_run_inverse  = prm%B/(tau_eff*prm%b_sl)
 
       dot_gamma_sl = sign(stt%rho_mob(:,en)*prm%b_sl/(v_wait_inverse+v_run_inverse),tau)