Merge branch 'development' into MiscImprovements

.gitlab-ci.yml

@@ -213,6 +213,14 @@ Post_OrientationConversion:
     - master
     - release
 
+Post_OrientationAverageMisorientation:
+  stage: postprocessing
+  script:
+    - OrientationAverageMisorientation/test.py
+  except:
+    - master
+    - release
+
 ###################################################################################################
 grid_mech_compile_Intel:
   stage: compilePETSc

VERSION

@@ -1 +1 @@
-v2.0.3-198-g1c762860
+v2.0.3-229-g90fa809f

python/damask/quaternion.py

@@ -0,0 +1,214 @@
+# -*- coding: UTF-8 no BOM -*-
+
+import numpy as np
+
+P = -1                                                                                              # convention (see DOI:10.1088/0965-0393/23/8/083501)
+
+####################################################################################################
+class Quaternion:
+    u"""
+    Quaternion with basic operations
+    
+    q is the real part, p = (x, y, z) are the imaginary parts.
+    Definition of multiplication depends on variable P, P ∈ {-1,1}.
+    """                                                                              
+
+    def __init__(self,
+                 q    = 0.0,
+                 p    = np.zeros(3,dtype=float)):
+      """Initializes to identity unless specified"""
+      self.q = float(q)
+      self.p = np.array(p,dtype=float)
+
+
+    def __copy__(self):
+      """Copy"""
+      return self.__class__(q=self.q,
+                            p=self.p.copy())
+
+    copy = __copy__
+
+
+    def __iter__(self):
+      """Components"""
+      return iter(self.asList())
+
+    def __repr__(self):
+      """Readable string"""
+      return 'Quaternion: (real={q:+.6f}, imag=<{p[0]:+.6f}, {p[1]:+.6f}, {p[2]:+.6f}>)'.format(q=self.q,p=self.p)
+      
+
+    def __add__(self, other):
+      """Addition"""
+      if isinstance(other, Quaternion):
+        return self.__class__(q=self.q + other.q,
+                              p=self.p + other.p)
+      else:
+        return NotImplemented
+
+    def __iadd__(self, other):
+      """In-place addition"""
+      if isinstance(other, Quaternion):
+        self.q += other.q
+        self.p += other.p
+        return self
+      else:
+        return NotImplemented
+      
+    def __pos__(self):                                                                      
+      """Unary positive operator"""                                                                       
+      return self
+      
+      
+    def __sub__(self, other):
+      """Subtraction"""
+      if isinstance(other, Quaternion):
+        return self.__class__(q=self.q - other.q,
+                              p=self.p - other.p)
+      else:
+        return NotImplemented
+
+    def __isub__(self, other):
+      """In-place subtraction"""
+      if isinstance(other, Quaternion):
+        self.q -= other.q
+        self.p -= other.p
+        return self
+      else:
+        return NotImplemented
+        
+    def __neg__(self):                                                                      
+      """Unary negative operator"""  
+      return self * -1.0
+      
+      
+    def __mul__(self, other):
+      """Multiplication with quaternion or scalar"""
+      if isinstance(other, Quaternion):
+        return self.__class__(q=self.q*other.q - np.dot(self.p,other.p),
+                              p=self.q*other.p + other.q*self.p + P * np.cross(self.p,other.p))
+      elif isinstance(other, (int, float)):
+        return self.__class__(q=self.q*other,
+                              p=self.p*other)
+      else:
+        return NotImplemented
+
+    def __imul__(self, other):
+      """In-place multiplication with quaternion or scalar"""
+      if isinstance(other, Quaternion):
+        self.q = self.q*other.q - np.dot(self.p,other.p)
+        self.p = self.q*other.p + other.q*self.p + P * np.cross(self.p,other.p)
+        return self
+      elif isinstance(other, (int, float)):
+        self.q *= other
+        self.p *= other
+        return self
+      else:
+        return NotImplemented
+
+
+    def __truediv__(self, other):
+      """Divsion with quaternion or scalar"""
+      if isinstance(other, Quaternion):
+        s = other.conjugate()/abs(other)**2.
+        return self.__class__(q=self.q * s,
+                              p=self.p * s)
+      elif isinstance(other, (int, float)):
+        return self.__class__(q=self.q / other,
+                              p=self.p / other)
+      else:
+        return NotImplemented
+
+    def __itruediv__(self, other):
+      """In-place divsion with quaternion or scalar"""
+      if isinstance(other, Quaternion):
+        s = other.conjugate()/abs(other)**2.
+        self *= s
+        return self
+      elif isinstance(other, (int, float)):
+        self.q /= other
+        self.p /= other
+        return self
+      else:
+        return NotImplemented
+
+
+    def __pow__(self, exponent):
+      """Power"""
+      if isinstance(exponent, (int, float)):
+        omega = np.acos(self.q)
+        return self.__class__(q=         np.cos(exponent*omega),
+                              p=self.p * np.sin(exponent*omega)/np.sin(omega))
+      else:
+        return NotImplemented
+
+    def __ipow__(self, exponent):
+      """In-place power"""
+      if isinstance(exponent, (int, float)):
+        omega = np.acos(self.q)
+        self.q  = np.cos(exponent*omega)
+        self.p *= np.sin(exponent*omega)/np.sin(omega)
+        return self
+      else:
+        return NotImplemented
+      
+
+    def __abs__(self):
+      """Norm"""
+      return np.sqrt(self.q ** 2 + np.dot(self.p,self.p))
+
+    magnitude = __abs__
+
+
+    def __eq__(self,other):
+      """Equal (sufficiently close) to each other"""
+      return np.isclose(( self-other).magnitude(),0.0) \
+          or np.isclose((-self-other).magnitude(),0.0)
+
+    def __ne__(self,other):
+      """Not equal (sufficiently close) to each other"""
+      return not self.__eq__(other)
+
+
+    def asM(self):
+      """Intermediate representation useful for quaternion averaging (see F. Landis Markley et al.)"""
+      return np.outer(self.asArray(),self.asArray())
+
+    def asArray(self):
+      """As numpy array"""
+      return np.array((self.q,self.p[0],self.p[1],self.p[2]))
+ 
+    def asList(self):
+      """As list"""
+      return [self.q]+list(self.p)
+
+
+    def normalize(self):
+      """Normalizes in-place"""
+      d = self.magnitude()
+      if d > 0.0: self /= d
+      return self
+      
+    def normalized(self):
+      """Returns normalized copy"""
+      return self.copy().normalize()
+
+
+    def conjugate(self):
+      """Conjugates in-place"""
+      self.p = -self.p
+      return self
+      
+    def conjugated(self):
+      """Returns conjugated copy"""
+      return self.copy().conjugate()
+
+
+    def homomorph(self):
+      """Homomorphs in-place"""
+      self *= -1.0
+      return self
+      
+    def homomorphed(self):
+      """Returns homomorphed copy"""
+      return self.copy().homomorph()

src/plastic_disloUCLA.f90

@@ -121,15 +121,7 @@ subroutine plastic_disloUCLA_init()
     math_expand
   use IO, only: &
     IO_error
-  use material, only: &
-    phase_plasticity, &
-    phase_plasticityInstance, &
-    phase_Noutput, &
-    material_allocatePlasticState, &
-    PLASTICITY_DISLOUCLA_label, &
-    PLASTICITY_DISLOUCLA_ID, &
-    material_phase, &
-    plasticState
+  use material
   use config, only: &
     config_phase
   use lattice

src/plastic_dislotwin.f90

@@ -185,15 +185,7 @@ subroutine plastic_dislotwin_init
    PI
  use IO, only: &
    IO_error
- use material, only: &
-   phase_plasticity, &
-   phase_plasticityInstance, &
-   phase_Noutput, &
-   material_allocatePlasticState, &
-   PLASTICITY_DISLOTWIN_label, &
-   PLASTICITY_DISLOTWIN_ID, &
-   material_phase, &  
-   plasticState
+ use material
  use config, only: &
    config_phase
  use lattice

src/plastic_isotropic.f90

@@ -91,19 +91,7 @@ subroutine plastic_isotropic_init
     debug_levelBasic
   use IO, only: &
     IO_error
-  use material, only: &
-    phase_plasticity, &
-    phase_plasticityInstance, &
-    phase_Noutput, &
-    material_allocatePlasticState, &
-    PLASTICITY_ISOTROPIC_label, &
-    PLASTICITY_ISOTROPIC_ID, &
-    material_phase, &
-    plasticState
-#ifdef DEBUG
-  use material, only: &
-    phasememberAt
-#endif
+  use material
   use config, only: &
     config_phase
   use lattice
@@ -464,7 +452,7 @@ function plastic_isotropic_postResults(Mp,instance,of) result(postResults)
         c = c + 1
       case (dot_gamma_ID)
         postResults(c+1) = prm%dot_gamma_0 &
-                         * (sqrt(1.5_pReal) * norm_Mp /(prm%M * stt%xi(of)))**prm%n
+                              * (sqrt(1.5_pReal) * norm_Mp /(prm%M * stt%xi(of)))**prm%n
         c = c + 1
  
     end select
@@ -479,11 +467,10 @@ end function plastic_isotropic_postResults
 !> @brief writes results to HDF5 output file
 !--------------------------------------------------------------------------------------------------
 subroutine plastic_isotropic_results(instance,group)
-#if defined(PETSc) || defined(DAMASK_HDF5)
-  use results, only: &
-    results_writeDataset
+#if defined(PETSc) || defined(DAMASKHDF5)
+  use results
 
-  integer,          intent(in) :: instance
+  integer, intent(in) :: instance
   character(len=*), intent(in) :: group
   
   integer :: o
@@ -496,10 +483,9 @@ subroutine plastic_isotropic_results(instance,group)
     end select
   enddo outputsLoop
   end associate
-  
 #else
-  integer,          intent(in) :: instance
-  character(len=*), intent(in) :: group
+  integer, intent(in) :: instance
+  character(len=*) :: group
 #endif
 
 end subroutine plastic_isotropic_results

src/plastic_kinematichardening.f90

@@ -112,19 +112,7 @@ subroutine plastic_kinehardening_init
    math_expand
  use IO, only: &
    IO_error
- use material, only: &
-   phase_plasticity, &
-   phase_plasticityInstance, &
-   phase_Noutput, &
-   material_allocatePlasticState, &
-   PLASTICITY_kinehardening_label, &
-   PLASTICITY_kinehardening_ID, &
-   material_phase, &
-   plasticState
-#ifdef DEBUG
- use material, only: &
-   phasememberAt
-#endif
+ use material
  use config, only: &
    config_phase
  use lattice

src/plastic_none.f90

@@ -23,13 +23,7 @@ subroutine plastic_none_init
    debug_level, &
    debug_constitutive, &
    debug_levelBasic
- use material, only: &
-   phase_plasticity, &
-   material_allocatePlasticState, &
-   PLASTICITY_NONE_label, &
-   PLASTICITY_NONE_ID, &
-   material_phase, &
-   plasticState
+ use material
 
  integer :: &
    Ninstance, &

src/plastic_nonlocal.f90

@@ -248,15 +248,7 @@ subroutine plastic_nonlocal_init
     debug_levelBasic
   use mesh, only: &
     theMesh
-  use material, only: &
-    phase_plasticity, &
-    phase_plasticityInstance, &
-    phase_Noutput, &
-    PLASTICITY_NONLOCAL_label, &
-    PLASTICITY_NONLOCAL_ID, &
-    plasticState, &
-    material_phase, &
-    material_allocatePlasticState
+  use material
   use config
   use lattice
 

src/plastic_phenopowerlaw.f90

@@ -116,15 +116,7 @@ subroutine plastic_phenopowerlaw_init
    math_expand
  use IO, only: &
    IO_error
- use material, only: &
-   phase_plasticity, &
-   phase_plasticityInstance, &
-   phase_Noutput, &
-   material_allocatePlasticState, &
-   PLASTICITY_PHENOPOWERLAW_LABEL, &
-   PLASTICITY_PHENOPOWERLAW_ID, &
-   material_phase, &
-   plasticState
+ use material
  use config, only: &
    config_phase
  use lattice