Merge branch 'misc-improvements' into Marc-2019.1

This commit is contained in:
Martin Diehl 2020-07-16 09:28:11 +02:00
commit 190c9a1b0d
42 changed files with 1105 additions and 1055 deletions

View File

@ -1 +1 @@
v2.0.3-2872-g0a2d3046
v2.0.3-2881-gc07efe84

2
env/DAMASK.csh vendored
View File

@ -9,7 +9,7 @@ source $ENV_ROOT/CONFIG
set path = ($DAMASK_ROOT/bin $path)
set SOLVER=`which DAMASK_spectral`
set SOLVER=`which DAMASK_grid`
if ( "x$DAMASK_NUM_THREADS" == "x" ) then
set DAMASK_NUM_THREADS=1
endif

2
env/DAMASK.sh vendored
View File

@ -35,7 +35,7 @@ cd $DAMASK_ROOT >/dev/null; BRANCH=$(git branch 2>/dev/null| grep -E '^\* '); cd
PATH=${DAMASK_ROOT}/bin:$PATH
SOLVER=$(type -p DAMASK_spectral || true 2>/dev/null)
SOLVER=$(type -p DAMASK_grid || true 2>/dev/null)
[ "x$SOLVER" == "x" ] && SOLVER=$(blink 'Not found!')
[ "x$DAMASK_NUM_THREADS" == "x" ] && DAMASK_NUM_THREADS=1

2
env/DAMASK.zsh vendored
View File

@ -27,7 +27,7 @@ cd $DAMASK_ROOT >/dev/null; BRANCH=$(git branch 2>/dev/null| grep -E '^\* '); cd
PATH=${DAMASK_ROOT}/bin:$PATH
SOLVER=$(which DAMASK_spectral || true 2>/dev/null)
SOLVER=$(which DAMASK_grid || true 2>/dev/null)
[[ "x$SOLVER" == "x" ]] && SOLVER=$(blink 'Not found!')
[[ "x$DAMASK_NUM_THREADS" == "x" ]] && DAMASK_NUM_THREADS=1

View File

@ -1,2 +0,0 @@
itmin 4
itmax 40

View File

@ -0,0 +1,3 @@
grid:
itmin: 4
itmax: 40

View File

@ -49,4 +49,5 @@ echo "program: $program"
exit 1
}
/bin/rm $userobj
/bin/rm $DIRJOB/*.mod
/bin/rm $DIRJOB/*.mod
/bin/rm $DIRJOB/*.smod

View File

@ -49,4 +49,5 @@ echo "program: $program"
exit 1
}
/bin/rm $userobj
/bin/rm $DIRJOB/*.mod
/bin/rm $DIRJOB/*.mod
/bin/rm $DIRJOB/*.smod

View File

@ -49,4 +49,5 @@ echo "program: $program"
exit 1
}
/bin/rm $userobj
/bin/rm $DIRJOB/*.mod
/bin/rm $DIRJOB/*.mod
/bin/rm $DIRJOB/*.smod

View File

@ -3267,6 +3267,7 @@ else
fi
/bin/rm $userobj 2>/dev/null
/bin/rm $DIRJOB/*.mod 2>/dev/null
/bin/rm $DIRJOB/*.smod 2>/dev/null
#
# run marc
@ -3599,6 +3600,7 @@ else # if test $link
fi # if test $link
/bin/rm $userobj 2>/dev/null
/bin/rm $DIRJOB/*.mod 2>/dev/null
/bin/rm $DIRJOB/*.smod 2>/dev/null
#
# run marc
@ -3933,6 +3935,7 @@ else # if test $link
fi # if test $link
/bin/rm $userobj 2>/dev/null
/bin/rm $DIRJOB/*.mod 2>/dev/null
/bin/rm $DIRJOB/*.smod 2>/dev/null
# done if no job id given
if test -z "$jid"
then

View File

@ -3267,6 +3267,7 @@ else
fi
/bin/rm $userobj 2>/dev/null
/bin/rm $DIRJOB/*.mod 2>/dev/null
/bin/rm $DIRJOB/*.smod 2>/dev/null
#
# run marc
@ -3599,6 +3600,7 @@ else # if test $link
fi # if test $link
/bin/rm $userobj 2>/dev/null
/bin/rm $DIRJOB/*.mod 2>/dev/null
/bin/rm $DIRJOB/*.smod 2>/dev/null
#
# run marc
@ -3933,6 +3935,7 @@ else # if test $link
fi # if test $link
/bin/rm $userobj 2>/dev/null
/bin/rm $DIRJOB/*.mod 2>/dev/null
/bin/rm $DIRJOB/*.smod 2>/dev/null
# done if no job id given
if test -z "$jid"
then

View File

@ -3267,6 +3267,7 @@ else
fi
/bin/rm $userobj 2>/dev/null
/bin/rm $DIRJOB/*.mod 2>/dev/null
/bin/rm $DIRJOB/*.smod 2>/dev/null
#
# run marc
@ -3599,6 +3600,7 @@ else # if test $link
fi # if test $link
/bin/rm $userobj 2>/dev/null
/bin/rm $DIRJOB/*.mod 2>/dev/null
/bin/rm $DIRJOB/*.smod 2>/dev/null
#
# run marc
@ -3933,6 +3935,7 @@ else # if test $link
fi # if test $link
/bin/rm $userobj 2>/dev/null
/bin/rm $DIRJOB/*.mod 2>/dev/null
/bin/rm $DIRJOB/*.smod 2>/dev/null
# done if no job id given
if test -z "$jid"
then

View File

@ -49,4 +49,5 @@ echo "program: $program"
exit 1
}
/bin/rm $userobj
/bin/rm $DIRJOB/*.mod
/bin/rm $DIRJOB/*.mod
/bin/rm $DIRJOB/*.smod

View File

@ -49,4 +49,5 @@ echo "program: $program"
exit 1
}
/bin/rm $userobj
/bin/rm $DIRJOB/*.mod
/bin/rm $DIRJOB/*.mod
/bin/rm $DIRJOB/*.smod

View File

@ -49,4 +49,5 @@ echo "program: $program"
exit 1
}
/bin/rm $userobj
/bin/rm $DIRJOB/*.mod
/bin/rm $DIRJOB/*.mod
/bin/rm $DIRJOB/*.smod

View File

@ -3267,6 +3267,7 @@ else
fi
/bin/rm $userobj 2>/dev/null
/bin/rm $DIRJOB/*.mod 2>/dev/null
/bin/rm $DIRJOB/*.smod 2>/dev/null
#
# run marc
@ -3599,6 +3600,7 @@ else # if test $link
fi # if test $link
/bin/rm $userobj 2>/dev/null
/bin/rm $DIRJOB/*.mod 2>/dev/null
/bin/rm $DIRJOB/*.smod 2>/dev/null
#
# run marc
@ -3933,6 +3935,7 @@ else # if test $link
fi # if test $link
/bin/rm $userobj 2>/dev/null
/bin/rm $DIRJOB/*.mod 2>/dev/null
/bin/rm $DIRJOB/*.smod 2>/dev/null
# done if no job id given
if test -z "$jid"
then

View File

@ -3267,6 +3267,7 @@ else
fi
/bin/rm $userobj 2>/dev/null
/bin/rm $DIRJOB/*.mod 2>/dev/null
/bin/rm $DIRJOB/*.smod 2>/dev/null
#
# run marc
@ -3599,6 +3600,7 @@ else # if test $link
fi # if test $link
/bin/rm $userobj 2>/dev/null
/bin/rm $DIRJOB/*.mod 2>/dev/null
/bin/rm $DIRJOB/*.smod 2>/dev/null
#
# run marc
@ -3933,6 +3935,7 @@ else # if test $link
fi # if test $link
/bin/rm $userobj 2>/dev/null
/bin/rm $DIRJOB/*.mod 2>/dev/null
/bin/rm $DIRJOB/*.smod 2>/dev/null
# done if no job id given
if test -z "$jid"
then

View File

@ -3267,6 +3267,7 @@ else
fi
/bin/rm $userobj 2>/dev/null
/bin/rm $DIRJOB/*.mod 2>/dev/null
/bin/rm $DIRJOB/*.smod 2>/dev/null
#
# run marc
@ -3599,6 +3600,7 @@ else # if test $link
fi # if test $link
/bin/rm $userobj 2>/dev/null
/bin/rm $DIRJOB/*.mod 2>/dev/null
/bin/rm $DIRJOB/*.smod 2>/dev/null
#
# run marc
@ -3933,6 +3935,7 @@ else # if test $link
fi # if test $link
/bin/rm $userobj 2>/dev/null
/bin/rm $DIRJOB/*.mod 2>/dev/null
/bin/rm $DIRJOB/*.smod 2>/dev/null
# done if no job id given
if test -z "$jid"
then

View File

@ -49,4 +49,5 @@ echo "program: $program"
exit 1
}
/bin/rm $userobj
/bin/rm $DIRJOB/*.mod
/bin/rm $DIRJOB/*.mod
/bin/rm $DIRJOB/*.smod

View File

@ -49,4 +49,5 @@ echo "program: $program"
exit 1
}
/bin/rm $userobj
/bin/rm $DIRJOB/*.mod
/bin/rm $DIRJOB/*.mod
/bin/rm $DIRJOB/*.smod

View File

@ -49,4 +49,5 @@ echo "program: $program"
exit 1
}
/bin/rm $userobj
/bin/rm $DIRJOB/*.mod
/bin/rm $DIRJOB/*.mod
/bin/rm $DIRJOB/*.smod

View File

@ -3241,6 +3241,7 @@ else
fi
/bin/rm $userobj 2>/dev/null
/bin/rm $DIRJOB/*.mod 2>/dev/null
/bin/rm $DIRJOB/*.smod 2>/dev/null
#
# run marc
@ -3573,6 +3574,7 @@ else # if test $link
fi # if test $link
/bin/rm $userobj 2>/dev/null
/bin/rm $DIRJOB/*.mod 2>/dev/null
/bin/rm $DIRJOB/*.smod 2>/dev/null
#
# run marc
@ -3907,6 +3909,7 @@ else # if test $link
fi # if test $link
/bin/rm $userobj 2>/dev/null
/bin/rm $DIRJOB/*.mod 2>/dev/null
/bin/rm $DIRJOB/*.smod 2>/dev/null
# done if no job id given
if test -z "$jid"
then

View File

@ -3241,6 +3241,7 @@ else
fi
/bin/rm $userobj 2>/dev/null
/bin/rm $DIRJOB/*.mod 2>/dev/null
/bin/rm $DIRJOB/*.smod 2>/dev/null
#
# run marc
@ -3573,6 +3574,7 @@ else # if test $link
fi # if test $link
/bin/rm $userobj 2>/dev/null
/bin/rm $DIRJOB/*.mod 2>/dev/null
/bin/rm $DIRJOB/*.smod 2>/dev/null
#
# run marc
@ -3907,6 +3909,7 @@ else # if test $link
fi # if test $link
/bin/rm $userobj 2>/dev/null
/bin/rm $DIRJOB/*.mod 2>/dev/null
/bin/rm $DIRJOB/*.smod 2>/dev/null
# done if no job id given
if test -z "$jid"
then

View File

@ -3241,6 +3241,7 @@ else
fi
/bin/rm $userobj 2>/dev/null
/bin/rm $DIRJOB/*.mod 2>/dev/null
/bin/rm $DIRJOB/*.smod 2>/dev/null
#
# run marc
@ -3573,6 +3574,7 @@ else # if test $link
fi # if test $link
/bin/rm $userobj 2>/dev/null
/bin/rm $DIRJOB/*.mod 2>/dev/null
/bin/rm $DIRJOB/*.smod 2>/dev/null
#
# run marc
@ -3907,6 +3909,7 @@ else # if test $link
fi # if test $link
/bin/rm $userobj 2>/dev/null
/bin/rm $DIRJOB/*.mod 2>/dev/null
/bin/rm $DIRJOB/*.smod 2>/dev/null
# done if no job id given
if test -z "$jid"
then

View File

@ -5,9 +5,11 @@ import re as _re
name = 'damask'
with open(_Path(__file__).parent/_Path('VERSION')) as _f:
version = _re.sub(r'^v','',_f.readline().strip())
__version__ = version
# make classes directly accessible as damask.Class
from ._environment import Environment # noqa
from ._environment import Environment as _ # noqa
environment = _()
from ._table import Table # noqa
from ._vtk import VTK # noqa
from ._colormap import Colormap # noqa
@ -19,6 +21,7 @@ from ._geom import Geom # noqa
from . import solver # noqa
# deprecated
Environment = _
from ._asciitable import ASCIItable # noqa
from ._test import Test # noqa
from .config import Material # noqa

View File

@ -1,8 +1,11 @@
import os
import json
import functools
import numpy as np
import matplotlib as mpl
if os.name == 'posix' and 'DISPLAY' not in os.environ:
mpl.use('Agg')
import matplotlib.pyplot as plt
from matplotlib import cm

View File

@ -6,9 +6,9 @@ from functools import partial
import numpy as np
from scipy import ndimage,spatial
import damask
from . import VTK
from . import util
from . import Environment
from . import grid_filters
@ -362,7 +362,7 @@ class Geom:
seeds_p = seeds
coords = grid_filters.cell_coord0(grid,size).reshape(-1,3)
pool = multiprocessing.Pool(processes = int(Environment().options['DAMASK_NUM_THREADS']))
pool = multiprocessing.Pool(processes = int(damask.environment.options['DAMASK_NUM_THREADS']))
result = pool.map_async(partial(Geom._find_closest_seed,seeds_p,weights_p), [coord for coord in coords])
pool.close()
pool.join()

View File

@ -18,7 +18,6 @@ from . import VTK
from . import Table
from . import Rotation
from . import Orientation
from . import Environment
from . import grid_filters
from . import mechanics
from . import util
@ -1060,7 +1059,7 @@ class Result:
Arguments parsed to func.
"""
num_threads = Environment().options['DAMASK_NUM_THREADS']
num_threads = damask.environment.options['DAMASK_NUM_THREADS']
pool = mp.Pool(int(num_threads) if num_threads is not None else None)
lock = mp.Manager().Lock()

View File

@ -805,7 +805,7 @@ class Rotation:
"""Rotation matrix to Bunge-Euler angles."""
with np.errstate(invalid='ignore',divide='ignore'):
zeta = 1.0/np.sqrt(1.0-om[...,2,2:3]**2)
eu = np.where(np.isclose(np.abs(om[...,2,2:3]),1.0,1e-9),
eu = np.where(np.isclose(np.abs(om[...,2,2:3]),1.0,0.0),
np.block([np.arctan2(om[...,0,1:2],om[...,0,0:1]),
np.pi*0.5*(1-om[...,2,2:3]),
np.zeros(om.shape[:-2]+(1,)),

View File

@ -180,7 +180,7 @@ class Test:
def fileInRoot(self,dir,file):
"""Path to a file in the root directory of DAMASK."""
return os.path.join(damask.Environment().rootDir(),dir,file)
return str(damask.environment.root_dir/dir/file)
def fileInReference(self,file):

View File

@ -9,7 +9,6 @@ from vtk.util.numpy_support import numpy_to_vtkIdTypeArray as np_to_vtkIdTypeArr
import damask
from . import Table
from . import Environment
class VTK:
@ -258,7 +257,7 @@ class VTK:
ren.AddActor(actor)
ren.SetBackground(0.2,0.2,0.2)
window.SetSize(Environment().screen_size[0],Environment().screen_size[1])
window.SetSize(damask.environment.screen_size[0],damask.environment.screen_size[1])
iren = vtk.vtkRenderWindowInteractor()
iren.SetRenderWindow(window)

View File

@ -3,12 +3,12 @@ import shlex
import string
from pathlib import Path
from .._environment import Environment
import damask
class Marc:
"""Wrapper to run DAMASK with MSCMarc."""
def __init__(self,version=Environment().options['MARC_VERSION']):
def __init__(self,version=damask.environment.options['MARC_VERSION']):
"""
Create a Marc solver object.
@ -24,7 +24,7 @@ class Marc:
@property
def library_path(self):
path_MSC = Environment().options['MSC_ROOT']
path_MSC = damask.environment.options['MSC_ROOT']
path_lib = Path(f'{path_MSC}/mentat{self.version}/shlib/linux64')
return path_lib if path_lib.is_dir() else None
@ -33,7 +33,7 @@ class Marc:
@property
def tools_path(self):
path_MSC = Environment().options['MSC_ROOT']
path_MSC = damask.environment.options['MSC_ROOT']
path_tools = Path(f'{path_MSC}/marc{self.version}/tools')
return path_tools if path_tools.is_dir() else None
@ -49,9 +49,7 @@ class Marc:
):
env = Environment()
usersub = env.root_dir/Path(f'src/DAMASK_marc{self.version}').with_suffix('.f90' if compile else '.marc')
usersub = damask.environment.root_dir/Path(f'src/DAMASK_marc{self.version}').with_suffix('.f90' if compile else '.marc')
if not usersub.is_file():
raise FileNotFoundError("DAMASK4Marc ({}) '{}' not found".format(('source' if compile else 'binary'),usersub))

View File

@ -1,8 +1,8 @@
import setuptools
import os
from pathlib import Path
import re
with open(os.path.join(os.path.dirname(__file__),'damask/VERSION')) as f:
with open(Path(__file__).parent/'damask/VERSION') as f:
version = re.sub(r'(-([^-]*)).*$',r'.\2',re.sub(r'^v(\d+\.\d+(\.\d+)?)',r'\1',f.readline().strip()))
setuptools.setup(

View File

@ -150,7 +150,7 @@ def om2qu(a):
def om2eu(om):
"""Rotation matrix to Bunge-Euler angles."""
if not np.isclose(np.abs(om[2,2]),1.0,1.e-9):
if not np.isclose(np.abs(om[2,2]),1.0,0.0):
zeta = 1.0/np.sqrt(1.0-om[2,2]**2)
eu = np.array([np.arctan2(om[2,0]*zeta,-om[2,1]*zeta),
np.arccos(om[2,2]),

View File

@ -248,7 +248,7 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS
endif validCalculation
if (debugCPFEM%extensive &
.and. (debugCPFEM%element == elCP .and. debugCPFEM%ip == ip) .or. .not. debugCPFEM%selective) then
.and. ((debugCPFEM%element == elCP .and. debugCPFEM%ip == ip) .or. .not. debugCPFEM%selective)) then
write(6,'(a,i8,1x,i2,/,12x,6(f10.3,1x)/)') &
'<< CPFEM >> stress/MPa at elFE ip ', elFE, ip, CPFEM_cs(1:6,ip,elCP)*1.0e-6_pReal
write(6,'(a,i8,1x,i2,/,6(12x,6(f10.3,1x)/))') &

View File

@ -529,6 +529,8 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg)
msg = 'unknown material parameter:'
case (211)
msg = 'material parameter out of bounds:'
case (212)
msg = 'nonlocal model not supported'
!--------------------------------------------------------------------------------------------------
! numerics error messages

View File

@ -15,7 +15,8 @@ module constitutive
use results
use lattice
use discretization
use geometry_plastic_nonlocal
use geometry_plastic_nonlocal, only: &
geometry_plastic_nonlocal_disable
use source_thermal_dissipation
use source_thermal_externalheat
use source_damage_isoBrittle

View File

@ -386,7 +386,11 @@ module subroutine plastic_nonlocal_init
call material_allocateState(plasticState(p),NipcMyPhase,sizeState,sizeDotState,sizeDeltaState)
plasticState(p)%nonlocal = config%KeyExists('/nonlocal/')
plasticState(p)%nonlocal = config%KeyExists('/nonlocal/')
if(plasticState(p)%nonlocal .and. .not. allocated(IPneighborhood)) &
call IO_error(212,ext_msg='IPneighborhood does not exist')
plasticState(p)%offsetDeltaState = 0 ! ToDo: state structure does not follow convention
st0%rho => plasticState(p)%state0 (0*prm%sum_N_sl+1:10*prm%sum_N_sl,:)

View File

@ -118,6 +118,7 @@ subroutine discretization_grid_init(restart)
call results_addAttribute('origin',origin, 'geometry')
call results_closeJobFile
endif
!--------------------------------------------------------------------------------------------------
! geometry information required by the nonlocal CP model
call geometry_plastic_nonlocal_setIPvolume(reshape([(product(mySize/real(myGrid,pReal)),j=1,product(myGrid))], &

View File

@ -119,6 +119,7 @@ subroutine discretization_marc_init
unscaledNormals = IPareaNormal(elem,nElems,connectivity_cell,node0_cell)
call geometry_plastic_nonlocal_setIParea(norm2(unscaledNormals,1))
call geometry_plastic_nonlocal_setIPareaNormal(unscaledNormals/spread(norm2(unscaledNormals,1),1,3))
!call geometry_plastic_nonlocal_setIPneighborhood ToDo: Support nonlocal
call geometry_plastic_nonlocal_results
end subroutine discretization_marc_init
@ -492,8 +493,8 @@ subroutine inputRead_mapNodes(FEM2DAMASK, &
chunkPos = IO_stringPos(fileContent(l))
if(chunkPos(1) < 1) cycle
if(IO_lc(IO_stringValue(fileContent(l),chunkPos,1)) == 'coordinates') then
chunkPos = [1,1,10]
do i = 1,nNodes
chunkPos = IO_stringPos(fileContent(l+1+i))
map_unsorted(:,i) = [IO_intValue(fileContent(l+1+i),chunkPos,1),i]
enddo
exit
@ -528,8 +529,8 @@ subroutine inputRead_elemNodes(nodes, &
chunkPos = IO_stringPos(fileContent(l))
if(chunkPos(1) < 1) cycle
if(IO_lc(IO_stringValue(fileContent(l),chunkPos,1)) == 'coordinates') then
chunkPos = [4,1,10,11,30,31,50,51,70]
do i=1,nNode
chunkPos = IO_stringPos(fileContent(l+1+i))
m = mesh_FEM2DAMASK_node(IO_intValue(fileContent(l+1+i),chunkPos,1))
do j = 1,3
nodes(j,m) = mesh_unitlength * IO_floatValue(fileContent(l+1+i),chunkPos,j+1)
@ -732,7 +733,7 @@ end subroutine inputRead_microstructureAndHomogenization
!--------------------------------------------------------------------------------------------------
!> @brief Calculates cell node coordinates from element node coordinates
!> @brief Determine cell connectivity and definition of cell nodes
!--------------------------------------------------------------------------------------------------
subroutine buildCells(connectivity_cell,cellNodeDefinition, &
elem,connectivity_elem)
@ -746,7 +747,8 @@ subroutine buildCells(connectivity_cell,cellNodeDefinition, &
integer,dimension(:), allocatable :: candidates_local
integer,dimension(:,:), allocatable :: parentsAndWeights,candidates_global
integer :: e, n, c, p, s,i,m,j,nParentNodes,nCellNode,Nelem,candidateID
integer :: e,n,c,p,s,i,m,j,&
nParentNodes,nCellNode,Nelem,candidateID
Nelem = size(connectivity_elem,2)
@ -760,9 +762,7 @@ subroutine buildCells(connectivity_cell,cellNodeDefinition, &
do e = 1, Nelem
do c = 1, elem%NcellNodes
realNode: if (count(elem%cellNodeParentNodeWeights(:,c) /= 0) == 1) then
where(connectivity_cell(:,:,e) == -c)
connectivity_cell(:,:,e) = connectivity_elem(c,e)
end where
where(connectivity_cell(:,:,e) == -c) connectivity_cell(:,:,e) = connectivity_elem(c,e)
endif realNode
enddo
enddo
@ -889,7 +889,7 @@ end subroutine buildCells
!--------------------------------------------------------------------------------------------------
!> @brief Calculates cell node coordinates from element node coordinates
!> @brief Calculate cell node coordinates from element node coordinates
!--------------------------------------------------------------------------------------------------
subroutine buildCellNodes(node_cell, &
definition,node_elem)
@ -919,7 +919,7 @@ end subroutine buildCellNodes
!--------------------------------------------------------------------------------------------------
!> @brief Calculates IP coordinates as center of cell
!> @brief Calculate IP coordinates as center of cell
!--------------------------------------------------------------------------------------------------
subroutine buildIPcoordinates(IPcoordinates, &
connectivity_cell,node_cell)

View File

@ -464,6 +464,8 @@ subroutine selfTest
real(pReal), dimension(4) :: qu
type(quaternion) :: q, q_2
if(dNeq(abs(P),1.0_pReal)) call IO_error(0,ext_msg='P not in {-1,+1}')
call random_number(qu)
qu = (qu-0.5_pReal) * 2.0_pReal
q = quaternion(qu)

View File

@ -574,6 +574,7 @@ end function om2qu
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief orientation matrix to Euler angles
!> @details Two step check for special cases to avoid invalid operations (not needed for python)
!---------------------------------------------------------------------------------------------------
pure function om2eu(om) result(eu)
@ -581,16 +582,16 @@ pure function om2eu(om) result(eu)
real(pReal), dimension(3) :: eu
real(pReal) :: zeta
if (abs(om(3,3)) < 1.0_pReal) then
zeta = 1.0_pReal/sqrt(1.0_pReal-om(3,3)**2.0_pReal)
if (dNeq(abs(om(3,3)),1.0_pReal,1.e-8_pReal)) then
zeta = 1.0_pReal/sqrt(math_clip(1.0_pReal-om(3,3)**2.0_pReal,1e-64_pReal,1.0_pReal))
eu = [atan2(om(3,1)*zeta,-om(3,2)*zeta), &
acos(om(3,3)), &
acos(math_clip(om(3,3),-1.0_pReal,1.0_pReal)), &
atan2(om(1,3)*zeta, om(2,3)*zeta)]
else
eu = [atan2(om(1,2),om(1,1)), 0.5_pReal*PI*(1.0_pReal-om(3,3)),0.0_pReal ]
end if
where(eu<0.0_pReal) eu = mod(eu+2.0_pReal*PI,[2.0_pReal*PI,PI,2.0_pReal*PI])
where(abs(eu) < 1.e-8_pReal) eu = 0.0_pReal
where(eu<0.0_pReal) eu = mod(eu+2.0_pReal*PI,[2.0_pReal*PI,PI,2.0_pReal*PI])
end function om2eu