From 1808b37357ed53180b1092dbd1a2e11e1d19b354 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 25 Nov 2018 11:17:12 +0100 Subject: [PATCH 01/70] can be easily computed during postprocessing --- src/plastic_phenopowerlaw.f90 | 17 +---------------- 1 file changed, 1 insertion(+), 16 deletions(-) diff --git a/src/plastic_phenopowerlaw.f90 b/src/plastic_phenopowerlaw.f90 index 57d48d109..0e85a65c9 100644 --- a/src/plastic_phenopowerlaw.f90 +++ b/src/plastic_phenopowerlaw.f90 @@ -24,12 +24,10 @@ module plastic_phenopowerlaw accumulatedshear_slip_ID, & shearrate_slip_ID, & resolvedstress_slip_ID, & - totalshear_ID, & resistance_twin_ID, & accumulatedshear_twin_ID, & shearrate_twin_ID, & - resolvedstress_twin_ID, & - totalvolfrac_twin_ID + resolvedstress_twin_ID end enum type, private :: tParameters @@ -338,12 +336,6 @@ subroutine plastic_phenopowerlaw_init outputID = merge(resolvedstress_twin_ID,undefined_ID,prm%totalNtwin>0_pInt) outputSize = prm%totalNtwin - case ('totalshear') - outputID = merge(totalshear_ID,undefined_ID,prm%totalNslip>0_pInt) - outputSize = 1_pInt - case ('totalvolfrac_twin') - outputID = merge(totalvolfrac_twin_ID,undefined_ID,prm%totalNtwin>0_pInt) - outputSize = 1_pInt end select if (outputID /= undefined_ID) then @@ -732,13 +724,6 @@ function plastic_phenopowerlaw_postResults(Mp,instance,of) result(postResults) enddo c = c + prm%totalNtwin - case (totalshear_ID) - postResults(c+1_pInt) = stt%sumGamma(of) - c = c + 1_pInt - case (totalvolfrac_twin_ID) - postResults(c+1_pInt) = stt%sumF(of) - c = c + 1_pInt - end select enddo outputsLoop end associate From 8a253856f1262b10bad16bd71b74c51f25c23407 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 25 Nov 2018 15:38:14 +0100 Subject: [PATCH 02/70] no need to perform the state integration --- src/plastic_phenopowerlaw.f90 | 15 ++++----------- 1 file changed, 4 insertions(+), 11 deletions(-) diff --git a/src/plastic_phenopowerlaw.f90 b/src/plastic_phenopowerlaw.f90 index 0e85a65c9..c777bff8c 100644 --- a/src/plastic_phenopowerlaw.f90 +++ b/src/plastic_phenopowerlaw.f90 @@ -79,7 +79,6 @@ module plastic_phenopowerlaw type, private :: tPhenopowerlawState real(pReal), pointer, dimension(:) :: & - sumGamma, & ! ToDo: why not make a dependent state? sumF ! ToDo: why not make a dependent state? real(pReal), pointer, dimension(:,:) :: & xi_slip, & @@ -351,7 +350,7 @@ subroutine plastic_phenopowerlaw_init NipcMyPhase = count(material_phase == p) ! number of IPCs containing my phase sizeState = size(['tau_slip ','gamma_slip']) * prm%TotalNslip & + size(['tau_twin ','gamma_twin']) * prm%TotalNtwin & - + size(['sum(gamma)','sum(f) ']) ! ToDo: only needed if either twin or slip active! + + size(['sum(f) ']) ! ToDo: only needed if either twin or slip active! sizeDotState = sizeState call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0_pInt, & @@ -375,12 +374,6 @@ subroutine plastic_phenopowerlaw_init dot%xi_twin => plasticState(p)%dotState(startIndex:endIndex,:) plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolResistance - startIndex = endIndex + 1_pInt - endIndex = endIndex + 1_pInt - stt%sumGamma => plasticState(p)%state (startIndex,:) - dot%sumGamma => plasticState(p)%dotState(startIndex,:) - plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolShear - startIndex = endIndex + 1_pInt endIndex = endIndex + 1_pInt stt%sumF=>plasticState(p)%state (startIndex,:) @@ -482,7 +475,7 @@ subroutine plastic_phenopowerlaw_dotState(Mp,instance,of) i real(pReal) :: & c_SlipSlip,c_TwinSlip,c_TwinTwin, & - xi_slip_sat_offset + xi_slip_sat_offset,sumGamma real(pReal), dimension(param(instance)%totalNslip) :: & left_SlipSlip,right_SlipSlip, & @@ -494,11 +487,12 @@ subroutine plastic_phenopowerlaw_dotState(Mp,instance,of) associate(prm => param(instance), stt => state(instance), dot => dotState(instance)) dot%whole(:,of) = 0.0_pReal + sumGamma = sum(stt%gamma_slip(:,of)) !-------------------------------------------------------------------------------------------------- ! system-independent (nonlinear) prefactors to M_Xx (X influenced by x) matrices c_SlipSlip = prm%h0_slipslip * (1.0_pReal + prm%twinC*stt%sumF(of)** prm%twinB) - c_TwinSlip = prm%h0_TwinSlip * stt%sumGamma(of)**prm%twinE + c_TwinSlip = prm%h0_TwinSlip * sumGamma**prm%twinE c_TwinTwin = prm%h0_TwinTwin * stt%sumF(of)**prm%twinD !-------------------------------------------------------------------------------------------------- @@ -512,7 +506,6 @@ subroutine plastic_phenopowerlaw_dotState(Mp,instance,of) ! shear rates call kinetics_slip(prm,stt,of,Mp,gdot_slip_pos,gdot_slip_neg) dot%gamma_slip(:,of) = abs(gdot_slip_pos+gdot_slip_neg) - dot%sumGamma(of) = sum(dot%gamma_slip(:,of)) call kinetics_twin(prm,stt,of,Mp,dot%gamma_twin(:,of)) if (prm%totalNtwin > 0_pInt) dot%sumF(of) = merge(sum(dot%gamma_twin(:,of)/prm%gamma_twin_char), & 0.0_pReal, & From 7cc2892e646eb25dc61953d3e85e7ec31ad9fab0 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 25 Nov 2018 16:14:46 +0100 Subject: [PATCH 03/70] no need to calculate twinned volume fraction as state --- src/plastic_phenopowerlaw.f90 | 33 ++++++++++++--------------------- 1 file changed, 12 insertions(+), 21 deletions(-) diff --git a/src/plastic_phenopowerlaw.f90 b/src/plastic_phenopowerlaw.f90 index c777bff8c..c2f5edb68 100644 --- a/src/plastic_phenopowerlaw.f90 +++ b/src/plastic_phenopowerlaw.f90 @@ -78,8 +78,6 @@ module plastic_phenopowerlaw type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance) type, private :: tPhenopowerlawState - real(pReal), pointer, dimension(:) :: & - sumF ! ToDo: why not make a dependent state? real(pReal), pointer, dimension(:,:) :: & xi_slip, & xi_twin, & @@ -349,8 +347,7 @@ subroutine plastic_phenopowerlaw_init ! allocate state arrays NipcMyPhase = count(material_phase == p) ! number of IPCs containing my phase sizeState = size(['tau_slip ','gamma_slip']) * prm%TotalNslip & - + size(['tau_twin ','gamma_twin']) * prm%TotalNtwin & - + size(['sum(f) ']) ! ToDo: only needed if either twin or slip active! + + size(['tau_twin ','gamma_twin']) * prm%TotalNtwin sizeDotState = sizeState call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0_pInt, & @@ -374,12 +371,6 @@ subroutine plastic_phenopowerlaw_init dot%xi_twin => plasticState(p)%dotState(startIndex:endIndex,:) plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolResistance - startIndex = endIndex + 1_pInt - endIndex = endIndex + 1_pInt - stt%sumF=>plasticState(p)%state (startIndex,:) - dot%sumF=>plasticState(p)%dotState(startIndex,:) - plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolTwinFrac - startIndex = endIndex + 1_pInt endIndex = endIndex + prm%totalNslip stt%gamma_slip => plasticState(p)%state (startIndex:endIndex,:) @@ -439,7 +430,8 @@ subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) call kinetics_slip(prm,stt,of,Mp,gdot_slip_pos,gdot_slip_neg,dgdot_dtauslip_pos,dgdot_dtauslip_neg) slipSystems: do i = 1_pInt, prm%totalNslip - Lp = Lp + (1.0_pReal-stt%sumF(of))*(gdot_slip_pos(i)+gdot_slip_neg(i))*prm%Schmid_slip(1:3,1:3,i) + Lp = Lp + (1.0_pReal- sum(stt%gamma_twin(:,of)/prm%gamma_twin_char)) & ! only shear in untwinned volume + * (gdot_slip_pos(i)+gdot_slip_neg(i))*prm%Schmid_slip(1:3,1:3,i) forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & + dgdot_dtauslip_pos(i) * prm%Schmid_slip(k,l,i) * prm%nonSchmid_pos(m,n,i) & @@ -475,8 +467,8 @@ subroutine plastic_phenopowerlaw_dotState(Mp,instance,of) i real(pReal) :: & c_SlipSlip,c_TwinSlip,c_TwinTwin, & - xi_slip_sat_offset,sumGamma - + xi_slip_sat_offset,& + sumGamma,sumF real(pReal), dimension(param(instance)%totalNslip) :: & left_SlipSlip,right_SlipSlip, & gdot_slip_pos,gdot_slip_neg @@ -488,17 +480,18 @@ subroutine plastic_phenopowerlaw_dotState(Mp,instance,of) dot%whole(:,of) = 0.0_pReal sumGamma = sum(stt%gamma_slip(:,of)) + sumF = sum(stt%gamma_twin(:,of)/prm%gamma_twin_char) !-------------------------------------------------------------------------------------------------- ! system-independent (nonlinear) prefactors to M_Xx (X influenced by x) matrices - c_SlipSlip = prm%h0_slipslip * (1.0_pReal + prm%twinC*stt%sumF(of)** prm%twinB) + c_SlipSlip = prm%h0_slipslip * (1.0_pReal + prm%twinC*sumF** prm%twinB) c_TwinSlip = prm%h0_TwinSlip * sumGamma**prm%twinE - c_TwinTwin = prm%h0_TwinTwin * stt%sumF(of)**prm%twinD + c_TwinTwin = prm%h0_TwinTwin * sumF**prm%twinD !-------------------------------------------------------------------------------------------------- ! calculate left and right vectors left_SlipSlip = 1.0_pReal + prm%H_int - xi_slip_sat_offset = prm%spr*sqrt(stt%sumF(of)) + xi_slip_sat_offset = prm%spr*sqrt(sumF) right_SlipSlip = abs(1.0_pReal-stt%xi_slip(:,of) / (prm%xi_slip_sat+xi_slip_sat_offset)) **prm%a_slip & * sign(1.0_pReal,1.0_pReal-stt%xi_slip(:,of) / (prm%xi_slip_sat+xi_slip_sat_offset)) @@ -507,9 +500,6 @@ subroutine plastic_phenopowerlaw_dotState(Mp,instance,of) call kinetics_slip(prm,stt,of,Mp,gdot_slip_pos,gdot_slip_neg) dot%gamma_slip(:,of) = abs(gdot_slip_pos+gdot_slip_neg) call kinetics_twin(prm,stt,of,Mp,dot%gamma_twin(:,of)) - if (prm%totalNtwin > 0_pInt) dot%sumF(of) = merge(sum(dot%gamma_twin(:,of)/prm%gamma_twin_char), & - 0.0_pReal, & - stt%sumF(of) < 0.98_pReal) !-------------------------------------------------------------------------------------------------- ! hardening @@ -605,7 +595,7 @@ end subroutine kinetics_slip !-------------------------------------------------------------------------------------------------- !> @brief calculates shear rates on twin systems and derivatives with respect to resolved stress -!> @details: Shear rates are calculated only optionally. NOTE: Against the common convention, the +!> @details: Derivates are calculated only optionally. NOTE: Against the common convention, the !> result (i.e. intent(out)) variables are the last to have the optional arguments at the end !-------------------------------------------------------------------------------------------------- pure subroutine kinetics_twin(prm,stt,of,Mp,gdot_twin,dgdot_dtau_twin) @@ -637,7 +627,8 @@ pure subroutine kinetics_twin(prm,stt,of,Mp,gdot_twin,dgdot_dtau_twin) enddo where(tau_twin > 0.0_pReal) - gdot_twin = (1.0_pReal-stt%sumF(of))*prm%gdot0_twin*(abs(tau_twin)/stt%xi_twin(:,of))**prm%n_twin + gdot_twin = (1.0_pReal-sum(stt%gamma_twin(:,of)/prm%gamma_twin_char)) & ! only twin in untwinned volume fraction + * prm%gdot0_twin*(abs(tau_twin)/stt%xi_twin(:,of))**prm%n_twin else where gdot_twin = 0.0_pReal end where From 638789b111d1540b6133b8c4449f9d4878d368c7 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 25 Nov 2018 16:16:46 +0100 Subject: [PATCH 04/70] totalvolfrac and totalshear don't exist anymore use postprocessing tools if needed --- examples/ConfigFiles/Phase_Phenopowerlaw_Aluminum.config | 1 - examples/ConfigFiles/Phase_Phenopowerlaw_BCC-Ferrite.config | 1 - examples/ConfigFiles/Phase_Phenopowerlaw_BCC-Martensite.config | 1 - examples/ConfigFiles/Phase_Phenopowerlaw_Gold.config | 2 -- examples/ConfigFiles/Phase_Phenopowerlaw_Magnesium.config | 2 -- examples/ConfigFiles/Phase_Phenopowerlaw_cpTi-alpha.config | 2 -- examples/ConfigFiles/Phase_Phenopowerlaw_multiField.config | 2 -- 7 files changed, 11 deletions(-) diff --git a/examples/ConfigFiles/Phase_Phenopowerlaw_Aluminum.config b/examples/ConfigFiles/Phase_Phenopowerlaw_Aluminum.config index 63ec8f3c8..938961d8e 100644 --- a/examples/ConfigFiles/Phase_Phenopowerlaw_Aluminum.config +++ b/examples/ConfigFiles/Phase_Phenopowerlaw_Aluminum.config @@ -6,7 +6,6 @@ plasticity phenopowerlaw (output) shearrate_slip (output) resolvedstress_slip (output) accumulated_shear_slip -(output) totalshear lattice_structure fcc Nslip 12 # per family diff --git a/examples/ConfigFiles/Phase_Phenopowerlaw_BCC-Ferrite.config b/examples/ConfigFiles/Phase_Phenopowerlaw_BCC-Ferrite.config index 594c5dc22..87b4e2312 100644 --- a/examples/ConfigFiles/Phase_Phenopowerlaw_BCC-Ferrite.config +++ b/examples/ConfigFiles/Phase_Phenopowerlaw_BCC-Ferrite.config @@ -19,4 +19,3 @@ tausat_slip 222.e6 412.7e6 # per family, optimization long h0_slipslip 1000.0e6 interaction_slipslip 1 1 1.4 1.4 1.4 1.4 w0_slip 2.0 -(output) totalshear diff --git a/examples/ConfigFiles/Phase_Phenopowerlaw_BCC-Martensite.config b/examples/ConfigFiles/Phase_Phenopowerlaw_BCC-Martensite.config index c86d516a9..a04f27e7f 100644 --- a/examples/ConfigFiles/Phase_Phenopowerlaw_BCC-Martensite.config +++ b/examples/ConfigFiles/Phase_Phenopowerlaw_BCC-Martensite.config @@ -19,4 +19,3 @@ tausat_slip 872.9e6 971.2e6 # per family h0_slipslip 563.0e9 interaction_slipslip 1 1 1.4 1.4 1.4 1.4 a_slip 2.0 -(output) totalshear diff --git a/examples/ConfigFiles/Phase_Phenopowerlaw_Gold.config b/examples/ConfigFiles/Phase_Phenopowerlaw_Gold.config index a2e06fc07..bb6f8fc83 100644 --- a/examples/ConfigFiles/Phase_Phenopowerlaw_Gold.config +++ b/examples/ConfigFiles/Phase_Phenopowerlaw_Gold.config @@ -14,11 +14,9 @@ plasticity phenopowerlaw (output) resistance_slip (output) shearrate_slip (output) resolvedstress_slip -(output) totalshear (output) resistance_twin (output) shearrate_twin (output) resolvedstress_twin -(output) totalvolfrac_twin lattice_structure fcc Nslip 12 # per family diff --git a/examples/ConfigFiles/Phase_Phenopowerlaw_Magnesium.config b/examples/ConfigFiles/Phase_Phenopowerlaw_Magnesium.config index ce6b06ff9..137606093 100644 --- a/examples/ConfigFiles/Phase_Phenopowerlaw_Magnesium.config +++ b/examples/ConfigFiles/Phase_Phenopowerlaw_Magnesium.config @@ -9,11 +9,9 @@ elasticity hooke (output) resistance_slip (output) shearrate_slip (output) resolvedstress_slip -(output) totalshear (output) resistance_twin (output) shearrate_twin (output) resolvedstress_twin -(output) totalvolfrac_twin lattice_structure hex covera_ratio 1.62350 # from Tromans 2011, Elastic Anisotropy of HCP Metal Crystals and Polycrystals diff --git a/examples/ConfigFiles/Phase_Phenopowerlaw_cpTi-alpha.config b/examples/ConfigFiles/Phase_Phenopowerlaw_cpTi-alpha.config index 64ecbca25..565cf02d9 100644 --- a/examples/ConfigFiles/Phase_Phenopowerlaw_cpTi-alpha.config +++ b/examples/ConfigFiles/Phase_Phenopowerlaw_cpTi-alpha.config @@ -5,11 +5,9 @@ elasticity hooke # (output) resistance_slip # (output) shearrate_slip # (output) resolvedstress_slip -# (output) totalshear # (output) resistance_twin # (output) shearrate_twin # (output) resolvedstress_twin -# (output) totalvolfrac_twin lattice_structure hex covera_ratio 1.587 diff --git a/examples/ConfigFiles/Phase_Phenopowerlaw_multiField.config b/examples/ConfigFiles/Phase_Phenopowerlaw_multiField.config index 05503a6e7..adad3710e 100644 --- a/examples/ConfigFiles/Phase_Phenopowerlaw_multiField.config +++ b/examples/ConfigFiles/Phase_Phenopowerlaw_multiField.config @@ -6,12 +6,10 @@ plasticity phenopowerlaw (output) shearrate_slip (output) resolvedstress_slip (output) accumulated_shear_slip -(output) totalshear (output) resistance_twin (output) shearrate_twin (output) resolvedstress_twin (output) accumulated_shear_twin -(output) totalvolfrac_twin lattice_structure fcc Nslip 12 # per family From cf4a0a69fd64b300fe161d314ee170d63386d45f Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 26 Nov 2018 09:15:44 +0100 Subject: [PATCH 05/70] Intel compiler detected use of unallocated variable bug was in for a while --- src/plastic_phenopowerlaw.f90 | 1 + 1 file changed, 1 insertion(+) diff --git a/src/plastic_phenopowerlaw.f90 b/src/plastic_phenopowerlaw.f90 index c2f5edb68..5fd8d8e3c 100644 --- a/src/plastic_phenopowerlaw.f90 +++ b/src/plastic_phenopowerlaw.f90 @@ -278,6 +278,7 @@ subroutine plastic_phenopowerlaw_init else twinActive allocate(prm%interaction_TwinTwin(0,0)) allocate(prm%xi_twin_0(0)) + allocate(prm%gamma_twin_char(0)) endif twinActive !-------------------------------------------------------------------------------------------------- From 7bb862f84ee92d19f93be7500eaffc96d6a9d416 Mon Sep 17 00:00:00 2001 From: Test User Date: Mon, 3 Dec 2018 11:31:52 +0100 Subject: [PATCH 06/70] [skip ci] updated version information after successful test of v2.0.2-1107-g17716b4f --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 4687d2404..85ef1a02a 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.2-1100-g65ff2157 +v2.0.2-1107-g17716b4f From 1d7172c971a5ca02a21d0db32a60241edf76bcb1 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Tue, 4 Dec 2018 17:05:35 -0500 Subject: [PATCH 07/70] adopted more intutitive alternative of P=-1 from Rowenhorst_etal2015 --- lib/damask/orientation.py | 102 +++++++++++++++++-------------- processing/pre/geom_fromTable.py | 8 +-- 2 files changed, 60 insertions(+), 50 deletions(-) diff --git a/lib/damask/orientation.py b/lib/damask/orientation.py index d0e549b4b..3516dab32 100644 --- a/lib/damask/orientation.py +++ b/lib/damask/orientation.py @@ -33,7 +33,7 @@ class Quaternion: All methods and naming conventions based on Rowenhorst_etal2015 Convention 1: coordinate frames are right-handed Convention 2: a rotation angle ω is taken to be positive for a counterclockwise rotation - when viewing from the end point of the rotation axis unit vector towards the origin + when viewing from the end point of the rotation axis towards the origin Convention 3: rotations will be interpreted in the passive sense Convention 4: Euler angle triplets are implemented using the Bunge convention, with the angular ranges as [0, 2π],[0, π],[0, 2π] @@ -49,11 +49,11 @@ class Quaternion: def __init__(self, quatArray = [1.0,0.0,0.0,0.0]): - """Initializes to identity if not given""" - self.w, \ - self.x, \ - self.y, \ - self.z = quatArray + """Initializes to identity unless specified""" + (self.w, + self.x, + self.y, + self.z ) = quatArray self.homomorph() def __iter__(self): @@ -61,16 +61,15 @@ class Quaternion: return iter([self.w,self.x,self.y,self.z]) def __copy__(self): - """Create copy""" + """Copy""" Q = Quaternion([self.w,self.x,self.y,self.z]) return Q copy = __copy__ def __repr__(self): - """Readbable string""" - return 'Quaternion(real=%+.6f, imag=<%+.6f, %+.6f, %+.6f>)' % \ - (self.w, self.x, self.y, self.z) + """Readable string""" + return 'Quaternion(real=%+.6f, imag=<%+.6f, %+.6f, %+.6f>)' % (self.w, self.x, self.y, self.z) def __pow__(self, exponent): """Power""" @@ -95,6 +94,8 @@ class Quaternion: def __mul__(self, other): """Multiplication""" + # Rowenhorst_etal2015 MSMSE: value of P is selected as -1 + P = -1.0 try: # quaternion Aw = self.w Ax = self.x @@ -106,9 +107,9 @@ class Quaternion: Bz = other.z Q = Quaternion() Q.w = - Ax * Bx - Ay * By - Az * Bz + Aw * Bw - Q.x = + Ax * Bw + Ay * Bz - Az * By + Aw * Bx - Q.y = - Ax * Bz + Ay * Bw + Az * Bx + Aw * By - Q.z = + Ax * By - Ay * Bx + Az * Bw + Aw * Bz + Q.x = + Ax * Bw + Aw * Bx + P * (Ay * Bz - Az * By) + Q.y = + Ay * Bw + Aw * By + P * (Az * Bx - Ax * Bz) + Q.z = + Az * Bw + Aw * Bz + P * (Ax * By - Ay * Bx) return Q except: pass try: # vector (perform active rotation, i.e. q*v*q.conjugated) @@ -120,16 +121,14 @@ class Quaternion: Vy = other[1] Vz = other[2] - return np.array([\ - w * w * Vx + 2 * y * w * Vz - 2 * z * w * Vy + \ - x * x * Vx + 2 * y * x * Vy + 2 * z * x * Vz - \ - z * z * Vx - y * y * Vx, - 2 * x * y * Vx + y * y * Vy + 2 * z * y * Vz + \ - 2 * w * z * Vx - z * z * Vy + w * w * Vy - \ - 2 * x * w * Vz - x * x * Vy, - 2 * x * z * Vx + 2 * y * z * Vy + \ - z * z * Vz - 2 * w * y * Vx - y * y * Vz + \ - 2 * w * x * Vy - x * x * Vz + w * w * Vz ]) + A = w**2 - x**2 - y**2 - z**2 + B = 2.0*(x*Vx + y*Vy + z*Vz) + + return np.array([ + A*Vx + B*x + 2*P*w * (y*Vz - z*Vy), + A*Vy + B*y + 2*P*w * (z*Vx - x*Vz), + A*Vz + B*z + 2*P*w * (x*Vy - y*Vx), + ]) except: pass try: # scalar Q = self.copy() @@ -143,6 +142,8 @@ class Quaternion: def __imul__(self, other): """In-place multiplication""" + # Rowenhorst_etal2015 MSMSE: value of P is selected as -1 + P = -1.0 try: # Quaternion Aw = self.w Ax = self.x @@ -153,9 +154,9 @@ class Quaternion: By = other.y Bz = other.z self.w = - Ax * Bx - Ay * By - Az * Bz + Aw * Bw - self.x = + Ax * Bw + Ay * Bz - Az * By + Aw * Bx - self.y = - Ax * Bz + Ay * Bw + Az * Bx + Aw * By - self.z = + Ax * By - Ay * Bx + Az * Bw + Aw * Bz + self.x = + Ax * Bw + Aw * Bx + P * (Ay * Bz - Az * By) + self.y = + Ay * Bw + Aw * By + P * (Az * Bx - Ax * Bz) + self.z = + Az * Bw + Aw * Bz + P * (Ax * By - Ay * Bx) except: pass return self @@ -314,13 +315,15 @@ class Quaternion: def asM(self): # to find Averaging Quaternions (see F. Landis Markley et al.) return np.outer([i for i in self],[i for i in self]) - + def asMatrix(self): + # Rowenhorst_etal2015 MSMSE: value of P is selected as -1 + P = -1.0 qbarhalf = 0.5*(self.w**2 - self.x**2 - self.y**2 - self.z**2) return 2.0*np.array( - [[ qbarhalf + self.x**2 , self.x*self.y - self.w*self.z, self.x*self.z + self.w*self.y], - [ self.x*self.y + self.w*self.z, qbarhalf + self.y**2 , self.y*self.z - self.w*self.x], - [ self.x*self.z - self.w*self.y, self.y*self.z + self.w*self.x, qbarhalf + self.z**2 ], + [[ qbarhalf + self.x**2 , self.x*self.y -P* self.w*self.z, self.x*self.z +P* self.w*self.y], + [ self.x*self.y +P* self.w*self.z, qbarhalf + self.y**2 , self.y*self.z -P* self.w*self.x], + [ self.x*self.z -P* self.w*self.y, self.y*self.z +P* self.w*self.x, qbarhalf + self.z**2 ], ]) def asAngleAxis(self, @@ -346,18 +349,20 @@ class Quaternion: def asEulers(self, degrees = False): """Orientation as Bunge-Euler angles.""" - q03 = self.w**2+self.z**2 - q12 = self.x**2+self.y**2 + # Rowenhorst_etal2015 MSMSE: value of P is selected as -1 + P = -1.0 + q03 = self.w**2 + self.z**2 + q12 = self.x**2 + self.y**2 chi = np.sqrt(q03*q12) if abs(chi) < 1e-10 and abs(q12) < 1e-10: - eulers = np.array([math.atan2(-2*self.w*self.z,self.w**2-self.z**2),0,0]) + eulers = np.array([math.atan2(-2*P*self.w*self.z,self.w**2-self.z**2),0,0]) elif abs(chi) < 1e-10 and abs(q03) < 1e-10: - eulers = np.array([math.atan2( 2*self.x*self.y,self.x**2-self.y**2),np.pi,0]) + eulers = np.array([math.atan2( 2 *self.x*self.y,self.x**2-self.y**2),np.pi,0]) else: - eulers = np.array([math.atan2((self.x*self.z-self.w*self.y)/chi,(-self.w*self.x-self.y*self.z)/chi), + eulers = np.array([math.atan2((self.x*self.z-P*self.w*self.y)/chi,(-P*self.w*self.x-self.y*self.z)/chi), math.atan2(2*chi,q03-q12), - math.atan2((self.w*self.y+self.x*self.z)/chi,( self.y*self.z-self.w*self.x)/chi), + math.atan2((P*self.w*self.y+self.x*self.z)/chi,( self.y*self.z-P*self.w*self.x)/chi), ]) return np.degrees(eulers) if degrees else eulers @@ -371,8 +376,9 @@ class Quaternion: @classmethod def fromRandom(cls,randomSeed = None): + import binascii if randomSeed is None: - randomSeed = int(os.urandom(4).encode('hex'), 16) + randomSeed = int(binascii.hexlify(os.urandom(4)),16) np.random.seed(randomSeed) r = np.random.random(3) w = math.cos(2.0*math.pi*r[0])*math.sqrt(r[2]) @@ -420,10 +426,12 @@ class Quaternion: c = np.cos(0.5*eulers[1]) s = np.sin(0.5*eulers[1]) - w = c * np.cos(sigma) - x = -s * np.cos(delta) - y = -s * np.sin(delta) - z = -c * np.sin(sigma) + # Rowenhorst_etal2015 MSMSE: value of P is selected as -1 + P = -1.0 + w = c * np.cos(sigma) + x = -P * s * np.cos(delta) + y = -P * s * np.sin(delta) + z = -P * c * np.sin(sigma) return cls([w,x,y,z]) @@ -435,16 +443,18 @@ class Quaternion: if m.shape != (3,3) and np.prod(m.shape) == 9: m = m.reshape(3,3) - w = 0.5*math.sqrt(1.+m[0,0]+m[1,1]+m[2,2]) - x = 0.5*math.sqrt(1.+m[0,0]-m[1,1]-m[2,2]) - y = 0.5*math.sqrt(1.-m[0,0]+m[1,1]-m[2,2]) - z = 0.5*math.sqrt(1.-m[0,0]-m[1,1]+m[2,2]) + # Rowenhorst_etal2015 MSMSE: value of P is selected as -1 + P = -1.0 + w = 0.5*math.sqrt(1.+m[0,0]+m[1,1]+m[2,2]) + x = P*0.5*math.sqrt(1.+m[0,0]-m[1,1]-m[2,2]) + y = P*0.5*math.sqrt(1.-m[0,0]+m[1,1]-m[2,2]) + z = P*0.5*math.sqrt(1.-m[0,0]-m[1,1]+m[2,2]) x *= -1 if m[2,1] < m[1,2] else 1 y *= -1 if m[0,2] < m[2,0] else 1 z *= -1 if m[1,0] < m[0,1] else 1 - return cls( np.array([w,x,y,z])/math.sqrt(w**2 + x**2 + y**2 + z**2)) + return cls(np.array([w,x,y,z])/math.sqrt(w**2 + x**2 + y**2 + z**2)) @classmethod diff --git a/processing/pre/geom_fromTable.py b/processing/pre/geom_fromTable.py index a2f71d0c9..3cdd1b2e6 100755 --- a/processing/pre/geom_fromTable.py +++ b/processing/pre/geom_fromTable.py @@ -43,7 +43,7 @@ parser.add_option('-e', '--eulers', parser.add_option('-d', '--degrees', dest = 'degrees', action = 'store_true', - help = 'angles are given in degrees [%default]') + help = 'all angles are in degrees') parser.add_option('-m', '--matrix', dest = 'matrix', type = 'string', metavar = 'string', @@ -71,7 +71,7 @@ parser.add_option('--axes', parser.add_option('-s', '--symmetry', dest = 'symmetry', action = 'extend', metavar = '', - help = 'crystal symmetry %default {{{}}} '.format(', '.join(damask.Symmetry.lattices[1:]))) + help = 'crystal symmetry of each phase %default {{{}}} '.format(', '.join(damask.Symmetry.lattices[1:]))) parser.add_option('--homogenization', dest = 'homogenization', type = 'int', metavar = 'int', @@ -234,7 +234,7 @@ for name in filenames: o = damask.Orientation(Eulers = myData[colOri:colOri+3]*toRadians, symmetry = mySym) elif inputtype == 'matrix': - o = damask.Orientation(matrix = myData[colOri:colOri+9].reshape(3,3).transpose(), + o = damask.Orientation(matrix = myData[colOri:colOri+9].reshape(3,3), symmetry = mySym) elif inputtype == 'frame': o = damask.Orientation(matrix = np.hstack((myData[colOri[0]:colOri[0]+3], @@ -246,7 +246,7 @@ for name in filenames: o = damask.Orientation(quaternion = myData[colOri:colOri+4], symmetry = mySym) - cos_disorientations = -np.ones(1,dtype='f') # largest possible disorientation + cos_disorientations = -np.ones(1,dtype=float) # largest possible disorientation closest_grain = -1 # invalid neighbor if options.tolerance > 0.0: # only try to compress orientations if asked to From a7554891a43676a8077e20d054c74b66178cf970 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Tue, 4 Dec 2018 19:20:24 -0500 Subject: [PATCH 08/70] changed internal quaternion representation to q,p and simplified math --- lib/damask/orientation.py | 289 ++++++++++++++------------------------ 1 file changed, 106 insertions(+), 183 deletions(-) diff --git a/lib/damask/orientation.py b/lib/damask/orientation.py index 3516dab32..9910b7a25 100644 --- a/lib/damask/orientation.py +++ b/lib/damask/orientation.py @@ -48,48 +48,42 @@ class Quaternion: """ def __init__(self, - quatArray = [1.0,0.0,0.0,0.0]): + quat = None, + q = 1.0, + p = np.zeros(3,dtype=float)): """Initializes to identity unless specified""" - (self.w, - self.x, - self.y, - self.z ) = quatArray + self.q = quat[0] if quat is not None else q + self.p = np.array(quat[1:4]) if quat is not None else p self.homomorph() def __iter__(self): """Components""" - return iter([self.w,self.x,self.y,self.z]) + return iter(self.asList()) def __copy__(self): """Copy""" - Q = Quaternion([self.w,self.x,self.y,self.z]) + Q = Quaternion(q=self.q,p=self.p) return Q copy = __copy__ def __repr__(self): """Readable string""" - return 'Quaternion(real=%+.6f, imag=<%+.6f, %+.6f, %+.6f>)' % (self.w, self.x, self.y, self.z) + return 'Quaternion(real={q:+.6f}, imag=<{p[0]:+.6f}, {p[1]:+.6f}, {p[2]:+.6f}>)'.format(q=self.q,p=self.p) def __pow__(self, exponent): """Power""" - omega = math.acos(self.w) - vRescale = math.sin(exponent*omega)/math.sin(omega) Q = Quaternion() - Q.w = math.cos(exponent*omega) - Q.x = self.x * vRescale - Q.y = self.y * vRescale - Q.z = self.z * vRescale + omega = math.acos(self.q) + Q.q = math.cos(exponent*omega) + Q.p = self.p * math.sin(exponent*omega)/math.sin(omega) return Q def __ipow__(self, exponent): """In-place power""" - omega = math.acos(self.w) - vRescale = math.sin(exponent*omega)/math.sin(omega) - self.w = np.cos(exponent*omega) - self.x *= vRescale - self.y *= vRescale - self.z *= vRescale + omega = math.acos(self.q[0]) + self.q = math.cos(exponent*omega) + self.p *= math.sin(exponent*omega)/math.sin(omega) return self def __mul__(self, other): @@ -97,45 +91,20 @@ class Quaternion: # Rowenhorst_etal2015 MSMSE: value of P is selected as -1 P = -1.0 try: # quaternion - Aw = self.w - Ax = self.x - Ay = self.y - Az = self.z - Bw = other.w - Bx = other.x - By = other.y - Bz = other.z Q = Quaternion() - Q.w = - Ax * Bx - Ay * By - Az * Bz + Aw * Bw - Q.x = + Ax * Bw + Aw * Bx + P * (Ay * Bz - Az * By) - Q.y = + Ay * Bw + Aw * By + P * (Az * Bx - Ax * Bz) - Q.z = + Az * Bw + Aw * Bz + P * (Ax * By - Ay * Bx) + Q.q = self.q*other.q - np.dot(self.p,other.p) + Q.p = self.q*other.p + other.q*self.p + P * np.cross(self.p,other.p) return Q except: pass - try: # vector (perform active rotation, i.e. q*v*q.conjugated) - w = self.w - x = self.x - y = self.y - z = self.z - Vx = other[0] - Vy = other[1] - Vz = other[2] - - A = w**2 - x**2 - y**2 - z**2 - B = 2.0*(x*Vx + y*Vy + z*Vz) - - return np.array([ - A*Vx + B*x + 2*P*w * (y*Vz - z*Vy), - A*Vy + B*y + 2*P*w * (z*Vx - x*Vz), - A*Vz + B*z + 2*P*w * (x*Vy - y*Vx), - ]) + try: # vector (perform passive rotation) + return (self.q*self.q - np.dot(self.p,self.p)) * np.array(other[:3]) \ + + 2.0*np.dot(self.p,other[:3]) * self.p \ + + 2.0*P*self.q * np.cross(self.p,other[:3]) except: pass try: # scalar Q = self.copy() - Q.w *= other - Q.x *= other - Q.y *= other - Q.z *= other + Q.q *= other + Q.p *= other return Q except: return self.copy() @@ -145,69 +114,49 @@ class Quaternion: # Rowenhorst_etal2015 MSMSE: value of P is selected as -1 P = -1.0 try: # Quaternion - Aw = self.w - Ax = self.x - Ay = self.y - Az = self.z - Bw = other.w - Bx = other.x - By = other.y - Bz = other.z - self.w = - Ax * Bx - Ay * By - Az * Bz + Aw * Bw - self.x = + Ax * Bw + Aw * Bx + P * (Ay * Bz - Az * By) - self.y = + Ay * Bw + Aw * By + P * (Az * Bx - Ax * Bz) - self.z = + Az * Bw + Aw * Bz + P * (Ax * By - Ay * Bx) + 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) except: pass return self def __div__(self, other): """Division""" if isinstance(other, (int,float)): - w = self.w / other - x = self.x / other - y = self.y / other - z = self.z / other - return self.__class__([w,x,y,z]) + q = self.q / other + p = self.p / other + return self.__class__(q=q,p=p) else: return NotImplemented def __idiv__(self, other): """In-place division""" if isinstance(other, (int,float)): - self.w /= other - self.x /= other - self.y /= other - self.z /= other + self.q /= other + self.p /= other return self def __add__(self, other): """Addition""" if isinstance(other, Quaternion): - w = self.w + other.w - x = self.x + other.x - y = self.y + other.y - z = self.z + other.z - return self.__class__([w,x,y,z]) + q = self.q + other.q + p = self.p + other.p + return self.__class__(q=q,p=p) else: return NotImplemented def __iadd__(self, other): """In-place addition""" if isinstance(other, Quaternion): - self.w += other.w - self.x += other.x - self.y += other.y - self.z += other.z + self.q += other.q + self.p += other.p return self def __sub__(self, other): """Subtraction""" if isinstance(other, Quaternion): Q = self.copy() - Q.w -= other.w - Q.x -= other.x - Q.y -= other.y - Q.z -= other.z + Q.q -= other.q + Q.p -= other.p return Q else: return self.copy() @@ -215,40 +164,25 @@ class Quaternion: def __isub__(self, other): """In-place subtraction""" if isinstance(other, Quaternion): - self.w -= other.w - self.x -= other.x - self.y -= other.y - self.z -= other.z + self.q -= other.q + self.p -= other.p return self def __neg__(self): """Additive inverse""" - self.w = -self.w - self.x = -self.x - self.y = -self.y - self.z = -self.z + self.q = -self.q + self.p = -self.p return self def __abs__(self): """Norm""" - return math.sqrt(self.w ** 2 + \ - self.x ** 2 + \ - self.y ** 2 + \ - self.z ** 2) + return math.sqrt(self.q ** 2 + np.dot(self.p,self.p)) magnitude = __abs__ def __eq__(self,other): """Equal at e-8 precision""" - return (abs(self.w-other.w) < 1e-8 and \ - abs(self.x-other.x) < 1e-8 and \ - abs(self.y-other.y) < 1e-8 and \ - abs(self.z-other.z) < 1e-8) \ - or \ - (abs(-self.w-other.w) < 1e-8 and \ - abs(-self.x-other.x) < 1e-8 and \ - abs(-self.y-other.y) < 1e-8 and \ - abs(-self.z-other.z) < 1e-8) + return (self-other).magnitude() < 1e-8 or (-self-other).magnitude() < 1e-8 def __ne__(self,other): """Not equal at e-8 precision""" @@ -259,16 +193,11 @@ class Quaternion: return (self.Rodrigues()>other.Rodrigues()) - (self.Rodrigues() 1: + if self.q > 1.: self.normalize() - s = math.sqrt(1. - self.w**2) - x = 2*self.w**2 - 1. - y = 2*self.w * s + s = math.sqrt(1. - self.q**2) + x = 2*self.q**2 - 1. + y = 2*self.q * s angle = math.atan2(y,x) if angle < 0.0: @@ -341,28 +272,28 @@ class Quaternion: s *= -1. return (np.degrees(angle) if degrees else angle, - np.array([1.0, 0.0, 0.0] if np.abs(angle) < 1e-6 else [self.x / s, self.y / s, self.z / s])) + np.array([1.0, 0.0, 0.0] if np.abs(angle) < 1e-6 else self.p / s)) def asRodrigues(self): - return np.inf*np.ones(3) if self.w == 0.0 else np.array([self.x, self.y, self.z])/self.w + return np.inf*np.ones(3) if self.q == 0.0 else self.p/self.q def asEulers(self, degrees = False): """Orientation as Bunge-Euler angles.""" # Rowenhorst_etal2015 MSMSE: value of P is selected as -1 P = -1.0 - q03 = self.w**2 + self.z**2 - q12 = self.x**2 + self.y**2 + q03 = self.q**2 + self.p[2]**2 + q12 = self.p[0]**2 + self.p[1]**2 chi = np.sqrt(q03*q12) if abs(chi) < 1e-10 and abs(q12) < 1e-10: - eulers = np.array([math.atan2(-2*P*self.w*self.z,self.w**2-self.z**2),0,0]) + eulers = np.array([math.atan2(-2*P*self.q*self.p[2],self.q**2-self.p[2]**2),0,0]) elif abs(chi) < 1e-10 and abs(q03) < 1e-10: - eulers = np.array([math.atan2( 2 *self.x*self.y,self.x**2-self.y**2),np.pi,0]) + eulers = np.array([math.atan2( 2 *self.p[0]*self.p[1],self.p[0]**2-self.p[1]**2),np.pi,0]) else: - eulers = np.array([math.atan2((self.x*self.z-P*self.w*self.y)/chi,(-P*self.w*self.x-self.y*self.z)/chi), + eulers = np.array([math.atan2((self.p[0]*self.p[2]-P*self.q*self.p[1])/chi,(-P*self.q*self.p[0]-self.p[1]*self.p[2])/chi), math.atan2(2*chi,q03-q12), - math.atan2((P*self.w*self.y+self.x*self.z)/chi,( self.y*self.z-P*self.w*self.x)/chi), + math.atan2((P*self.q*self.p[1]+self.p[0]*self.p[2])/chi,( self.p[1]*self.p[2]-P*self.q*self.p[0])/chi), ]) return np.degrees(eulers) if degrees else eulers @@ -385,7 +316,7 @@ class Quaternion: x = math.sin(2.0*math.pi*r[1])*math.sqrt(1.0-r[2]) y = math.cos(2.0*math.pi*r[1])*math.sqrt(1.0-r[2]) z = math.sin(2.0*math.pi*r[0])*math.sqrt(r[2]) - return cls([w,x,y,z]) + return cls(quat=[w,x,y,z]) @classmethod @@ -393,9 +324,7 @@ class Quaternion: if not isinstance(rodrigues, np.ndarray): rodrigues = np.array(rodrigues) halfangle = math.atan(np.linalg.norm(rodrigues)) c = math.cos(halfangle) - w = c - x,y,z = rodrigues/c - return cls([w,x,y,z]) + return cls(q=c,p=rodrigues/c) @classmethod @@ -403,22 +332,19 @@ class Quaternion: angle, axis, degrees = False): - if not isinstance(axis, np.ndarray): axis = np.array(axis,dtype='d') + if not isinstance(axis, np.ndarray): axis = np.array(axis,dtype=float) axis = axis.astype(float)/np.linalg.norm(axis) angle = np.radians(angle) if degrees else angle s = math.sin(0.5 * angle) - w = math.cos(0.5 * angle) - x = axis[0] * s - y = axis[1] * s - z = axis[2] * s - return cls([w,x,y,z]) + c = math.cos(0.5 * angle) + return cls(q=c,p=axis*s) @classmethod def fromEulers(cls, eulers, degrees = False): - if not isinstance(eulers, np.ndarray): eulers = np.array(eulers,dtype='d') + if not isinstance(eulers, np.ndarray): eulers = np.array(eulers,dtype=float) eulers = np.radians(eulers) if degrees else eulers sigma = 0.5*(eulers[0]+eulers[2]) @@ -432,7 +358,7 @@ class Quaternion: x = -P * s * np.cos(delta) y = -P * s * np.sin(delta) z = -P * c * np.sin(sigma) - return cls([w,x,y,z]) + return cls(quat=[w,x,y,z]) # Modified Method to calculate Quaternion from Orientation Matrix, @@ -454,7 +380,7 @@ class Quaternion: y *= -1 if m[0,2] < m[2,0] else 1 z *= -1 if m[1,0] < m[0,1] else 1 - return cls(np.array([w,x,y,z])/math.sqrt(w**2 + x**2 + y**2 + z**2)) + return cls(quat=np.array([w,x,y,z])/math.sqrt(w**2 + x**2 + y**2 + z**2)) @classmethod @@ -468,36 +394,30 @@ class Quaternion: assert isinstance(q1, Quaternion) and isinstance(q2, Quaternion) Q = cls() - costheta = q1.w * q2.w + q1.x * q2.x + q1.y * q2.y + q1.z * q2.z + costheta = q1.q*q2.q + np.dot(q1.p,q2.p) if costheta < 0.: costheta = -costheta q1 = q1.conjugated() - elif costheta > 1: - costheta = 1 + elif costheta > 1.: + costheta = 1. theta = math.acos(costheta) if abs(theta) < 0.01: - Q.w = q2.w - Q.x = q2.x - Q.y = q2.y - Q.z = q2.z + Q.q = q2.q + Q.p = q2.p return Q sintheta = math.sqrt(1.0 - costheta * costheta) if abs(sintheta) < 0.01: - Q.w = (q1.w + q2.w) * 0.5 - Q.x = (q1.x + q2.x) * 0.5 - Q.y = (q1.y + q2.y) * 0.5 - Q.z = (q1.z + q2.z) * 0.5 + Q.q = (q1.q + q2.q) * 0.5 + Q.p = (q1.p + q2.p) * 0.5 return Q - ratio1 = math.sin((1 - t) * theta) / sintheta - ratio2 = math.sin(t * theta) / sintheta + ratio1 = math.sin((1.0 - t) * theta) / sintheta + ratio2 = math.sin( t * theta) / sintheta - Q.w = q1.w * ratio1 + q2.w * ratio2 - Q.x = q1.x * ratio1 + q2.x * ratio2 - Q.y = q1.y * ratio1 + q2.y * ratio2 - Q.z = q1.z * ratio1 + q2.z * ratio2 + Q.q = q1.q * ratio1 + q2.q * ratio2 + Q.p = q1.p * ratio1 + q2.p * ratio2 return Q @@ -523,7 +443,7 @@ class Symmetry: def __repr__(self): """Readbable string""" - return '%s' % (self.lattice) + return '{}'.format(self.lattice) def __eq__(self, other): @@ -536,7 +456,7 @@ class Symmetry: def __cmp__(self,other): """Linear ordering""" - myOrder = Symmetry.lattices.index(self.lattice) + myOrder = Symmetry.lattices.index(self.lattice) otherOrder = Symmetry.lattices.index(other.lattice) return (myOrder > otherOrder) - (myOrder < otherOrder) @@ -732,7 +652,7 @@ class Symmetry: else: return True - v = np.array(vector,dtype = float) + v = np.array(vector,dtype=float) if proper: # check both improper ... theComponents = np.dot(basis['improper'],v) inSST = np.all(theComponents >= 0.0) @@ -747,10 +667,10 @@ class Symmetry: if color: # have to return color array if inSST: rgb = np.power(theComponents/np.linalg.norm(theComponents),0.5) # smoothen color ramps - rgb = np.minimum(np.ones(3,'d'),rgb) # limit to maximum intensity + rgb = np.minimum(np.ones(3,dtype=float),rgb) # limit to maximum intensity rgb /= max(rgb) # normalize to (HS)V = 1 else: - rgb = np.zeros(3,'d') + rgb = np.zeros(3,dtype=float) return (inSST,rgb) else: return inSST @@ -790,8 +710,9 @@ class Orientation: self.quaternion = Quaternion.fromRodrigues(Rodrigues) elif isinstance(quaternion, Quaternion): # based on given quaternion self.quaternion = quaternion.homomorphed() - elif isinstance(quaternion, np.ndarray) and quaternion.shape == (4,): # based on given quaternion-like array - self.quaternion = Quaternion(quaternion).homomorphed() + elif (isinstance(quaternion, np.ndarray) and quaternion.shape == (4,)) or \ + (isinstance(quaternion, list) and len(quaternion) == 4 ): # based on given quaternion-like array + self.quaternion = Quaternion(quat=quaternion).homomorphed() self.symmetry = Symmetry(symmetry) @@ -804,10 +725,12 @@ class Orientation: def __repr__(self): """Value as all implemented representations""" - return 'Symmetry: %s\n' % (self.symmetry) + \ - 'Quaternion: %s\n' % (self.quaternion) + \ - 'Matrix:\n%s\n' % ( '\n'.join(['\t'.join(map(str,self.asMatrix()[i,:])) for i in range(3)]) ) + \ - 'Bunge Eulers / deg: %s' % ('\t'.join(map(str,self.asEulers(degrees=True))) ) + return '\n'.join([ + 'Symmetry: {}'.format(self.symmetry), + 'Quaternion: {}'.format(self.quaternion), + 'Matrix:\n{}'.format( '\n'.join(['\t'.join(list(map(str,self.asMatrix()[i,:]))) for i in range(3)]) ), + 'Bunge Eulers / deg: {}'.format('\t'.join(list(map(str,self.asEulers(degrees=True)))) ), + ]) def asQuaternion(self): return self.quaternion.asList() From 49af6c70a7c4187ea7b4175cc8b96a007b7bfdb9 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 5 Dec 2018 08:56:26 +0100 Subject: [PATCH 09/70] PGI compiler complained about kind mismatch --- src/spectral_damage.f90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/spectral_damage.f90 b/src/spectral_damage.f90 index c6f2f3822..13ddf0e74 100644 --- a/src/spectral_damage.f90 +++ b/src/spectral_damage.f90 @@ -120,8 +120,8 @@ subroutine spectral_damage_init() trim(snes_type) == 'vinewtonssls') then call DMGetGlobalVector(damage_grid,lBound,ierr); CHKERRQ(ierr) call DMGetGlobalVector(damage_grid,uBound,ierr); CHKERRQ(ierr) - call VecSet(lBound,0.0,ierr); CHKERRQ(ierr) - call VecSet(uBound,1.0,ierr); CHKERRQ(ierr) + call VecSet(lBound,0.0_pReal,ierr); CHKERRQ(ierr) + call VecSet(uBound,1.0_pReal,ierr); CHKERRQ(ierr) call SNESVISetVariableBounds(damage_snes,lBound,uBound,ierr) !< variable bounds for variational inequalities like contact mechanics, damage etc. call DMRestoreGlobalVector(damage_grid,lBound,ierr); CHKERRQ(ierr) call DMRestoreGlobalVector(damage_grid,uBound,ierr); CHKERRQ(ierr) @@ -134,7 +134,7 @@ subroutine spectral_damage_init() xend = xstart + xend - 1 yend = ystart + yend - 1 zend = zstart + zend - 1 - call VecSet(solution,1.0,ierr); CHKERRQ(ierr) + call VecSet(solution,1.0_pReal,ierr); CHKERRQ(ierr) allocate(damage_current(grid(1),grid(2),grid3), source=1.0_pReal) allocate(damage_lastInc(grid(1),grid(2),grid3), source=1.0_pReal) allocate(damage_stagInc(grid(1),grid(2),grid3), source=1.0_pReal) From deedbc4fda591578356a413713df09fdc42a92b7 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 5 Dec 2018 08:59:32 +0100 Subject: [PATCH 10/70] standard conforming line continuation --- src/kinematics_thermal_expansion.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/kinematics_thermal_expansion.f90 b/src/kinematics_thermal_expansion.f90 index 3cec1da4c..3d1de3d0a 100644 --- a/src/kinematics_thermal_expansion.f90 +++ b/src/kinematics_thermal_expansion.f90 @@ -225,7 +225,7 @@ subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar3333, ipc, + lattice_thermalExpansion33(1:3,1:3,2,phase)*(T - TRef)**1 & ! linear coefficient + lattice_thermalExpansion33(1:3,1:3,3,phase)*(T - TRef)**2 & ! quadratic coefficient ) / & - (1.0_pReal \ + (1.0_pReal & + lattice_thermalExpansion33(1:3,1:3,1,phase)*(T - TRef)**1 / 1. & + lattice_thermalExpansion33(1:3,1:3,2,phase)*(T - TRef)**2 / 2. & + lattice_thermalExpansion33(1:3,1:3,3,phase)*(T - TRef)**3 / 3. & From e48fd338bcce4eb02e25ea398b45a34ee1515af8 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Wed, 5 Dec 2018 08:55:26 -0500 Subject: [PATCH 11/70] fixed shallow copy problem of quaternion.p object --- lib/damask/orientation.py | 45 +++++++++++++++++---------------------- 1 file changed, 20 insertions(+), 25 deletions(-) diff --git a/lib/damask/orientation.py b/lib/damask/orientation.py index 9910b7a25..d3cd4f94a 100644 --- a/lib/damask/orientation.py +++ b/lib/damask/orientation.py @@ -62,8 +62,7 @@ class Quaternion: def __copy__(self): """Copy""" - Q = Quaternion(q=self.q,p=self.p) - return Q + return self.__class__(q=self.q,p=self.p.copy()) copy = __copy__ @@ -91,10 +90,8 @@ class Quaternion: # Rowenhorst_etal2015 MSMSE: value of P is selected as -1 P = -1.0 try: # quaternion - Q = Quaternion() - Q.q = self.q*other.q - np.dot(self.p,other.p) - Q.p = self.q*other.p + other.q*self.p + P * np.cross(self.p,other.p) - return Q + 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)) except: pass try: # vector (perform passive rotation) return (self.q*self.q - np.dot(self.p,self.p)) * np.array(other[:3]) \ @@ -102,10 +99,8 @@ class Quaternion: + 2.0*P*self.q * np.cross(self.p,other[:3]) except: pass try: # scalar - Q = self.copy() - Q.q *= other - Q.p *= other - return Q + return self.__class__(q=self.q*other, + p=self.p*other) except: return self.copy() @@ -122,9 +117,8 @@ class Quaternion: def __div__(self, other): """Division""" if isinstance(other, (int,float)): - q = self.q / other - p = self.p / other - return self.__class__(q=q,p=p) + return self.__class__(q=self.q / other, + p=self.p / other) else: return NotImplemented @@ -138,9 +132,8 @@ class Quaternion: def __add__(self, other): """Addition""" if isinstance(other, Quaternion): - q = self.q + other.q - p = self.p + other.p - return self.__class__(q=q,p=p) + return self.__class__(q=self.q + other.q, + p=self.p + other.p) else: return NotImplemented @@ -154,12 +147,10 @@ class Quaternion: def __sub__(self, other): """Subtraction""" if isinstance(other, Quaternion): - Q = self.copy() - Q.q -= other.q - Q.p -= other.p - return Q + return self.__class__(q=self.q - other.q, + p=self.p - other.p) else: - return self.copy() + return NotImplemented def __isub__(self, other): """In-place subtraction""" @@ -190,7 +181,8 @@ class Quaternion: def __cmp__(self,other): """Linear ordering""" - return (self.Rodrigues()>other.Rodrigues()) - (self.Rodrigues() np.linalg.norm(other.asRodrigues()) else 0) \ + - (1 if np.linalg.norm(self.asRodrigues()) < np.linalg.norm(other.asRodrigues()) else 0) def magnitude_squared(self): return self.q ** 2 + np.dot(self.p,self.p) @@ -203,7 +195,8 @@ class Quaternion: def normalize(self): d = self.magnitude() if d > 0.0: - self /= d + self.q /= d + self.p /= d return self def conjugate(self): @@ -322,9 +315,11 @@ class Quaternion: @classmethod def fromRodrigues(cls, rodrigues): if not isinstance(rodrigues, np.ndarray): rodrigues = np.array(rodrigues) - halfangle = math.atan(np.linalg.norm(rodrigues)) + norm = np.linalg.norm(rodrigues) + halfangle = math.atan(norm) + s = math.sin(halfangle) c = math.cos(halfangle) - return cls(q=c,p=rodrigues/c) + return cls(q=c,p=s*rodrigues/norm) @classmethod From c0f7ae27983ed6c75e68234113e9857fce77f9e4 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Wed, 5 Dec 2018 09:07:29 -0500 Subject: [PATCH 12/70] updated scripts that still used formerly valid object properties of quaternions --- lib/damask/orientation.py | 2 +- processing/post/addGrainID.py | 4 ++-- processing/pre/geom_fromTable.py | 4 ++-- 3 files changed, 5 insertions(+), 5 deletions(-) diff --git a/lib/damask/orientation.py b/lib/damask/orientation.py index d3cd4f94a..63fb124a6 100644 --- a/lib/damask/orientation.py +++ b/lib/damask/orientation.py @@ -855,7 +855,7 @@ class Orientation: M = closest.quaternion.asM() * n if i == 0 else M + closest.quaternion.asM() * n # noqa add (multiples) of this orientation to average noqa eig, vec = np.linalg.eig(M/N) - return Orientation(quaternion = Quaternion(quatArray = np.real(vec.T[eig.argmax()])), + return Orientation(quaternion = Quaternion(quat = np.real(vec.T[eig.argmax()])), symmetry = reference.symmetry.lattice) diff --git a/processing/post/addGrainID.py b/processing/post/addGrainID.py index 45034034b..c3b98f4e6 100755 --- a/processing/post/addGrainID.py +++ b/processing/post/addGrainID.py @@ -200,9 +200,9 @@ for name in filenames: if gID != -1 and gID not in alreadyChecked: # indexed point belonging to a grain not yet tested? alreadyChecked[gID] = True # remember not to check again disorientation = o.disorientation(orientations[gID],SST = False)[0] # compare against other orientation - if disorientation.quaternion.w > cos_disorientation: # within threshold ... + if disorientation.quaternion.q > cos_disorientation: # within threshold ... candidates.append(gID) # remember as potential candidate - if disorientation.quaternion.w >= bestDisorientation.w: # ... and better than current best? + if disorientation.quaternion.q >= bestDisorientation.q: # ... and better than current best? matched = True matchedID = gID # remember that grain bestDisorientation = disorientation.quaternion diff --git a/processing/pre/geom_fromTable.py b/processing/pre/geom_fromTable.py index 3cdd1b2e6..e1157d325 100755 --- a/processing/pre/geom_fromTable.py +++ b/processing/pre/geom_fromTable.py @@ -258,7 +258,7 @@ for name in filenames: if len(grains) > 0: # check immediate neighborhood first cos_disorientations = np.array([o.disorientation(orientations[grainID], - SST = False)[0].quaternion.w \ + SST = False)[0].quaternion.q \ for grainID in grains]) # store disorientation per grainID closest_grain = np.argmax(cos_disorientations) # grain among grains with closest orientation to myself match = 'local' @@ -269,7 +269,7 @@ for name in filenames: if len(grains) > 0: cos_disorientations = np.array([o.disorientation(orientations[grainID], - SST = False)[0].quaternion.w \ + SST = False)[0].quaternion.q \ for grainID in grains]) # store disorientation per grainID closest_grain = np.argmax(cos_disorientations) # grain among grains with closest orientation to myself match = 'global' From a6d4c73de0e559f1484930f29e52777c15111b06 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Wed, 5 Dec 2018 10:20:05 -0500 Subject: [PATCH 13/70] added list of map and introduced "quat" keyword in quaternion init --- lib/damask/orientation.py | 22 ++++++++++++++-------- processing/pre/geom_addPrimitive.py | 8 ++++---- 2 files changed, 18 insertions(+), 12 deletions(-) diff --git a/lib/damask/orientation.py b/lib/damask/orientation.py index 63fb124a6..c1862e13c 100644 --- a/lib/damask/orientation.py +++ b/lib/damask/orientation.py @@ -72,15 +72,13 @@ class Quaternion: def __pow__(self, exponent): """Power""" - Q = Quaternion() omega = math.acos(self.q) - Q.q = math.cos(exponent*omega) - Q.p = self.p * math.sin(exponent*omega)/math.sin(omega) - return Q + return self.__class__(q= math.cos(exponent*omega), + p=self.p * math.sin(exponent*omega)/math.sin(omega)) def __ipow__(self, exponent): """In-place power""" - omega = math.acos(self.q[0]) + omega = math.acos(self.q) self.q = math.cos(exponent*omega) self.p *= math.sin(exponent*omega)/math.sin(omega) return self @@ -94,9 +92,17 @@ class Quaternion: p=self.q*other.p + other.q*self.p + P * np.cross(self.p,other.p)) except: pass try: # vector (perform passive rotation) - return (self.q*self.q - np.dot(self.p,self.p)) * np.array(other[:3]) \ - + 2.0*np.dot(self.p,other[:3]) * self.p \ - + 2.0*P*self.q * np.cross(self.p,other[:3]) + ( x, y, z) = self.p + (Vx,Vy,Vz) = other[0:3] + A = self.q*self.q - np.dot(self.p,self.p) + B = 2.0 * (x*Vx + y*Vy + z*Vz) + C = 2.0 * P*self.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), + ]) except: pass try: # scalar return self.__class__(q=self.q*other, diff --git a/processing/pre/geom_addPrimitive.py b/processing/pre/geom_addPrimitive.py index 35000d8bf..54de558f7 100755 --- a/processing/pre/geom_addPrimitive.py +++ b/processing/pre/geom_addPrimitive.py @@ -64,11 +64,11 @@ if options.dimension is None: parser.error('no dimension specified.') if options.angleaxis is not None: options.angleaxis = list(map(float,options.angleaxis)) - rotation = damask.Quaternion().fromAngleAxis(np.radians(options.angleaxis[0]) if options.degrees else options.angleaxis[0], - options.angleaxis[1:4]) + rotation = damask.Quaternion.fromAngleAxis(np.radians(options.angleaxis[0]) if options.degrees else options.angleaxis[0], + options.angleaxis[1:4]) elif options.quaternion is not None: - options.quaternion = map(float,options.quaternion) - rotation = damask.Quaternion(options.quaternion) + options.quaternion = list(map(float,options.quaternion)) + rotation = damask.Quaternion(quat=options.quaternion) else: rotation = damask.Quaternion() From 72304638f9ac2852a6e0206f3dce6fb3499dde64 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 6 Dec 2018 00:42:28 +0100 Subject: [PATCH 14/70] PGI fails if optional argument is not present --- src/IO.f90 | 6 +++++- src/config.f90 | 12 ++++++++++-- 2 files changed, 15 insertions(+), 3 deletions(-) diff --git a/src/IO.f90 b/src/IO.f90 index a76959a44..4b962bd1e 100644 --- a/src/IO.f90 +++ b/src/IO.f90 @@ -684,7 +684,11 @@ function IO_stringValue(string,chunkPos,myChunk,silent) logical :: warn - warn = merge(silent,.false.,present(silent)) + if (present(silent)) then + warn = silent + else + warn = .false. + endif IO_stringValue = '' valuePresent: if (myChunk > chunkPos(1) .or. myChunk < 1_pInt) then diff --git a/src/config.f90 b/src/config.f90 index 4d5a76432..7ae800f30 100644 --- a/src/config.f90 +++ b/src/config.f90 @@ -513,8 +513,12 @@ character(len=65536) function getString(this,key,defaultVal,raw) type(tPartitionedStringList), pointer :: item logical :: found, & whole + if (present(raw)) then + whole = raw + else + whole = .false. + endif - whole = merge(raw,.false.,present(raw)) ! whole string or white space splitting found = present(defaultVal) if (found) then getString = trim(defaultVal) @@ -661,7 +665,11 @@ function getStrings(this,key,defaultVal,requiredShape,raw) cumulative cumulative = (key(1:1) == '(' .and. key(len_trim(key):len_trim(key)) == ')') - whole = merge(raw,.false.,present(raw)) + if (present(raw)) then + whole = raw + else + whole = .false. + endif found = .false. item => this From 57683566423bf2c5a3b1e3324620639fd0e2b8f1 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 6 Dec 2018 01:11:41 +0100 Subject: [PATCH 15/70] PGI compatible --- src/IO.f90 | 4 +++- src/math.f90 | 11 ++++------- 2 files changed, 7 insertions(+), 8 deletions(-) diff --git a/src/IO.f90 b/src/IO.f90 index 4b962bd1e..af59b11b9 100644 --- a/src/IO.f90 +++ b/src/IO.f90 @@ -191,7 +191,9 @@ recursive function IO_recursiveRead(fileName,cnt) result(fileContent) l,i, & myStat - if (merge(cnt,0_pInt,present(cnt))>10_pInt) call IO_error(106_pInt,ext_msg=trim(fileName)) + if (present(cnt)) then + if (cnt>10_pInt) call IO_error(106_pInt,ext_msg=trim(fileName)) + endif !-------------------------------------------------------------------------------------------------- ! read data as stream diff --git a/src/math.f90 b/src/math.f90 index 440ee5303..725c0446e 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -302,7 +302,7 @@ subroutine math_check endif end subroutine math_check - + !-------------------------------------------------------------------------------------------------- !> @brief Quicksort algorithm for two-dimensional integer arrays @@ -2625,12 +2625,9 @@ real(pReal) pure function math_clip(a, left, right) real(pReal), intent(in) :: a real(pReal), intent(in), optional :: left, right - - math_clip = min ( & - max (merge(left, -huge(a), present(left)), a), & - merge(right, huge(a), present(right)) & - ) - + math_clip = a + if (present(left)) math_clip = max(left,math_clip) + if (present(right)) math_clip = min(right,math_clip) if (present(left) .and. present(right)) & math_clip = merge (IEEE_value(1.0_pReal,IEEE_quiet_NaN),math_clip, left>right) From 95e6a41807a433e8abbc40b7f6f34d1cca544e89 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 6 Dec 2018 06:25:29 +0100 Subject: [PATCH 16/70] disabled test for post processing since orientation conversion give unexpected results at the moment --- .gitlab-ci.yml | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 06f4356b0..368ca94db 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -158,12 +158,12 @@ Post_AverageDown: - master - release -Post_General: - stage: postprocessing - script: PostProcessing/test.py - except: - - master - - release +#Post_General: +# stage: postprocessing +# script: PostProcessing/test.py +# except: +# - master +# - release Post_GeometryReconstruction: stage: postprocessing From ee7e4ee0d983ee6b988cc822e658360b82dd4a82 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 6 Dec 2018 06:27:07 +0100 Subject: [PATCH 17/70] python3 has different encoding syntax --- processing/pre/hybridIA_linODFsampling.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/processing/pre/hybridIA_linODFsampling.py b/processing/pre/hybridIA_linODFsampling.py index acedc3cca..cf1a473cf 100755 --- a/processing/pre/hybridIA_linODFsampling.py +++ b/processing/pre/hybridIA_linODFsampling.py @@ -244,7 +244,7 @@ for name in filenames: continue damask.util.report(scriptName,name) - randomSeed = int(os.urandom(4).encode('hex'), 16) if options.randomSeed is None else options.randomSeed # random seed per file for second phase + randomSeed = int(os.urandom(4).hex(), 16) if options.randomSeed is None else options.randomSeed # random seed per file for second phase random.seed(randomSeed) # ------------------------------------------ read header and data --------------------------------- From 2349341e25e9b2767e1028215b526db4b27bd1c5 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 6 Dec 2018 22:23:34 +0100 Subject: [PATCH 18/70] texture component test also needs a closer look --- .gitlab-ci.yml | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 368ca94db..99daef5e5 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -364,12 +364,12 @@ Phenopowerlaw_singleSlip: - master - release -TextureComponents: - stage: spectral - script: TextureComponents/test.py - except: - - master - - release +#TextureComponents: +# stage: spectral +# script: TextureComponents/test.py +# except: +# - master +# - release ################################################################################################### From 590f83a944a2ad5dac811d0c2f554ab19c507f72 Mon Sep 17 00:00:00 2001 From: Test User Date: Fri, 7 Dec 2018 03:32:36 +0100 Subject: [PATCH 19/70] [skip ci] updated version information after successful test of v2.0.2-1122-g2349341e --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 85ef1a02a..3ad5cdbde 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.2-1107-g17716b4f +v2.0.2-1122-g2349341e From d7f505f24e3369f2e09c28b089b9fce0a8f41e79 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Fri, 7 Dec 2018 10:34:45 -0500 Subject: [PATCH 20/70] making sure that roots are taken from strictly non-negative summation results --- lib/damask/orientation.py | 20 +++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/lib/damask/orientation.py b/lib/damask/orientation.py index c1862e13c..7b90c5c09 100644 --- a/lib/damask/orientation.py +++ b/lib/damask/orientation.py @@ -311,10 +311,12 @@ class Quaternion: randomSeed = int(binascii.hexlify(os.urandom(4)),16) np.random.seed(randomSeed) r = np.random.random(3) - w = math.cos(2.0*math.pi*r[0])*math.sqrt(r[2]) - x = math.sin(2.0*math.pi*r[1])*math.sqrt(1.0-r[2]) - y = math.cos(2.0*math.pi*r[1])*math.sqrt(1.0-r[2]) - z = math.sin(2.0*math.pi*r[0])*math.sqrt(r[2]) + A = math.sqrt(max(0.0,r[2])) + B = math.sqrt(max(0.0,1.0-r[2])) + w = math.cos(2.0*math.pi*r[0])*A + x = math.sin(2.0*math.pi*r[1])*B + y = math.cos(2.0*math.pi*r[1])*B + z = math.sin(2.0*math.pi*r[0])*A return cls(quat=[w,x,y,z]) @@ -372,10 +374,10 @@ class Quaternion: # Rowenhorst_etal2015 MSMSE: value of P is selected as -1 P = -1.0 - w = 0.5*math.sqrt(1.+m[0,0]+m[1,1]+m[2,2]) - x = P*0.5*math.sqrt(1.+m[0,0]-m[1,1]-m[2,2]) - y = P*0.5*math.sqrt(1.-m[0,0]+m[1,1]-m[2,2]) - z = P*0.5*math.sqrt(1.-m[0,0]-m[1,1]+m[2,2]) + w = 0.5*math.sqrt(max(0.0,1.0+m[0,0]+m[1,1]+m[2,2])) + x = P*0.5*math.sqrt(max(0.0,1.0+m[0,0]-m[1,1]-m[2,2])) + y = P*0.5*math.sqrt(max(0.0,1.0-m[0,0]+m[1,1]-m[2,2])) + z = P*0.5*math.sqrt(max(0.0,1.0-m[0,0]-m[1,1]+m[2,2])) x *= -1 if m[2,1] < m[1,2] else 1 y *= -1 if m[0,2] < m[2,0] else 1 @@ -668,7 +670,7 @@ class Symmetry: if color: # have to return color array if inSST: rgb = np.power(theComponents/np.linalg.norm(theComponents),0.5) # smoothen color ramps - rgb = np.minimum(np.ones(3,dtype=float),rgb) # limit to maximum intensity + rgb = np.minimum(np.ones(3,dtype=float),rgb) # limit to maximum intensity rgb /= max(rgb) # normalize to (HS)V = 1 else: rgb = np.zeros(3,dtype=float) From 3d7a73a0dee6ec30d57a7e078a05d523b9f06088 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Fri, 7 Dec 2018 22:02:30 -0500 Subject: [PATCH 21/70] added option to return angleAxis as flat list --- lib/damask/orientation.py | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/lib/damask/orientation.py b/lib/damask/orientation.py index 7b90c5c09..6fd0993ea 100644 --- a/lib/damask/orientation.py +++ b/lib/damask/orientation.py @@ -257,7 +257,8 @@ class Quaternion: ]) def asAngleAxis(self, - degrees = False): + degrees = False, + flat = False): if self.q > 1.: self.normalize() @@ -270,8 +271,12 @@ class Quaternion: angle *= -1. s *= -1. - return (np.degrees(angle) if degrees else angle, - np.array([1.0, 0.0, 0.0] if np.abs(angle) < 1e-6 else self.p / s)) + if flat: + return np.hstack((np.degrees(angle) if degrees else angle, + np.array([1.0, 0.0, 0.0] if np.abs(angle) < 1e-6 else self.p / s))) + else: + return (np.degrees(angle) if degrees else angle, + np.array([1.0, 0.0, 0.0] if np.abs(angle) < 1e-6 else self.p / s)) def asRodrigues(self): return np.inf*np.ones(3) if self.q == 0.0 else self.p/self.q @@ -749,8 +754,9 @@ class Orientation: rodrigues = property(asRodrigues) def asAngleAxis(self, - degrees = False): - return self.quaternion.asAngleAxis(degrees) + degrees = False, + flat = False): + return self.quaternion.asAngleAxis(degrees,flat) angleAxis = property(asAngleAxis) def asMatrix(self): From ac45a9e2d33c2e10ea037c19632c695df92c6bff Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Fri, 7 Dec 2018 22:03:14 -0500 Subject: [PATCH 22/70] more output options, frame input changed to "lab expressed in crystal" --- processing/post/addOrientations.py | 76 +++++++++++++++--------------- 1 file changed, 38 insertions(+), 38 deletions(-) diff --git a/processing/post/addOrientations.py b/processing/post/addOrientations.py index b4bda4770..b4a8ed07c 100755 --- a/processing/post/addOrientations.py +++ b/processing/post/addOrientations.py @@ -21,19 +21,22 @@ Additional (globally fixed) rotations of the lab frame and/or crystal frame can """, version = scriptID) -outputChoices = ['quaternion','rodrigues','eulers'] +outputChoices = { + 'quaternion': ['quat',4], + 'rodrigues': ['rodr',3], + 'eulers': ['eulr',3], + 'matrix': ['mtrx',9], + 'angleaxis': ['aaxs',4], + } + parser.add_option('-o', '--output', dest = 'output', action = 'extend', metavar = '', help = 'output orientation formats {{{}}}'.format(', '.join(outputChoices))) -parser.add_option('-s', '--symmetry', - dest = 'symmetry', - type = 'choice', choices = damask.Symmetry.lattices[1:], metavar='string', - help = 'crystal symmetry [%default] {{{}}} '.format(', '.join(damask.Symmetry.lattices[1:]))) parser.add_option('-d', '--degrees', dest = 'degrees', action = 'store_true', - help = 'angles are given in degrees [%default]') + help = 'all angles in degrees') parser.add_option('-R', '--labrotation', dest='labrotation', type = 'float', nargs = 4, metavar = ' '.join(['float']*4), @@ -58,21 +61,20 @@ parser.add_option( '--quaternion', dest = 'quaternion', type = 'string', metavar = 'string', help = 'quaternion label') -parser.add_option('-a', - dest = 'a', +parser.add_option('-x', + dest = 'x', type = 'string', metavar = 'string', - help = 'crystal frame a vector label') -parser.add_option('-b', - dest = 'b', + help = 'label of lab x vector (expressed in crystal coords)') +parser.add_option('-y', + dest = 'y', type = 'string', metavar = 'string', - help = 'crystal frame b vector label') -parser.add_option('-c', - dest = 'c', + help = 'label of lab y vector (expressed in crystal coords)') +parser.add_option('-z', + dest = 'z', type = 'string', metavar = 'string', - help = 'crystal frame c vector label') + help = 'label of lab z vector (expressed in crystal coords)') parser.set_defaults(output = [], - symmetry = damask.Symmetry.lattices[-1], labrotation = (0.,1.,1.,1.), # no rotation about 1,1,1 crystalrotation = (0.,1.,1.,1.), # no rotation about 1,1,1 degrees = False, @@ -86,9 +88,9 @@ if options.output == [] or (not set(options.output).issubset(set(outputChoices)) input = [options.eulers is not None, options.rodrigues is not None, - options.a is not None and \ - options.b is not None and \ - options.c is not None, + options.x is not None and \ + options.y is not None and \ + options.z is not None, options.matrix is not None, options.quaternion is not None, ] @@ -97,13 +99,13 @@ if np.sum(input) != 1: parser.error('needs exactly one input format.') (label,dim,inputtype) = [(options.eulers,3,'eulers'), (options.rodrigues,3,'rodrigues'), - ([options.a,options.b,options.c],[3,3,3],'frame'), + ([options.x,options.y,options.z],[3,3,3],'frame'), (options.matrix,9,'matrix'), (options.quaternion,4,'quaternion'), ][np.where(input)[0][0]] # select input label that was requested toRadians = math.pi/180.0 if options.degrees else 1.0 # rescale degrees to radians -r = damask.Quaternion().fromAngleAxis(toRadians*options.crystalrotation[0],options.crystalrotation[1:]) # crystal frame rotation -R = damask.Quaternion().fromAngleAxis(toRadians*options. labrotation[0],options. labrotation[1:]) # lab frame rotation +r = damask.Quaternion.fromAngleAxis(toRadians*options.crystalrotation[0],options.crystalrotation[1:]) # crystal frame rotation +R = damask.Quaternion.fromAngleAxis(toRadians*options. labrotation[0],options. labrotation[1:]) # lab frame rotation # --- loop over input files ------------------------------------------------------------------------ @@ -137,32 +139,28 @@ for name in filenames: table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:])) for output in options.output: - if output == 'quaternion': table.labels_append(['{}_{}_{}({})'.format(i+1,'quat',options.symmetry,label) for i in range(4)]) - elif output == 'rodrigues': table.labels_append(['{}_{}_{}({})'.format(i+1,'rodr',options.symmetry,label) for i in range(3)]) - elif output == 'eulers': table.labels_append(['{}_{}_{}({})'.format(i+1,'eulr',options.symmetry,label) for i in range(3)]) + if output in outputChoices: + table.labels_append(['{}_{}({})'.format(i+1,outputChoices[output][0],label) \ + for i in range(outputChoices[output][1])]) table.head_write() # ------------------------------------------ process data ------------------------------------------ outputAlive = True while outputAlive and table.data_read(): # read next data line of ASCII table - if inputtype == 'eulers': - o = damask.Orientation(Eulers = np.array(list(map(float,table.data[column:column+3])))*toRadians, - symmetry = options.symmetry).reduced() + if inputtype == 'eulers': + o = damask.Orientation(Eulers = np.array(list(map(float,table.data[column:column+3])))*toRadians) elif inputtype == 'rodrigues': - o = damask.Orientation(Rodrigues= np.array(list(map(float,table.data[column:column+3]))), - symmetry = options.symmetry).reduced() + o = damask.Orientation(Rodrigues= np.array(list(map(float,table.data[column:column+3])))) elif inputtype == 'matrix': - o = damask.Orientation(matrix = np.array(list(map(float,table.data[column:column+9]))).reshape(3,3).transpose(), - symmetry = options.symmetry).reduced() + o = damask.Orientation(matrix = np.array(list(map(float,table.data[column:column+9]))).reshape(3,3)) elif inputtype == 'frame': - o = damask.Orientation(matrix = np.array(list(map(float,table.data[column[0]:column[0]+3] + \ - table.data[column[1]:column[1]+3] + \ - table.data[column[2]:column[2]+3]))).reshape(3,3), - symmetry = options.symmetry).reduced() + M = np.array(list(map(float,table.data[column[0]:column[0]+3] + \ + table.data[column[1]:column[1]+3] + \ + table.data[column[2]:column[2]+3]))).reshape(3,3).T + o = damask.Orientation(matrix = M/np.linalg.norm(M,axis=0)) elif inputtype == 'quaternion': - o = damask.Orientation(quaternion = np.array(list(map(float,table.data[column:column+4]))), - symmetry = options.symmetry).reduced() + o = damask.Orientation(quaternion = np.array(list(map(float,table.data[column:column+4])))) o.quaternion = r*o.quaternion*R # apply additional lab and crystal frame rotations @@ -170,6 +168,8 @@ for name in filenames: if output == 'quaternion': table.data_append(o.asQuaternion()) elif output == 'rodrigues': table.data_append(o.asRodrigues()) elif output == 'eulers': table.data_append(o.asEulers(degrees=options.degrees)) + elif output == 'matrix': table.data_append(o.asMatrix()) + elif output == 'angleaxis': table.data_append(o.asAngleAxis(degrees=options.degrees,flat=True)) outputAlive = table.data_write() # output processed line # ------------------------------------------ output finalization ----------------------------------- From 18b791bb63ce946720d4cd5d9cbffb0beeebc096 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Fri, 7 Dec 2018 22:06:47 -0500 Subject: [PATCH 23/70] fixed bug (not clear how it escaped up to now) --- lib/damask/asciitable.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/lib/damask/asciitable.py b/lib/damask/asciitable.py index 839adb944..448587db0 100644 --- a/lib/damask/asciitable.py +++ b/lib/damask/asciitable.py @@ -493,8 +493,8 @@ class ASCIItable(): (d if str(c) != str(labels[present[i]]) else 1))) use = np.array(columns) if len(columns) > 0 else None - - self.tags = list(np.array(self.tags)[use]) # update labels with valid subset + + self.tags = list(np.array(self.__IO__['tags'])[use]) # update labels with valid subset self.data = np.loadtxt(self.__IO__['in'],usecols=use,ndmin=2) # self.data = np.genfromtxt(self.__IO__['in'],dtype=None,names=self.tags,usecols=use) From 02c2e902a7a89d3e0bffbc2dcb7ff0d86527959b Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 6 Dec 2018 22:23:34 +0100 Subject: [PATCH 24/70] texture component test also needs a closer look --- .gitlab-ci.yml | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 368ca94db..99daef5e5 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -364,12 +364,12 @@ Phenopowerlaw_singleSlip: - master - release -TextureComponents: - stage: spectral - script: TextureComponents/test.py - except: - - master - - release +#TextureComponents: +# stage: spectral +# script: TextureComponents/test.py +# except: +# - master +# - release ################################################################################################### From bb81731bc7d19e85f1e58a961cdbfe4de9fcc5ed Mon Sep 17 00:00:00 2001 From: Test User Date: Fri, 7 Dec 2018 03:32:36 +0100 Subject: [PATCH 25/70] [skip ci] updated version information after successful test of v2.0.2-1122-g2349341e --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 85ef1a02a..3ad5cdbde 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.2-1107-g17716b4f +v2.0.2-1122-g2349341e From 0b99c92fb575344269cee3f4c4b31449eb3f278f Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 8 Dec 2018 07:25:39 +0100 Subject: [PATCH 26/70] updated tests --- PRIVATE | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PRIVATE b/PRIVATE index e9f93abae..64c39e902 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit e9f93abaecafbfbf11072ae70bca213a7201ed38 +Subproject commit 64c39e90200a8dd90dff3665e27050d9db81ac0d From b777253dd64997a0f97717636fb221ff47fa5d83 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Sat, 8 Dec 2018 11:52:55 -0500 Subject: [PATCH 27/70] quat->Eulers now returns always positive angles --- lib/damask/orientation.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/lib/damask/orientation.py b/lib/damask/orientation.py index 6fd0993ea..7131d39e3 100644 --- a/lib/damask/orientation.py +++ b/lib/damask/orientation.py @@ -263,8 +263,8 @@ class Quaternion: self.normalize() s = math.sqrt(1. - self.q**2) - x = 2*self.q**2 - 1. - y = 2*self.q * s + x = 2.*self.q**2 - 1. + y = 2.*self.q * s angle = math.atan2(y,x) if angle < 0.0: @@ -300,6 +300,7 @@ class Quaternion: math.atan2((P*self.q*self.p[1]+self.p[0]*self.p[2])/chi,( self.p[1]*self.p[2]-P*self.q*self.p[0])/chi), ]) + eulers %= 2.0*math.pi return np.degrees(eulers) if degrees else eulers From 1de37b9397b3c294beecc474d32694395ccff0b9 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Sat, 8 Dec 2018 12:01:46 -0500 Subject: [PATCH 28/70] polished formatting of comments --- lib/damask/orientation.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/lib/damask/orientation.py b/lib/damask/orientation.py index 7131d39e3..db2baad92 100644 --- a/lib/damask/orientation.py +++ b/lib/damask/orientation.py @@ -87,11 +87,11 @@ class Quaternion: """Multiplication""" # Rowenhorst_etal2015 MSMSE: value of P is selected as -1 P = -1.0 - try: # quaternion + try: # 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)) except: pass - try: # vector (perform passive rotation) + try: # vector (perform passive rotation) ( x, y, z) = self.p (Vx,Vy,Vz) = other[0:3] A = self.q*self.q - np.dot(self.p,self.p) @@ -104,7 +104,7 @@ class Quaternion: A*Vz + B*z + C*(x*Vy - y*Vx), ]) except: pass - try: # scalar + try: # scalar return self.__class__(q=self.q*other, p=self.p*other) except: @@ -114,7 +114,7 @@ class Quaternion: """In-place multiplication""" # Rowenhorst_etal2015 MSMSE: value of P is selected as -1 P = -1.0 - try: # Quaternion + try: # 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) except: pass @@ -237,7 +237,7 @@ class Quaternion: def asList(self): return [self.q]+list(self.p) - def asM(self): # to find Averaging Quaternions (see F. Landis Markley et al.) + def asM(self): # to find Averaging Quaternions (see F. Landis Markley et al.) return np.outer(self.asList(),self.asList()) def asMatrix(self): @@ -300,7 +300,7 @@ class Quaternion: math.atan2((P*self.q*self.p[1]+self.p[0]*self.p[2])/chi,( self.p[1]*self.p[2]-P*self.q*self.p[0])/chi), ]) - eulers %= 2.0*math.pi + eulers %= 2.0*math.pi # enforce positive angles return np.degrees(eulers) if degrees else eulers From 0ff9e9fa067ab8034634153b2722f1d5c5e91387 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 8 Dec 2018 18:13:16 +0100 Subject: [PATCH 29/70] report CMake should be updated on the homepage also --- DAMASK_prerequisites.sh | 3 +++ 1 file changed, 3 insertions(+) diff --git a/DAMASK_prerequisites.sh b/DAMASK_prerequisites.sh index c9e165c3c..f0a525c34 100755 --- a/DAMASK_prerequisites.sh +++ b/DAMASK_prerequisites.sh @@ -119,6 +119,9 @@ for executable in mpirun mpiexec; do getDetails $executable '--version' done +firstLevel "CMake" +getDetails cmake --version + firstLevel "Abaqus" cd installation/mods_Abaqus # to have the right environment file for executable in abaqus abq2017 abq2018; do From 3f8726464d20e3bdc37961a08a58309524ef4e70 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Sat, 8 Dec 2018 13:22:22 -0500 Subject: [PATCH 30/70] addOrientations now checks its input for convention conformity --- processing/post/addOrientations.py | 37 ++++++++++++++++++++++++++---- 1 file changed, 32 insertions(+), 5 deletions(-) diff --git a/processing/post/addOrientations.py b/processing/post/addOrientations.py index b4a8ed07c..efd4dac24 100755 --- a/processing/post/addOrientations.py +++ b/processing/post/addOrientations.py @@ -9,6 +9,29 @@ import damask scriptName = os.path.splitext(os.path.basename(__file__))[0] scriptID = ' '.join([scriptName,damask.version]) +# -------------------------------------------------------------------- +# convention conformity checks +# -------------------------------------------------------------------- + +def check_Eulers(eulers): + if np.any(eulers < 0.0) or np.any(eulers > 2.0*math.pi) or eulers[1] > math.pi: # Euler angles within valid range? + raise ValueError('Euler angles outside of [0..2π],[0..π],[0..2π].\n{} {} {}.'.format(*eulers)) + return eulers + +def check_quaternion(q): + if q.q < 0.0: # positive first quaternion component? + raise ValueError('quaternion has negative first component.\n{}'.format(q)) + return q + +def check_matrix(M): + if abs(1.0-np.linalg.det(M)) > 1e-8: # proper rotation? + raise ValueError('matrix is not a proper rotation.\n{}'.format(M)) + if abs(np.dot(M[0],M[1])) > 1e-8 \ + or abs(np.dot(M[1],M[2])) > 1e-8 \ + or abs(np.dot(M[2],M[0])) > 1e-8: # all orthogonal? + raise ValueError('matrix is not orthogonal.\n{}'.format(M)) + return M + # -------------------------------------------------------------------- # MAIN # -------------------------------------------------------------------- @@ -103,6 +126,7 @@ if np.sum(input) != 1: parser.error('needs exactly one input format.') (options.matrix,9,'matrix'), (options.quaternion,4,'quaternion'), ][np.where(input)[0][0]] # select input label that was requested + toRadians = math.pi/180.0 if options.degrees else 1.0 # rescale degrees to radians r = damask.Quaternion.fromAngleAxis(toRadians*options.crystalrotation[0],options.crystalrotation[1:]) # crystal frame rotation R = damask.Quaternion.fromAngleAxis(toRadians*options. labrotation[0],options. labrotation[1:]) # lab frame rotation @@ -149,18 +173,21 @@ for name in filenames: outputAlive = True while outputAlive and table.data_read(): # read next data line of ASCII table if inputtype == 'eulers': - o = damask.Orientation(Eulers = np.array(list(map(float,table.data[column:column+3])))*toRadians) + + o = damask.Orientation(Eulers = check_Eulers(np.array(list(map(float,table.data[column:column+3])))*toRadians)) elif inputtype == 'rodrigues': - o = damask.Orientation(Rodrigues= np.array(list(map(float,table.data[column:column+3])))) + o = damask.Orientation(Rodrigues = np.array(list(map(float,table.data[column:column+3])))) elif inputtype == 'matrix': - o = damask.Orientation(matrix = np.array(list(map(float,table.data[column:column+9]))).reshape(3,3)) + + o = damask.Orientation(matrix = check_matrix(np.array(list(map(float,table.data[column:column+9]))).reshape(3,3))) elif inputtype == 'frame': M = np.array(list(map(float,table.data[column[0]:column[0]+3] + \ table.data[column[1]:column[1]+3] + \ table.data[column[2]:column[2]+3]))).reshape(3,3).T - o = damask.Orientation(matrix = M/np.linalg.norm(M,axis=0)) + o = damask.Orientation(matrix = check_matrix(M/np.linalg.norm(M,axis=0))) elif inputtype == 'quaternion': - o = damask.Orientation(quaternion = np.array(list(map(float,table.data[column:column+4])))) + + o = damask.Orientation(quaternion = check_quaternion(np.array(list(map(float,table.data[column:column+4]))))) o.quaternion = r*o.quaternion*R # apply additional lab and crystal frame rotations From 4130cbcffa7705c360530a1b4ec7ecc3bc5e3637 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 9 Dec 2018 07:10:31 +0100 Subject: [PATCH 31/70] simplified - always use quaternion as input for orientation - addAPS34IDEstrainCoords uses table functionality - rotateData can figure out whether input is vector or tensor --- processing/post/addAPS34IDEstrainCoords.py | 30 +++++---- processing/post/addGrainID.py | 76 +++------------------- processing/post/addIPFcolor.py | 64 ++---------------- processing/post/addPole.py | 70 +++----------------- processing/post/rotateData.py | 65 ++++++++---------- 5 files changed, 67 insertions(+), 238 deletions(-) diff --git a/processing/post/addAPS34IDEstrainCoords.py b/processing/post/addAPS34IDEstrainCoords.py index 2e753bfa4..0b54637f6 100755 --- a/processing/post/addAPS34IDEstrainCoords.py +++ b/processing/post/addAPS34IDEstrainCoords.py @@ -35,21 +35,27 @@ datainfo = {'len':3, datainfo['label'] += options.frame -# --- loop over input files ------------------------------------------------------------------------- -if filenames == []: - filenames = ['STDIN'] +# --- loop over input files ------------------------------------------------------------------------ + +if filenames == []: filenames = [None] for name in filenames: - if name == 'STDIN': - file = {'name':'STDIN', 'input':sys.stdin, 'output':sys.stdout, 'croak':sys.stderr} - file['croak'].write('\033[1m'+scriptName+'\033[0m\n') - else: - if not os.path.exists(name): continue - file = {'name':name, 'input':open(name), 'output':open(name+'_tmp','w'), 'croak':sys.stderr} - file['croak'].write('\033[1m'+scriptName+'\033[0m: '+file['name']+'\n') + try: table = damask.ASCIItable(name = name, + buffered = False) + except: continue + damask.util.report(scriptName,name) + +# ------------------------------------------ read header ------------------------------------------ + + table.head_read() + + if not table.label_dimension(options.quaternion) == 4: + damask.util.croak('input {} does not have dimension 4.'.format(options.quaternion)) + table.close(dismiss = True) # close ASCIItable and remove empty file + continue + + - table = damask.ASCIItable(file['input'],file['output'],buffered=False) # make unbuffered ASCII_table - table.head_read() # read ASCII header info table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:])) # --------------- figure out columns to process --------------------------------------------------- diff --git a/processing/post/addGrainID.py b/processing/post/addGrainID.py index c3b98f4e6..dedf5b278 100755 --- a/processing/post/addGrainID.py +++ b/processing/post/addGrainID.py @@ -31,32 +31,6 @@ parser.add_option('-s', dest = 'symmetry', type = 'string', metavar = 'string', help = 'crystal symmetry [%default]') -parser.add_option('-e', - '--eulers', - dest = 'eulers', - type = 'string', metavar = 'string', - help = 'label of Euler angles') -parser.add_option('--degrees', - dest = 'degrees', - action = 'store_true', - help = 'Euler angles are given in degrees [%default]') -parser.add_option('-m', - '--matrix', - dest = 'matrix', - type = 'string', metavar = 'string', - help = 'label of orientation matrix') -parser.add_option('-a', - dest = 'a', - type = 'string', metavar = 'string', - help = 'label of crystal frame a vector') -parser.add_option('-b', - dest = 'b', - type = 'string', metavar = 'string', - help = 'label of crystal frame b vector') -parser.add_option('-c', - dest = 'c', - type = 'string', metavar = 'string', - help = 'label of crystal frame c vector') parser.add_option('-q', '--quaternion', dest = 'quaternion', @@ -69,9 +43,9 @@ parser.add_option('-p', help = 'label of coordinates [%default]') parser.set_defaults(disorientation = 5, + quaternion = 'orientation', symmetry = 'cubic', pos = 'pos', - degrees = False, ) (options, filenames) = parser.parse_args() @@ -79,22 +53,6 @@ parser.set_defaults(disorientation = 5, if options.radius is None: parser.error('no radius specified.') -input = [options.eulers is not None, - options.a is not None and \ - options.b is not None and \ - options.c is not None, - options.matrix is not None, - options.quaternion is not None, - ] - -if np.sum(input) != 1: parser.error('needs exactly one input format.') - -(label,dim,inputtype) = [(options.eulers,3,'eulers'), - ([options.a,options.b,options.c],[3,3,3],'frame'), - (options.matrix,9,'matrix'), - (options.quaternion,4,'quaternion'), - ][np.where(input)[0][0]] # select input label that was requested -toRadians = np.pi/180.0 if options.degrees else 1.0 # rescale degrees to radians cos_disorientation = np.cos(np.radians(options.disorientation/2.)) # cos of half the disorientation angle # --- loop over input files ------------------------------------------------------------------------- @@ -118,8 +76,8 @@ for name in filenames: if not 3 >= table.label_dimension(options.pos) >= 1: errors.append('coordinates "{}" need to have one, two, or three dimensions.'.format(options.pos)) - if not np.all(table.label_dimension(label) == dim): - errors.append('input "{}" does not have dimension {}.'.format(label,dim)) + if not np.all(table.label_dimension(options.quaternion) == 4): + errors.append('input "{}" does not have dimension 4.'.format(options.quaternion)) else: column = table.label_index(label) if remarks != []: damask.util.croak(remarks) @@ -131,9 +89,7 @@ for name in filenames: # ------------------------------------------ assemble header --------------------------------------- table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:])) - table.labels_append('grainID_{}@{:g}'.format('+'.join(label) - if isinstance(label, (list,tuple)) - else label, + table.labels_append('grainID_{}@{:g}'.format('+'.join(options.quaternion), options.disorientation)) # report orientation source and disorientation table.head_write() @@ -141,16 +97,14 @@ for name in filenames: # ------------------------------------------ build KD tree ----------------------------------------- + table.data_readArray(options.pos) # read position vectors + grainID = -np.ones(len(table.data),dtype=int) + # --- start background messaging bg = damask.util.backgroundMessage() bg.start() - bg.set_message('reading positions...') - - table.data_readArray(options.pos) # read position vectors - grainID = -np.ones(len(table.data),dtype=int) - start = tick = time.clock() bg.set_message('building KD tree...') kdtree = spatial.KDTree(copy.deepcopy(table.data)) @@ -175,20 +129,8 @@ for name in filenames: bg.set_message('(%02i:%02i:%02i) processing point %i of %i (grain count %i)...'\ %(time_delta//3600,time_delta%3600//60,time_delta%60,p,len(grainID),np.count_nonzero(memberCounts))) - if inputtype == 'eulers': - o = damask.Orientation(Eulers = np.array(map(float,table.data[column:column+3]))*toRadians, - symmetry = options.symmetry).reduced() - elif inputtype == 'matrix': - o = damask.Orientation(matrix = np.array(map(float,table.data[column:column+9])).reshape(3,3).transpose(), - symmetry = options.symmetry).reduced() - elif inputtype == 'frame': - o = damask.Orientation(matrix = np.array(map(float,table.data[column[0]:column[0]+3] + \ - table.data[column[1]:column[1]+3] + \ - table.data[column[2]:column[2]+3])).reshape(3,3), - symmetry = options.symmetry).reduced() - elif inputtype == 'quaternion': - o = damask.Orientation(quaternion = np.array(map(float,table.data[column:column+4])), - symmetry = options.symmetry).reduced() + o = damask.Orientation(quaternion = np.array(map(float,table.data[column:column+4])), + symmetry = options.symmetry).reduced() matched = False alreadyChecked = {} diff --git a/processing/post/addIPFcolor.py b/processing/post/addIPFcolor.py index ffe3dac3f..1e5d8c2f3 100755 --- a/processing/post/addIPFcolor.py +++ b/processing/post/addIPFcolor.py @@ -26,58 +26,18 @@ parser.add_option('-s', '--symmetry', dest = 'symmetry', type = 'choice', choices = damask.Symmetry.lattices[1:], metavar='string', help = 'crystal symmetry [%default] {{{}}} '.format(', '.join(damask.Symmetry.lattices[1:]))) -parser.add_option('-e', '--eulers', - dest = 'eulers', - type = 'string', metavar = 'string', - help = 'Euler angles label') -parser.add_option('-d', '--degrees', - dest = 'degrees', - action = 'store_true', - help = 'Euler angles are given in degrees [%default]') -parser.add_option('-m', '--matrix', - dest = 'matrix', - type = 'string', metavar = 'string', - help = 'orientation matrix label') -parser.add_option('-a', - dest = 'a', - type = 'string', metavar = 'string', - help = 'crystal frame a vector label') -parser.add_option('-b', - dest = 'b', - type = 'string', metavar = 'string', - help = 'crystal frame b vector label') -parser.add_option('-c', - dest = 'c', - type = 'string', metavar = 'string', - help = 'crystal frame c vector label') parser.add_option('-q', '--quaternion', dest = 'quaternion', type = 'string', metavar = 'string', help = 'quaternion label') parser.set_defaults(pole = (0.0,0.0,1.0), + quaternion = 'orientation', symmetry = damask.Symmetry.lattices[-1], - degrees = False, ) (options, filenames) = parser.parse_args() -input = [options.eulers is not None, - options.a is not None and \ - options.b is not None and \ - options.c is not None, - options.matrix is not None, - options.quaternion is not None, - ] - -if np.sum(input) != 1: parser.error('needs exactly one input format.') - -(label,dim,inputtype) = [(options.eulers,3,'eulers'), - ([options.a,options.b,options.c],[3,3,3],'frame'), - (options.matrix,9,'matrix'), - (options.quaternion,4,'quaternion'), - ][np.where(input)[0][0]] # select input label that was requested -toRadians = math.pi/180.0 if options.degrees else 1.0 # rescale degrees to radians pole = np.array(options.pole) pole /= np.linalg.norm(pole) @@ -98,12 +58,12 @@ for name in filenames: # ------------------------------------------ sanity checks ---------------------------------------- - if not np.all(table.label_dimension(label) == dim): - damask.util.croak('input {} does not have dimension {}.'.format(label,dim)) + if not table.label_dimension(options.quaternion) == 4: + damask.util.croak('input {} does not have dimension 4.'.format(options.quaternion)) table.close(dismiss = True) # close ASCIItable and remove empty file continue - column = table.label_index(label) + column = table.label_index(options.quaternion) # ------------------------------------------ assemble header --------------------------------------- @@ -115,20 +75,8 @@ for name in filenames: outputAlive = True while outputAlive and table.data_read(): # read next data line of ASCII table - if inputtype == 'eulers': - o = damask.Orientation(Eulers = np.array(list(map(float,table.data[column:column+3])))*toRadians, - symmetry = options.symmetry).reduced() - elif inputtype == 'matrix': - o = damask.Orientation(matrix = np.array(list(map(float,table.data[column:column+9]))).reshape(3,3).transpose(), - symmetry = options.symmetry).reduced() - elif inputtype == 'frame': - o = damask.Orientation(matrix = np.array(list(map(float,table.data[column[0]:column[0]+3] + \ - table.data[column[1]:column[1]+3] + \ - table.data[column[2]:column[2]+3]))).reshape(3,3), - symmetry = options.symmetry).reduced() - elif inputtype == 'quaternion': - o = damask.Orientation(quaternion = np.array(list(map(float,table.data[column:column+4]))), - symmetry = options.symmetry).reduced() + o = damask.Orientation(quaternion = np.array(list(map(float,table.data[column:column+4]))), + symmetry = options.symmetry).reduced() table.data_append(o.IPFcolor(pole)) outputAlive = table.data_write() # output processed line diff --git a/processing/post/addPole.py b/processing/post/addPole.py index 95bc87637..83a031d46 100755 --- a/processing/post/addPole.py +++ b/processing/post/addPole.py @@ -1,4 +1,4 @@ -#!/usr/bin/env python2.7 +#!/usr/bin/env python3 # -*- coding: UTF-8 no BOM -*- import os,sys,math @@ -26,58 +26,18 @@ parser.add_option('--polar', dest = 'polar', action = 'store_true', help = 'output polar coordinates r,phi [%default]') -parser.add_option('-e', '--eulers', - dest = 'eulers', - type = 'string', metavar = 'string', - help = 'Euler angles label') -parser.add_option('-d', '--degrees', - dest = 'degrees', - action = 'store_true', - help = 'Euler angles are given in degrees [%default]') -parser.add_option('-m', '--matrix', - dest = 'matrix', - type = 'string', metavar = 'string', - help = 'orientation matrix label') -parser.add_option('-a', - dest = 'a', - type = 'string', metavar = 'string', - help = 'crystal frame a vector label') -parser.add_option('-b', - dest = 'b', - type = 'string', metavar = 'string', - help = 'crystal frame b vector label') -parser.add_option('-c', - dest = 'c', - type = 'string', metavar = 'string', - help = 'crystal frame c vector label') parser.add_option('-q', '--quaternion', dest = 'quaternion', type = 'string', metavar = 'string', help = 'quaternion label') parser.set_defaults(pole = (1.0,0.0,0.0), - degrees = False, + quaternion = 'orientation', polar = False, ) (options, filenames) = parser.parse_args() -input = [options.eulers is not None, - options.a is not None and \ - options.b is not None and \ - options.c is not None, - options.matrix is not None, - options.quaternion is not None, - ] - -if np.sum(input) != 1: parser.error('needs exactly one input format.') - -(label,dim,inputtype) = [(options.eulers,3,'eulers'), - ([options.a,options.b,options.c],[3,3,3],'frame'), - (options.matrix,9,'matrix'), - (options.quaternion,4,'quaternion'), - ][np.where(input)[0][0]] # select input label that was requested -toRadians = math.pi/180.0 if options.degrees else 1.0 # rescale degrees to radians pole = np.array(options.pole) pole /= np.linalg.norm(pole) @@ -98,18 +58,13 @@ for name in filenames: # ------------------------------------------ sanity checks ---------------------------------------- - errors = [] - remarks = [] - - if not np.all(table.label_dimension(label) == dim): errors.append('input {} does not have dimension {}.'.format(label,dim)) - else: column = table.label_index(label) - - if remarks != []: damask.util.croak(remarks) - if errors != []: - damask.util.croak(errors) - table.close(dismiss = True) + if not table.label_dimension(options.quaternion) == 4: + damask.util.croak('input {} does not have dimension 4.'.format(options.quaternion)) + table.close(dismiss = True) # close ASCIItable and remove empty file continue + column = table.label_index(options.quaternion) + # ------------------------------------------ assemble header --------------------------------------- table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:])) @@ -119,16 +74,7 @@ for name in filenames: # ------------------------------------------ process data ------------------------------------------ outputAlive = True while outputAlive and table.data_read(): # read next data line of ASCII table - if inputtype == 'eulers': - o = damask.Orientation(Eulers = np.array(list(map(float,table.data[column:column+3])))*toRadians) - elif inputtype == 'matrix': - o = damask.Orientation(matrix = np.array(list(map(float,table.data[column:column+9]))).reshape(3,3).transpose()) - elif inputtype == 'frame': - o = damask.Orientation(matrix = np.array(list(map(float,table.data[column[0]:column[0]+3] + \ - table.data[column[1]:column[1]+3] + \ - table.data[column[2]:column[2]+3]))).reshape(3,3)) - elif inputtype == 'quaternion': - o = damask.Orientation(quaternion = np.array(list(map(float,table.data[column:column+4])))) + o = damask.Orientation(quaternion = np.array(list(map(float,table.data[column:column+4])))) rotatedPole = o.quaternion*pole # rotate pole according to crystal orientation (x,y) = rotatedPole[0:2]/(1.+abs(pole[2])) # stereographic projection diff --git a/processing/post/rotateData.py b/processing/post/rotateData.py index 84d4cd3e1..683fbf863 100755 --- a/processing/post/rotateData.py +++ b/processing/post/rotateData.py @@ -18,19 +18,15 @@ Rotate vector and/or tensor column data by given angle around given axis. """, version = scriptID) -parser.add_option('-v','--vector', - dest = 'vector', +parser.add_option('-d', '--data', + dest = 'data', action = 'extend', metavar = '', - help = 'column heading of vector(s) to rotate') -parser.add_option('-t','--tensor', - dest = 'tensor', - action = 'extend', metavar = '', - help = 'column heading of tensor(s) to rotate') + help = 'vector/tensor value(s) label(s)') parser.add_option('-r', '--rotation', dest = 'rotation', type = 'float', nargs = 4, metavar = ' '.join(['float']*4), help = 'angle and axis to rotate data [%default]') -parser.add_option('-d', '--degrees', +parser.add_option('--degrees', dest = 'degrees', action = 'store_true', help = 'angles are given in degrees [%default]') @@ -41,7 +37,7 @@ parser.set_defaults(rotation = (0.,1.,1.,1.), (options,filenames) = parser.parse_args() -if options.vector is None and options.tensor is None: +if options.data is None: parser.error('no data column specified.') toRadians = math.pi/180.0 if options.degrees else 1.0 # rescale degrees to radians @@ -59,27 +55,24 @@ for name in filenames: except: continue damask.util.report(scriptName,name) -# ------------------------------------------ read header ------------------------------------------ +# --- interpret header ---------------------------------------------------------------------------- table.head_read() -# ------------------------------------------ sanity checks ---------------------------------------- - - items = { - 'tensor': {'dim': 9, 'shape': [3,3], 'labels':options.tensor, 'active':[], 'column': []}, - 'vector': {'dim': 3, 'shape': [3], 'labels':options.vector, 'active':[], 'column': []}, - } errors = [] remarks = [] - column = {} - - for type, data in items.items(): - for what in data['labels']: - dim = table.label_dimension(what) - if dim != data['dim']: remarks.append('column {} is not a {}.'.format(what,type)) - else: - items[type]['active'].append(what) - items[type]['column'].append(table.label_index(what)) + active = {'vector':[],'tensor':[]} + + for i,dim in enumerate(table.label_dimension(options.data)): + label = options.data[i] + if dim == -1: + remarks.append('{} "{}" not found...'.format(datatype,me)) + elif dim == 3: + remarks.append('adding vector "{}"...'.format(label)) + active['vector'].append(label) + elif dim == 9: + remarks.append('adding tensor "{}"...'.format(label)) + active['tensor'].append(label) if remarks != []: damask.util.croak(remarks) if errors != []: @@ -95,20 +88,14 @@ for name in filenames: # ------------------------------------------ process data ------------------------------------------ outputAlive = True while outputAlive and table.data_read(): # read next data line of ASCII table - - datatype = 'vector' - - for column in items[datatype]['column']: # loop over all requested labels - table.data[column:column+items[datatype]['dim']] = \ - q * np.array(list(map(float,table.data[column:column+items[datatype]['dim']]))) - - datatype = 'tensor' - - for column in items[datatype]['column']: # loop over all requested labels - table.data[column:column+items[datatype]['dim']] = \ - np.dot(R,np.dot(np.array(list(map(float,table.data[column:column+items[datatype]['dim']]))).\ - reshape(items[datatype]['shape']),R.transpose())).reshape(items[datatype]['dim']) - + for v in active['vector']: + column = table.label_index(v) + table.data[column:column+3] = q * np.array(list(map(float,table.data[column:column+3]))) + for t in active['tensor']: + column = table.label_index(v) + table.data[column:column+9] = \ + np.dot(R,np.dot(np.array(list(map(float,table.data[column:column+9]]))).\ + reshape((3,3)),R.transpose())).reshape((9)) outputAlive = table.data_write() # output processed line # ------------------------------------------ output finalization ----------------------------------- From d145b257aaf433f255e9b739130be894af039fb3 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 9 Dec 2018 07:57:05 +0100 Subject: [PATCH 32/70] new version of numpy complain about overlong range argument --- processing/post/binXY.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/processing/post/binXY.py b/processing/post/binXY.py index aac9aaffa..726fca39f 100755 --- a/processing/post/binXY.py +++ b/processing/post/binXY.py @@ -1,4 +1,4 @@ -#!/usr/bin/env python2.7 +#!/usr/bin/env python3 # -*- coding: UTF-8 no BOM -*- import os,sys @@ -118,10 +118,9 @@ for name in filenames: minmax[c] = np.log(minmax[c]) # change minmax to log, too delta = minmax[:,1]-minmax[:,0] - (grid,xedges,yedges) = np.histogram2d(table.data[:,0],table.data[:,1], bins=options.bins, - range=minmax, + range=minmax[0:2,0:2], weights=None if options.weight is None else table.data[:,2]) if options.normCol: From 83a3628765d4a876610f2e3afb73f921106faefe Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 9 Dec 2018 07:58:42 +0100 Subject: [PATCH 33/70] bug fixes --- processing/post/rotateData.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/processing/post/rotateData.py b/processing/post/rotateData.py index 683fbf863..7ed9542cf 100755 --- a/processing/post/rotateData.py +++ b/processing/post/rotateData.py @@ -92,10 +92,10 @@ for name in filenames: column = table.label_index(v) table.data[column:column+3] = q * np.array(list(map(float,table.data[column:column+3]))) for t in active['tensor']: - column = table.label_index(v) + column = table.label_index(t) table.data[column:column+9] = \ - np.dot(R,np.dot(np.array(list(map(float,table.data[column:column+9]]))).\ - reshape((3,3)),R.transpose())).reshape((9)) + np.dot(R,np.dot(np.array(list(map(float,table.data[column:column+9]))).reshape((3,3)), + R.transpose())).reshape((9)) outputAlive = table.data_write() # output processed line # ------------------------------------------ output finalization ----------------------------------- From 6b73840853b1418d9176de22a8fd14d27a9d1b2b Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 9 Dec 2018 08:07:40 +0100 Subject: [PATCH 34/70] works with python3 --- processing/post/addInfo.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/processing/post/addInfo.py b/processing/post/addInfo.py index 50f003d6b..59efcd973 100755 --- a/processing/post/addInfo.py +++ b/processing/post/addInfo.py @@ -1,4 +1,4 @@ -#!/usr/bin/env python2.7 +#!/usr/bin/env python3 # -*- coding: UTF-8 no BOM -*- import os From 7eaea64d091ba06add88bee16ed5d759edeb5286 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 9 Dec 2018 09:08:33 +0100 Subject: [PATCH 35/70] python3 compatible hope that the new bar causes less trouble than the background message --- lib/damask/util.py | 24 ++++++++++++++++++++++++ processing/post/addGrainID.py | 28 ++++++---------------------- 2 files changed, 30 insertions(+), 22 deletions(-) diff --git a/lib/damask/util.py b/lib/damask/util.py index 8727a1473..39715730c 100644 --- a/lib/damask/util.py +++ b/lib/damask/util.py @@ -132,6 +132,30 @@ class extendableOption(Option): else: Option.take_action(self, action, dest, opt, value, values, parser) +# Print iterations progress +# from https://gist.github.com/aubricus/f91fb55dc6ba5557fbab06119420dd6a +def print_progress(iteration, total, prefix='', suffix='', decimals=1, bar_length=100): + """ + Call in a loop to create terminal progress bar + @params: + iteration - Required : current iteration (Int) + total - Required : total iterations (Int) + prefix - Optional : prefix string (Str) + suffix - Optional : suffix string (Str) + decimals - Optional : positive number of decimals in percent complete (Int) + bar_length - Optional : character length of bar (Int) + """ + str_format = "{0:." + str(decimals) + "f}" + percents = str_format.format(100 * (iteration / float(total))) + filled_length = int(round(bar_length * iteration / float(total))) + bar = '█' * filled_length + '-' * (bar_length - filled_length) + + sys.stdout.write('\r%s |%s| %s%s %s' % (prefix, bar, percents, '%', suffix)), + + if iteration == total: + sys.stdout.write('\n') +sys.stdout.flush() + # ----------------------------- class backgroundMessage(threading.Thread): """Reporting with animation to indicate progress""" diff --git a/processing/post/addGrainID.py b/processing/post/addGrainID.py index dedf5b278..84bdaa30d 100755 --- a/processing/post/addGrainID.py +++ b/processing/post/addGrainID.py @@ -1,4 +1,4 @@ -#!/usr/bin/env python2.7 +#!/usr/bin/env python3 # -*- coding: UTF-8 no BOM -*- import os,sys,time,copy @@ -78,7 +78,7 @@ for name in filenames: errors.append('coordinates "{}" need to have one, two, or three dimensions.'.format(options.pos)) if not np.all(table.label_dimension(options.quaternion) == 4): errors.append('input "{}" does not have dimension 4.'.format(options.quaternion)) - else: column = table.label_index(label) + else: column = table.label_index(options.quaternion) if remarks != []: damask.util.croak(remarks) if errors != []: @@ -93,26 +93,15 @@ for name in filenames: options.disorientation)) # report orientation source and disorientation table.head_write() -# ------------------------------------------ process data ------------------------------------------ - # ------------------------------------------ build KD tree ----------------------------------------- table.data_readArray(options.pos) # read position vectors grainID = -np.ones(len(table.data),dtype=int) - -# --- start background messaging - - bg = damask.util.backgroundMessage() - bg.start() - - start = tick = time.clock() - bg.set_message('building KD tree...') + Npoints = table.data.shape[0] kdtree = spatial.KDTree(copy.deepcopy(table.data)) # ------------------------------------------ assign grain IDs -------------------------------------- - tick = time.clock() - orientations = [] # quaternions found for grain memberCounts = [] # number of voxels in grain p = 0 # point counter @@ -123,13 +112,10 @@ for name in filenames: table.data_rewind() while table.data_read(): # read next data line of ASCII table - if p > 0 and p % 1000 == 0: + if Npoints > 100 and p%(Npoints//100) == 0: # report in 1% steps if possible and avoid modulo by zero + damask.util.print_progress(iteration=p,total=Npoints) - time_delta = (time.clock()-tick) * (len(grainID) - p) / p - bg.set_message('(%02i:%02i:%02i) processing point %i of %i (grain count %i)...'\ - %(time_delta//3600,time_delta%3600//60,time_delta%60,p,len(grainID),np.count_nonzero(memberCounts))) - - o = damask.Orientation(quaternion = np.array(map(float,table.data[column:column+4])), + o = damask.Orientation(quaternion = np.array(list(map(float,table.data[column:column+4]))), symmetry = options.symmetry).reduced() matched = False @@ -180,8 +166,6 @@ for name in filenames: outputAlive = table.data_write() # output processed line p += 1 - bg.set_message('done after {} seconds'.format(time.clock()-start)) - # ------------------------------------------ output finalization ----------------------------------- table.close() # close ASCII tables From d9b47f09bc3f5a311194d793bcf4e0eb36c36e9c Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 9 Dec 2018 10:52:37 +0100 Subject: [PATCH 36/70] modernized split arguments logically, assume 3-vector for coordinates --- processing/post/addAPS34IDEstrainCoords.py | 71 ++++++++-------------- 1 file changed, 27 insertions(+), 44 deletions(-) diff --git a/processing/post/addAPS34IDEstrainCoords.py b/processing/post/addAPS34IDEstrainCoords.py index 0b54637f6..1071baa91 100755 --- a/processing/post/addAPS34IDEstrainCoords.py +++ b/processing/post/addAPS34IDEstrainCoords.py @@ -1,4 +1,4 @@ -#!/usr/bin/env python2.7 +#!/usr/bin/env python3 # -*- coding: UTF-8 no BOM -*- import os,sys @@ -19,21 +19,22 @@ Transform X,Y,Z,F APS BeamLine 34 coordinates to x,y,z APS strain coordinates. """, version = scriptID) -parser.add_option('-f','--frame', dest='frame', nargs=4, type='string', metavar='string string string string', - help='APS X,Y,Z coords, and depth F') -parser.set_defaults(frame = None) +parser.add_option('-f', + '--frame', + dest='frame', + metavar='string', + help='APS X,Y,Z coords') +parser.add_option('--depth', + dest='depth', + metavar='string', + help='depth') (options,filenames) = parser.parse_args() if options.frame is None: - parser.error('no data column specified...') - - -datainfo = {'len':3, - 'label':[] - } - -datainfo['label'] += options.frame + parser.error('frame not specified') +if options.depth is None: + parser.error('depth not specified') # --- loop over input files ------------------------------------------------------------------------ @@ -49,31 +50,19 @@ for name in filenames: table.head_read() - if not table.label_dimension(options.quaternion) == 4: - damask.util.croak('input {} does not have dimension 4.'.format(options.quaternion)) - table.close(dismiss = True) # close ASCIItable and remove empty file +# ------------------------------------------ sanity checks ----------------------------------------- + errors = [] + if table.label_dimension(options.frame) != 3: + errors.append('input {} does not have dimension 3.'.format(options.frame)) + if table.label_dimension(options.depth) != 1: + errors.append('input {} does not have dimension 1.'.format(options.depth)) + if errors != []: + damask.util.croak(errors) + table.close(dismiss = True) continue - - table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:])) -# --------------- figure out columns to process --------------------------------------------------- - active = [] - column = {} - columnMissing = False - - for label in datainfo['label']: - key = label - if key in table.labels(raw = True): - active.append(label) - column[label] = table.labels.index(key) # remember columns of requested data - else: - file['croak'].write('column %s not found...\n'%label) - columnMissing = True - - if columnMissing: continue - # ------------------------------------------ assemble header --------------------------------------- table.labels_append(['%i_coord'%(i+1) for i in range(3)]) # extend ASCII header with new labels table.head_write() @@ -83,21 +72,15 @@ for name in filenames: RotMat2TSL=np.array([[1., 0., 0.], [0., np.cos(theta), np.sin(theta)], # Orientation to account for -135 deg [0., -np.sin(theta), np.cos(theta)]]) # rotation for TSL convention - vec = np.zeros(4) - outputAlive = True while outputAlive and table.data_read(): # read next data line of ASCII table - for i,label in enumerate(active): - vec[i] = table.data[column[label]] + coord = list(map(float,table.data[table.label_index(options.frame):table.label_index(options.frame)+3])) + depth = float(table.data[table.label_index(options.depth)]) - table.data_append(np.dot(RotMat2TSL,np.array([-vec[0], -vec[1],-vec[2]+vec[3]]))) + table.data_append(np.dot(RotMat2TSL,np.array([-coord[0],-coord[1],-coord[2]+depth]))) outputAlive = table.data_write() # output processed line -# ------------------------------------------ output result ----------------------------------------- - outputAlive and table.output_flush() # just in case of buffered ASCII table +# ------------------------------------------ output finalization ----------------------------------- - table.input_close() # close input ASCII table (works for stdin) - table.output_close() # close output ASCII table (works for stdout) - if file['name'] != 'STDIN': - os.rename(file['name']+'_tmp',file['name']) # overwrite old one with tmp new + table.close() # close ASCII tables From 213e31ff8700279df8264f8d135b9d030c63cba8 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 9 Dec 2018 10:53:25 +0100 Subject: [PATCH 37/70] string was split into letters --- processing/post/addGrainID.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/processing/post/addGrainID.py b/processing/post/addGrainID.py index 84bdaa30d..ac2c53a02 100755 --- a/processing/post/addGrainID.py +++ b/processing/post/addGrainID.py @@ -89,8 +89,7 @@ for name in filenames: # ------------------------------------------ assemble header --------------------------------------- table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:])) - table.labels_append('grainID_{}@{:g}'.format('+'.join(options.quaternion), - options.disorientation)) # report orientation source and disorientation + table.labels_append('grainID_{}@{:g}'.format(options.quaternion,options.disorientation)) # report orientation source and disorientation table.head_write() # ------------------------------------------ build KD tree ----------------------------------------- From aefa006d090fe2ebecd7c1b8b7e808ec756fdc15 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 9 Dec 2018 10:53:53 +0100 Subject: [PATCH 38/70] simplified assume quaternions as input --- processing/post/addSchmidfactors.py | 64 +++-------------------------- 1 file changed, 5 insertions(+), 59 deletions(-) diff --git a/processing/post/addSchmidfactors.py b/processing/post/addSchmidfactors.py index 9c5e9b0b4..a64708e1e 100755 --- a/processing/post/addSchmidfactors.py +++ b/processing/post/addSchmidfactors.py @@ -123,46 +123,19 @@ parser.add_option('-n', '--normal', dest = 'normal', type = 'float', nargs = 3, metavar = 'float float float', help = 'stress plane normal in lab frame [%default]') -parser.add_option('-e', '--eulers', - dest = 'eulers', - type = 'string', metavar = 'string', - help = 'Euler angles label') -parser.add_option('-d', '--degrees', - dest = 'degrees', - action = 'store_true', - help = 'Euler angles are given in degrees [%default]') -parser.add_option('-m', '--matrix', - dest = 'matrix', - type = 'string', metavar = 'string', - help = 'orientation matrix label') -parser.add_option('-a', - dest = 'a', - type = 'string', metavar = 'string', - help = 'crystal frame a vector label') -parser.add_option('-b', - dest = 'b', - type = 'string', metavar = 'string', - help = 'crystal frame b vector label') -parser.add_option('-c', - dest = 'c', - type = 'string', metavar = 'string', - help = 'crystal frame c vector label') parser.add_option('-q', '--quaternion', dest = 'quaternion', - type = 'string', metavar = 'string', + metavar = 'string', help = 'quaternion label') parser.set_defaults(force = (0.0,0.0,1.0), normal = None, lattice = latticeChoices[0], CoverA = math.sqrt(8./3.), - degrees = False, ) (options, filenames) = parser.parse_args() -toRadians = math.pi/180.0 if options.degrees else 1.0 # rescale degrees to radians - force = np.array(options.force) force /= np.linalg.norm(force) @@ -174,22 +147,6 @@ if options.normal: else: normal = force -input = [options.eulers is not None, - options.a is not None and \ - options.b is not None and \ - options.c is not None, - options.matrix is not None, - options.quaternion is not None, - ] - -if np.sum(input) != 1: parser.error('needs exactly one input format.') - -(label,dim,inputtype) = [(options.eulers,3,'eulers'), - ([options.a,options.b,options.c],[3,3,3],'frame'), - (options.matrix,9,'matrix'), - (options.quaternion,4,'quaternion'), - ][np.where(input)[0][0]] # select input label that was requested - slip_direction = np.zeros((len(slipSystems[options.lattice]),3),'f') slip_normal = np.zeros_like(slip_direction) @@ -227,13 +184,12 @@ for name in filenames: table.head_read() # ------------------------------------------ sanity checks ---------------------------------------- - - if not np.all(table.label_dimension(label) == dim): - damask.util.croak('input {} does not have dimension {}.'.format(label,dim)) + if not table.label_dimension(options.quaternion) == 4: + damask.util.croak('input {} does not have dimension 4.'.format(options.quaternion)) table.close(dismiss = True) # close ASCIItable and remove empty file continue - column = table.label_index(label) + column = table.label_index(options.quaternion) # ------------------------------------------ assemble header --------------------------------------- @@ -251,17 +207,7 @@ for name in filenames: outputAlive = True while outputAlive and table.data_read(): # read next data line of ASCII table - if inputtype == 'eulers': - o = damask.Orientation(Eulers = np.array(list(map(float,table.data[column:column+3])))*toRadians,) - elif inputtype == 'matrix': - o = damask.Orientation(matrix = np.array(list(map(float,table.data[column:column+9]))).reshape(3,3).transpose(),) - elif inputtype == 'frame': - o = damask.Orientation(matrix = np.array(list(map(float,table.data[column[0]:column[0]+3] + \ - table.data[column[1]:column[1]+3] + \ - table.data[column[2]:column[2]+3]))).reshape(3,3),) - elif inputtype == 'quaternion': - o = damask.Orientation(quaternion = np.array(list(map(float,table.data[column:column+4]))),) - + o = damask.Orientation(quaternion = np.array(list(map(float,table.data[column:column+4])))) table.data_append( np.abs( np.sum(slip_direction * (o.quaternion * force) ,axis=1) \ * np.sum(slip_normal * (o.quaternion * normal),axis=1))) From f7af4eba9cb6ec8719311080aa4ed0c090fa2a1b Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 9 Dec 2018 10:59:48 +0100 Subject: [PATCH 39/70] updated + new tests all tests python3 compatible and with quaternions for orientations --- PRIVATE | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PRIVATE b/PRIVATE index 64c39e902..efa9fe4e9 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 64c39e90200a8dd90dff3665e27050d9db81ac0d +Subproject commit efa9fe4e97cf2dfebff97f3a07503fc277dabf8a From 4a58fd98f96a6da70e0dec442ccf5e5c15cafda5 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 9 Dec 2018 11:02:28 +0100 Subject: [PATCH 40/70] pep257 compatible --- lib/damask/util.py | 1 + 1 file changed, 1 insertion(+) diff --git a/lib/damask/util.py b/lib/damask/util.py index 39715730c..754dd7c5d 100644 --- a/lib/damask/util.py +++ b/lib/damask/util.py @@ -137,6 +137,7 @@ class extendableOption(Option): def print_progress(iteration, total, prefix='', suffix='', decimals=1, bar_length=100): """ Call in a loop to create terminal progress bar + @params: iteration - Required : current iteration (Int) total - Required : total iterations (Int) From 741032339e43bd880fdeb6e72e87d267f84d2d51 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 9 Dec 2018 11:03:18 +0100 Subject: [PATCH 41/70] time/math/datatype are not used/defined any more --- processing/post/addGrainID.py | 2 +- processing/post/addIPFcolor.py | 2 +- processing/post/addPole.py | 2 +- processing/post/rotateData.py | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/processing/post/addGrainID.py b/processing/post/addGrainID.py index ac2c53a02..2eb3aec3a 100755 --- a/processing/post/addGrainID.py +++ b/processing/post/addGrainID.py @@ -1,7 +1,7 @@ #!/usr/bin/env python3 # -*- coding: UTF-8 no BOM -*- -import os,sys,time,copy +import os,sys,copy import numpy as np import damask from optparse import OptionParser diff --git a/processing/post/addIPFcolor.py b/processing/post/addIPFcolor.py index 1e5d8c2f3..2916df703 100755 --- a/processing/post/addIPFcolor.py +++ b/processing/post/addIPFcolor.py @@ -1,7 +1,7 @@ #!/usr/bin/env python3 # -*- coding: UTF-8 no BOM -*- -import os,sys,math +import os,sys import numpy as np from optparse import OptionParser import damask diff --git a/processing/post/addPole.py b/processing/post/addPole.py index 83a031d46..f2e9b978d 100755 --- a/processing/post/addPole.py +++ b/processing/post/addPole.py @@ -1,7 +1,7 @@ #!/usr/bin/env python3 # -*- coding: UTF-8 no BOM -*- -import os,sys,math +import os,sys import numpy as np from optparse import OptionParser import damask diff --git a/processing/post/rotateData.py b/processing/post/rotateData.py index 7ed9542cf..95102345b 100755 --- a/processing/post/rotateData.py +++ b/processing/post/rotateData.py @@ -66,7 +66,7 @@ for name in filenames: for i,dim in enumerate(table.label_dimension(options.data)): label = options.data[i] if dim == -1: - remarks.append('{} "{}" not found...'.format(datatype,me)) + remarks.append('"{}" not found...'.format(label)) elif dim == 3: remarks.append('adding vector "{}"...'.format(label)) active['vector'].append(label) From 50542915579ef549381a22d863c03d39fc2e29ac Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 9 Dec 2018 12:25:31 +0100 Subject: [PATCH 42/70] python3 compatible status message --- processing/post/postResults.py | 51 +++++++++++++++------------------- 1 file changed, 22 insertions(+), 29 deletions(-) diff --git a/processing/post/postResults.py b/processing/post/postResults.py index c6ed25d4b..74477d6a8 100755 --- a/processing/post/postResults.py +++ b/processing/post/postResults.py @@ -803,12 +803,6 @@ if not options.constitutiveResult: options.constitutiveResult = [] options.sort.reverse() options.sep.reverse() -# --- start background messaging - -if options.verbose: - bg = damask.util.backgroundMessage() - bg.start() - # --- parse .output and .t16 files if os.path.splitext(files[0])[1] == '': @@ -825,18 +819,13 @@ me = { 'Constitutive': options.phase, } -if options.verbose: bg.set_message('parsing .output files...') - for what in me: outputFormat[what] = ParseOutputFormat(filename, what, me[what]) if '_id' not in outputFormat[what]['specials']: print("\nsection '{}' not found in <{}>".format(me[what], what)) print('\n'.join(map(lambda x:' [%s]'%x, outputFormat[what]['specials']['brothers']))) -if options.verbose: bg.set_message('opening result file...') - p = OpenPostfile(filename+extension,options.filetype,options.nodal) -if options.verbose: bg.set_message('parsing result file...') stat = ParsePostfile(p, filename, outputFormat) if options.filetype == 'marc': stat['NumberOfIncrements'] -= 1 # t16 contains one "virtual" increment (at 0) @@ -879,8 +868,10 @@ if options.info: # --- build connectivity maps elementsOfNode = {} -for e in range(stat['NumberOfElements']): - if options.verbose and e%1000 == 0: bg.set_message('connect elem %i...'%e) +Nelems = stat['NumberOfElements'] +for e in range(Nelems): + if Nelems > 100 and e%(Nelems//100) == 0: # report in 1% steps if possible and avoid modulo by zero + damask.util.print_progress(iteration=e,total=Nelems,prefix='1/3: connecting elements') for n in map(p.node_sequence,p.element(e).items): if n not in elementsOfNode: elementsOfNode[n] = [p.element_id(e)] @@ -899,10 +890,12 @@ index = {} groups = [] groupCount = 0 memberCount = 0 - +print('\n') if options.nodalScalar: - for n in range(stat['NumberOfNodes']): - if options.verbose and n%1000 == 0: bg.set_message('scan node %i...'%n) + Npoints = stat['NumberOfNodes'] + for n in range(Npoints): + if Npoints > 100 and e%(Npoints//100) == 0: # report in 1% steps if possible and avoid modulo by zero + damask.util.print_progress(iteration=n,total=Npoints,prefix='2/3: scanning nodes ') myNodeID = p.node_id(n) myNodeCoordinates = [p.node(n).x, p.node(n).y, p.node(n).z] myElemID = 0 @@ -933,10 +926,13 @@ if options.nodalScalar: myNodeCoordinates) # incrementally update average location groups[index[grp]].append([myElemID,myNodeID,myIpID,myGrainID,0]) # append a new list defining each group member memberCount += 1 + print('\n') else: - for e in range(stat['NumberOfElements']): - if options.verbose and e%1000 == 0: bg.set_message('scan elem %i...'%e) + Nelems = stat['NumberOfElements'] + for e in range(Nelems): + if Nelems > 100 and e%(Nelems//100) == 0: # report in 1% steps if possible and avoid modulo by zero + damask.util.print_progress(iteration=e,total=Nelems,prefix='2/3: scanning elements ') myElemID = p.element_id(e) myIpCoordinates = ipCoords(p.element(e).type, list(map(lambda node: [node.x, node.y, node.z], list(map(p.node, map(p.node_sequence, p.element(e).items)))))) @@ -976,6 +972,7 @@ else: myIpCoordinates[n]) # incrementally update average location groups[index[grp]].append([myElemID,myNodeID,myIpID,myGrainID,n]) # append a new list defining each group member memberCount += 1 + print('\n') # --------------------------- sort groups -------------------------------- @@ -1002,7 +999,6 @@ if 'none' not in map(str.lower, options.sort): theKeys.append('x[0][%i]'%where[criterium]) sortKeys = eval('lambda x:(%s)'%(','.join(theKeys))) -if options.verbose: bg.set_message('sorting groups...') groups.sort(key = sortKeys) # in-place sorting to save mem @@ -1021,8 +1017,6 @@ standard = ['inc'] + \ # --------------------------- loop over positions -------------------------------- -if options.verbose: bg.set_message('getting map between positions and increments...') - incAtPosition = {} positionOfInc = {} @@ -1048,8 +1042,8 @@ increments = [incAtPosition[x] for x in locations] # build list of increments to time_start = time.time() +Nincs = len([i for i in locations]) for incCount,position in enumerate(locations): # walk through locations - p.moveto(position+offset_pos) # wind to correct position # --------------------------- file management -------------------------------- @@ -1075,16 +1069,14 @@ for incCount,position in enumerate(locations): # walk through locations # --------------------------- read and map data per group -------------------------------- member = 0 - for group in groups: - + Ngroups = len(groups) + for j,group in enumerate(groups): + f = incCount*Ngroups + j + if (Ngroups*Nincs) > 100 and f%((Ngroups*Nincs)//100) == 0: # report in 1% steps if possible and avoid modulo by zero + damask.util.print_progress(iteration=f,total=Ngroups*Nincs,prefix='3/3: processing points ') N = 0 # group member counter for (e,n,i,g,n_local) in group[1:]: # loop over group members member += 1 - if member%1000 == 0: - time_delta = ((len(locations)*memberCount)/float(member+incCount*memberCount)-1.0)*(time.time()-time_start) - if options.verbose: bg.set_message('(%02i:%02i:%02i) processing point %i of %i from increment %i (position %i)...' - %(time_delta//3600,time_delta%3600//60,time_delta%60,member,memberCount,increments[incCount],position)) - newby = [] # current member's data if options.nodalScalar: @@ -1172,6 +1164,7 @@ for incCount,position in enumerate(locations): # walk through locations group[0] + \ mappedResult) )) + '\n') +print('') if fileOpen: file.close() From 697d97cd3839d21f8e298c4a501fdb2fdf937960 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 9 Dec 2018 12:30:07 +0100 Subject: [PATCH 43/70] cleaned legacy format was for pre-MPI output where the extra bytes at start and end of a chunk (Fortran out) had to be handled manually --- processing/post/postResults.py | 66 +++++++--------------------------- 1 file changed, 12 insertions(+), 54 deletions(-) diff --git a/processing/post/postResults.py b/processing/post/postResults.py index 74477d6a8..66142be43 100755 --- a/processing/post/postResults.py +++ b/processing/post/postResults.py @@ -121,12 +121,8 @@ class MPIEspectral_result: # mimic py_post result object self._logscales = self._keyedPackedArray('logscales',count=self.N_loadcases,type='i') self.size = self._keyedPackedArray('size:',count=3,type='d') - if self.size == [None,None,None]: # no 'size' found, try legacy alias 'dimension' - self.size = self._keyedPackedArray('dimension',count=3,type='d') self.grid = self._keyedPackedArray('grid:',count=3,type='i') - if self.grid == [None,None,None]: # no 'grid' found, try legacy alias 'resolution' - self.grid = self._keyedPackedArray('resolution',count=3,type='i') self.N_nodes = (self.grid[0]+1)*(self.grid[1]+1)*(self.grid[2]+1) self.N_elements = self.grid[0] * self.grid[1] * self.grid[2] @@ -142,13 +138,8 @@ class MPIEspectral_result: # mimic py_post result object # parameters for file handling depending on output format - if options.legacy: - self.tagLen=8 - self.fourByteLimit = 2**31 -1 -8 - else: - self.tagLen=0 + self.tagLen=0 self.expectedFileSize = self.dataOffset+self.N_increments*(self.tagLen+self.N_elements*self.N_element_scalars*8) - if options.legacy: self.expectedFileSize+=self.expectedFileSize//self.fourByteLimit*8 # add extra 8 bytes for additional headers at 4 GB limits if self.expectedFileSize != self.filesize: print('\n**\n* Unexpected file size. Incomplete simulation or file corrupted!\n**') @@ -280,42 +271,16 @@ class MPIEspectral_result: # mimic py_post result object return self.N_element_scalars def element_scalar(self,e,idx): - if not options.legacy: - incStart = self.dataOffset \ - + self.position*8*self.N_elements*self.N_element_scalars - where = (e*self.N_element_scalars + idx)*8 - try: - self.file.seek(incStart+where) - value = struct.unpack('d',self.file.read(8))[0] - except: - print('seeking {}'.format(incStart+where)) - print('e {} idx {}'.format(e,idx)) - sys.exit(1) - - else: - self.fourByteLimit = 2**31 -1 -8 -# header & footer + extra header and footer for 4 byte int range (Fortran) -# values - incStart = self.dataOffset \ - + self.position*8*( 1 + self.N_elements*self.N_element_scalars*8//self.fourByteLimit \ - + self.N_elements*self.N_element_scalars) - - where = (e*self.N_element_scalars + idx)*8 - try: - if where%self.fourByteLimit + 8 >= self.fourByteLimit: # danger of reading into fortran record footer at 4 byte limit - data='' - for i in range(8): - self.file.seek(incStart+where+(where//self.fourByteLimit)*8+4) - data += self.file.read(1) - where += 1 - value = struct.unpack('d',data)[0] - else: - self.file.seek(incStart+where+(where//self.fourByteLimit)*8+4) - value = struct.unpack('d',self.file.read(8))[0] - except: - print('seeking {}'.format(incStart+where+(where//self.fourByteLimit)*8+4)) - print('e {} idx {}'.format(e,idx)) - sys.exit(1) + incStart = self.dataOffset \ + + self.position*8*self.N_elements*self.N_element_scalars + where = (e*self.N_element_scalars + idx)*8 + try: + self.file.seek(incStart+where) + value = struct.unpack('d',self.file.read(8))[0] + except: + print('seeking {}'.format(incStart+where)) + print('e {} idx {}'.format(e,idx)) + sys.exit(1) return [elemental_scalar(node,value) for node in self.element(e).items] @@ -645,8 +610,6 @@ of already processed data points for evaluation. parser.add_option('-i','--info', action='store_true', dest='info', help='list contents of resultfile') -parser.add_option('-l','--legacy', action='store_true', dest='legacy', - help='data format of spectral solver is in legacy format (no MPI out)') parser.add_option('-n','--nodal', action='store_true', dest='nodal', help='data is extrapolated to nodal value') parser.add_option( '--prefix', dest='prefix', @@ -673,10 +636,7 @@ parser.add_option('-p','--type', dest='filetype', help = 'type of result file [auto]') parser.add_option('-q','--quiet', dest='verbose', action = 'store_false', - help = 'suppress verbose output') -parser.add_option('--verbose', dest='verbose', - action = 'store_true', - help = 'enable verbose output') + help = 'legacy switch, no effect') group_material = OptionGroup(parser,'Material identifier') @@ -718,8 +678,6 @@ parser.add_option_group(group_general) parser.add_option_group(group_special) parser.set_defaults(info = False, - verbose = False, - legacy = False, nodal = False, prefix = '', suffix = '', From 65165ffc5ed31fd9159ad2b9f304d792eac2fae6 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 9 Dec 2018 12:35:32 +0100 Subject: [PATCH 44/70] assumed wrong type quaternion should have lenght one --- processing/post/addOrientations.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/processing/post/addOrientations.py b/processing/post/addOrientations.py index efd4dac24..27cc92b1f 100755 --- a/processing/post/addOrientations.py +++ b/processing/post/addOrientations.py @@ -14,21 +14,23 @@ scriptID = ' '.join([scriptName,damask.version]) # -------------------------------------------------------------------- def check_Eulers(eulers): - if np.any(eulers < 0.0) or np.any(eulers > 2.0*math.pi) or eulers[1] > math.pi: # Euler angles within valid range? + if np.any(eulers < 0.0) or np.any(eulers > 2.0*math.pi) or eulers[1] > math.pi: # Euler angles within valid range? raise ValueError('Euler angles outside of [0..2π],[0..π],[0..2π].\n{} {} {}.'.format(*eulers)) return eulers def check_quaternion(q): - if q.q < 0.0: # positive first quaternion component? + if q[0] < 0.0: # positive first quaternion component? raise ValueError('quaternion has negative first component.\n{}'.format(q)) + if not(np.isclose(np.linalg.norm(q), 1.0)): # unit quaternion? + raise ValueError('quaternion is not of unit length.\n{} {} {} {}'.format(*q)) return q def check_matrix(M): - if abs(1.0-np.linalg.det(M)) > 1e-8: # proper rotation? + if abs(1.0-np.linalg.det(M)) > 1e-8: # proper rotation? raise ValueError('matrix is not a proper rotation.\n{}'.format(M)) if abs(np.dot(M[0],M[1])) > 1e-8 \ or abs(np.dot(M[1],M[2])) > 1e-8 \ - or abs(np.dot(M[2],M[0])) > 1e-8: # all orthogonal? + or abs(np.dot(M[2],M[0])) > 1e-8: # all orthogonal? raise ValueError('matrix is not orthogonal.\n{}'.format(M)) return M From d43dfa569b0fc08569b17ba1db0b260d18706588 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 9 Dec 2018 12:37:27 +0100 Subject: [PATCH 45/70] lines too long --- processing/post/postResults.py | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/processing/post/postResults.py b/processing/post/postResults.py index 66142be43..e075acdf5 100755 --- a/processing/post/postResults.py +++ b/processing/post/postResults.py @@ -852,7 +852,7 @@ print('\n') if options.nodalScalar: Npoints = stat['NumberOfNodes'] for n in range(Npoints): - if Npoints > 100 and e%(Npoints//100) == 0: # report in 1% steps if possible and avoid modulo by zero + if Npoints > 100 and e%(Npoints//100) == 0: # report in 1% steps if possible and avoid modulo by zero damask.util.print_progress(iteration=n,total=Npoints,prefix='2/3: scanning nodes ') myNodeID = p.node_id(n) myNodeCoordinates = [p.node(n).x, p.node(n).y, p.node(n).z] @@ -862,27 +862,27 @@ if options.nodalScalar: # generate an expression that is only true for the locations specified by options.filter filter = substituteLocation(options.filter, [myElemID,myNodeID,myIpID,myGrainID], myNodeCoordinates) - if filter != '' and not eval(filter): # for all filter expressions that are not true:... - continue # ... ignore this data point and continue with next + if filter != '' and not eval(filter): # for all filter expressions that are not true:... + continue # ... ignore this data point and continue with next # --- group data locations # generate a unique key for a group of separated data based on the separation criterium for the location grp = substituteLocation('#'.join(options.sep), [myElemID,myNodeID,myIpID,myGrainID], myNodeCoordinates) - if grp not in index: # create a new group if not yet present + if grp not in index: # create a new group if not yet present index[grp] = groupCount - groups.append([[0,0,0,0,0.0,0.0,0.0]]) # initialize with avg location + groups.append([[0,0,0,0,0.0,0.0,0.0]]) # initialize with avg location groupCount += 1 groups[index[grp]][0][:4] = mapIncremental('','unique', len(groups[index[grp]])-1, groups[index[grp]][0][:4], - [myElemID,myNodeID,myIpID,myGrainID]) # keep only if unique average location + [myElemID,myNodeID,myIpID,myGrainID]) # keep only if unique average location groups[index[grp]][0][4:] = mapIncremental('','avg', len(groups[index[grp]])-1, groups[index[grp]][0][4:], - myNodeCoordinates) # incrementally update average location - groups[index[grp]].append([myElemID,myNodeID,myIpID,myGrainID,0]) # append a new list defining each group member + myNodeCoordinates) # incrementally update average location + groups[index[grp]].append([myElemID,myNodeID,myIpID,myGrainID,0]) # append a new list defining each group member memberCount += 1 print('\n') @@ -1030,7 +1030,7 @@ for incCount,position in enumerate(locations): # walk through locations Ngroups = len(groups) for j,group in enumerate(groups): f = incCount*Ngroups + j - if (Ngroups*Nincs) > 100 and f%((Ngroups*Nincs)//100) == 0: # report in 1% steps if possible and avoid modulo by zero + if (Ngroups*Nincs) > 100 and f%((Ngroups*Nincs)//100) == 0: # report in 1% steps if possible and avoid modulo by zero damask.util.print_progress(iteration=f,total=Ngroups*Nincs,prefix='3/3: processing points ') N = 0 # group member counter for (e,n,i,g,n_local) in group[1:]: # loop over group members From 130fac21c6dd676072d5b06d5a15ee5c6f03b801 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 9 Dec 2018 13:02:07 +0100 Subject: [PATCH 46/70] cleaning --- lib/damask/util.py | 35 +++++++++++++++++----------------- processing/post/addGrainID.py | 14 ++++++++++---- processing/post/postResults.py | 22 ++++++++++++--------- 3 files changed, 40 insertions(+), 31 deletions(-) diff --git a/lib/damask/util.py b/lib/damask/util.py index 754dd7c5d..a701a1288 100644 --- a/lib/damask/util.py +++ b/lib/damask/util.py @@ -135,27 +135,26 @@ class extendableOption(Option): # Print iterations progress # from https://gist.github.com/aubricus/f91fb55dc6ba5557fbab06119420dd6a def print_progress(iteration, total, prefix='', suffix='', decimals=1, bar_length=100): - """ - Call in a loop to create terminal progress bar + """ + Call in a loop to create terminal progress bar - @params: - iteration - Required : current iteration (Int) - total - Required : total iterations (Int) - prefix - Optional : prefix string (Str) - suffix - Optional : suffix string (Str) - decimals - Optional : positive number of decimals in percent complete (Int) - bar_length - Optional : character length of bar (Int) - """ - str_format = "{0:." + str(decimals) + "f}" - percents = str_format.format(100 * (iteration / float(total))) - filled_length = int(round(bar_length * iteration / float(total))) - bar = '█' * filled_length + '-' * (bar_length - filled_length) + @params: + iteration - Required : current iteration (Int) + total - Required : total iterations (Int) + prefix - Optional : prefix string (Str) + suffix - Optional : suffix string (Str) + decimals - Optional : positive number of decimals in percent complete (Int) + bar_length - Optional : character length of bar (Int) + """ + str_format = "{0:." + str(decimals) + "f}" + percents = str_format.format(100 * (iteration / float(total))) + filled_length = int(round(bar_length * iteration / float(total))) + bar = '█' * filled_length + '-' * (bar_length - filled_length) - sys.stdout.write('\r%s |%s| %s%s %s' % (prefix, bar, percents, '%', suffix)), + sys.stderr.write('\r%s |%s| %s%s %s' % (prefix, bar, percents, '%', suffix)), - if iteration == total: - sys.stdout.write('\n') -sys.stdout.flush() + if iteration == total: sys.stderr.write('\n\n') + sys.stderr.flush() # ----------------------------- class backgroundMessage(threading.Thread): diff --git a/processing/post/addGrainID.py b/processing/post/addGrainID.py index 2eb3aec3a..2c4cfe5d6 100755 --- a/processing/post/addGrainID.py +++ b/processing/post/addGrainID.py @@ -29,20 +29,25 @@ parser.add_option('-d', parser.add_option('-s', '--symmetry', dest = 'symmetry', - type = 'string', metavar = 'string', + metavar = 'string', help = 'crystal symmetry [%default]') parser.add_option('-q', '--quaternion', dest = 'quaternion', - type = 'string', metavar = 'string', + metavar = 'string', help = 'label of quaternion') parser.add_option('-p', '--pos', '--position', dest = 'pos', - type = 'string', metavar = 'string', + metavar = 'string', help = 'label of coordinates [%default]') +parser.add_option('--quiet', + dest='verbose', + action = 'store_false', + help = 'hide status bar (useful when piping to file)') parser.set_defaults(disorientation = 5, + verbose = True, quaternion = 'orientation', symmetry = 'cubic', pos = 'pos', @@ -111,7 +116,7 @@ for name in filenames: table.data_rewind() while table.data_read(): # read next data line of ASCII table - if Npoints > 100 and p%(Npoints//100) == 0: # report in 1% steps if possible and avoid modulo by zero + if options.verbose and Npoints > 100 and p%(Npoints//100) == 0: # report in 1% steps if possible and avoid modulo by zero damask.util.print_progress(iteration=p,total=Npoints) o = damask.Orientation(quaternion = np.array(list(map(float,table.data[column:column+4]))), @@ -160,6 +165,7 @@ for name in filenames: outputAlive = True p = 0 + damask.util.print_progress(iteration=1,total=1) while outputAlive and table.data_read(): # read next data line of ASCII table table.data_append(1+packingMap[grainID[p]]) # add (condensed) grain ID outputAlive = table.data_write() # output processed line diff --git a/processing/post/postResults.py b/processing/post/postResults.py index e075acdf5..6ccdd7c3f 100755 --- a/processing/post/postResults.py +++ b/processing/post/postResults.py @@ -636,7 +636,7 @@ parser.add_option('-p','--type', dest='filetype', help = 'type of result file [auto]') parser.add_option('-q','--quiet', dest='verbose', action = 'store_false', - help = 'legacy switch, no effect') + help = 'hide status bar (useful when piping to file)') group_material = OptionGroup(parser,'Material identifier') @@ -679,6 +679,7 @@ parser.add_option_group(group_special) parser.set_defaults(info = False, nodal = False, + verbose = True, prefix = '', suffix = '', dir = 'postProc', @@ -705,6 +706,8 @@ if files == []: parser.print_help() parser.error('no file specified...') +damask.util.report(scriptName,files[0]) + if not os.path.exists(files[0]): parser.print_help() parser.error('invalid file "%s" specified...'%files[0]) @@ -828,7 +831,7 @@ if options.info: elementsOfNode = {} Nelems = stat['NumberOfElements'] for e in range(Nelems): - if Nelems > 100 and e%(Nelems//100) == 0: # report in 1% steps if possible and avoid modulo by zero + if options.verbose and Nelems > 100 and e%(Nelems//100) == 0: # report in 1% steps if possible and avoid modulo by zero damask.util.print_progress(iteration=e,total=Nelems,prefix='1/3: connecting elements') for n in map(p.node_sequence,p.element(e).items): if n not in elementsOfNode: @@ -848,11 +851,12 @@ index = {} groups = [] groupCount = 0 memberCount = 0 -print('\n') +damask.util.print_progress(iteration=1,total=1,prefix='1/3: connecting elements') + if options.nodalScalar: Npoints = stat['NumberOfNodes'] for n in range(Npoints): - if Npoints > 100 and e%(Npoints//100) == 0: # report in 1% steps if possible and avoid modulo by zero + if options.verbose and Npoints > 100 and e%(Npoints//100) == 0: # report in 1% steps if possible and avoid modulo by zero damask.util.print_progress(iteration=n,total=Npoints,prefix='2/3: scanning nodes ') myNodeID = p.node_id(n) myNodeCoordinates = [p.node(n).x, p.node(n).y, p.node(n).z] @@ -884,12 +888,12 @@ if options.nodalScalar: myNodeCoordinates) # incrementally update average location groups[index[grp]].append([myElemID,myNodeID,myIpID,myGrainID,0]) # append a new list defining each group member memberCount += 1 - print('\n') + damask.util.print_progress(iteration=1,total=1,prefix='2/3: scanning nodes ') else: Nelems = stat['NumberOfElements'] for e in range(Nelems): - if Nelems > 100 and e%(Nelems//100) == 0: # report in 1% steps if possible and avoid modulo by zero + if options.verbose and Nelems > 100 and e%(Nelems//100) == 0: # report in 1% steps if possible and avoid modulo by zero damask.util.print_progress(iteration=e,total=Nelems,prefix='2/3: scanning elements ') myElemID = p.element_id(e) myIpCoordinates = ipCoords(p.element(e).type, list(map(lambda node: [node.x, node.y, node.z], @@ -930,7 +934,7 @@ else: myIpCoordinates[n]) # incrementally update average location groups[index[grp]].append([myElemID,myNodeID,myIpID,myGrainID,n]) # append a new list defining each group member memberCount += 1 - print('\n') + damask.util.print_progress(iteration=1,total=1,prefix='2/3: scanning elements ') # --------------------------- sort groups -------------------------------- @@ -1030,7 +1034,7 @@ for incCount,position in enumerate(locations): # walk through locations Ngroups = len(groups) for j,group in enumerate(groups): f = incCount*Ngroups + j - if (Ngroups*Nincs) > 100 and f%((Ngroups*Nincs)//100) == 0: # report in 1% steps if possible and avoid modulo by zero + if options.verbose and (Ngroups*Nincs) > 100 and f%((Ngroups*Nincs)//100) == 0: # report in 1% steps if possible and avoid modulo by zero damask.util.print_progress(iteration=f,total=Ngroups*Nincs,prefix='3/3: processing points ') N = 0 # group member counter for (e,n,i,g,n_local) in group[1:]: # loop over group members @@ -1122,7 +1126,7 @@ for incCount,position in enumerate(locations): # walk through locations group[0] + \ mappedResult) )) + '\n') -print('') +damask.util.print_progress(iteration=1,total=1,prefix='3/3: processing points ') if fileOpen: file.close() From 92c356a2c8fc9345a165850282a38f84d19a4436 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 9 Dec 2018 13:18:03 +0100 Subject: [PATCH 47/70] updated tests --- PRIVATE | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PRIVATE b/PRIVATE index efa9fe4e9..d15faafa8 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit efa9fe4e97cf2dfebff97f3a07503fc277dabf8a +Subproject commit d15faafa81e7133977dcc6a160ed73a58de69ccb From 08ac49fac026e91fd45d695de5f32ebfa852a839 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 9 Dec 2018 13:21:48 +0100 Subject: [PATCH 48/70] line too long --- processing/post/postResults.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/processing/post/postResults.py b/processing/post/postResults.py index 6ccdd7c3f..de78180c1 100755 --- a/processing/post/postResults.py +++ b/processing/post/postResults.py @@ -1034,7 +1034,7 @@ for incCount,position in enumerate(locations): # walk through locations Ngroups = len(groups) for j,group in enumerate(groups): f = incCount*Ngroups + j - if options.verbose and (Ngroups*Nincs) > 100 and f%((Ngroups*Nincs)//100) == 0: # report in 1% steps if possible and avoid modulo by zero + if options.verbose and (Ngroups*Nincs) > 100 and f%((Ngroups*Nincs)//100) == 0: # report in 1% steps if possible and avoid modulo by zero damask.util.print_progress(iteration=f,total=Ngroups*Nincs,prefix='3/3: processing points ') N = 0 # group member counter for (e,n,i,g,n_local) in group[1:]: # loop over group members From bc04b83d8ad5082545cdef8d4a5cf5c571dec7e0 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 9 Dec 2018 14:09:26 +0100 Subject: [PATCH 49/70] no reason for using python2.7 --- processing/post/addMises.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/processing/post/addMises.py b/processing/post/addMises.py index 4719c2e35..55cf6552e 100755 --- a/processing/post/addMises.py +++ b/processing/post/addMises.py @@ -1,4 +1,4 @@ -#!/usr/bin/env python2.7 +#!/usr/bin/env python3 # -*- coding: UTF-8 no BOM -*- import os,sys,math From 60686fb72c3191defbce5523f6e5c710f37ca7e6 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 9 Dec 2018 18:05:37 +0100 Subject: [PATCH 50/70] [skip ci] improved reporting --- processing/post/addOrientations.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/processing/post/addOrientations.py b/processing/post/addOrientations.py index 27cc92b1f..141083036 100755 --- a/processing/post/addOrientations.py +++ b/processing/post/addOrientations.py @@ -1,7 +1,7 @@ #!/usr/bin/env python3 # -*- coding: UTF-8 no BOM -*- -import os,sys,math +import os,sys import numpy as np from optparse import OptionParser import damask @@ -14,19 +14,19 @@ scriptID = ' '.join([scriptName,damask.version]) # -------------------------------------------------------------------- def check_Eulers(eulers): - if np.any(eulers < 0.0) or np.any(eulers > 2.0*math.pi) or eulers[1] > math.pi: # Euler angles within valid range? + if np.any(eulers < 0.0) or np.any(eulers > 2.0*np.pi) or eulers[1] > np.pi: # Euler angles within valid range? raise ValueError('Euler angles outside of [0..2π],[0..π],[0..2π].\n{} {} {}.'.format(*eulers)) return eulers def check_quaternion(q): if q[0] < 0.0: # positive first quaternion component? - raise ValueError('quaternion has negative first component.\n{}'.format(q)) - if not(np.isclose(np.linalg.norm(q), 1.0)): # unit quaternion? + raise ValueError('quaternion has negative first component.\n{}'.format(q[0])) + if not np.isclose(np.linalg.norm(q), 1.0): # unit quaternion? raise ValueError('quaternion is not of unit length.\n{} {} {} {}'.format(*q)) return q def check_matrix(M): - if abs(1.0-np.linalg.det(M)) > 1e-8: # proper rotation? + if not np.isclose(np.linalg.det(M),1.0): # proper rotation? raise ValueError('matrix is not a proper rotation.\n{}'.format(M)) if abs(np.dot(M[0],M[1])) > 1e-8 \ or abs(np.dot(M[1],M[2])) > 1e-8 \ @@ -129,7 +129,7 @@ if np.sum(input) != 1: parser.error('needs exactly one input format.') (options.quaternion,4,'quaternion'), ][np.where(input)[0][0]] # select input label that was requested -toRadians = math.pi/180.0 if options.degrees else 1.0 # rescale degrees to radians +toRadians = np.pi/180.0 if options.degrees else 1.0 # rescale degrees to radians r = damask.Quaternion.fromAngleAxis(toRadians*options.crystalrotation[0],options.crystalrotation[1:]) # crystal frame rotation R = damask.Quaternion.fromAngleAxis(toRadians*options. labrotation[0],options. labrotation[1:]) # lab frame rotation From d219842ad86b93392c06b3d2ee1cf674ae1e5fed Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 9 Dec 2018 20:19:17 +0100 Subject: [PATCH 51/70] [skip ci] consistent tolerances --- processing/post/addOrientations.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/processing/post/addOrientations.py b/processing/post/addOrientations.py index 141083036..803d585e3 100755 --- a/processing/post/addOrientations.py +++ b/processing/post/addOrientations.py @@ -28,9 +28,9 @@ def check_quaternion(q): def check_matrix(M): if not np.isclose(np.linalg.det(M),1.0): # proper rotation? raise ValueError('matrix is not a proper rotation.\n{}'.format(M)) - if abs(np.dot(M[0],M[1])) > 1e-8 \ - or abs(np.dot(M[1],M[2])) > 1e-8 \ - or abs(np.dot(M[2],M[0])) > 1e-8: # all orthogonal? + if not np.isclose(np.dot(M[0],M[1]), 0.0) \ + or not np.isclose(np.dot(M[1],M[2]), 0.0) \ + or not np.isclose(np.dot(M[2],M[0]), 0.0): # all orthogonal? raise ValueError('matrix is not orthogonal.\n{}'.format(M)) return M From c038dc9f93cdfc69b8ceccb80c62e737e3a9b5ed Mon Sep 17 00:00:00 2001 From: Test User Date: Mon, 10 Dec 2018 06:53:20 +0100 Subject: [PATCH 52/70] [skip ci] updated version information after successful test of v2.0.2-1124-g0ff9e9fa --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 3ad5cdbde..8817d8cc0 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.2-1122-g2349341e +v2.0.2-1124-g0ff9e9fa From 767650e0023176ac6b47720a9d5c2645db3f098a Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 10 Dec 2018 09:27:39 +0100 Subject: [PATCH 53/70] more general option + improved descriptions hopefully more user friendly --- PRIVATE | 2 +- processing/post/addGrainID.py | 6 ++--- processing/post/addIPFcolor.py | 13 ++++++----- processing/post/addOrientations.py | 34 ++++++++++++++++------------- processing/post/addPole.py | 14 +++++++----- processing/post/addSchmidfactors.py | 19 ++++++++++------ 6 files changed, 51 insertions(+), 37 deletions(-) diff --git a/PRIVATE b/PRIVATE index d15faafa8..e3dac27b7 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit d15faafa81e7133977dcc6a160ed73a58de69ccb +Subproject commit e3dac27b709d7fb3630bbd75271b220827221492 diff --git a/processing/post/addGrainID.py b/processing/post/addGrainID.py index 2c4cfe5d6..ed707f9ea 100755 --- a/processing/post/addGrainID.py +++ b/processing/post/addGrainID.py @@ -31,11 +31,11 @@ parser.add_option('-s', dest = 'symmetry', metavar = 'string', help = 'crystal symmetry [%default]') -parser.add_option('-q', - '--quaternion', +parser.add_option('-o', + '--orientation', dest = 'quaternion', metavar = 'string', - help = 'label of quaternion') + help = 'label of crystal orientation given as unit quaternion [%default]') parser.add_option('-p', '--pos', '--position', dest = 'pos', diff --git a/processing/post/addIPFcolor.py b/processing/post/addIPFcolor.py index 2916df703..c5a59a63a 100755 --- a/processing/post/addIPFcolor.py +++ b/processing/post/addIPFcolor.py @@ -18,18 +18,21 @@ Add RGB color value corresponding to TSL-OIM scheme for inverse pole figures. """, version = scriptID) -parser.add_option('-p', '--pole', +parser.add_option('-p', + '--pole', dest = 'pole', type = 'float', nargs = 3, metavar = 'float float float', help = 'lab frame direction for inverse pole figure [%default]') -parser.add_option('-s', '--symmetry', +parser.add_option('-s', + '--symmetry', dest = 'symmetry', type = 'choice', choices = damask.Symmetry.lattices[1:], metavar='string', help = 'crystal symmetry [%default] {{{}}} '.format(', '.join(damask.Symmetry.lattices[1:]))) -parser.add_option('-q', '--quaternion', +parser.add_option('-o', + '--orientation', dest = 'quaternion', - type = 'string', metavar = 'string', - help = 'quaternion label') + metavar = 'string', + help = 'label of crystal orientation given as unit quaternion [%default]') parser.set_defaults(pole = (0.0,0.0,1.0), quaternion = 'orientation', diff --git a/processing/post/addOrientations.py b/processing/post/addOrientations.py index 803d585e3..a33f96b91 100755 --- a/processing/post/addOrientations.py +++ b/processing/post/addOrientations.py @@ -54,49 +54,53 @@ outputChoices = { 'angleaxis': ['aaxs',4], } -parser.add_option('-o', '--output', +parser.add_option('-o', + '--output', dest = 'output', action = 'extend', metavar = '', help = 'output orientation formats {{{}}}'.format(', '.join(outputChoices))) -parser.add_option('-d', '--degrees', +parser.add_option('-d', + '--degrees', dest = 'degrees', action = 'store_true', help = 'all angles in degrees') -parser.add_option('-R', '--labrotation', +parser.add_option('-R', + '--labrotation', dest='labrotation', type = 'float', nargs = 4, metavar = ' '.join(['float']*4), help = 'angle and axis of additional lab frame rotation') -parser.add_option('-r', '--crystalrotation', +parser.add_option('-r', + '--crystalrotation', dest='crystalrotation', type = 'float', nargs = 4, metavar = ' '.join(['float']*4), help = 'angle and axis of additional crystal frame rotation') -parser.add_option( '--eulers', +parser.add_option('--eulers', dest = 'eulers', - type = 'string', metavar = 'string', + metavar = 'string', help = 'Euler angles label') -parser.add_option( '--rodrigues', +parser.add_option('--rodrigues', dest = 'rodrigues', - type = 'string', metavar = 'string', + metavar = 'string', help = 'Rodrigues vector label') -parser.add_option( '--matrix', +parser.add_option('--matrix', dest = 'matrix', - type = 'string', metavar = 'string', + metavar = 'string', help = 'orientation matrix label') -parser.add_option( '--quaternion', +parser.add_option('--quaternion', dest = 'quaternion', - type = 'string', metavar = 'string', + metavar = 'string', help = 'quaternion label') parser.add_option('-x', dest = 'x', - type = 'string', metavar = 'string', + metavar = 'string', help = 'label of lab x vector (expressed in crystal coords)') parser.add_option('-y', dest = 'y', - type = 'string', metavar = 'string', + metavar = 'string', help = 'label of lab y vector (expressed in crystal coords)') parser.add_option('-z', dest = 'z', - type = 'string', metavar = 'string', + metavar = 'string', help = 'label of lab z vector (expressed in crystal coords)') parser.set_defaults(output = [], diff --git a/processing/post/addPole.py b/processing/post/addPole.py index f2e9b978d..3098effc7 100755 --- a/processing/post/addPole.py +++ b/processing/post/addPole.py @@ -14,22 +14,24 @@ scriptID = ' '.join([scriptName,damask.version]) # -------------------------------------------------------------------- parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """ -Add x,y coordinates of stereographic projection of given direction (pole) in crystal frame. +Add coordinates of stereographic projection of given direction (pole) in crystal frame. """, version = scriptID) -parser.add_option('-p', '--pole', +parser.add_option('-p', + '--pole', dest = 'pole', type = 'float', nargs = 3, metavar = 'float float float', help = 'crystal frame direction for pole figure [%default]') parser.add_option('--polar', dest = 'polar', action = 'store_true', - help = 'output polar coordinates r,phi [%default]') -parser.add_option('-q', '--quaternion', + help = 'output polar coordinates (r,φ) instead of Cartesian coordinates (x,y)') +parser.add_option('-o', + '--orientation', dest = 'quaternion', - type = 'string', metavar = 'string', - help = 'quaternion label') + metavar = 'string', + help = 'label of crystal orientation given as unit quaternion [%default]') parser.set_defaults(pole = (1.0,0.0,0.0), quaternion = 'orientation', diff --git a/processing/post/addSchmidfactors.py b/processing/post/addSchmidfactors.py index a64708e1e..6335b419e 100755 --- a/processing/post/addSchmidfactors.py +++ b/processing/post/addSchmidfactors.py @@ -109,26 +109,31 @@ Add columns listing Schmid factors (and optional trace vector of selected system """, version = scriptID) latticeChoices = ('fcc','bcc','hex') -parser.add_option('-l','--lattice', +parser.add_option('-l', + '--lattice', dest = 'lattice', type = 'choice', choices = latticeChoices, metavar='string', help = 'type of lattice structure [%default] {}'.format(latticeChoices)) parser.add_option('--covera', dest = 'CoverA', type = 'float', metavar = 'float', help = 'C over A ratio for hexagonal systems') -parser.add_option('-f', '--force', +parser.add_option('-f', + '--force', dest = 'force', type = 'float', nargs = 3, metavar = 'float float float', help = 'force direction in lab frame [%default]') -parser.add_option('-n', '--normal', +parser.add_option('-n', + '--normal', dest = 'normal', type = 'float', nargs = 3, metavar = 'float float float', - help = 'stress plane normal in lab frame [%default]') -parser.add_option('-q', '--quaternion', + help = 'stress plane normal in lab frame, per default perpendicular to the force') +parser.add_option('-o', + '--orientation', dest = 'quaternion', metavar = 'string', - help = 'quaternion label') + help = 'label of crystal orientation given as unit quaternion [%default]') parser.set_defaults(force = (0.0,0.0,1.0), + quaternion='orientation', normal = None, lattice = latticeChoices[0], CoverA = math.sqrt(8./3.), @@ -139,7 +144,7 @@ parser.set_defaults(force = (0.0,0.0,1.0), force = np.array(options.force) force /= np.linalg.norm(force) -if options.normal: +if options.normal is not None: normal = np.array(options.normal) normal /= np.linalg.norm(normal) if abs(np.dot(force,normal)) > 1e-3: From 1f2fbbee21ff50c6fcf746655fd2a8ce0a3b234f Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 10 Dec 2018 09:38:57 +0100 Subject: [PATCH 54/70] consistent name --- lib/damask/util.py | 2 +- processing/post/addGrainID.py | 4 ++-- processing/post/postResults.py | 16 ++++++++-------- 3 files changed, 11 insertions(+), 11 deletions(-) diff --git a/lib/damask/util.py b/lib/damask/util.py index a701a1288..e7d4a16bc 100644 --- a/lib/damask/util.py +++ b/lib/damask/util.py @@ -134,7 +134,7 @@ class extendableOption(Option): # Print iterations progress # from https://gist.github.com/aubricus/f91fb55dc6ba5557fbab06119420dd6a -def print_progress(iteration, total, prefix='', suffix='', decimals=1, bar_length=100): +def progressBar(iteration, total, prefix='', suffix='', decimals=1, bar_length=100): """ Call in a loop to create terminal progress bar diff --git a/processing/post/addGrainID.py b/processing/post/addGrainID.py index ed707f9ea..3c4eaf4fa 100755 --- a/processing/post/addGrainID.py +++ b/processing/post/addGrainID.py @@ -117,7 +117,7 @@ for name in filenames: while table.data_read(): # read next data line of ASCII table if options.verbose and Npoints > 100 and p%(Npoints//100) == 0: # report in 1% steps if possible and avoid modulo by zero - damask.util.print_progress(iteration=p,total=Npoints) + damask.util.progressBar(iteration=p,total=Npoints) o = damask.Orientation(quaternion = np.array(list(map(float,table.data[column:column+4]))), symmetry = options.symmetry).reduced() @@ -165,7 +165,7 @@ for name in filenames: outputAlive = True p = 0 - damask.util.print_progress(iteration=1,total=1) + damask.util.progressBar(iteration=1,total=1) while outputAlive and table.data_read(): # read next data line of ASCII table table.data_append(1+packingMap[grainID[p]]) # add (condensed) grain ID outputAlive = table.data_write() # output processed line diff --git a/processing/post/postResults.py b/processing/post/postResults.py index de78180c1..a5a2669d7 100755 --- a/processing/post/postResults.py +++ b/processing/post/postResults.py @@ -832,7 +832,7 @@ elementsOfNode = {} Nelems = stat['NumberOfElements'] for e in range(Nelems): if options.verbose and Nelems > 100 and e%(Nelems//100) == 0: # report in 1% steps if possible and avoid modulo by zero - damask.util.print_progress(iteration=e,total=Nelems,prefix='1/3: connecting elements') + damask.util.progressBar(iteration=e,total=Nelems,prefix='1/3: connecting elements') for n in map(p.node_sequence,p.element(e).items): if n not in elementsOfNode: elementsOfNode[n] = [p.element_id(e)] @@ -851,13 +851,13 @@ index = {} groups = [] groupCount = 0 memberCount = 0 -damask.util.print_progress(iteration=1,total=1,prefix='1/3: connecting elements') +damask.util.progressBar(iteration=1,total=1,prefix='1/3: connecting elements') if options.nodalScalar: Npoints = stat['NumberOfNodes'] for n in range(Npoints): if options.verbose and Npoints > 100 and e%(Npoints//100) == 0: # report in 1% steps if possible and avoid modulo by zero - damask.util.print_progress(iteration=n,total=Npoints,prefix='2/3: scanning nodes ') + damask.util.progressBar(iteration=n,total=Npoints,prefix='2/3: scanning nodes ') myNodeID = p.node_id(n) myNodeCoordinates = [p.node(n).x, p.node(n).y, p.node(n).z] myElemID = 0 @@ -888,13 +888,13 @@ if options.nodalScalar: myNodeCoordinates) # incrementally update average location groups[index[grp]].append([myElemID,myNodeID,myIpID,myGrainID,0]) # append a new list defining each group member memberCount += 1 - damask.util.print_progress(iteration=1,total=1,prefix='2/3: scanning nodes ') + damask.util.progressBar(iteration=1,total=1,prefix='2/3: scanning nodes ') else: Nelems = stat['NumberOfElements'] for e in range(Nelems): if options.verbose and Nelems > 100 and e%(Nelems//100) == 0: # report in 1% steps if possible and avoid modulo by zero - damask.util.print_progress(iteration=e,total=Nelems,prefix='2/3: scanning elements ') + damask.util.progressBar(iteration=e,total=Nelems,prefix='2/3: scanning elements ') myElemID = p.element_id(e) myIpCoordinates = ipCoords(p.element(e).type, list(map(lambda node: [node.x, node.y, node.z], list(map(p.node, map(p.node_sequence, p.element(e).items)))))) @@ -934,7 +934,7 @@ else: myIpCoordinates[n]) # incrementally update average location groups[index[grp]].append([myElemID,myNodeID,myIpID,myGrainID,n]) # append a new list defining each group member memberCount += 1 - damask.util.print_progress(iteration=1,total=1,prefix='2/3: scanning elements ') + damask.util.progressBar(iteration=1,total=1,prefix='2/3: scanning elements ') # --------------------------- sort groups -------------------------------- @@ -1035,7 +1035,7 @@ for incCount,position in enumerate(locations): # walk through locations for j,group in enumerate(groups): f = incCount*Ngroups + j if options.verbose and (Ngroups*Nincs) > 100 and f%((Ngroups*Nincs)//100) == 0: # report in 1% steps if possible and avoid modulo by zero - damask.util.print_progress(iteration=f,total=Ngroups*Nincs,prefix='3/3: processing points ') + damask.util.progressBar(iteration=f,total=Ngroups*Nincs,prefix='3/3: processing points ') N = 0 # group member counter for (e,n,i,g,n_local) in group[1:]: # loop over group members member += 1 @@ -1126,7 +1126,7 @@ for incCount,position in enumerate(locations): # walk through locations group[0] + \ mappedResult) )) + '\n') -damask.util.print_progress(iteration=1,total=1,prefix='3/3: processing points ') +damask.util.progressBar(iteration=1,total=1,prefix='3/3: processing points ') if fileOpen: file.close() From 2aa6b12126b62c9c2d9e282f971877e2f71f36d1 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 12 Dec 2018 16:13:57 +0100 Subject: [PATCH 55/70] IMPORTANT Behavior change: Slip (Lp) happens in twinned volume fraction aliases for associate do not have to be defined --- src/plastic_phenopowerlaw.f90 | 52 ++++++++++++++--------------------- 1 file changed, 21 insertions(+), 31 deletions(-) diff --git a/src/plastic_phenopowerlaw.f90 b/src/plastic_phenopowerlaw.f90 index 5fd8d8e3c..136fe07e6 100644 --- a/src/plastic_phenopowerlaw.f90 +++ b/src/plastic_phenopowerlaw.f90 @@ -148,12 +148,6 @@ subroutine plastic_phenopowerlaw_init real(pReal), dimension(0), parameter :: emptyRealArray = [real(pReal)::] character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::] - type(tParameters) :: & - prm - type(tPhenopowerlawState) :: & - stt, & - dot - integer(kind(undefined_ID)) :: & outputID !< ID of each post result output @@ -161,7 +155,7 @@ subroutine plastic_phenopowerlaw_init structure = '',& extmsg = '' character(len=65536), dimension(:), allocatable :: & - outputs + outputs write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_PHENOPOWERLAW_label//' init -+>>>' write(6,'(a15,a)') ' Current time: ',IO_timeStamp() @@ -398,6 +392,8 @@ end subroutine plastic_phenopowerlaw_init !-------------------------------------------------------------------------------------------------- !> @brief calculates plastic velocity gradient and its tangent +!> @details asumme that deformation by dislocation glide affects twinned and untwinned volume +! equally (Taylor assumption). Twinning happens only in untwinned volume ( !-------------------------------------------------------------------------------------------------- subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) @@ -420,19 +416,15 @@ subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) gdot_slip_pos,gdot_slip_neg real(pReal), dimension(param(instance)%totalNtwin) :: & gdot_twin,dgdot_dtautwin - - type(tParameters) :: prm - type(tPhenopowerlawState) :: stt - - associate(prm => param(instance), stt => state(instance)) - + Lp = 0.0_pReal dLp_dMp = 0.0_pReal - + + associate(prm => param(instance), stt => state(instance)) + call kinetics_slip(prm,stt,of,Mp,gdot_slip_pos,gdot_slip_neg,dgdot_dtauslip_pos,dgdot_dtauslip_neg) slipSystems: do i = 1_pInt, prm%totalNslip - Lp = Lp + (1.0_pReal- sum(stt%gamma_twin(:,of)/prm%gamma_twin_char)) & ! only shear in untwinned volume - * (gdot_slip_pos(i)+gdot_slip_neg(i))*prm%Schmid_slip(1:3,1:3,i) + Lp = Lp + (gdot_slip_pos(i)+gdot_slip_neg(i))*prm%Schmid_slip(1:3,1:3,i) forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & + dgdot_dtauslip_pos(i) * prm%Schmid_slip(k,l,i) * prm%nonSchmid_pos(m,n,i) & @@ -446,9 +438,9 @@ subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & + dgdot_dtautwin(i)*prm%Schmid_twin(k,l,i)*prm%Schmid_twin(m,n,i) enddo twinSystems + end associate - end subroutine plastic_phenopowerlaw_LpAndItsTangent @@ -474,9 +466,6 @@ subroutine plastic_phenopowerlaw_dotState(Mp,instance,of) left_SlipSlip,right_SlipSlip, & gdot_slip_pos,gdot_slip_neg - type(tParameters) :: prm - type(tPhenopowerlawState) :: dot,stt - associate(prm => param(instance), stt => state(instance), dot => dotState(instance)) dot%whole(:,of) = 0.0_pReal @@ -506,7 +495,7 @@ subroutine plastic_phenopowerlaw_dotState(Mp,instance,of) ! hardening hardeningSlip: do i = 1_pInt, prm%totalNslip dot%xi_slip(i,of) = dot_product(prm%interaction_SlipSlip(i,:),right_SlipSlip*dot%gamma_slip(:,of)) & - * c_SlipSlip * left_SlipSlip(i) & + * c_SlipSlip * left_SlipSlip(i) & + dot_product(prm%interaction_SlipTwin(i,:),dot%gamma_twin(:,of)) enddo hardeningSlip @@ -522,8 +511,9 @@ end subroutine plastic_phenopowerlaw_dotState !-------------------------------------------------------------------------------------------------- !> @brief calculates shear rates on slip systems and derivatives with respect to resolved stress -!> @details: Shear rates are calculated only optionally. NOTE: Against the common convention, the -!> result (i.e. intent(out)) variables are the last to have the optional arguments at the end +!> @details Shear rates are calculated only optionally. +! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to +! have the optional arguments at the end !-------------------------------------------------------------------------------------------------- pure subroutine kinetics_slip(prm,stt,of,Mp,gdot_slip_pos,gdot_slip_neg, & dgdot_dtau_slip_pos,dgdot_dtau_slip_neg) @@ -595,9 +585,11 @@ end subroutine kinetics_slip !-------------------------------------------------------------------------------------------------- -!> @brief calculates shear rates on twin systems and derivatives with respect to resolved stress -!> @details: Derivates are calculated only optionally. NOTE: Against the common convention, the -!> result (i.e. intent(out)) variables are the last to have the optional arguments at the end +!> @brief calculates shear rates on twin systems and derivatives with respect to resolved stress. +! twinning is assumed to take place only in untwinned volume. +!> @details Derivates are calculated only optionally. +! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to +! have the optional arguments at the end. !-------------------------------------------------------------------------------------------------- pure subroutine kinetics_twin(prm,stt,of,Mp,gdot_twin,dgdot_dtau_twin) use prec, only: & @@ -667,13 +659,10 @@ function plastic_phenopowerlaw_postResults(Mp,instance,of) result(postResults) real(pReal), dimension(param(instance)%totalNslip) :: & gdot_slip_pos,gdot_slip_neg - type(tParameters) :: prm - type(tPhenopowerlawState) :: stt - - associate( prm => param(instance), stt => state(instance)) - postResults = 0.0_pReal c = 0_pInt + + associate( prm => param(instance), stt => state(instance)) outputsLoop: do o = 1_pInt,size(prm%outputID) select case(prm%outputID(o)) @@ -711,6 +700,7 @@ function plastic_phenopowerlaw_postResults(Mp,instance,of) result(postResults) end select enddo outputsLoop + end associate end function plastic_phenopowerlaw_postResults From 15f1badb60b915f8eeea60085608a88c4d809b75 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 13 Dec 2018 17:27:45 +0100 Subject: [PATCH 56/70] test with reference solution for modified phenopowerlaw --- PRIVATE | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PRIVATE b/PRIVATE index 2e0ab3e25..58137906b 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 2e0ab3e2564fea6ac3779623065b47dfc3261ef5 +Subproject commit 58137906b84b6cf0e273dfdde623a2986d03f98e From 4dd064a2753868b2a8bec0d9522ec892a9a9992e Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 14 Dec 2018 09:23:52 +0100 Subject: [PATCH 57/70] polishing --- DAMASK_prerequisites.sh | 2 +- src/plastic_phenopowerlaw.f90 | 12 ++++++------ 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/DAMASK_prerequisites.sh b/DAMASK_prerequisites.sh index f0a525c34..5e96f4aee 100755 --- a/DAMASK_prerequisites.sh +++ b/DAMASK_prerequisites.sh @@ -79,7 +79,7 @@ ls $PETSC_DIR/lib firstLevel "Python" DEFAULT_PYTHON=python3 -for executable in python python2 python3 python2.7; do +for executable in python python3; do getDetails $executable '--version' done secondLevel "Details on $DEFAULT_PYTHON:" diff --git a/src/plastic_phenopowerlaw.f90 b/src/plastic_phenopowerlaw.f90 index 136fe07e6..053fe958b 100644 --- a/src/plastic_phenopowerlaw.f90 +++ b/src/plastic_phenopowerlaw.f90 @@ -53,7 +53,7 @@ module plastic_phenopowerlaw xi_twin_0, & !< initial critical shear stress for twin xi_slip_sat, & !< maximum critical shear stress for slip nonSchmidCoeff, & - H_int, & !< per family hardening activity (optional) !ToDo: Better name! + H_int, & !< per family hardening activity (optional) gamma_twin_char !< characteristic shear for twins real(pReal), allocatable, dimension(:,:) :: & interaction_SlipSlip, & !< slip resistance from slip activity @@ -234,8 +234,8 @@ subroutine plastic_phenopowerlaw_init ! sanity checks if (prm%gdot0_slip <= 0.0_pReal) extmsg = trim(extmsg)//'gdot0_slip ' - if (dEq0(prm%a_slip)) extmsg = trim(extmsg)//'a_slip ' ! ToDo: negative values ok? - if (dEq0(prm%n_slip)) extmsg = trim(extmsg)//'n_slip ' ! ToDo: negative values ok? + if (dEq0(prm%a_slip)) extmsg = trim(extmsg)//'a_slip ' ! ToDo: negative values ok? + if (dEq0(prm%n_slip)) extmsg = trim(extmsg)//'n_slip ' ! ToDo: negative values ok? if (any(prm%xi_slip_0 <= 0.0_pReal)) extmsg = trim(extmsg)//'xi_slip_0 ' if (any(prm%xi_slip_sat < prm%xi_slip_0)) extmsg = trim(extmsg)//'xi_slip_sat ' else slipActive @@ -268,7 +268,7 @@ subroutine plastic_phenopowerlaw_init ! sanity checks if (prm%gdot0_twin <= 0.0_pReal) extmsg = trim(extmsg)//'gdot0_twin ' - if (dEq0(prm%n_twin)) extmsg = trim(extmsg)//'n_twin ' ! ToDo: negative values ok? + if (dEq0(prm%n_twin)) extmsg = trim(extmsg)//'n_twin ' ! ToDo: negative values ok? else twinActive allocate(prm%interaction_TwinTwin(0,0)) allocate(prm%xi_twin_0(0)) @@ -285,8 +285,8 @@ subroutine plastic_phenopowerlaw_init config_phase(p)%getFloats('interaction_twinslip'), & structure(1:3)) else slipAndTwinActive - allocate(prm%interaction_SlipTwin(prm%totalNslip,prm%TotalNtwin)) ! at least one dimension is 0 - allocate(prm%interaction_TwinSlip(prm%totalNtwin,prm%TotalNslip)) ! at least one dimension is 0 + allocate(prm%interaction_SlipTwin(prm%totalNslip,prm%TotalNtwin)) ! at least one dimension is 0 + allocate(prm%interaction_TwinSlip(prm%totalNtwin,prm%TotalNslip)) ! at least one dimension is 0 prm%h0_TwinSlip = 0.0_pReal endif slipAndTwinActive From afe34feff55e1d0b984ba178045e894e4f6e82fb Mon Sep 17 00:00:00 2001 From: Test User Date: Sat, 15 Dec 2018 00:55:09 +0100 Subject: [PATCH 58/70] [skip ci] updated version information after successful test of v2.0.2-1137-g4dd064a2 --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 8817d8cc0..99b9812a2 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.2-1124-g0ff9e9fa +v2.0.2-1137-g4dd064a2 From e82d3557d1d13327612944ab24ab79d5a1ae65e6 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Mon, 17 Dec 2018 09:55:02 -0500 Subject: [PATCH 59/70] removed Rodrigues class; removed quaternion "inverse" method; adopted Rowenhorst algorithm for angle--axis; use math.isclose(...,0.0) to check for ==0.0 --- lib/damask/__init__.py | 2 +- lib/damask/orientation.py | 81 +++++++++++---------------------------- 2 files changed, 24 insertions(+), 59 deletions(-) diff --git a/lib/damask/__init__.py b/lib/damask/__init__.py index 9809ce5b2..c8981069d 100644 --- a/lib/damask/__init__.py +++ b/lib/damask/__init__.py @@ -13,7 +13,7 @@ from .asciitable import ASCIItable # noqa from .config import Material # noqa from .colormaps import Colormap, Color # noqa -from .orientation import Quaternion, Rodrigues, Symmetry, Orientation # noqa +from .orientation import Quaternion, Symmetry, Orientation # noqa #from .block import Block # only one class from .result import Result # noqa diff --git a/lib/damask/orientation.py b/lib/damask/orientation.py index db2baad92..fbf4dd8a8 100644 --- a/lib/damask/orientation.py +++ b/lib/damask/orientation.py @@ -7,24 +7,6 @@ import math,os import numpy as np -# ****************************************************************************************** -class Rodrigues: - - def __init__(self, vector = np.zeros(3)): - self.vector = vector - - def asQuaternion(self): - norm = np.linalg.norm(self.vector) - halfAngle = np.arctan(norm) - return Quaternion(np.cos(halfAngle),np.sin(halfAngle)*self.vector/norm) - - def asAngleAxis(self): - norm = np.linalg.norm(self.vector) - halfAngle = np.arctan(norm) - return (2.0*halfAngle,self.vector/norm) - - - # ****************************************************************************************** class Quaternion: u""" @@ -178,12 +160,13 @@ class Quaternion: magnitude = __abs__ def __eq__(self,other): - """Equal at e-8 precision""" - return (self-other).magnitude() < 1e-8 or (-self-other).magnitude() < 1e-8 + """Equal (sufficiently close) to each other""" + return math.isclose(( self-other).magnitude(),0.0) \ + or math.isclose((-self-other).magnitude(),0.0) def __ne__(self,other): - """Not equal at e-8 precision""" - return not self.__eq__(self,other) + """Not equal (sufficiently close) to each other""" + return not self.__eq__(other) def __cmp__(self,other): """Linear ordering""" @@ -193,11 +176,6 @@ class Quaternion: def magnitude_squared(self): return self.q ** 2 + np.dot(self.p,self.p) - def identity(self): - self.q = 1. - self.p = np.zeros(3,dtype=float) - return self - def normalize(self): d = self.magnitude() if d > 0.0: @@ -209,13 +187,6 @@ class Quaternion: self.p = -self.p return self - def inverse(self): - d = self.magnitude() - if d > 0.0: - self.conjugate() - self /= d - return self - def homomorph(self): if self.q < 0.0: self.q = -self.q @@ -228,9 +199,6 @@ class Quaternion: def conjugated(self): return self.copy().conjugate() - def inversed(self): - return self.copy().inverse() - def homomorphed(self): return self.copy().homomorph() @@ -259,27 +227,24 @@ class Quaternion: def asAngleAxis(self, degrees = False, flat = False): - if self.q > 1.: - self.normalize() - s = math.sqrt(1. - self.q**2) - x = 2.*self.q**2 - 1. - y = 2.*self.q * s + angle = 2.0*math.acos(self.q) - angle = math.atan2(y,x) - if angle < 0.0: - angle *= -1. - s *= -1. - - if flat: - return np.hstack((np.degrees(angle) if degrees else angle, - np.array([1.0, 0.0, 0.0] if np.abs(angle) < 1e-6 else self.p / s))) + if math.isclose(angle,0.0): + angle = 0.0 + axis = np.array([0.0,0.0,1.0]) + elif math.isclose(self.q,0.0): + angle = math.pi + axis = self.p else: - return (np.degrees(angle) if degrees else angle, - np.array([1.0, 0.0, 0.0] if np.abs(angle) < 1e-6 else self.p / s)) + axis = np.sign(self.q)*self.p/np.linalg.norm(self.p) + + angle = np.degrees(angle) if degrees else angle + + return np.hstack((angle,axis)) if flat else (angle,axis) def asRodrigues(self): - return np.inf*np.ones(3) if self.q == 0.0 else self.p/self.q + return np.inf*np.ones(3) if math.isclose(self.q,0.0) else self.p/self.q def asEulers(self, degrees = False): @@ -290,9 +255,9 @@ class Quaternion: q12 = self.p[0]**2 + self.p[1]**2 chi = np.sqrt(q03*q12) - if abs(chi) < 1e-10 and abs(q12) < 1e-10: + if math.isclose(chi,0.0) and math.isclose(q12,0.0): eulers = np.array([math.atan2(-2*P*self.q*self.p[2],self.q**2-self.p[2]**2),0,0]) - elif abs(chi) < 1e-10 and abs(q03) < 1e-10: + elif math.isclose(chi,0.0) and math.isclose(q03,0.0): eulers = np.array([math.atan2( 2 *self.p[0]*self.p[1],self.p[0]**2-self.p[1]**2),np.pi,0]) else: eulers = np.array([math.atan2((self.p[0]*self.p[2]-P*self.q*self.p[1])/chi,(-P*self.q*self.p[0]-self.p[1]*self.p[2])/chi), @@ -456,11 +421,11 @@ class Symmetry: def __eq__(self, other): - """Equal""" + """Equal to other""" return self.lattice == other.lattice def __neq__(self, other): - """Not equal""" + """Not equal to other""" return not self.__eq__(other) def __cmp__(self,other): @@ -549,7 +514,7 @@ class Symmetry: def inFZ(self,R): """Check whether given Rodrigues vector falls into fundamental zone of own symmetry.""" - if isinstance(R, Quaternion): R = R.asRodrigues() # translate accidentially passed quaternion + if isinstance(R, Quaternion): R = R.asRodrigues() # translate accidentally passed quaternion # fundamental zone in Rodrigues space is point symmetric around origin R = abs(R) if self.lattice == 'cubic': From 5d7f9305920c1e76b59986ceb7eab8c23a1c9d8a Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 17 Dec 2018 20:04:16 +0100 Subject: [PATCH 60/70] math.isclose does not exist for older versions np.isclose does the same --- lib/damask/orientation.py | 14 +++++++------- processing/pre/geom_rotate.py | 2 -- 2 files changed, 7 insertions(+), 9 deletions(-) diff --git a/lib/damask/orientation.py b/lib/damask/orientation.py index fbf4dd8a8..946c97d1e 100644 --- a/lib/damask/orientation.py +++ b/lib/damask/orientation.py @@ -161,8 +161,8 @@ class Quaternion: def __eq__(self,other): """Equal (sufficiently close) to each other""" - return math.isclose(( self-other).magnitude(),0.0) \ - or math.isclose((-self-other).magnitude(),0.0) + 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""" @@ -230,10 +230,10 @@ class Quaternion: angle = 2.0*math.acos(self.q) - if math.isclose(angle,0.0): + if np.isclose(angle,0.0): angle = 0.0 axis = np.array([0.0,0.0,1.0]) - elif math.isclose(self.q,0.0): + elif np.isclose(self.q,0.0): angle = math.pi axis = self.p else: @@ -244,7 +244,7 @@ class Quaternion: return np.hstack((angle,axis)) if flat else (angle,axis) def asRodrigues(self): - return np.inf*np.ones(3) if math.isclose(self.q,0.0) else self.p/self.q + return np.inf*np.ones(3) if np.isclose(self.q,0.0) else self.p/self.q def asEulers(self, degrees = False): @@ -255,9 +255,9 @@ class Quaternion: q12 = self.p[0]**2 + self.p[1]**2 chi = np.sqrt(q03*q12) - if math.isclose(chi,0.0) and math.isclose(q12,0.0): + if np.isclose(chi,0.0) and np.isclose(q12,0.0): eulers = np.array([math.atan2(-2*P*self.q*self.p[2],self.q**2-self.p[2]**2),0,0]) - elif math.isclose(chi,0.0) and math.isclose(q03,0.0): + elif np.isclose(chi,0.0) and np.isclose(q03,0.0): eulers = np.array([math.atan2( 2 *self.p[0]*self.p[1],self.p[0]**2-self.p[1]**2),np.pi,0]) else: eulers = np.array([math.atan2((self.p[0]*self.p[2]-P*self.q*self.p[1])/chi,(-P*self.q*self.p[0]-self.p[1]*self.p[2])/chi), diff --git a/processing/pre/geom_rotate.py b/processing/pre/geom_rotate.py index 939f2b943..eb70f7137 100755 --- a/processing/pre/geom_rotate.py +++ b/processing/pre/geom_rotate.py @@ -60,8 +60,6 @@ eulers = np.array(damask.orientation.Orientation( degrees = options.degrees, ).asEulers(degrees=True)) -damask.util.croak('{} {} {}'.format(*eulers)) - # --- loop over input files ------------------------------------------------------------------------- if filenames == []: filenames = [None] From d794e152491fcfb63de2b768360f310272410070 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Mon, 17 Dec 2018 16:07:51 -0500 Subject: [PATCH 61/70] [skip ci] small polish --- lib/damask/orientation.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/lib/damask/orientation.py b/lib/damask/orientation.py index 946c97d1e..63880a3e6 100644 --- a/lib/damask/orientation.py +++ b/lib/damask/orientation.py @@ -416,7 +416,7 @@ class Symmetry: def __repr__(self): - """Readbable string""" + """Readable string""" return '{}'.format(self.lattice) @@ -502,7 +502,7 @@ class Symmetry: ] return list(map(Quaternion, - np.array(symQuats)[np.atleast_1d(np.array(who)) if who != [] else range(len(symQuats))])) + np.array(symQuats)[np.atleast_1d(np.array(who)) if who != [] else range(len(symQuats))])) def equivalentQuaternions(self, From f34111cf31992d2507333526de9b88d4ef4ec7f2 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Tue, 18 Dec 2018 09:11:44 -0500 Subject: [PATCH 62/70] [skip ci] env/DAMASK.xxx now reports current repository branch of DAMASK_ROOT --- env/DAMASK.csh | 7 ++++++- env/DAMASK.sh | 7 +++++-- env/DAMASK.zsh | 7 +++++-- 3 files changed, 16 insertions(+), 5 deletions(-) diff --git a/env/DAMASK.csh b/env/DAMASK.csh index 2e1dc979c..6b6c58d9d 100644 --- a/env/DAMASK.csh +++ b/env/DAMASK.csh @@ -7,6 +7,11 @@ set DAMASK_ROOT=`python -c "import os,sys; print(os.path.realpath(os.path.expand source $DAMASK_ROOT/CONFIG +# add BRANCH if DAMASK_ROOT is a git repository +cd $DAMASK_ROOT >/dev/null +set BRANCH = `git branch 2>/dev/null| grep -E '^\* ')` +cd - >/dev/null + # if DAMASK_BIN is present if ( $?DAMASK_BIN) then set path = ($DAMASK_BIN $path) @@ -41,7 +46,7 @@ if ( $?prompt ) then echo https://damask.mpie.de echo echo Using environment with ... - echo "DAMASK $DAMASK_ROOT" + echo "DAMASK $DAMASK_ROOT $BRANCH" echo "Spectral Solver $SOLVER" echo "Post Processing $PROCESSING" if ( $?PETSC_DIR) then diff --git a/env/DAMASK.sh b/env/DAMASK.sh index e790ae3cc..bd26a3ebb 100644 --- a/env/DAMASK.sh +++ b/env/DAMASK.sh @@ -30,6 +30,9 @@ set() { source $DAMASK_ROOT/CONFIG unset -f set +# add BRANCH if DAMASK_ROOT is a git repository +cd $DAMASK_ROOT >/dev/null; BRANCH=$(git branch 2>/dev/null| grep -E '^\* '); cd - >/dev/null + # add DAMASK_BIN if present [ "x$DAMASK_BIN" != "x" ] && PATH=$DAMASK_BIN:$PATH @@ -59,7 +62,7 @@ if [ ! -z "$PS1" ]; then echo https://damask.mpie.de echo echo Using environment with ... - echo "DAMASK $DAMASK_ROOT" + echo "DAMASK $DAMASK_ROOT $BRANCH" echo "Spectral Solver $SOLVER" echo "Post Processing $PROCESSING" if [ "x$PETSC_DIR" != "x" ]; then @@ -94,7 +97,7 @@ fi export DAMASK_NUM_THREADS export PYTHONPATH=$DAMASK_ROOT/lib:$PYTHONPATH -for var in BASE STAT SOLVER PROCESSING FREE DAMASK_BIN; do +for var in BASE STAT SOLVER PROCESSING FREE DAMASK_BIN BRANCH; do unset "${var}" done for var in DAMASK MSC; do diff --git a/env/DAMASK.zsh b/env/DAMASK.zsh index dbfde767d..ae34fdfd2 100644 --- a/env/DAMASK.zsh +++ b/env/DAMASK.zsh @@ -21,6 +21,9 @@ set() { source $DAMASK_ROOT/CONFIG unset -f set +# add BRANCH if DAMASK_ROOT is a git repository +cd $DAMASK_ROOT >/dev/null; BRANCH=$(git branch 2>/dev/null| grep -E '^\* '); cd - >/dev/null + # add DAMASK_BIN if present [ "x$DAMASK_BIN != x" ] && PATH=$DAMASK_BIN:$PATH @@ -50,7 +53,7 @@ if [ ! -z "$PS1" ]; then echo https://damask.mpie.de echo echo "Using environment with ..." - echo "DAMASK $DAMASK_ROOT" + echo "DAMASK $DAMASK_ROOT $BRANCH" echo "Spectral Solver $SOLVER" echo "Post Processing $PROCESSING" if [ "x$PETSC_DIR" != "x" ]; then @@ -87,7 +90,7 @@ fi export DAMASK_NUM_THREADS export PYTHONPATH=$DAMASK_ROOT/lib:$PYTHONPATH -for var in BASE STAT SOLVER PROCESSING FREE DAMASK_BIN; do +for var in BASE STAT SOLVER PROCESSING FREE DAMASK_BIN BRANCH; do unset "${var}" done for var in DAMASK MSC; do From 576c4492d49c1456188ca7c2cf290bf821d20b0d Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Tue, 18 Dec 2018 09:24:25 -0500 Subject: [PATCH 63/70] [skip ci] bugfix for DAMASK_BIN==empty check; [ logical ] --> [[ logical ]] to align with bash syntax --- env/DAMASK.zsh | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/env/DAMASK.zsh b/env/DAMASK.zsh index ae34fdfd2..4d5a1e47d 100644 --- a/env/DAMASK.zsh +++ b/env/DAMASK.zsh @@ -25,15 +25,15 @@ unset -f set cd $DAMASK_ROOT >/dev/null; BRANCH=$(git branch 2>/dev/null| grep -E '^\* '); cd - >/dev/null # add DAMASK_BIN if present -[ "x$DAMASK_BIN != x" ] && PATH=$DAMASK_BIN:$PATH +[[ "x$DAMASK_BIN" != "x" ]] && PATH=$DAMASK_BIN:$PATH SOLVER=$(which DAMASK_spectral || true 2>/dev/null) -[ "x$SOLVER" = "x" ] && SOLVER=$(blink 'Not found!') +[[ "x$SOLVER" == "x" ]] && SOLVER=$(blink 'Not found!') PROCESSING=$(which postResults || true 2>/dev/null) -[ "x$PROCESSING" = "x" ] && PROCESSING=$(blink 'Not found!') +[[ "x$PROCESSING" == "x" ]] && PROCESSING=$(blink 'Not found!') -[ "x$DAMASK_NUM_THREADS" = "x" ] && DAMASK_NUM_THREADS=1 +[[ "x$DAMASK_NUM_THREADS" == "x" ]] && DAMASK_NUM_THREADS=1 # currently, there is no information that unlimited causes problems # still, http://software.intel.com/en-us/forums/topic/501500 suggest to fix it @@ -56,13 +56,13 @@ if [ ! -z "$PS1" ]; then echo "DAMASK $DAMASK_ROOT $BRANCH" echo "Spectral Solver $SOLVER" echo "Post Processing $PROCESSING" - if [ "x$PETSC_DIR" != "x" ]; then + if [ "x$PETSC_DIR" != "x" ]; then echo -n "PETSc location " [ -d $PETSC_DIR ] && echo $PETSC_DIR || blink $PETSC_DIR [[ $(canonicalPath "$PETSC_DIR") == $PETSC_DIR ]] \ || echo " ~~> "$(canonicalPath "$PETSC_DIR") fi - [[ "x$PETSC_ARCH" == "x" ]] \ + [[ "x$PETSC_ARCH" == "x" ]] \ || echo "PETSc architecture $PETSC_ARCH" echo -n "MSC.Marc/Mentat " [ -d $MSC_ROOT ] && echo $MSC_ROOT || blink $MSC_ROOT From 7df6dfbf3fdb07bc55e3fc7518d7068a14d22852 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Tue, 18 Dec 2018 15:10:26 -0500 Subject: [PATCH 64/70] renamed "showTable" to "viewTable" to avoid clashing with astropy.table module --- processing/post/{showTable.py => viewTable.py} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename processing/post/{showTable.py => viewTable.py} (100%) diff --git a/processing/post/showTable.py b/processing/post/viewTable.py similarity index 100% rename from processing/post/showTable.py rename to processing/post/viewTable.py From 6a0a0c37542e82f6f745a3ae553427e2b2d5256c Mon Sep 17 00:00:00 2001 From: Test User Date: Wed, 19 Dec 2018 04:20:41 +0100 Subject: [PATCH 65/70] [skip ci] updated version information after successful test of v2.0.2-1187-gcd8ee450 --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 99b9812a2..4fd318a98 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.2-1137-g4dd064a2 +v2.0.2-1187-gcd8ee450 From 10266bad5c1eef069dd75859568eef91b088f931 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 19 Dec 2018 16:57:40 +0100 Subject: [PATCH 66/70] using updated branch in PRIVATE due to changes in the behavior, testing needs a separate branch which is also called orientationClass_with_negative_P --- PRIVATE | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PRIVATE b/PRIVATE index 58137906b..b9a52a85c 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 58137906b84b6cf0e273dfdde623a2986d03f98e +Subproject commit b9a52a85cd65cc27a8e863302bd984abdcad1455 From 6fe17e05f6fb9b8e8e6644ff8670dc079c22ebe4 Mon Sep 17 00:00:00 2001 From: Test User Date: Thu, 20 Dec 2018 07:37:39 +0100 Subject: [PATCH 67/70] [skip ci] updated version information after successful test of v2.0.2-1190-g7df6dfbf --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 99b9812a2..23cc93ef2 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.2-1137-g4dd064a2 +v2.0.2-1190-g7df6dfbf From 25a80723c0c5b477779bb82beec4d47311b1c175 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Thu, 20 Dec 2018 17:09:53 -0500 Subject: [PATCH 68/70] numpy interface changed to be strict about 2D array shape in histogram --- processing/post/binXY.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/processing/post/binXY.py b/processing/post/binXY.py index 726fca39f..ea73d13b9 100755 --- a/processing/post/binXY.py +++ b/processing/post/binXY.py @@ -120,7 +120,7 @@ for name in filenames: delta = minmax[:,1]-minmax[:,0] (grid,xedges,yedges) = np.histogram2d(table.data[:,0],table.data[:,1], bins=options.bins, - range=minmax[0:2,0:2], + range=minmax[:2], weights=None if options.weight is None else table.data[:,2]) if options.normCol: From a1ff380ef4b6cecdf89e17950a0fec624a022627 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Thu, 20 Dec 2018 17:10:43 -0500 Subject: [PATCH 69/70] slightly better aligned debug output --- src/config.f90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/config.f90 b/src/config.f90 index 7ae800f30..441dd953c 100644 --- a/src/config.f90 +++ b/src/config.f90 @@ -142,7 +142,7 @@ subroutine config_init() case (trim(material_partPhase)) call parseFile(phase_name,config_phase,line,fileContent(i+1:)) - if (iand(myDebug,debug_levelBasic) /= 0_pInt) write(6,'(a)') ' Phase parsed'; flush(6) + if (iand(myDebug,debug_levelBasic) /= 0_pInt) write(6,'(a)') ' Phase parsed'; flush(6) case (trim(material_partMicrostructure)) call parseFile(microstructure_name,config_microstructure,line,fileContent(i+1:)) @@ -150,7 +150,7 @@ subroutine config_init() case (trim(material_partCrystallite)) call parseFile(crystallite_name,config_crystallite,line,fileContent(i+1:)) - if (iand(myDebug,debug_levelBasic) /= 0_pInt) write(6,'(a)') ' Crystallite parsed'; flush(6) + if (iand(myDebug,debug_levelBasic) /= 0_pInt) write(6,'(a)') ' Crystallite parsed'; flush(6) case (trim(material_partHomogenization)) call parseFile(homogenization_name,config_homogenization,line,fileContent(i+1:)) @@ -158,7 +158,7 @@ subroutine config_init() case (trim(material_partTexture)) call parseFile(texture_name,config_texture,line,fileContent(i+1:)) - if (iand(myDebug,debug_levelBasic) /= 0_pInt) write(6,'(a)') ' Texture parsed'; flush(6) + if (iand(myDebug,debug_levelBasic) /= 0_pInt) write(6,'(a)') ' Texture parsed'; flush(6) end select From 9bdd4d1d684fc026832127a794aab5547da799a1 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Thu, 20 Dec 2018 17:12:05 -0500 Subject: [PATCH 70/70] added line to stop config parser from trying to interpret rest of geom description---huge speedup!! --- processing/pre/geom_fromVoronoiTessellation.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/processing/pre/geom_fromVoronoiTessellation.py b/processing/pre/geom_fromVoronoiTessellation.py index 359acb366..5610c939a 100755 --- a/processing/pre/geom_fromVoronoiTessellation.py +++ b/processing/pre/geom_fromVoronoiTessellation.py @@ -323,12 +323,13 @@ for name in filenames: ] if hasEulers: config_header += [''] + theAxes = [] if options.axes is None else ['axes\t{} {} {}'.format(*options.axes)] for ID in grainIDs: eulerID = np.nonzero(grains == ID)[0][0] # find first occurrence of this grain id config_header += ['[Grain{}]'.format(str(ID).zfill(formatwidth)), '(gauss)\tphi1 {:g}\tPhi {:g}\tphi2 {:g}\tscatter 0.0\tfraction 1.0'.format(*eulers[eulerID]) - ] - if options.axes is not None: config_header.append('axes\t{} {} {}'.format(*options.axes)) + ] + theAxes + config_header += [''] table.labels_clear() table.info_clear()