diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index d82c7c8ac..e2f95c82d 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -506,7 +506,7 @@ Processing: - rm abq_addUserOutput.py marc_addUserOutput.py - $DAMASKROOT/PRIVATE/documenting/scriptHelpToWiki.py --debug *.py - cd $DAMASKROOT/processing/post - - rm marc_to_vtk.py vtk2ang.py DAD*.py + - rm vtk2ang.py DAD*.py - $DAMASKROOT/PRIVATE/documenting/scriptHelpToWiki.py --debug *.py except: - master diff --git a/VERSION b/VERSION index 4b8c1b29a..50693fa6c 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.3-252-g0e77be25 +v2.0.3-261-g99878952 diff --git a/processing/post/marc_to_vtk.py b/processing/post/marc_to_vtk.py deleted file mode 100755 index 05f7a6908..000000000 --- a/processing/post/marc_to_vtk.py +++ /dev/null @@ -1,199 +0,0 @@ -#!/usr/bin/env python2.7 -# -*- coding: UTF-8 no BOM -*- - -import os,sys,re -import argparse -import damask -import vtk, numpy as np - -scriptName = os.path.splitext(os.path.basename(__file__))[0] -scriptID = ' '.join([scriptName, damask.version]) - -parser = argparse.ArgumentParser(description='Convert from Marc input file format (.dat) to VTK format (.vtu)', version = scriptID) -parser.add_argument('filename', type=str, help='file to convert') -parser.add_argument('-t', '--table', type=str, help='ASCIItable file containing nodal data to subdivide and interpolate') - -args = parser.parse_args() - -with open(args.filename, 'r') as marcfile: - marctext = marcfile.read(); - -# Load table (if any) -if args.table is not None: - try: - table = damask.ASCIItable( - name=args.table, - outname='subdivided_{}'.format(args.table), - buffered=True - ) - - table.head_read() - table.data_readArray() - - # Python list is faster for appending - nodal_data = list(table.data) - except: args.table = None - -# Extract connectivity chunk from file... -connectivity_text = re.findall(r'connectivity[\n\r]+(.*?)[\n\r]+[a-zA-Z]', marctext, flags=(re.MULTILINE | re.DOTALL))[0] -connectivity_lines = re.split(r'[\n\r]+', connectivity_text, flags=(re.MULTILINE | re.DOTALL)) -connectivity_header = connectivity_lines[0] -connectivity_lines = connectivity_lines[1:] - -# Construct element map -elements = dict(map(lambda line: - ( - int(line[0:10]), # index - { - 'type': int(line[10:20]), - 'verts': list(map(int, re.split(r' +', line[20:].strip()))) - } - ), connectivity_lines)) - -# Extract coordinate chunk from file -coordinates_text = re.findall(r'coordinates[\n\r]+(.*?)[\n\r]+[a-zA-Z]', marctext, flags=(re.MULTILINE | re.DOTALL))[0] -coordinates_lines = re.split(r'[\n\r]+', coordinates_text, flags=(re.MULTILINE | re.DOTALL)) -coordinates_header = coordinates_lines[0] -coordinates_lines = coordinates_lines[1:] - -# marc input file does not use "e" in scientific notation, this adds it and converts -fl_format = lambda string: float(re.sub(r'(\d)([\+\-])', r'\1e\2', string)) -# Construct coordinate map -coordinates = dict(map(lambda line: - ( - int(line[0:10]), - np.array([ - fl_format(line[10:30]), - fl_format(line[30:50]), - fl_format(line[50:70]) - ]) - ), coordinates_lines)) - -# Subdivide volumes -grid = vtk.vtkUnstructuredGrid() -vertex_count = len(coordinates) -edge_to_vert = dict() # when edges are subdivided, a new vertex in the middle is produced and placed in here -ordered_pair = lambda a, b: (a, b) if a < b else (b, a) # edges are bidirectional - -def subdivide_edge(vert1, vert2): - edge = ordered_pair(vert1, vert2) - - if edge in edge_to_vert: - return edge_to_vert[edge] - - # Vertex does not exist, create it - newvert = len(coordinates) + 1 - coordinates[newvert] = 0.5 * (coordinates[vert1] + coordinates[vert2]) # Average - edge_to_vert[edge] = newvert; - - # Interpolate nodal data - if args.table is not None: - nodal_data.append(0.5 * (nodal_data[vert1 - 1] + nodal_data[vert2 - 1])) - return newvert; - -for el_id in range(1, len(elements) + 1): # Marc starts counting at 1 - el = elements[el_id] - if el['type'] == 7: - # Hexahedron, subdivided - - # There may be a better way to iterate over these, but this is consistent - # with the ordering scheme provided at https://damask.mpie.de/pub/Documentation/ElementType - - subverts = np.zeros((3,3,3), dtype=int) - # Get corners - subverts[0, 0, 0] = el['verts'][0] - subverts[2, 0, 0] = el['verts'][1] - subverts[2, 2, 0] = el['verts'][2] - subverts[0, 2, 0] = el['verts'][3] - subverts[0, 0, 2] = el['verts'][4] - subverts[2, 0, 2] = el['verts'][5] - subverts[2, 2, 2] = el['verts'][6] - subverts[0, 2, 2] = el['verts'][7] - - # lower edges - subverts[1, 0, 0] = subdivide_edge(subverts[0, 0, 0], subverts[2, 0, 0]) - subverts[2, 1, 0] = subdivide_edge(subverts[2, 0, 0], subverts[2, 2, 0]) - subverts[1, 2, 0] = subdivide_edge(subverts[2, 2, 0], subverts[0, 2, 0]) - subverts[0, 1, 0] = subdivide_edge(subverts[0, 2, 0], subverts[0, 0, 0]) - - # middle edges - subverts[0, 0, 1] = subdivide_edge(subverts[0, 0, 0], subverts[0, 0, 2]) - subverts[2, 0, 1] = subdivide_edge(subverts[2, 0, 0], subverts[2, 0, 2]) - subverts[2, 2, 1] = subdivide_edge(subverts[2, 2, 0], subverts[2, 2, 2]) - subverts[0, 2, 1] = subdivide_edge(subverts[0, 2, 0], subverts[0, 2, 2]) - - # top edges - subverts[1, 0, 2] = subdivide_edge(subverts[0, 0, 2], subverts[2, 0, 2]) - subverts[2, 1, 2] = subdivide_edge(subverts[2, 0, 2], subverts[2, 2, 2]) - subverts[1, 2, 2] = subdivide_edge(subverts[2, 2, 2], subverts[0, 2, 2]) - subverts[0, 1, 2] = subdivide_edge(subverts[0, 2, 2], subverts[0, 0, 2]) - - # then faces... The edge_to_vert addition is due to there being two ways - # to calculate a face vertex, depending on which opposite vertices are used to subdivide. - # This way, we avoid creating duplicate vertices. - subverts[1, 1, 0] = subdivide_edge(subverts[1, 0, 0], subverts[1, 2, 0]) - edge_to_vert[ordered_pair(subverts[0, 1, 0], subverts[2, 1, 0])] = subverts[1, 1, 0] - - subverts[1, 0, 1] = subdivide_edge(subverts[1, 0, 0], subverts[1, 0, 2]) - edge_to_vert[ordered_pair(subverts[0, 0, 1], subverts[2, 0, 1])] = subverts[1, 0, 1] - - subverts[2, 1, 1] = subdivide_edge(subverts[2, 1, 0], subverts[2, 1, 2]) - edge_to_vert[ordered_pair(subverts[2, 0, 1], subverts[2, 2, 1])] = subverts[2, 1, 1] - - subverts[1, 2, 1] = subdivide_edge(subverts[1, 2, 0], subverts[1, 2, 2]) - edge_to_vert[ordered_pair(subverts[0, 2, 1], subverts[2, 2, 1])] = subverts[1, 2, 1] - - subverts[0, 1, 1] = subdivide_edge(subverts[0, 1, 0], subverts[0, 1, 2]) - edge_to_vert[ordered_pair(subverts[0, 0, 1], subverts[0, 2, 1])] = subverts[0, 1, 1] - - subverts[1, 1, 2] = subdivide_edge(subverts[1, 0, 2], subverts[1, 2, 2]) - edge_to_vert[ordered_pair(subverts[0, 1, 2], subverts[2, 1, 2])] = subverts[1, 1, 2] - - # and finally the center. There are three ways to calculate, but elements should - # not intersect, so the edge_to_vert part isn't needed here. - subverts[1, 1, 1] = subdivide_edge(subverts[1, 1, 0], subverts[1, 1, 2]) - - - # Now make the hexahedron subelements - # order in which vtk expects vertices for a hexahedron - order = np.array([(0,0,0),(1,0,0),(1,1,0),(0,1,0),(0,0,1),(1,0,1),(1,1,1),(0,1,1)]) - for z in range(2): - for y in range(2): - for x in range(2): - hex_ = vtk.vtkHexahedron() - for vert_id in range(8): - coord = order[vert_id] + (x, y, z) - # minus one, since vtk starts at zero but marc starts at one - hex_.GetPointIds().SetId(vert_id, subverts[coord[0], coord[1], coord[2]] - 1) - grid.InsertNextCell(hex_.GetCellType(), hex_.GetPointIds()) - - else: - damask.util.croak('Unsupported Marc element type: {} (skipping)'.format(el['type'])) - -# Load all points -points = vtk.vtkPoints() -for point in range(1, len(coordinates) + 1): # marc indices start at 1 - points.InsertNextPoint(coordinates[point].tolist()) - -grid.SetPoints(points) - -# grid now contains the elements from the given marc file -writer = vtk.vtkXMLUnstructuredGridWriter() -writer.SetFileName(re.sub(r'\..+', ".vtu", args.filename)) # *.vtk extension does not work in paraview - -if vtk.VTK_MAJOR_VERSION <= 5: writer.SetInput(grid) -else: writer.SetInputData(grid) -writer.Write() - -if args.table is not None: - table.info_append([ - scriptID + ' ' + ' '.join(sys.argv[1:]), - ]) - table.head_write() - table.output_flush() - - table.data = np.array(nodal_data) - - table.data_writeArray() - -table.close() diff --git a/src/CPFEM2.f90 b/src/CPFEM2.f90 index bf8deeb78..38229d050 100644 --- a/src/CPFEM2.f90 +++ b/src/CPFEM2.f90 @@ -299,6 +299,8 @@ subroutine CPFEM_results(inc,time) pInt use results use HDF5_utilities + use homogenization, only: & + homogenization_results use constitutive, only: & constitutive_results use crystallite, only: & @@ -312,6 +314,7 @@ subroutine CPFEM_results(inc,time) call results_addIncrement(inc,time) call constitutive_results call crystallite_results + call homogenization_results call results_removeLink('current') ! ToDo: put this into closeJobFile call results_closeJobFile diff --git a/src/DAMASK_interface.f90 b/src/DAMASK_interface.f90 index ebe625ef6..b76813fe6 100644 --- a/src/DAMASK_interface.f90 +++ b/src/DAMASK_interface.f90 @@ -15,6 +15,9 @@ #define PETSC_MINOR_MIN 10 #define PETSC_MINOR_MAX 11 module DAMASK_interface + use prec + use system_routines + implicit none private logical, public, protected :: & @@ -39,6 +42,7 @@ module DAMASK_interface getLoadCaseFile, & rectifyPath, & makeRelativePath + contains !-------------------------------------------------------------------------------------------------- @@ -46,18 +50,8 @@ contains !! information on computation to screen !-------------------------------------------------------------------------------------------------- subroutine DAMASK_interface_init - use, intrinsic :: & - iso_fortran_env - use, intrinsic :: & - iso_c_binding + use, intrinsic :: iso_fortran_env use PETScSys - use prec, only: & - pReal - use system_routines, only: & - signalusr1_C, & - signalusr2_C, & - getHostName, & - getCWD #include #if defined(__GFORTRAN__) && __GNUC__ @brief basename of geometry file with extension from command line arguments !-------------------------------------------------------------------------------------------------- character(len=1024) function getGeometryFile(geometryParameter) - use system_routines, only: & - getCWD - implicit none character(len=1024), intent(in) :: geometryParameter logical :: file_exists external :: quit @@ -385,10 +366,7 @@ end function getGeometryFile !> @brief relative path of loadcase from command line arguments !-------------------------------------------------------------------------------------------------- character(len=1024) function getLoadCaseFile(loadCaseParameter) - use system_routines, only: & - getCWD - implicit none character(len=1024), intent(in) :: loadCaseParameter logical :: file_exists external :: quit @@ -412,7 +390,6 @@ end function getLoadCaseFile !-------------------------------------------------------------------------------------------------- function rectifyPath(path) - implicit none character(len=*) :: path character(len=1024) :: rectifyPath integer :: i,j,k,l @@ -457,7 +434,6 @@ end function rectifyPath !-------------------------------------------------------------------------------------------------- character(len=1024) function makeRelativePath(a,b) - implicit none character (len=*), intent(in) :: a,b character (len=1024) :: a_cleaned,b_cleaned integer :: i,posLastCommonSlash,remainingSlashes @@ -484,9 +460,7 @@ end function makeRelativePath !> @brief sets global variable SIGTERM to .true. !-------------------------------------------------------------------------------------------------- subroutine catchSIGTERM(signal) bind(C) - use :: iso_c_binding - implicit none integer(C_INT), value :: signal SIGTERM = .true. @@ -500,7 +474,6 @@ end subroutine catchSIGTERM !-------------------------------------------------------------------------------------------------- subroutine setSIGTERM(state) - implicit none logical, intent(in) :: state SIGTERM = state @@ -511,9 +484,7 @@ end subroutine setSIGTERM !> @brief sets global variable SIGUSR1 to .true. !-------------------------------------------------------------------------------------------------- subroutine catchSIGUSR1(signal) bind(C) - use :: iso_c_binding - implicit none integer(C_INT), value :: signal SIGUSR1 = .true. @@ -527,7 +498,6 @@ end subroutine catchSIGUSR1 !-------------------------------------------------------------------------------------------------- subroutine setSIGUSR1(state) - implicit none logical, intent(in) :: state SIGUSR1 = state @@ -538,9 +508,7 @@ end subroutine setSIGUSR1 !> @brief sets global variable SIGUSR2 to .true. if program receives SIGUSR2 !-------------------------------------------------------------------------------------------------- subroutine catchSIGUSR2(signal) bind(C) - use :: iso_c_binding - implicit none integer(C_INT), value :: signal SIGUSR2 = .true. @@ -554,7 +522,6 @@ end subroutine catchSIGUSR2 !-------------------------------------------------------------------------------------------------- subroutine setSIGUSR2(state) - implicit none logical, intent(in) :: state SIGUSR2 = state diff --git a/src/Lambert.f90 b/src/Lambert.f90 index 40bd04a01..c7b2c0d49 100644 --- a/src/Lambert.f90 +++ b/src/Lambert.f90 @@ -38,11 +38,7 @@ !> Modeling and Simulations in Materials Science and Engineering 22, 075013 (2014). !-------------------------------------------------------------------------- module Lambert - use prec, only: & - pReal - use math, only: & - PI - use future + use math implicit none private @@ -73,11 +69,7 @@ contains !> @brief map from 3D cubic grid to 3D ball !-------------------------------------------------------------------------- function LambertCubeToBall(cube) result(ball) - use, intrinsic :: IEEE_ARITHMETIC - use prec, only: & - dEq0 - implicit none real(pReal), intent(in), dimension(3) :: cube real(pReal), dimension(3) :: ball, LamXYZ, XYZ real(pReal), dimension(2) :: T @@ -133,15 +125,7 @@ end function LambertCubeToBall !> @brief map from 3D ball to 3D cubic grid !-------------------------------------------------------------------------- pure function LambertBallToCube(xyz) result(cube) - use, intrinsic :: IEEE_ARITHMETIC, only:& - IEEE_positive_inf, & - IEEE_value - use prec, only: & - dEq0 - use math, only: & - math_clip - implicit none real(pReal), intent(in), dimension(3) :: xyz real(pReal), dimension(3) :: cube, xyz1, xyz3 real(pReal), dimension(2) :: Tinv, xyz2 @@ -196,7 +180,6 @@ end function LambertBallToCube !-------------------------------------------------------------------------- pure function GetPyramidOrder(xyz) - implicit none real(pReal),intent(in),dimension(3) :: xyz integer, dimension(3) :: GetPyramidOrder diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 6428e19d1..4ac88797b 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -4,8 +4,7 @@ !> @brief elasticity, plasticity, internal microstructure state !-------------------------------------------------------------------------------------------------- module constitutive - use prec, only: & - pInt + use math implicit none private @@ -38,8 +37,6 @@ contains !> @brief allocates arrays pointing to array of the various constitutive modules !-------------------------------------------------------------------------------------------------- subroutine constitutive_init - use prec, only: & - pReal use debug, only: & debug_constitutive, & debug_levelBasic @@ -109,7 +106,6 @@ subroutine constitutive_init use kinematics_slipplane_opening use kinematics_thermal_expansion - implicit none integer, parameter :: FILEUNIT = 204 integer :: & o, & !< counter in output loop @@ -285,8 +281,6 @@ end subroutine constitutive_init !> ToDo: homogenizedC66 would be more consistent !-------------------------------------------------------------------------------------------------- function constitutive_homogenizedC(ipc,ip,el) - use prec, only: & - pReal use material, only: & phase_plasticity, & material_phase, & @@ -297,7 +291,6 @@ function constitutive_homogenizedC(ipc,ip,el) use lattice, only: & lattice_C66 - implicit none real(pReal), dimension(6,6) :: constitutive_homogenizedC integer, intent(in) :: & ipc, & !< component-ID of integration point @@ -317,8 +310,6 @@ end function constitutive_homogenizedC !> @brief calls microstructure function of the different constitutive models !-------------------------------------------------------------------------------------------------- subroutine constitutive_microstructure(Fe, Fp, ipc, ip, el) - use prec, only: & - pReal use material, only: & phasememberAt, & phase_plasticity, & @@ -337,7 +328,6 @@ subroutine constitutive_microstructure(Fe, Fp, ipc, ip, el) use plastic_disloUCLA, only: & plastic_disloUCLA_dependentState - implicit none integer, intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point @@ -376,8 +366,6 @@ end subroutine constitutive_microstructure !-------------------------------------------------------------------------------------------------- subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & S, Fi, ipc, ip, el) - use prec, only: & - pReal use material, only: & phasememberAt, & phase_plasticity, & @@ -408,7 +396,6 @@ subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & use plastic_nonlocal, only: & plastic_nonlocal_LpAndItsTangent - implicit none integer, intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point @@ -488,12 +475,6 @@ end subroutine constitutive_LpAndItsTangents !-------------------------------------------------------------------------------------------------- subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & S, Fi, ipc, ip, el) - use prec, only: & - pReal - use math, only: & - math_I3, & - math_inv33, & - math_det33 use material, only: & phasememberAt, & phase_plasticity, & @@ -515,7 +496,6 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & use kinematics_thermal_expansion, only: & kinematics_thermal_expansion_LiAndItsTangent - implicit none integer, intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point @@ -593,10 +573,6 @@ end subroutine constitutive_LiAndItsTangents !> @brief collects initial intermediate deformation gradient !-------------------------------------------------------------------------------------------------- pure function constitutive_initialFi(ipc, ip, el) - use prec, only: & - pReal - use math, only: & - math_I3 use material, only: & material_phase, & material_homogenizationAt, & @@ -608,7 +584,6 @@ pure function constitutive_initialFi(ipc, ip, el) use kinematics_thermal_expansion, only: & kinematics_thermal_expansion_initialStrain - implicit none integer, intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point @@ -643,10 +618,7 @@ end function constitutive_initialFi !! (so far no case switch because only Hooke is implemented) !-------------------------------------------------------------------------------------------------- subroutine constitutive_SandItsTangents(S, dS_dFe, dS_dFi, Fe, Fi, ipc, ip, el) - use prec, only: & - pReal - implicit none integer, intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point @@ -672,12 +644,6 @@ end subroutine constitutive_SandItsTangents !-------------------------------------------------------------------------------------------------- subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, & Fe, Fi, ipc, ip, el) - use prec, only: & - pReal - use math, only : & - math_mul3333xx33, & - math_66toSym3333, & - math_I3 use material, only: & material_phase, & material_homogenizationAt, & @@ -687,7 +653,6 @@ subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, & damageMapping, & STIFFNESS_DEGRADATION_damage_ID - implicit none integer, intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point @@ -735,8 +700,6 @@ end subroutine constitutive_hooke_SandItsTangents !> @brief contains the constitutive equation for calculating the rate of change of microstructure !-------------------------------------------------------------------------------------------------- subroutine constitutive_collectDotState(S, FeArray, Fi, FpArray, subdt, ipc, ip, el) - use prec, only: & - pReal use debug, only: & debug_level, & debug_constitutive, & @@ -786,7 +749,6 @@ subroutine constitutive_collectDotState(S, FeArray, Fi, FpArray, subdt, ipc, ip, use source_thermal_externalheat, only: & source_thermal_externalheat_dotState - implicit none integer, intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point @@ -873,8 +835,6 @@ end subroutine constitutive_collectDotState !> will return false if delta state is not needed/supported by the constitutive model !-------------------------------------------------------------------------------------------------- subroutine constitutive_collectDeltaState(S, Fe, Fi, ipc, ip, el) - use prec, only: & - pReal use debug, only: & debug_level, & debug_constitutive, & @@ -896,7 +856,6 @@ subroutine constitutive_collectDeltaState(S, Fe, Fi, ipc, ip, el) use source_damage_isoBrittle, only: & source_damage_isoBrittle_deltaState - implicit none integer, intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point @@ -944,8 +903,6 @@ end subroutine constitutive_collectDeltaState !> @brief returns array of constitutive results !-------------------------------------------------------------------------------------------------- function constitutive_postResults(S, Fi, ipc, ip, el) - use prec, only: & - pReal use material, only: & phasememberAt, & phase_plasticityInstance, & @@ -990,7 +947,6 @@ function constitutive_postResults(S, Fi, ipc, ip, el) use source_damage_anisoDuctile, only: & source_damage_anisoDuctile_postResults - implicit none integer, intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point diff --git a/src/crystallite.f90 b/src/crystallite.f90 index e972ff162..07a8d068a 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -937,8 +937,6 @@ end function crystallite_push33ToRef !-------------------------------------------------------------------------------------------------- function crystallite_postResults(ipc, ip, el) use math, only: & - math_qToEuler, & - math_qToEulerAxisAngle, & math_det33, & math_I3, & inDeg diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 06da6ab2e..fe15364b8 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -7,6 +7,7 @@ module homogenization use prec, only: & pReal + use material !-------------------------------------------------------------------------------------------------- ! General variables for the homogenization at a material point @@ -66,7 +67,8 @@ module homogenization public :: & homogenization_init, & materialpoint_stressAndItsTangent, & - materialpoint_postResults + materialpoint_postResults, & + homogenization_results private :: & partitionDeformation, & updateState, & @@ -100,7 +102,6 @@ subroutine homogenization_init config_deallocate, & config_homogenization, & homogenization_name - use material use homogenization_mech_RGC use thermal_isothermal use thermal_adiabatic @@ -291,7 +292,6 @@ subroutine homogenization_init write(6,'(a32,1x,7(i8,1x))') 'materialpoint_requested: ', shape(materialpoint_requested) write(6,'(a32,1x,7(i8,1x))') 'materialpoint_converged: ', shape(materialpoint_converged) write(6,'(a32,1x,7(i8,1x),/)') 'materialpoint_doneAndHappy: ', shape(materialpoint_doneAndHappy) - write(6,'(a32,1x,7(i8,1x))') 'maxSizePostResults: ', homogenization_maxSizePostResults endif flush(6) @@ -316,17 +316,6 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt) terminallyIll use mesh, only: & mesh_element - use material, only: & - plasticState, & - sourceState, & - homogState, & - thermalState, & - damageState, & - phase_Nsources, & - material_homogenizationAt, & - mappingHomogenization, & - phaseAt, phasememberAt, & - homogenization_Ngrains use crystallite, only: & crystallite_F0, & crystallite_Fp0, & @@ -658,16 +647,6 @@ subroutine materialpoint_postResults FEsolving_execIP use mesh, only: & mesh_element - use material, only: & - material_homogenizationAt, & - homogState, & - thermalState, & - damageState, & - plasticState, & - sourceState, & - material_phase, & - homogenization_Ngrains, & - microstructure_crystallite use crystallite, only: & crystallite_sizePostResults, & crystallite_postResults @@ -723,12 +702,6 @@ end subroutine materialpoint_postResults subroutine partitionDeformation(ip,el) use mesh, only: & mesh_element - use material, only: & - homogenization_type, & - homogenization_Ngrains, & - HOMOGENIZATION_NONE_ID, & - HOMOGENIZATION_ISOSTRAIN_ID, & - HOMOGENIZATION_RGC_ID use crystallite, only: & crystallite_partionedF use homogenization_mech_RGC, only: & @@ -767,14 +740,6 @@ end subroutine partitionDeformation function updateState(ip,el) use mesh, only: & mesh_element - use material, only: & - homogenization_type, & - thermal_type, & - damage_type, & - homogenization_Ngrains, & - HOMOGENIZATION_RGC_ID, & - THERMAL_adiabatic_ID, & - DAMAGE_local_ID use crystallite, only: & crystallite_P, & crystallite_dPdF, & @@ -835,13 +800,6 @@ end function updateState subroutine averageStressAndItsTangent(ip,el) use mesh, only: & mesh_element - use material, only: & - homogenization_type, & - homogenization_typeInstance, & - homogenization_Ngrains, & - HOMOGENIZATION_NONE_ID, & - HOMOGENIZATION_ISOSTRAIN_ID, & - HOMOGENIZATION_RGC_ID use crystallite, only: & crystallite_P,crystallite_dPdF use homogenization_mech_RGC, only: & @@ -884,27 +842,6 @@ end subroutine averageStressAndItsTangent function postResults(ip,el) use mesh, only: & mesh_element - use material, only: & - thermalMapping, & - thermal_typeInstance, & - material_homogenizationAt, & - homogenization_typeInstance,& - mappingHomogenization, & - homogState, & - thermalState, & - damageState, & - homogenization_type, & - thermal_type, & - damage_type, & - HOMOGENIZATION_NONE_ID, & - HOMOGENIZATION_ISOSTRAIN_ID, & - HOMOGENIZATION_RGC_ID, & - THERMAL_isothermal_ID, & - THERMAL_adiabatic_ID, & - THERMAL_conduction_ID, & - DAMAGE_none_ID, & - DAMAGE_local_ID, & - DAMAGE_nonlocal_ID use homogenization_mech_RGC, only: & homogenization_RGC_postResults use thermal_adiabatic, only: & @@ -969,4 +906,41 @@ function postResults(ip,el) end function postResults + +!-------------------------------------------------------------------------------------------------- +!> @brief writes homogenization results to HDF5 output file +!-------------------------------------------------------------------------------------------------- +subroutine homogenization_results +#if defined(PETSc) || defined(DAMASK_HDF5) + use results + use homogenization_mech_RGC + use HDF5_utilities + use config, only: & + config_name_homogenization => homogenization_name ! anticipate logical name + + use material, only: & + homogenization_typeInstance, & + material_homogenization_type => homogenization_type + + integer :: p + character(len=256) :: group + + do p=1,size(config_name_homogenization) + group = trim('current/materialpoint')//'/'//trim(config_name_homogenization(p)) + call HDF5_closeGroup(results_addGroup(group)) + + group = trim(group)//'/mech' + + call HDF5_closeGroup(results_addGroup(group)) + select case(material_homogenization_type(p)) + case(HOMOGENIZATION_rgc_ID) + call mech_RGC_results(homogenization_typeInstance(p),group) + end select + + enddo +#endif + + +end subroutine homogenization_results + end module homogenization diff --git a/src/homogenization_mech_RGC.f90 b/src/homogenization_mech_RGC.f90 index 012ad6086..27a52432e 100644 --- a/src/homogenization_mech_RGC.f90 +++ b/src/homogenization_mech_RGC.f90 @@ -9,6 +9,7 @@ module homogenization_mech_RGC use prec, only: & pReal + use material implicit none private @@ -75,7 +76,8 @@ module homogenization_mech_RGC homogenization_RGC_partitionDeformation, & homogenization_RGC_averageStressAndItsTangent, & homogenization_RGC_updateState, & - homogenization_RGC_postResults + homogenization_RGC_postResults, & + mech_RGC_results ! name suited for planned submodule situation private :: & relaxationVector, & interfaceNormal, & @@ -104,18 +106,6 @@ subroutine homogenization_RGC_init() INRAD use IO, only: & IO_error - use material, only: & -#ifdef DEBUG - mappingHomogenization, & -#endif - homogenization_type, & - material_homogenizationAt, & - homogState, & - HOMOGENIZATION_RGC_ID, & - HOMOGENIZATION_RGC_LABEL, & - homogenization_typeInstance, & - homogenization_Noutput, & - homogenization_Ngrains use config, only: & config_homogenization @@ -322,10 +312,6 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) #endif use math, only: & math_invert2 - use material, only: & - material_homogenizationAt, & - homogenization_typeInstance, & - mappingHomogenization use numerics, only: & absTol_RGC, & relTol_RGC, & @@ -1109,6 +1095,54 @@ pure function homogenization_RGC_postResults(instance,of) result(postResults) end function homogenization_RGC_postResults +!-------------------------------------------------------------------------------------------------- +!> @brief writes results to HDF5 output file +! ToDo: check wheter units are correct +!-------------------------------------------------------------------------------------------------- +subroutine mech_RGC_results(instance,group) +#if defined(PETSc) || defined(DAMASK_HDF5) + use results, only: & + results_writeDataset + + integer, intent(in) :: instance + character(len=*) :: group + integer :: o + + associate(stt => state(instance), dst => dependentState(instance), prm => param(instance)) + + outputsLoop: do o = 1,size(prm%outputID) + select case(prm%outputID(o)) + + case (constitutivework_ID) + call results_writeDataset(group,stt%work,'W',& + 'work density','J/m³') + case (magnitudemismatch_ID) + call results_writeDataset(group,dst%mismatch,'N',& + 'average mismatch tensor','1') + case (penaltyenergy_ID) + call results_writeDataset(group,stt%penaltyEnergy,'R',& + 'mismatch penalty density','J/m³') + case (volumediscrepancy_ID) + call results_writeDataset(group,dst%volumeDiscrepancy,'Delta_V',& + 'volume discrepancy','m³') + case (maximumrelaxrate_ID) + call results_writeDataset(group,dst%relaxationrate_max,'max_alpha_dot',& + 'maximum relaxation rate','m/s') + case (averagerelaxrate_ID) + call results_writeDataset(group,dst%relaxationrate_avg,'avg_alpha_dot',& + 'average relaxation rate','m/s') + end select + enddo outputsLoop + end associate + +#else + integer, intent(in) :: instance + character(len=*) :: group +#endif + +end subroutine mech_RGC_results + + !-------------------------------------------------------------------------------------------------- !> @brief collect relaxation vectors of an interface !-------------------------------------------------------------------------------------------------- diff --git a/src/homogenization_mech_isostrain.f90 b/src/homogenization_mech_isostrain.f90 index 1fdf5435c..39ed3a8c6 100644 --- a/src/homogenization_mech_isostrain.f90 +++ b/src/homogenization_mech_isostrain.f90 @@ -36,13 +36,6 @@ module subroutine mech_isostrain_init debug_levelBasic use IO, only: & IO_error - use material, only: & - homogenization_type, & - material_homogenizationAt, & - homogState, & - HOMOGENIZATION_ISOSTRAIN_ID, & - HOMOGENIZATION_ISOSTRAIN_LABEL, & - homogenization_typeInstance use config, only: & config_homogenization diff --git a/src/homogenization_mech_none.f90 b/src/homogenization_mech_none.f90 index 4ac509363..8300acce6 100644 --- a/src/homogenization_mech_none.f90 +++ b/src/homogenization_mech_none.f90 @@ -20,13 +20,7 @@ module subroutine mech_none_init debug_levelBasic use config, only: & config_homogenization - use material, only: & - homogenization_type, & - material_homogenizationAt, & - homogState, & - HOMOGENIZATION_NONE_LABEL, & - HOMOGENIZATION_NONE_ID - + implicit none integer :: & Ninstance, & diff --git a/src/list.f90 b/src/list.f90 index 1c4c4a8a0..be80b151d 100644 --- a/src/list.f90 +++ b/src/list.f90 @@ -70,7 +70,6 @@ subroutine add(this,string) IO_lc, & IO_stringPos - implicit none class(tPartitionedStringList), target, intent(in) :: this character(len=*), intent(in) :: string type(tPartitionedStringList), pointer :: new, temp @@ -95,7 +94,6 @@ end subroutine add !-------------------------------------------------------------------------------------------------- subroutine show(this) - implicit none class(tPartitionedStringList), target, intent(in) :: this type(tPartitionedStringList), pointer :: item @@ -114,7 +112,6 @@ end subroutine show !-------------------------------------------------------------------------------------------------- subroutine free(this) - implicit none class(tPartitionedStringList), intent(inout) :: this if(associated(this%next)) deallocate(this%next) @@ -128,7 +125,6 @@ end subroutine free !-------------------------------------------------------------------------------------------------- recursive subroutine finalize(this) - implicit none type(tPartitionedStringList), intent(inout) :: this if(associated(this%next)) deallocate(this%next) @@ -142,7 +138,6 @@ end subroutine finalize !-------------------------------------------------------------------------------------------------- subroutine finalizeArray(this) - implicit none integer :: i type(tPartitionedStringList), intent(inout), dimension(:) :: this type(tPartitionedStringList), pointer :: temp ! bug in Gfortran? @@ -165,7 +160,6 @@ logical function keyExists(this,key) use IO, only: & IO_stringValue - implicit none class(tPartitionedStringList), target, intent(in) :: this character(len=*), intent(in) :: key type(tPartitionedStringList), pointer :: item @@ -189,8 +183,6 @@ integer function countKeys(this,key) use IO, only: & IO_stringValue - implicit none - class(tPartitionedStringList), target, intent(in) :: this character(len=*), intent(in) :: key type(tPartitionedStringList), pointer :: item @@ -218,7 +210,6 @@ real(pReal) function getFloat(this,key,defaultVal) IO_stringValue, & IO_FloatValue - implicit none class(tPartitionedStringList), target, intent(in) :: this character(len=*), intent(in) :: key real(pReal), intent(in), optional :: defaultVal @@ -255,7 +246,6 @@ integer function getInt(this,key,defaultVal) IO_stringValue, & IO_IntValue - implicit none class(tPartitionedStringList), target, intent(in) :: this character(len=*), intent(in) :: key integer, intent(in), optional :: defaultVal @@ -292,7 +282,6 @@ character(len=65536) function getString(this,key,defaultVal,raw) IO_error, & IO_stringValue - implicit none class(tPartitionedStringList), target, intent(in) :: this character(len=*), intent(in) :: key character(len=*), intent(in), optional :: defaultVal @@ -343,7 +332,6 @@ function getFloats(this,key,defaultVal,requiredSize) IO_stringValue, & IO_FloatValue - implicit none real(pReal), dimension(:), allocatable :: getFloats class(tPartitionedStringList), target, intent(in) :: this character(len=*), intent(in) :: key @@ -393,7 +381,6 @@ function getInts(this,key,defaultVal,requiredSize) IO_stringValue, & IO_IntValue - implicit none integer, dimension(:), allocatable :: getInts class(tPartitionedStringList), target, intent(in) :: this character(len=*), intent(in) :: key @@ -443,7 +430,6 @@ function getStrings(this,key,defaultVal,raw) IO_error, & IO_StringValue - implicit none character(len=65536),dimension(:), allocatable :: getStrings class(tPartitionedStringList),target, intent(in) :: this character(len=*), intent(in) :: key diff --git a/src/material.f90 b/src/material.f90 index 64c3c2529..bb8d6dbff 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -8,15 +8,9 @@ !! 'phase', 'texture', and 'microstucture' !-------------------------------------------------------------------------------------------------- module material - use prec, only: & - pReal, & - pInt, & - tState, & - tPlasticState, & - tSourceState, & - tHomogMapping, & - group_float, & - group_int + use prec + use math + use config implicit none private @@ -53,55 +47,35 @@ module material enum, bind(c) enumerator :: ELASTICITY_undefined_ID, & - ELASTICITY_hooke_ID - end enum - enum, bind(c) - enumerator :: PLASTICITY_undefined_ID, & + ELASTICITY_hooke_ID, & + PLASTICITY_undefined_ID, & PLASTICITY_none_ID, & PLASTICITY_isotropic_ID, & PLASTICITY_phenopowerlaw_ID, & PLASTICITY_kinehardening_ID, & PLASTICITY_dislotwin_ID, & PLASTICITY_disloucla_ID, & - PLASTICITY_nonlocal_ID - end enum - - enum, bind(c) - enumerator :: SOURCE_undefined_ID, & + PLASTICITY_nonlocal_ID, & + SOURCE_undefined_ID, & SOURCE_thermal_dissipation_ID, & SOURCE_thermal_externalheat_ID, & SOURCE_damage_isoBrittle_ID, & SOURCE_damage_isoDuctile_ID, & SOURCE_damage_anisoBrittle_ID, & - SOURCE_damage_anisoDuctile_ID - end enum - - enum, bind(c) - enumerator :: KINEMATICS_undefined_ID, & + SOURCE_damage_anisoDuctile_ID, & + KINEMATICS_undefined_ID, & KINEMATICS_cleavage_opening_ID, & KINEMATICS_slipplane_opening_ID, & - KINEMATICS_thermal_expansion_ID - end enum - - enum, bind(c) - enumerator :: STIFFNESS_DEGRADATION_undefined_ID, & - STIFFNESS_DEGRADATION_damage_ID - end enum - - enum, bind(c) - enumerator :: THERMAL_isothermal_ID, & + KINEMATICS_thermal_expansion_ID, & + STIFFNESS_DEGRADATION_undefined_ID, & + STIFFNESS_DEGRADATION_damage_ID, & + THERMAL_isothermal_ID, & THERMAL_adiabatic_ID, & - THERMAL_conduction_ID - end enum - - enum, bind(c) - enumerator :: DAMAGE_none_ID, & + THERMAL_conduction_ID, & + DAMAGE_none_ID, & DAMAGE_local_ID, & - DAMAGE_nonlocal_ID - end enum - - enum, bind(c) - enumerator :: HOMOGENIZATION_undefined_ID, & + DAMAGE_nonlocal_ID, & + HOMOGENIZATION_undefined_ID, & HOMOGENIZATION_none_ID, & HOMOGENIZATION_isostrain_ID, & HOMOGENIZATION_rgc_ID @@ -279,20 +253,9 @@ subroutine material_init debug_material, & debug_levelBasic, & debug_levelExtensive - use config, only: & - config_crystallite, & - config_homogenization, & - config_microstructure, & - config_phase, & - config_texture, & - homogenization_name, & - microstructure_name, & - phase_name, & - texture_name use mesh, only: & theMesh - implicit none integer(pInt), parameter :: FILEUNIT = 210_pInt integer(pInt) :: m,c,h, myDebug, myPhase, myHomog integer(pInt) :: & @@ -461,14 +424,11 @@ end subroutine material_init !> @brief parses the homogenization part from the material configuration !-------------------------------------------------------------------------------------------------- subroutine material_parseHomogenization - use config, only : & - config_homogenization use mesh, only: & theMesh use IO, only: & IO_error - implicit none integer(pInt) :: h character(len=65536) :: tag @@ -559,21 +519,15 @@ end subroutine material_parseHomogenization !> @brief parses the microstructure part in the material configuration file !-------------------------------------------------------------------------------------------------- subroutine material_parseMicrostructure - use prec, only: & - dNeq use IO, only: & IO_floatValue, & IO_intValue, & IO_stringValue, & IO_stringPos, & IO_error - use config, only: & - config_microstructure, & - microstructure_name use mesh, only: & theMesh - implicit none character(len=65536), dimension(:), allocatable :: & strings integer(pInt), allocatable, dimension(:) :: chunkPos @@ -637,10 +591,7 @@ end subroutine material_parseMicrostructure !> @brief parses the crystallite part in the material configuration file !-------------------------------------------------------------------------------------------------- subroutine material_parseCrystallite - use config, only: & - config_crystallite - implicit none integer(pInt) :: c allocate(crystallite_Noutput(size(config_crystallite)),source=0_pInt) @@ -659,10 +610,7 @@ subroutine material_parsePhase IO_error, & IO_getTag, & IO_stringValue - use config, only: & - config_phase - implicit none integer(pInt) :: sourceCtr, kinematicsCtr, stiffDegradationCtr, p character(len=65536), dimension(:), allocatable :: str @@ -785,23 +733,12 @@ end subroutine material_parsePhase !> @brief parses the texture part in the material configuration file !-------------------------------------------------------------------------------------------------- subroutine material_parseTexture - use prec, only: & - dNeq use IO, only: & IO_error, & IO_stringPos, & IO_floatValue, & IO_stringValue - use config, only: & - config_deallocate, & - config_texture - use math, only: & - inRad, & - math_sampleRandomOri, & - math_I3, & - math_det33 - implicit none integer(pInt) :: section, gauss, j, t, i character(len=65536), dimension(:), allocatable :: strings ! Values for given key in material config integer(pInt), dimension(:), allocatable :: chunkPos @@ -880,7 +817,6 @@ subroutine material_allocatePlasticState(phase,NofMyPhase,& use numerics, only: & numerics_integrator - implicit none integer(pInt), intent(in) :: & phase, & NofMyPhase, & @@ -928,7 +864,6 @@ subroutine material_allocateSourceState(phase,of,NofMyPhase,& use numerics, only: & numerics_integrator - implicit none integer(pInt), intent(in) :: & phase, & of, & @@ -967,47 +902,15 @@ end subroutine material_allocateSourceState !! calculates the volume of the grains and deals with texture components !-------------------------------------------------------------------------------------------------- subroutine material_populateGrains - use math, only: & - math_EulertoR, & - math_RtoEuler, & - math_mul33x33, & - math_range use mesh, only: & theMesh - use config, only: & - config_homogenization, & - config_microstructure, & - config_deallocate - use IO, only: & - IO_error - - implicit none - integer(pInt), dimension (:,:), allocatable :: Ngrains - integer(pInt), dimension (microstructure_maxNconstituents) :: & - NgrainsOfConstituent, & - currentGrainOfConstituent, & - randomOrder - real(pReal), dimension (microstructure_maxNconstituents) :: & - rndArray - real(pReal), dimension (:,:), allocatable :: orientationOfGrain - real(pReal), dimension (3) :: orientation - integer(pInt), dimension (:), allocatable :: phaseOfGrain, textureOfGrain - integer(pInt) :: t,e,i,g,j,m,c,r,homog,micro,sgn,hme, & - phaseID,textureID,dGrains,myNgrains,myNorientations,myNconstituents, & - grain,constituentGrain,ipGrain,ip - real(pReal) :: deviation,extreme,rnd - integer(pInt), dimension (:,:), allocatable :: Nelems ! counts number of elements in homog, micro array - type(group_int), dimension (:,:), allocatable :: elemsOfHomogMicro ! lists element number in homog, micro array + integer(pInt) :: e,i,c,homog,micro allocate(material_phase(homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%Nelems), source=0_pInt) allocate(material_texture(homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%Nelems), source=0_pInt) allocate(material_EulerAngles(3,homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%Nelems),source=0.0_pReal) - allocate(Ngrains(size(config_homogenization),size(config_microstructure)), source=0_pInt) - allocate(Nelems (size(config_homogenization),size(config_microstructure)), source=0_pInt) - - do e = 1, theMesh%Nelems do i = 1, theMesh%elem%nIPs homog = theMesh%homogenizationAt(e) @@ -1017,7 +920,7 @@ subroutine material_populateGrains material_texture(c,i,e) = microstructure_texture(c,micro) material_EulerAngles(1:3,c,i,e) = texture_Gauss(1:3,1,material_texture(c,i,e)) material_EulerAngles(1:3,c,i,e) = math_RtoEuler( & ! translate back to Euler angles - math_mul33x33( & ! pre-multiply + matmul( & ! pre-multiply math_EulertoR(material_EulerAngles(1:3,c,i,e)), & ! face-value orientation texture_transformation(1:3,1:3,material_texture(c,i,e)) & ! and transformation matrix ) & @@ -1026,209 +929,6 @@ subroutine material_populateGrains enddo enddo - return -!-------------------------------------------------------------------------------------------------- -! precounting of elements for each homog/micro pair - do e = 1_pInt, theMesh%Nelems - homog = theMesh%homogenizationAt(e) - micro = theMesh%microstructureAt(e) - Nelems(homog,micro) = Nelems(homog,micro) + 1_pInt - enddo - allocate(elemsOfHomogMicro(size(config_homogenization),size(config_microstructure))) - do homog = 1,size(config_homogenization) - do micro = 1,size(config_microstructure) - if (Nelems(homog,micro) > 0_pInt) then - allocate(elemsOfHomogMicro(homog,micro)%p(Nelems(homog,micro))) - elemsOfHomogMicro(homog,micro)%p = 0_pInt - endif - enddo - enddo - -!-------------------------------------------------------------------------------------------------- -! identify maximum grain count per IP (from element) and find grains per homog/micro pair - Nelems = 0_pInt ! reuse as counter - elementLooping: do e = 1_pInt,theMesh%Nelems - homog = theMesh%homogenizationAt(e) - micro = theMesh%microstructureAt(e) - if (homog < 1_pInt .or. homog > size(config_homogenization)) & ! out of bounds - call IO_error(154_pInt,e,0_pInt,0_pInt) - if (micro < 1_pInt .or. micro > size(config_microstructure)) & ! out of bounds - call IO_error(155_pInt,e,0_pInt,0_pInt) - if (microstructure_elemhomo(micro)) then ! how many grains are needed at this element? - dGrains = homogenization_Ngrains(homog) ! only one set of Ngrains (other IPs are plain copies) - else - dGrains = homogenization_Ngrains(homog) * theMesh%elem%nIPs ! each IP has Ngrains - endif - Ngrains(homog,micro) = Ngrains(homog,micro) + dGrains ! total grain count - Nelems(homog,micro) = Nelems(homog,micro) + 1_pInt ! total element count - elemsOfHomogMicro(homog,micro)%p(Nelems(homog,micro)) = e ! remember elements active in this homog/micro pair - enddo elementLooping - - allocate(phaseOfGrain(maxval(Ngrains)), source=0_pInt) ! reserve memory for maximum case - allocate(textureOfGrain(maxval(Ngrains)), source=0_pInt) ! reserve memory for maximum case - allocate(orientationOfGrain(3,maxval(Ngrains)),source=0.0_pReal) ! reserve memory for maximum case - - - homogenizationLoop: do homog = 1_pInt,size(config_homogenization) - dGrains = homogenization_Ngrains(homog) ! grain number per material point - microstructureLoop: do micro = 1_pInt,size(config_microstructure) ! all pairs of homog and micro - activePair: if (Ngrains(homog,micro) > 0_pInt) then - myNgrains = Ngrains(homog,micro) ! assign short name for total number of grains to populate - myNconstituents = microstructure_Nconstituents(micro) ! assign short name for number of constituents - -!-------------------------------------------------------------------------------------------------- -! divide myNgrains as best over constituents -! -! example: three constituents with fractions of 0.25, 0.25, and 0.5 distributed over 20 (microstructure) grains -! -! ***** ***** ********** -! NgrainsOfConstituent: 5, 5, 10 -! counters: -! |-----> grain (if constituent == 2) -! |--> constituentGrain (of constituent 2) -! - - NgrainsOfConstituent = 0_pInt ! reset counter of grains per constituent - forall (i = 1_pInt:myNconstituents) & - NgrainsOfConstituent(i) = nint(microstructure_fraction(i,micro)*real(myNgrains,pReal),pInt)! do rounding integer conversion - do while (sum(NgrainsOfConstituent) /= myNgrains) ! total grain count over constituents wrong? - sgn = sign(1_pInt, myNgrains - sum(NgrainsOfConstituent)) ! direction of required change - extreme = 0.0_pReal - t = 0_pInt - do i = 1_pInt,myNconstituents ! find largest deviator - deviation = real(sgn,pReal)*log( microstructure_fraction(i,micro) / & - !-------------------------------- & - (real(NgrainsOfConstituent(i),pReal)/real(myNgrains,pReal) ) ) - if (deviation > extreme) then - extreme = deviation - t = i - endif - enddo - NgrainsOfConstituent(t) = NgrainsOfConstituent(t) + sgn ! change that by one - enddo - -!-------------------------------------------------------------------------------------------------- -! assign phase and texture info - - phaseOfGrain = 0_pInt - textureOfGrain = 0_pInt - orientationOfGrain = 0.0_pReal - - texture: do i = 1_pInt,myNconstituents ! loop over constituents - grain = sum(NgrainsOfConstituent(1_pInt:i-1_pInt)) ! set microstructure grain index of current constituent - ! "grain" points to start of this constituent's grain population - constituentGrain = 0_pInt ! constituent grain index - - phaseID = microstructure_phase(i,micro) - textureID = microstructure_texture(i,micro) - phaseOfGrain (grain+1_pInt:grain+NgrainsOfConstituent(i)) = phaseID ! assign resp. phase - textureOfGrain(grain+1_pInt:grain+NgrainsOfConstituent(i)) = textureID ! assign resp. texture - - myNorientations = ceiling(real(NgrainsOfConstituent(i),pReal)/1.0,pInt) ! max number of unique orientations (excl. symmetry) - -!-------------------------------------------------------------------------------------------------- -! has texture components - gauss: do t = 1_pInt,texture_Ngauss(textureID) ! loop over Gauss components - do g = 1_pInt,int(real(myNorientations,pReal)*texture_Gauss(5,t,textureID),pInt) ! loop over required grain count - orientationOfGrain(:,grain+constituentGrain+g) = texture_Gauss(1:3,t,textureID) - enddo - constituentGrain = & - constituentGrain + int(real(myNorientations,pReal)*texture_Gauss(5,t,textureID)) ! advance counter for grains of current constituent - enddo gauss - - -!-------------------------------------------------------------------------------------------------- -! ...texture transformation - - do j = 1_pInt,myNorientations ! loop over each "real" orientation - orientationOfGrain(1:3,grain+j) = math_RtoEuler( & ! translate back to Euler angles - math_mul33x33( & ! pre-multiply - math_EulertoR(orientationOfGrain(1:3,grain+j)), & ! face-value orientation - texture_transformation(1:3,1:3,textureID) & ! and transformation matrix - ) & - ) - enddo - -!-------------------------------------------------------------------------------------------------- -! shuffle grains within current constituent - - do j = 1_pInt,NgrainsOfConstituent(i)-1_pInt ! walk thru grains of current constituent - call random_number(rnd) - t = nint(rnd*real(NgrainsOfConstituent(i)-j,pReal)+real(j,pReal)+0.5_pReal,pInt) ! select a grain in remaining list - m = phaseOfGrain(grain+t) ! exchange current with random - phaseOfGrain(grain+t) = phaseOfGrain(grain+j) - phaseOfGrain(grain+j) = m - m = textureOfGrain(grain+t) ! exchange current with random - textureOfGrain(grain+t) = textureOfGrain(grain+j) - textureOfGrain(grain+j) = m - orientation = orientationOfGrain(1:3,grain+t) ! exchange current with random - orientationOfGrain(1:3,grain+t) = orientationOfGrain(1:3,grain+j) - orientationOfGrain(1:3,grain+j) = orientation - enddo - - enddo texture -!< @todo calc fraction after weighing with volumePerGrain, exchange in MC steps to improve result (humbug at the moment) - - - -!-------------------------------------------------------------------------------------------------- -! distribute grains of all constituents as accurately as possible to given constituent fractions - - ip = 0_pInt - currentGrainOfConstituent = 0_pInt - - do hme = 1_pInt, Nelems(homog,micro) - e = elemsOfHomogMicro(homog,micro)%p(hme) ! only perform calculations for elements with homog, micro combinations which is indexed in cpElemsindex - if (microstructure_elemhomo(micro)) then ! homogeneous distribution of grains over each element's IPs - m = 1_pInt ! process only first IP - else - m = theMesh%elem%nIPs - endif - - do i = 1_pInt, m ! loop over necessary IPs - ip = ip + 1_pInt ! keep track of total ip count - ipGrain = 0_pInt ! count number of grains assigned at this IP - randomOrder = math_range(microstructure_maxNconstituents) ! start out with ordered sequence of constituents - call random_number(rndArray) ! as many rnd numbers as (max) constituents - do j = 1_pInt, myNconstituents - 1_pInt ! loop over constituents ... - r = nint(rndArray(j)*real(myNconstituents-j,pReal)+real(j,pReal)+0.5_pReal,pInt) ! ... select one in remaining list - c = randomOrder(r) ! ... call it "c" - randomOrder(r) = randomOrder(j) ! ... and exchange with present position in constituent list - grain = sum(NgrainsOfConstituent(1:c-1_pInt)) ! figure out actual starting index in overall/consecutive grain population - do g = 1_pInt, min(dGrains-ipGrain, & ! leftover number of grains at this IP - max(0_pInt, & ! no negative values - nint(real(ip * dGrains * NgrainsOfConstituent(c)) / & ! fraction of grains scaled to this constituent... - real(myNgrains),pInt) - & ! ...minus those already distributed - currentGrainOfConstituent(c))) - ipGrain = ipGrain + 1_pInt ! advance IP grain counter - currentGrainOfConstituent(c) = currentGrainOfConstituent(c) + 1_pInt ! advance index of grain population for constituent c - material_phase(ipGrain,i,e) = phaseOfGrain(grain+currentGrainOfConstituent(c)) - material_texture(ipGrain,i,e) = textureOfGrain(grain+currentGrainOfConstituent(c)) - material_EulerAngles(1:3,ipGrain,i,e) = orientationOfGrain(1:3,grain+currentGrainOfConstituent(c)) - enddo; enddo - - c = randomOrder(microstructure_Nconstituents(micro)) ! look up constituent remaining after random shuffling - grain = sum(NgrainsOfConstituent(1:c-1_pInt)) ! figure out actual starting index in overall/consecutive grain population - do ipGrain = ipGrain + 1_pInt, dGrains ! ensure last constituent fills up to dGrains - currentGrainOfConstituent(c) = currentGrainOfConstituent(c) + 1_pInt - material_phase(ipGrain,i,e) = phaseOfGrain(grain+currentGrainOfConstituent(c)) - material_texture(ipGrain,i,e) = textureOfGrain(grain+currentGrainOfConstituent(c)) - material_EulerAngles(1:3,ipGrain,i,e) = orientationOfGrain(1:3,grain+currentGrainOfConstituent(c)) - enddo - - enddo - - do i = i, theMesh%elem%nIPs ! loop over IPs to (possibly) distribute copies from first IP - material_phase (1_pInt:dGrains,i,e) = material_phase (1_pInt:dGrains,1,e) - material_texture(1_pInt:dGrains,i,e) = material_texture(1_pInt:dGrains,1,e) - material_EulerAngles(1:3,1_pInt:dGrains,i,e) = material_EulerAngles(1:3,1_pInt:dGrains,1,e) - enddo - - enddo - endif activePair - enddo microstructureLoop - enddo homogenizationLoop - deallocate(texture_transformation) call config_deallocate('material.config/microstructure') diff --git a/src/math.f90 b/src/math.f90 index 3338b99e3..1740ebdb7 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -6,212 +6,117 @@ !> @brief Mathematical library, including random number generation and tensor representations !-------------------------------------------------------------------------------------------------- module math - use prec, only: & - pReal - use future + use prec + use future - implicit none - private - real(pReal), parameter, public :: PI = acos(-1.0_pReal) !< ratio of a circle's circumference to its diameter - real(pReal), parameter, public :: INDEG = 180.0_pReal/PI !< conversion from radian into degree - real(pReal), parameter, public :: INRAD = PI/180.0_pReal !< conversion from degree into radian - complex(pReal), parameter, public :: TWOPIIMG = (0.0_pReal,2.0_pReal)*(PI,0.0_pReal) !< Re(0.0), Im(2xPi) + implicit none + real(pReal), parameter, public :: PI = acos(-1.0_pReal) !< ratio of a circle's circumference to its diameter + real(pReal), parameter, public :: INDEG = 180.0_pReal/PI !< conversion from radian into degree + real(pReal), parameter, public :: INRAD = PI/180.0_pReal !< conversion from degree into radian + complex(pReal), parameter, public :: TWOPIIMG = cmplx(0.0_pReal,2.0_pReal*PI) !< Re(0.0), Im(2xPi) - real(pReal), dimension(3,3), parameter, public :: & - MATH_I3 = reshape([& - 1.0_pReal,0.0_pReal,0.0_pReal, & - 0.0_pReal,1.0_pReal,0.0_pReal, & - 0.0_pReal,0.0_pReal,1.0_pReal & - ],[3,3]) !< 3x3 Identity + real(pReal), dimension(3,3), parameter, public :: & + MATH_I3 = reshape([& + 1.0_pReal,0.0_pReal,0.0_pReal, & + 0.0_pReal,1.0_pReal,0.0_pReal, & + 0.0_pReal,0.0_pReal,1.0_pReal & + ],[3,3]) !< 3x3 Identity - real(pReal), dimension(6), parameter, private :: & - nrmMandel = [& - 1.0_pReal, 1.0_pReal, 1.0_pReal, & - sqrt(2.0_pReal), sqrt(2.0_pReal), sqrt(2.0_pReal) ] !< weighting for Mandel notation (forward) + real(pReal), dimension(6), parameter, private :: & + nrmMandel = [& + 1.0_pReal, 1.0_pReal, 1.0_pReal, & + sqrt(2.0_pReal), sqrt(2.0_pReal), sqrt(2.0_pReal) ] !< weighting for Mandel notation (forward) - real(pReal), dimension(6), parameter , private :: & - invnrmMandel = [& - 1.0_pReal, 1.0_pReal, 1.0_pReal, & - 1.0_pReal/sqrt(2.0_pReal), 1.0_pReal/sqrt(2.0_pReal), 1.0_pReal/sqrt(2.0_pReal) ] !< weighting for Mandel notation (backward) + real(pReal), dimension(6), parameter , private :: & + invnrmMandel = [& + 1.0_pReal, 1.0_pReal, 1.0_pReal, & + 1.0_pReal/sqrt(2.0_pReal), 1.0_pReal/sqrt(2.0_pReal), 1.0_pReal/sqrt(2.0_pReal) ] !< weighting for Mandel notation (backward) - integer, dimension (2,6), parameter, private :: & - mapNye = reshape([& - 1,1, & - 2,2, & - 3,3, & - 1,2, & - 2,3, & - 1,3 & - ],[2,6]) !< arrangement in Nye notation. + integer, dimension (2,6), parameter, private :: & + mapNye = reshape([& + 1,1, & + 2,2, & + 3,3, & + 1,2, & + 2,3, & + 1,3 & + ],[2,6]) !< arrangement in Nye notation. - integer, dimension (2,6), parameter, private :: & - mapVoigt = reshape([& - 1,1, & - 2,2, & - 3,3, & - 2,3, & - 1,3, & - 1,2 & - ],[2,6]) !< arrangement in Voigt notation + integer, dimension (2,6), parameter, private :: & + mapVoigt = reshape([& + 1,1, & + 2,2, & + 3,3, & + 2,3, & + 1,3, & + 1,2 & + ],[2,6]) !< arrangement in Voigt notation - integer, dimension (2,9), parameter, private :: & - mapPlain = reshape([& - 1,1, & - 1,2, & - 1,3, & - 2,1, & - 2,2, & - 2,3, & - 3,1, & - 3,2, & - 3,3 & - ],[2,9]) !< arrangement in Plain notation + integer, dimension (2,9), parameter, private :: & + mapPlain = reshape([& + 1,1, & + 1,2, & + 1,3, & + 2,1, & + 2,2, & + 2,3, & + 3,1, & + 3,2, & + 3,3 & + ],[2,9]) !< arrangement in Plain notation !-------------------------------------------------------------------------------------------------- ! Provide deprecated name for compatibility - interface math_crossproduct - module procedure math_cross - end interface math_crossproduct - interface math_mul3x3 - module procedure math_inner - end interface math_mul3x3 - public :: & - math_mul3x3, & - math_crossproduct + interface math_mul3x3 + module procedure math_inner + end interface math_mul3x3 + public :: & + math_mul3x3 !--------------------------------------------------------------------------------------------------- - - public :: & - math_init, & - math_qsort, & - math_expand, & - math_range, & - math_identity2nd, & - math_identity4th, & - math_civita, & - math_delta, & - math_cross, & - math_outer, & - math_inner, & - math_mul33xx33, & - math_mul3333xx33, & - math_mul3333xx3333, & - math_mul33x33, & - math_mul33x3, & - math_mul33x3_complex, & - math_mul66x6 , & - math_exp33 , & - math_transpose33, & - math_inv33, & - math_invert33, & - math_invSym3333, & - math_invert, & - math_invert2, & - math_symmetric33, & - math_symmetric66, & - math_skew33, & - math_spherical33, & - math_deviatoric33, & - math_equivStrain33, & - math_equivStress33, & - math_trace33, & - math_det33, & - math_33to9, & - math_9to33, & - math_sym33to6, & - math_6toSym33, & - math_3333to99, & - math_99to3333, & - math_sym3333to66, & - math_66toSym3333, & - math_Voigt66to3333, & - math_qRand, & - math_qMul, & - math_qDot, & - math_qConj, & - math_qInv, & - math_qRot, & - math_RtoEuler, & - math_RtoQ, & - math_EulerToR, & - math_EulerToQ, & - math_EulerAxisAngleToR, & - math_axisAngleToR, & - math_EulerAxisAngleToQ, & - math_axisAngleToQ, & - math_qToRodrig, & - math_qToEuler, & - math_qToEulerAxisAngle, & - math_qToAxisAngle, & - math_qToR, & - math_EulerMisorientation, & - math_sampleRandomOri, & - math_sampleGaussOri, & - math_sampleFiberOri, & - math_sampleGaussVar, & - math_symmetricEulers, & - math_eigenvectorBasisSym33, & - math_eigenvectorBasisSym33_log, & - math_eigenvectorBasisSym, & - math_eigenValuesVectorsSym33, & - math_eigenValuesVectorsSym, & - math_rotationalPart33, & - math_invariantsSym33, & - math_eigenvaluesSym33, & - math_factorial, & - math_binomial, & - math_multinomial, & - math_volTetrahedron, & - math_areaTriangle, & - math_rotate_forward33, & - math_rotate_backward33, & - math_rotate_forward3333, & - math_clip private :: & - math_check, & - halton + math_check contains !-------------------------------------------------------------------------------------------------- -!> @brief initialization of random seed generator +!> @brief initialization of random seed generator and internal checks !-------------------------------------------------------------------------------------------------- subroutine math_init - use numerics, only: & - randomSeed + use numerics, only: & + randomSeed - implicit none - integer :: i - real(pReal), dimension(4) :: randTest - integer :: randSize - integer, dimension(:), allocatable :: randInit - - write(6,'(/,a)') ' <<<+- math init -+>>>' + integer :: i + real(pReal), dimension(4) :: randTest + integer :: randSize + integer, dimension(:), allocatable :: randInit + + write(6,'(/,a)') ' <<<+- math init -+>>>' - call random_seed(size=randSize) - allocate(randInit(randSize)) - if (randomSeed > 0) then - randInit = randomSeed - call random_seed(put=randInit) - else - call random_seed() - call random_seed(get = randInit) - randInit(2:randSize) = randInit(1) - call random_seed(put = randInit) - endif + call random_seed(size=randSize) + allocate(randInit(randSize)) + if (randomSeed > 0) then + randInit = randomSeed + call random_seed(put=randInit) + else + call random_seed() + call random_seed(get = randInit) + randInit(2:randSize) = randInit(1) + call random_seed(put = randInit) + endif - do i = 1, 4 - call random_number(randTest(i)) - enddo + do i = 1, 4 + call random_number(randTest(i)) + enddo - write(6,'(a,I2)') ' size of random seed: ', randSize - do i = 1,randSize - write(6,'(a,I2,I14)') ' value of random seed: ', i, randInit(i) - enddo - write(6,'(a,4(/,26x,f17.14),/)') ' start of random sequence: ', randTest + write(6,'(a,I2)') ' size of random seed: ', randSize + do i = 1,randSize + write(6,'(a,I2,I14)') ' value of random seed: ', i, randInit(i) + enddo + write(6,'(a,4(/,26x,f17.14),/)') ' start of random sequence: ', randTest - call random_seed(put = randInit) - call math_check() + call random_seed(put = randInit) + call math_check end subroutine math_init @@ -219,83 +124,26 @@ end subroutine math_init !> @brief check correctness of (some) math functions !-------------------------------------------------------------------------------------------------- subroutine math_check - use prec, only: tol_math_check - use IO, only: IO_error + use IO, only: IO_error - implicit none - character(len=64) :: error_msg + character(len=64) :: error_msg - real(pReal), dimension(3,3) :: R,R2 - real(pReal), dimension(3) :: Eulers,v - real(pReal), dimension(4) :: q,q2,axisangle - - ! --- check rotation dictionary --- - - q = math_qRand() ! random quaternion - - ! +++ q -> a -> q +++ - axisangle = math_qToAxisAngle(q) - q2 = math_axisAngleToQ(axisangle(1:3),axisangle(4)) - if ( any(abs( q-q2) > tol_math_check) .and. & - any(abs(-q-q2) > tol_math_check) ) then - write (error_msg, '(a,e14.6)' ) & - 'quat -> axisAngle -> quat maximum deviation ',min(maxval(abs( q-q2)),maxval(abs(-q-q2))) - call IO_error(401,ext_msg=error_msg) - endif - - ! +++ q -> R -> q +++ - R = math_qToR(q) - q2 = math_RtoQ(R) - if ( any(abs( q-q2) > tol_math_check) .and. & - any(abs(-q-q2) > tol_math_check) ) then - write (error_msg, '(a,e14.6)' ) & - 'quat -> R -> quat maximum deviation ',min(maxval(abs( q-q2)),maxval(abs(-q-q2))) - call IO_error(401,ext_msg=error_msg) - endif - - ! +++ q -> euler -> q +++ - Eulers = math_qToEuler(q) - q2 = math_EulerToQ(Eulers) - if ( any(abs( q-q2) > tol_math_check) .and. & - any(abs(-q-q2) > tol_math_check) ) then - write (error_msg, '(a,e14.6)' ) & - 'quat -> euler -> quat maximum deviation ',min(maxval(abs( q-q2)),maxval(abs(-q-q2))) - call IO_error(401,ext_msg=error_msg) - endif - - ! +++ R -> euler -> R +++ - Eulers = math_RtoEuler(R) - R2 = math_EulerToR(Eulers) - if ( any(abs( R-R2) > tol_math_check) ) then - write (error_msg, '(a,e14.6)' ) & - 'R -> euler -> R maximum deviation ',maxval(abs( R-R2)) - call IO_error(401,ext_msg=error_msg) - endif - - ! +++ check rotation sense of q and R +++ - v = halton([2,8,5]) ! random vector - R = math_qToR(q) - if (any(abs(matmul(R,v) - math_qRot(q,v)) > tol_math_check)) then - write (error_msg, '(a)' ) 'R(q)*v has different sense than q*v' - call IO_error(401,ext_msg=error_msg) - endif - - ! +++ check vector expansion +++ - if (any(abs([1.0_pReal,2.0_pReal,2.0_pReal,3.0_pReal,3.0_pReal,3.0_pReal] - & - math_expand([1.0_pReal,2.0_pReal,3.0_pReal],[1,2,3,0])) > tol_math_check)) then - write (error_msg, '(a)' ) 'math_expand [1,2,3] by [1,2,3,0] => [1,2,2,3,3,3]' - call IO_error(401,ext_msg=error_msg) - endif - if (any(abs([1.0_pReal,2.0_pReal,2.0_pReal] - & - math_expand([1.0_pReal,2.0_pReal,3.0_pReal],[1,2])) > tol_math_check)) then - write (error_msg, '(a)' ) 'math_expand [1,2,3] by [1,2] => [1,2,2]' - call IO_error(401,ext_msg=error_msg) - endif - if (any(abs([1.0_pReal,2.0_pReal,2.0_pReal,1.0_pReal,1.0_pReal,1.0_pReal] - & - math_expand([1.0_pReal,2.0_pReal],[1,2,3])) > tol_math_check)) then - write (error_msg, '(a)' ) 'math_expand [1,2] by [1,2,3] => [1,2,2,1,1,1]' - call IO_error(401,ext_msg=error_msg) - endif + ! +++ check vector expansion +++ + if (any(abs([1.0_pReal,2.0_pReal,2.0_pReal,3.0_pReal,3.0_pReal,3.0_pReal] - & + math_expand([1.0_pReal,2.0_pReal,3.0_pReal],[1,2,3,0])) > tol_math_check)) then + write (error_msg, '(a)' ) 'math_expand [1,2,3] by [1,2,3,0] => [1,2,2,3,3,3]' + call IO_error(401,ext_msg=error_msg) + endif + if (any(abs([1.0_pReal,2.0_pReal,2.0_pReal] - & + math_expand([1.0_pReal,2.0_pReal,3.0_pReal],[1,2])) > tol_math_check)) then + write (error_msg, '(a)' ) 'math_expand [1,2,3] by [1,2] => [1,2,2]' + call IO_error(401,ext_msg=error_msg) + endif + if (any(abs([1.0_pReal,2.0_pReal,2.0_pReal,1.0_pReal,1.0_pReal,1.0_pReal] - & + math_expand([1.0_pReal,2.0_pReal],[1,2,3])) > tol_math_check)) then + write (error_msg, '(a)' ) 'math_expand [1,2] by [1,2,3] => [1,2,2,1,1,1]' + call IO_error(401,ext_msg=error_msg) + endif end subroutine math_check @@ -305,76 +153,74 @@ end subroutine math_check ! Sorting is done with respect to array(sort,:) and keeps array(/=sort,:) linked to it. ! default: sort=1 !-------------------------------------------------------------------------------------------------- -recursive subroutine math_qsort(a, istart, iend, sortDim) +recursive subroutine math_sort(a, istart, iend, sortDim) - implicit none - integer, dimension(:,:), intent(inout) :: a - integer, intent(in),optional :: istart,iend, sortDim - integer :: ipivot,s,e,d - - if(present(istart)) then - s = istart - else - s = lbound(a,2) - endif - - if(present(iend)) then - e = iend - else - e = ubound(a,2) - endif - - if(present(sortDim)) then - d = sortDim - else - d = 1 - endif + integer, dimension(:,:), intent(inout) :: a + integer, intent(in),optional :: istart,iend, sortDim + integer :: ipivot,s,e,d - if (s < e) then - ipivot = qsort_partition(a,s, e, d) - call math_qsort(a, s, ipivot-1, d) - call math_qsort(a, ipivot+1, e, d) - endif - -!-------------------------------------------------------------------------------------------------- - contains - - !------------------------------------------------------------------------------------------------- - !> @brief Partitioning required for quicksort - !------------------------------------------------------------------------------------------------- - integer function qsort_partition(a, istart, iend, sort) - - implicit none - integer, dimension(:,:), intent(inout) :: a - integer, intent(in) :: istart,iend,sort - integer, dimension(size(a,1)) :: tmp - integer :: i,j + if(present(istart)) then + s = istart + else + s = lbound(a,2) + endif - do - ! find the first element on the right side less than or equal to the pivot point - do j = iend, istart, -1 - if (a(sort,j) <= a(sort,istart)) exit - enddo - ! find the first element on the left side greater than the pivot point - do i = istart, iend - if (a(sort,i) > a(sort,istart)) exit - enddo - cross: if (i >= j) then ! if the indices cross, exchange left value with pivot and return with the partition index - tmp = a(:,istart) - a(:,istart) = a(:,j) - a(:,j) = tmp - qsort_partition = j - return - else cross ! if they do not cross, exchange values - tmp = a(:,i) - a(:,i) = a(:,j) - a(:,j) = tmp - endif cross - enddo - - end function qsort_partition + if(present(iend)) then + e = iend + else + e = ubound(a,2) + endif + + if(present(sortDim)) then + d = sortDim + else + d = 1 + endif + + if (s < e) then + ipivot = qsort_partition(a,s, e, d) + call math_sort(a, s, ipivot-1, d) + call math_sort(a, ipivot+1, e, d) + endif -end subroutine math_qsort + + contains + + !------------------------------------------------------------------------------------------------- + !> @brief Partitioning required for quicksort + !------------------------------------------------------------------------------------------------- + integer function qsort_partition(a, istart, iend, sort) + + integer, dimension(:,:), intent(inout) :: a + integer, intent(in) :: istart,iend,sort + integer, dimension(size(a,1)) :: tmp + integer :: i,j + + do + ! find the first element on the right side less than or equal to the pivot point + do j = iend, istart, -1 + if (a(sort,j) <= a(sort,istart)) exit + enddo + ! find the first element on the left side greater than the pivot point + do i = istart, iend + if (a(sort,i) > a(sort,istart)) exit + enddo + cross: if (i >= j) then ! if the indices cross, exchange left value with pivot and return with the partition index + tmp = a(:,istart) + a(:,istart) = a(:,j) + a(:,j) = tmp + qsort_partition = j + return + else cross ! if they do not cross, exchange values + tmp = a(:,i) + a(:,i) = a(:,j) + a(:,j) = tmp + endif cross + enddo + + end function qsort_partition + +end subroutine math_sort !-------------------------------------------------------------------------------------------------- @@ -384,18 +230,16 @@ end subroutine math_qsort !-------------------------------------------------------------------------------------------------- pure function math_expand(what,how) - implicit none - real(pReal), dimension(:), intent(in) :: what - integer, dimension(:), intent(in) :: how - real(pReal), dimension(sum(how)) :: math_expand - integer :: i - - if (sum(how) == 0) & - return - - do i = 1, size(how) - math_expand(sum(how(1:i-1))+1:sum(how(1:i))) = what(mod(i-1,size(what))+1) - enddo + real(pReal), dimension(:), intent(in) :: what + integer, dimension(:), intent(in) :: how + real(pReal), dimension(sum(how)) :: math_expand + integer :: i + + if (sum(how) == 0) return + + do i = 1, size(how) + math_expand(sum(how(1:i-1))+1:sum(how(1:i))) = what(mod(i-1,size(what))+1) + enddo end function math_expand @@ -405,12 +249,11 @@ end function math_expand !-------------------------------------------------------------------------------------------------- pure function math_range(N) - implicit none - integer, intent(in) :: N !< length of range - integer :: i - integer, dimension(N) :: math_range + integer, intent(in) :: N !< length of range + integer :: i + integer, dimension(N) :: math_range - math_range = [(i,i=1,N)] + math_range = [(i,i=1,N)] end function math_range @@ -420,13 +263,14 @@ end function math_range !-------------------------------------------------------------------------------------------------- pure function math_identity2nd(dimen) - implicit none - integer, intent(in) :: dimen !< tensor dimension - integer :: i - real(pReal), dimension(dimen,dimen) :: math_identity2nd + integer, intent(in) :: dimen !< tensor dimension + integer :: i + real(pReal), dimension(dimen,dimen) :: math_identity2nd - math_identity2nd = 0.0_pReal - forall(i=1:dimen) math_identity2nd(i,i) = 1.0_pReal + math_identity2nd = 0.0_pReal + do i=1, dimen + math_identity2nd(i,i) = 1.0_pReal + enddo end function math_identity2nd @@ -436,15 +280,14 @@ end function math_identity2nd !-------------------------------------------------------------------------------------------------- pure function math_identity4th(dimen) - implicit none - integer, intent(in) :: dimen !< tensor dimension - integer :: i,j,k,l - real(pReal), dimension(dimen,dimen,dimen,dimen) :: math_identity4th - real(pReal), dimension(dimen,dimen) :: identity2nd + integer, intent(in) :: dimen !< tensor dimension + integer :: i,j,k,l + real(pReal), dimension(dimen,dimen,dimen,dimen) :: math_identity4th + real(pReal), dimension(dimen,dimen) :: identity2nd - identity2nd = math_identity2nd(dimen) - forall(i=1:dimen,j=1:dimen,k=1:dimen,l=1:dimen) & - math_identity4th(i,j,k,l) = 0.5_pReal*(identity2nd(i,k)*identity2nd(j,l)+identity2nd(i,l)*identity2nd(j,k)) + identity2nd = math_identity2nd(dimen) + forall(i=1:dimen,j=1:dimen,k=1:dimen,l=1:dimen) & + math_identity4th(i,j,k,l) = 0.5_pReal*(identity2nd(i,k)*identity2nd(j,l)+identity2nd(i,l)*identity2nd(j,k)) end function math_identity4th @@ -457,16 +300,15 @@ end function math_identity4th !-------------------------------------------------------------------------------------------------- real(pReal) pure function math_civita(i,j,k) - implicit none - integer, intent(in) :: i,j,k + integer, intent(in) :: i,j,k - math_civita = 0.0_pReal - if (((i == 1).and.(j == 2).and.(k == 3)) .or. & - ((i == 2).and.(j == 3).and.(k == 1)) .or. & - ((i == 3).and.(j == 1).and.(k == 2))) math_civita = 1.0_pReal - if (((i == 1).and.(j == 3).and.(k == 2)) .or. & - ((i == 2).and.(j == 1).and.(k == 3)) .or. & - ((i == 3).and.(j == 2).and.(k == 1))) math_civita = -1.0_pReal + math_civita = 0.0_pReal + if (((i == 1).and.(j == 2).and.(k == 3)) .or. & + ((i == 2).and.(j == 3).and.(k == 1)) .or. & + ((i == 3).and.(j == 1).and.(k == 2))) math_civita = 1.0_pReal + if (((i == 1).and.(j == 3).and.(k == 2)) .or. & + ((i == 2).and.(j == 1).and.(k == 3)) .or. & + ((i == 3).and.(j == 2).and.(k == 1))) math_civita = -1.0_pReal end function math_civita @@ -479,10 +321,9 @@ end function math_civita !-------------------------------------------------------------------------------------------------- real(pReal) pure function math_delta(i,j) - implicit none - integer, intent (in) :: i,j + integer, intent (in) :: i,j - math_delta = merge(0.0_pReal, 1.0_pReal, i /= j) + math_delta = merge(0.0_pReal, 1.0_pReal, i /= j) end function math_delta @@ -492,13 +333,12 @@ end function math_delta !-------------------------------------------------------------------------------------------------- pure function math_cross(A,B) - implicit none - real(pReal), dimension(3), intent(in) :: A,B - real(pReal), dimension(3) :: math_cross + real(pReal), dimension(3), intent(in) :: A,B + real(pReal), dimension(3) :: math_cross - math_cross = [ A(2)*B(3) -A(3)*B(2), & - A(3)*B(1) -A(1)*B(3), & - A(1)*B(2) -A(2)*B(1) ] + math_cross = [ A(2)*B(3) -A(3)*B(2), & + A(3)*B(1) -A(1)*B(3), & + A(1)*B(2) -A(2)*B(1) ] end function math_cross @@ -508,12 +348,11 @@ end function math_cross !-------------------------------------------------------------------------------------------------- pure function math_outer(A,B) - implicit none - real(pReal), dimension(:), intent(in) :: A,B - real(pReal), dimension(size(A,1),size(B,1)) :: math_outer - integer :: i,j + real(pReal), dimension(:), intent(in) :: A,B + real(pReal), dimension(size(A,1),size(B,1)) :: math_outer + integer :: i,j - forall(i=1:size(A,1),j=1:size(B,1)) math_outer(i,j) = A(i)*B(j) + forall(i=1:size(A,1),j=1:size(B,1)) math_outer(i,j) = A(i)*B(j) end function math_outer @@ -523,11 +362,10 @@ end function math_outer !-------------------------------------------------------------------------------------------------- real(pReal) pure function math_inner(A,B) - implicit none - real(pReal), dimension(:), intent(in) :: A - real(pReal), dimension(size(A,1)), intent(in) :: B + real(pReal), dimension(:), intent(in) :: A + real(pReal), dimension(size(A,1)), intent(in) :: B - math_inner = sum(A*B) + math_inner = sum(A*B) end function math_inner @@ -537,13 +375,12 @@ end function math_inner !-------------------------------------------------------------------------------------------------- real(pReal) pure function math_mul33xx33(A,B) - implicit none - real(pReal), dimension(3,3), intent(in) :: A,B - integer :: i,j - real(pReal), dimension(3,3) :: C + real(pReal), dimension(3,3), intent(in) :: A,B + integer :: i,j + real(pReal), dimension(3,3) :: C - forall(i=1:3,j=1:3) C(i,j) = A(i,j) * B(i,j) - math_mul33xx33 = sum(C) + forall(i=1:3,j=1:3) C(i,j) = A(i,j) * B(i,j) + math_mul33xx33 = sum(C) end function math_mul33xx33 @@ -553,13 +390,12 @@ end function math_mul33xx33 !-------------------------------------------------------------------------------------------------- pure function math_mul3333xx33(A,B) - implicit none - real(pReal), dimension(3,3) :: math_mul3333xx33 - real(pReal), dimension(3,3,3,3), intent(in) :: A - real(pReal), dimension(3,3), intent(in) :: B - integer :: i,j + real(pReal), dimension(3,3) :: math_mul3333xx33 + real(pReal), dimension(3,3,3,3), intent(in) :: A + real(pReal), dimension(3,3), intent(in) :: B + integer :: i,j - forall(i = 1:3,j = 1:3) math_mul3333xx33(i,j) = sum(A(i,j,1:3,1:3)*B(1:3,1:3)) + forall(i = 1:3,j = 1:3) math_mul3333xx33(i,j) = sum(A(i,j,1:3,1:3)*B(1:3,1:3)) end function math_mul3333xx33 @@ -569,190 +405,71 @@ end function math_mul3333xx33 !-------------------------------------------------------------------------------------------------- pure function math_mul3333xx3333(A,B) - implicit none - integer :: i,j,k,l - real(pReal), dimension(3,3,3,3), intent(in) :: A - real(pReal), dimension(3,3,3,3), intent(in) :: B - real(pReal), dimension(3,3,3,3) :: math_mul3333xx3333 + integer :: i,j,k,l + real(pReal), dimension(3,3,3,3), intent(in) :: A + real(pReal), dimension(3,3,3,3), intent(in) :: B + real(pReal), dimension(3,3,3,3) :: math_mul3333xx3333 - forall(i = 1:3,j = 1:3, k = 1:3, l= 1:3) & - math_mul3333xx3333(i,j,k,l) = sum(A(i,j,1:3,1:3)*B(1:3,1:3,k,l)) + forall(i = 1:3,j = 1:3, k = 1:3, l= 1:3) & + math_mul3333xx3333(i,j,k,l) = sum(A(i,j,1:3,1:3)*B(1:3,1:3,k,l)) end function math_mul3333xx3333 -!-------------------------------------------------------------------------------------------------- -!> @brief matrix multiplication 33x33 = 33 -!-------------------------------------------------------------------------------------------------- -pure function math_mul33x33(A,B) - - implicit none - real(pReal), dimension(3,3) :: math_mul33x33 - real(pReal), dimension(3,3), intent(in) :: A,B - integer :: i,j - - forall(i=1:3,j=1:3) math_mul33x33(i,j) = A(i,1)*B(1,j) + A(i,2)*B(2,j) + A(i,3)*B(3,j) - -end function math_mul33x33 - - -!-------------------------------------------------------------------------------------------------- -!> @brief matrix multiplication 66x66 = 66 -!-------------------------------------------------------------------------------------------------- -pure function math_mul66x66(A,B) - - implicit none - real(pReal), dimension(6,6) :: math_mul66x66 - real(pReal), dimension(6,6), intent(in) :: A,B - integer :: i,j - - forall(i=1:6,j=1:6) & - math_mul66x66(i,j) = A(i,1)*B(1,j) + A(i,2)*B(2,j) + A(i,3)*B(3,j) & - + A(i,4)*B(4,j) + A(i,5)*B(5,j) + A(i,6)*B(6,j) - -end function math_mul66x66 - - -!-------------------------------------------------------------------------------------------------- -!> @brief matrix multiplication 99x99 = 99 -!-------------------------------------------------------------------------------------------------- -pure function math_mul99x99(A,B) - - implicit none - real(pReal), dimension(9,9) :: math_mul99x99 - real(pReal), dimension(9,9), intent(in) :: A,B - integer i,j - - forall(i=1:9,j=1:9) & - math_mul99x99(i,j) = A(i,1)*B(1,j) + A(i,2)*B(2,j) + A(i,3)*B(3,j) & - + A(i,4)*B(4,j) + A(i,5)*B(5,j) + A(i,6)*B(6,j) & - + A(i,7)*B(7,j) + A(i,8)*B(8,j) + A(i,9)*B(9,j) - -end function math_mul99x99 - - -!-------------------------------------------------------------------------------------------------- -!> @brief matrix multiplication 33x3 = 3 -!-------------------------------------------------------------------------------------------------- -pure function math_mul33x3(A,B) - - implicit none - real(pReal), dimension(3) :: math_mul33x3 - real(pReal), dimension(3,3), intent(in) :: A - real(pReal), dimension(3), intent(in) :: B - integer :: i - - forall (i=1:3) math_mul33x3(i) = sum(A(i,1:3)*B) - -end function math_mul33x3 - - -!-------------------------------------------------------------------------------------------------- -!> @brief matrix multiplication complex(33) x real(3) = complex(3) -!-------------------------------------------------------------------------------------------------- -pure function math_mul33x3_complex(A,B) - - implicit none - complex(pReal), dimension(3) :: math_mul33x3_complex - complex(pReal), dimension(3,3), intent(in) :: A - real(pReal), dimension(3), intent(in) :: B - integer :: i - - forall (i=1:3) math_mul33x3_complex(i) = sum(A(i,1:3)*cmplx(B,0.0_pReal,pReal)) - -end function math_mul33x3_complex - - -!-------------------------------------------------------------------------------------------------- -!> @brief matrix multiplication 66x6 = 6 -!-------------------------------------------------------------------------------------------------- -pure function math_mul66x6(A,B) - - implicit none - real(pReal), dimension(6) :: math_mul66x6 - real(pReal), dimension(6,6), intent(in) :: A - real(pReal), dimension(6), intent(in) :: B - integer :: i - - forall (i=1:6) math_mul66x6(i) = A(i,1)*B(1) + A(i,2)*B(2) + A(i,3)*B(3) & - + A(i,4)*B(4) + A(i,5)*B(5) + A(i,6)*B(6) - -end function math_mul66x6 - - !-------------------------------------------------------------------------------------------------- !> @brief 3x3 matrix exponential up to series approximation order n (default 5) !-------------------------------------------------------------------------------------------------- pure function math_exp33(A,n) - implicit none - integer :: i - integer, intent(in), optional :: n - real(pReal), dimension(3,3), intent(in) :: A - real(pReal), dimension(3,3) :: B, math_exp33 - real(pReal) :: invFac + integer :: i + integer, intent(in), optional :: n + real(pReal), dimension(3,3), intent(in) :: A + real(pReal), dimension(3,3) :: B, math_exp33 + real(pReal) :: invFac - B = math_I3 ! init - invFac = 1.0_pReal ! 0! - math_exp33 = B ! A^0 = eye2 + B = math_I3 ! init + invFac = 1.0_pReal ! 0! + math_exp33 = B ! A^0 = eye2 - do i = 1, merge(n,5,present(n)) - invFac = invFac/real(i,pReal) ! invfac = 1/i! - B = matmul(B,A) - math_exp33 = math_exp33 + invFac*B ! exp = SUM (A^i)/i! - enddo + do i = 1, merge(n,5,present(n)) + invFac = invFac/real(i,pReal) ! invfac = 1/i! + B = matmul(B,A) + math_exp33 = math_exp33 + invFac*B ! exp = SUM (A^i)/i! + enddo end function math_exp33 -!-------------------------------------------------------------------------------------------------- -!> @brief transposition of a 33 matrix -!-------------------------------------------------------------------------------------------------- -pure function math_transpose33(A) - - implicit none - real(pReal),dimension(3,3) :: math_transpose33 - real(pReal),dimension(3,3),intent(in) :: A - integer :: i,j - - forall(i=1:3, j=1:3) math_transpose33(i,j) = A(j,i) - -end function math_transpose33 - - !-------------------------------------------------------------------------------------------------- !> @brief Cramer inversion of 33 matrix (function) !> @details Direct Cramer inversion of matrix A. Returns all zeroes if not possible, i.e. ! if determinant is close to zero !-------------------------------------------------------------------------------------------------- pure function math_inv33(A) - use prec, only: & - dNeq0 - implicit none - real(pReal),dimension(3,3),intent(in) :: A - real(pReal) :: DetA - real(pReal),dimension(3,3) :: math_inv33 + real(pReal),dimension(3,3),intent(in) :: A + real(pReal) :: DetA + real(pReal),dimension(3,3) :: math_inv33 - math_inv33(1,1) = A(2,2) * A(3,3) - A(2,3) * A(3,2) - math_inv33(2,1) = -A(2,1) * A(3,3) + A(2,3) * A(3,1) - math_inv33(3,1) = A(2,1) * A(3,2) - A(2,2) * A(3,1) + math_inv33(1,1) = A(2,2) * A(3,3) - A(2,3) * A(3,2) + math_inv33(2,1) = -A(2,1) * A(3,3) + A(2,3) * A(3,1) + math_inv33(3,1) = A(2,1) * A(3,2) - A(2,2) * A(3,1) - DetA = A(1,1) * math_inv33(1,1) + A(1,2) * math_inv33(2,1) + A(1,3) * math_inv33(3,1) + DetA = A(1,1) * math_inv33(1,1) + A(1,2) * math_inv33(2,1) + A(1,3) * math_inv33(3,1) - if (dNeq0(DetA)) then - math_inv33(1,2) = -A(1,2) * A(3,3) + A(1,3) * A(3,2) - math_inv33(2,2) = A(1,1) * A(3,3) - A(1,3) * A(3,1) - math_inv33(3,2) = -A(1,1) * A(3,2) + A(1,2) * A(3,1) + if (dNeq0(DetA)) then + math_inv33(1,2) = -A(1,2) * A(3,3) + A(1,3) * A(3,2) + math_inv33(2,2) = A(1,1) * A(3,3) - A(1,3) * A(3,1) + math_inv33(3,2) = -A(1,1) * A(3,2) + A(1,2) * A(3,1) - math_inv33(1,3) = A(1,2) * A(2,3) - A(1,3) * A(2,2) - math_inv33(2,3) = -A(1,1) * A(2,3) + A(1,3) * A(2,1) - math_inv33(3,3) = A(1,1) * A(2,2) - A(1,2) * A(2,1) - - math_inv33 = math_inv33/DetA - else - math_inv33 = 0.0_pReal - endif + math_inv33(1,3) = A(1,2) * A(2,3) - A(1,3) * A(2,2) + math_inv33(2,3) = -A(1,1) * A(2,3) + A(1,3) * A(2,1) + math_inv33(3,3) = A(1,1) * A(2,2) - A(1,2) * A(2,1) + + math_inv33 = math_inv33/DetA + else + math_inv33 = 0.0_pReal + endif end function math_inv33 @@ -764,36 +481,33 @@ end function math_inv33 ! ToDo: Output arguments should be first !-------------------------------------------------------------------------------------------------- pure subroutine math_invert33(A, InvA, DetA, error) - use prec, only: & - dEq0 - implicit none - logical, intent(out) :: error - real(pReal),dimension(3,3),intent(in) :: A - real(pReal),dimension(3,3),intent(out) :: InvA - real(pReal), intent(out) :: DetA + logical, intent(out) :: error + real(pReal),dimension(3,3),intent(in) :: A + real(pReal),dimension(3,3),intent(out) :: InvA + real(pReal), intent(out) :: DetA - InvA(1,1) = A(2,2) * A(3,3) - A(2,3) * A(3,2) - InvA(2,1) = -A(2,1) * A(3,3) + A(2,3) * A(3,1) - InvA(3,1) = A(2,1) * A(3,2) - A(2,2) * A(3,1) + InvA(1,1) = A(2,2) * A(3,3) - A(2,3) * A(3,2) + InvA(2,1) = -A(2,1) * A(3,3) + A(2,3) * A(3,1) + InvA(3,1) = A(2,1) * A(3,2) - A(2,2) * A(3,1) - DetA = A(1,1) * InvA(1,1) + A(1,2) * InvA(2,1) + A(1,3) * InvA(3,1) + DetA = A(1,1) * InvA(1,1) + A(1,2) * InvA(2,1) + A(1,3) * InvA(3,1) - if (dEq0(DetA)) then - InvA = 0.0_pReal - error = .true. - else - InvA(1,2) = -A(1,2) * A(3,3) + A(1,3) * A(3,2) - InvA(2,2) = A(1,1) * A(3,3) - A(1,3) * A(3,1) - InvA(3,2) = -A(1,1) * A(3,2) + A(1,2) * A(3,1) + if (dEq0(DetA)) then + InvA = 0.0_pReal + error = .true. + else + InvA(1,2) = -A(1,2) * A(3,3) + A(1,3) * A(3,2) + InvA(2,2) = A(1,1) * A(3,3) - A(1,3) * A(3,1) + InvA(3,2) = -A(1,1) * A(3,2) + A(1,2) * A(3,1) - InvA(1,3) = A(1,2) * A(2,3) - A(1,3) * A(2,2) - InvA(2,3) = -A(1,1) * A(2,3) + A(1,3) * A(2,1) - InvA(3,3) = A(1,1) * A(2,2) - A(1,2) * A(2,1) + InvA(1,3) = A(1,2) * A(2,3) - A(1,3) * A(2,2) + InvA(2,3) = -A(1,1) * A(2,3) + A(1,3) * A(2,1) + InvA(3,3) = A(1,1) * A(2,2) - A(1,2) * A(2,1) - InvA = InvA/DetA - error = .false. - endif + InvA = InvA/DetA + error = .false. + endif end subroutine math_invert33 @@ -802,30 +516,29 @@ end subroutine math_invert33 !> @brief Inversion of symmetriced 3x3x3x3 tensor. !-------------------------------------------------------------------------------------------------- function math_invSym3333(A) - use IO, only: & - IO_error + use IO, only: & + IO_error - implicit none - real(pReal),dimension(3,3,3,3) :: math_invSym3333 + real(pReal),dimension(3,3,3,3) :: math_invSym3333 - real(pReal),dimension(3,3,3,3),intent(in) :: A + real(pReal),dimension(3,3,3,3),intent(in) :: A - integer :: ierr - integer, dimension(6) :: ipiv6 - real(pReal), dimension(6,6) :: temp66_Real - real(pReal), dimension(6) :: work6 - external :: & - dgetrf, & - dgetri + integer :: ierr + integer, dimension(6) :: ipiv6 + real(pReal), dimension(6,6) :: temp66_Real + real(pReal), dimension(6) :: work6 + external :: & + dgetrf, & + dgetri - temp66_real = math_sym3333to66(A) - call dgetrf(6,6,temp66_real,6,ipiv6,ierr) - call dgetri(6,temp66_real,6,ipiv6,work6,6,ierr) - if (ierr == 0) then - math_invSym3333 = math_66toSym3333(temp66_real) - else - call IO_error(400, ext_msg = 'math_invSym3333') - endif + temp66_real = math_sym3333to66(A) + call dgetrf(6,6,temp66_real,6,ipiv6,ierr) + call dgetri(6,temp66_real,6,ipiv6,work6,6,ierr) + if (ierr == 0) then + math_invSym3333 = math_66toSym3333(temp66_real) + else + call IO_error(400, ext_msg = 'math_invSym3333') + endif end function math_invSym3333 @@ -836,13 +549,12 @@ end function math_invSym3333 !-------------------------------------------------------------------------------------------------- subroutine math_invert2(InvA, error, A) - implicit none - real(pReal), dimension(:,:), intent(in) :: A - - real(pReal), dimension(size(A,1),size(A,1)), intent(out) :: invA - logical, intent(out) :: error + real(pReal), dimension(:,:), intent(in) :: A + + real(pReal), dimension(size(A,1),size(A,1)), intent(out) :: invA + logical, intent(out) :: error - call math_invert(size(A,1), A, InvA, error) + call math_invert(size(A,1), A, InvA, error) end subroutine math_invert2 @@ -854,25 +566,24 @@ end subroutine math_invert2 !-------------------------------------------------------------------------------------------------- subroutine math_invert(myDim,A, InvA, error) - implicit none - integer, intent(in) :: myDim - real(pReal), dimension(myDim,myDim), intent(in) :: A + integer, intent(in) :: myDim + real(pReal), dimension(myDim,myDim), intent(in) :: A - integer :: ierr - integer, dimension(myDim) :: ipiv - real(pReal), dimension(myDim) :: work - - real(pReal), dimension(myDim,myDim), intent(out) :: invA - logical, intent(out) :: error - external :: & - dgetrf, & - dgetri - - invA = A - call dgetrf(myDim,myDim,invA,myDim,ipiv,ierr) - call dgetri(myDim,InvA,myDim,ipiv,work,myDim,ierr) - error = merge(.true.,.false., ierr /= 0) + integer :: ierr + integer, dimension(myDim) :: ipiv + real(pReal), dimension(myDim) :: work + + real(pReal), dimension(myDim,myDim), intent(out) :: invA + logical, intent(out) :: error + external :: & + dgetrf, & + dgetri + + invA = A + call dgetrf(myDim,myDim,invA,myDim,ipiv,ierr) + call dgetri(myDim,InvA,myDim,ipiv,work,myDim,ierr) + error = merge(.true.,.false., ierr /= 0) end subroutine math_invert @@ -882,11 +593,10 @@ end subroutine math_invert !-------------------------------------------------------------------------------------------------- pure function math_symmetric33(m) - implicit none - real(pReal), dimension(3,3) :: math_symmetric33 - real(pReal), dimension(3,3), intent(in) :: m - - math_symmetric33 = 0.5_pReal * (m + transpose(m)) + real(pReal), dimension(3,3) :: math_symmetric33 + real(pReal), dimension(3,3), intent(in) :: m + + math_symmetric33 = 0.5_pReal * (m + transpose(m)) end function math_symmetric33 @@ -896,11 +606,10 @@ end function math_symmetric33 !-------------------------------------------------------------------------------------------------- pure function math_symmetric66(m) - implicit none - real(pReal), dimension(6,6) :: math_symmetric66 - real(pReal), dimension(6,6), intent(in) :: m + real(pReal), dimension(6,6) :: math_symmetric66 + real(pReal), dimension(6,6), intent(in) :: m - math_symmetric66 = 0.5_pReal * (m + transpose(m)) + math_symmetric66 = 0.5_pReal * (m + transpose(m)) end function math_symmetric66 @@ -910,11 +619,10 @@ end function math_symmetric66 !-------------------------------------------------------------------------------------------------- pure function math_skew33(m) - implicit none - real(pReal), dimension(3,3) :: math_skew33 - real(pReal), dimension(3,3), intent(in) :: m + real(pReal), dimension(3,3) :: math_skew33 + real(pReal), dimension(3,3), intent(in) :: m - math_skew33 = m - math_symmetric33(m) + math_skew33 = m - math_symmetric33(m) end function math_skew33 @@ -923,11 +631,10 @@ end function math_skew33 !-------------------------------------------------------------------------------------------------- pure function math_spherical33(m) - implicit none - real(pReal), dimension(3,3) :: math_spherical33 - real(pReal), dimension(3,3), intent(in) :: m + real(pReal), dimension(3,3) :: math_spherical33 + real(pReal), dimension(3,3), intent(in) :: m - math_spherical33 = math_I3 * math_trace33(m)/3.0_pReal + math_spherical33 = math_I3 * math_trace33(m)/3.0_pReal end function math_spherical33 @@ -937,11 +644,10 @@ end function math_spherical33 !-------------------------------------------------------------------------------------------------- pure function math_deviatoric33(m) - implicit none - real(pReal), dimension(3,3) :: math_deviatoric33 - real(pReal), dimension(3,3), intent(in) :: m + real(pReal), dimension(3,3) :: math_deviatoric33 + real(pReal), dimension(3,3), intent(in) :: m - math_deviatoric33 = m - math_spherical33(m) + math_deviatoric33 = m - math_spherical33(m) end function math_deviatoric33 @@ -951,18 +657,17 @@ end function math_deviatoric33 !-------------------------------------------------------------------------------------------------- pure function math_equivStrain33(m) - implicit none - real(pReal), dimension(3,3), intent(in) :: m - real(pReal), dimension(3) :: e,s - real(pReal) :: math_equivStrain33 + real(pReal), dimension(3,3), intent(in) :: m + real(pReal), dimension(3) :: e,s + real(pReal) :: math_equivStrain33 - e = [2.0_pReal*m(1,1)-m(2,2)-m(3,3), & - 2.0_pReal*m(2,2)-m(3,3)-m(1,1), & - 2.0_pReal*m(3,3)-m(1,1)-m(2,2)]/3.0_pReal - s = [m(1,2),m(2,3),m(1,3)]*2.0_pReal + e = [2.0_pReal*m(1,1)-m(2,2)-m(3,3), & + 2.0_pReal*m(2,2)-m(3,3)-m(1,1), & + 2.0_pReal*m(3,3)-m(1,1)-m(2,2)]/3.0_pReal + s = [m(1,2),m(2,3),m(1,3)]*2.0_pReal - math_equivStrain33 = 2.0_pReal/3.0_pReal & - * (1.50_pReal*(sum(e**2.0_pReal))+ 0.75_pReal*(sum(s**2.0_pReal)))**(0.5_pReal) + math_equivStrain33 = 2.0_pReal/3.0_pReal & + * (1.50_pReal*(sum(e**2.0_pReal))+ 0.75_pReal*(sum(s**2.0_pReal)))**(0.5_pReal) end function math_equivStrain33 @@ -972,19 +677,18 @@ end function math_equivStrain33 !-------------------------------------------------------------------------------------------------- pure function math_equivStress33(m) - implicit none - real(pReal), dimension(3,3), intent(in) :: m - real(pReal) :: math_equivStress33 + real(pReal), dimension(3,3), intent(in) :: m + real(pReal) :: math_equivStress33 - math_equivStress33 =( ( (m(1,1)-m(2,2))**2.0_pReal + & - (m(2,2)-m(3,3))**2.0_pReal + & - (m(3,3)-m(1,1))**2.0_pReal + & - 6.0_pReal*( m(1,2)**2.0_pReal + & - m(2,3)**2.0_pReal + & - m(1,3)**2.0_pReal & - ) & - )**0.5_pReal & - )/sqrt(2.0_pReal) + math_equivStress33 =( ( (m(1,1)-m(2,2))**2.0_pReal + & + (m(2,2)-m(3,3))**2.0_pReal + & + (m(3,3)-m(1,1))**2.0_pReal + & + 6.0_pReal*( m(1,2)**2.0_pReal + & + m(2,3)**2.0_pReal + & + m(1,3)**2.0_pReal & + ) & + )**0.5_pReal & + )/sqrt(2.0_pReal) end function math_equivStress33 @@ -994,10 +698,9 @@ end function math_equivStress33 !-------------------------------------------------------------------------------------------------- real(pReal) pure function math_trace33(m) - implicit none - real(pReal), dimension(3,3), intent(in) :: m + real(pReal), dimension(3,3), intent(in) :: m - math_trace33 = m(1,1) + m(2,2) + m(3,3) + math_trace33 = m(1,1) + m(2,2) + m(3,3) end function math_trace33 @@ -1007,12 +710,11 @@ end function math_trace33 !-------------------------------------------------------------------------------------------------- real(pReal) pure function math_det33(m) - implicit none - real(pReal), dimension(3,3), intent(in) :: m - - math_det33 = m(1,1)* (m(2,2)*m(3,3)-m(2,3)*m(3,2)) & - - m(1,2)* (m(2,1)*m(3,3)-m(2,3)*m(3,1)) & - + m(1,3)* (m(2,1)*m(3,2)-m(2,2)*m(3,1)) + real(pReal), dimension(3,3), intent(in) :: m + + math_det33 = m(1,1)* (m(2,2)*m(3,3)-m(2,3)*m(3,2)) & + - m(1,2)* (m(2,1)*m(3,3)-m(2,3)*m(3,1)) & + + m(1,3)* (m(2,1)*m(3,2)-m(2,2)*m(3,1)) end function math_det33 @@ -1022,8 +724,7 @@ end function math_det33 !-------------------------------------------------------------------------------------------------- real(pReal) pure function math_detSym33(m) - implicit none - real(pReal), dimension(3,3), intent(in) :: m + real(pReal), dimension(3,3), intent(in) :: m math_detSym33 = -(m(1,1)*m(2,3)**2 + m(2,2)*m(1,3)**2 + m(3,3)*m(1,2)**2) & + m(1,1)*m(2,2)*m(3,3) + 2.0_pReal * m(1,2)*m(1,3)*m(2,3) @@ -1036,13 +737,14 @@ end function math_detSym33 !-------------------------------------------------------------------------------------------------- pure function math_33to9(m33) - implicit none - real(pReal), dimension(9) :: math_33to9 - real(pReal), dimension(3,3), intent(in) :: m33 - - integer :: i - - forall(i=1:9) math_33to9(i) = m33(mapPlain(1,i),mapPlain(2,i)) + real(pReal), dimension(9) :: math_33to9 + real(pReal), dimension(3,3), intent(in) :: m33 + + integer :: i + + do i = 1, 9 + math_33to9(i) = m33(mapPlain(1,i),mapPlain(2,i)) + enddo end function math_33to9 @@ -1052,13 +754,14 @@ end function math_33to9 !-------------------------------------------------------------------------------------------------- pure function math_9to33(v9) - implicit none - real(pReal), dimension(3,3) :: math_9to33 - real(pReal), dimension(9), intent(in) :: v9 - - integer :: i + real(pReal), dimension(3,3) :: math_9to33 + real(pReal), dimension(9), intent(in) :: v9 + + integer :: i - forall(i=1:9) math_9to33(mapPlain(1,i),mapPlain(2,i)) = v9(i) + do i = 1, 9 + math_9to33(mapPlain(1,i),mapPlain(2,i)) = v9(i) + enddo end function math_9to33 @@ -1071,21 +774,22 @@ end function math_9to33 !-------------------------------------------------------------------------------------------------- pure function math_sym33to6(m33,weighted) - implicit none - real(pReal), dimension(6) :: math_sym33to6 - real(pReal), dimension(3,3), intent(in) :: m33 - logical, optional, intent(in) :: weighted - - real(pReal), dimension(6) :: w - integer :: i - - if(present(weighted)) then - w = merge(nrmMandel,1.0_pReal,weighted) - else - w = nrmMandel - endif + real(pReal), dimension(6) :: math_sym33to6 + real(pReal), dimension(3,3), intent(in) :: m33 + logical, optional, intent(in) :: weighted + + real(pReal), dimension(6) :: w + integer :: i + + if(present(weighted)) then + w = merge(nrmMandel,1.0_pReal,weighted) + else + w = nrmMandel + endif - forall(i=1:6) math_sym33to6(i) = w(i)*m33(mapNye(1,i),mapNye(2,i)) + do i = 1, 6 + math_sym33to6(i) = w(i)*m33(mapNye(1,i),mapNye(2,i)) + enddo end function math_sym33to6 @@ -1098,24 +802,23 @@ end function math_sym33to6 !-------------------------------------------------------------------------------------------------- pure function math_6toSym33(v6,weighted) - implicit none - real(pReal), dimension(3,3) :: math_6toSym33 - real(pReal), dimension(6), intent(in) :: v6 - logical, optional, intent(in) :: weighted + real(pReal), dimension(3,3) :: math_6toSym33 + real(pReal), dimension(6), intent(in) :: v6 + logical, optional, intent(in) :: weighted - real(pReal), dimension(6) :: w - integer :: i - - if(present(weighted)) then - w = merge(invnrmMandel,1.0_pReal,weighted) - else - w = invnrmMandel - endif + real(pReal), dimension(6) :: w + integer :: i + + if(present(weighted)) then + w = merge(invnrmMandel,1.0_pReal,weighted) + else + w = invnrmMandel + endif - do i=1,6 - math_6toSym33(mapNye(1,i),mapNye(2,i)) = w(i)*v6(i) - math_6toSym33(mapNye(2,i),mapNye(1,i)) = w(i)*v6(i) - enddo + do i=1,6 + math_6toSym33(mapNye(1,i),mapNye(2,i)) = w(i)*v6(i) + math_6toSym33(mapNye(2,i),mapNye(1,i)) = w(i)*v6(i) + enddo end function math_6toSym33 @@ -1125,14 +828,14 @@ end function math_6toSym33 !-------------------------------------------------------------------------------------------------- pure function math_3333to99(m3333) - implicit none - real(pReal), dimension(9,9) :: math_3333to99 - real(pReal), dimension(3,3,3,3), intent(in) :: m3333 + real(pReal), dimension(9,9) :: math_3333to99 + real(pReal), dimension(3,3,3,3), intent(in) :: m3333 - integer :: i,j + integer :: i,j - forall(i=1:9,j=1:9) & - math_3333to99(i,j) = m3333(mapPlain(1,i),mapPlain(2,i),mapPlain(1,j),mapPlain(2,j)) + do i=1,9; do j=1,9 + math_3333to99(i,j) = m3333(mapPlain(1,i),mapPlain(2,i),mapPlain(1,j),mapPlain(2,j)) + enddo; enddo end function math_3333to99 @@ -1142,14 +845,14 @@ end function math_3333to99 !-------------------------------------------------------------------------------------------------- pure function math_99to3333(m99) - implicit none - real(pReal), dimension(3,3,3,3) :: math_99to3333 - real(pReal), dimension(9,9), intent(in) :: m99 + real(pReal), dimension(3,3,3,3) :: math_99to3333 + real(pReal), dimension(9,9), intent(in) :: m99 - integer :: i,j + integer :: i,j - forall(i=1:9,j=1:9) & - math_99to3333(mapPlain(1,i),mapPlain(2,i),mapPlain(1,j),mapPlain(2,j)) = m99(i,j) + do i=1,9; do j=1,9 + math_99to3333(mapPlain(1,i),mapPlain(2,i),mapPlain(1,j),mapPlain(2,j)) = m99(i,j) + enddo; enddo end function math_99to3333 @@ -1162,22 +865,22 @@ end function math_99to3333 !-------------------------------------------------------------------------------------------------- pure function math_sym3333to66(m3333,weighted) - implicit none - real(pReal), dimension(6,6) :: math_sym3333to66 - real(pReal), dimension(3,3,3,3), intent(in) :: m3333 - logical, optional, intent(in) :: weighted + real(pReal), dimension(6,6) :: math_sym3333to66 + real(pReal), dimension(3,3,3,3), intent(in) :: m3333 + logical, optional, intent(in) :: weighted - real(pReal), dimension(6) :: w - integer :: i,j - - if(present(weighted)) then - w = merge(nrmMandel,1.0_pReal,weighted) - else - w = nrmMandel - endif + real(pReal), dimension(6) :: w + integer :: i,j + + if(present(weighted)) then + w = merge(nrmMandel,1.0_pReal,weighted) + else + w = nrmMandel + endif - forall(i=1:6,j=1:6) & - math_sym3333to66(i,j) = w(i)*w(j)*m3333(mapNye(1,i),mapNye(2,i),mapNye(1,j),mapNye(2,j)) + do i=1,6; do j=1,6 + math_sym3333to66(i,j) = w(i)*w(j)*m3333(mapNye(1,i),mapNye(2,i),mapNye(1,j),mapNye(2,j)) + enddo; enddo end function math_sym3333to66 @@ -1190,26 +893,25 @@ end function math_sym3333to66 !-------------------------------------------------------------------------------------------------- pure function math_66toSym3333(m66,weighted) - implicit none - real(pReal), dimension(3,3,3,3) :: math_66toSym3333 - real(pReal), dimension(6,6), intent(in) :: m66 - logical, optional, intent(in) :: weighted - - real(pReal), dimension(6) :: w - integer :: i,j - - if(present(weighted)) then - w = merge(invnrmMandel,1.0_pReal,weighted) - else - w = invnrmMandel - endif + real(pReal), dimension(3,3,3,3) :: math_66toSym3333 + real(pReal), dimension(6,6), intent(in) :: m66 + logical, optional, intent(in) :: weighted + + real(pReal), dimension(6) :: w + integer :: i,j + + if(present(weighted)) then + w = merge(invnrmMandel,1.0_pReal,weighted) + else + w = invnrmMandel + endif - do i=1,6; do j=1, 6 - math_66toSym3333(mapNye(1,i),mapNye(2,i),mapNye(1,j),mapNye(2,j)) = w(i)*w(j)*m66(i,j) - math_66toSym3333(mapNye(2,i),mapNye(1,i),mapNye(1,j),mapNye(2,j)) = w(i)*w(j)*m66(i,j) - math_66toSym3333(mapNye(1,i),mapNye(2,i),mapNye(2,j),mapNye(1,j)) = w(i)*w(j)*m66(i,j) - math_66toSym3333(mapNye(2,i),mapNye(1,i),mapNye(2,j),mapNye(1,j)) = w(i)*w(j)*m66(i,j) - enddo; enddo + do i=1,6; do j=1,6 + math_66toSym3333(mapNye(1,i),mapNye(2,i),mapNye(1,j),mapNye(2,j)) = w(i)*w(j)*m66(i,j) + math_66toSym3333(mapNye(2,i),mapNye(1,i),mapNye(1,j),mapNye(2,j)) = w(i)*w(j)*m66(i,j) + math_66toSym3333(mapNye(1,i),mapNye(2,i),mapNye(2,j),mapNye(1,j)) = w(i)*w(j)*m66(i,j) + math_66toSym3333(mapNye(2,i),mapNye(1,i),mapNye(2,j),mapNye(1,j)) = w(i)*w(j)*m66(i,j) + enddo; enddo end function math_66toSym3333 @@ -1219,142 +921,71 @@ end function math_66toSym3333 !-------------------------------------------------------------------------------------------------- pure function math_Voigt66to3333(m66) - implicit none - real(pReal), dimension(3,3,3,3) :: math_Voigt66to3333 - real(pReal), dimension(6,6), intent(in) :: m66 - integer :: i,j + real(pReal), dimension(3,3,3,3) :: math_Voigt66to3333 + real(pReal), dimension(6,6), intent(in) :: m66 + integer :: i,j - do i=1,6; do j=1, 6 - math_Voigt66to3333(mapVoigt(1,i),mapVoigt(2,i),mapVoigt(1,j),mapVoigt(2,j)) = m66(i,j) - math_Voigt66to3333(mapVoigt(2,i),mapVoigt(1,i),mapVoigt(1,j),mapVoigt(2,j)) = m66(i,j) - math_Voigt66to3333(mapVoigt(1,i),mapVoigt(2,i),mapVoigt(2,j),mapVoigt(1,j)) = m66(i,j) - math_Voigt66to3333(mapVoigt(2,i),mapVoigt(1,i),mapVoigt(2,j),mapVoigt(1,j)) = m66(i,j) - enddo; enddo + do i=1,6; do j=1, 6 + math_Voigt66to3333(mapVoigt(1,i),mapVoigt(2,i),mapVoigt(1,j),mapVoigt(2,j)) = m66(i,j) + math_Voigt66to3333(mapVoigt(2,i),mapVoigt(1,i),mapVoigt(1,j),mapVoigt(2,j)) = m66(i,j) + math_Voigt66to3333(mapVoigt(1,i),mapVoigt(2,i),mapVoigt(2,j),mapVoigt(1,j)) = m66(i,j) + math_Voigt66to3333(mapVoigt(2,i),mapVoigt(1,i),mapVoigt(2,j),mapVoigt(1,j)) = m66(i,j) + enddo; enddo end function math_Voigt66to3333 -!-------------------------------------------------------------------------------------------------- -!> @brief random quaternion -! http://math.stackexchange.com/questions/131336/uniform-random-quaternion-in-a-restricted-angle-range -! K. Shoemake. Uniform random rotations. In D. Kirk, editor, Graphics Gems III, pages 124-132. -! Academic, New York, 1992. -!-------------------------------------------------------------------------------------------------- -function math_qRand() - - implicit none - real(pReal), dimension(4) :: math_qRand - real(pReal), dimension(3) :: rnd - - rnd = halton([8,4,9]) - math_qRand = [cos(2.0_pReal*PI*rnd(1))*sqrt(rnd(3)), & - sin(2.0_pReal*PI*rnd(2))*sqrt(1.0_pReal-rnd(3)), & - cos(2.0_pReal*PI*rnd(2))*sqrt(1.0_pReal-rnd(3)), & - sin(2.0_pReal*PI*rnd(1))*sqrt(rnd(3))] - -end function math_qRand - - !-------------------------------------------------------------------------------------------------- !> @brief quaternion multiplication q1xq2 = q12 !-------------------------------------------------------------------------------------------------- pure function math_qMul(A,B) - implicit none - real(pReal), dimension(4) :: math_qMul - real(pReal), dimension(4), intent(in) :: A, B + real(pReal), dimension(4) :: math_qMul + real(pReal), dimension(4), intent(in) :: A, B - math_qMul = [ A(1)*B(1) - A(2)*B(2) - A(3)*B(3) - A(4)*B(4), & - A(1)*B(2) + A(2)*B(1) + A(3)*B(4) - A(4)*B(3), & - A(1)*B(3) - A(2)*B(4) + A(3)*B(1) + A(4)*B(2), & - A(1)*B(4) + A(2)*B(3) - A(3)*B(2) + A(4)*B(1) ] + math_qMul = [ A(1)*B(1) - A(2)*B(2) - A(3)*B(3) - A(4)*B(4), & + A(1)*B(2) + A(2)*B(1) + A(3)*B(4) - A(4)*B(3), & + A(1)*B(3) - A(2)*B(4) + A(3)*B(1) + A(4)*B(2), & + A(1)*B(4) + A(2)*B(3) - A(3)*B(2) + A(4)*B(1) ] end function math_qMul -!-------------------------------------------------------------------------------------------------- -!> @brief quaternion dotproduct -!-------------------------------------------------------------------------------------------------- -real(pReal) pure function math_qDot(A,B) - - implicit none - real(pReal), dimension(4), intent(in) :: A, B - - math_qDot = sum(A*B) - -end function math_qDot - - !-------------------------------------------------------------------------------------------------- !> @brief quaternion conjugation !-------------------------------------------------------------------------------------------------- pure function math_qConj(Q) - implicit none - real(pReal), dimension(4) :: math_qConj - real(pReal), dimension(4), intent(in) :: Q + real(pReal), dimension(4) :: math_qConj + real(pReal), dimension(4), intent(in) :: Q - math_qConj = [Q(1), -Q(2:4)] + math_qConj = [Q(1), -Q(2:4)] end function math_qConj -!-------------------------------------------------------------------------------------------------- -!> @brief quaternion norm -!-------------------------------------------------------------------------------------------------- -real(pReal) pure function math_qNorm(Q) - - implicit none - real(pReal), dimension(4), intent(in) :: Q - - math_qNorm = norm2(Q) - -end function math_qNorm - - -!-------------------------------------------------------------------------------------------------- -!> @brief quaternion inversion -!-------------------------------------------------------------------------------------------------- -pure function math_qInv(Q) - use prec, only: & - dNeq0 - - implicit none - real(pReal), dimension(4), intent(in) :: Q - real(pReal), dimension(4) :: math_qInv - real(pReal) :: squareNorm - - math_qInv = 0.0_pReal - - squareNorm = math_qDot(Q,Q) - if (dNeq0(squareNorm)) math_qInv = math_qConj(Q) / squareNorm - -end function math_qInv - - !-------------------------------------------------------------------------------------------------- !> @brief action of a quaternion on a vector (rotate vector v with Q) !-------------------------------------------------------------------------------------------------- pure function math_qRot(Q,v) - implicit none - real(pReal), dimension(4), intent(in) :: Q - real(pReal), dimension(3), intent(in) :: v - real(pReal), dimension(3) :: math_qRot - real(pReal), dimension(4,4) :: T - integer :: i, j + real(pReal), dimension(4), intent(in) :: Q + real(pReal), dimension(3), intent(in) :: v + real(pReal), dimension(3) :: math_qRot + real(pReal), dimension(4,4) :: T + integer :: i, j - do i = 1,4 - do j = 1,i - T(i,j) = Q(i) * Q(j) - enddo - enddo + do i = 1,4 + do j = 1,i + T(i,j) = Q(i) * Q(j) + enddo + enddo - math_qRot = [-v(1)*(T(3,3)+T(4,4)) + v(2)*(T(3,2)-T(4,1)) + v(3)*(T(4,2)+T(3,1)), & - v(1)*(T(3,2)+T(4,1)) - v(2)*(T(2,2)+T(4,4)) + v(3)*(T(4,3)-T(2,1)), & - v(1)*(T(4,2)-T(3,1)) + v(2)*(T(4,3)+T(2,1)) - v(3)*(T(2,2)+T(3,3))] + math_qRot = [-v(1)*(T(3,3)+T(4,4)) + v(2)*(T(3,2)-T(4,1)) + v(3)*(T(4,2)+T(3,1)), & + v(1)*(T(3,2)+T(4,1)) - v(2)*(T(2,2)+T(4,4)) + v(3)*(T(4,3)-T(2,1)), & + v(1)*(T(4,2)-T(3,1)) + v(2)*(T(4,3)+T(2,1)) - v(3)*(T(2,2)+T(3,3))] - math_qRot = 2.0_pReal * math_qRot + v + math_qRot = 2.0_pReal * math_qRot + v end function math_qRot @@ -1368,28 +999,26 @@ end function math_qRot !-------------------------------------------------------------------------------------------------- pure function math_RtoEuler(R) - implicit none - real(pReal), dimension (3,3), intent(in) :: R - real(pReal), dimension(3) :: math_RtoEuler - real(pReal) :: sqhkl, squvw, sqhk + real(pReal), dimension (3,3), intent(in) :: R + real(pReal), dimension(3) :: math_RtoEuler + real(pReal) :: sqhkl, squvw, sqhk - sqhkl=sqrt(R(1,3)*R(1,3)+R(2,3)*R(2,3)+R(3,3)*R(3,3)) - squvw=sqrt(R(1,1)*R(1,1)+R(2,1)*R(2,1)+R(3,1)*R(3,1)) - sqhk =sqrt(R(1,3)*R(1,3)+R(2,3)*R(2,3)) + sqhkl=sqrt(R(1,3)*R(1,3)+R(2,3)*R(2,3)+R(3,3)*R(3,3)) + squvw=sqrt(R(1,1)*R(1,1)+R(2,1)*R(2,1)+R(3,1)*R(3,1)) + sqhk =sqrt(R(1,3)*R(1,3)+R(2,3)*R(2,3)) -! calculate PHI - math_RtoEuler(2) = acos(math_clip(R(3,3)/sqhkl,-1.0_pReal, 1.0_pReal)) + math_RtoEuler(2) = acos(math_clip(R(3,3)/sqhkl,-1.0_pReal, 1.0_pReal)) - if((math_RtoEuler(2) < 1.0e-8_pReal) .or. (pi-math_RtoEuler(2) < 1.0e-8_pReal)) then - math_RtoEuler(3) = 0.0_pReal - math_RtoEuler(1) = acos(math_clip(R(1,1)/squvw, -1.0_pReal, 1.0_pReal)) - if(R(2,1) > 0.0_pReal) math_RtoEuler(1) = 2.0_pReal*pi-math_RtoEuler(1) - else - math_RtoEuler(3) = acos(math_clip(R(2,3)/sqhk, -1.0_pReal, 1.0_pReal)) - if(R(1,3) < 0.0) math_RtoEuler(3) = 2.0_pReal*pi-math_RtoEuler(3) - math_RtoEuler(1) = acos(math_clip(-R(3,2)/sin(math_RtoEuler(2)), -1.0_pReal, 1.0_pReal)) - if(R(3,1) < 0.0) math_RtoEuler(1) = 2.0_pReal*pi-math_RtoEuler(1) - end if + if((math_RtoEuler(2) < 1.0e-8_pReal) .or. (pi-math_RtoEuler(2) < 1.0e-8_pReal)) then + math_RtoEuler(3) = 0.0_pReal + math_RtoEuler(1) = acos(math_clip(R(1,1)/squvw, -1.0_pReal, 1.0_pReal)) + if(R(2,1) > 0.0_pReal) math_RtoEuler(1) = 2.0_pReal*pi-math_RtoEuler(1) + else + math_RtoEuler(3) = acos(math_clip(R(2,3)/sqhk, -1.0_pReal, 1.0_pReal)) + if(R(1,3) < 0.0) math_RtoEuler(3) = 2.0_pReal*pi-math_RtoEuler(3) + math_RtoEuler(1) = acos(math_clip(-R(3,2)/sin(math_RtoEuler(2)), -1.0_pReal, 1.0_pReal)) + if(R(3,1) < 0.0) math_RtoEuler(1) = 2.0_pReal*pi-math_RtoEuler(1) + end if end function math_RtoEuler @@ -1400,50 +1029,49 @@ end function math_RtoEuler !-------------------------------------------------------------------------------------------------- pure function math_RtoQ(R) - implicit none - real(pReal), dimension(3,3), intent(in) :: R - real(pReal), dimension(4) :: absQ, math_RtoQ - real(pReal) :: max_absQ - integer, dimension(1) :: largest + real(pReal), dimension(3,3), intent(in) :: R + real(pReal), dimension(4) :: absQ, math_RtoQ + real(pReal) :: max_absQ + integer, dimension(1) :: largest - math_RtoQ = 0.0_pReal + math_RtoQ = 0.0_pReal - absQ = [+ R(1,1) + R(2,2) + R(3,3), & - + R(1,1) - R(2,2) - R(3,3), & - - R(1,1) + R(2,2) - R(3,3), & - - R(1,1) - R(2,2) + R(3,3)] + 1.0_pReal + absQ = [+ R(1,1) + R(2,2) + R(3,3), & + + R(1,1) - R(2,2) - R(3,3), & + - R(1,1) + R(2,2) - R(3,3), & + - R(1,1) - R(2,2) + R(3,3)] + 1.0_pReal - largest = maxloc(absQ) + largest = maxloc(absQ) - largestComponent: select case(largest(1)) - case (1) largestComponent - !1---------------------------------- - math_RtoQ(2) = R(3,2) - R(2,3) - math_RtoQ(3) = R(1,3) - R(3,1) - math_RtoQ(4) = R(2,1) - R(1,2) - - case (2) largestComponent - math_RtoQ(1) = R(3,2) - R(2,3) - !2---------------------------------- - math_RtoQ(3) = R(2,1) + R(1,2) - math_RtoQ(4) = R(1,3) + R(3,1) - - case (3) largestComponent - math_RtoQ(1) = R(1,3) - R(3,1) - math_RtoQ(2) = R(2,1) + R(1,2) - !3---------------------------------- - math_RtoQ(4) = R(3,2) + R(2,3) - - case (4) largestComponent - math_RtoQ(1) = R(2,1) - R(1,2) - math_RtoQ(2) = R(1,3) + R(3,1) - math_RtoQ(3) = R(2,3) + R(3,2) - !4---------------------------------- - end select largestComponent + largestComponent: select case(largest(1)) + case (1) largestComponent + !1---------------------------------- + math_RtoQ(2) = R(3,2) - R(2,3) + math_RtoQ(3) = R(1,3) - R(3,1) + math_RtoQ(4) = R(2,1) - R(1,2) + + case (2) largestComponent + math_RtoQ(1) = R(3,2) - R(2,3) + !2---------------------------------- + math_RtoQ(3) = R(2,1) + R(1,2) + math_RtoQ(4) = R(1,3) + R(3,1) + + case (3) largestComponent + math_RtoQ(1) = R(1,3) - R(3,1) + math_RtoQ(2) = R(2,1) + R(1,2) + !3---------------------------------- + math_RtoQ(4) = R(3,2) + R(2,3) + + case (4) largestComponent + math_RtoQ(1) = R(2,1) - R(1,2) + math_RtoQ(2) = R(1,3) + R(3,1) + math_RtoQ(3) = R(2,3) + R(3,2) + !4---------------------------------- + end select largestComponent - max_absQ = 0.5_pReal * sqrt(absQ(largest(1))) - math_RtoQ = math_RtoQ * 0.25_pReal / max_absQ - math_RtoQ(largest(1)) = max_absQ + max_absQ = 0.5_pReal * sqrt(absQ(largest(1))) + math_RtoQ = math_RtoQ * 0.25_pReal / max_absQ + math_RtoQ(largest(1)) = max_absQ end function math_RtoQ @@ -1457,64 +1085,34 @@ end function math_RtoQ !-------------------------------------------------------------------------------------------------- pure function math_EulerToR(Euler) - implicit none - real(pReal), dimension(3), intent(in) :: Euler - real(pReal), dimension(3,3) :: math_EulerToR - real(pReal) :: c1, C, c2, s1, S, s2 + real(pReal), dimension(3), intent(in) :: Euler + real(pReal), dimension(3,3) :: math_EulerToR + real(pReal) :: c1, C, c2, s1, S, s2 - c1 = cos(Euler(1)) - C = cos(Euler(2)) - c2 = cos(Euler(3)) - s1 = sin(Euler(1)) - S = sin(Euler(2)) - s2 = sin(Euler(3)) + c1 = cos(Euler(1)) + C = cos(Euler(2)) + c2 = cos(Euler(3)) + s1 = sin(Euler(1)) + S = sin(Euler(2)) + s2 = sin(Euler(3)) - math_EulerToR(1,1) = c1*c2 -s1*C*s2 - math_EulerToR(1,2) = -c1*s2 -s1*C*c2 - math_EulerToR(1,3) = s1*S + math_EulerToR(1,1) = c1*c2 -s1*C*s2 + math_EulerToR(1,2) = -c1*s2 -s1*C*c2 + math_EulerToR(1,3) = s1*S - math_EulerToR(2,1) = s1*c2 +c1*C*s2 - math_EulerToR(2,2) = -s1*s2 +c1*C*c2 - math_EulerToR(2,3) = -c1*S + math_EulerToR(2,1) = s1*c2 +c1*C*s2 + math_EulerToR(2,2) = -s1*s2 +c1*C*c2 + math_EulerToR(2,3) = -c1*S - math_EulerToR(3,1) = S*s2 - math_EulerToR(3,2) = S*c2 - math_EulerToR(3,3) = C - - math_EulerToR = transpose(math_EulerToR) ! convert to passive rotation + math_EulerToR(3,1) = S*s2 + math_EulerToR(3,2) = S*c2 + math_EulerToR(3,3) = C + + math_EulerToR = transpose(math_EulerToR) ! convert to passive rotation end function math_EulerToR -!-------------------------------------------------------------------------------------------------- -!> @brief quaternion (w+ix+jy+kz) from Bunge-Euler (3-1-3) angles (in radians) -!> @details rotation matrix is meant to represent a PASSIVE rotation, composed of INTRINSIC -!> @details rotations around the axes of the details rotating reference frame. -!> @details similar to eu2qu from "D Rowenhorst et al. Consistent representations of and -!> @details conversions between 3D rotations, Model. Simul. Mater. Sci. Eng. 23-8 (2015)", but -!> @details Q is conjucated and Q is not reversed for Q(0) < 0. -!-------------------------------------------------------------------------------------------------- -pure function math_EulerToQ(eulerangles) - - implicit none - real(pReal), dimension(3), intent(in) :: eulerangles - real(pReal), dimension(4) :: math_EulerToQ - real(pReal) :: c, s, sigma, delta - - c = cos(0.5_pReal * eulerangles(2)) - s = sin(0.5_pReal * eulerangles(2)) - sigma = 0.5_pReal * (eulerangles(1)+eulerangles(3)) - delta = 0.5_pReal * (eulerangles(1)-eulerangles(3)) - - math_EulerToQ= [c * cos(sigma), & - s * cos(delta), & - s * sin(delta), & - c * sin(sigma) ] - math_EulerToQ = math_qConj(math_EulerToQ) ! convert to passive rotation - -end function math_EulerToQ - - !-------------------------------------------------------------------------------------------------- !> @brief rotation matrix from axis and angle (in radians) !> @details rotation matrix is meant to represent a ACTIVE rotation @@ -1525,434 +1123,98 @@ end function math_EulerToQ !-------------------------------------------------------------------------------------------------- pure function math_axisAngleToR(axis,omega) - implicit none - real(pReal), dimension(3,3) :: math_axisAngleToR - real(pReal), dimension(3), intent(in) :: axis - real(pReal), intent(in) :: omega - real(pReal), dimension(3) :: n - real(pReal) :: norm,s,c,c1 + real(pReal), dimension(3,3) :: math_axisAngleToR + real(pReal), dimension(3), intent(in) :: axis + real(pReal), intent(in) :: omega + real(pReal), dimension(3) :: n + real(pReal) :: norm,s,c,c1 - norm = norm2(axis) - wellDefined: if (norm > 1.0e-8_pReal) then - n = axis/norm ! normalize axis to be sure + norm = norm2(axis) + wellDefined: if (norm > 1.0e-8_pReal) then + n = axis/norm ! normalize axis to be sure - s = sin(omega) - c = cos(omega) - c1 = 1.0_pReal - c + s = sin(omega) + c = cos(omega) + c1 = 1.0_pReal - c - math_axisAngleToR(1,1) = c + c1*n(1)**2.0_pReal - math_axisAngleToR(1,2) = c1*n(1)*n(2) - s*n(3) - math_axisAngleToR(1,3) = c1*n(1)*n(3) + s*n(2) - - math_axisAngleToR(2,1) = c1*n(1)*n(2) + s*n(3) - math_axisAngleToR(2,2) = c + c1*n(2)**2.0_pReal - math_axisAngleToR(2,3) = c1*n(2)*n(3) - s*n(1) - - math_axisAngleToR(3,1) = c1*n(1)*n(3) - s*n(2) - math_axisAngleToR(3,2) = c1*n(2)*n(3) + s*n(1) - math_axisAngleToR(3,3) = c + c1*n(3)**2.0_pReal - else wellDefined - math_axisAngleToR = math_I3 - endif wellDefined + math_axisAngleToR(1,1) = c + c1*n(1)**2.0_pReal + math_axisAngleToR(1,2) = c1*n(1)*n(2) - s*n(3) + math_axisAngleToR(1,3) = c1*n(1)*n(3) + s*n(2) + + math_axisAngleToR(2,1) = c1*n(1)*n(2) + s*n(3) + math_axisAngleToR(2,2) = c + c1*n(2)**2.0_pReal + math_axisAngleToR(2,3) = c1*n(2)*n(3) - s*n(1) + + math_axisAngleToR(3,1) = c1*n(1)*n(3) - s*n(2) + math_axisAngleToR(3,2) = c1*n(2)*n(3) + s*n(1) + math_axisAngleToR(3,3) = c + c1*n(3)**2.0_pReal + else wellDefined + math_axisAngleToR = math_I3 + endif wellDefined end function math_axisAngleToR -!-------------------------------------------------------------------------------------------------- -!> @brief rotation matrix from axis and angle (in radians) -!> @details rotation matrix is meant to represent a PASSIVE rotation -!> @details (see http://en.wikipedia.org/wiki/Euler_angles for definitions) -!> @details eq-uivalent to eu2qu (P=+1) from "D Rowenhorst et al. Consistent representations of and -!> @details conversions between 3D rotations, Model. Simul. Mater. Sci. Eng. 23-8 (2015)" -!-------------------------------------------------------------------------------------------------- -pure function math_EulerAxisAngleToR(axis,omega) - - implicit none - real(pReal), dimension(3,3) :: math_EulerAxisAngleToR - real(pReal), dimension(3), intent(in) :: axis - real(pReal), intent(in) :: omega - - math_EulerAxisAngleToR = transpose(math_axisAngleToR(axis,omega)) ! convert to passive rotation - -end function math_EulerAxisAngleToR - - -!-------------------------------------------------------------------------------------------------- -!> @brief quaternion (w+ix+jy+kz) from Euler axis and angle (in radians) -!> @details quaternion is meant to represent a PASSIVE rotation -!> @details (see http://en.wikipedia.org/wiki/Euler_angles for definitions) -!> @details formula for active rotation taken from -!> @details http://en.wikipedia.org/wiki/Rotation_representation_%28mathematics%29#Rodrigues_parameters -!-------------------------------------------------------------------------------------------------- -pure function math_EulerAxisAngleToQ(axis,omega) - - implicit none - real(pReal), dimension(4) :: math_EulerAxisAngleToQ - real(pReal), dimension(3), intent(in) :: axis - real(pReal), intent(in) :: omega - - math_EulerAxisAngleToQ = math_qConj(math_axisAngleToQ(axis,omega)) ! convert to passive rotation - -end function math_EulerAxisAngleToQ - - -!-------------------------------------------------------------------------------------------------- -!> @brief quaternion (w+ix+jy+kz) from axis and angle (in radians) -!> @details quaternion is meant to represent an ACTIVE rotation -!> @details (see http://en.wikipedia.org/wiki/Euler_angles for definitions) -!> @details formula for active rotation taken from -!> @details http://en.wikipedia.org/wiki/Rotation_representation_%28mathematics%29#Rodrigues_parameters -!> @details equivalent to eu2qu (P=+1) from "D Rowenhorst et al. Consistent representations of and -!> @details conversions between 3D rotations, Model. Simul. Mater. Sci. Eng. 23-8 (2015)" -!-------------------------------------------------------------------------------------------------- -pure function math_axisAngleToQ(axis,omega) - - implicit none - real(pReal), dimension(4) :: math_axisAngleToQ - real(pReal), dimension(3), intent(in) :: axis - real(pReal), intent(in) :: omega - real(pReal), dimension(3) :: axisNrm - real(pReal) :: norm - - norm = norm2(axis) - wellDefined: if (norm > 1.0e-8_pReal) then - axisNrm = axis/norm ! normalize axis to be sure - math_axisAngleToQ = [cos(0.5_pReal*omega), sin(0.5_pReal*omega) * axisNrm(1:3)] - else wellDefined - math_axisAngleToQ = [1.0_pReal,0.0_pReal,0.0_pReal,0.0_pReal] - endif wellDefined - -end function math_axisAngleToQ - - -!-------------------------------------------------------------------------------------------------- -!> @brief orientation matrix from quaternion (w+ix+jy+kz) -!> @details taken from http://arxiv.org/pdf/math/0701759v1.pdf -!> @details see also http://en.wikipedia.org/wiki/Rotation_formalisms_in_three_dimensions -!-------------------------------------------------------------------------------------------------- -pure function math_qToR(q) - - implicit none - real(pReal), dimension(4), intent(in) :: q - real(pReal), dimension(3,3) :: math_qToR, T,S - integer :: i, j - - forall(i = 1:3, j = 1:3) T(i,j) = q(i+1) * q(j+1) - - S = reshape( [0.0_pReal, -q(4), q(3), & - q(4), 0.0_pReal, -q(2), & - -q(3), q(2), 0.0_pReal],[3,3]) ! notation is transposed - - math_qToR = (2.0_pReal * q(1)*q(1) - 1.0_pReal) * math_I3 & - + 2.0_pReal * T - 2.0_pReal * q(1) * S - -end function math_qToR - - -!-------------------------------------------------------------------------------------------------- -!> @brief 3-1-3 Euler angles (in radians) from quaternion (w+ix+jy+kz) -!> @details quaternion is meant to represent a PASSIVE rotation, -!> @details composed of INTRINSIC rotations around the axes of the -!> @details rotating reference frame -!> @details (see http://en.wikipedia.org/wiki/Euler_angles for definitions) -!-------------------------------------------------------------------------------------------------- -pure function math_qToEuler(qPassive) - - implicit none - real(pReal), dimension(4), intent(in) :: qPassive - real(pReal), dimension(4) :: q - real(pReal), dimension(3) :: math_qToEuler - - q = math_qConj(qPassive) ! convert to active rotation, since formulas are defined for active rotations - - math_qToEuler(2) = acos(1.0_pReal-2.0_pReal*(q(2)**2+q(3)**2)) - - if (abs(math_qToEuler(2)) < 1.0e-6_pReal) then - math_qToEuler(1) = sign(2.0_pReal*acos(math_clip(q(1),-1.0_pReal, 1.0_pReal)),q(4)) - math_qToEuler(3) = 0.0_pReal - else - math_qToEuler(1) = atan2(+q(1)*q(3)+q(2)*q(4), q(1)*q(2)-q(3)*q(4)) - math_qToEuler(3) = atan2(-q(1)*q(3)+q(2)*q(4), q(1)*q(2)+q(3)*q(4)) - endif - - math_qToEuler = merge(math_qToEuler + [2.0_pReal*PI, PI, 2.0_pReal*PI], & ! ensure correct range - math_qToEuler, math_qToEuler<0.0_pReal) - -end function math_qToEuler - - -!-------------------------------------------------------------------------------------------------- -!> @brief axis-angle (x, y, z, ang in radians) from quaternion (w+ix+jy+kz) -!> @details quaternion is meant to represent an ACTIVE rotation -!> @details (see http://en.wikipedia.org/wiki/Euler_angles for definitions) -!> @details formula for active rotation taken from -!> @details http://en.wikipedia.org/wiki/Rotation_representation_%28mathematics%29#Rodrigues_parameters -!-------------------------------------------------------------------------------------------------- -pure function math_qToAxisAngle(Q) - - implicit none - real(pReal), dimension(4), intent(in) :: Q - real(pReal) :: halfAngle, sinHalfAngle - real(pReal), dimension(4) :: math_qToAxisAngle - - halfAngle = acos(math_clip(Q(1),-1.0_pReal,1.0_pReal)) - sinHalfAngle = sin(halfAngle) - - smallRotation: if (sinHalfAngle <= 1.0e-4_pReal) then - math_qToAxisAngle = 0.0_pReal - else smallRotation - math_qToAxisAngle= [ Q(2:4)/sinHalfAngle, halfAngle*2.0_pReal] - endif smallRotation - -end function math_qToAxisAngle - - -!-------------------------------------------------------------------------------------------------- -!> @brief Euler axis-angle (x, y, z, ang in radians) from quaternion (w+ix+jy+kz) -!> @details quaternion is meant to represent a PASSIVE rotation -!> @details (see http://en.wikipedia.org/wiki/Euler_angles for definitions) -!-------------------------------------------------------------------------------------------------- -pure function math_qToEulerAxisAngle(qPassive) - - implicit none - real(pReal), dimension(4), intent(in) :: qPassive - real(pReal), dimension(4) :: q - real(pReal), dimension(4) :: math_qToEulerAxisAngle - - q = math_qConj(qPassive) ! convert to active rotation - math_qToEulerAxisAngle = math_qToAxisAngle(q) - -end function math_qToEulerAxisAngle - - !-------------------------------------------------------------------------------------------------- !> @brief Rodrigues vector (x, y, z) from unit quaternion (w+ix+jy+kz) !-------------------------------------------------------------------------------------------------- pure function math_qToRodrig(Q) - use, intrinsic :: & - IEEE_arithmetic - use prec, only: & - tol_math_check - implicit none - real(pReal), dimension(4), intent(in) :: Q - real(pReal), dimension(3) :: math_qToRodrig + real(pReal), dimension(4), intent(in) :: Q + real(pReal), dimension(3) :: math_qToRodrig - math_qToRodrig = merge(Q(2:4)/Q(1),IEEE_value(1.0_pReal,IEEE_quiet_NaN),abs(Q(1)) > tol_math_check)! NaN for 180 deg since Rodrig is unbound + math_qToRodrig = merge(Q(2:4)/Q(1),IEEE_value(1.0_pReal,IEEE_quiet_NaN),abs(Q(1)) > tol_math_check)! NaN for 180 deg since Rodrig is unbound end function math_qToRodrig -!-------------------------------------------------------------------------------------------------- -!> @brief misorientation angle between two sets of Euler angles -!-------------------------------------------------------------------------------------------------- -real(pReal) pure function math_EulerMisorientation(EulerA,EulerB) - - implicit none - real(pReal), dimension(3), intent(in) :: EulerA,EulerB - real(pReal) :: cosTheta - - cosTheta = (math_trace33(matmul(math_EulerToR(EulerB), & - transpose(math_EulerToR(EulerA)))) - 1.0_pReal) * 0.5_pReal - - math_EulerMisorientation = acos(math_clip(cosTheta,-1.0_pReal,1.0_pReal)) - -end function math_EulerMisorientation - - -!-------------------------------------------------------------------------------------------------- -!> @brief draw a random sample from Euler space -!-------------------------------------------------------------------------------------------------- -function math_sampleRandomOri() - - implicit none - real(pReal), dimension(3) :: math_sampleRandomOri, rnd - - rnd = halton([1,7,3]) - math_sampleRandomOri = [rnd(1)*2.0_pReal*PI, & - acos(2.0_pReal*rnd(2)-1.0_pReal), & - rnd(3)*2.0_pReal*PI] - -end function math_sampleRandomOri - - -!-------------------------------------------------------------------------------------------------- -!> @brief draw a sample from an Gaussian distribution around given orientation and Full Width -! at Half Maximum (FWHM) -!> @details: A uniform misorientation (limited to 2*FWHM) is sampled followed by convolution with -! a Gausian distribution -!-------------------------------------------------------------------------------------------------- -function math_sampleGaussOri(center,FWHM) - - implicit none - real(pReal), intent(in) :: FWHM - real(pReal), dimension(3), intent(in) :: center - real(pReal) :: angle - real(pReal), dimension(3) :: math_sampleGaussOri, axis - real(pReal), dimension(4) :: rnd - real(pReal), dimension(3,3) :: R - - if (FWHM < 0.1_pReal*INRAD) then - math_sampleGaussOri = center - else - GaussConvolution: do - rnd = halton([8,3,6,11]) - axis(1) = rnd(1)*2.0_pReal-1.0_pReal ! uniform on [-1,1] - axis(2:3) = [sqrt(1.0-axis(1)**2.0_pReal)*cos(rnd(2)*2.0*PI),& - sqrt(1.0-axis(1)**2.0_pReal)*sin(rnd(2)*2.0*PI)] ! random axis - angle = (rnd(3)-0.5_pReal)*4.0_pReal*FWHM ! rotation by [0, +-2 FWHM] - R = math_axisAngleToR(axis,angle) - angle = math_EulerMisorientation([0.0_pReal,0.0_pReal,0.0_pReal],math_RtoEuler(R)) - if (rnd(4) <= exp(-4.0_pReal*log(2.0_pReal)*(angle/FWHM)**2_pReal)) exit ! rejection sampling (Gaussian) - enddo GaussConvolution - math_sampleGaussOri = math_RtoEuler(matmul(R,math_EulerToR(center))) - endif - -end function math_sampleGaussOri - - -!-------------------------------------------------------------------------------------------------- -!> @brief draw a sample from an Gaussian distribution around given fiber texture and Full Width -! at Half Maximum (FWHM) -!------------------------------------------------------------------------------------------------- -function math_sampleFiberOri(alpha,beta,FWHM) - - implicit none - real(pReal), dimension(2), intent(in) :: alpha,beta - real(pReal), intent(in) :: FWHM - real(pReal), dimension(3) :: math_sampleFiberOri, & - fInC,& !< fiber axis in crystal coordinate system - fInS,& !< fiber axis in sample coordinate system - u - real(pReal), dimension(3) :: rnd - real(pReal), dimension(:),allocatable :: a !< 2D vector to tilt - integer, dimension(:),allocatable :: idx !< components of 2D vector - real(pReal), dimension(3,3) :: R !< Rotation matrix (composed of three components) - real(pReal):: angle,c - integer:: j,& !< index of smallest component - i - - allocate(a(0)) - allocate(idx(0)) - fInC = [sin(alpha(1))*cos(alpha(2)), sin(alpha(1))*sin(alpha(2)), cos(alpha(1))] - fInS = [sin(beta(1))*cos(beta(2)), sin(beta(1))*sin(beta(2)), cos(beta(1))] - - R = math_EulerAxisAngleToR(math_crossproduct(fInC,fInS),-acos(dot_product(fInC,fInS))) !< rotation to align fiber axis in crystal and sample system - - rnd = halton([7,10,3]) - R = matmul(R,math_EulerAxisAngleToR(fInS,rnd(1)*2.0_pReal*PI)) !< additional rotation (0..360deg) perpendicular to fiber axis - - if (FWHM > 0.1_pReal*INRAD) then - reducedTo2D: do i=1,3 - if (i /= minloc(abs(fInS),1)) then - a=[a,fInS(i)] - idx=[idx,i] - else - j = i - endif - enddo reducedTo2D - GaussConvolution: do - angle = (rnd(2)-0.5_pReal)*4.0_pReal*FWHM ! rotation by [0, +-2 FWHM] - ! solve cos(angle) = dot_product(fInS,u) under the assumption that their smallest component is the same - c = cos(angle)-fInS(j)**2 - u(idx(2)) = -(2.0_pReal*c*a(2) + sqrt(4*((c*a(2))**2-sum(a**2)*(c**2-a(1)**2*(1-fInS(j)**2)))))/& - (2*sum(a**2)) - u(idx(1)) = sqrt(1-u(idx(2))**2-fInS(j)**2) - u(j) = fInS(j) - - rejectionSampling: if (rnd(3) <= exp(-4.0_pReal*log(2.0_pReal)*(angle/FWHM)**2_pReal)) then - R = matmul(R,math_EulerAxisAngleToR(math_crossproduct(u,fInS),angle)) ! tilt around direction of smallest component - exit - endif rejectionSampling - rnd = halton([7,10,3]) - enddo GaussConvolution - endif - math_sampleFiberOri = math_RtoEuler(R) - -end function math_sampleFiberOri - - !-------------------------------------------------------------------------------------------------- !> @brief draw a random sample from Gauss variable !-------------------------------------------------------------------------------------------------- real(pReal) function math_sampleGaussVar(meanvalue, stddev, width) - use prec, only: & - tol_math_check - implicit none - real(pReal), intent(in) :: meanvalue, & ! meanvalue of gauss distribution - stddev ! standard deviation of gauss distribution - real(pReal), intent(in), optional :: width ! width of considered values as multiples of standard deviation - real(pReal), dimension(2) :: rnd ! random numbers - real(pReal) :: scatter, & ! normalized scatter around meanvalue - myWidth + real(pReal), intent(in) :: meanvalue, & ! meanvalue of gauss distribution + stddev ! standard deviation of gauss distribution + real(pReal), intent(in), optional :: width ! width of considered values as multiples of standard deviation + real(pReal), dimension(2) :: rnd ! random numbers + real(pReal) :: scatter, & ! normalized scatter around meanvalue + myWidth - if (abs(stddev) < tol_math_check) then - math_sampleGaussVar = meanvalue - else - myWidth = merge(width,3.0_pReal,present(width)) ! use +-3*sigma as default value for scatter if not given - - do - call random_number(rnd) - scatter = myWidth * (2.0_pReal * rnd(1) - 1.0_pReal) - if (rnd(2) <= exp(-0.5_pReal * scatter ** 2.0_pReal)) exit ! test if scattered value is drawn - enddo + if (abs(stddev) < tol_math_check) then + math_sampleGaussVar = meanvalue + else + myWidth = merge(width,3.0_pReal,present(width)) ! use +-3*sigma as default value for scatter if not given + + do + call random_number(rnd) + scatter = myWidth * (2.0_pReal * rnd(1) - 1.0_pReal) + if (rnd(2) <= exp(-0.5_pReal * scatter ** 2.0_pReal)) exit ! test if scattered value is drawn + enddo - math_sampleGaussVar = scatter * stddev - endif + math_sampleGaussVar = scatter * stddev + endif end function math_sampleGaussVar -!-------------------------------------------------------------------------------------------------- -!> @brief symmetrically equivalent Euler angles for given sample symmetry -!> @detail 1 (equivalent to != 2 and !=4):triclinic, 2:monoclinic, 4:orthotropic -!-------------------------------------------------------------------------------------------------- -pure function math_symmetricEulers(sym,Euler) - - implicit none - integer, intent(in) :: sym !< symmetry Class - real(pReal), dimension(3), intent(in) :: Euler - real(pReal), dimension(3,3) :: math_symmetricEulers - - math_symmetricEulers = transpose(reshape([PI+Euler(1), PI-Euler(1), 2.0_pReal*PI-Euler(1), & - Euler(2), PI-Euler(2), PI -Euler(2), & - Euler(3), PI+Euler(3), PI +Euler(3)],[3,3])) ! transpose is needed to have symbolic notation instead of column-major - - math_symmetricEulers = modulo(math_symmetricEulers,2.0_pReal*pi) - - select case (sym) - case (4) ! orthotropic: all done - - case (2) ! monoclinic: return only first - math_symmetricEulers(1:3,2:3) = 0.0_pReal - - case default ! triclinic: return blank - math_symmetricEulers = 0.0_pReal - end select - -end function math_symmetricEulers - - !-------------------------------------------------------------------------------------------------- !> @brief eigenvalues and eigenvectors of symmetric matrix m ! ToDo: has wrong oder of arguments !-------------------------------------------------------------------------------------------------- subroutine math_eigenValuesVectorsSym(m,values,vectors,error) - implicit none - real(pReal), dimension(:,:), intent(in) :: m - real(pReal), dimension(size(m,1)), intent(out) :: values - real(pReal), dimension(size(m,1),size(m,1)), intent(out) :: vectors - logical, intent(out) :: error - integer :: info - real(pReal), dimension((64+2)*size(m,1)) :: work ! block size of 64 taken from http://www.netlib.org/lapack/double/dsyev.f - external :: & - dsyev + real(pReal), dimension(:,:), intent(in) :: m + real(pReal), dimension(size(m,1)), intent(out) :: values + real(pReal), dimension(size(m,1),size(m,1)), intent(out) :: vectors + logical, intent(out) :: error + integer :: info + real(pReal), dimension((64+2)*size(m,1)) :: work ! block size of 64 taken from http://www.netlib.org/lapack/double/dsyev.f + external :: & + dsyev - vectors = m ! copy matrix to input (doubles as output) array - call dsyev('V','U',size(m,1),vectors,size(m,1),values,work,(64+2)*size(m,1),info) - error = (info == 0) + vectors = m ! copy matrix to input (doubles as output) array + call dsyev('V','U',size(m,1),vectors,size(m,1),values,work,(64+2)*size(m,1),info) + error = (info == 0) end subroutine math_eigenValuesVectorsSym @@ -1967,50 +1229,49 @@ end subroutine math_eigenValuesVectorsSym !-------------------------------------------------------------------------------------------------- subroutine math_eigenValuesVectorsSym33(m,values,vectors) - implicit none - real(pReal), dimension(3,3),intent(in) :: m - real(pReal), dimension(3), intent(out) :: values - real(pReal), dimension(3,3),intent(out) :: vectors - real(pReal) :: T, U, norm, threshold - logical :: error - - values = math_eigenvaluesSym33(m) - - vectors(1:3,2) = [ m(1, 2) * m(2, 3) - m(1, 3) * m(2, 2), & - m(1, 3) * m(1, 2) - m(2, 3) * m(1, 1), & - m(1, 2)**2] - - T = maxval(abs(values)) - U = max(T, T**2) - threshold = sqrt(5.68e-14_pReal * U**2) - -! Calculate first eigenvector by the formula v[0] = (m - lambda[0]).e1 x (m - lambda[0]).e2 - vectors(1:3,1) = [ vectors(1,2) + m(1, 3) * values(1), & - vectors(2,2) + m(2, 3) * values(1), & - (m(1,1) - values(1)) * (m(2,2) - values(1)) - vectors(3,2)] - norm = norm2(vectors(1:3, 1)) - - fallback1: if(norm < threshold) then - call math_eigenValuesVectorsSym(m,values,vectors,error) - return - endif fallback1 - - vectors(1:3,1) = vectors(1:3, 1) / norm + real(pReal), dimension(3,3),intent(in) :: m + real(pReal), dimension(3), intent(out) :: values + real(pReal), dimension(3,3),intent(out) :: vectors + real(pReal) :: T, U, norm, threshold + logical :: error -! Calculate second eigenvector by the formula v[1] = (m - lambda[1]).e1 x (m - lambda[1]).e2 - vectors(1:3,2) = [ vectors(1,2) + m(1, 3) * values(2), & - vectors(2,2) + m(2, 3) * values(2), & - (m(1,1) - values(2)) * (m(2,2) - values(2)) - vectors(3,2)] - norm = norm2(vectors(1:3, 2)) - - fallback2: if(norm < threshold) then - call math_eigenValuesVectorsSym(m,values,vectors,error) - return - endif fallback2 - vectors(1:3,2) = vectors(1:3, 2) / norm - -! Calculate third eigenvector according to v[2] = v[0] x v[1] - vectors(1:3,3) = math_crossproduct(vectors(1:3,1),vectors(1:3,2)) + values = math_eigenvaluesSym33(m) + + vectors(1:3,2) = [ m(1, 2) * m(2, 3) - m(1, 3) * m(2, 2), & + m(1, 3) * m(1, 2) - m(2, 3) * m(1, 1), & + m(1, 2)**2] + + T = maxval(abs(values)) + U = max(T, T**2) + threshold = sqrt(5.68e-14_pReal * U**2) + + ! Calculate first eigenvector by the formula v[0] = (m - lambda[0]).e1 x (m - lambda[0]).e2 + vectors(1:3,1) = [ vectors(1,2) + m(1, 3) * values(1), & + vectors(2,2) + m(2, 3) * values(1), & + (m(1,1) - values(1)) * (m(2,2) - values(1)) - vectors(3,2)] + norm = norm2(vectors(1:3, 1)) + + fallback1: if(norm < threshold) then + call math_eigenValuesVectorsSym(m,values,vectors,error) + return + endif fallback1 + + vectors(1:3,1) = vectors(1:3, 1) / norm + + ! Calculate second eigenvector by the formula v[1] = (m - lambda[1]).e1 x (m - lambda[1]).e2 + vectors(1:3,2) = [ vectors(1,2) + m(1, 3) * values(2), & + vectors(2,2) + m(2, 3) * values(2), & + (m(1,1) - values(2)) * (m(2,2) - values(2)) - vectors(3,2)] + norm = norm2(vectors(1:3, 2)) + + fallback2: if(norm < threshold) then + call math_eigenValuesVectorsSym(m,values,vectors,error) + return + endif fallback2 + vectors(1:3,2) = vectors(1:3, 2) / norm + + ! Calculate third eigenvector according to v[2] = v[0] x v[1] + vectors(1:3,3) = math_cross(vectors(1:3,1),vectors(1:3,2)) end subroutine math_eigenValuesVectorsSym33 @@ -2020,22 +1281,21 @@ end subroutine math_eigenValuesVectorsSym33 !-------------------------------------------------------------------------------------------------- function math_eigenvectorBasisSym(m) - implicit none - real(pReal), dimension(:,:), intent(in) :: m - real(pReal), dimension(size(m,1)) :: values - real(pReal), dimension(size(m,1),size(m,1)) :: vectors - real(pReal), dimension(size(m,1),size(m,1)) :: math_eigenvectorBasisSym - logical :: error - integer :: i - - math_eigenvectorBasisSym = 0.0_pReal - call math_eigenValuesVectorsSym(m,values,vectors,error) - if(error) return + real(pReal), dimension(:,:), intent(in) :: m + real(pReal), dimension(size(m,1)) :: values + real(pReal), dimension(size(m,1),size(m,1)) :: vectors + real(pReal), dimension(size(m,1),size(m,1)) :: math_eigenvectorBasisSym + logical :: error + integer :: i - do i=1, size(m,1) - math_eigenvectorBasisSym = math_eigenvectorBasisSym & - + sqrt(values(i)) * math_outer(vectors(:,i),vectors(:,i)) - enddo + math_eigenvectorBasisSym = 0.0_pReal + call math_eigenValuesVectorsSym(m,values,vectors,error) + if(error) return + + do i=1, size(m,1) + math_eigenvectorBasisSym = math_eigenvectorBasisSym & + + sqrt(values(i)) * math_outer(vectors(:,i),vectors(:,i)) + enddo end function math_eigenvectorBasisSym @@ -2045,58 +1305,57 @@ end function math_eigenvectorBasisSym !-------------------------------------------------------------------------------------------------- pure function math_eigenvectorBasisSym33(m) - implicit none - real(pReal), dimension(3,3) :: math_eigenvectorBasisSym33 - real(pReal), dimension(3) :: invariants, values - real(pReal), dimension(3,3), intent(in) :: m - real(pReal) :: P, Q, rho, phi - real(pReal), parameter :: TOL=1.e-14_pReal - real(pReal), dimension(3,3,3) :: N, EB - - invariants = math_invariantsSym33(m) - EB = 0.0_pReal - - P = invariants(2)-invariants(1)**2.0_pReal/3.0_pReal - Q = -2.0_pReal/27.0_pReal*invariants(1)**3.0_pReal+product(invariants(1:2))/3.0_pReal-invariants(3) - - threeSimilarEigenvalues: if(all(abs([P,Q]) < TOL)) then - values = invariants(1)/3.0_pReal -! this is not really correct, but at least the basis is correct - EB(1,1,1)=1.0_pReal - EB(2,2,2)=1.0_pReal - EB(3,3,3)=1.0_pReal - else threeSimilarEigenvalues - rho=sqrt(-3.0_pReal*P**3.0_pReal)/9.0_pReal - phi=acos(math_clip(-Q/rho*0.5_pReal,-1.0_pReal,1.0_pReal)) - values = 2.0_pReal*rho**(1.0_pReal/3.0_pReal)* & - [cos(phi/3.0_pReal), & - cos((phi+2.0_pReal*PI)/3.0_pReal), & - cos((phi+4.0_pReal*PI)/3.0_pReal) & - ] + invariants(1)/3.0_pReal - N(1:3,1:3,1) = m-values(1)*math_I3 - N(1:3,1:3,2) = m-values(2)*math_I3 - N(1:3,1:3,3) = m-values(3)*math_I3 - twoSimilarEigenvalues: if(abs(values(1)-values(2)) < TOL) then - EB(1:3,1:3,3)=matmul(N(1:3,1:3,1),N(1:3,1:3,2))/ & - ((values(3)-values(1))*(values(3)-values(2))) - EB(1:3,1:3,1)=math_I3-EB(1:3,1:3,3) - elseif(abs(values(2)-values(3)) < TOL) then twoSimilarEigenvalues - EB(1:3,1:3,1)=matmul(N(1:3,1:3,2),N(1:3,1:3,3))/ & - ((values(1)-values(2))*(values(1)-values(3))) - EB(1:3,1:3,2)=math_I3-EB(1:3,1:3,1) - elseif(abs(values(3)-values(1)) < TOL) then twoSimilarEigenvalues - EB(1:3,1:3,2)=matmul(N(1:3,1:3,1),N(1:3,1:3,3))/ & - ((values(2)-values(1))*(values(2)-values(3))) - EB(1:3,1:3,1)=math_I3-EB(1:3,1:3,2) - else twoSimilarEigenvalues - EB(1:3,1:3,1)=matmul(N(1:3,1:3,2),N(1:3,1:3,3))/ & - ((values(1)-values(2))*(values(1)-values(3))) - EB(1:3,1:3,2)=matmul(N(1:3,1:3,1),N(1:3,1:3,3))/ & - ((values(2)-values(1))*(values(2)-values(3))) - EB(1:3,1:3,3)=matmul(N(1:3,1:3,1),N(1:3,1:3,2))/ & - ((values(3)-values(1))*(values(3)-values(2))) - endif twoSimilarEigenvalues - endif threeSimilarEigenvalues + real(pReal), dimension(3,3) :: math_eigenvectorBasisSym33 + real(pReal), dimension(3) :: invariants, values + real(pReal), dimension(3,3), intent(in) :: m + real(pReal) :: P, Q, rho, phi + real(pReal), parameter :: TOL=1.e-14_pReal + real(pReal), dimension(3,3,3) :: N, EB + + invariants = math_invariantsSym33(m) + EB = 0.0_pReal + + P = invariants(2)-invariants(1)**2.0_pReal/3.0_pReal + Q = -2.0_pReal/27.0_pReal*invariants(1)**3.0_pReal+product(invariants(1:2))/3.0_pReal-invariants(3) + + threeSimilarEigenvalues: if(all(abs([P,Q]) < TOL)) then + values = invariants(1)/3.0_pReal + ! this is not really correct, but at least the basis is correct + EB(1,1,1)=1.0_pReal + EB(2,2,2)=1.0_pReal + EB(3,3,3)=1.0_pReal + else threeSimilarEigenvalues + rho=sqrt(-3.0_pReal*P**3.0_pReal)/9.0_pReal + phi=acos(math_clip(-Q/rho*0.5_pReal,-1.0_pReal,1.0_pReal)) + values = 2.0_pReal*rho**(1.0_pReal/3.0_pReal)* & + [cos(phi/3.0_pReal), & + cos((phi+2.0_pReal*PI)/3.0_pReal), & + cos((phi+4.0_pReal*PI)/3.0_pReal) & + ] + invariants(1)/3.0_pReal + N(1:3,1:3,1) = m-values(1)*math_I3 + N(1:3,1:3,2) = m-values(2)*math_I3 + N(1:3,1:3,3) = m-values(3)*math_I3 + twoSimilarEigenvalues: if(abs(values(1)-values(2)) < TOL) then + EB(1:3,1:3,3)=matmul(N(1:3,1:3,1),N(1:3,1:3,2))/ & + ((values(3)-values(1))*(values(3)-values(2))) + EB(1:3,1:3,1)=math_I3-EB(1:3,1:3,3) + elseif(abs(values(2)-values(3)) < TOL) then twoSimilarEigenvalues + EB(1:3,1:3,1)=matmul(N(1:3,1:3,2),N(1:3,1:3,3))/ & + ((values(1)-values(2))*(values(1)-values(3))) + EB(1:3,1:3,2)=math_I3-EB(1:3,1:3,1) + elseif(abs(values(3)-values(1)) < TOL) then twoSimilarEigenvalues + EB(1:3,1:3,2)=matmul(N(1:3,1:3,1),N(1:3,1:3,3))/ & + ((values(2)-values(1))*(values(2)-values(3))) + EB(1:3,1:3,1)=math_I3-EB(1:3,1:3,2) + else twoSimilarEigenvalues + EB(1:3,1:3,1)=matmul(N(1:3,1:3,2),N(1:3,1:3,3))/ & + ((values(1)-values(2))*(values(1)-values(3))) + EB(1:3,1:3,2)=matmul(N(1:3,1:3,1),N(1:3,1:3,3))/ & + ((values(2)-values(1))*(values(2)-values(3))) + EB(1:3,1:3,3)=matmul(N(1:3,1:3,1),N(1:3,1:3,2))/ & + ((values(3)-values(1))*(values(3)-values(2))) + endif twoSimilarEigenvalues + endif threeSimilarEigenvalues math_eigenvectorBasisSym33 = sqrt(values(1)) * EB(1:3,1:3,1) & + sqrt(values(2)) * EB(1:3,1:3,2) & @@ -2110,62 +1369,61 @@ end function math_eigenvectorBasisSym33 !-------------------------------------------------------------------------------------------------- pure function math_eigenvectorBasisSym33_log(m) - implicit none - real(pReal), dimension(3,3) :: math_eigenvectorBasisSym33_log - real(pReal), dimension(3) :: invariants, values - real(pReal), dimension(3,3), intent(in) :: m - real(pReal) :: P, Q, rho, phi - real(pReal), parameter :: TOL=1.e-14_pReal - real(pReal), dimension(3,3,3) :: N, EB + real(pReal), dimension(3,3) :: math_eigenvectorBasisSym33_log + real(pReal), dimension(3) :: invariants, values + real(pReal), dimension(3,3), intent(in) :: m + real(pReal) :: P, Q, rho, phi + real(pReal), parameter :: TOL=1.e-14_pReal + real(pReal), dimension(3,3,3) :: N, EB - invariants = math_invariantsSym33(m) - EB = 0.0_pReal + invariants = math_invariantsSym33(m) + EB = 0.0_pReal - P = invariants(2)-invariants(1)**2.0_pReal/3.0_pReal - Q = -2.0_pReal/27.0_pReal*invariants(1)**3.0_pReal+product(invariants(1:2))/3.0_pReal-invariants(3) + P = invariants(2)-invariants(1)**2.0_pReal/3.0_pReal + Q = -2.0_pReal/27.0_pReal*invariants(1)**3.0_pReal+product(invariants(1:2))/3.0_pReal-invariants(3) - threeSimilarEigenvalues: if(all(abs([P,Q]) < TOL)) then - values = invariants(1)/3.0_pReal -! this is not really correct, but at least the basis is correct - EB(1,1,1)=1.0_pReal - EB(2,2,2)=1.0_pReal - EB(3,3,3)=1.0_pReal - else threeSimilarEigenvalues - rho=sqrt(-3.0_pReal*P**3.0_pReal)/9.0_pReal - phi=acos(math_clip(-Q/rho*0.5_pReal,-1.0_pReal,1.0_pReal)) - values = 2.0_pReal*rho**(1.0_pReal/3.0_pReal)* & - [cos(phi/3.0_pReal), & - cos((phi+2.0_pReal*PI)/3.0_pReal), & - cos((phi+4.0_pReal*PI)/3.0_pReal) & - ] + invariants(1)/3.0_pReal - N(1:3,1:3,1) = m-values(1)*math_I3 - N(1:3,1:3,2) = m-values(2)*math_I3 - N(1:3,1:3,3) = m-values(3)*math_I3 - twoSimilarEigenvalues: if(abs(values(1)-values(2)) < TOL) then - EB(1:3,1:3,3)=matmul(N(1:3,1:3,1),N(1:3,1:3,2))/ & - ((values(3)-values(1))*(values(3)-values(2))) - EB(1:3,1:3,1)=math_I3-EB(1:3,1:3,3) - elseif(abs(values(2)-values(3)) < TOL) then twoSimilarEigenvalues - EB(1:3,1:3,1)=matmul(N(1:3,1:3,2),N(1:3,1:3,3))/ & - ((values(1)-values(2))*(values(1)-values(3))) - EB(1:3,1:3,2)=math_I3-EB(1:3,1:3,1) - elseif(abs(values(3)-values(1)) < TOL) then twoSimilarEigenvalues - EB(1:3,1:3,2)=matmul(N(1:3,1:3,1),N(1:3,1:3,3))/ & - ((values(2)-values(1))*(values(2)-values(3))) - EB(1:3,1:3,1)=math_I3-EB(1:3,1:3,2) - else twoSimilarEigenvalues - EB(1:3,1:3,1)=matmul(N(1:3,1:3,2),N(1:3,1:3,3))/ & - ((values(1)-values(2))*(values(1)-values(3))) - EB(1:3,1:3,2)=matmul(N(1:3,1:3,1),N(1:3,1:3,3))/ & - ((values(2)-values(1))*(values(2)-values(3))) - EB(1:3,1:3,3)=matmul(N(1:3,1:3,1),N(1:3,1:3,2))/ & - ((values(3)-values(1))*(values(3)-values(2))) - endif twoSimilarEigenvalues - endif threeSimilarEigenvalues + threeSimilarEigenvalues: if(all(abs([P,Q]) < TOL)) then + values = invariants(1)/3.0_pReal + ! this is not really correct, but at least the basis is correct + EB(1,1,1)=1.0_pReal + EB(2,2,2)=1.0_pReal + EB(3,3,3)=1.0_pReal + else threeSimilarEigenvalues + rho=sqrt(-3.0_pReal*P**3.0_pReal)/9.0_pReal + phi=acos(math_clip(-Q/rho*0.5_pReal,-1.0_pReal,1.0_pReal)) + values = 2.0_pReal*rho**(1.0_pReal/3.0_pReal)* & + [cos(phi/3.0_pReal), & + cos((phi+2.0_pReal*PI)/3.0_pReal), & + cos((phi+4.0_pReal*PI)/3.0_pReal) & + ] + invariants(1)/3.0_pReal + N(1:3,1:3,1) = m-values(1)*math_I3 + N(1:3,1:3,2) = m-values(2)*math_I3 + N(1:3,1:3,3) = m-values(3)*math_I3 + twoSimilarEigenvalues: if(abs(values(1)-values(2)) < TOL) then + EB(1:3,1:3,3)=matmul(N(1:3,1:3,1),N(1:3,1:3,2))/ & + ((values(3)-values(1))*(values(3)-values(2))) + EB(1:3,1:3,1)=math_I3-EB(1:3,1:3,3) + elseif(abs(values(2)-values(3)) < TOL) then twoSimilarEigenvalues + EB(1:3,1:3,1)=matmul(N(1:3,1:3,2),N(1:3,1:3,3))/ & + ((values(1)-values(2))*(values(1)-values(3))) + EB(1:3,1:3,2)=math_I3-EB(1:3,1:3,1) + elseif(abs(values(3)-values(1)) < TOL) then twoSimilarEigenvalues + EB(1:3,1:3,2)=matmul(N(1:3,1:3,1),N(1:3,1:3,3))/ & + ((values(2)-values(1))*(values(2)-values(3))) + EB(1:3,1:3,1)=math_I3-EB(1:3,1:3,2) + else twoSimilarEigenvalues + EB(1:3,1:3,1)=matmul(N(1:3,1:3,2),N(1:3,1:3,3))/ & + ((values(1)-values(2))*(values(1)-values(3))) + EB(1:3,1:3,2)=matmul(N(1:3,1:3,1),N(1:3,1:3,3))/ & + ((values(2)-values(1))*(values(2)-values(3))) + EB(1:3,1:3,3)=matmul(N(1:3,1:3,1),N(1:3,1:3,2))/ & + ((values(3)-values(1))*(values(3)-values(2))) + endif twoSimilarEigenvalues + endif threeSimilarEigenvalues - math_eigenvectorBasisSym33_log = log(sqrt(values(1))) * EB(1:3,1:3,1) & - + log(sqrt(values(2))) * EB(1:3,1:3,2) & - + log(sqrt(values(3))) * EB(1:3,1:3,3) + math_eigenvectorBasisSym33_log = log(sqrt(values(1))) * EB(1:3,1:3,1) & + + log(sqrt(values(2))) * EB(1:3,1:3,2) & + + log(sqrt(values(3))) * EB(1:3,1:3,3) end function math_eigenvectorBasisSym33_log @@ -2174,25 +1432,22 @@ end function math_eigenvectorBasisSym33_log !> @brief rotational part from polar decomposition of 33 tensor m !-------------------------------------------------------------------------------------------------- function math_rotationalPart33(m) - use prec, only: & - dEq0 - use IO, only: & - IO_warning + use IO, only: & + IO_warning - implicit none - real(pReal), intent(in), dimension(3,3) :: m - real(pReal), dimension(3,3) :: math_rotationalPart33 - real(pReal), dimension(3,3) :: U , Uinv + real(pReal), intent(in), dimension(3,3) :: m + real(pReal), dimension(3,3) :: math_rotationalPart33 + real(pReal), dimension(3,3) :: U , Uinv - U = math_eigenvectorBasisSym33(matmul(transpose(m),m)) - Uinv = math_inv33(U) + U = math_eigenvectorBasisSym33(matmul(transpose(m),m)) + Uinv = math_inv33(U) - inversionFailed: if (all(dEq0(Uinv))) then - math_rotationalPart33 = math_I3 - call IO_warning(650) - else inversionFailed - math_rotationalPart33 = matmul(m,Uinv) - endif inversionFailed + inversionFailed: if (all(dEq0(Uinv))) then + math_rotationalPart33 = math_I3 + call IO_warning(650) + else inversionFailed + math_rotationalPart33 = matmul(m,Uinv) + endif inversionFailed end function math_rotationalPart33 @@ -2202,21 +1457,18 @@ end function math_rotationalPart33 ! will return NaN on error !-------------------------------------------------------------------------------------------------- function math_eigenvaluesSym(m) - use, intrinsic :: & - IEEE_arithmetic - implicit none - real(pReal), dimension(:,:), intent(in) :: m - real(pReal), dimension(size(m,1)) :: math_eigenvaluesSym - real(pReal), dimension(size(m,1),size(m,1)) :: vectors - integer :: info - real(pReal), dimension((64+2)*size(m,1)) :: work ! block size of 64 taken from http://www.netlib.org/lapack/double/dsyev.f - external :: & - dsyev + real(pReal), dimension(:,:), intent(in) :: m + real(pReal), dimension(size(m,1)) :: math_eigenvaluesSym + real(pReal), dimension(size(m,1),size(m,1)) :: vectors + integer :: info + real(pReal), dimension((64+2)*size(m,1)) :: work ! block size of 64 taken from http://www.netlib.org/lapack/double/dsyev.f + external :: & + dsyev - vectors = m ! copy matrix to input (doubles as output) array - call dsyev('N','U',size(m,1),vectors,size(m,1),math_eigenvaluesSym,work,(64+2)*size(m,1),info) - if (info /= 0) math_eigenvaluesSym = IEEE_value(1.0_pReal,IEEE_quiet_NaN) + vectors = m ! copy matrix to input (doubles as output) array + call dsyev('N','U',size(m,1),vectors,size(m,1),math_eigenvaluesSym,work,(64+2)*size(m,1),info) + if (info /= 0) math_eigenvaluesSym = IEEE_value(1.0_pReal,IEEE_quiet_NaN) end function math_eigenvaluesSym @@ -2230,28 +1482,27 @@ end function math_eigenvaluesSym !-------------------------------------------------------------------------------------------------- function math_eigenvaluesSym33(m) - implicit none - real(pReal), intent(in), dimension(3,3) :: m - real(pReal), dimension(3) :: math_eigenvaluesSym33,invariants - real(pReal) :: P, Q, rho, phi - real(pReal), parameter :: TOL=1.e-14_pReal + real(pReal), intent(in), dimension(3,3) :: m + real(pReal), dimension(3) :: math_eigenvaluesSym33,invariants + real(pReal) :: P, Q, rho, phi + real(pReal), parameter :: TOL=1.e-14_pReal - invariants = math_invariantsSym33(m) ! invariants are coefficients in characteristic polynomial apart for the sign of c0 and c2 in http://arxiv.org/abs/physics/0610206 + invariants = math_invariantsSym33(m) ! invariants are coefficients in characteristic polynomial apart for the sign of c0 and c2 in http://arxiv.org/abs/physics/0610206 - P = invariants(2)-invariants(1)**2.0_pReal/3.0_pReal ! different from http://arxiv.org/abs/physics/0610206 (this formulation was in DAMASK) - Q = -2.0_pReal/27.0_pReal*invariants(1)**3.0_pReal+product(invariants(1:2))/3.0_pReal-invariants(3)! different from http://arxiv.org/abs/physics/0610206 (this formulation was in DAMASK) + P = invariants(2)-invariants(1)**2.0_pReal/3.0_pReal ! different from http://arxiv.org/abs/physics/0610206 (this formulation was in DAMASK) + Q = -2.0_pReal/27.0_pReal*invariants(1)**3.0_pReal+product(invariants(1:2))/3.0_pReal-invariants(3)! different from http://arxiv.org/abs/physics/0610206 (this formulation was in DAMASK) - if(all(abs([P,Q]) < TOL)) then - math_eigenvaluesSym33 = math_eigenvaluesSym(m) - else - rho=sqrt(-3.0_pReal*P**3.0_pReal)/9.0_pReal - phi=acos(math_clip(-Q/rho*0.5_pReal,-1.0_pReal,1.0_pReal)) - math_eigenvaluesSym33 = 2.0_pReal*rho**(1.0_pReal/3.0_pReal)* & - [cos(phi/3.0_pReal), & - cos((phi+2.0_pReal*PI)/3.0_pReal), & - cos((phi+4.0_pReal*PI)/3.0_pReal) & - ] + invariants(1)/3.0_pReal - endif + if(all(abs([P,Q]) < TOL)) then + math_eigenvaluesSym33 = math_eigenvaluesSym(m) + else + rho=sqrt(-3.0_pReal*P**3.0_pReal)/9.0_pReal + phi=acos(math_clip(-Q/rho*0.5_pReal,-1.0_pReal,1.0_pReal)) + math_eigenvaluesSym33 = 2.0_pReal*rho**(1.0_pReal/3.0_pReal)* & + [cos(phi/3.0_pReal), & + cos((phi+2.0_pReal*PI)/3.0_pReal), & + cos((phi+4.0_pReal*PI)/3.0_pReal) & + ] + invariants(1)/3.0_pReal + endif end function math_eigenvaluesSym33 @@ -2261,252 +1512,26 @@ end function math_eigenvaluesSym33 !-------------------------------------------------------------------------------------------------- pure function math_invariantsSym33(m) - implicit none - real(pReal), dimension(3,3), intent(in) :: m - real(pReal), dimension(3) :: math_invariantsSym33 + real(pReal), dimension(3,3), intent(in) :: m + real(pReal), dimension(3) :: math_invariantsSym33 - math_invariantsSym33(1) = math_trace33(m) - math_invariantsSym33(2) = m(1,1)*m(2,2) + m(1,1)*m(3,3) + m(2,2)*m(3,3) & - -(m(1,2)**2 + m(1,3)**2 + m(2,3)**2) - math_invariantsSym33(3) = math_detSym33(m) + math_invariantsSym33(1) = math_trace33(m) + math_invariantsSym33(2) = m(1,1)*m(2,2) + m(1,1)*m(3,3) + m(2,2)*m(3,3) & + -(m(1,2)**2 + m(1,3)**2 + m(2,3)**2) + math_invariantsSym33(3) = math_detSym33(m) end function math_invariantsSym33 -!------------------------------------------------------------------------------------------------- -!> @brief computes an element of a Halton sequence. -!> @author John Burkardt -!> @author Martin Diehl -!> @details Incrementally increasing elements of the Halton sequence for given bases (> 0) -!> @details Reference: -!> @details J.H. Halton: On the efficiency of certain quasi-random sequences of points in evaluating -!> @details multi-dimensional integrals, Numerische Mathematik, Volume 2, pages 84-90, 1960. -!> @details Reference for prime numbers: -!> @details Milton Abramowitz and Irene Stegun: Handbook of Mathematical Functions, -!> @details US Department of Commerce, 1964, pages 870-873. -!> @details Daniel Zwillinger: CRC Standard Mathematical Tables and Formulae, -!> @details 30th Edition, CRC Press, 1996, pages 95-98. -!------------------------------------------------------------------------------------------------- -function halton(bases) - - implicit none - integer, intent(in), dimension(:):: & - bases !< bases (prime number ID) - real(pReal), dimension(size(bases)) :: & - halton - integer, save :: & - current = 1 - real(pReal), dimension(size(bases)) :: & - base_inv - integer, dimension(size(bases)) :: & - base, & - t - integer, dimension(0:1600), parameter :: & - prime = [& - 1, & - 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, & - 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, & - 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, & - 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, & - 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, & - 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, & - 283, 293, 307, 311, 313, 317, 331, 337, 347, 349, & - 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, & - 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, & - 467, 479, 487, 491, 499, 503, 509, 521, 523, 541, & - ! 101:200 - 547, 557, 563, 569, 571, 577, 587, 593, 599, 601, & - 607, 613, 617, 619, 631, 641, 643, 647, 653, 659, & - 661, 673, 677, 683, 691, 701, 709, 719, 727, 733, & - 739, 743, 751, 757, 761, 769, 773, 787, 797, 809, & - 811, 821, 823, 827, 829, 839, 853, 857, 859, 863, & - 877, 881, 883, 887, 907, 911, 919, 929, 937, 941, & - 947, 953, 967, 971, 977, 983, 991, 997, 1009, 1013, & - 1019, 1021, 1031, 1033, 1039, 1049, 1051, 1061, 1063, 1069, & - 1087, 1091, 1093, 1097, 1103, 1109, 1117, 1123, 1129, 1151, & - 1153, 1163, 1171, 1181, 1187, 1193, 1201, 1213, 1217, 1223, & - ! 201:300 - 1229, 1231, 1237, 1249, 1259, 1277, 1279, 1283, 1289, 1291, & - 1297, 1301, 1303, 1307, 1319, 1321, 1327, 1361, 1367, 1373, & - 1381, 1399, 1409, 1423, 1427, 1429, 1433, 1439, 1447, 1451, & - 1453, 1459, 1471, 1481, 1483, 1487, 1489, 1493, 1499, 1511, & - 1523, 1531, 1543, 1549, 1553, 1559, 1567, 1571, 1579, 1583, & - 1597, 1601, 1607, 1609, 1613, 1619, 1621, 1627, 1637, 1657, & - 1663, 1667, 1669, 1693, 1697, 1699, 1709, 1721, 1723, 1733, & - 1741, 1747, 1753, 1759, 1777, 1783, 1787, 1789, 1801, 1811, & - 1823, 1831, 1847, 1861, 1867, 1871, 1873, 1877, 1879, 1889, & - 1901, 1907, 1913, 1931, 1933, 1949, 1951, 1973, 1979, 1987, & - ! 301:400 - 1993, 1997, 1999, 2003, 2011, 2017, 2027, 2029, 2039, 2053, & - 2063, 2069, 2081, 2083, 2087, 2089, 2099, 2111, 2113, 2129, & - 2131, 2137, 2141, 2143, 2153, 2161, 2179, 2203, 2207, 2213, & - 2221, 2237, 2239, 2243, 2251, 2267, 2269, 2273, 2281, 2287, & - 2293, 2297, 2309, 2311, 2333, 2339, 2341, 2347, 2351, 2357, & - 2371, 2377, 2381, 2383, 2389, 2393, 2399, 2411, 2417, 2423, & - 2437, 2441, 2447, 2459, 2467, 2473, 2477, 2503, 2521, 2531, & - 2539, 2543, 2549, 2551, 2557, 2579, 2591, 2593, 2609, 2617, & - 2621, 2633, 2647, 2657, 2659, 2663, 2671, 2677, 2683, 2687, & - 2689, 2693, 2699, 2707, 2711, 2713, 2719, 2729, 2731, 2741, & - ! 401:500 - 2749, 2753, 2767, 2777, 2789, 2791, 2797, 2801, 2803, 2819, & - 2833, 2837, 2843, 2851, 2857, 2861, 2879, 2887, 2897, 2903, & - 2909, 2917, 2927, 2939, 2953, 2957, 2963, 2969, 2971, 2999, & - 3001, 3011, 3019, 3023, 3037, 3041, 3049, 3061, 3067, 3079, & - 3083, 3089, 3109, 3119, 3121, 3137, 3163, 3167, 3169, 3181, & - 3187, 3191, 3203, 3209, 3217, 3221, 3229, 3251, 3253, 3257, & - 3259, 3271, 3299, 3301, 3307, 3313, 3319, 3323, 3329, 3331, & - 3343, 3347, 3359, 3361, 3371, 3373, 3389, 3391, 3407, 3413, & - 3433, 3449, 3457, 3461, 3463, 3467, 3469, 3491, 3499, 3511, & - 3517, 3527, 3529, 3533, 3539, 3541, 3547, 3557, 3559, 3571, & - ! 501:600 - 3581, 3583, 3593, 3607, 3613, 3617, 3623, 3631, 3637, 3643, & - 3659, 3671, 3673, 3677, 3691, 3697, 3701, 3709, 3719, 3727, & - 3733, 3739, 3761, 3767, 3769, 3779, 3793, 3797, 3803, 3821, & - 3823, 3833, 3847, 3851, 3853, 3863, 3877, 3881, 3889, 3907, & - 3911, 3917, 3919, 3923, 3929, 3931, 3943, 3947, 3967, 3989, & - 4001, 4003, 4007, 4013, 4019, 4021, 4027, 4049, 4051, 4057, & - 4073, 4079, 4091, 4093, 4099, 4111, 4127, 4129, 4133, 4139, & - 4153, 4157, 4159, 4177, 4201, 4211, 4217, 4219, 4229, 4231, & - 4241, 4243, 4253, 4259, 4261, 4271, 4273, 4283, 4289, 4297, & - 4327, 4337, 4339, 4349, 4357, 4363, 4373, 4391, 4397, 4409, & - ! 601:700 - 4421, 4423, 4441, 4447, 4451, 4457, 4463, 4481, 4483, 4493, & - 4507, 4513, 4517, 4519, 4523, 4547, 4549, 4561, 4567, 4583, & - 4591, 4597, 4603, 4621, 4637, 4639, 4643, 4649, 4651, 4657, & - 4663, 4673, 4679, 4691, 4703, 4721, 4723, 4729, 4733, 4751, & - 4759, 4783, 4787, 4789, 4793, 4799, 4801, 4813, 4817, 4831, & - 4861, 4871, 4877, 4889, 4903, 4909, 4919, 4931, 4933, 4937, & - 4943, 4951, 4957, 4967, 4969, 4973, 4987, 4993, 4999, 5003, & - 5009, 5011, 5021, 5023, 5039, 5051, 5059, 5077, 5081, 5087, & - 5099, 5101, 5107, 5113, 5119, 5147, 5153, 5167, 5171, 5179, & - 5189, 5197, 5209, 5227, 5231, 5233, 5237, 5261, 5273, 5279, & - ! 701:800 - 5281, 5297, 5303, 5309, 5323, 5333, 5347, 5351, 5381, 5387, & - 5393, 5399, 5407, 5413, 5417, 5419, 5431, 5437, 5441, 5443, & - 5449, 5471, 5477, 5479, 5483, 5501, 5503, 5507, 5519, 5521, & - 5527, 5531, 5557, 5563, 5569, 5573, 5581, 5591, 5623, 5639, & - 5641, 5647, 5651, 5653, 5657, 5659, 5669, 5683, 5689, 5693, & - 5701, 5711, 5717, 5737, 5741, 5743, 5749, 5779, 5783, 5791, & - 5801, 5807, 5813, 5821, 5827, 5839, 5843, 5849, 5851, 5857, & - 5861, 5867, 5869, 5879, 5881, 5897, 5903, 5923, 5927, 5939, & - 5953, 5981, 5987, 6007, 6011, 6029, 6037, 6043, 6047, 6053, & - 6067, 6073, 6079, 6089, 6091, 6101, 6113, 6121, 6131, 6133, & - ! 801:900 - 6143, 6151, 6163, 6173, 6197, 6199, 6203, 6211, 6217, 6221, & - 6229, 6247, 6257, 6263, 6269, 6271, 6277, 6287, 6299, 6301, & - 6311, 6317, 6323, 6329, 6337, 6343, 6353, 6359, 6361, 6367, & - 6373, 6379, 6389, 6397, 6421, 6427, 6449, 6451, 6469, 6473, & - 6481, 6491, 6521, 6529, 6547, 6551, 6553, 6563, 6569, 6571, & - 6577, 6581, 6599, 6607, 6619, 6637, 6653, 6659, 6661, 6673, & - 6679, 6689, 6691, 6701, 6703, 6709, 6719, 6733, 6737, 6761, & - 6763, 6779, 6781, 6791, 6793, 6803, 6823, 6827, 6829, 6833, & - 6841, 6857, 6863, 6869, 6871, 6883, 6899, 6907, 6911, 6917, & - 6947, 6949, 6959, 6961, 6967, 6971, 6977, 6983, 6991, 6997, & - ! 901:1000 - 7001, 7013, 7019, 7027, 7039, 7043, 7057, 7069, 7079, 7103, & - 7109, 7121, 7127, 7129, 7151, 7159, 7177, 7187, 7193, 7207, & - 7211, 7213, 7219, 7229, 7237, 7243, 7247, 7253, 7283, 7297, & - 7307, 7309, 7321, 7331, 7333, 7349, 7351, 7369, 7393, 7411, & - 7417, 7433, 7451, 7457, 7459, 7477, 7481, 7487, 7489, 7499, & - 7507, 7517, 7523, 7529, 7537, 7541, 7547, 7549, 7559, 7561, & - 7573, 7577, 7583, 7589, 7591, 7603, 7607, 7621, 7639, 7643, & - 7649, 7669, 7673, 7681, 7687, 7691, 7699, 7703, 7717, 7723, & - 7727, 7741, 7753, 7757, 7759, 7789, 7793, 7817, 7823, 7829, & - 7841, 7853, 7867, 7873, 7877, 7879, 7883, 7901, 7907, 7919, & - ! 1001:1100 - 7927, 7933, 7937, 7949, 7951, 7963, 7993, 8009, 8011, 8017, & - 8039, 8053, 8059, 8069, 8081, 8087, 8089, 8093, 8101, 8111, & - 8117, 8123, 8147, 8161, 8167, 8171, 8179, 8191, 8209, 8219, & - 8221, 8231, 8233, 8237, 8243, 8263, 8269, 8273, 8287, 8291, & - 8293, 8297, 8311, 8317, 8329, 8353, 8363, 8369, 8377, 8387, & - 8389, 8419, 8423, 8429, 8431, 8443, 8447, 8461, 8467, 8501, & - 8513, 8521, 8527, 8537, 8539, 8543, 8563, 8573, 8581, 8597, & - 8599, 8609, 8623, 8627, 8629, 8641, 8647, 8663, 8669, 8677, & - 8681, 8689, 8693, 8699, 8707, 8713, 8719, 8731, 8737, 8741, & - 8747, 8753, 8761, 8779, 8783, 8803, 8807, 8819, 8821, 8831, & - ! 1101:1200 - 8837, 8839, 8849, 8861, 8863, 8867, 8887, 8893, 8923, 8929, & - 8933, 8941, 8951, 8963, 8969, 8971, 8999, 9001, 9007, 9011, & - 9013, 9029, 9041, 9043, 9049, 9059, 9067, 9091, 9103, 9109, & - 9127, 9133, 9137, 9151, 9157, 9161, 9173, 9181, 9187, 9199, & - 9203, 9209, 9221, 9227, 9239, 9241, 9257, 9277, 9281, 9283, & - 9293, 9311, 9319, 9323, 9337, 9341, 9343, 9349, 9371, 9377, & - 9391, 9397, 9403, 9413, 9419, 9421, 9431, 9433, 9437, 9439, & - 9461, 9463, 9467, 9473, 9479, 9491, 9497, 9511, 9521, 9533, & - 9539, 9547, 9551, 9587, 9601, 9613, 9619, 9623, 9629, 9631, & - 9643, 9649, 9661, 9677, 9679, 9689, 9697, 9719, 9721, 9733, & - ! 1201:1300 - 9739, 9743, 9749, 9767, 9769, 9781, 9787, 9791, 9803, 9811, & - 9817, 9829, 9833, 9839, 9851, 9857, 9859, 9871, 9883, 9887, & - 9901, 9907, 9923, 9929, 9931, 9941, 9949, 9967, 9973, 10007, & - 10009, 10037, 10039, 10061, 10067, 10069, 10079, 10091, 10093, 10099, & - 10103, 10111, 10133, 10139, 10141, 10151, 10159, 10163, 10169, 10177, & - 10181, 10193, 10211, 10223, 10243, 10247, 10253, 10259, 10267, 10271, & - 10273, 10289, 10301, 10303, 10313, 10321, 10331, 10333, 10337, 10343, & - 10357, 10369, 10391, 10399, 10427, 10429, 10433, 10453, 10457, 10459, & - 10463, 10477, 10487, 10499, 10501, 10513, 10529, 10531, 10559, 10567, & - 10589, 10597, 10601, 10607, 10613, 10627, 10631, 10639, 10651, 10657, & - ! 1301:1400 - 10663, 10667, 10687, 10691, 10709, 10711, 10723, 10729, 10733, 10739, & - 10753, 10771, 10781, 10789, 10799, 10831, 10837, 10847, 10853, 10859, & - 10861, 10867, 10883, 10889, 10891, 10903, 10909, 10937, 10939, 10949, & - 10957, 10973, 10979, 10987, 10993, 11003, 11027, 11047, 11057, 11059, & - 11069, 11071, 11083, 11087, 11093, 11113, 11117, 11119, 11131, 11149, & - 11159, 11161, 11171, 11173, 11177, 11197, 11213, 11239, 11243, 11251, & - 11257, 11261, 11273, 11279, 11287, 11299, 11311, 11317, 11321, 11329, & - 11351, 11353, 11369, 11383, 11393, 11399, 11411, 11423, 11437, 11443, & - 11447, 11467, 11471, 11483, 11489, 11491, 11497, 11503, 11519, 11527, & - 11549, 11551, 11579, 11587, 11593, 11597, 11617, 11621, 11633, 11657, & - ! 1401:1500 - 11677, 11681, 11689, 11699, 11701, 11717, 11719, 11731, 11743, 11777, & - 11779, 11783, 11789, 11801, 11807, 11813, 11821, 11827, 11831, 11833, & - 11839, 11863, 11867, 11887, 11897, 11903, 11909, 11923, 11927, 11933, & - 11939, 11941, 11953, 11959, 11969, 11971, 11981, 11987, 12007, 12011, & - 12037, 12041, 12043, 12049, 12071, 12073, 12097, 12101, 12107, 12109, & - 12113, 12119, 12143, 12149, 12157, 12161, 12163, 12197, 12203, 12211, & - 12227, 12239, 12241, 12251, 12253, 12263, 12269, 12277, 12281, 12289, & - 12301, 12323, 12329, 12343, 12347, 12373, 12377, 12379, 12391, 12401, & - 12409, 12413, 12421, 12433, 12437, 12451, 12457, 12473, 12479, 12487, & - 12491, 12497, 12503, 12511, 12517, 12527, 12539, 12541, 12547, 12553, & - ! 1501:1600 - 12569, 12577, 12583, 12589, 12601, 12611, 12613, 12619, 12637, 12641, & - 12647, 12653, 12659, 12671, 12689, 12697, 12703, 12713, 12721, 12739, & - 12743, 12757, 12763, 12781, 12791, 12799, 12809, 12821, 12823, 12829, & - 12841, 12853, 12889, 12893, 12899, 12907, 12911, 12917, 12919, 12923, & - 12941, 12953, 12959, 12967, 12973, 12979, 12983, 13001, 13003, 13007, & - 13009, 13033, 13037, 13043, 13049, 13063, 13093, 13099, 13103, 13109, & - 13121, 13127, 13147, 13151, 13159, 13163, 13171, 13177, 13183, 13187, & - 13217, 13219, 13229, 13241, 13249, 13259, 13267, 13291, 13297, 13309, & - 13313, 13327, 13331, 13337, 13339, 13367, 13381, 13397, 13399, 13411, & - 13417, 13421, 13441, 13451, 13457, 13463, 13469, 13477, 13487, 13499] - - current = current + 1 - - base = prime(bases) - base_inv = 1.0_pReal/real(base,pReal) - - halton = 0.0_pReal - t = current - - do while (any( t /= 0) ) - halton = halton + real(mod(t,base), pReal) * base_inv - base_inv = base_inv / real(base, pReal) - t = t / base - enddo - -end function halton - - !-------------------------------------------------------------------------------------------------- !> @brief factorial !-------------------------------------------------------------------------------------------------- integer pure function math_factorial(n) - implicit none - integer, intent(in) :: n - integer :: i - - math_factorial = product([(i, i=1,n)]) + integer, intent(in) :: n + integer :: i + + math_factorial = product([(i, i=1,n)]) end function math_factorial @@ -2516,12 +1541,11 @@ end function math_factorial !-------------------------------------------------------------------------------------------------- integer pure function math_binomial(n,k) - implicit none - integer, intent(in) :: n, k - integer :: i, j - - j = min(k,n-k) - math_binomial = product([(i, i=n, n-j+1, -1)])/math_factorial(j) + integer, intent(in) :: n, k + integer :: i, j + + j = min(k,n-k) + math_binomial = product([(i, i=n, n-j+1, -1)])/math_factorial(j) end function math_binomial @@ -2531,14 +1555,13 @@ end function math_binomial !-------------------------------------------------------------------------------------------------- integer pure function math_multinomial(alpha) - implicit none - integer, intent(in), dimension(:) :: alpha - integer :: i + integer, intent(in), dimension(:) :: alpha + integer :: i - math_multinomial = 1 - do i = 1, size(alpha) - math_multinomial = math_multinomial*math_binomial(sum(alpha(1:i)),alpha(i)) - enddo + math_multinomial = 1 + do i = 1, size(alpha) + math_multinomial = math_multinomial*math_binomial(sum(alpha(1:i)),alpha(i)) + enddo end function math_multinomial @@ -2548,15 +1571,14 @@ end function math_multinomial !-------------------------------------------------------------------------------------------------- real(pReal) pure function math_volTetrahedron(v1,v2,v3,v4) - implicit none - real(pReal), dimension (3), intent(in) :: v1,v2,v3,v4 - real(pReal), dimension (3,3) :: m + real(pReal), dimension (3), intent(in) :: v1,v2,v3,v4 + real(pReal), dimension (3,3) :: m - m(1:3,1) = v1-v2 - m(1:3,2) = v2-v3 - m(1:3,3) = v3-v4 + m(1:3,1) = v1-v2 + m(1:3,2) = v2-v3 + m(1:3,3) = v3-v4 - math_volTetrahedron = math_det33(m)/6.0_pReal + math_volTetrahedron = math_det33(m)/6.0_pReal end function math_volTetrahedron @@ -2566,10 +1588,9 @@ end function math_volTetrahedron !-------------------------------------------------------------------------------------------------- real(pReal) pure function math_areaTriangle(v1,v2,v3) - implicit none - real(pReal), dimension (3), intent(in) :: v1,v2,v3 + real(pReal), dimension (3), intent(in) :: v1,v2,v3 - math_areaTriangle = 0.5_pReal * norm2(math_crossproduct(v1-v2,v1-v3)) + math_areaTriangle = 0.5_pReal * norm2(math_cross(v1-v2,v1-v3)) end function math_areaTriangle @@ -2579,12 +1600,10 @@ end function math_areaTriangle !-------------------------------------------------------------------------------------------------- pure function math_rotate_forward33(tensor,rot_tensor) - implicit none + real(pReal), dimension(3,3) :: math_rotate_forward33 + real(pReal), dimension(3,3), intent(in) :: tensor, rot_tensor - real(pReal), dimension(3,3) :: math_rotate_forward33 - real(pReal), dimension(3,3), intent(in) :: tensor, rot_tensor - - math_rotate_forward33 = matmul(rot_tensor,matmul(tensor,transpose(rot_tensor))) + math_rotate_forward33 = matmul(rot_tensor,matmul(tensor,transpose(rot_tensor))) end function math_rotate_forward33 @@ -2594,11 +1613,10 @@ end function math_rotate_forward33 !-------------------------------------------------------------------------------------------------- pure function math_rotate_backward33(tensor,rot_tensor) - implicit none - real(pReal), dimension(3,3) :: math_rotate_backward33 - real(pReal), dimension(3,3), intent(in) :: tensor, rot_tensor + real(pReal), dimension(3,3) :: math_rotate_backward33 + real(pReal), dimension(3,3), intent(in) :: tensor, rot_tensor - math_rotate_backward33 = matmul(transpose(rot_tensor),matmul(tensor,rot_tensor)) + math_rotate_backward33 = matmul(transpose(rot_tensor),matmul(tensor,rot_tensor)) end function math_rotate_backward33 @@ -2608,19 +1626,18 @@ end function math_rotate_backward33 !-------------------------------------------------------------------------------------------------- pure function math_rotate_forward3333(tensor,rot_tensor) - implicit none - real(pReal), dimension(3,3,3,3) :: math_rotate_forward3333 - real(pReal), dimension(3,3), intent(in) :: rot_tensor - real(pReal), dimension(3,3,3,3), intent(in) :: tensor - integer :: i,j,k,l,m,n,o,p + real(pReal), dimension(3,3,3,3) :: math_rotate_forward3333 + real(pReal), dimension(3,3), intent(in) :: rot_tensor + real(pReal), dimension(3,3,3,3), intent(in) :: tensor + integer :: i,j,k,l,m,n,o,p - math_rotate_forward3333 = 0.0_pReal - do i = 1,3;do j = 1,3;do k = 1,3;do l = 1,3 - do m = 1,3;do n = 1,3;do o = 1,3;do p = 1,3 - math_rotate_forward3333(i,j,k,l) & - = math_rotate_forward3333(i,j,k,l) & - + rot_tensor(i,m) * rot_tensor(j,n) * rot_tensor(k,o) * rot_tensor(l,p) * tensor(m,n,o,p) - enddo; enddo; enddo; enddo; enddo; enddo; enddo; enddo + math_rotate_forward3333 = 0.0_pReal + do i = 1,3;do j = 1,3;do k = 1,3;do l = 1,3 + do m = 1,3;do n = 1,3;do o = 1,3;do p = 1,3 + math_rotate_forward3333(i,j,k,l) & + = math_rotate_forward3333(i,j,k,l) & + + rot_tensor(i,m) * rot_tensor(j,n) * rot_tensor(k,o) * rot_tensor(l,p) * tensor(m,n,o,p) + enddo; enddo; enddo; enddo; enddo; enddo; enddo; enddo end function math_rotate_forward3333 @@ -2630,18 +1647,15 @@ end function math_rotate_forward3333 ! Will return NaN if left > right !-------------------------------------------------------------------------------------------------- real(pReal) pure elemental function math_clip(a, left, right) - use, intrinsic :: & - IEEE_arithmetic - implicit none - real(pReal), intent(in) :: a - real(pReal), intent(in), optional :: left, right + real(pReal), intent(in) :: a + real(pReal), intent(in), optional :: left, right - math_clip = a - if (present(left)) math_clip = max(left,math_clip) - if (present(right)) math_clip = min(right,math_clip) - if (present(left) .and. present(right)) & - math_clip = merge (IEEE_value(1.0_pReal,IEEE_quiet_NaN),math_clip, left>right) + math_clip = a + if (present(left)) math_clip = max(left,math_clip) + if (present(right)) math_clip = min(right,math_clip) + if (present(left) .and. present(right)) & + math_clip = merge (IEEE_value(1.0_pReal,IEEE_quiet_NaN),math_clip, left>right) end function math_clip diff --git a/src/mesh/FEM_mech.f90 b/src/mesh/FEM_mech.f90 index 68a24d548..dd9872110 100644 --- a/src/mesh/FEM_mech.f90 +++ b/src/mesh/FEM_mech.f90 @@ -55,8 +55,7 @@ module FEM_mech public :: & FEM_mech_init, & FEM_mech_solution ,& - FEM_mech_forward, & - FEM_mech_destroy + FEM_mech_forward contains !-------------------------------------------------------------------------------------------------- @@ -583,6 +582,7 @@ subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr) end subroutine FEM_mech_formJacobian + !-------------------------------------------------------------------------------------------------- !> @brief forwarding routine !-------------------------------------------------------------------------------------------------- @@ -655,7 +655,6 @@ end subroutine FEM_mech_forward !-------------------------------------------------------------------------------------------------- subroutine FEM_mech_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr) use numerics, only: & - worldrank, & err_struct_tolAbs, & err_struct_tolRel use IO, only: & @@ -677,30 +676,13 @@ subroutine FEM_mech_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dumm call SNESConvergedDefault(snes_local,PETScIter,xnorm,snorm,fnorm/divTol,reason,dummy,ierr) CHKERRQ(ierr) if (terminallyIll) reason = SNES_DIVERGED_FUNCTION_DOMAIN - if (worldrank == 0) then - write(6,'(1/,1x,a,a,i0,a,i0,f0.3)') trim(incInfo), & - ' @ Iteration ',PETScIter,' mechanical residual norm = ', & - int(fnorm/divTol),fnorm/divTol-int(fnorm/divTol) - write(6,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress / MPa =',& + write(6,'(1/,1x,a,a,i0,a,i0,f0.3)') trim(incInfo), & + ' @ Iteration ',PETScIter,' mechanical residual norm = ', & + int(fnorm/divTol),fnorm/divTol-int(fnorm/divTol) + write(6,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress / MPa =',& transpose(P_av)*1.e-6_pReal - flush(6) - endif + flush(6) end subroutine FEM_mech_converged - -!-------------------------------------------------------------------------------------------------- -!> @brief destroy routine -!-------------------------------------------------------------------------------------------------- -subroutine FEM_mech_destroy() - - implicit none - PetscErrorCode :: ierr - - call VecDestroy(solution,ierr); CHKERRQ(ierr) - call VecDestroy(solution_rate,ierr); CHKERRQ(ierr) - call SNESDestroy(mech_snes,ierr); CHKERRQ(ierr) - -end subroutine FEM_mech_destroy - end module FEM_mech diff --git a/src/mesh/FEM_utilities.f90 b/src/mesh/FEM_utilities.f90 index c0e37001e..a2ba2d345 100644 --- a/src/mesh/FEM_utilities.f90 +++ b/src/mesh/FEM_utilities.f90 @@ -23,7 +23,6 @@ use PETScis !-------------------------------------------------------------------------------------------------- ! grid related information information real(pReal), public :: wgt !< weighting factor 1/Nelems - real(pReal), public :: wgtDof !< weighting factor 1/Nelems !-------------------------------------------------------------------------------------------------- ! output data @@ -35,33 +34,18 @@ use PETScis enum, bind(c) enumerator :: FIELD_UNDEFINED_ID, & - FIELD_MECH_ID, & - FIELD_THERMAL_ID, & - FIELD_DAMAGE_ID, & - FIELD_SOLUTE_ID, & - FIELD_MGTWIN_ID + FIELD_MECH_ID end enum enum, bind(c) enumerator :: COMPONENT_UNDEFINED_ID, & COMPONENT_MECH_X_ID, & COMPONENT_MECH_Y_ID, & - COMPONENT_MECH_Z_ID, & - COMPONENT_THERMAL_T_ID, & - COMPONENT_DAMAGE_PHI_ID, & - COMPONENT_SOLUTE_CV_ID, & - COMPONENT_SOLUTE_CVPOT_ID, & - COMPONENT_SOLUTE_CH_ID, & - COMPONENT_SOLUTE_CHPOT_ID, & - COMPONENT_SOLUTE_CVaH_ID, & - COMPONENT_SOLUTE_CVaHPOT_ID, & - COMPONENT_MGTWIN_PHI_ID + COMPONENT_MECH_Z_ID end enum !-------------------------------------------------------------------------------------------------- ! variables controlling debugging logical, private :: & - debugGeneral, & !< general debugging of FEM solver - debugRotation, & !< also printing out results in lab frame debugPETSc !< use some in debug defined options for more verbose PETSc solution !-------------------------------------------------------------------------------------------------- @@ -111,36 +95,26 @@ use PETScis public :: & utilities_init, & utilities_constitutiveResponse, & - utilities_indexBoundaryDofs, & utilities_projectBCValues, & - utilities_indexActiveSet, & - utilities_destroy, & FIELD_MECH_ID, & COMPONENT_MECH_X_ID, & COMPONENT_MECH_Y_ID, & - COMPONENT_MECH_Z_ID, & - COMPONENT_THERMAL_T_ID + COMPONENT_MECH_Z_ID contains !-------------------------------------------------------------------------------------------------- !> @brief allocates all neccessary fields, sets debug flags !-------------------------------------------------------------------------------------------------- -subroutine utilities_init() +subroutine utilities_init use numerics, only: & structOrder, & - integrationOrder, & - worldsize, & - worldrank, & petsc_defaultOptions, & petsc_options use debug, only: & debug_level, & debug_SPECTRAL, & - debug_LEVELBASIC, & - debug_SPECTRALPETSC, & - debug_SPECTRALROTATION - use debug, only: & + debug_SPECTRALPETSC,& PETSCDEBUG use math ! must use the whole module for use of FFTW use mesh, only: & @@ -151,16 +125,12 @@ subroutine utilities_init() implicit none character(len=1024) :: petsc_optionsPhysics - integer(pInt) :: dimPlex - PetscInt :: dim PetscErrorCode :: ierr write(6,'(/,a)') ' <<<+- DAMASK_FEM_utilities init -+>>>' !-------------------------------------------------------------------------------------------------- ! set debugging parameters - debugGeneral = iand(debug_level(debug_SPECTRAL),debug_LEVELBASIC) /= 0 - debugRotation = iand(debug_level(debug_SPECTRAL),debug_SPECTRALROTATION) /= 0 debugPETSc = iand(debug_level(debug_SPECTRAL),debug_SPECTRALPETSC) /= 0 if(debugPETSc) write(6,'(3(/,a),/)') & ' Initializing PETSc with debug options: ', & @@ -180,7 +150,6 @@ subroutine utilities_init() wgt = 1.0/real(mesh_maxNips*mesh_NcpElemsGlobal,pReal) - call DMGetDimension(geomMesh,dimPlex,ierr); CHKERRQ(ierr) end subroutine utilities_init @@ -188,21 +157,13 @@ end subroutine utilities_init !> @brief calculates constitutive response !-------------------------------------------------------------------------------------------------- subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData) - use debug, only: & - debug_reset, & - debug_info use math, only: & - math_transpose33, & - math_rotate_forward33, & math_det33 use FEsolving, only: & restartWrite use homogenization, only: & - materialpoint_F, & materialpoint_P, & materialpoint_stressAndItsTangent - use mesh, only: & - mesh_NcpElems implicit none real(pReal), intent(in) :: timeinc !< loading time @@ -213,9 +174,6 @@ subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData) logical :: & age - integer(pInt) :: & - j - real(pReal) :: defgradDetMin, defgradDetMax, defgradDet PetscErrorCode :: ierr write(6,'(/,a)') ' ... evaluating constitutive response ......................................' @@ -227,27 +185,9 @@ subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData) if (cutBack) then ! restore saved variables age = .False. endif - call debug_reset() - -!-------------------------------------------------------------------------------------------------- -! calculate bounds of det(F) and report - if(debugGeneral) then - defgradDetMax = -huge(1.0_pReal) - defgradDetMin = +huge(1.0_pReal) - do j = 1_pInt, mesh_NcpElems - defgradDet = math_det33(materialpoint_F(1:3,1:3,1,j)) - defgradDetMax = max(defgradDetMax,defgradDet) - defgradDetMin = min(defgradDetMin,defgradDet) - end do - write(6,'(a,1x,es11.4)') ' max determinant of deformation =', defgradDetMax - write(6,'(a,1x,es11.4)') ' min determinant of deformation =', defgradDetMin - flush(6) - endif call materialpoint_stressAndItsTangent(.true.,timeinc) ! calculate P field - call debug_info() - restartWrite = .false. ! reset restartWrite status cutBack = .false. ! reset cutBack status @@ -257,97 +197,6 @@ subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData) end subroutine utilities_constitutiveResponse -!-------------------------------------------------------------------------------------------------- -!> @brief Create index sets of boundary dofs (in local and global numbering) -!-------------------------------------------------------------------------------------------------- -subroutine utilities_indexBoundaryDofs(dm_local,nFaceSets,numFields,local2global,section,localIS,globalIS) - - implicit none - - DM :: dm_local - ISLocalToGlobalMapping :: local2global - PetscSection :: section - PetscInt :: nFaceSets, numFields, nDof - IS, dimension(nFaceSets,numFields) :: localIS, globalIS - PetscInt :: field, faceSet, point, dof, offset - PetscInt :: localSize, storageSize, ISSize - PetscInt, dimension(:) , allocatable :: localIndices - IS :: faceSetIS, BC_IS, dummyIS - PetscInt, dimension(:) , pointer :: pFaceSets, pBCvertex, pBCvertexlc - DMLabel :: BCLabel - PetscErrorCode :: ierr - - call DMGetLabel(dm_local,'Face Sets',BCLabel,ierr); CHKERRQ(ierr) - call DMPlexLabelComplete(dm_local,BCLabel,ierr); CHKERRQ(ierr) - call PetscSectionGetStorageSize(section,storageSize,ierr); CHKERRQ(ierr) - call DMGetLabelIdIS(dm_local,'Face Sets',faceSetIS,ierr); CHKERRQ(ierr) - call ISGetIndicesF90(faceSetIS,pFaceSets,ierr); CHKERRQ(ierr) - allocate(localIndices (storageSize)) - do faceSet = 1, nFaceSets - call DMGetStratumSize(dm_local,'Face Sets',pFaceSets(faceSet),ISSize,ierr) - CHKERRQ(ierr) - call DMGetStratumIS(dm_local,'Face Sets',pFaceSets(faceSet),BC_IS,ierr) - CHKERRQ(ierr) - if (ISSize > 0) call ISGetIndicesF90(BC_IS,pBCvertex,ierr) - do field = 1, numFields - localSize = 0 - do point = 1, ISSize - call PetscSectionGetFieldDof(section,pBCvertex(point),field-1,nDof,ierr) - CHKERRQ(ierr) - call PetscSectionGetFieldOffset(section,pBCvertex(point),field-1,offset,ierr) - CHKERRQ(ierr) - do dof = 1, nDof - localSize = localSize + 1 - localIndices(localSize) = offset + dof - 1 - enddo - enddo - call ISCreateGeneral(PETSC_COMM_SELF,localSize,localIndices,PETSC_COPY_VALUES, & - localIS(faceSet,field),ierr) - CHKERRQ(ierr) - call ISLocalToGlobalMappingApplyIS(local2global,localIS(faceSet,field), & - globalIS(faceSet,field),ierr) - CHKERRQ(ierr) - enddo - if (ISSize > 0) call ISRestoreIndicesF90(BC_IS,pBCvertex,ierr) - call ISDestroy(BC_IS,ierr); CHKERRQ(ierr) - enddo - call ISRestoreIndicesF90(faceSetIS,pFaceSets,ierr); CHKERRQ(ierr) - call ISDestroy(faceSetIS,ierr); CHKERRQ(ierr) - - do faceSet = 1, nFaceSets; do field = 1, numFields - call ISGetSize(globalIS(faceSet,field),ISSize,ierr); CHKERRQ(ierr) - if (ISSize > 0) then - call ISGetIndicesF90(localIS(faceSet,field),pBCvertexlc,ierr); CHKERRQ(ierr) - call ISGetIndicesF90(globalIS(faceSet,field),pBCvertex,ierr); CHKERRQ(ierr) - endif - localSize = 0 - do point = 1, ISSize - if (pBCvertex(point) >= 0) then - localSize = localSize + 1 - localIndices(localSize) = pBCvertexlc(point) - endif - enddo - if (ISSize > 0) then - call ISRestoreIndicesF90(localIS(faceSet,field),pBCvertexlc,ierr); CHKERRQ(ierr) - call ISRestoreIndicesF90(globalIS(faceSet,field),pBCvertex,ierr); CHKERRQ(ierr) - endif - call ISDestroy(globalIS(faceSet,field),ierr); CHKERRQ(ierr) - call ISCreateGeneral(PETSC_COMM_SELF,localSize,localIndices,PETSC_COPY_VALUES, & - globalIS(faceSet,field),ierr) - CHKERRQ(ierr) - if (ISSize > 0) then - call ISDuplicate(localIS(faceSet,field),dummyIS,ierr); CHKERRQ(ierr) - call ISDestroy(localIS(faceSet,field),ierr); CHKERRQ(ierr) - call ISDifference(dummyIS,globalIS(faceSet,field),localIS(faceSet,field),ierr) - CHKERRQ(ierr) - call ISDestroy(dummyIS,ierr); CHKERRQ(ierr) - endif - enddo; enddo - deallocate(localIndices) - -end subroutine utilities_indexBoundaryDofs - - !-------------------------------------------------------------------------------------------------- !> @brief Project BC values to local vector !-------------------------------------------------------------------------------------------------- @@ -384,104 +233,4 @@ subroutine utilities_projectBCValues(localVec,section,field,comp,bcPointsIS,BCVa end subroutine utilities_projectBCValues -!-------------------------------------------------------------------------------------------------- -!> @brief Create index sets of boundary dofs (in local and global numbering) -!-------------------------------------------------------------------------------------------------- -subroutine utilities_indexActiveSet(field,section,x_local,f_local,localIS,globalIS) - use mesh, only: & - geomMesh - - implicit none - - ISLocalToGlobalMapping :: local2global - PetscSection :: section - Vec :: x_local, f_local - PetscInt :: field - IS :: localIS, globalIS, dummyIS - PetscScalar, dimension(:) , pointer :: x_scal, f_scal - PetscInt :: ISSize - PetscInt :: chart, chartStart, chartEnd, nDof, dof, offset - PetscInt :: localSize - PetscInt, dimension(:) , allocatable :: localIndices - PetscInt, dimension(:) , pointer :: pBCvertex, pBCvertexlc - PetscErrorCode :: ierr - - call DMGetLocalToGlobalMapping(geomMesh,local2global,ierr) - CHKERRQ(ierr) - call DMPlexGetChart(geomMesh,chartStart,chartEnd,ierr) - CHKERRQ(ierr) - call VecGetArrayF90(x_local,x_scal,ierr); CHKERRQ(ierr) - call VecGetArrayF90(f_local,f_scal,ierr); CHKERRQ(ierr) - localSize = 0 - do chart = chartStart, chartEnd-1 - call PetscSectionGetFieldDof(section,chart,field-1,nDof,ierr); CHKERRQ(ierr) - call PetscSectionGetFieldOffset(section,chart,field-1,offset,ierr); CHKERRQ(ierr) - do dof = offset+1, offset+nDof - if (((x_scal(dof) < 1.0e-8) .and. (f_scal(dof) > 0.0)) .or. & - ((x_scal(dof) > 1.0 - 1.0e-8) .and. (f_scal(dof) < 0.0))) localSize = localSize + 1 - enddo - enddo - allocate(localIndices(localSize)) - localSize = 0 - do chart = chartStart, chartEnd-1 - call PetscSectionGetFieldDof(section,chart,field-1,nDof,ierr); CHKERRQ(ierr) - call PetscSectionGetFieldOffset(section,chart,field-1,offset,ierr); CHKERRQ(ierr) - do dof = offset+1, offset+nDof - if (((x_scal(dof) < 1.0e-8) .and. (f_scal(dof) > 0.0)) .or. & - ((x_scal(dof) > 1.0 - 1.0e-8) .and. (f_scal(dof) < 0.0))) then - localSize = localSize + 1 - localIndices(localSize) = dof-1 - endif - enddo - enddo - call VecRestoreArrayF90(x_local,x_scal,ierr); CHKERRQ(ierr) - call VecRestoreArrayF90(f_local,f_scal,ierr); CHKERRQ(ierr) - call ISCreateGeneral(PETSC_COMM_SELF,localSize,localIndices,PETSC_COPY_VALUES,localIS,ierr) - CHKERRQ(ierr) - call ISLocalToGlobalMappingApplyIS(local2global,localIS,globalIS,ierr) - CHKERRQ(ierr) - call ISGetSize(globalIS,ISSize,ierr); CHKERRQ(ierr) - if (ISSize > 0) then - call ISGetIndicesF90(localIS,pBCvertexlc,ierr); CHKERRQ(ierr) - call ISGetIndicesF90(globalIS,pBCvertex,ierr); CHKERRQ(ierr) - endif - localSize = 0 - do chart = 1, ISSize - if (pBCvertex(chart) >= 0) then - localSize = localSize + 1 - localIndices(localSize) = pBCvertexlc(chart) - endif - enddo - if (ISSize > 0) then - call ISRestoreIndicesF90(localIS,pBCvertexlc,ierr); CHKERRQ(ierr) - call ISRestoreIndicesF90(globalIS,pBCvertex,ierr); CHKERRQ(ierr) - endif - call ISDestroy(globalIS,ierr); CHKERRQ(ierr) - call ISCreateGeneral(PETSC_COMM_SELF,localSize,localIndices,PETSC_COPY_VALUES,globalIS,ierr) - CHKERRQ(ierr) - if (ISSize > 0) then - call ISDuplicate(localIS,dummyIS,ierr); CHKERRQ(ierr) - call ISDestroy(localIS,ierr); CHKERRQ(ierr) - call ISDifference(dummyIS,globalIS,localIS,ierr) - CHKERRQ(ierr) - call ISDestroy(dummyIS,ierr); CHKERRQ(ierr) - endif - -end subroutine utilities_indexActiveSet - -!-------------------------------------------------------------------------------------------------- -!> @brief cleans up -!-------------------------------------------------------------------------------------------------- -subroutine utilities_destroy() - - !implicit none - !PetscInt :: homog, cryst, grain, phase - !PetscErrorCode :: ierr - - !call VecDestroy(coordinatesVec,ierr); CHKERRQ(ierr) - !call PetscViewerDestroy(resUnit, ierr); CHKERRQ(ierr) - -end subroutine utilities_destroy - - end module FEM_utilities diff --git a/src/mesh_abaqus.f90 b/src/mesh_abaqus.f90 index adcaa8513..0467d09aa 100644 --- a/src/mesh_abaqus.f90 +++ b/src/mesh_abaqus.f90 @@ -884,7 +884,7 @@ end subroutine mesh_abaqus_count_cpElements !-------------------------------------------------------------------------------------------------- subroutine mesh_abaqus_map_elements(fileUnit) - use math, only: math_qsort + use math, only: math_sort use IO, only: IO_lc, & IO_stringValue, & IO_stringPos, & @@ -935,7 +935,7 @@ subroutine mesh_abaqus_map_elements(fileUnit) endselect enddo - call math_qsort(mesh_mapFEtoCPelem,1_pInt,int(size(mesh_mapFEtoCPelem,2_pInt),pInt)) ! should be mesh_NcpElems + call math_sort(mesh_mapFEtoCPelem,1_pInt,int(size(mesh_mapFEtoCPelem,2_pInt),pInt)) ! should be mesh_NcpElems if (int(size(mesh_mapFEtoCPelem),pInt) < 2_pInt) call IO_error(error_ID=907_pInt) @@ -948,7 +948,7 @@ end subroutine mesh_abaqus_map_elements !-------------------------------------------------------------------------------------------------- subroutine mesh_abaqus_map_nodes(fileUnit) - use math, only: math_qsort + use math, only: math_sort use IO, only: IO_lc, & IO_stringValue, & IO_stringPos, & @@ -999,7 +999,7 @@ subroutine mesh_abaqus_map_nodes(fileUnit) endif enddo - call math_qsort(mesh_mapFEtoCPnode,1_pInt,int(size(mesh_mapFEtoCPnode,2_pInt),pInt)) + call math_sort(mesh_mapFEtoCPnode,1_pInt,int(size(mesh_mapFEtoCPnode,2_pInt),pInt)) if (int(size(mesh_mapFEtoCPnode),pInt) == 0_pInt) call IO_error(error_ID=908_pInt) @@ -1541,7 +1541,7 @@ pure function mesh_cellCenterCoordinates(ip,el) !-------------------------------------------------------------------------------------------------- subroutine mesh_build_ipAreas use math, only: & - math_crossproduct + math_cross implicit none integer(pInt) :: e,t,g,c,i,f,n,m @@ -1576,7 +1576,7 @@ subroutine mesh_build_ipAreas do f = 1_pInt,FE_NipNeighbors(c) ! loop over cell faces forall(n = 1_pInt:FE_NcellnodesPerCellface(c)) & nodePos(1:3,n) = mesh_cellnode(1:3,mesh_cell(FE_cellface(n,f,c),i,e)) - normal = math_crossproduct(nodePos(1:3,2) - nodePos(1:3,1), & + normal = math_cross(nodePos(1:3,2) - nodePos(1:3,1), & nodePos(1:3,3) - nodePos(1:3,1)) mesh_ipArea(f,i,e) = norm2(normal) mesh_ipAreaNormal(1:3,f,i,e) = normal / norm2(normal) ! ensure unit length of area normal @@ -1595,7 +1595,7 @@ subroutine mesh_build_ipAreas nodePos(1:3,n) = mesh_cellnode(1:3,mesh_cell(FE_cellface(n,f,c),i,e)) forall(n = 1_pInt:FE_NcellnodesPerCellface(c)) & normals(1:3,n) = 0.5_pReal & - * math_crossproduct(nodePos(1:3,1+mod(n ,m)) - nodePos(1:3,n), & + * math_cross(nodePos(1:3,1+mod(n ,m)) - nodePos(1:3,n), & nodePos(1:3,1+mod(n+1,m)) - nodePos(1:3,n)) normal = 0.5_pReal * sum(normals,2) mesh_ipArea(f,i,e) = norm2(normal) diff --git a/src/mesh_marc.f90 b/src/mesh_marc.f90 index c999e32d5..79718c37f 100644 --- a/src/mesh_marc.f90 +++ b/src/mesh_marc.f90 @@ -640,7 +640,7 @@ subroutine mesh_marc_map_elementSets(nameElemSet,mapElemSet,fileUnit) !-------------------------------------------------------------------------------------------------- subroutine mesh_marc_map_elements(tableStyle,nameElemSet,mapElemSet,nElems,fileUnit) - use math, only: math_qsort + use math, only: math_sort use IO, only: IO_lc, & IO_intValue, & IO_stringValue, & @@ -701,7 +701,7 @@ subroutine mesh_marc_map_elements(tableStyle,nameElemSet,mapElemSet,nElems,fileU mesh_mapFEtoCPelem(2,cpElem) = cpElem enddo -call math_qsort(mesh_mapFEtoCPelem,1_pInt,int(size(mesh_mapFEtoCPelem,2_pInt),pInt)) ! should be mesh_NcpElems +call math_sort(mesh_mapFEtoCPelem,1_pInt,int(size(mesh_mapFEtoCPelem,2_pInt),pInt)) ! should be mesh_NcpElems end subroutine mesh_marc_map_elements @@ -711,7 +711,7 @@ end subroutine mesh_marc_map_elements !-------------------------------------------------------------------------------------------------- subroutine mesh_marc_map_nodes(nNodes,fileUnit) - use math, only: math_qsort + use math, only: math_sort use IO, only: IO_lc, & IO_stringValue, & IO_stringPos, & @@ -743,7 +743,7 @@ subroutine mesh_marc_map_nodes(nNodes,fileUnit) endif enddo -650 call math_qsort(mesh_mapFEtoCPnode,1_pInt,int(size(mesh_mapFEtoCPnode,2_pInt),pInt)) +650 call math_sort(mesh_mapFEtoCPnode,1_pInt,int(size(mesh_mapFEtoCPnode,2_pInt),pInt)) end subroutine mesh_marc_map_nodes @@ -1230,7 +1230,7 @@ pure function mesh_cellCenterCoordinates(ip,el) !-------------------------------------------------------------------------------------------------- subroutine mesh_build_ipAreas use math, only: & - math_crossproduct + math_cross implicit none integer(pInt) :: e,t,g,c,i,f,n,m @@ -1265,7 +1265,7 @@ subroutine mesh_build_ipAreas do f = 1_pInt,FE_NipNeighbors(c) ! loop over cell faces forall(n = 1_pInt:FE_NcellnodesPerCellface(c)) & nodePos(1:3,n) = mesh_cellnode(1:3,mesh_cell(FE_cellface(n,f,c),i,e)) - normal = math_crossproduct(nodePos(1:3,2) - nodePos(1:3,1), & + normal = math_cross(nodePos(1:3,2) - nodePos(1:3,1), & nodePos(1:3,3) - nodePos(1:3,1)) mesh_ipArea(f,i,e) = norm2(normal) mesh_ipAreaNormal(1:3,f,i,e) = normal / norm2(normal) ! ensure unit length of area normal @@ -1284,7 +1284,7 @@ subroutine mesh_build_ipAreas nodePos(1:3,n) = mesh_cellnode(1:3,mesh_cell(FE_cellface(n,f,c),i,e)) forall(n = 1_pInt:FE_NcellnodesPerCellface(c)) & normals(1:3,n) = 0.5_pReal & - * math_crossproduct(nodePos(1:3,1+mod(n ,m)) - nodePos(1:3,n), & + * math_cross(nodePos(1:3,1+mod(n ,m)) - nodePos(1:3,n), & nodePos(1:3,1+mod(n+1,m)) - nodePos(1:3,n)) normal = 0.5_pReal * sum(normals,2) mesh_ipArea(f,i,e) = norm2(normal) diff --git a/src/prec.f90 b/src/prec.f90 index cba8a68ef..2a1c2f14a 100644 --- a/src/prec.f90 +++ b/src/prec.f90 @@ -7,11 +7,9 @@ !> @brief setting precision for real and int type !-------------------------------------------------------------------------------------------------- module prec - use, intrinsic :: IEEE_arithmetic, only:& - IEEE_selected_real_kind + use, intrinsic :: IEEE_arithmetic implicit none - private ! https://software.intel.com/en-us/blogs/2017/03/27/doctor-fortran-in-it-takes-all-kinds #ifdef Abaqus integer, parameter, public :: pReal = selected_real_kind(15,307) !< number with 15 significant digits, up to 1e+-307 (typically 64 bit) @@ -102,7 +100,6 @@ contains !-------------------------------------------------------------------------------------------------- subroutine prec_init - implicit none integer, allocatable, dimension(:) :: realloc_lhs_test external :: & @@ -131,7 +128,6 @@ end subroutine prec_init !-------------------------------------------------------------------------------------------------- logical elemental pure function dEq(a,b,tol) - implicit none real(pReal), intent(in) :: a,b real(pReal), intent(in), optional :: tol real(pReal) :: eps @@ -155,7 +151,6 @@ end function dEq !-------------------------------------------------------------------------------------------------- logical elemental pure function dNeq(a,b,tol) - implicit none real(pReal), intent(in) :: a,b real(pReal), intent(in), optional :: tol real(pReal) :: eps @@ -179,7 +174,6 @@ end function dNeq !-------------------------------------------------------------------------------------------------- logical elemental pure function dEq0(a,tol) - implicit none real(pReal), intent(in) :: a real(pReal), intent(in), optional :: tol real(pReal) :: eps @@ -203,7 +197,6 @@ end function dEq0 !-------------------------------------------------------------------------------------------------- logical elemental pure function dNeq0(a,tol) - implicit none real(pReal), intent(in) :: a real(pReal), intent(in), optional :: tol real(pReal) :: eps @@ -228,7 +221,6 @@ end function dNeq0 !-------------------------------------------------------------------------------------------------- logical elemental pure function cEq(a,b,tol) - implicit none complex(pReal), intent(in) :: a,b real(pReal), intent(in), optional :: tol real(pReal) :: eps @@ -253,7 +245,6 @@ end function cEq !-------------------------------------------------------------------------------------------------- logical elemental pure function cNeq(a,b,tol) - implicit none complex(pReal), intent(in) :: a,b real(pReal), intent(in), optional :: tol real(pReal) :: eps diff --git a/src/source_damage_anisoBrittle.f90 b/src/source_damage_anisoBrittle.f90 index 717ea9fb4..da00db2b2 100644 --- a/src/source_damage_anisoBrittle.f90 +++ b/src/source_damage_anisoBrittle.f90 @@ -91,8 +91,6 @@ subroutine source_damage_anisoBrittle_init lattice_SchmidMatrix_cleavage, & lattice_maxNcleavageFamily - implicit none - integer(pInt) :: Ninstance,phase,instance,source,sourceOffset integer(pInt) :: NofMyPhase,p ,i integer(pInt), dimension(0), parameter :: emptyIntArray = [integer(pInt)::] @@ -219,7 +217,6 @@ subroutine source_damage_anisoBrittle_dotState(S, ipc, ip, el) lattice_maxNcleavageFamily, & lattice_NcleavageSystem - implicit none integer(pInt), intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point @@ -279,7 +276,6 @@ subroutine source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalph use material, only: & sourceState - implicit none integer(pInt), intent(in) :: & phase, & constituent @@ -307,7 +303,6 @@ function source_damage_anisoBrittle_postResults(phase, constituent) use material, only: & sourceState - implicit none integer(pInt), intent(in) :: & phase, & constituent diff --git a/src/source_damage_anisoDuctile.f90 b/src/source_damage_anisoDuctile.f90 index f0774df62..8f6d68a88 100644 --- a/src/source_damage_anisoDuctile.f90 +++ b/src/source_damage_anisoDuctile.f90 @@ -82,7 +82,6 @@ subroutine source_damage_anisoDuctile_init config_phase - implicit none integer(pInt) :: Ninstance,phase,instance,source,sourceOffset integer(pInt) :: NofMyPhase,p ,i @@ -194,7 +193,6 @@ subroutine source_damage_anisoDuctile_dotState(ipc, ip, el) damage, & damageMapping - implicit none integer(pInt), intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point @@ -231,7 +229,6 @@ subroutine source_damage_anisoDuctile_getRateAndItsTangent(localphiDot, dLocalph use material, only: & sourceState - implicit none integer(pInt), intent(in) :: & phase, & constituent @@ -259,7 +256,6 @@ function source_damage_anisoDuctile_postResults(phase, constituent) use material, only: & sourceState - implicit none integer(pInt), intent(in) :: & phase, & constituent diff --git a/src/source_damage_isoBrittle.f90 b/src/source_damage_isoBrittle.f90 index 2d13277c7..b60218458 100644 --- a/src/source_damage_isoBrittle.f90 +++ b/src/source_damage_isoBrittle.f90 @@ -74,7 +74,6 @@ subroutine source_damage_isoBrittle_init config_phase, & material_Nphase - implicit none integer(pInt) :: Ninstance,phase,instance,source,sourceOffset integer(pInt) :: NofMyPhase,p,i @@ -176,7 +175,6 @@ subroutine source_damage_isoBrittle_deltaState(C, Fe, ipc, ip, el) math_sym33to6, & math_I3 - implicit none integer(pInt), intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point @@ -221,7 +219,6 @@ subroutine source_damage_isoBrittle_getRateAndItsTangent(localphiDot, dLocalphiD use material, only: & sourceState - implicit none integer(pInt), intent(in) :: & phase, & constituent @@ -251,7 +248,6 @@ function source_damage_isoBrittle_postResults(phase, constituent) use material, only: & sourceState - implicit none integer(pInt), intent(in) :: & phase, & constituent diff --git a/src/source_damage_isoDuctile.f90 b/src/source_damage_isoDuctile.f90 index 279d436c9..149803693 100644 --- a/src/source_damage_isoDuctile.f90 +++ b/src/source_damage_isoDuctile.f90 @@ -74,7 +74,6 @@ subroutine source_damage_isoDuctile_init config_phase, & material_Nphase - implicit none integer(pInt) :: Ninstance,phase,instance,source,sourceOffset integer(pInt) :: NofMyPhase,p,i @@ -177,7 +176,6 @@ subroutine source_damage_isoDuctile_dotState(ipc, ip, el) damage, & damageMapping - implicit none integer(pInt), intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point @@ -206,7 +204,6 @@ subroutine source_damage_isoDuctile_getRateAndItsTangent(localphiDot, dLocalphiD use material, only: & sourceState - implicit none integer(pInt), intent(in) :: & phase, & constituent @@ -234,7 +231,6 @@ function source_damage_isoDuctile_postResults(phase, constituent) use material, only: & sourceState - implicit none integer(pInt), intent(in) :: & phase, & constituent diff --git a/src/source_thermal_dissipation.f90 b/src/source_thermal_dissipation.f90 index 23d6e8b58..94452eb47 100644 --- a/src/source_thermal_dissipation.f90 +++ b/src/source_thermal_dissipation.f90 @@ -56,7 +56,6 @@ subroutine source_thermal_dissipation_init config_phase, & material_Nphase - implicit none integer :: Ninstance,instance,source,sourceOffset integer :: NofMyPhase,p @@ -103,7 +102,6 @@ end subroutine source_thermal_dissipation_init !-------------------------------------------------------------------------------------------------- subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDOT_dT, Tstar, Lp, phase) - implicit none integer, intent(in) :: & phase real(pReal), intent(in), dimension(3,3) :: & diff --git a/src/source_thermal_externalheat.f90 b/src/source_thermal_externalheat.f90 index 377513c57..699902ad3 100644 --- a/src/source_thermal_externalheat.f90 +++ b/src/source_thermal_externalheat.f90 @@ -63,7 +63,6 @@ subroutine source_thermal_externalheat_init config_phase, & material_Nphase - implicit none integer :: maxNinstance,instance,source,sourceOffset,NofMyPhase,p @@ -120,7 +119,6 @@ subroutine source_thermal_externalheat_dotState(phase, of) use material, only: & sourceState - implicit none integer, intent(in) :: & phase, & of @@ -140,7 +138,6 @@ subroutine source_thermal_externalheat_getRateAndItsTangent(TDot, dTDot_dT, phas use material, only: & sourceState - implicit none integer, intent(in) :: & phase, & of diff --git a/src/system_routines.f90 b/src/system_routines.f90 index d7a27a4f9..0611c96db 100644 --- a/src/system_routines.f90 +++ b/src/system_routines.f90 @@ -3,13 +3,9 @@ !> @brief provides wrappers to C routines !-------------------------------------------------------------------------------------------------- module system_routines - use, intrinsic :: ISO_C_Binding, only: & - C_INT, & - C_CHAR, & - C_NULL_CHAR + use, intrinsic :: ISO_C_Binding implicit none - private public :: & signalterm_C, & @@ -81,16 +77,15 @@ contains !-------------------------------------------------------------------------------------------------- logical function isDirectory(path) - implicit none - character(len=*), intent(in) :: path - character(kind=C_CHAR), dimension(1024) :: strFixedLength ! C string as array - integer :: i - - strFixedLength = repeat(C_NULL_CHAR,len(strFixedLength)) - do i=1,len(path) ! copy array components - strFixedLength(i)=path(i:i) - enddo - isDirectory=merge(.True.,.False.,isDirectory_C(strFixedLength) /= 0_C_INT) + character(len=*), intent(in) :: path + character(kind=C_CHAR), dimension(1024) :: strFixedLength ! C string as array + integer :: i + + strFixedLength = repeat(C_NULL_CHAR,len(strFixedLength)) + do i=1,len(path) ! copy array components + strFixedLength(i)=path(i:i) + enddo + isDirectory=merge(.True.,.False.,isDirectory_C(strFixedLength) /= 0_C_INT) end function isDirectory @@ -100,24 +95,23 @@ end function isDirectory !-------------------------------------------------------------------------------------------------- character(len=1024) function getCWD() - implicit none - character(kind=C_CHAR), dimension(1024) :: charArray ! C string is an array - integer(C_INT) :: stat - integer :: i - - call getCurrentWorkDir_C(charArray,stat) - if (stat /= 0_C_INT) then - getCWD = 'Error occured when getting currend working directory' - else - getCWD = repeat('',len(getCWD)) - arrayToString: do i=1,len(getCWD) - if (charArray(i) /= C_NULL_CHAR) then - getCWD(i:i)=charArray(i) - else - exit - endif - enddo arrayToString - endif + character(kind=C_CHAR), dimension(1024) :: charArray ! C string is an array + integer(C_INT) :: stat + integer :: i + + call getCurrentWorkDir_C(charArray,stat) + if (stat /= 0_C_INT) then + getCWD = 'Error occured when getting currend working directory' + else + getCWD = repeat('',len(getCWD)) + arrayToString: do i=1,len(getCWD) + if (charArray(i) /= C_NULL_CHAR) then + getCWD(i:i)=charArray(i) + else + exit + endif + enddo arrayToString + endif end function getCWD @@ -126,24 +120,24 @@ end function getCWD !> @brief gets the current host name !-------------------------------------------------------------------------------------------------- character(len=1024) function getHostName() - implicit none - character(kind=C_CHAR), dimension(1024) :: charArray ! C string is an array - integer(C_INT) :: stat - integer :: i - call getHostName_C(charArray,stat) - if (stat /= 0_C_INT) then - getHostName = 'Error occured when getting host name' - else - getHostName = repeat('',len(getHostName)) - arrayToString: do i=1,len(getHostName) - if (charArray(i) /= C_NULL_CHAR) then - getHostName(i:i)=charArray(i) - else - exit - endif - enddo arrayToString - endif + character(kind=C_CHAR), dimension(1024) :: charArray ! C string is an array + integer(C_INT) :: stat + integer :: i + + call getHostName_C(charArray,stat) + if (stat /= 0_C_INT) then + getHostName = 'Error occured when getting host name' + else + getHostName = repeat('',len(getHostName)) + arrayToString: do i=1,len(getHostName) + if (charArray(i) /= C_NULL_CHAR) then + getHostName(i:i)=charArray(i) + else + exit + endif + enddo arrayToString + endif end function getHostName @@ -152,16 +146,16 @@ end function getHostName !> @brief changes the current working directory !-------------------------------------------------------------------------------------------------- logical function setCWD(path) - implicit none - character(len=*), intent(in) :: path - character(kind=C_CHAR), dimension(1024) :: strFixedLength ! C string is an array - integer :: i - strFixedLength = repeat(C_NULL_CHAR,len(strFixedLength)) - do i=1,len(path) ! copy array components - strFixedLength(i)=path(i:i) - enddo - setCWD=merge(.True.,.False.,chdir_C(strFixedLength) /= 0_C_INT) + character(len=*), intent(in) :: path + character(kind=C_CHAR), dimension(1024) :: strFixedLength ! C string is an array + integer :: i + + strFixedLength = repeat(C_NULL_CHAR,len(strFixedLength)) + do i=1,len(path) ! copy array components + strFixedLength(i)=path(i:i) + enddo + setCWD=merge(.True.,.False.,chdir_C(strFixedLength) /= 0_C_INT) end function setCWD diff --git a/src/thermal_adiabatic.f90 b/src/thermal_adiabatic.f90 index 8fdf55c97..bfc34d1c4 100644 --- a/src/thermal_adiabatic.f90 +++ b/src/thermal_adiabatic.f90 @@ -57,7 +57,6 @@ subroutine thermal_adiabatic_init use config, only: & config_homogenization - implicit none integer :: maxNinstance,section,instance,i,sizeState,NofMyHomog character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::] character(len=65536), dimension(:), allocatable :: outputs @@ -124,7 +123,6 @@ function thermal_adiabatic_updateState(subdt, ip, el) temperatureRate, & thermalMapping - implicit none integer, intent(in) :: & ip, & !< integration point number el !< element number @@ -181,7 +179,6 @@ subroutine thermal_adiabatic_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el) crystallite_S, & crystallite_Lp - implicit none integer, intent(in) :: & ip, & !< integration point number el !< element number @@ -246,7 +243,6 @@ function thermal_adiabatic_getSpecificHeat(ip,el) use mesh, only: & mesh_element - implicit none integer, intent(in) :: & ip, & !< integration point number el !< element number @@ -282,7 +278,6 @@ function thermal_adiabatic_getMassDensity(ip,el) use mesh, only: & mesh_element - implicit none integer, intent(in) :: & ip, & !< integration point number el !< element number @@ -312,7 +307,6 @@ function thermal_adiabatic_postResults(homog,instance,of) result(postResults) use material, only: & temperature - implicit none integer, intent(in) :: & homog, & instance, & diff --git a/src/thermal_conduction.f90 b/src/thermal_conduction.f90 index 8dd352357..292d8a3f6 100644 --- a/src/thermal_conduction.f90 +++ b/src/thermal_conduction.f90 @@ -58,7 +58,6 @@ subroutine thermal_conduction_init use config, only: & config_homogenization - implicit none integer :: maxNinstance,section,instance,i integer :: sizeState integer :: NofMyHomog @@ -135,7 +134,6 @@ subroutine thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el) crystallite_S, & crystallite_Lp - implicit none integer, intent(in) :: & ip, & !< integration point number el !< element number @@ -205,7 +203,6 @@ function thermal_conduction_getConductivity33(ip,el) use crystallite, only: & crystallite_push33ToRef - implicit none integer, intent(in) :: & ip, & !< integration point number el !< element number @@ -239,7 +236,6 @@ function thermal_conduction_getSpecificHeat(ip,el) use mesh, only: & mesh_element - implicit none integer, intent(in) :: & ip, & !< integration point number el !< element number @@ -273,7 +269,6 @@ function thermal_conduction_getMassDensity(ip,el) use mesh, only: & mesh_element - implicit none integer, intent(in) :: & ip, & !< integration point number el !< element number @@ -306,7 +301,6 @@ subroutine thermal_conduction_putTemperatureAndItsRate(T,Tdot,ip,el) temperatureRate, & thermalMapping - implicit none integer, intent(in) :: & ip, & !< integration point number el !< element number @@ -332,7 +326,6 @@ function thermal_conduction_postResults(homog,instance,of) result(postResults) use material, only: & temperature - implicit none integer, intent(in) :: & homog, & instance, & diff --git a/src/thermal_isothermal.f90 b/src/thermal_isothermal.f90 index a161094e7..c15dfa8f1 100644 --- a/src/thermal_isothermal.f90 +++ b/src/thermal_isothermal.f90 @@ -22,7 +22,6 @@ subroutine thermal_isothermal_init() material_Nhomogenization use material - implicit none integer :: & homog, & NofMyHomog