223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

Merge remote-tracking branch 'origin/development' into typehints_table

This commit is contained in:
Martin Diehl 2022-01-23 12:45:06 +01:00
parent 4d3e0e6a51 efd9b37235
commit 803c85c2ef
26 changed files with 550 additions and 380 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

22
.gitlab-ci.yml
View File

@ -46,7 +46,7 @@ variables:
# ++++++++++++ PETSc ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
PETSC_GNU: "Libraries/PETSc/3.16.1/GNU-10-OpenMPI-4.1.1"
PETSC_INTELLLVM: "Libraries/PETSc/3.16.3/oneAPI-2022.0.1-IntelMPI-2021.5.0"
PETSC_INTEL: "Libraries/PETSc/3.16.2/Intel-2022.0.1-IntelMPI-2021.5.0"
PETSC_INTEL: "Libraries/PETSc/3.16.3/Intel-2022.0.1-IntelMPI-2021.5.0"
# ++++++++++++ MSC Marc +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
MSC: "FEM/MSC/2021.3.1"
IntelMarc: "Compiler/Intel/19.1.2 Libraries/IMKL/2020"
@ -79,63 +79,63 @@ mypy:
###################################################################################################
test_grid_GNU:
grid_GNU:
stage: compile
script:
- module load ${COMPILER_GNU} ${MPI_GNU} ${PETSC_GNU}
- cd PRIVATE/testing/pytest
- pytest -k 'compile and grid' --basetemp ${TESTROOT}/compile_grid_GNU
test_mesh_GNU:
mesh_GNU:
stage: compile
script:
- module load ${COMPILER_GNU} ${MPI_GNU} ${PETSC_GNU}
- cd PRIVATE/testing/pytest
- pytest -k 'compile and mesh' --basetemp ${TESTROOT}/compile_mesh_GNU
test_grid_GNU-64bit:
grid_GNU-64bit:
stage: compile
script:
- module load Compiler/GNU/10 Libraries/PETSc/3.16.2/64bit
- cd PRIVATE/testing/pytest
- pytest -k 'compile and grid' --basetemp ${TESTROOT}/compile_grid_GNU-64bit
test_mesh_GNU-64bit:
mesh_GNU-64bit:
stage: compile
script:
- module load Compiler/GNU/10 Libraries/PETSc/3.16.2/64bit
- cd PRIVATE/testing/pytest
- pytest -k 'compile and mesh' --basetemp ${TESTROOT}/compile_mesh_GNU-64bit
test_grid_IntelLLVM:
grid_IntelLLVM:
stage: compile
script:
- module load ${COMPILER_INTELLLVM} ${MPI_INTELLLVM} ${PETSC_INTELLLVM}
- cd PRIVATE/testing/pytest
- pytest -k 'compile and grid' --basetemp ${TESTROOT}/compile_grid_IntelLLVM
test_mesh_IntelLLVM:
mesh_IntelLLVM:
stage: compile
script:
- module load ${COMPILER_INTELLLVM} ${MPI_INTELLLVM} ${PETSC_INTELLLVM}
- cd PRIVATE/testing/pytest
- pytest -k 'compile and mesh' --basetemp ${TESTROOT}/compile_mesh_IntelLLVM
test_grid_Intel:
grid_Intel:
stage: compile
script:
- module load ${COMPILER_INTEL} ${MPI_INTEL} ${PETSC_INTEL}
- cd PRIVATE/testing/pytest
- pytest -k 'compile and grid' --basetemp ${TESTROOT}/compile_grid_Intel
test_mesh_Intel:
mesh_Intel:
stage: compile
script:
- module load ${COMPILER_INTEL} ${MPI_INTEL} ${PETSC_INTEL}
- cd PRIVATE/testing/pytest
- pytest -k 'compile and mesh' --basetemp ${TESTROOT}/compile_mesh_Intel
test_Marc:
Marc_Intel:
stage: compile
script:
- module load $IntelMarc $HDF5Marc $MSC
@ -158,7 +158,7 @@ setup_mesh:
- cmake -DDAMASK_SOLVER=MESH -DCMAKE_INSTALL_PREFIX=${TESTROOT} ${CI_PROJECT_DIR}
- make -j2 all install
compile_Marc:
setup_Marc:
stage: compile
script:
- module load $IntelMarc $HDF5Marc $MSC

2
PRIVATE

@ -1 +1 @@
Subproject commit b898a8b5552bd9d1c555edc3d8134564dd32fe53
Subproject commit d570b4a1fac5b9b99d026d3ff9bf593615e22ce5

12
examples/config/phase/mechanical/plastic/phenopowerlaw_Mg.yaml
View File

@ -6,7 +6,7 @@ references:
output: [xi_sl, xi_tw]
N_sl: [3, 3, 6, 0, 6] # basal, prism, -, 1. pyr<a>, -, 2. pyr<c+a>
N_sl: [3, 3, 6, 0, 6] # basal, prism, 1. pyr<a>, -, 2. pyr<c+a>
N_tw: [6, 0, 6] # tension, -, compression
xi_0_sl: [10.e+6, 55.e+6, 60.e+6, 0., 60.e+6]
@ -23,9 +23,13 @@ f_sat_sl-tw: 10.0
h_0_sl-sl: 500.0e+6
h_0_tw-tw: 50.0e+6
h_0_tw-sl: 150.0e+6
h_sl-sl: [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
1.0, 1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0,
1.0, 1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, +1.0, 1.0] # unused entries are indicated by -1.0
h_sl-sl: [+1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, -1.0,
-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0,
-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0,
-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 1.0, 1.0, 1.0, 1.0,
+1.0, 1.0, -1.0, -1.0, -1.0, -1.0, 1.0, 1.0, 1.0, 1.0,
-1.0, -1.0, -1.0, -1.0, 1.0, 1.0, 1.0, 1.0, 1.0] # unused entries are indicated by -1.0
h_tw-tw: [+1.0, 1.0, -1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0,
-1.0, 1.0] # unused entries are indicated by -1.0
h_tw-sl: [1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0,

9
examples/config/phase/mechanical/plastic/phenopowerlaw_Ti.yaml
View File

@ -8,7 +8,7 @@ references:
https://doi.org/10.1016/j.actamat.2017.05.015
output: [gamma_sl]
N_sl: [3, 3, 0, 12] # basal, 1. prism, -, 1. pyr<c+a>
N_sl: [3, 3, 0, 12] # basal, prism, -, 1. pyr<c+a>
n_sl: 20
a_sl: 2.0
dot_gamma_0_sl: 0.001
@ -20,5 +20,8 @@ xi_inf_sl: [568.e+6, 150.e+7, 0.0, 3420.e+6]
# L. Wang et al. :
# xi_0_sl: [127.e+6, 96.e+6, 0.0, 240.e+6]
h_sl-sl: [+1.0, 1.0, 1.0, 1.0, 1.0, 1.0, -1.0, -1.0, -1.0, -1.0,
-1.0, -1.0, +1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0] # unused entries are indicated by -1.0
h_sl-sl: [+1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, -1.0, -1.0,
-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 1.0,
+1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+1.0, 1.0, 1.0, 1.0, 1.0, 1.0] # unused entries are indicated by -1.0

2
python/damask/VERSION
View File

@ -1 +1 @@
v3.0.0-alpha5-389-ga000e477c
v3.0.0-alpha5-475-g160eb1c60

2
python/damask/__init__.py
View File

@ -21,8 +21,8 @@ from ._rotation import Rotation # noqa
from ._crystal import Crystal # noqa
from ._orientation import Orientation # noqa
from ._table import Table # noqa
from ._vtk import VTK # noqa
from ._colormap import Colormap # noqa
from ._vtk import VTK # noqa
from ._config import Config # noqa
from ._configmaterial import ConfigMaterial # noqa
from ._grid import Grid # noqa

2
python/damask/_colormap.py
View File

@ -250,7 +250,7 @@ class Colormap(mpl.colors.ListedColormap):
np.logical_or (np.isnan(field), field == gap)) # mask NaN (and gap if present)
l,r = (field[mask].min(),field[mask].max()) if bounds is None else \
np.array(bounds,float)[:2]
(bounds[0],bounds[1])
delta,avg = r-l,0.5*abs(r+l)

71
python/damask/_grid.py
View File

@ -4,7 +4,7 @@ import warnings
import multiprocessing as mp
from functools import partial
import typing
from typing import Union, Optional, TextIO, List, Sequence
from typing import Union, Optional, TextIO, List, Sequence, Literal
from pathlib import Path
import numpy as np
@ -70,7 +70,7 @@ class Grid:
])
def __copy__(self) -> "Grid":
def __copy__(self) -> 'Grid':
"""Create deep copy."""
return copy.deepcopy(self)
@ -161,7 +161,7 @@ class Grid:
@staticmethod
def load(fname: Union[str, Path]) -> "Grid":
def load(fname: Union[str, Path]) -> 'Grid':
"""
Load from VTK image data file.
@ -190,7 +190,7 @@ class Grid:
@typing. no_type_check
@staticmethod
def load_ASCII(fname)-> "Grid":
def load_ASCII(fname)-> 'Grid':
"""
Load from geom file.
@ -264,7 +264,7 @@ class Grid:
@staticmethod
def load_Neper(fname: Union[str, Path]) -> "Grid":
def load_Neper(fname: Union[str, Path]) -> 'Grid':
"""
Load from Neper VTK file.
@ -279,7 +279,7 @@ class Grid:
Grid-based geometry from file.
"""
v = VTK.load(fname,'vtkImageData')
v = VTK.load(fname,'ImageData')
cells = np.array(v.vtk_data.GetDimensions())-1
bbox = np.array(v.vtk_data.GetBounds()).reshape(3,2).T
@ -292,7 +292,7 @@ class Grid:
def load_DREAM3D(fname: Union[str, Path],
feature_IDs: str = None, cell_data: str = None,
phases: str = 'Phases', Euler_angles: str = 'EulerAngles',
base_group: str = None) -> "Grid":
base_group: str = None) -> 'Grid':
"""
Load DREAM.3D (HDF5) file.
@ -354,7 +354,7 @@ class Grid:
@staticmethod
def from_table(table: Table,
coordinates: str,
labels: Union[str, Sequence[str]]) -> "Grid":
labels: Union[str, Sequence[str]]) -> 'Grid':
"""
Create grid from ASCII table.
@ -422,13 +422,14 @@ class Grid:
Grid-based geometry from tessellation.
"""
weights_p: FloatSequence
if periodic:
weights_p = np.tile(weights,27) # Laguerre weights (1,2,3,1,2,3,...,1,2,3)
seeds_p = np.vstack((seeds -np.array([size[0],0.,0.]),seeds, seeds +np.array([size[0],0.,0.])))
seeds_p = np.vstack((seeds_p-np.array([0.,size[1],0.]),seeds_p,seeds_p+np.array([0.,size[1],0.])))
seeds_p = np.vstack((seeds_p-np.array([0.,0.,size[2]]),seeds_p,seeds_p+np.array([0.,0.,size[2]])))
else:
weights_p = weights
weights_p = np.array(weights,float)
seeds_p = seeds
coords = grid_filters.coordinates0_point(cells,size).reshape(-1,3)
@ -452,7 +453,7 @@ class Grid:
size: FloatSequence,
seeds: np.ndarray,
material: IntSequence = None,
periodic: bool = True) -> "Grid":
periodic: bool = True) -> 'Grid':
"""
Create grid from Voronoi tessellation.
@ -538,7 +539,7 @@ class Grid:
surface: str,
threshold: float = 0.0,
periods: int = 1,
materials: IntSequence = (0,1)) -> "Grid":
materials: IntSequence = (0,1)) -> 'Grid':
"""
Create grid from definition of triply periodic minimal surface.
@ -674,7 +675,7 @@ class Grid:
def show(self) -> None:
"""Show on screen."""
VTK.from_rectilinear_grid(self.cells,self.size,self.origin).show()
VTK.from_image_data(self.cells,self.size,self.origin).show()
def add_primitive(self,
@ -684,7 +685,7 @@ class Grid:
fill: int = None,
R: Rotation = Rotation(),
inverse: bool = False,
periodic: bool = True) -> "Grid":
periodic: bool = True) -> 'Grid':
"""
Insert a primitive geometric object at a given position.
@ -769,7 +770,7 @@ class Grid:
)
def mirror(self, directions: Sequence[str], reflect: bool = False) -> "Grid":
def mirror(self, directions: Sequence[str], reflect: bool = False) -> 'Grid':
"""
Mirror grid along given directions.
@ -821,7 +822,7 @@ class Grid:
)
def flip(self, directions: Sequence[str]) -> "Grid":
def flip(self, directions: Union[Literal['x', 'y', 'z'], Sequence[Literal['x', 'y', 'z']]]) -> 'Grid':
"""
Flip grid along given directions.
@ -841,7 +842,8 @@ class Grid:
if not set(directions).issubset(valid):
raise ValueError(f'invalid direction {set(directions).difference(valid)} specified')
mat = np.flip(self.material, (valid.index(d) for d in directions if d in valid))
mat = np.flip(self.material, [valid.index(d) for d in directions if d in valid])
return Grid(material = mat,
size = self.size,
@ -850,7 +852,7 @@ class Grid:
)
def scale(self, cells: IntSequence, periodic: bool = True) -> "Grid":
def scale(self, cells: IntSequence, periodic: bool = True) -> 'Grid':
"""
Scale grid to new cells.
@ -897,7 +899,7 @@ class Grid:
def clean(self,
stencil: int = 3,
selection: IntSequence = None,
periodic: bool = True) -> "Grid":
periodic: bool = True) -> 'Grid':
"""
Smooth grid by selecting most frequent material index within given stencil at each location.
@ -937,7 +939,7 @@ class Grid:
)
def renumber(self) -> "Grid":
def renumber(self) -> 'Grid':
"""
Renumber sorted material indices as 0,...,N-1.
@ -956,7 +958,7 @@ class Grid:
)
def rotate(self, R: Rotation, fill: int = None) -> "Grid":
def rotate(self, R: Rotation, fill: int = None) -> 'Grid':
"""
Rotate grid (pad if required).
@ -973,14 +975,13 @@ class Grid:
Updated grid-based geometry.
"""
if fill is None: fill = np.nanmax(self.material) + 1
dtype = float if isinstance(fill,float) or self.material.dtype in np.sctypes['float'] else int
material = self.material
# These rotations are always applied in the reference coordinate system, i.e. (z,x,z) not (z,x',z'')
# see https://www.cs.utexas.edu/~theshark/courses/cs354/lectures/cs354-14.pdf
for angle,axes in zip(R.as_Euler_angles(degrees=True)[::-1], [(0,1),(1,2),(0,1)]):
material_temp = ndimage.rotate(material,angle,axes,order=0,prefilter=False,output=dtype,cval=fill)
material_temp = ndimage.rotate(material,angle,axes,order=0,prefilter=False,
output=self.material.dtype,
cval=np.nanmax(self.material) + 1 if fill is None else fill)
# avoid scipy interpolation errors for rotations close to multiples of 90°
material = material_temp if np.prod(material_temp.shape) != np.prod(material.shape) else \
np.rot90(material,k=np.rint(angle/90.).astype(int),axes=axes)
@ -997,7 +998,7 @@ class Grid:
def canvas(self,
cells: IntSequence = None,
offset: IntSequence = None,
fill: int = None) -> "Grid":
fill: int = None) -> 'Grid':
"""
Crop or enlarge/pad grid.
@ -1031,10 +1032,8 @@ class Grid:
"""
offset_ = np.array(offset,int) if offset is not None else np.zeros(3,int)
cells_ = np.array(cells,int) if cells is not None else self.cells
if fill is None: fill = np.nanmax(self.material) + 1
dtype = float if int(fill) != fill or self.material.dtype in np.sctypes['float'] else int
canvas = np.full(cells_,fill,dtype)
canvas = np.full(cells_,np.nanmax(self.material) + 1 if fill is None else fill,self.material.dtype)
LL = np.clip( offset_, 0,np.minimum(self.cells, cells_+offset_))
UR = np.clip( offset_+cells_, 0,np.minimum(self.cells, cells_+offset_))
@ -1050,7 +1049,7 @@ class Grid:
)
def substitute(self, from_material: IntSequence, to_material: IntSequence) -> "Grid":
def substitute(self, from_material: IntSequence, to_material: IntSequence) -> 'Grid':
"""
Substitute material indices.
@ -1067,20 +1066,18 @@ class Grid:
Updated grid-based geometry.
"""
def mp(entry, mapper):
return mapper[entry] if entry in mapper else entry
material = self.material.copy()
for f,t in zip(from_material,to_material): # ToDo Python 3.10 has strict mode for zip
material[self.material==f] = t
mp = np.vectorize(mp)
mapper = dict(zip(from_material,to_material))
return Grid(material = mp(self.material,mapper).reshape(self.cells),
return Grid(material = material,
size = self.size,
origin = self.origin,
comments = self.comments+[util.execution_stamp('Grid','substitute')],
)
def sort(self) -> "Grid":
def sort(self) -> 'Grid':
"""
Sort material indices such that min(material) is located at (0,0,0).
@ -1106,7 +1103,7 @@ class Grid:
vicinity: int = 1,
offset: int = None,
trigger: IntSequence = [],
periodic: bool = True) -> "Grid":
periodic: bool = True) -> 'Grid':
"""
Offset material index of points in the vicinity of xxx.

2
python/damask/_result.py
View File

@ -981,7 +981,7 @@ class Result:
t = 'tensor'
if o is None: o = 'fro'
else:
raise ValueError(f'invalid norm order {ord}')
raise ValueError(f'invalid shape of {x["label"]}')
return {
'data': np.linalg.norm(x['data'],ord=o,axis=axis,keepdims=True),

6
python/damask/_rotation.py
View File

@ -812,8 +812,7 @@ class Rotation:
return Rotation.from_basis(R)
@staticmethod
def from_parallel(a,b,
**kwargs):
def from_parallel(a,b):
"""
Initialize from pairs of two orthogonal lattice basis vectors.
@ -963,8 +962,7 @@ class Rotation:
N = 500,
degrees = True,
fractions = True,
rng_seed = None,
**kwargs):
rng_seed = None):
"""
Sample discrete values from a binned orientation distribution function (ODF).

56
python/damask/_vtk.py
View File

@ -2,6 +2,7 @@ import os
import warnings
import multiprocessing as mp
from pathlib import Path
from typing import Union, Literal, List
import numpy as np
import vtk
@ -9,6 +10,7 @@ from vtk.util.numpy_support import numpy_to_vtk as np_to_vtk
from vtk.util.numpy_support import numpy_to_vtkIdTypeArray as np_to_vtkIdTypeArray
from vtk.util.numpy_support import vtk_to_numpy as vtk_to_np
from ._typehints import FloatSequence, IntSequence
from . import util
from . import Table
@ -20,7 +22,7 @@ class VTK:
High-level interface to VTK.
"""
def __init__(self,vtk_data):
def __init__(self, vtk_data: vtk.vtkDataSet):
"""
New spatial visualization.
@ -36,7 +38,7 @@ class VTK:
@staticmethod
def from_image_data(cells,size,origin=np.zeros(3)):
def from_image_data(cells: IntSequence, size: FloatSequence, origin: FloatSequence = np.zeros(3)) -> 'VTK':
"""
Create VTK of type vtk.vtkImageData.
@ -60,13 +62,13 @@ class VTK:
vtk_data = vtk.vtkImageData()
vtk_data.SetDimensions(*(np.array(cells)+1))
vtk_data.SetOrigin(*(np.array(origin)))
vtk_data.SetSpacing(*(size/cells))
vtk_data.SetSpacing(*(np.array(size)/np.array(cells)))
return VTK(vtk_data)
@staticmethod
def from_rectilinear_grid(grid,size,origin=np.zeros(3)):
def from_rectilinear_grid(grid: np.ndarray, size: FloatSequence, origin: FloatSequence = np.zeros(3)) -> 'VTK':
"""
Create VTK of type vtk.vtkRectilinearGrid.
@ -98,7 +100,7 @@ class VTK:
@staticmethod
def from_unstructured_grid(nodes,connectivity,cell_type):
def from_unstructured_grid(nodes: np.ndarray, connectivity: np.ndarray, cell_type: str) -> 'VTK':
"""
Create VTK of type vtk.vtkUnstructuredGrid.
@ -138,7 +140,7 @@ class VTK:
@staticmethod
def from_poly_data(points):
def from_poly_data(points: np.ndarray) -> 'VTK':
"""
Create VTK of type vtk.polyData.
@ -172,15 +174,17 @@ class VTK:
@staticmethod
def load(fname,dataset_type=None):
def load(fname: Union[str, Path],
dataset_type: Literal['ImageData', 'UnstructuredGrid', 'PolyData'] = None) -> 'VTK':
"""
Load from VTK file.
Parameters
----------
fname : str or pathlib.Path
Filename for reading. Valid extensions are .vti, .vtr, .vtu, .vtp, and .vtk.
dataset_type : {'vtkImageData', ''vtkRectilinearGrid', 'vtkUnstructuredGrid', 'vtkPolyData'}, optional
Filename for reading.
Valid extensions are .vti, .vtr, .vtu, .vtp, and .vtk.
dataset_type : {'ImageData', 'UnstructuredGrid', 'PolyData'}, optional
Name of the vtk.vtkDataSet subclass when opening a .vtk file.
Returns
@ -234,7 +238,7 @@ class VTK:
def _write(writer):
"""Wrapper for parallel writing."""
writer.Write()
def save(self,fname,parallel=True,compress=True):
def save(self, fname: Union[str, Path], parallel: bool = True, compress: bool = True):
"""
Save as VTK file.
@ -280,7 +284,7 @@ class VTK:
# Check https://blog.kitware.com/ghost-and-blanking-visibility-changes/ for missing data
# Needs support for damask.Table
def add(self,data,label=None):
def add(self, data: Union[np.ndarray, np.ma.MaskedArray], label: str = None):
"""
Add data to either cells or points.
@ -327,7 +331,7 @@ class VTK:
raise TypeError
def get(self,label):
def get(self, label: str) -> np.ndarray:
"""
Get either cell or point data.
@ -369,7 +373,7 @@ class VTK:
raise ValueError(f'Array "{label}" not found.')
def get_comments(self):
def get_comments(self) -> List[str]:
"""Return the comments."""
fielddata = self.vtk_data.GetFieldData()
for a in range(fielddata.GetNumberOfArrays()):
@ -379,7 +383,7 @@ class VTK:
return []
def set_comments(self,comments):
def set_comments(self, comments: Union[str, List[str]]):
"""
Set comments.
@ -396,7 +400,7 @@ class VTK:
self.vtk_data.GetFieldData().AddArray(s)
def add_comments(self,comments):
def add_comments(self, comments: Union[str, List[str]]):
"""
Add comments.
@ -409,7 +413,7 @@ class VTK:
self.set_comments(self.get_comments() + ([comments] if isinstance(comments,str) else comments))
def __repr__(self):
def __repr__(self) -> str:
"""ASCII representation of the VTK data."""
writer = vtk.vtkDataSetWriter()
writer.SetHeader(f'# {util.execution_stamp("VTK")}')
@ -419,7 +423,7 @@ class VTK:
return writer.GetOutputString()
def show(self) -> None:
def show(self):
"""
Render.
@ -442,22 +446,24 @@ class VTK:
mapper = vtk.vtkDataSetMapper()
mapper.SetInputData(self.vtk_data)
actor = vtk.vtkActor()
actor.SetMapper(mapper)
actor.GetProperty().SetColor(230/255,150/255,68/255)
ren = vtk.vtkRenderer()
ren.AddActor(actor)
ren.SetBackground(67/255,128/255,208/255)
window = vtk.vtkRenderWindow()
window.AddRenderer(ren)
ren.AddActor(actor)
ren.SetBackground(0.2,0.2,0.2)
window.SetSize(width,height)
window.SetWindowName(util.execution_stamp('VTK','show'))
iren = vtk.vtkRenderWindowInteractor()
iren.SetRenderWindow(window)
iren.Initialize()
window.Render()
iren.Start()
if os.name == 'posix' and 'DISPLAY' not in os.environ:
print('Found no rendering device')
else:
window.Render()
iren.Start()

2
python/damask/seeds.py
View File

@ -79,7 +79,7 @@ def from_Poisson_disc(size: _FloatSequence, N_seeds: int, N_candidates: int, dis
s = 1
i = 0
progress = _util._ProgressBar(N_seeds+1,'',50)
progress = _util.ProgressBar(N_seeds+1,'',50)
while s < N_seeds:
i += 1
candidates = rng.random((N_candidates,3))*_np.broadcast_to(size,(N_candidates,3))

167
python/damask/util.py
View File

@ -7,12 +7,16 @@ import subprocess
import shlex
import re
import fractions
import collections.abc as abc
from functools import reduce
from typing import Union, Tuple, Iterable, Callable, Dict, List, Any, Literal, Optional
from pathlib import Path
import numpy as np
import h5py
from . import version
from ._typehints import FloatSequence
# limit visibility
__all__=[
@ -50,16 +54,16 @@ _colors = {
####################################################################################################
# Functions
####################################################################################################
def srepr(arg,glue = '\n'):
def srepr(msg, glue: str = '\n') -> str:
r"""
Join items with glue string.
Parameters
----------
arg : iterable
msg : object with __repr__ or sequence of objects with __repr__
Items to join.
glue : str, optional
Glue used for joining operation. Defaults to \n.
Glue used for joining operation. Defaults to '\n'.
Returns
-------
@ -67,21 +71,21 @@ def srepr(arg,glue = '\n'):
String representation of the joined items.
"""
if (not hasattr(arg, 'strip') and
(hasattr(arg, '__getitem__') or
hasattr(arg, '__iter__'))):
return glue.join(str(x) for x in arg)
if (not hasattr(msg, 'strip') and
(hasattr(msg, '__getitem__') or
hasattr(msg, '__iter__'))):
return glue.join(str(x) for x in msg)
else:
return arg if isinstance(arg,str) else repr(arg)
return msg if isinstance(msg,str) else repr(msg)
def emph(what):
def emph(msg) -> str:
"""
Format with emphasis.
Parameters
----------
what : object with __repr__ or iterable of objects with __repr__.
msg : object with __repr__ or sequence of objects with __repr__
Message to format.
Returns
@ -90,15 +94,15 @@ def emph(what):
Formatted string representation of the joined items.
"""
return _colors['bold']+srepr(what)+_colors['end_color']
return _colors['bold']+srepr(msg)+_colors['end_color']
def deemph(what):
def deemph(msg) -> str:
"""
Format with deemphasis.
Parameters
----------
what : object with __repr__ or iterable of objects with __repr__.
msg : object with __repr__ or sequence of objects with __repr__
Message to format.
Returns
@ -107,15 +111,15 @@ def deemph(what):
Formatted string representation of the joined items.
"""
return _colors['dim']+srepr(what)+_colors['end_color']
return _colors['dim']+srepr(msg)+_colors['end_color']
def warn(what):
def warn(msg) -> str:
"""
Format for warning.
Parameters
----------
what : object with __repr__ or iterable of objects with __repr__.
msg : object with __repr__ or sequence of objects with __repr__
Message to format.
Returns
@ -124,15 +128,15 @@ def warn(what):
Formatted string representation of the joined items.
"""
return _colors['warning']+emph(what)+_colors['end_color']
return _colors['warning']+emph(msg)+_colors['end_color']
def strikeout(what):
def strikeout(msg) -> str:
"""
Format as strikeout.
Parameters
----------
what : object with __repr__ or iterable of objects with __repr__.
msg : object with __repr__ or iterable of objects with __repr__
Message to format.
Returns
@ -141,10 +145,10 @@ def strikeout(what):
Formatted string representation of the joined items.
"""
return _colors['crossout']+srepr(what)+_colors['end_color']
return _colors['crossout']+srepr(msg)+_colors['end_color']
def run(cmd,wd='./',env=None,timeout=None):
def run(cmd: str, wd: str = './', env: Dict[str, str] = None, timeout: int = None) -> Tuple[str, str]:
"""
Run a command.
@ -153,7 +157,7 @@ def run(cmd,wd='./',env=None,timeout=None):
cmd : str
Command to be executed.
wd : str, optional
Working directory of process. Defaults to ./ .
Working directory of process. Defaults to './'.
env : dict, optional
Environment for execution.
timeout : integer, optional
@ -185,7 +189,7 @@ def run(cmd,wd='./',env=None,timeout=None):
execute = run
def natural_sort(key):
def natural_sort(key: str) -> List[Union[int, str]]:
"""
Natural sort.
@ -200,7 +204,10 @@ def natural_sort(key):
return [ convert(c) for c in re.split('([0-9]+)', key) ]
def show_progress(iterable,N_iter=None,prefix='',bar_length=50):
def show_progress(iterable: Iterable,
N_iter: int = None,
prefix: str = '',
bar_length: int = 50) -> Any:
"""
Decorate a loop with a progress bar.
@ -208,39 +215,49 @@ def show_progress(iterable,N_iter=None,prefix='',bar_length=50):
Parameters
----------
iterable : iterable or function with yield statement
Iterable (or function with yield statement) to be decorated.
iterable : iterable
Iterable to be decorated.
N_iter : int, optional
Total number of iterations. Required unless obtainable as len(iterable).
Total number of iterations. Required if iterable is not a sequence.
prefix : str, optional
Prefix string.
bar_length : int, optional
Length of progress bar in characters. Defaults to 50.
"""
if N_iter in [0,1] or (hasattr(iterable,'__len__') and len(iterable) <= 1):
if isinstance(iterable,abc.Sequence):
if N_iter is None:
N = len(iterable)
else:
raise ValueError('N_iter given for sequence')
else:
if N_iter is None:
raise ValueError('N_iter not given')
else:
N = N_iter
if N <= 1:
for item in iterable:
yield item
else:
status = _ProgressBar(N_iter if N_iter is not None else len(iterable),prefix,bar_length)
status = ProgressBar(N,prefix,bar_length)
for i,item in enumerate(iterable):
yield item
status.update(i)
def scale_to_coprime(v):
def scale_to_coprime(v: FloatSequence) -> np.ndarray:
"""
Scale vector to co-prime (relatively prime) integers.
Parameters
----------
v : numpy.ndarray of shape (:)
v : sequence of float, len (:)
Vector to scale.
Returns
-------
m : numpy.ndarray of shape (:)
m : numpy.ndarray, shape (:)
Vector scaled to co-prime numbers.
"""
@ -257,17 +274,21 @@ def scale_to_coprime(v):
except AttributeError:
return a * b // np.gcd(a, b)
m = (np.array(v) * reduce(lcm, map(lambda x: int(get_square_denominator(x)),v)) ** 0.5).astype(int)
v_ = np.array(v)
m = (v_ * reduce(lcm, map(lambda x: int(get_square_denominator(x)),v_))**0.5).astype(int)
m = m//reduce(np.gcd,m)
with np.errstate(invalid='ignore'):
if not np.allclose(np.ma.masked_invalid(v/m),v[np.argmax(abs(v))]/m[np.argmax(abs(v))]):
raise ValueError(f'Invalid result {m} for input {v}. Insufficient precision?')
if not np.allclose(np.ma.masked_invalid(v_/m),v_[np.argmax(abs(v_))]/m[np.argmax(abs(v_))]):
raise ValueError(f'Invalid result {m} for input {v_}. Insufficient precision?')
return m
def project_equal_angle(vector,direction='z',normalize=True,keepdims=False):
def project_equal_angle(vector: np.ndarray,
direction: Literal['x', 'y', 'z'] = 'z',
normalize: bool = True,
keepdims: bool = False) -> np.ndarray:
"""
Apply equal-angle projection to vector.
@ -275,9 +296,8 @@ def project_equal_angle(vector,direction='z',normalize=True,keepdims=False):
----------
vector : numpy.ndarray, shape (...,3)
Vector coordinates to be projected.
direction : str
Projection direction 'x', 'y', or 'z'.
Defaults to 'z'.
direction : {'x', 'y', 'z'}
Projection direction. Defaults to 'z'.
normalize : bool
Ensure unit length of input vector. Defaults to True.
keepdims : bool
@ -309,7 +329,10 @@ def project_equal_angle(vector,direction='z',normalize=True,keepdims=False):
return np.roll(np.block([v[...,:2]/(1.0+np.abs(v[...,2:3])),np.zeros_like(v[...,2:3])]),
-shift if keepdims else 0,axis=-1)[...,:3 if keepdims else 2]
def project_equal_area(vector,direction='z',normalize=True,keepdims=False):
def project_equal_area(vector: np.ndarray,
direction: Literal['x', 'y', 'z'] = 'z',
normalize: bool = True,
keepdims: bool = False) -> np.ndarray:
"""
Apply equal-area projection to vector.
@ -317,9 +340,8 @@ def project_equal_area(vector,direction='z',normalize=True,keepdims=False):
----------
vector : numpy.ndarray, shape (...,3)
Vector coordinates to be projected.
direction : str
Projection direction 'x', 'y', or 'z'.
Defaults to 'z'.
direction : {'x', 'y', 'z'}
Projection direction. Defaults to 'z'.
normalize : bool
Ensure unit length of input vector. Defaults to True.
keepdims : bool
@ -351,15 +373,14 @@ def project_equal_area(vector,direction='z',normalize=True,keepdims=False):
return np.roll(np.block([v[...,:2]/np.sqrt(1.0+np.abs(v[...,2:3])),np.zeros_like(v[...,2:3])]),
-shift if keepdims else 0,axis=-1)[...,:3 if keepdims else 2]
def execution_stamp(class_name,function_name=None):
def execution_stamp(class_name: str, function_name: str = None) -> str:
"""Timestamp the execution of a (function within a) class."""
now = datetime.datetime.now().astimezone().strftime('%Y-%m-%d %H:%M:%S%z')
_function_name = '' if function_name is None else f'.{function_name}'
return f'damask.{class_name}{_function_name} v{version} ({now})'
def hybrid_IA(dist,N,rng_seed=None):
def hybrid_IA(dist: np.ndarray, N: int, rng_seed = None) -> np.ndarray:
"""
Hybrid integer approximation.
@ -387,7 +408,10 @@ def hybrid_IA(dist,N,rng_seed=None):
return np.repeat(np.arange(len(dist)),repeats)[np.random.default_rng(rng_seed).permutation(N_inv_samples)[:N]]
def shapeshifter(fro,to,mode='left',keep_ones=False):
def shapeshifter(fro: Tuple[int, ...],
to: Tuple[int, ...],
mode: Literal['left','right'] = 'left',
keep_ones: bool = False) -> Tuple[Optional[int], ...]:
"""
Return dimensions that reshape 'fro' to become broadcastable to 'to'.
@ -398,9 +422,9 @@ def shapeshifter(fro,to,mode='left',keep_ones=False):
to : tuple
Target shape of array after broadcasting.
len(to) cannot be less than len(fro).
mode : str, optional
mode : {'left', 'right'}, optional
Indicates whether new axes are preferably added to
either 'left' or 'right' of the original shape.
either left or right of the original shape.
Defaults to 'left'.
keep_ones : bool, optional
Treat '1' in fro as literal value instead of dimensional placeholder.
@ -434,21 +458,22 @@ def shapeshifter(fro,to,mode='left',keep_ones=False):
fro = (1,) if not len(fro) else fro
to = (1,) if not len(to) else to
try:
grp = re.match(beg[mode]
match = re.match(beg[mode]
+f',{sep[mode]}'.join(map(lambda x: f'{x}'
if x>1 or (keep_ones and len(fro)>1) else
'\\d+',fro))
+f',{end[mode]}',
','.join(map(str,to))+',').groups()
except AttributeError:
+f',{end[mode]}',','.join(map(str,to))+',')
assert match
grp = match.groups()
except AssertionError:
raise ValueError(f'Shapes can not be shifted {fro} --> {to}')
fill = ()
fill: Tuple[Optional[int], ...] = ()
for g,d in zip(grp,fro+(None,)):
fill += (1,)*g.count(',')+(d,)
return fill[:-1]
def shapeblender(a,b):
def shapeblender(a: Tuple[int, ...], b: Tuple[int, ...]) -> Tuple[int, ...]:
"""
Return a shape that overlaps the rightmost entries of 'a' with the leftmost of 'b'.
@ -476,7 +501,7 @@ def shapeblender(a,b):
return a + b[i:]
def extend_docstring(extra_docstring):
def extend_docstring(extra_docstring: str) -> Callable:
"""
Decorator: Append to function's docstring.
@ -492,7 +517,7 @@ def extend_docstring(extra_docstring):
return _decorator
def extended_docstring(f,extra_docstring):
def extended_docstring(f: Callable, extra_docstring: str) -> Callable:
"""
Decorator: Combine another function's docstring with a given docstring.
@ -510,7 +535,7 @@ def extended_docstring(f,extra_docstring):
return _decorator
def DREAM3D_base_group(fname):
def DREAM3D_base_group(fname: Union[str, Path]) -> str:
"""
Determine the base group of a DREAM.3D file.
@ -536,7 +561,7 @@ def DREAM3D_base_group(fname):
return base_group
def DREAM3D_cell_data_group(fname):
def DREAM3D_cell_data_group(fname: Union[str, Path]) -> str:
"""
Determine the cell data group of a DREAM.3D file.
@ -568,18 +593,18 @@ def DREAM3D_cell_data_group(fname):
return cell_data_group
def Bravais_to_Miller(*,uvtw=None,hkil=None):
def Bravais_to_Miller(*, uvtw: np.ndarray = None, hkil: np.ndarray = None) -> np.ndarray:
"""
Transform 4 Miller–Bravais indices to 3 Miller indices of crystal direction [uvw] or plane normal (hkl).
Parameters
----------
uvtw|hkil : numpy.ndarray of shape (...,4)
uvtw|hkil : numpy.ndarray, shape (...,4)
Miller–Bravais indices of crystallographic direction [uvtw] or plane normal (hkil).
Returns
-------
uvw|hkl : numpy.ndarray of shape (...,3)
uvw|hkl : numpy.ndarray, shape (...,3)
Miller indices of [uvw] direction or (hkl) plane normal.
"""
@ -595,18 +620,18 @@ def Bravais_to_Miller(*,uvtw=None,hkil=None):
return np.einsum('il,...l',basis,axis)
def Miller_to_Bravais(*,uvw=None,hkl=None):
def Miller_to_Bravais(*, uvw: np.ndarray = None, hkl: np.ndarray = None) -> np.ndarray:
"""
Transform 3 Miller indices to 4 Miller–Bravais indices of crystal direction [uvtw] or plane normal (hkil).
Parameters
----------
uvw|hkl : numpy.ndarray of shape (...,3)
uvw|hkl : numpy.ndarray, shape (...,3)
Miller indices of crystallographic direction [uvw] or plane normal (hkl).
Returns
-------
uvtw|hkil : numpy.ndarray of shape (...,4)
uvtw|hkil : numpy.ndarray, shape (...,4)
Miller–Bravais indices of [uvtw] direction or (hkil) plane normal.
"""
@ -624,7 +649,7 @@ def Miller_to_Bravais(*,uvw=None,hkl=None):
return np.einsum('il,...l',basis,axis)
def dict_prune(d):
def dict_prune(d: Dict) -> Dict:
"""
Recursively remove empty dictionaries.
@ -650,7 +675,7 @@ def dict_prune(d):
return new
def dict_flatten(d):
def dict_flatten(d: Dict) -> Dict:
"""
Recursively remove keys of single-entry dictionaries.
@ -678,14 +703,14 @@ def dict_flatten(d):
####################################################################################################
# Classes
####################################################################################################
class _ProgressBar:
class ProgressBar:
"""
Report progress of an interation as a status bar.
Works for 0-based loops, ETA is estimated by linear extrapolation.
"""
def __init__(self,total,prefix,bar_length):
def __init__(self, total: int, prefix: str, bar_length: int):
"""
Set current time as basis for ETA estimation.
@ -708,7 +733,7 @@ class _ProgressBar:
sys.stderr.write(f"{self.prefix} {'░'*self.bar_length} 0% ETA n/a")
sys.stderr.flush()
def update(self,iteration):
def update(self, iteration: int) -> None:
fraction = (iteration+1) / self.total
filled_length = int(self.bar_length * fraction)

1
python/tests/test_Colormap.py
View File

@ -81,6 +81,7 @@ class TestColormap:
assert Colormap.from_predefined('strain') == Colormap.from_predefined('strain')
assert Colormap.from_predefined('strain') != Colormap.from_predefined('stress')
assert Colormap.from_predefined('strain',N=128) != Colormap.from_predefined('strain',N=64)
assert not Colormap.from_predefined('strain',N=128) == 1
@pytest.mark.parametrize('low,high',[((0,0,0),(1,1,1)),
([0,0,0],[1,1,1])])

3
python/tests/test_Crystal.py
View File

@ -40,6 +40,9 @@ class TestCrystal:
alpha=alpha,beta=beta,gamma=gamma)
assert np.allclose(np.eye(3),np.einsum('ik,jk',c.basis_real,c.basis_reciprocal))
def test_basis_invalid(self):
with pytest.raises(KeyError):
Crystal(family='cubic').basis_real
@pytest.mark.parametrize('keyFrame,keyLattice',[('uvw','direction'),('hkl','plane'),])
@pytest.mark.parametrize('vector',np.array([

1
python/tests/test_Grid.py
View File

@ -44,6 +44,7 @@ class TestGrid:
def test_equal(self,default):
assert default == default
assert not default == 42
def test_repr(self,default):
print(default)

9
python/tests/test_Orientation.py
View File

@ -364,6 +364,11 @@ class TestOrientation:
table.save(reference)
assert np.allclose(P,Table.load(reference).get('Schmid'))
def test_Schmid_invalid(self):
with pytest.raises(KeyError):
Orientation(lattice='fcc').Schmid()
### vectorization tests ###
@pytest.mark.parametrize('lattice',['hP','cI','cF']) # tI not included yet
@ -505,3 +510,7 @@ class TestOrientation:
for loc in np.random.randint(0,blend,(10,len(blend))):
assert np.allclose(o[tuple(loc[:len(o.shape)])].to_pole(uvw=v[tuple(loc[-len(v.shape[:-1]):])]),
o.to_pole(uvw=v)[tuple(loc)])
def test_mul_invalid(self):
with pytest.raises(TypeError):
Orientation.from_random(lattice='cF')*np.ones(3)

3
python/tests/test_Result.py
View File

@ -102,6 +102,9 @@ class TestResult:
with pytest.raises(AttributeError):
default.view('invalid',True)
def test_add_invalid(self,default):
default.add_absolute('xxxx')
def test_add_absolute(self,default):
default.add_absolute('F_e')
in_memory = np.abs(default.place('F_e'))

9
python/tests/test_Rotation.py
View File

@ -792,6 +792,11 @@ class TestRotation:
R = Rotation.from_random(shape,rng_seed=1)
assert R == R if shape is None else (R == R).all()
@pytest.mark.parametrize('shape',[None,5,(4,6)])
def test_allclose(self,shape):
R = Rotation.from_random(shape,rng_seed=1)
assert R.allclose(R)
@pytest.mark.parametrize('shape',[None,5,(4,6)])
def test_unequal(self,shape):
R = Rotation.from_random(shape,rng_seed=1)
@ -1124,3 +1129,7 @@ class TestRotation:
weights_r = np.histogramdd(Eulers_r,steps,rng)[0].flatten(order='F')/N * np.sum(weights)
assert np.sqrt(((weights_r - weights) ** 2).mean()) < 5
def test_mul_invalid(self):
with pytest.raises(TypeError):
Rotation.from_random()*np.ones(3)

4
python/tests/test_VTK.py
View File

@ -28,6 +28,10 @@ class TestVTK:
def _patch_execution_stamp(self, patch_execution_stamp):
print('patched damask.util.execution_stamp')
def test_show(sef,default,monkeypatch):
monkeypatch.delenv('DISPLAY',raising=False)
default.show()
def test_rectilinearGrid(self,tmp_path):
cells = np.random.randint(5,10,3)*2
size = np.random.random(3) + 1.0

286
src/YAML_types.f90
View File

@ -62,6 +62,8 @@ module YAML_types
tNode_get_byKey_as1dString => tNode_get_byKey_as1dString
procedure :: &
getKey => tNode_get_byIndex_asKey
procedure :: &
Keys => tNode_getKeys
procedure :: &
getIndex => tNode_get_byKey_asIndex
procedure :: &
@ -117,7 +119,7 @@ module YAML_types
type, extends(tNode), public :: tList
class(tItem), pointer :: first => null()
class(tItem), pointer :: first => NULL()
contains
procedure :: asFormattedString => tList_asFormattedString
@ -144,8 +146,8 @@ module YAML_types
type :: tItem
character(len=:), allocatable :: key
class(tNode), pointer :: node => null()
class(tItem), pointer :: next => null()
class(tNode), pointer :: node => NULL()
class(tItem), pointer :: next => NULL()
contains
final :: tItem_finalize
@ -221,22 +223,22 @@ subroutine selfTest
select type(s1)
class is(tScalar)
s1 = '2'
endselect
end select
select type(s2)
class is(tScalar)
s2 = '3'
endselect
end select
select type(s3)
class is(tScalar)
s3 = '4'
endselect
end select
select type(s4)
class is(tScalar)
s4 = '5'
endselect
end select
allocate(tList::l1)
@ -249,14 +251,14 @@ subroutine selfTest
if (any(dNeq(l1%as1dFloat(),[2.0_pReal,3.0_pReal]))) error stop 'tList_as1dFloat'
if (n%get_asInt(1) /= 2) error stop 'byIndex_asInt'
if (dNeq(n%get_asFloat(2),3.0_pReal)) error stop 'byIndex_asFloat'
endselect
end select
allocate(tList::l3)
select type(l3)
class is(tList)
call l3%append(s3)
call l3%append(s4)
endselect
end select
allocate(tList::l2)
select type(l2)
@ -332,9 +334,12 @@ function tNode_asScalar(self) result(scalar)
class(tNode), intent(in), target :: self
class(tScalar), pointer :: scalar
select type(self)
class is(tScalar)
scalar => self
class default
nullify(scalar)
end select
end function tNode_asScalar
@ -348,9 +353,12 @@ function tNode_asList(self) result(list)
class(tNode), intent(in), target :: self
class(tList), pointer :: list
select type(self)
class is(tList)
list => self
class default
nullify(list)
end select
end function tNode_asList
@ -364,9 +372,12 @@ function tNode_asDict(self) result(dict)
class(tNode), intent(in), target :: self
class(tDict), pointer :: dict
select type(self)
class is(tDict)
dict => self
class default
nullify(dict)
end select
end function tNode_asDict
@ -385,12 +396,13 @@ function tNode_get_byIndex(self,i) result(node)
class(tItem), pointer :: item
integer :: j
select type(self)
class is(tList)
self_ => self%asList()
class default
call IO_error(706,ext_msg='Expected list')
endselect
end select
item => self_%first
@ -409,15 +421,14 @@ end function tNode_get_byIndex
!--------------------------------------------------------------------------------------------------
function tNode_get_byIndex_asFloat(self,i) result(nodeAsFloat)
class(tNode), intent(in), target :: self
integer, intent(in) :: i
class(tNode), intent(in) :: self
integer, intent(in) :: i
real(pReal) :: nodeAsFloat
class(tNode), pointer :: node
type(tScalar), pointer :: scalar
node => self%get(i)
select type(node)
select type(node => self%get(i))
class is(tScalar)
scalar => node%asScalar()
nodeAsFloat = scalar%asFloat()
@ -433,15 +444,15 @@ end function tNode_get_byIndex_asFloat
!--------------------------------------------------------------------------------------------------
function tNode_get_byIndex_asInt(self,i) result(nodeAsInt)
class(tNode), intent(in), target :: self
integer, intent(in) :: i
class(tNode), intent(in) :: self
integer, intent(in) :: i
integer :: nodeAsInt
class(tNode), pointer :: node
type(tScalar), pointer :: scalar
node => self%get(i)
select type(node)
select type(node => self%get(i))
class is(tScalar)
scalar => node%asScalar()
nodeAsInt = scalar%asInt()
@ -457,21 +468,20 @@ end function tNode_get_byIndex_asInt
!--------------------------------------------------------------------------------------------------
function tNode_get_byIndex_asBool(self,i) result(nodeAsBool)
class(tNode), intent(in), target :: self
integer, intent(in) :: i
class(tNode), intent(in) :: self
integer, intent(in) :: i
logical :: nodeAsBool
class(tNode), pointer :: node
type(tScalar), pointer :: scalar
node => self%get(i)
select type(node)
select type(node => self%get(i))
class is(tScalar)
scalar => node%asScalar()
nodeAsBool = scalar%asBool()
class default
call IO_error(706,ext_msg='Expected scalar Boolean')
endselect
end select
end function tNode_get_byIndex_asBool
@ -481,21 +491,20 @@ end function tNode_get_byIndex_asBool
!--------------------------------------------------------------------------------------------------
function tNode_get_byIndex_asString(self,i) result(nodeAsString)
class(tNode), intent(in), target :: self
integer, intent(in) :: i
class(tNode), intent(in) :: self
integer, intent(in) :: i
character(len=:), allocatable :: nodeAsString
class(tNode), pointer :: node
type(tScalar), pointer :: scalar
node => self%get(i)
select type(node)
select type(node => self%get(i))
class is(tScalar)
scalar => node%asScalar()
nodeAsString = scalar%asString()
class default
call IO_error(706,ext_msg='Expected scalar string')
endselect
end select
end function tNode_get_byIndex_asString
@ -505,21 +514,20 @@ end function tNode_get_byIndex_asString
!--------------------------------------------------------------------------------------------------
function tNode_get_byIndex_as1dFloat(self,i) result(nodeAs1dFloat)
class(tNode), intent(in), target :: self
integer, intent(in) :: i
class(tNode), intent(in) :: self
integer, intent(in) :: i
real(pReal), dimension(:), allocatable :: nodeAs1dFloat
class(tNode), pointer :: node
class(tList), pointer :: list
node => self%get(i)
select type(node)
select type(node => self%get(i))
class is(tList)
list => node%asList()
nodeAs1dFloat = list%as1dFloat()
class default
call IO_error(706,ext_msg='Expected list of floats')
endselect
end select
end function tNode_get_byIndex_as1dFloat
@ -529,21 +537,20 @@ end function tNode_get_byIndex_as1dFloat
!--------------------------------------------------------------------------------------------------
function tNode_get_byIndex_as1dInt(self,i) result(nodeAs1dInt)
class(tNode), intent(in), target :: self
integer, intent(in) :: i
class(tNode), intent(in) :: self
integer, intent(in) :: i
integer, dimension(:), allocatable :: nodeAs1dInt
class(tNode), pointer :: node
class(tList), pointer :: list
node => self%get(i)
select type(node)
select type(node => self%get(i))
class is(tList)
list => node%asList()
nodeAs1dInt = list%as1dInt()
class default
call IO_error(706,ext_msg='Expected list of integers')
endselect
end select
end function tNode_get_byIndex_as1dInt
@ -553,21 +560,20 @@ end function tNode_get_byIndex_as1dInt
!--------------------------------------------------------------------------------------------------
function tNode_get_byIndex_as1dBool(self,i) result(nodeAs1dBool)
class(tNode), intent(in), target :: self
integer, intent(in) :: i
class(tNode), intent(in) :: self
integer, intent(in) :: i
logical, dimension(:), allocatable :: nodeAs1dBool
class(tNode), pointer :: node
class(tList), pointer :: list
node => self%get(i)
select type(node)
select type(node => self%get(i))
class is(tList)
list => node%asList()
nodeAs1dBool = list%as1dBool()
class default
call IO_error(706,ext_msg='Expected list of Booleans')
endselect
end select
end function tNode_get_byIndex_as1dBool
@ -577,21 +583,20 @@ end function tNode_get_byIndex_as1dBool
!--------------------------------------------------------------------------------------------------
function tNode_get_byIndex_as1dString(self,i) result(nodeAs1dString)
class(tNode), intent(in), target :: self
integer, intent(in) :: i
class(tNode), intent(in) :: self
integer, intent(in) :: i
character(len=:), allocatable, dimension(:) :: nodeAs1dString
class(tNode), pointer :: node
type(tList), pointer :: list
node => self%get(i)
select type(node)
select type(node => self%get(i))
class is(tList)
list => node%asList()
nodeAs1dString = list%as1dString()
class default
call IO_error(706,ext_msg='Expected list of strings')
endselect
end select
end function tNode_get_byIndex_as1dString
@ -609,22 +614,50 @@ function tNode_get_byIndex_asKey(self,i) result(key)
type(tDict), pointer :: dict
type(tItem), pointer :: item
select type(self)
class is(tDict)
dict => self%asDict()
item => dict%first
do j = 1, min(i,dict%length)-1
item => item%next
enddo
end do
class default
call IO_error(706,ext_msg='Expected dict')
endselect
end select
key = item%key
end function tNode_get_byIndex_asKey
!--------------------------------------------------------------------------------------------------
!> @brief Get all keys from a dictionary
!--------------------------------------------------------------------------------------------------
function tNode_getKeys(self) result(keys)
class(tNode), intent(in) :: self
character(len=:), dimension(:), allocatable :: keys
character(len=pStringLen), dimension(:), allocatable :: temp
integer :: j, l
allocate(temp(self%length))
l = 0
do j = 1, self%length
temp(j) = self%getKey(j)
l = max(len_trim(temp(j)),l)
end do
allocate(character(l)::keys(self%length))
do j = 1, self%length
keys(j) = trim(temp(j))
end do
end function tNode_getKeys
!-------------------------------------------------------------------------------------------------
!> @brief Checks if a given key/item is present in the dict/list
!-------------------------------------------------------------------------------------------------
@ -658,7 +691,7 @@ function tNode_contains(self,k) result(exists)
enddo
class default
call IO_error(706,ext_msg='Expected list or dict')
endselect
end select
end function tNode_contains
@ -686,7 +719,7 @@ function tNode_get_byKey(self,k,defaultVal) result(node)
self_ => self%asDict()
class default
call IO_error(706,ext_msg='Expected dict for key '//k)
endselect
end select
j = 1
item => self_%first
@ -713,23 +746,22 @@ end function tNode_get_byKey
!--------------------------------------------------------------------------------------------------
function tNode_get_byKey_asFloat(self,k,defaultVal) result(nodeAsFloat)
class(tNode), intent(in), target :: self
character(len=*), intent(in) :: k
real(pReal), intent(in),optional :: defaultVal
class(tNode), intent(in) :: self
character(len=*), intent(in) :: k
real(pReal), intent(in), optional :: defaultVal
real(pReal) :: nodeAsFloat
class(tNode), pointer :: node
type(tScalar), pointer :: scalar
if (self%contains(k)) then
node => self%get(k)
select type(node)
select type(node => self%get(k))
class is(tScalar)
scalar => node%asScalar()
nodeAsFloat = scalar%asFloat()
class default
call IO_error(706,ext_msg='Expected scalar float for key '//k)
endselect
end select
elseif (present(defaultVal)) then
nodeAsFloat = defaultVal
else
@ -744,23 +776,22 @@ end function tNode_get_byKey_asFloat
!--------------------------------------------------------------------------------------------------
function tNode_get_byKey_asInt(self,k,defaultVal) result(nodeAsInt)
class(tNode), intent(in), target :: self
character(len=*), intent(in) :: k
integer, intent(in),optional :: defaultVal
class(tNode), intent(in) :: self
character(len=*), intent(in) :: k
integer, intent(in), optional :: defaultVal
integer :: nodeAsInt
class(tNode), pointer :: node
type(tScalar), pointer :: scalar
if (self%contains(k)) then
node => self%get(k)
select type(node)
select type(node => self%get(k))
class is(tScalar)
scalar => node%asScalar()
nodeAsInt = scalar%asInt()
class default
call IO_error(706,ext_msg='Expected scalar integer for key '//k)
endselect
end select
elseif (present(defaultVal)) then
nodeAsInt = defaultVal
else
@ -775,23 +806,22 @@ end function tNode_get_byKey_asInt
!--------------------------------------------------------------------------------------------------
function tNode_get_byKey_asBool(self,k,defaultVal) result(nodeAsBool)
class(tNode), intent(in), target :: self
character(len=*), intent(in) :: k
logical, intent(in),optional :: defaultVal
class(tNode), intent(in) :: self
character(len=*), intent(in) :: k
logical, intent(in), optional :: defaultVal
logical :: nodeAsBool
class(tNode), pointer :: node
type(tScalar), pointer :: scalar
if (self%contains(k)) then
node => self%get(k)
select type(node)
select type(node => self%get(k))
class is(tScalar)
scalar => node%asScalar()
nodeAsBool = scalar%asBool()
class default
call IO_error(706,ext_msg='Expected scalar Boolean for key '//k)
endselect
end select
elseif (present(defaultVal)) then
nodeAsBool = defaultVal
else
@ -806,23 +836,22 @@ end function tNode_get_byKey_asBool
!--------------------------------------------------------------------------------------------------
function tNode_get_byKey_asString(self,k,defaultVal) result(nodeAsString)
class(tNode), intent(in), target :: self
character(len=*), intent(in) :: k
character(len=*), intent(in),optional :: defaultVal
class(tNode), intent(in) :: self
character(len=*), intent(in) :: k
character(len=*), intent(in), optional :: defaultVal
character(len=:), allocatable :: nodeAsString
class(tNode), pointer :: node
type(tScalar), pointer :: scalar
if (self%contains(k)) then
node => self%get(k)
select type(node)
select type(node => self%get(k))
class is(tScalar)
scalar => node%asScalar()
nodeAsString = scalar%asString()
class default
call IO_error(706,ext_msg='Expected scalar string for key '//k)
endselect
end select
elseif (present(defaultVal)) then
nodeAsString = defaultVal
else
@ -837,25 +866,24 @@ end function tNode_get_byKey_asString
!--------------------------------------------------------------------------------------------------
function tNode_get_byKey_as1dFloat(self,k,defaultVal,requiredSize) result(nodeAs1dFloat)
class(tNode), intent(in), target :: self
character(len=*), intent(in) :: k
class(tNode), intent(in) :: self
character(len=*), intent(in) :: k
real(pReal), intent(in), dimension(:), optional :: defaultVal
integer, intent(in), optional :: requiredSize
real(pReal), dimension(:), allocatable :: nodeAs1dFloat
class(tNode), pointer :: node
type(tList), pointer :: list
if (self%contains(k)) then
node => self%get(k)
select type(node)
select type(node => self%get(k))
class is(tList)
list => node%asList()
nodeAs1dFloat = list%as1dFloat()
class default
call IO_error(706,ext_msg='Expected 1D float array for key '//k)
endselect
end select
elseif (present(defaultVal)) then
nodeAs1dFloat = defaultVal
else
@ -874,25 +902,24 @@ end function tNode_get_byKey_as1dFloat
!--------------------------------------------------------------------------------------------------
function tNode_get_byKey_as2dFloat(self,k,defaultVal,requiredShape) result(nodeAs2dFloat)
class(tNode), intent(in), target :: self
character(len=*), intent(in) :: k
class(tNode), intent(in) :: self
character(len=*), intent(in) :: k
real(pReal), intent(in), dimension(:,:), optional :: defaultVal
integer, intent(in), dimension(2), optional :: requiredShape
real(pReal), dimension(:,:), allocatable :: nodeAs2dFloat
class(tNode), pointer :: node
type(tList), pointer :: rows
if(self%contains(k)) then
node => self%get(k)
select type(node)
select type(node => self%get(k))
class is(tList)
rows => node%asList()
nodeAs2dFloat = rows%as2dFloat()
class default
call IO_error(706,ext_msg='Expected 2D float array for key '//k)
endselect
end select
elseif(present(defaultVal)) then
nodeAs2dFloat = defaultVal
else
@ -911,24 +938,22 @@ end function tNode_get_byKey_as2dFloat
!--------------------------------------------------------------------------------------------------
function tNode_get_byKey_as1dInt(self,k,defaultVal,requiredSize) result(nodeAs1dInt)
class(tNode), intent(in), target :: self
character(len=*), intent(in) :: k
class(tNode), intent(in) :: self
character(len=*), intent(in) :: k
integer, dimension(:), intent(in), optional :: defaultVal
integer, intent(in), optional :: requiredSize
integer, dimension(:), allocatable :: nodeAs1dInt
class(tNode), pointer :: node
type(tList), pointer :: list
if (self%contains(k)) then
node => self%get(k)
select type(node)
select type(node => self%get(k))
class is(tList)
list => node%asList()
nodeAs1dInt = list%as1dInt()
class default
call IO_error(706,ext_msg='Expected 1D integer array for key '//k)
endselect
end select
elseif (present(defaultVal)) then
nodeAs1dInt = defaultVal
else
@ -947,23 +972,22 @@ end function tNode_get_byKey_as1dInt
!--------------------------------------------------------------------------------------------------
function tNode_get_byKey_as1dBool(self,k,defaultVal) result(nodeAs1dBool)
class(tNode), intent(in), target :: self
character(len=*), intent(in) :: k
class(tNode), intent(in) :: self
character(len=*), intent(in) :: k
logical, dimension(:), intent(in), optional :: defaultVal
logical, dimension(:), allocatable :: nodeAs1dBool
class(tNode), pointer :: node
type(tList), pointer :: list
if (self%contains(k)) then
node => self%get(k)
select type(node)
select type(node => self%get(k))
class is(tList)
list => node%asList()
nodeAs1dBool = list%as1dBool()
class default
call IO_error(706,ext_msg='Expected 1D Boolean array for key '//k)
endselect
end select
elseif (present(defaultVal)) then
nodeAs1dBool = defaultVal
else
@ -978,23 +1002,22 @@ end function tNode_get_byKey_as1dBool
!--------------------------------------------------------------------------------------------------
function tNode_get_byKey_as1dString(self,k,defaultVal) result(nodeAs1dString)
class(tNode), intent(in), target :: self
character(len=*), intent(in) :: k
class(tNode), intent(in) :: self
character(len=*), intent(in) :: k
character(len=*), intent(in), dimension(:), optional :: defaultVal
character(len=:), allocatable, dimension(:) :: nodeAs1dString
class(tNode), pointer :: node
type(tList), pointer :: list
if (self%contains(k)) then
node => self%get(k)
select type(node)
select type(node => self%get(k))
class is(tList)
list => node%asList()
nodeAs1dString = list%as1dString()
class default
call IO_error(706,ext_msg='Expected 1D string array for key '//k)
endselect
end select
elseif (present(defaultVal)) then
nodeAs1dString = defaultVal
else
@ -1015,11 +1038,11 @@ function output_as1dString(self) result(output) !ToDo: SR: Re
class(tNode), pointer :: output_list
integer :: o
output_list => self%get('output',defaultVal=emptyList)
output_list => self%get('output',defaultVal=emptyList)
allocate(output(output_list%length))
do o = 1, output_list%length
output(o) = output_list%get_asString(o)
enddo
end do
end function output_as1dString
@ -1042,7 +1065,7 @@ function tNode_get_byKey_asIndex(self,key) result(keyIndex)
do while (associated(item%next) .and. item%key /= key)
item => item%next
keyIndex = keyIndex+1
enddo
end do
if (item%key /= key) call IO_error(140,ext_msg=key)
@ -1087,7 +1110,7 @@ recursive function tList_asFormattedString(self,indent) result(str)
if (i /= 1) str = str//repeat(' ',indent_)
str = str//'- '//item%node%asFormattedString(indent_+2)
item => item%next
enddo
end do
end function tList_asFormattedString
@ -1119,9 +1142,9 @@ recursive function tDict_asFormattedString(self,indent) result(str)
str = str//trim(item%key)//': '//item%node%asFormattedString(indent_+len_trim(item%key)+2)
class default
str = str//trim(item%key)//':'//IO_EOL//repeat(' ',indent_+2)//item%node%asFormattedString(indent_+2)
endselect
end select
item => item%next
enddo
end do
end function tDict_asFormattedString
@ -1190,13 +1213,14 @@ function tList_as1dFloat(self)
type(tItem), pointer :: item
type(tScalar), pointer :: scalar
allocate(tList_as1dFloat(self%length))
item => self%first
do i = 1, self%length
scalar => item%node%asScalar()
tList_as1dFloat(i) = scalar%asFloat()
item => item%next
enddo
end do
end function tList_as1dFloat
@ -1213,6 +1237,7 @@ function tList_as2dFloat(self)
class(tNode), pointer :: row
type(tList), pointer :: row_data
row => self%get(1)
row_data => row%asList()
allocate(tList_as2dFloat(self%length,row_data%length))
@ -1220,9 +1245,9 @@ function tList_as2dFloat(self)
do i=1,self%length
row => self%get(i)
row_data => row%asList()
if(row_data%length /= size(tList_as2dFloat,2)) call IO_error(709,ext_msg='Varying number of columns')
if (row_data%length /= size(tList_as2dFloat,2)) call IO_error(709,ext_msg='Varying number of columns')
tList_as2dFloat(i,:) = self%get_as1dFloat(i)
enddo
end do
end function tList_as2dFloat
@ -1239,13 +1264,14 @@ function tList_as1dInt(self)
type(tItem), pointer :: item
type(tScalar), pointer :: scalar
allocate(tList_as1dInt(self%length))
item => self%first
do i = 1, self%length
scalar => item%node%asScalar()
tList_as1dInt(i) = scalar%asInt()
item => item%next
enddo
end do
end function tList_as1dInt
@ -1262,13 +1288,14 @@ function tList_as1dBool(self)
type(tItem), pointer :: item
type(tScalar), pointer :: scalar
allocate(tList_as1dBool(self%length))
item => self%first
do i = 1, self%length
scalar => item%node%asScalar()
tList_as1dBool(i) = scalar%asBool()
item => item%next
enddo
end do
end function tList_as1dBool
@ -1285,13 +1312,14 @@ function tList_as1dString(self)
type(tItem), pointer :: item
type(tScalar), pointer :: scalar
len_max = 0
item => self%first
do i = 1, self%length
scalar => item%node%asScalar()
len_max = max(len_max, len_trim(scalar%asString()))
item => item%next
enddo
end do
allocate(character(len=len_max) :: tList_as1dString(self%length))
item => self%first

60
src/grid/spectral_utilities.f90
View File

@ -156,7 +156,7 @@ subroutine spectral_utilities_init
integer(C_INTPTR_T) :: alloc_local, local_K, local_K_offset
integer(C_INTPTR_T), parameter :: &
scalarSize = 1_C_INTPTR_T, &
vecSize = 3_C_INTPTR_T, &
vectorSize = 3_C_INTPTR_T, &
tensorSize = 9_C_INTPTR_T
character(len=*), parameter :: &
PETSCDEBUG = ' -snes_view -snes_monitor '
@ -274,7 +274,7 @@ subroutine spectral_utilities_init
call c_f_pointer(tensorField, tensorField_fourier, [3_C_INTPTR_T,3_C_INTPTR_T, &
gridFFTW(1)/2_C_INTPTR_T + 1_C_INTPTR_T , gridFFTW(2),local_K]) ! place a pointer for a fourier tensor representation
vectorField = fftw_alloc_complex(vecSize*alloc_local)
vectorField = fftw_alloc_complex(vectorSize*alloc_local)
call c_f_pointer(vectorField, vectorField_real, [3_C_INTPTR_T,&
2_C_INTPTR_T*(gridFFTW(1)/2_C_INTPTR_T + 1_C_INTPTR_T),gridFFTW(2),local_K]) ! place a pointer for a real vector representation
call c_f_pointer(vectorField, vectorField_fourier,[3_C_INTPTR_T,&
@ -288,42 +288,42 @@ subroutine spectral_utilities_init
!--------------------------------------------------------------------------------------------------
! tensor MPI fftw plans
planTensorForth = fftw_mpi_plan_many_dft_r2c(3, [gridFFTW(3),gridFFTW(2),gridFFTW(1)], & ! dimension, logical length in each dimension in reversed order
tensorSize, FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, &! no. of transforms, default iblock and oblock
tensorField_real, tensorField_fourier, & ! input data, output data
PETSC_COMM_WORLD, FFTW_planner_flag) ! use all processors, planer precision
planTensorForth = fftw_mpi_plan_many_dft_r2c(3,gridFFTW(3:1:-1),tensorSize, &
FFTW_MPI_DEFAULT_BLOCK,FFTW_MPI_DEFAULT_BLOCK, &
tensorField_real,tensorField_fourier, &
PETSC_COMM_WORLD,FFTW_planner_flag)
if (.not. c_associated(planTensorForth)) error stop 'FFTW error'
planTensorBack = fftw_mpi_plan_many_dft_c2r(3, [gridFFTW(3),gridFFTW(2),gridFFTW(1)], & ! dimension, logical length in each dimension in reversed order
tensorSize, FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, &! no. of transforms, default iblock and oblock
tensorField_fourier,tensorField_real, & ! input data, output data
PETSC_COMM_WORLD, FFTW_planner_flag) ! all processors, planer precision
planTensorBack = fftw_mpi_plan_many_dft_c2r(3,gridFFTW(3:1:-1),tensorSize, &
FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, &
tensorField_fourier,tensorField_real, &
PETSC_COMM_WORLD, FFTW_planner_flag)
if (.not. c_associated(planTensorBack)) error stop 'FFTW error'
!--------------------------------------------------------------------------------------------------
! vector MPI fftw plans
planVectorForth = fftw_mpi_plan_many_dft_r2c(3, [gridFFTW(3),gridFFTW(2),gridFFTW(1)], & ! dimension, logical length in each dimension in reversed order
vecSize, FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK,&! no. of transforms, default iblock and oblock
vectorField_real, vectorField_fourier, & ! input data, output data
PETSC_COMM_WORLD, FFTW_planner_flag) ! use all processors, planer precision
if (.not. C_ASSOCIATED(planVectorForth)) error stop 'FFTW error'
planVectorBack = fftw_mpi_plan_many_dft_c2r(3, [gridFFTW(3),gridFFTW(2),gridFFTW(1)], & ! dimension, logical length in each dimension in reversed order
vecSize, FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, & ! no. of transforms, default iblock and oblock
vectorField_fourier,vectorField_real, & ! input data, output data
PETSC_COMM_WORLD, FFTW_planner_flag) ! all processors, planer precision
if (.not. C_ASSOCIATED(planVectorBack)) error stop 'FFTW error'
planVectorForth = fftw_mpi_plan_many_dft_r2c(3,gridFFTW(3:1:-1),vectorSize, &
FFTW_MPI_DEFAULT_BLOCK,FFTW_MPI_DEFAULT_BLOCK, &
vectorField_real,vectorField_fourier, &
PETSC_COMM_WORLD,FFTW_planner_flag)
if (.not. c_associated(planVectorForth)) error stop 'FFTW error'
planVectorBack = fftw_mpi_plan_many_dft_c2r(3,gridFFTW(3:1:-1),vectorSize, &
FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, &
vectorField_fourier,vectorField_real, &
PETSC_COMM_WORLD, FFTW_planner_flag)
if (.not. c_associated(planVectorBack)) error stop 'FFTW error'
!--------------------------------------------------------------------------------------------------
! scalar MPI fftw plans
planScalarForth = fftw_mpi_plan_many_dft_r2c(3, [gridFFTW(3),gridFFTW(2),gridFFTW(1)], & ! dimension, logical length in each dimension in reversed order
scalarSize, FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, &! no. of transforms, default iblock and oblock
scalarField_real, scalarField_fourier, & ! input data, output data
PETSC_COMM_WORLD, FFTW_planner_flag) ! use all processors, planer precision
if (.not. C_ASSOCIATED(planScalarForth)) error stop 'FFTW error'
planScalarBack = fftw_mpi_plan_many_dft_c2r(3, [gridFFTW(3),gridFFTW(2),gridFFTW(1)], & ! dimension, logical length in each dimension in reversed order, no. of transforms
scalarSize, FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, &! no. of transforms, default iblock and oblock
scalarField_fourier,scalarField_real, & ! input data, output data
PETSC_COMM_WORLD, FFTW_planner_flag) ! use all processors, planer precision
if (.not. C_ASSOCIATED(planScalarBack)) error stop 'FFTW error'
planScalarForth = fftw_mpi_plan_many_dft_r2c(3,gridFFTW(3:1:-1),scalarSize, &
FFTW_MPI_DEFAULT_BLOCK,FFTW_MPI_DEFAULT_BLOCK, &
scalarField_real,scalarField_fourier, &
PETSC_COMM_WORLD,FFTW_planner_flag)
if (.not. c_associated(planScalarForth)) error stop 'FFTW error'
planScalarBack = fftw_mpi_plan_many_dft_c2r(3,gridFFTW(3:1:-1),scalarSize, &
FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, &
scalarField_fourier,scalarField_real, &
PETSC_COMM_WORLD, FFTW_planner_flag)
if (.not. c_associated(planScalarBack)) error stop 'FFTW error'
!--------------------------------------------------------------------------------------------------
! calculation of discrete angular frequencies, ordered as in FFTW (wrap around)

142
src/lattice.f90
View File

@ -587,8 +587,8 @@ function lattice_C66_trans(Ntrans,C_parent66,lattice_target, &
!--------------------------------------------------------------------------------------------------
!> @brief Non-schmid projections for bcc with up to 6 coefficients
! Koester et al. 2012, Acta Materialia 60 (2012) 3894–3901, eq. (17)
! Gröger et al. 2008, Acta Materialia 56 (2008) 5412–5425, table 1
! https://doi.org/10.1016/j.actamat.2012.03.053, eq. (17)
! https://doi.org/10.1016/j.actamat.2008.07.037, table 1
!--------------------------------------------------------------------------------------------------
function lattice_nonSchmidMatrix(Nslip,nonSchmidCoefficients,sense) result(nonSchmidMatrix)
@ -602,6 +602,7 @@ function lattice_nonSchmidMatrix(Nslip,nonSchmidCoefficients,sense) result(nonSc
type(rotation) :: R
integer :: i
if (abs(sense) /= 1) error stop 'Sense in lattice_nonSchmidMatrix'
coordinateSystem = buildCoordinateSystem(Nslip,BCC_NSLIPSYSTEM,BCC_SYSTEMSLIP,'cI',0.0_pReal)
@ -634,7 +635,9 @@ end function lattice_nonSchmidMatrix
!--------------------------------------------------------------------------------------------------
!> @brief Slip-slip interaction matrix
!> details only active slip systems are considered
!> @details only active slip systems are considered
!> @details https://doi.org/10.1016/j.actamat.2016.12.040 (fcc: Tab S4-1, bcc: Tab S5-1)
!> @details https://doi.org/10.1016/j.ijplas.2014.06.010 (hex: Tab 3b)
!--------------------------------------------------------------------------------------------------
function lattice_interaction_SlipBySlip(Nslip,interactionValues,lattice) result(interactionMatrix)
@ -646,6 +649,7 @@ function lattice_interaction_SlipBySlip(Nslip,interactionValues,lattice) result(
integer, dimension(:), allocatable :: NslipMax
integer, dimension(:,:), allocatable :: interactionTypes
integer, dimension(FCC_NSLIP,FCC_NSLIP), parameter :: &
FCC_INTERACTIONSLIPSLIP = reshape( [&
1, 2, 2, 4, 7, 5, 3, 5, 5, 4, 6, 7, 10,11,10,11,12,13, & ! -----> acting (forest)
@ -750,41 +754,113 @@ function lattice_interaction_SlipBySlip(Nslip,interactionValues,lattice) result(
integer, dimension(HEX_NSLIP,HEX_NSLIP), parameter :: &
HEX_INTERACTIONSLIPSLIP = reshape( [&
! basal prism 1. pyr<a> 1. pyr<c+a> 2. pyr<c+a>
1, 2, 2, 3, 3, 3, 7, 7, 7, 7, 7, 7, 13,13,13,13,13,13,13,13,13,13,13,13, 21,21,21,21,21,21, & ! -----> acting (forest)
2, 1, 2, 3, 3, 3, 7, 7, 7, 7, 7, 7, 13,13,13,13,13,13,13,13,13,13,13,13, 21,21,21,21,21,21, & ! | basal
2, 2, 1, 3, 3, 3, 7, 7, 7, 7, 7, 7, 13,13,13,13,13,13,13,13,13,13,13,13, 21,21,21,21,21,21, & ! |
1, 2, 2, 3, 4, 4, 9,10, 9, 9,10, 9, 20,21,22,22,21,20,20,21,22,22,21,20, 47,47,48,47,47,48, & ! -----> acting (forest)
2, 1, 2, 4, 3, 4, 10, 9, 9,10, 9, 9, 22,22,21,20,20,21,22,22,21,20,20,21, 47,48,47,47,48,47, & ! | basal
2, 2, 1, 4, 4, 3, 9, 9,10, 9, 9,10, 21,20,20,21,22,22,21,20,20,21,22,22, 48,47,47,48,47,47, & ! |
! v
6, 6, 6, 4, 5, 5, 8, 8, 8, 8, 8, 8, 14,14,14,14,14,14,14,14,14,14,14,14, 22,22,22,22,22,22, & ! reacting (primary)
6, 6, 6, 5, 4, 5, 8, 8, 8, 8, 8, 8, 14,14,14,14,14,14,14,14,14,14,14,14, 22,22,22,22,22,22, & ! prism
6, 6, 6, 5, 5, 4, 8, 8, 8, 8, 8, 8, 14,14,14,14,14,14,14,14,14,14,14,14, 22,22,22,22,22,22, &
7, 8, 8, 5, 6, 6, 11,12,11,11,12,11, 23,24,25,25,24,23,23,24,25,25,24,23, 49,49,50,49,49,50, & ! reacting (primary)
8, 7, 8, 6, 5, 6, 12,11,11,12,11,11, 25,25,24,23,23,24,25,25,24,23,23,24, 49,50,49,49,50,49, & ! prism
8, 8, 7, 6, 6, 5, 11,11,12,11,11,12, 24,23,23,24,25,25,24,23,23,24,25,25, 50,49,49,50,49,49, &
12,12,12, 11,11,11, 9,10,10,10,10,10, 15,15,15,15,15,15,15,15,15,15,15,15, 23,23,23,23,23,23, &
12,12,12, 11,11,11, 10, 9,10,10,10,10, 15,15,15,15,15,15,15,15,15,15,15,15, 23,23,23,23,23,23, &
12,12,12, 11,11,11, 10,10, 9,10,10,10, 15,15,15,15,15,15,15,15,15,15,15,15, 23,23,23,23,23,23, &
12,12,12, 11,11,11, 10,10,10, 9,10,10, 15,15,15,15,15,15,15,15,15,15,15,15, 23,23,23,23,23,23, & ! 1. pyr<a>
12,12,12, 11,11,11, 10,10,10,10, 9,10, 15,15,15,15,15,15,15,15,15,15,15,15, 23,23,23,23,23,23, &
12,12,12, 11,11,11, 10,10,10,10,10, 9, 15,15,15,15,15,15,15,15,15,15,15,15, 23,23,23,23,23,23, &
18,19,18, 16,17,16, 13,14,14,15,14,14, 26,26,27,28,28,27,29,29,27,28,28,27, 51,52,51,51,52,51, &
19,18,18, 17,16,16, 14,13,14,14,15,14, 28,27,26,26,27,28,28,27,29,29,27,28, 51,51,52,51,51,52, &
18,18,19, 16,16,17, 14,14,13,14,14,15, 27,28,28,27,26,26,27,28,28,27,29,29, 52,51,51,52,51,51, &
18,19,18, 16,17,16, 15,14,14,13,14,14, 29,29,27,28,28,27,26,26,27,28,28,27, 51,52,51,51,52,51, & ! 1. pyr<a>
19,18,18, 17,16,16, 14,15,14,14,13,14, 28,27,29,29,27,28,28,27,26,26,27,28, 51,51,52,51,51,52, &
18,18,19, 16,16,17, 14,14,15,14,14,13, 27,28,28,27,29,29,27,28,28,27,26,26, 52,51,51,52,51,51, &
20,20,20, 19,19,19, 18,18,18,18,18,18, 16,17,17,17,17,17,17,17,17,17,17,17, 24,24,24,24,24,24, &
20,20,20, 19,19,19, 18,18,18,18,18,18, 17,16,17,17,17,17,17,17,17,17,17,17, 24,24,24,24,24,24, &
20,20,20, 19,19,19, 18,18,18,18,18,18, 17,17,16,17,17,17,17,17,17,17,17,17, 24,24,24,24,24,24, &
20,20,20, 19,19,19, 18,18,18,18,18,18, 17,17,17,16,17,17,17,17,17,17,17,17, 24,24,24,24,24,24, &
20,20,20, 19,19,19, 18,18,18,18,18,18, 17,17,17,17,16,17,17,17,17,17,17,17, 24,24,24,24,24,24, &
20,20,20, 19,19,19, 18,18,18,18,18,18, 17,17,17,17,17,16,17,17,17,17,17,17, 24,24,24,24,24,24, &
20,20,20, 19,19,19, 18,18,18,18,18,18, 17,17,17,17,17,17,16,17,17,17,17,17, 24,24,24,24,24,24, & ! 1. pyr<c+a>
20,20,20, 19,19,19, 18,18,18,18,18,18, 17,17,17,17,17,17,17,16,17,17,17,17, 24,24,24,24,24,24, &
20,20,20, 19,19,19, 18,18,18,18,18,18, 17,17,17,17,17,17,17,17,16,17,17,17, 24,24,24,24,24,24, &
20,20,20, 19,19,19, 18,18,18,18,18,18, 17,17,17,17,17,17,17,17,17,16,17,17, 24,24,24,24,24,24, &
20,20,20, 19,19,19, 18,18,18,18,18,18, 17,17,17,17,17,17,17,17,17,17,16,17, 24,24,24,24,24,24, &
20,20,20, 19,19,19, 18,18,18,18,18,18, 17,17,17,17,17,17,17,17,17,17,17,16, 24,24,24,24,24,24, &
44,45,46, 41,42,43, 37,38,39,40,38,39, 30,31,32,32,32,33,34,35,32,32,32,36, 53,54,55,53,54,56, &
46,45,44, 43,42,41, 37,39,38,40,39,38, 31,30,36,32,32,32,35,34,33,32,32,32, 56,54,53,55,54,53, &
45,46,44, 42,43,41, 39,37,38,39,40,38, 32,36,30,31,32,32,32,33,34,35,32,32, 56,53,54,55,53,54, &
45,44,46, 42,41,43, 38,37,39,38,40,39, 32,32,31,30,36,32,32,32,35,34,33,32, 53,56,54,53,55,54, &
46,44,45, 43,41,42, 38,39,37,38,39,40, 32,32,32,36,30,31,32,32,32,33,34,35, 54,56,53,54,55,53, &
44,46,45, 41,43,42, 39,38,37,39,38,40, 33,32,32,32,31,30,36,32,32,32,35,34, 54,53,56,54,53,55, &
44,45,46, 41,42,43, 40,38,39,37,38,39, 34,35,32,32,32,36,30,31,32,32,32,33, 53,54,56,53,54,55, & ! 1. pyr<c+a>
46,45,44, 43,42,41, 40,39,38,37,39,38, 35,34,33,32,32,32,31,30,36,32,32,32, 55,54,53,56,54,53, &
45,46,44, 42,43,41, 39,40,38,39,37,38, 32,33,34,35,32,32,32,36,30,31,32,32, 55,53,54,56,53,54, &
45,44,46, 42,41,43, 38,40,39,38,37,39, 32,32,35,34,33,32,32,32,31,30,36,32, 53,55,54,53,56,54, &
46,44,45, 43,41,42, 38,39,40,38,39,37, 32,32,32,33,34,35,32,32,32,36,30,31, 54,55,53,54,56,53, &
44,46,45, 41,43,42, 39,38,40,39,38,37, 36,32,32,32,35,34,33,32,32,32,31,30, 54,53,55,54,53,56, &
30,30,30, 29,29,29, 28,28,28,28,28,28, 27,27,27,27,27,27,27,27,27,27,27,27, 25,26,26,26,26,26, &
30,30,30, 29,29,29, 28,28,28,28,28,28, 27,27,27,27,27,27,27,27,27,27,27,27, 26,25,26,26,26,26, &
30,30,30, 29,29,29, 28,28,28,28,28,28, 27,27,27,27,27,27,27,27,27,27,27,27, 26,26,25,26,26,26, &
30,30,30, 29,29,29, 28,28,28,28,28,28, 27,27,27,27,27,27,27,27,27,27,27,27, 26,26,26,25,26,26, & ! 2. pyr<c+a>
30,30,30, 29,29,29, 28,28,28,28,28,28, 27,27,27,27,27,27,27,27,27,27,27,27, 26,26,26,26,25,26, &
30,30,30, 29,29,29, 28,28,28,28,28,28, 27,27,27,27,27,27,27,27,27,27,27,27, 26,26,26,26,26,25 &
68,68,69, 66,66,67, 64,64,65,64,65,65, 60,61,61,60,62,62,60,63,63,60,62,62, 57,58,58,59,58,58, &
68,69,68, 66,67,66, 65,64,64,65,64,64, 62,62,60,61,61,60,62,62,60,63,63,60, 58,57,58,58,59,58, &
69,68,68, 67,66,66, 64,65,64,64,65,64, 63,60,62,62,60,61,61,60,62,62,60,63, 58,58,57,58,58,59, &
68,68,69, 66,66,67, 64,64,65,64,64,65, 60,63,63,60,62,62,60,61,61,60,62,62, 59,58,58,57,58,58, & ! 2. pyr<c+a>
68,69,68, 66,67,66, 65,64,64,65,64,64, 62,62,60,63,63,60,62,62,60,61,61,60, 58,59,58,58,57,58, &
69,68,68, 67,66,66, 64,65,64,64,65,64, 61,60,62,62,60,63,63,60,62,62,60,61, 58,58,59,58,58,57 &
],shape(HEX_INTERACTIONSLIPSLIP)) !< Slip-slip interaction types for hex (onion peel naming scheme)
!< 10.1016/j.ijplas.2014.06.010 table 3
!< 10.1080/14786435.2012.699689 table 2 and 3
!< index & label & description
!< 1 & S1 & basal self-interaction
!< 2 & 1 & basal/basal coplanar
!< 3 & 3 & basal/prismatic collinear
!< 4 & 4 & basal/prismatic non-collinear
!< 5 & S2 & prismatic self-interaction
!< 6 & 2 & prismatic/prismatic
!< 7 & 5 & prismatic/basal collinear
!< 8 & 6 & prismatic/basal non-collinear
!< 9 & - & basal/pyramidal <a> non-collinear
!< 10 & - & basal/pyramidal <a> collinear
!< 11 & - & prismatic/pyramidal <a> non-collinear
!< 12 & - & prismatic/pyramidal <a> collinear
!< 13 & - & pyramidal <a> self-interaction
!< 14 & - & pyramidal <a> non-collinear
!< 15 & - & pyramidal <a> collinear
!< 16 & - & pyramidal <a>/prismatic non-collinear
!< 17 & - & pyramidal <a>/prismatic collinear
!< 18 & - & pyramidal <a>/basal non-collinear
!< 19 & - & pyramidal <a>/basal collinear
!< 20 & - & basal/1. order pyramidal <c+a> semi-collinear
!< 21 & - & basal/1. order pyramidal <c+a>
!< 22 & - & basal/1. order pyramidal <c+a>
!< 23 & - & prismatic/1. order pyramidal <c+a> semi-collinear
!< 24 & - & prismatic/1. order pyramidal <c+a>
!< 25 & - & prismatic/1. order pyramidal <c+a> semi-coplanar?
!< 26 & - & pyramidal <a>/1. order pyramidal <c+a> coplanar
!< 27 & - & pyramidal <a>/1. order pyramidal <c+a>
!< 28 & - & pyramidal <a>/1. order pyramidal <c+a> semi-collinear
!< 29 & - & pyramidal <a>/1. order pyramidal <c+a> semi-coplanar
!< 30 & - & 1. order pyramidal <c+a> self-interaction
!< 31 & - & 1. order pyramidal <c+a> coplanar
!< 32 & - & 1. order pyramidal <c+a>
!< 33 & - & 1. order pyramidal <c+a>
!< 34 & - & 1. order pyramidal <c+a> semi-coplanar
!< 35 & - & 1. order pyramidal <c+a> semi-coplanar
!< 36 & - & 1. order pyramidal <c+a> collinear
!< 37 & - & 1. order pyramidal <c+a>/pyramidal <a> coplanar
!< 38 & - & 1. order pyramidal <c+a>/pyramidal <a> semi-collinear
!< 39 & - & 1. order pyramidal <c+a>/pyramidal <a>
!< 40 & - & 1. order pyramidal <c+a>/pyramidal <a> semi-coplanar
!< 41 & - & 1. order pyramidal <c+a>/prismatic semi-collinear
!< 42 & - & 1. order pyramidal <c+a>/prismatic semi-coplanar
!< 43 & - & 1. order pyramidal <c+a>/prismatic
!< 44 & - & 1. order pyramidal <c+a>/basal semi-collinear
!< 45 & - & 1. order pyramidal <c+a>/basal
!< 46 & - & 1. order pyramidal <c+a>/basal
!< 47 & 8 & basal/2. order pyramidal <c+a> non-collinear
!< 48 & 7 & basal/2. order pyramidal <c+a> semi-collinear
!< 49 & 10 & prismatic/2. order pyramidal <c+a>
!< 50 & 9 & prismatic/2. order pyramidal <c+a> semi-collinear
!< 51 & - & pyramidal <a>/2. order pyramidal <c+a>
!< 52 & - & pyramidal <a>/2. order pyramidal <c+a> semi collinear
!< 53 & - & 1. order pyramidal <c+a>/2. order pyramidal <c+a>
!< 54 & - & 1. order pyramidal <c+a>/2. order pyramidal <c+a>
!< 55 & - & 1. order pyramidal <c+a>/2. order pyramidal <c+a>
!< 56 & - & 1. order pyramidal <c+a>/2. order pyramidal <c+a> collinear
!< 57 & S3 & 2. order pyramidal <c+a> self-interaction
!< 58 & 16 & 2. order pyramidal <c+a> non-collinear
!< 59 & 15 & 2. order pyramidal <c+a> semi-collinear
!< 60 & - & 2. order pyramidal <c+a>/1. order pyramidal <c+a>
!< 61 & - & 2. order pyramidal <c+a>/1. order pyramidal <c+a> collinear
!< 62 & - & 2. order pyramidal <c+a>/1. order pyramidal <c+a>
!< 63 & - & 2. order pyramidal <c+a>/1. order pyramidal <c+a>
!< 64 & - & 2. order pyramidal <c+a>/pyramidal <a> non-collinear
!< 65 & - & 2. order pyramidal <c+a>/pyramidal <a> semi-collinear
!< 66 & 14 & 2. order pyramidal <c+a>/prismatic non-collinear
!< 67 & 13 & 2. order pyramidal <c+a>/prismatic semi-collinear
!< 68 & 12 & 2. order pyramidal <c+a>/basal non-collinear
!< 69 & 11 & 2. order pyramidal <c+a>/basal semi-collinear
integer, dimension(BCT_NSLIP,BCT_NSLIP), parameter :: &
BCT_INTERACTIONSLIPSLIP = reshape( [&

11
src/material.f90
View File

@ -108,8 +108,14 @@ subroutine parse()
homogenizations => config_material%get('homogenization')
call sanityCheck(materials, homogenizations)
#if defined (__GFORTRAN__)
material_name_phase = getKeys(phases)
material_name_homogenization = getKeys(homogenizations)
#else
material_name_phase = phases%Keys()
material_name_homogenization = homogenizations%Keys()
#endif
allocate(homogenization_Nconstituents(homogenizations%length))
do h=1, homogenizations%length
@ -203,9 +209,9 @@ subroutine sanityCheck(materials,homogenizations)
end subroutine sanityCheck
#if defined (__GFORTRAN__)
!--------------------------------------------------------------------------------------------------
!> @brief Get all keys from a dictionary
!> @brief %keys() is broken on gfortran
!--------------------------------------------------------------------------------------------------
function getKeys(dict)
@ -228,5 +234,6 @@ function getKeys(dict)
end do
end function getKeys
#endif
end module material

32
src/math.f90
View File

@ -84,38 +84,34 @@ contains
subroutine math_init
real(pReal), dimension(4) :: randTest
integer :: &
randSize, &
randomSeed !< fixed seeding for pseudo-random number generator, Default 0: use random seed
integer, dimension(:), allocatable :: randInit
integer :: randSize
integer, dimension(:), allocatable :: seed
class(tNode), pointer :: &
num_generic
print'(/,1x,a)', '<<<+- math init -+>>>'; flush(IO_STDOUT)
num_generic => config_numerics%get('generic',defaultVal=emptyDict)
randomSeed = num_generic%get_asInt('random_seed', defaultVal = 0)
call random_seed(size=randSize)
allocate(randInit(randSize))
if (randomSeed > 0) then
randInit = randomSeed
allocate(seed(randSize))
if (num_generic%contains('random_seed')) then
seed = num_generic%get_as1dInt('random_seed',requiredSize=randSize)
else
call random_seed()
call random_seed(get = randInit)
randInit(2:randSize) = randInit(1)
endif
call random_seed(get = seed)
end if
call random_seed(put = randInit)
call random_seed(put = seed)
call random_number(randTest)
print'(/,a,i2)', ' size of random seed: ', randSize
print'( a,i0)', ' value of random seed: ', randInit(1)
print'( a,4(/,26x,f17.14),/)', ' start of random sequence: ', randTest
print'(/,a,i2)', ' size of random seed: ', randSize
print*, 'value of random seed: ', seed
print'( a,4(/,26x,f17.14))', ' start of random sequence: ', randTest
call random_seed(put = randInit)
call selfTest
call selfTest()
end subroutine math_init

14
src/rotations.f90
View File

@ -270,7 +270,7 @@ pure elemental subroutine standardize(self)
class(rotation), intent(inout) :: self
if (self%q(1) < 0.0_pReal) self%q = - self%q
if (sign(1.0_pReal,self%q(1)) < 0.0_pReal) self%q = - self%q
end subroutine standardize
@ -450,7 +450,7 @@ pure function qu2om(qu) result(om)
om(3,2) = 2.0_pReal*(qu(4)*qu(3)+qu(1)*qu(2))
om(1,3) = 2.0_pReal*(qu(2)*qu(4)+qu(1)*qu(3))
if (P < 0.0_pReal) om = transpose(om)
if (sign(1.0_pReal,P) < 0.0_pReal) om = transpose(om)
end function qu2om
@ -480,7 +480,7 @@ pure function qu2eu(qu) result(eu)
atan2( 2.0_pReal*chi, q03-q12 ), &
atan2(( P*qu(1)*qu(3)+qu(2)*qu(4))*chi, (-P*qu(1)*qu(2)+qu(3)*qu(4))*chi )]
endif degenerated
where(eu<0.0_pReal) eu = mod(eu+2.0_pReal*PI,[2.0_pReal*PI,PI,2.0_pReal*PI])
where(sign(1.0_pReal,eu)<0.0_pReal) eu = mod(eu+2.0_pReal*PI,[2.0_pReal*PI,PI,2.0_pReal*PI])
end function qu2eu
@ -602,7 +602,7 @@ pure function om2qu(om) result(qu)
qu = [ (om(2,1) - om(1,2)) /s,(om(1,3) + om(3,1)) / s,(om(2,3) + om(3,2)) / s,0.25_pReal * s]
endif
endif
if(qu(1)<0._pReal) qu =-1.0_pReal * qu
if(sign(1.0_pReal,qu(1))<0.0_pReal) qu =-1.0_pReal * qu
qu = qu*[1.0_pReal,P,P,P]
end function om2qu
@ -628,7 +628,7 @@ pure function om2eu(om) result(eu)
eu = [atan2(om(1,2),om(1,1)), 0.5_pReal*PI*(1.0_pReal-om(3,3)),0.0_pReal ]
end if
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])
where(sign(1.0_pReal,eu)<0.0_pReal) eu = mod(eu+2.0_pReal*PI,[2.0_pReal*PI,PI,2.0_pReal*PI])
end function om2eu
@ -735,7 +735,7 @@ pure function eu2qu(eu) result(qu)
-P*sPhi*cos(ee(1)-ee(3)), &
-P*sPhi*sin(ee(1)-ee(3)), &
-P*cPhi*sin(ee(1)+ee(3))]
if(qu(1) < 0.0_pReal) qu = qu * (-1.0_pReal)
if(sign(1.0_pReal,qu(1)) < 0.0_pReal) qu = qu * (-1.0_pReal)
end function eu2qu
@ -792,7 +792,7 @@ pure function eu2ax(eu) result(ax)
else
ax(1:3) = -P/tau * [ t*cos(delta), t*sin(delta), sin(sigma) ] ! passive axis-angle pair so a minus sign in front
ax(4) = alpha
if (alpha < 0.0_pReal) ax = -ax ! ensure alpha is positive
if (sign(1.0_pReal,alpha) < 0.0_pReal) ax = -ax ! ensure alpha is positive
end if
end function eu2ax