From e11be7e6007581f349de3172b5584a4574921a06 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Thu, 17 Dec 2020 10:47:56 -0500 Subject: [PATCH 01/19] preinitialize a ConfigMaterial object with 'constituents','homogenization','phase' keys --- python/damask/_configmaterial.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/python/damask/_configmaterial.py b/python/damask/_configmaterial.py index 6de2283f4..17d8c796d 100644 --- a/python/damask/_configmaterial.py +++ b/python/damask/_configmaterial.py @@ -9,6 +9,13 @@ from . import Orientation class ConfigMaterial(Config): """Material configuration.""" + def __init__(self): + """Initialize object with all required dictionary keys.""" + super().__init__() + self['material'] = [] + self['homogenization'] = {} + self['phase'] = {} + def save(self,fname='material.yaml',**kwargs): """ Save to yaml file. @@ -274,6 +281,7 @@ class ConfigMaterial(Config): c = [{} for _ in range(length)] if constituents is None else \ [{'constituents':u} for u in ConfigMaterial._constituents(**constituents)] + if len(c) == 1: c = [copy.deepcopy(c[0]) for _ in range(length)] if length != 1 and length != len(c): From 403ac693dab87d752961bd8732ea77c8f9213204 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Thu, 17 Dec 2020 18:08:55 -0500 Subject: [PATCH 02/19] need to pass init argument to dict superclass --- python/damask/_configmaterial.py | 15 +++++++++------ .../tests/reference/ConfigMaterial/material.yaml | 10 +++++----- 2 files changed, 14 insertions(+), 11 deletions(-) diff --git a/python/damask/_configmaterial.py b/python/damask/_configmaterial.py index 17d8c796d..ed72aa27d 100644 --- a/python/damask/_configmaterial.py +++ b/python/damask/_configmaterial.py @@ -9,12 +9,15 @@ from . import Orientation class ConfigMaterial(Config): """Material configuration.""" - def __init__(self): - """Initialize object with all required dictionary keys.""" - super().__init__() - self['material'] = [] - self['homogenization'] = {} - self['phase'] = {} + _defaults = {'material': [], + 'homogenization': {}, + 'phase': {}} + + def __init__(self,d={}): + """Initialize object with default dictionary keys.""" + super().__init__(d) + for k,v in self._defaults.items(): + if k not in self: self[k] = v def save(self,fname='material.yaml',**kwargs): """ diff --git a/python/tests/reference/ConfigMaterial/material.yaml b/python/tests/reference/ConfigMaterial/material.yaml index 933e295b3..fbba6a631 100644 --- a/python/tests/reference/ConfigMaterial/material.yaml +++ b/python/tests/reference/ConfigMaterial/material.yaml @@ -1,10 +1,10 @@ homogenization: SX: - N_constituents: 2 - mech: {type: none} + N_constituents: 1 + mechanics: {type: none} Taylor: N_constituents: 2 - mech: {type: isostrain} + mechanics: {type: isostrain} material: - constituents: @@ -34,11 +34,11 @@ material: phase: Aluminum: lattice: cF - mech: + mechanics: output: [F, P, F_e, F_p, L_p] elasticity: {C_11: 106.75e9, C_12: 60.41e9, C_44: 28.34e9, type: hooke} Steel: lattice: cI - mech: + mechanics: output: [F, P, F_e, F_p, L_p] elasticity: {C_11: 233.3e9, C_12: 135.5e9, C_44: 118.0e9, type: hooke} From 5fb0e4908b9b553b33a2e5858eb2fa4afcc2709e Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 18 Dec 2020 07:09:05 +0100 Subject: [PATCH 03/19] Examples reflect actual behavior --- python/damask/_configmaterial.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/python/damask/_configmaterial.py b/python/damask/_configmaterial.py index ed72aa27d..b94e9897a 100644 --- a/python/damask/_configmaterial.py +++ b/python/damask/_configmaterial.py @@ -85,6 +85,8 @@ class ConfigMaterial(Config): fraction: 1.0 phase: Steel homogenization: SX + homogenization: {} + phase: {} """ constituents_ = {k:table.get(v) for k,v in constituents.items()} @@ -271,6 +273,8 @@ class ConfigMaterial(Config): fraction: 1.0 phase: Aluminum homogenization: SX + homogenization: {} + phase: {} """ length = -1 From 3010d11c8e919ff4fa62cf277270d45e39f62732 Mon Sep 17 00:00:00 2001 From: Test User Date: Fri, 18 Dec 2020 10:51:46 +0100 Subject: [PATCH 04/19] [skip ci] updated version information after successful test of v3.0.0-alpha2-39-g5fb0e4908 --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index b7128c9af..d50c5b204 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v3.0.0-alpha2-35-g1ebd10745 +v3.0.0-alpha2-39-g5fb0e4908 From 35f9861818bbb3b3e6fdfd1d6f520c5fddc531c6 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 18 Dec 2020 15:19:04 +0100 Subject: [PATCH 05/19] Fortran standard is 2018 will not work for older compilers --- CMakeLists.txt | 76 +++------ PRIVATE | 2 +- cmake/Compiler-Intel.cmake | 2 +- python/damask/_config.py | 3 + python/tests/test_Table.py | 2 +- src/DAMASK_interface.f90 | 8 +- src/grid/grid_damage_spectral.f90 | 3 +- src/grid/grid_mech_FEM.f90 | 9 +- src/grid/grid_mech_spectral_basic.f90 | 9 +- src/grid/grid_mech_spectral_polarisation.f90 | 9 +- src/grid/grid_thermal_spectral.f90 | 3 +- src/grid/spectral_utilities.f90 | 23 ++- src/material.f90 | 2 +- src/mesh/mesh_mech_FEM.f90 | 33 +--- src/results.f90 | 169 ++++++++++--------- 15 files changed, 147 insertions(+), 206 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index dd2348fd1..8db6dd0c0 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,6 +1,18 @@ -######################################################################################## -# Compiler options for building DAMASK -cmake_minimum_required (VERSION 3.10.0 FATAL_ERROR) +cmake_minimum_required (VERSION 3.10.0) +include (FindPkgConfig REQUIRED) + +# Dummy project to determine compiler names and version +project (Prerequisites LANGUAGES) +set(ENV{PKG_CONFIG_PATH} "$ENV{PETSC_DIR}/$ENV{PETSC_ARCH}/lib/pkgconfig") +pkg_search_module (PETSC REQUIRED PETSc>3.12.0) +pkg_get_variable (CMAKE_Fortran_COMPILER PETSc fcompiler) +pkg_get_variable (CMAKE_C_COMPILER PETSc ccompiler) + +find_program (CAT_EXECUTABLE NAMES cat) +execute_process (COMMAND ${CAT_EXECUTABLE} ${PROJECT_SOURCE_DIR}/VERSION + RESULT_VARIABLE DAMASK_VERSION_RETURN + OUTPUT_VARIABLE DAMASK_VERSION + OUTPUT_STRIP_TRAILING_WHITESPACE) #--------------------------------------------------------------------------------------- # Find PETSc from system environment @@ -28,19 +40,10 @@ include ${petsc_conf_rules} include ${petsc_conf_variables} INCLUDE_DIRS := \${PETSC_FC_INCLUDES} LIBRARIES := \${PETSC_WITH_EXTERNAL_LIB} -COMPILERF := \${FC} -COMPILERC := \${CC} -LINKERNAME := \${FLINKER} includes: \t@echo \${INCLUDE_DIRS} extlibs: \t@echo \${LIBRARIES} -compilerf: -\t@echo \${COMPILERF} -compilerc: -\t@echo \${COMPILERC} -linker: -\t@echo \${LINKERNAME} ") # CMake will execute each target in the ${petsc_config_makefile} @@ -52,26 +55,10 @@ execute_process (COMMAND ${MAKE_EXECUTABLE} --no-print-directory -f ${petsc_conf OUTPUT_VARIABLE petsc_includes OUTPUT_STRIP_TRAILING_WHITESPACE) # Find the PETSc external linking directory settings -# required for final linking, must be appended after the executable execute_process (COMMAND ${MAKE_EXECUTABLE} --no-print-directory -f ${petsc_config_makefile} "extlibs" RESULT_VARIABLE PETSC_EXTERNAL_LIB_RETURN OUTPUT_VARIABLE petsc_external_lib OUTPUT_STRIP_TRAILING_WHITESPACE) -# PETSc specified fortran compiler -execute_process (COMMAND ${MAKE_EXECUTABLE} --no-print-directory -f ${petsc_config_makefile} "compilerf" - RESULT_VARIABLE PETSC_MPIFC_RETURN - OUTPUT_VARIABLE PETSC_MPIFC - OUTPUT_STRIP_TRAILING_WHITESPACE) -# PETSc specified C compiler -execute_process (COMMAND ${MAKE_EXECUTABLE} --no-print-directory -f ${petsc_config_makefile} "compilerc" - RESULT_VARIABLE PETSC_MPICC_RETURN - OUTPUT_VARIABLE PETSC_MPICC - OUTPUT_STRIP_TRAILING_WHITESPACE) -# PETSc specified linker (Fortran compiler + PETSc linking flags) -execute_process (COMMAND ${MAKE_EXECUTABLE} --no-print-directory -f ${petsc_config_makefile} "linker" - RESULT_VARIABLE PETSC_LINKER_RETURN - OUTPUT_VARIABLE PETSC_LINKER - OUTPUT_STRIP_TRAILING_WHITESPACE) # Remove temporary makefile, no need to keep it anymore. file (REMOVE_RECURSE ${TEMPDIR}) @@ -90,14 +77,6 @@ endforeach (exlib) message ("Found PETSC_DIR:\n${PETSC_DIR}\n" ) message ("Found PETSC_INCLUDES:\n${PETSC_INCLUDES}\n" ) message ("Found PETSC_EXTERNAL_LIB:\n${PETSC_EXTERNAL_LIB}\n") -message ("Found PETSC_LINKER:\n${PETSC_LINKER}\n" ) -message ("Found MPI Fortran Compiler:\n${PETSC_MPIFC}\n" ) -message ("Found MPI C Compiler:\n${PETSC_MPICC}\n" ) - -# set compiler commands to match PETSc (needs to be done before defining the project) -# https://cmake.org/Wiki/CMake_FAQ#How_do_I_use_a_different_compiler.3F -set (CMAKE_Fortran_COMPILER "${PETSC_MPIFC}") -set (CMAKE_C_COMPILER "${PETSC_MPICC}") #--------------------------------------------------------------------------------------- # Now start to care about DAMASK @@ -105,17 +84,18 @@ set (CMAKE_C_COMPILER "${PETSC_MPICC}") # DAMASK solver defines project to build string(TOLOWER ${DAMASK_SOLVER} DAMASK_SOLVER) if (DAMASK_SOLVER STREQUAL "grid") - project (damask-grid Fortran C) + project (damask-grid HOMEPAGE_URL https://damask.mpie.de LANGUAGES Fortran C) add_definitions (-DGrid) - message ("Building Grid Solver\n") elseif (DAMASK_SOLVER STREQUAL "mesh") - project (damask-mesh Fortran C) + project (damask-mesh HOMEPAGE_URL https://damask.mpie.de LANGUAGES Fortran C) add_definitions (-DMesh) - message ("Building Mesh Solver\n") else () message (FATAL_ERROR "Build target (DAMASK_SOLVER) is not defined") endif () -list(APPEND CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/cmake) +add_definitions (-DDAMASKVERSION="${DAMASK_VERSION}") +add_definitions (-DPETSc) + +message ("\nBuilding ${CMAKE_PROJECT_NAME}\n") if (CMAKE_BUILD_TYPE STREQUAL "") set (CMAKE_BUILD_TYPE "RELEASE") @@ -153,17 +133,8 @@ if (CMAKE_BUILD_TYPE STREQUAL "SYNTAXONLY") set (BUILDCMD_POST "${BUILDCMD_POST} -fsyntax-only") endif () -# Parse DAMASK version from VERSION file -find_program (CAT_EXECUTABLE NAMES cat) -execute_process (COMMAND ${CAT_EXECUTABLE} ${PROJECT_SOURCE_DIR}/VERSION - RESULT_VARIABLE DAMASK_VERSION_RETURN - OUTPUT_VARIABLE DAMASK_V - OUTPUT_STRIP_TRAILING_WHITESPACE) -add_definitions (-DDAMASKVERSION="${DAMASK_V}") - -# definition of other macros -add_definitions (-DPETSc) +list(APPEND CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/cmake) if (CMAKE_Fortran_COMPILER_ID STREQUAL "Intel") include (Compiler-Intel) elseif(CMAKE_Fortran_COMPILER_ID STREQUAL "GNU") @@ -174,9 +145,8 @@ else () message (FATAL_ERROR "Compiler type (CMAKE_Fortran_COMPILER_ID) not recognized") endif () - set (CMAKE_Fortran_FLAGS_${CMAKE_BUILD_TYPE} "${BUILDCMD_PRE} ${OPENMP_FLAGS} ${STANDARD_CHECK} ${OPTIMIZATION_FLAGS} ${COMPILE_FLAGS} ${PRECISION_FLAGS}") -set (CMAKE_Fortran_LINK_EXECUTABLE "${BUILDCMD_PRE} ${PETSC_LINKER} ${OPENMP_FLAGS} ${OPTIMIZATION_FLAGS} ${LINKER_FLAGS}") +set (CMAKE_Fortran_LINK_EXECUTABLE "${BUILDCMD_PRE} ${CMAKE_Fortran_COMPILER} ${OPENMP_FLAGS} ${OPTIMIZATION_FLAGS} ${LINKER_FLAGS}") if (CMAKE_BUILD_TYPE STREQUAL "DEBUG") set (CMAKE_Fortran_FLAGS_${CMAKE_BUILD_TYPE} "${CMAKE_Fortran_FLAGS_${CMAKE_BUILD_TYPE}} ${DEBUG_FLAGS}") diff --git a/PRIVATE b/PRIVATE index de65e1df5..e1a1048e1 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit de65e1df5a76362de93667e9820dbf330b56f96d +Subproject commit e1a1048e1f593683b4b432d41455bd236008c3ad diff --git a/cmake/Compiler-Intel.cmake b/cmake/Compiler-Intel.cmake index 719ed885b..5b551069e 100644 --- a/cmake/Compiler-Intel.cmake +++ b/cmake/Compiler-Intel.cmake @@ -20,7 +20,7 @@ endif () # -assume std_mod_proc_name (included in -standard-semantics) causes problems if other modules # (PETSc, HDF5) are not compiled with this option (https://software.intel.com/en-us/forums/intel-fortran-compiler-for-linux-and-mac-os-x/topic/62172) -set (STANDARD_CHECK "-stand f15 -standard-semantics -assume nostd_mod_proc_name") +set (STANDARD_CHECK "-stand f18 -standard-semantics -assume nostd_mod_proc_name") set (LINKER_FLAGS "${LINKER_FLAGS} -shared-intel") # Link against shared Intel libraries instead of static ones diff --git a/python/damask/_config.py b/python/damask/_config.py index 24245f4bd..76955588f 100644 --- a/python/damask/_config.py +++ b/python/damask/_config.py @@ -21,6 +21,9 @@ class NiceDumper(yaml.SafeDumper): return self.represent_data(dict(data)) if isinstance(data, dict) and type(data) != dict else \ super().represent_data(data) + def ignore_aliases(self, data): + """No references.""" + return True class Config(dict): """YAML-based configuration.""" diff --git a/python/tests/test_Table.py b/python/tests/test_Table.py index ac5859ecb..8f617aff5 100644 --- a/python/tests/test_Table.py +++ b/python/tests/test_Table.py @@ -22,7 +22,7 @@ class TestTable: @pytest.mark.parametrize('N',[10,40]) def test_len(self,N): - len(Table(np.random.rand(N,3),{'X':3})) == N + assert len(Table(np.random.rand(N,3),{'X':3})) == N def test_get_scalar(self,default): d = default.get('s') diff --git a/src/DAMASK_interface.f90 b/src/DAMASK_interface.f90 index 41f421eb8..d38020225 100644 --- a/src/DAMASK_interface.f90 +++ b/src/DAMASK_interface.f90 @@ -10,7 +10,7 @@ !> and working directory. !-------------------------------------------------------------------------------------------------- #define PETSC_MAJOR 3 -#define PETSC_MINOR_MIN 10 +#define PETSC_MINOR_MIN 12 #define PETSC_MINOR_MAX 14 module DAMASK_interface @@ -392,7 +392,7 @@ end function makeRelativePath subroutine catchSIGTERM(signal) bind(C) integer(C_INT), value :: signal - interface_SIGTERM = .true. + call interface_setSIGTERM(.true.) print'(a,i0,a)', ' received signal ',signal, ', set SIGTERM=TRUE' @@ -417,7 +417,7 @@ end subroutine interface_setSIGTERM subroutine catchSIGUSR1(signal) bind(C) integer(C_INT), value :: signal - interface_SIGUSR1 = .true. + call interface_setSIGUSR1(.true.) print'(a,i0,a)', ' received signal ',signal, ', set SIGUSR1=TRUE' @@ -442,7 +442,7 @@ end subroutine interface_setSIGUSR1 subroutine catchSIGUSR2(signal) bind(C) integer(C_INT), value :: signal - interface_SIGUSR2 = .true. + call interface_setSIGUSR2(.true.) print'(a,i0,a)', ' received signal ',signal, ', set SIGUSR2=TRUE' diff --git a/src/grid/grid_damage_spectral.f90 b/src/grid/grid_damage_spectral.f90 index 4c014f3c0..79437945b 100644 --- a/src/grid/grid_damage_spectral.f90 +++ b/src/grid/grid_damage_spectral.f90 @@ -203,8 +203,7 @@ function grid_damage_spectral_solution(timeinc) result(solution) call VecMax(solution_vec,devNull,phi_max,ierr); CHKERRQ(ierr) if (solution%converged) & print'(/,a)', ' ... nonlocal damage converged .....................................' - write(IO_STDOUT,'(/,a,f8.6,2x,f8.6,2x,e11.4,/)',advance='no') ' Minimum|Maximum|Delta Damage = ',& - phi_min, phi_max, stagNorm + print'(/,a,f8.6,2x,f8.6,2x,e11.4)', ' Minimum|Maximum|Delta Damage = ', phi_min, phi_max, stagNorm print'(/,a)', ' ===========================================================================' flush(IO_STDOUT) diff --git a/src/grid/grid_mech_FEM.f90 b/src/grid/grid_mech_FEM.f90 index 146f28567..741ce404a 100644 --- a/src/grid/grid_mech_FEM.f90 +++ b/src/grid/grid_mech_FEM.f90 @@ -509,11 +509,10 @@ subroutine formResidual(da_local,x_local, & newIteration: if (totalIter <= PETScIter) then totalIter = totalIter + 1 print'(1x,a,3(a,i0))', trim(incInfo), ' @ Iteration ', num%itmin, '≤',totalIter+1, '≤', num%itmax - if (debugRotation) & - write(IO_STDOUT,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & - ' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.)) - write(IO_STDOUT,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & - ' deformation gradient aim =', transpose(F_aim) + if (debugRotation) print'(/,a,/,2(3(f12.7,1x)/),3(f12.7,1x))', & + ' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.)) + print'(/,a,/,2(3(f12.7,1x)/),3(f12.7,1x))', & + ' deformation gradient aim =', transpose(F_aim) flush(IO_STDOUT) endif newIteration diff --git a/src/grid/grid_mech_spectral_basic.f90 b/src/grid/grid_mech_spectral_basic.f90 index 4f5ceff61..d87a22fb7 100644 --- a/src/grid/grid_mech_spectral_basic.f90 +++ b/src/grid/grid_mech_spectral_basic.f90 @@ -471,11 +471,10 @@ subroutine formResidual(in, F, & newIteration: if (totalIter <= PETScIter) then totalIter = totalIter + 1 print'(1x,a,3(a,i0))', trim(incInfo), ' @ Iteration ', num%itmin, '≤',totalIter, '≤', num%itmax - if(debugRotation) & - write(IO_STDOUT,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & - ' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.)) - write(IO_STDOUT,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & - ' deformation gradient aim =', transpose(F_aim) + if (debugRotation) print'(/,a,/,2(3(f12.7,1x)/),3(f12.7,1x))', & + ' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.)) + print'(/,a,/,2(3(f12.7,1x)/),3(f12.7,1x))', & + ' deformation gradient aim =', transpose(F_aim) flush(IO_STDOUT) endif newIteration diff --git a/src/grid/grid_mech_spectral_polarisation.f90 b/src/grid/grid_mech_spectral_polarisation.f90 index d09b7fcb2..9bbb40a53 100644 --- a/src/grid/grid_mech_spectral_polarisation.f90 +++ b/src/grid/grid_mech_spectral_polarisation.f90 @@ -552,11 +552,10 @@ subroutine formResidual(in, FandF_tau, & newIteration: if (totalIter <= PETScIter) then totalIter = totalIter + 1 print'(1x,a,3(a,i0))', trim(incInfo), ' @ Iteration ', num%itmin, '≤',totalIter, '≤', num%itmax - if (debugRotation) & - write(IO_STDOUT,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & - ' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.)) - write(IO_STDOUT,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & - ' deformation gradient aim =', transpose(F_aim) + if (debugRotation) print'(/,a,/,2(3(f12.7,1x)/),3(f12.7,1x))', & + ' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.)) + print'(/,a,/,2(3(f12.7,1x)/),3(f12.7,1x))', & + ' deformation gradient aim =', transpose(F_aim) flush(IO_STDOUT) endif newIteration diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90 index 68a1c5ed1..f5d1a33bc 100644 --- a/src/grid/grid_thermal_spectral.f90 +++ b/src/grid/grid_thermal_spectral.f90 @@ -197,8 +197,7 @@ function grid_thermal_spectral_solution(timeinc) result(solution) call VecMax(solution_vec,devNull,T_max,ierr); CHKERRQ(ierr) if (solution%converged) & print'(/,a)', ' ... thermal conduction converged ..................................' - write(IO_STDOUT,'(/,a,f8.4,2x,f8.4,2x,f8.4,/)',advance='no') ' Minimum|Maximum|Delta Temperature / K = ',& - T_min, T_max, stagNorm + print'(/,a,f8.4,2x,f8.4,2x,f8.4)', ' Minimum|Maximum|Delta Temperature / K = ', T_min, T_max, stagNorm print'(/,a)', ' ===========================================================================' flush(IO_STDOUT) diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index 6d6e26cae..989448dc3 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -688,8 +688,8 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C) if(debugGeneral) then print'(/,a)', ' ... updating masked compliance ............................................' - write(IO_STDOUT,'(/,a,/,9(9(2x,f12.7,1x)/))',advance='no') ' Stiffness C (load) / GPa =',& - transpose(temp99_Real)*1.0e-9_pReal + print'(/,a,/,8(9(2x,f12.7,1x)/),9(2x,f12.7,1x))', & + ' Stiffness C (load) / GPa =', transpose(temp99_Real)*1.0e-9_pReal flush(IO_STDOUT) endif @@ -709,9 +709,8 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C) if (debugGeneral .or. errmatinv) then write(formatString, '(i2)') size_reduced formatString = '(/,a,/,'//trim(formatString)//'('//trim(formatString)//'(2x,es9.2,1x)/))' - write(IO_STDOUT,trim(formatString),advance='no') ' C * S (load) ', & - transpose(matmul(c_reduced,s_reduced)) - write(IO_STDOUT,trim(formatString),advance='no') ' S (load) ', transpose(s_reduced) + print trim(formatString), ' C * S (load) ', transpose(matmul(c_reduced,s_reduced)) + print trim(formatString), ' S (load) ', transpose(s_reduced) if(errmatinv) error stop 'matrix inversion error' endif temp99_real = reshape(unpack(reshape(s_reduced,[size_reduced**2]),reshape(mask,[81]),0.0_pReal),[9,9]) @@ -722,7 +721,7 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C) utilities_maskedCompliance = math_99to3333(temp99_Real) if(debugGeneral) then - write(IO_STDOUT,'(/,a,/,9(9(2x,f10.5,1x)/),/)',advance='no') & + print'(/,a,/,9(9(2x,f10.5,1x)/),9(2x,f10.5,1x))', & ' Masked Compliance (load) * GPa =', transpose(temp99_Real)*1.0e9_pReal flush(IO_STDOUT) endif @@ -818,13 +817,11 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,& P = reshape(homogenization_P, [3,3,grid(1),grid(2),grid3]) P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt ! average of P call MPI_Allreduce(MPI_IN_PLACE,P_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) - if (debugRotation) & - write(IO_STDOUT,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress (lab) / MPa =',& - transpose(P_av)*1.e-6_pReal - if(present(rotation_BC)) & - P_av = rotation_BC%rotate(P_av) - write(IO_STDOUT,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress / MPa =',& - transpose(P_av)*1.e-6_pReal + if (debugRotation) print'(/,a,/,2(3(2x,f12.4,1x)/),3(2x,f12.4,1x))', & + ' Piola--Kirchhoff stress (lab) / MPa =', transpose(P_av)*1.e-6_pReal + if(present(rotation_BC)) P_av = rotation_BC%rotate(P_av) + print'(/,a,/,2(3(2x,f12.4,1x)/),3(2x,f12.4,1x))', & + ' Piola--Kirchhoff stress / MPa =', transpose(P_av)*1.e-6_pReal flush(IO_STDOUT) dPdF_max = 0.0_pReal diff --git a/src/material.f90 b/src/material.f90 index 223ea6ed8..b05979298 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -176,7 +176,7 @@ subroutine material_init(restart) if (.not. restart) then call results_openJobFile - call results_mapping_constituent(material_phaseAt,material_phaseMemberAt,material_name_phase) + call results_mapping_phase(material_phaseAt,material_phaseMemberAt,material_name_phase) call results_mapping_homogenization(material_homogenizationAt,material_homogenizationMemberAt,material_name_homogenization) call results_closeJobFile endif diff --git a/src/mesh/mesh_mech_FEM.f90 b/src/mesh/mesh_mech_FEM.f90 index b6ce1e175..eb5f862c2 100644 --- a/src/mesh/mesh_mech_FEM.f90 +++ b/src/mesh/mesh_mech_FEM.f90 @@ -146,14 +146,9 @@ subroutine FEM_mech_init(fieldBC) call PetscFESetQuadrature(mechFE,mechQuad,ierr); CHKERRQ(ierr) call PetscFEGetDimension(mechFE,nBasis,ierr); CHKERRQ(ierr) nBasis = nBasis/nc -#if (PETSC_VERSION_MINOR > 10) call DMAddField(mech_mesh,PETSC_NULL_DMLABEL,mechFE,ierr); CHKERRQ(ierr) call DMCreateDS(mech_mesh,ierr); CHKERRQ(ierr) -#endif call DMGetDS(mech_mesh,mechDS,ierr); CHKERRQ(ierr) -#if (PETSC_VERSION_MINOR < 11) - call PetscDSAddDiscretization(mechDS,mechFE,ierr); CHKERRQ(ierr) -#endif call PetscDSGetTotalDimension(mechDS,cellDof,ierr); CHKERRQ(ierr) call PetscFEDestroy(mechFE,ierr); CHKERRQ(ierr) call PetscQuadratureDestroy(mechQuad,ierr); CHKERRQ(ierr) @@ -162,11 +157,7 @@ subroutine FEM_mech_init(fieldBC) ! Setup FEM mech boundary conditions call DMGetLabel(mech_mesh,'Face Sets',BCLabel,ierr); CHKERRQ(ierr) call DMPlexLabelComplete(mech_mesh,BCLabel,ierr); CHKERRQ(ierr) -#if (PETSC_VERSION_MINOR < 12) - call DMGetSection(mech_mesh,section,ierr); CHKERRQ(ierr) -#else call DMGetLocalSection(mech_mesh,section,ierr); CHKERRQ(ierr) -#endif allocate(pnumComp(1), source=dimPlex) allocate(pnumDof(0:dimPlex), source = 0) do topologDim = 0, dimPlex @@ -204,14 +195,8 @@ subroutine FEM_mech_init(fieldBC) endif endif enddo; enddo -#if (PETSC_VERSION_MINOR < 11) - call DMPlexCreateSection(mech_mesh,dimPlex,1,pNumComp,pNumDof, & - numBC,pBcField,pBcComps,pBcPoints,PETSC_NULL_IS,section,ierr) -#else call DMPlexCreateSection(mech_mesh,nolabel,pNumComp,pNumDof, & numBC,pBcField,pBcComps,pBcPoints,PETSC_NULL_IS,section,ierr) - -#endif CHKERRQ(ierr) call DMSetSection(mech_mesh,section,ierr); CHKERRQ(ierr) do faceSet = 1, numBC @@ -266,11 +251,7 @@ subroutine FEM_mech_init(fieldBC) x_scal(basis+1:basis+dimPlex) = pV0 + matmul(transpose(cellJMat),nodalPointsP + 1.0_pReal) enddo px_scal => x_scal -#if (PETSC_VERSION_MINOR < 11) - call DMPlexVecSetClosure(mech_mesh,section,solution_local,cell,px_scal,INSERT_ALL_VALUES,ierr) -#else - call DMPlexVecSetClosure(mech_mesh,section,solution_local,cell,px_scal,5,ierr) ! PETSc: cbee0a90b60958e5c50c89b1e41f4451dfa6008c -#endif + call DMPlexVecSetClosure(mech_mesh,section,solution_local,cell,px_scal,5,ierr) CHKERRQ(ierr) enddo @@ -353,11 +334,7 @@ subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr) allocate(pinvcellJ(dimPlex**2)) allocate(x_scal(cellDof)) -#if (PETSC_VERSION_MINOR < 12) - call DMGetSection(dm_local,section,ierr); CHKERRQ(ierr) -#else call DMGetLocalSection(dm_local,section,ierr); CHKERRQ(ierr) -#endif call DMGetDS(dm_local,prob,ierr); CHKERRQ(ierr) call PetscDSGetTabulation(prob,0,basisField,basisFieldDer,ierr) CHKERRQ(ierr) @@ -500,11 +477,7 @@ subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr) call MatZeroEntries(Jac,ierr); CHKERRQ(ierr) call DMGetDS(dm_local,prob,ierr); CHKERRQ(ierr) call PetscDSGetTabulation(prob,0,basisField,basisFieldDer,ierr) -#if (PETSC_VERSION_MINOR < 12) - call DMGetSection(dm_local,section,ierr); CHKERRQ(ierr) -#else call DMGetLocalSection(dm_local,section,ierr); CHKERRQ(ierr) -#endif call DMGetGlobalSection(dm_local,gSection,ierr); CHKERRQ(ierr) call DMGetLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr) @@ -684,8 +657,8 @@ subroutine FEM_mech_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dumm print'(/,1x,a,a,i0,a,i0,f0.3)', trim(incInfo), & ' @ Iteration ',PETScIter,' mechanical residual norm = ', & int(fnorm/divTol),fnorm/divTol-int(fnorm/divTol) - write(IO_STDOUT,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress / MPa =',& - transpose(P_av)*1.e-6_pReal + print'(/,a,/,2(3(2x,f12.4,1x)/),3(2x,f12.4,1x))', & + ' Piola--Kirchhoff stress / MPa =',transpose(P_av)*1.e-6_pReal flush(IO_STDOUT) end subroutine FEM_mech_converged diff --git a/src/results.f90 b/src/results.f90 index f15ad4e4a..ea9fd62d4 100644 --- a/src/results.f90 +++ b/src/results.f90 @@ -49,7 +49,7 @@ module results results_setLink, & results_addAttribute, & results_removeLink, & - results_mapping_constituent, & + results_mapping_phase, & results_mapping_homogenization contains @@ -461,7 +461,7 @@ end subroutine results_writeTensorDataset_int !-------------------------------------------------------------------------------------------------- !> @brief adds the unique mapping from spatial position and constituent ID to results !-------------------------------------------------------------------------------------------------- -subroutine results_mapping_constituent(phaseAt,memberAtLocal,label) +subroutine results_mapping_phase(phaseAt,memberAtLocal,label) integer, dimension(:,:), intent(in) :: phaseAt !< phase section at (constituent,element) integer, dimension(:,:,:), intent(in) :: memberAtLocal !< phase member at (constituent,IP,element) @@ -491,6 +491,47 @@ subroutine results_mapping_constituent(phaseAt,memberAtLocal,label) integer(SIZE_T) :: type_size_string, type_size_int integer :: hdferr, ierr, i +!-------------------------------------------------------------------------------------------------- +! prepare MPI communication (transparent for non-MPI runs) + call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr) + if(hdferr < 0) error stop 'HDF5 error' + memberOffset = 0 + do i=1, size(label) + memberOffset(i,worldrank) = count(phaseAt == i)*size(memberAtLocal,2) ! number of points/instance of this process + enddo + writeSize = 0 + writeSize(worldrank) = size(memberAtLocal(1,:,:)) ! total number of points by this process + +!-------------------------------------------------------------------------------------------------- +! MPI settings and communication +#ifdef PETSc + call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, hdferr) + if(hdferr < 0) error stop 'HDF5 error' + + call MPI_allreduce(MPI_IN_PLACE,writeSize,worldsize,MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr) ! get output at each process + if(ierr /= 0) error stop 'MPI error' + + call MPI_allreduce(MPI_IN_PLACE,memberOffset,size(memberOffset),MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr)! get offset at each process + if(ierr /= 0) error stop 'MPI error' +#endif + + myShape = int([size(phaseAt,1),writeSize(worldrank)], HSIZE_T) + myOffset = int([0,sum(writeSize(0:worldrank-1))], HSIZE_T) + totalShape = int([size(phaseAt,1),sum(writeSize)], HSIZE_T) + + +!--------------------------------------------------------------------------------------------------- +! expand phaseAt to consider IPs (is not stored per IP) + do i = 1, size(phaseAtMaterialpoint,2) + phaseAtMaterialpoint(:,i,:) = phaseAt + enddo + +!--------------------------------------------------------------------------------------------------- +! renumber member from my process to all processes + do i = 1, size(label) + where(phaseAtMaterialpoint == i) memberAtGlobal = memberAtLocal + sum(memberOffset(i,0:worldrank-1)) -1 ! convert to 0-based + enddo + !--------------------------------------------------------------------------------------------------- ! compound type: name of phase section + position/index within results array call h5tcopy_f(H5T_NATIVE_CHARACTER, dt_id, hdferr) @@ -525,34 +566,6 @@ subroutine results_mapping_constituent(phaseAt,memberAtLocal,label) call h5tclose_f(dt_id, hdferr) if(hdferr < 0) error stop 'HDF5 error' -!-------------------------------------------------------------------------------------------------- -! prepare MPI communication (transparent for non-MPI runs) - call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr) - if(hdferr < 0) error stop 'HDF5 error' - memberOffset = 0 - do i=1, size(label) - memberOffset(i,worldrank) = count(phaseAt == i)*size(memberAtLocal,2) ! number of points/instance of this process - enddo - writeSize = 0 - writeSize(worldrank) = size(memberAtLocal(1,:,:)) ! total number of points by this process - -!-------------------------------------------------------------------------------------------------- -! MPI settings and communication -#ifdef PETSc - call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, hdferr) - if(hdferr < 0) error stop 'HDF5 error' - - call MPI_allreduce(MPI_IN_PLACE,writeSize,worldsize,MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr) ! get output at each process - if(ierr /= 0) error stop 'MPI error' - - call MPI_allreduce(MPI_IN_PLACE,memberOffset,size(memberOffset),MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr)! get offset at each process - if(ierr /= 0) error stop 'MPI error' -#endif - - myShape = int([size(phaseAt,1),writeSize(worldrank)], HSIZE_T) - myOffset = int([0,sum(writeSize(0:worldrank-1))], HSIZE_T) - totalShape = int([size(phaseAt,1),sum(writeSize)], HSIZE_T) - !-------------------------------------------------------------------------------------------------- ! create dataspace in memory (local shape = hyperslab) and in file (global shape) call h5screate_simple_f(2,myShape,memspace_id,hdferr,myShape) @@ -564,18 +577,6 @@ subroutine results_mapping_constituent(phaseAt,memberAtLocal,label) call h5sselect_hyperslab_f(filespace_id, H5S_SELECT_SET_F, myOffset, myShape, hdferr) if(hdferr < 0) error stop 'HDF5 error' -!--------------------------------------------------------------------------------------------------- -! expand phaseAt to consider IPs (is not stored per IP) - do i = 1, size(phaseAtMaterialpoint,2) - phaseAtMaterialpoint(:,i,:) = phaseAt - enddo - -!--------------------------------------------------------------------------------------------------- -! renumber member from my process to all processes - do i = 1, size(label) - where(phaseAtMaterialpoint == i) memberAtGlobal = memberAtLocal + sum(memberOffset(i,0:worldrank-1)) -1 ! convert to 0-based - enddo - !-------------------------------------------------------------------------------------------------- ! write the components of the compound type individually call h5pset_preserve_f(plist_id, .TRUE., hdferr) @@ -609,7 +610,7 @@ subroutine results_mapping_constituent(phaseAt,memberAtLocal,label) if(hdferr < 0) error stop 'HDF5 error' call h5tclose_f(position_id, hdferr) -end subroutine results_mapping_constituent +end subroutine results_mapping_phase !-------------------------------------------------------------------------------------------------- @@ -645,6 +646,48 @@ subroutine results_mapping_homogenization(homogenizationAt,memberAtLocal,label) integer(SIZE_T) :: type_size_string, type_size_int integer :: hdferr, ierr, i + +!-------------------------------------------------------------------------------------------------- +! prepare MPI communication (transparent for non-MPI runs) + call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr) + if(hdferr < 0) error stop 'HDF5 error' + memberOffset = 0 + do i=1, size(label) + memberOffset(i,worldrank) = count(homogenizationAt == i)*size(memberAtLocal,1) ! number of points/instance of this process + enddo + writeSize = 0 + writeSize(worldrank) = size(memberAtLocal) ! total number of points by this process + +!-------------------------------------------------------------------------------------------------- +! MPI settings and communication +#ifdef PETSc + call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, hdferr) + if(hdferr < 0) error stop 'HDF5 error' + + call MPI_allreduce(MPI_IN_PLACE,writeSize,worldsize,MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr) ! get output at each process + if(ierr /= 0) error stop 'MPI error' + + call MPI_allreduce(MPI_IN_PLACE,memberOffset,size(memberOffset),MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr)! get offset at each process + if(ierr /= 0) error stop 'MPI error' +#endif + + myShape = int([writeSize(worldrank)], HSIZE_T) + myOffset = int([sum(writeSize(0:worldrank-1))], HSIZE_T) + totalShape = int([sum(writeSize)], HSIZE_T) + + +!--------------------------------------------------------------------------------------------------- +! expand phaseAt to consider IPs (is not stored per IP) + do i = 1, size(homogenizationAtMaterialpoint,1) + homogenizationAtMaterialpoint(i,:) = homogenizationAt + enddo + +!--------------------------------------------------------------------------------------------------- +! renumber member from my process to all processes + do i = 1, size(label) + where(homogenizationAtMaterialpoint == i) memberAtGlobal = memberAtLocal + sum(memberOffset(i,0:worldrank-1)) - 1 ! convert to 0-based + enddo + !--------------------------------------------------------------------------------------------------- ! compound type: name of phase section + position/index within results array call h5tcopy_f(H5T_NATIVE_CHARACTER, dt_id, hdferr) @@ -679,34 +722,6 @@ subroutine results_mapping_homogenization(homogenizationAt,memberAtLocal,label) call h5tclose_f(dt_id, hdferr) if(hdferr < 0) error stop 'HDF5 error' -!-------------------------------------------------------------------------------------------------- -! prepare MPI communication (transparent for non-MPI runs) - call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr) - if(hdferr < 0) error stop 'HDF5 error' - memberOffset = 0 - do i=1, size(label) - memberOffset(i,worldrank) = count(homogenizationAt == i)*size(memberAtLocal,1) ! number of points/instance of this process - enddo - writeSize = 0 - writeSize(worldrank) = size(memberAtLocal) ! total number of points by this process - -!-------------------------------------------------------------------------------------------------- -! MPI settings and communication -#ifdef PETSc - call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, hdferr) - if(hdferr < 0) error stop 'HDF5 error' - - call MPI_allreduce(MPI_IN_PLACE,writeSize,worldsize,MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr) ! get output at each process - if(ierr /= 0) error stop 'MPI error' - - call MPI_allreduce(MPI_IN_PLACE,memberOffset,size(memberOffset),MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr)! get offset at each process - if(ierr /= 0) error stop 'MPI error' -#endif - - myShape = int([writeSize(worldrank)], HSIZE_T) - myOffset = int([sum(writeSize(0:worldrank-1))], HSIZE_T) - totalShape = int([sum(writeSize)], HSIZE_T) - !-------------------------------------------------------------------------------------------------- ! create dataspace in memory (local shape = hyperslab) and in file (global shape) call h5screate_simple_f(1,myShape,memspace_id,hdferr,myShape) @@ -718,18 +733,6 @@ subroutine results_mapping_homogenization(homogenizationAt,memberAtLocal,label) call h5sselect_hyperslab_f(filespace_id, H5S_SELECT_SET_F, myOffset, myShape, hdferr) if(hdferr < 0) error stop 'HDF5 error' -!--------------------------------------------------------------------------------------------------- -! expand phaseAt to consider IPs (is not stored per IP) - do i = 1, size(homogenizationAtMaterialpoint,1) - homogenizationAtMaterialpoint(i,:) = homogenizationAt - enddo - -!--------------------------------------------------------------------------------------------------- -! renumber member from my process to all processes - do i = 1, size(label) - where(homogenizationAtMaterialpoint == i) memberAtGlobal = memberAtLocal + sum(memberOffset(i,0:worldrank-1)) - 1 ! convert to 0-based - enddo - !-------------------------------------------------------------------------------------------------- ! write the components of the compound type individually call h5pset_preserve_f(plist_id, .TRUE., hdferr) From 35b833e2adddf19b96c616a54ae217e49b16f2cb Mon Sep 17 00:00:00 2001 From: Test User Date: Fri, 18 Dec 2020 20:08:35 +0100 Subject: [PATCH 06/19] [skip ci] updated version information after successful test of v3.0.0-alpha2-42-g6cc78cb41 --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index d50c5b204..85027fcd3 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v3.0.0-alpha2-39-g5fb0e4908 +v3.0.0-alpha2-42-g6cc78cb41 From c913a577d01e0b2f5887a9748fc1cd459426f94c Mon Sep 17 00:00:00 2001 From: Test User Date: Mon, 21 Dec 2020 19:15:06 +0100 Subject: [PATCH 07/19] [skip ci] updated version information after successful test of v3.0.0-alpha2-62-gd1c7ec068 --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 85027fcd3..5c2eab7c5 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v3.0.0-alpha2-42-g6cc78cb41 +v3.0.0-alpha2-62-gd1c7ec068 From 3ea05cabadb680c16e0dead27f5f30142bd58e24 Mon Sep 17 00:00:00 2001 From: Test User Date: Mon, 21 Dec 2020 21:32:13 +0100 Subject: [PATCH 08/19] [skip ci] updated version information after successful test of v3.0.0-alpha2-68-gbc32797fa --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 85027fcd3..72db22a20 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v3.0.0-alpha2-42-g6cc78cb41 +v3.0.0-alpha2-68-gbc32797fa From d8b57680ec0c43520f045954c7edc01c81fcb7df Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Mon, 21 Dec 2020 15:46:07 -0500 Subject: [PATCH 09/19] raise NotImplemented when using R*b instead of R@b --- python/damask/_rotation.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index f4d2cf248..780e81891 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -144,6 +144,11 @@ class Rotation: return self.copy(rotation=Rotation(np.block([np.cos(pwr*phi),np.sin(pwr*phi)*p]))._standardize()) + def __mul__(self,other): + """Standard multiplication is not implemented.""" + raise NotImplementedError('Use "R@b", i.e. matmul, to apply rotation "R" to object "b"') + + def __matmul__(self,other): """ Rotation of vector, second or fourth order tensor, or rotation object. From 8b2f75b99b02a6f9844d4697890a9c2f59857da7 Mon Sep 17 00:00:00 2001 From: Test User Date: Tue, 22 Dec 2020 11:22:57 +0100 Subject: [PATCH 10/19] [skip ci] updated version information after successful test of v3.0.0-alpha2-75-gac45427e9 --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 5c2eab7c5..3e5a7aa8c 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v3.0.0-alpha2-62-gd1c7ec068 +v3.0.0-alpha2-75-gac45427e9 From 6d5c3a5d12587975c2da3eba9c95f787d5dd62de Mon Sep 17 00:00:00 2001 From: Test User Date: Thu, 24 Dec 2020 20:02:59 +0100 Subject: [PATCH 11/19] [skip ci] updated version information after successful test of v3.0.0-alpha2-97-g10bbeb561 --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 3e5a7aa8c..616042312 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v3.0.0-alpha2-75-gac45427e9 +v3.0.0-alpha2-97-g10bbeb561 From 4796afdd92f7b576c36bf62017fd8ba1b036899f Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Mon, 28 Dec 2020 12:10:21 -0500 Subject: [PATCH 12/19] fix for broken representation of no-rotation orientations and averaging weights --- python/damask/_orientation.py | 6 +++--- python/damask/_rotation.py | 14 +++++++++++--- 2 files changed, 14 insertions(+), 6 deletions(-) diff --git a/python/damask/_orientation.py b/python/damask/_orientation.py index 05301561f..e6434813e 100644 --- a/python/damask/_orientation.py +++ b/python/damask/_orientation.py @@ -226,9 +226,9 @@ class Orientation(Rotation): """ return super().__eq__(other) \ - and self.family == other.family \ - and self.lattice == other.lattice \ - and self.parameters == other.parameters + and hasattr(other, 'family') and self.family == other.family \ + and hasattr(other, 'lattice') and self.lattice == other.lattice \ + and hasattr(other, 'parameters') and self.parameters == other.parameters def __matmul__(self,other): diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index 780e81891..8d3050ab8 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -204,8 +204,16 @@ class Rotation: def append(self,other): - """Extend rotation array along first dimension with other array.""" - return self.copy(rotation=np.vstack((self.quaternion,other.quaternion))) + """ + Extend rotation array along first dimension with other array(s). + + Parameters + ---------- + other : Rotation or list of Rotations. + + """ + return self.copy(rotation=np.vstack(tuple(map(lambda x:x.quaternion, + [self]+other if type(other) == list else [self,other])))) def flatten(self,order = 'C'): @@ -263,7 +271,7 @@ class Rotation: """Intermediate representation supporting quaternion averaging.""" return np.einsum('...i,...j',quat,quat) - if not weights: + if weights is None: weights = np.ones(self.shape,dtype=float) eig, vec = np.linalg.eig(np.sum(_M(self.quaternion) * weights[...,np.newaxis,np.newaxis],axis=-3) \ From da62daf15d37b6a1d42a1b23d171b3f82e4f6120 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Mon, 28 Dec 2020 12:26:09 -0500 Subject: [PATCH 13/19] added test for appending rotation lists; better check for type==list --- python/damask/_rotation.py | 2 +- python/tests/test_Rotation.py | 8 ++++++++ 2 files changed, 9 insertions(+), 1 deletion(-) diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index 8d3050ab8..78029b59a 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -213,7 +213,7 @@ class Rotation: """ return self.copy(rotation=np.vstack(tuple(map(lambda x:x.quaternion, - [self]+other if type(other) == list else [self,other])))) + [self]+other if isinstance(other,list) else [self,other])))) def flatten(self,order = 'C'): diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py index c60029046..36e3a3ac9 100644 --- a/python/tests/test_Rotation.py +++ b/python/tests/test_Rotation.py @@ -800,6 +800,14 @@ class TestRotation: print(f'append 2x {shape} --> {s.shape}') assert s[0,...] == r[0,...] and s[-1,...] == p[-1,...] + @pytest.mark.parametrize('shape',[None,1,(1,),(4,2),(3,3,2)]) + def test_append_list(self,shape): + r = Rotation.from_random(shape=shape) + p = Rotation.from_random(shape=shape) + s = r.append([r,p]) + print(f'append 3x {shape} --> {s.shape}') + assert s[0,...] == r[0,...] and s[-1,...] == p[-1,...] + @pytest.mark.parametrize('quat,standardized',[ ([-1,0,0,0],[1,0,0,0]), ([-0.5,-0.5,-0.5,-0.5],[0.5,0.5,0.5,0.5]), From 2791f3d192da8179f3f980aa56bbcbd9ddb3beb6 Mon Sep 17 00:00:00 2001 From: Test User Date: Tue, 29 Dec 2020 00:33:12 +0100 Subject: [PATCH 14/19] [skip ci] updated version information after successful test of v3.0.0-alpha2-101-gab2a4a987 --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 616042312..cb10a9c4f 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v3.0.0-alpha2-97-g10bbeb561 +v3.0.0-alpha2-101-gab2a4a987 From 0f28d8048b1bc7316f3f381fef96b9a9c9bd0196 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 2 Jan 2021 17:57:11 +0100 Subject: [PATCH 15/19] KISS --- src/commercialFEM_fileList.f90 | 1 - src/quaternions.f90 | 534 --------------------------------- src/rotations.f90 | 83 +++-- 3 files changed, 57 insertions(+), 561 deletions(-) delete mode 100644 src/quaternions.f90 diff --git a/src/commercialFEM_fileList.f90 b/src/commercialFEM_fileList.f90 index 08e7b9c1c..3e3e017eb 100644 --- a/src/commercialFEM_fileList.f90 +++ b/src/commercialFEM_fileList.f90 @@ -11,7 +11,6 @@ #include "config.f90" #include "LAPACK_interface.f90" #include "math.f90" -#include "quaternions.f90" #include "rotations.f90" #include "FEsolving.f90" #include "element.f90" diff --git a/src/quaternions.f90 b/src/quaternions.f90 deleted file mode 100644 index c5c43e3c1..000000000 --- a/src/quaternions.f90 +++ /dev/null @@ -1,534 +0,0 @@ -!--------------------------------------------------------------------------------------------------- -!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH -!> @author Philip Eisenlohr, Michigan State University -!> @brief general quaternion math, not limited to unit quaternions -!> @details w is the real part, (x, y, z) are the imaginary parts. -!> @details https://en.wikipedia.org/wiki/Quaternion -!--------------------------------------------------------------------------------------------------- -module quaternions - use prec - - implicit none - private - - real(pReal), parameter, public :: P = -1.0_pReal !< parameter for orientation conversion. - - type, public :: quaternion - real(pReal), private :: w = 0.0_pReal - real(pReal), private :: x = 0.0_pReal - real(pReal), private :: y = 0.0_pReal - real(pReal), private :: z = 0.0_pReal - - - contains - procedure, private :: add__ - procedure, private :: pos__ - generic, public :: operator(+) => add__,pos__ - - procedure, private :: sub__ - procedure, private :: neg__ - generic, public :: operator(-) => sub__,neg__ - - procedure, private :: mul_quat__ - procedure, private :: mul_scal__ - generic, public :: operator(*) => mul_quat__, mul_scal__ - - procedure, private :: div_quat__ - procedure, private :: div_scal__ - generic, public :: operator(/) => div_quat__, div_scal__ - - procedure, private :: eq__ - generic, public :: operator(==) => eq__ - - procedure, private :: neq__ - generic, public :: operator(/=) => neq__ - - procedure, private :: pow_quat__ - procedure, private :: pow_scal__ - generic, public :: operator(**) => pow_quat__, pow_scal__ - - procedure, public :: abs => abs__ - procedure, public :: conjg => conjg__ - procedure, public :: real => real__ - procedure, public :: aimag => aimag__ - - procedure, public :: homomorphed - procedure, public :: asArray - procedure, public :: inverse - - end type - - interface assignment (=) - module procedure assign_quat__ - module procedure assign_vec__ - end interface assignment (=) - - interface quaternion - module procedure init__ - end interface quaternion - - interface abs - procedure abs__ - end interface abs - - interface dot_product - procedure dot_product__ - end interface dot_product - - interface conjg - module procedure conjg__ - end interface conjg - - interface exp - module procedure exp__ - end interface exp - - interface log - module procedure log__ - end interface log - - interface real - module procedure real__ - end interface real - - interface aimag - module procedure aimag__ - end interface aimag - - public :: & - quaternions_init, & - assignment(=), & - conjg, aimag, & - log, exp, & - abs, dot_product, & - inverse, & - real - -contains - - -!-------------------------------------------------------------------------------------------------- -!> @brief Do self test. -!-------------------------------------------------------------------------------------------------- -subroutine quaternions_init - - print'(/,a)', ' <<<+- quaternions init -+>>>'; flush(6) - - call selfTest - -end subroutine quaternions_init - - -!--------------------------------------------------------------------------------------------------- -!> @brief construct a quaternion from a 4-vector -!--------------------------------------------------------------------------------------------------- -type(quaternion) pure function init__(array) - - real(pReal), intent(in), dimension(4) :: array - - init__%w = array(1) - init__%x = array(2) - init__%y = array(3) - init__%z = array(4) - -end function init__ - - -!--------------------------------------------------------------------------------------------------- -!> @brief assign a quaternion -!--------------------------------------------------------------------------------------------------- -elemental pure subroutine assign_quat__(self,other) - - type(quaternion), intent(out) :: self - type(quaternion), intent(in) :: other - - self = [other%w,other%x,other%y,other%z] - -end subroutine assign_quat__ - - -!--------------------------------------------------------------------------------------------------- -!> @brief assign a 4-vector -!--------------------------------------------------------------------------------------------------- -pure subroutine assign_vec__(self,other) - - type(quaternion), intent(out) :: self - real(pReal), intent(in), dimension(4) :: other - - self%w = other(1) - self%x = other(2) - self%y = other(3) - self%z = other(4) - -end subroutine assign_vec__ - - -!--------------------------------------------------------------------------------------------------- -!> @brief add a quaternion -!--------------------------------------------------------------------------------------------------- -type(quaternion) elemental pure function add__(self,other) - - class(quaternion), intent(in) :: self,other - - add__ = [ self%w, self%x, self%y ,self%z] & - + [other%w, other%x, other%y,other%z] - -end function add__ - - -!--------------------------------------------------------------------------------------------------- -!> @brief return (unary positive operator) -!--------------------------------------------------------------------------------------------------- -type(quaternion) elemental pure function pos__(self) - - class(quaternion), intent(in) :: self - - pos__ = self * (+1.0_pReal) - -end function pos__ - - -!--------------------------------------------------------------------------------------------------- -!> @brief subtract a quaternion -!--------------------------------------------------------------------------------------------------- -type(quaternion) elemental pure function sub__(self,other) - - class(quaternion), intent(in) :: self,other - - sub__ = [ self%w, self%x, self%y ,self%z] & - - [other%w, other%x, other%y,other%z] - -end function sub__ - - -!--------------------------------------------------------------------------------------------------- -!> @brief negate (unary negative operator) -!--------------------------------------------------------------------------------------------------- -type(quaternion) elemental pure function neg__(self) - - class(quaternion), intent(in) :: self - - neg__ = self * (-1.0_pReal) - -end function neg__ - - -!--------------------------------------------------------------------------------------------------- -!> @brief multiply with a quaternion -!--------------------------------------------------------------------------------------------------- -type(quaternion) elemental pure function mul_quat__(self,other) - - class(quaternion), intent(in) :: self, other - - mul_quat__%w = self%w*other%w - self%x*other%x - self%y*other%y - self%z*other%z - mul_quat__%x = self%w*other%x + self%x*other%w + P * (self%y*other%z - self%z*other%y) - mul_quat__%y = self%w*other%y + self%y*other%w + P * (self%z*other%x - self%x*other%z) - mul_quat__%z = self%w*other%z + self%z*other%w + P * (self%x*other%y - self%y*other%x) - -end function mul_quat__ - - -!--------------------------------------------------------------------------------------------------- -!> @brief multiply with a scalar -!--------------------------------------------------------------------------------------------------- -type(quaternion) elemental pure function mul_scal__(self,scal) - - class(quaternion), intent(in) :: self - real(pReal), intent(in) :: scal - - mul_scal__ = [self%w,self%x,self%y,self%z]*scal - -end function mul_scal__ - - -!--------------------------------------------------------------------------------------------------- -!> @brief divide by a quaternion -!--------------------------------------------------------------------------------------------------- -type(quaternion) elemental pure function div_quat__(self,other) - - class(quaternion), intent(in) :: self, other - - div_quat__ = self * (conjg(other)/(abs(other)**2.0_pReal)) - -end function div_quat__ - - -!--------------------------------------------------------------------------------------------------- -!> @brief divide by a scalar -!--------------------------------------------------------------------------------------------------- -type(quaternion) elemental pure function div_scal__(self,scal) - - class(quaternion), intent(in) :: self - real(pReal), intent(in) :: scal - - div_scal__ = [self%w,self%x,self%y,self%z]/scal - -end function div_scal__ - - -!--------------------------------------------------------------------------------------------------- -!> @brief test equality -!--------------------------------------------------------------------------------------------------- -logical elemental pure function eq__(self,other) - - class(quaternion), intent(in) :: self,other - - eq__ = all(dEq([ self%w, self%x, self%y, self%z], & - [other%w,other%x,other%y,other%z])) - -end function eq__ - - -!--------------------------------------------------------------------------------------------------- -!> @brief test inequality -!--------------------------------------------------------------------------------------------------- -logical elemental pure function neq__(self,other) - - class(quaternion), intent(in) :: self,other - - neq__ = .not. self%eq__(other) - -end function neq__ - - -!--------------------------------------------------------------------------------------------------- -!> @brief raise to the power of a quaternion -!--------------------------------------------------------------------------------------------------- -type(quaternion) elemental pure function pow_quat__(self,expon) - - class(quaternion), intent(in) :: self - type(quaternion), intent(in) :: expon - - pow_quat__ = exp(log(self)*expon) - -end function pow_quat__ - - -!--------------------------------------------------------------------------------------------------- -!> @brief raise to the power of a scalar -!--------------------------------------------------------------------------------------------------- -type(quaternion) elemental pure function pow_scal__(self,expon) - - class(quaternion), intent(in) :: self - real(pReal), intent(in) :: expon - - pow_scal__ = exp(log(self)*expon) - -end function pow_scal__ - - -!--------------------------------------------------------------------------------------------------- -!> @brief take exponential -!--------------------------------------------------------------------------------------------------- -type(quaternion) elemental pure function exp__(a) - - class(quaternion), intent(in) :: a - real(pReal) :: absImag - - absImag = norm2(aimag(a)) - - exp__ = merge(exp(a%w) * [ cos(absImag), & - a%x/absImag * sin(absImag), & - a%y/absImag * sin(absImag), & - a%z/absImag * sin(absImag)], & - IEEE_value(1.0_pReal,IEEE_SIGNALING_NAN), & - dNeq0(absImag)) - -end function exp__ - - -!--------------------------------------------------------------------------------------------------- -!> @brief take logarithm -!--------------------------------------------------------------------------------------------------- -type(quaternion) elemental pure function log__(a) - - class(quaternion), intent(in) :: a - real(pReal) :: absImag - - absImag = norm2(aimag(a)) - - log__ = merge([log(abs(a)), & - a%x/absImag * acos(a%w/abs(a)), & - a%y/absImag * acos(a%w/abs(a)), & - a%z/absImag * acos(a%w/abs(a))], & - IEEE_value(1.0_pReal,IEEE_SIGNALING_NAN), & - dNeq0(absImag)) - -end function log__ - - -!--------------------------------------------------------------------------------------------------- -!> @brief return norm -!--------------------------------------------------------------------------------------------------- -real(pReal) elemental pure function abs__(self) - - class(quaternion), intent(in) :: self - - abs__ = norm2([self%w,self%x,self%y,self%z]) - -end function abs__ - - -!--------------------------------------------------------------------------------------------------- -!> @brief calculate dot product -!--------------------------------------------------------------------------------------------------- -real(pReal) elemental pure function dot_product__(a,b) - - class(quaternion), intent(in) :: a,b - - dot_product__ = a%w*b%w + a%x*b%x + a%y*b%y + a%z*b%z - -end function dot_product__ - - -!--------------------------------------------------------------------------------------------------- -!> @brief take conjugate complex -!--------------------------------------------------------------------------------------------------- -type(quaternion) elemental pure function conjg__(self) - - class(quaternion), intent(in) :: self - - conjg__ = [self%w,-self%x,-self%y,-self%z] - -end function conjg__ - - -!--------------------------------------------------------------------------------------------------- -!> @brief homomorph -!--------------------------------------------------------------------------------------------------- -type(quaternion) elemental pure function homomorphed(self) - - class(quaternion), intent(in) :: self - - homomorphed = - self - -end function homomorphed - - -!--------------------------------------------------------------------------------------------------- -!> @brief return as plain array -!--------------------------------------------------------------------------------------------------- -pure function asArray(self) - - real(pReal), dimension(4) :: asArray - class(quaternion), intent(in) :: self - - asArray = [self%w,self%x,self%y,self%z] - -end function asArray - - -!--------------------------------------------------------------------------------------------------- -!> @brief real part (scalar) -!--------------------------------------------------------------------------------------------------- -pure function real__(self) - - real(pReal) :: real__ - class(quaternion), intent(in) :: self - - real__ = self%w - -end function real__ - - -!--------------------------------------------------------------------------------------------------- -!> @brief imaginary part (3-vector) -!--------------------------------------------------------------------------------------------------- -pure function aimag__(self) - - real(pReal), dimension(3) :: aimag__ - class(quaternion), intent(in) :: self - - aimag__ = [self%x,self%y,self%z] - -end function aimag__ - - -!--------------------------------------------------------------------------------------------------- -!> @brief inverse -!--------------------------------------------------------------------------------------------------- -type(quaternion) elemental pure function inverse(self) - - class(quaternion), intent(in) :: self - - inverse = conjg(self)/abs(self)**2.0_pReal - -end function inverse - - -!-------------------------------------------------------------------------------------------------- -!> @brief check correctness of some quaternions functions -!-------------------------------------------------------------------------------------------------- -subroutine selfTest - - real(pReal), dimension(4) :: qu - type(quaternion) :: q, q_2 - - if(dNeq(abs(P),1.0_pReal)) error stop 'P not in {-1,+1}' - - call random_number(qu) - qu = (qu-0.5_pReal) * 2.0_pReal - q = quaternion(qu) - - q_2= qu - if(any(dNeq(q%asArray(),q_2%asArray()))) error stop 'assign_vec__' - - q_2 = q + q - if(any(dNeq(q_2%asArray(),2.0_pReal*qu))) error stop 'add__' - - q_2 = q - q - if(any(dNeq0(q_2%asArray()))) error stop 'sub__' - - q_2 = q * 5.0_pReal - if(any(dNeq(q_2%asArray(),5.0_pReal*qu))) error stop 'mul__' - - q_2 = q / 0.5_pReal - if(any(dNeq(q_2%asArray(),2.0_pReal*qu))) error stop 'div__' - - q_2 = q * 0.3_pReal - if(dNeq0(abs(q)) .and. q_2 == q) error stop 'eq__' - - q_2 = q - if(q_2 /= q) error stop 'neq__' - - if(dNeq(abs(q),norm2(qu))) error stop 'abs__' - if(dNeq(abs(q)**2.0_pReal, real(q*q%conjg()),1.0e-14_pReal)) & - error stop 'abs__/*conjg' - - if(any(dNeq(q%asArray(),qu))) error stop 'eq__' - if(dNeq(q%real(), qu(1))) error stop 'real()' - if(any(dNeq(q%aimag(), qu(2:4)))) error stop 'aimag()' - - q_2 = q%homomorphed() - if(q /= q_2* (-1.0_pReal)) error stop 'homomorphed' - if(dNeq(q_2%real(), qu(1)* (-1.0_pReal))) error stop 'homomorphed/real' - if(any(dNeq(q_2%aimag(),qu(2:4)*(-1.0_pReal)))) error stop 'homomorphed/aimag' - - q_2 = conjg(q) - if(dNeq(abs(q),abs(q_2))) error stop 'conjg/abs' - if(q /= conjg(q_2)) error stop 'conjg/involution' - if(dNeq(q_2%real(), q%real())) error stop 'conjg/real' - if(any(dNeq(q_2%aimag(),q%aimag()*(-1.0_pReal)))) error stop 'conjg/aimag' - - if(abs(q) > 0.0_pReal) then - q_2 = q * q%inverse() - if( dNeq(real(q_2), 1.0_pReal,1.0e-15_pReal)) error stop 'inverse/real' - if(any(dNeq0(aimag(q_2), 1.0e-15_pReal))) error stop 'inverse/aimag' - - q_2 = q/abs(q) - q_2 = conjg(q_2) - inverse(q_2) - if(any(dNeq0(q_2%asArray(),1.0e-15_pReal))) error stop 'inverse/conjg' - endif - if(dNeq(dot_product(qu,qu),dot_product(q,q))) error stop 'dot_product' - -#if !(defined(__GFORTRAN__) && __GNUC__ < 9) - if (norm2(aimag(q)) > 0.0_pReal) then - if (dNeq0(abs(q-exp(log(q))),1.0e-13_pReal)) error stop 'exp/log' - if (dNeq0(abs(q-log(exp(q))),1.0e-13_pReal)) error stop 'log/exp' - endif -#endif - -end subroutine selfTest - - -end module quaternions diff --git a/src/rotations.f90 b/src/rotations.f90 index 888e73762..73f8a16a1 100644 --- a/src/rotations.f90 +++ b/src/rotations.f90 @@ -47,16 +47,16 @@ !--------------------------------------------------------------------------------------------------- module rotations - use prec use IO use math - use quaternions implicit none private + real(pReal), parameter :: P = -1.0_pReal !< parameter for orientation conversion. + type, public :: rotation - type(quaternion) :: q + real(pReal), dimension(4) :: q contains procedure, public :: asQuaternion procedure, public :: asEulers @@ -103,7 +103,6 @@ contains !-------------------------------------------------------------------------------------------------- subroutine rotations_init - call quaternions_init print'(/,a)', ' <<<+- rotations init -+>>>'; flush(IO_STDOUT) print*, 'Rowenhorst et al., Modelling and Simulation in Materials Science and Engineering 23:083501, 2015' @@ -122,7 +121,7 @@ pure function asQuaternion(self) class(rotation), intent(in) :: self real(pReal), dimension(4) :: asQuaternion - asQuaternion = self%q%asArray() + asQuaternion = self%q end function asQuaternion !--------------------------------------------------------------------------------------------------- @@ -131,7 +130,7 @@ pure function asEulers(self) class(rotation), intent(in) :: self real(pReal), dimension(3) :: asEulers - asEulers = qu2eu(self%q%asArray()) + asEulers = qu2eu(self%q) end function asEulers !--------------------------------------------------------------------------------------------------- @@ -140,7 +139,7 @@ pure function asAxisAngle(self) class(rotation), intent(in) :: self real(pReal), dimension(4) :: asAxisAngle - asAxisAngle = qu2ax(self%q%asArray()) + asAxisAngle = qu2ax(self%q) end function asAxisAngle !--------------------------------------------------------------------------------------------------- @@ -149,7 +148,7 @@ pure function asMatrix(self) class(rotation), intent(in) :: self real(pReal), dimension(3,3) :: asMatrix - asMatrix = qu2om(self%q%asArray()) + asMatrix = qu2om(self%q) end function asMatrix !--------------------------------------------------------------------------------------------------- @@ -158,7 +157,7 @@ pure function asRodrigues(self) class(rotation), intent(in) :: self real(pReal), dimension(4) :: asRodrigues - asRodrigues = qu2ro(self%q%asArray()) + asRodrigues = qu2ro(self%q) end function asRodrigues !--------------------------------------------------------------------------------------------------- @@ -167,7 +166,7 @@ pure function asHomochoric(self) class(rotation), intent(in) :: self real(pReal), dimension(3) :: asHomochoric - asHomochoric = qu2ho(self%q%asArray()) + asHomochoric = qu2ho(self%q) end function asHomochoric @@ -259,7 +258,7 @@ pure elemental function rotRot__(self,R) result(rRot) type(rotation) :: rRot class(rotation), intent(in) :: self,R - rRot = rotation(self%q*R%q) + rRot = rotation(multiply_quaternion(self%q,R%q)) call rRot%standardize() end function rotRot__ @@ -272,14 +271,14 @@ pure elemental subroutine standardize(self) class(rotation), intent(inout) :: self - if (real(self%q) < 0.0_pReal) self%q = self%q%homomorphed() + if (self%q(1) < 0.0_pReal) self%q = - self%q end subroutine standardize !--------------------------------------------------------------------------------------------------- !> @author Marc De Graef, Carnegie Mellon University -!> @brief rotate a vector passively (default) or actively +!> @brief Rotate a vector passively (default) or actively. !--------------------------------------------------------------------------------------------------- pure function rotVector(self,v,active) result(vRot) @@ -288,9 +287,8 @@ pure function rotVector(self,v,active) result(vRot) real(pReal), intent(in), dimension(3) :: v logical, intent(in), optional :: active - real(pReal), dimension(3) :: v_normed - type(quaternion) :: q - logical :: passive + real(pReal), dimension(4) :: v_normed, q + logical :: passive if (present(active)) then passive = .not. active @@ -301,13 +299,13 @@ pure function rotVector(self,v,active) result(vRot) if (dEq0(norm2(v))) then vRot = v else - v_normed = v/norm2(v) + v_normed = [0.0_pReal,v]/norm2(v) if (passive) then - q = self%q * (quaternion([0.0_pReal, v_normed(1), v_normed(2), v_normed(3)]) * conjg(self%q) ) + q = multiply_quaternion(self%q, multiply_quaternion(v_normed, conjugate_quaternion(self%q))) else - q = conjg(self%q) * (quaternion([0.0_pReal, v_normed(1), v_normed(2), v_normed(3)]) * self%q ) + q = multiply_quaternion(conjugate_quaternion(self%q), multiply_quaternion(v_normed, self%q)) endif - vRot = q%aimag()*norm2(v) + vRot = q(2:4)*norm2(v) endif end function rotVector @@ -315,8 +313,8 @@ end function rotVector !--------------------------------------------------------------------------------------------------- !> @author Marc De Graef, Carnegie Mellon University -!> @brief rotate a rank-2 tensor passively (default) or actively -!> @details: rotation is based on rotation matrix +!> @brief Rotate a rank-2 tensor passively (default) or actively. +!> @details: Rotation is based on rotation matrix !--------------------------------------------------------------------------------------------------- pure function rotTensor2(self,T,active) result(tRot) @@ -403,7 +401,7 @@ pure elemental function misorientation(self,other) type(rotation) :: misorientation class(rotation), intent(in) :: self, other - misorientation%q = other%q * conjg(self%q) + misorientation%q = multiply_quaternion(other%q, [self%q(1),-self%q(2:4)]) end function misorientation @@ -1338,7 +1336,7 @@ end function cu2ho !-------------------------------------------------------------------------- !> @author Marc De Graef, Carnegie Mellon University !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH -!> @brief determine to which pyramid a point in a cubic grid belongs +!> @brief Determine to which pyramid a point in a cubic grid belongs. !-------------------------------------------------------------------------- pure function GetPyramidOrder(xyz) @@ -1362,7 +1360,39 @@ end function GetPyramidOrder !-------------------------------------------------------------------------------------------------- -!> @brief check correctness of some rotations functions +!> @brief Multiply two quaternions. +!-------------------------------------------------------------------------------------------------- +pure function multiply_quaternion(qu1,qu2) + + real(pReal), dimension(4), intent(in) :: qu1, qu2 + real(pReal), dimension(4) :: multiply_quaternion + + + multiply_quaternion(1) = qu1(1)*qu2(1) - qu1(2)*qu2(2) - qu1(3)*qu2(3) - qu1(4)*qu2(4) + multiply_quaternion(2) = qu1(1)*qu2(2) + qu1(2)*qu2(1) + P * (qu1(3)*qu2(4) - qu1(4)*qu2(3)) + multiply_quaternion(3) = qu1(1)*qu2(3) + qu1(3)*qu2(1) + P * (qu1(4)*qu2(2) - qu1(2)*qu2(4)) + multiply_quaternion(4) = qu1(1)*qu2(4) + qu1(4)*qu2(1) + P * (qu1(2)*qu2(3) - qu1(3)*qu2(2)) + +end function multiply_quaternion + + +!-------------------------------------------------------------------------------------------------- +!> @brief Calculate conjugate complex of a quaternion. +!-------------------------------------------------------------------------------------------------- +pure function conjugate_quaternion(qu) + + real(pReal), dimension(4), intent(in) :: qu + real(pReal), dimension(4) :: conjugate_quaternion + + + conjugate_quaternion = [qu(1), -qu(2), -qu(3), -qu(4)] + + +end function conjugate_quaternion + + +!-------------------------------------------------------------------------------------------------- +!> @brief Check correctness of some rotations functions. !-------------------------------------------------------------------------------------------------- subroutine selfTest @@ -1374,7 +1404,8 @@ subroutine selfTest real :: A,B integer :: i - do i=1,10 + + do i = 1, 10 #if defined(__GFORTRAN__) && __GNUC__<9 if(i<7) cycle From 69c11383cfd4e5f7d91a82b1001087bf67c21743 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 6 Jan 2021 13:37:45 +0100 Subject: [PATCH 16/19] better use function --- src/rotations.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/rotations.f90 b/src/rotations.f90 index 73f8a16a1..57dd16b53 100644 --- a/src/rotations.f90 +++ b/src/rotations.f90 @@ -401,7 +401,7 @@ pure elemental function misorientation(self,other) type(rotation) :: misorientation class(rotation), intent(in) :: self, other - misorientation%q = multiply_quaternion(other%q, [self%q(1),-self%q(2:4)]) + misorientation%q = multiply_quaternion(other%q, conjugate_quaternion(self%q)) end function misorientation From 437d91495bd54e29af1ba1b4e555a513844ea1f1 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 6 Jan 2021 14:24:02 +0100 Subject: [PATCH 17/19] No FEsolving --- src/CPFEM.f90 | 5 +- src/CPFEM2.f90 | 1 - src/DAMASK_marc.f90 | 1 - src/FEsolving.f90 | 15 ---- src/commercialFEM_fileList.f90 | 1 - src/constitutive.f90 | 55 ++++++-------- src/grid/discretization_grid.f90 | 4 -- src/grid/grid_mech_FEM.f90 | 1 - src/grid/grid_mech_spectral_basic.f90 | 1 - src/grid/grid_mech_spectral_polarisation.f90 | 1 - src/grid/spectral_utilities.f90 | 4 +- src/homogenization.f90 | 76 +++++++++----------- src/homogenization_mech.f90 | 24 +++---- src/marc/discretization_marc.f90 | 4 -- src/mesh/DAMASK_mesh.f90 | 1 - src/mesh/FEM_utilities.f90 | 2 +- src/mesh/discretization_mesh.f90 | 6 +- 17 files changed, 72 insertions(+), 130 deletions(-) delete mode 100644 src/FEsolving.f90 diff --git a/src/CPFEM.f90 b/src/CPFEM.f90 index abbcce04a..240688a8c 100644 --- a/src/CPFEM.f90 +++ b/src/CPFEM.f90 @@ -5,7 +5,6 @@ !-------------------------------------------------------------------------------------------------- module CPFEM use prec - use FEsolving use math use rotations use YAML_types @@ -197,11 +196,9 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS CPFEM_dcsde(1:6,1:6,ip,elCP) = ODD_JACOBIAN * math_eye(6) else validCalculation - FEsolving_execElem = elCP - FEsolving_execIP = ip if (debugCPFEM%extensive) & print'(a,i8,1x,i2)', '<< CPFEM >> calculation for elFE ip ',elFE,ip - call materialpoint_stressAndItsTangent(dt) + call materialpoint_stressAndItsTangent(dt,[ip,ip],[elCP,elCP]) terminalIllness: if (terminallyIll) then diff --git a/src/CPFEM2.f90 b/src/CPFEM2.f90 index 44b93d1cb..5a500875d 100644 --- a/src/CPFEM2.f90 +++ b/src/CPFEM2.f90 @@ -6,7 +6,6 @@ module CPFEM2 use prec use config - use FEsolving use math use rotations use YAML_types diff --git a/src/DAMASK_marc.f90 b/src/DAMASK_marc.f90 index ea7430c6b..0ad68445c 100644 --- a/src/DAMASK_marc.f90 +++ b/src/DAMASK_marc.f90 @@ -176,7 +176,6 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, & use DAMASK_interface use config use YAML_types - use FEsolving use discretization_marc use homogenization use CPFEM diff --git a/src/FEsolving.f90 b/src/FEsolving.f90 deleted file mode 100644 index 3fc1482d3..000000000 --- a/src/FEsolving.f90 +++ /dev/null @@ -1,15 +0,0 @@ -!-------------------------------------------------------------------------------------------------- -!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH -!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH -!> @brief global variables for flow control -!-------------------------------------------------------------------------------------------------- -module FEsolving - - implicit none - public - - integer, dimension(2) :: & - FEsolving_execElem, & !< for ping-pong scheme always whole range, otherwise one specific element - FEsolving_execIP !< for ping-pong scheme always range to max IP, otherwise one specific IP - -end module FEsolving diff --git a/src/commercialFEM_fileList.f90 b/src/commercialFEM_fileList.f90 index 08e7b9c1c..d8ab6390d 100644 --- a/src/commercialFEM_fileList.f90 +++ b/src/commercialFEM_fileList.f90 @@ -13,7 +13,6 @@ #include "math.f90" #include "quaternions.f90" #include "rotations.f90" -#include "FEsolving.f90" #include "element.f90" #include "HDF5_utilities.f90" #include "results.f90" diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 7e380f8cd..b3fb0b246 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -16,7 +16,6 @@ module constitutive use parallelization use HDF5_utilities use DAMASK_interface - use FEsolving use results implicit none @@ -940,8 +939,8 @@ subroutine crystallite_init flush(IO_STDOUT) !$OMP PARALLEL DO PRIVATE(ph,me) - do el = FEsolving_execElem(1),FEsolving_execElem(2) - do ip = FEsolving_execIP(1), FEsolving_execIP(2); do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) + do el = 1, size(material_phaseMemberAt,3) + do ip = 1, size(material_phaseMemberAt,2); do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) ph = material_phaseAt(co,el) me = material_phaseMemberAt(co,ip,el) @@ -967,14 +966,14 @@ subroutine crystallite_init crystallite_partitionedF0 = crystallite_F0 crystallite_partitionedF = crystallite_F0 - call crystallite_orientations() !$OMP PARALLEL DO PRIVATE(ph,me) - do el = FEsolving_execElem(1),FEsolving_execElem(2) - do ip = FEsolving_execIP(1),FEsolving_execIP(2) + do el = 1, size(material_phaseMemberAt,3) + do ip = 1, size(material_phaseMemberAt,2) do co = 1,homogenization_Nconstituents(material_homogenizationAt(el)) ph = material_phaseAt(co,el) me = material_phaseMemberAt(co,ip,el) + call crystallite_orientations(co,ip,el) call constitutive_plastic_dependentState(crystallite_partitionedF0(1:3,1:3,co,ip,el),co,ip,el) ! update dependent state variables to be consistent with basic states enddo enddo @@ -1089,7 +1088,7 @@ function crystallite_stressTangent(co,ip,el) result(dPdF) el !< counter in element loop integer :: & o, & - p, pp, m + p, ph, me real(pReal), dimension(3,3) :: devNull, & invSubFp0,invSubFi0,invFp,invFi, & @@ -1109,19 +1108,19 @@ function crystallite_stressTangent(co,ip,el) result(dPdF) real(pReal), dimension(9,9):: temp_99 logical :: error - pp = material_phaseAt(co,el) - m = material_phaseMemberAt(co,ip,el) + ph = material_phaseAt(co,el) + me = material_phaseMemberAt(co,ip,el) call constitutive_hooke_SandItsTangents(devNull,dSdFe,dSdFi, & crystallite_Fe(1:3,1:3,co,ip,el), & - constitutive_mech_Fi(pp)%data(1:3,1:3,m),co,ip,el) + constitutive_mech_Fi(ph)%data(1:3,1:3,me),co,ip,el) call constitutive_LiAndItsTangents(devNull,dLidS,dLidFi, & crystallite_S (1:3,1:3,co,ip,el), & - constitutive_mech_Fi(pp)%data(1:3,1:3,m), & + constitutive_mech_Fi(ph)%data(1:3,1:3,me), & co,ip,el) - invFp = math_inv33(constitutive_mech_Fp(pp)%data(1:3,1:3,m)) - invFi = math_inv33(constitutive_mech_Fi(pp)%data(1:3,1:3,m)) + invFp = math_inv33(constitutive_mech_Fp(ph)%data(1:3,1:3,me)) + invFi = math_inv33(constitutive_mech_Fi(ph)%data(1:3,1:3,me)) invSubFp0 = math_inv33(crystallite_subFp0(1:3,1:3,co,ip,el)) invSubFi0 = math_inv33(crystallite_subFi0(1:3,1:3,co,ip,el)) @@ -1150,7 +1149,7 @@ function crystallite_stressTangent(co,ip,el) result(dPdF) call constitutive_plastic_LpAndItsTangents(devNull,dLpdS,dLpdFi, & crystallite_S (1:3,1:3,co,ip,el), & - constitutive_mech_Fi(pp)%data(1:3,1:3,m),co,ip,el) + constitutive_mech_Fi(ph)%data(1:3,1:3,me),co,ip,el) dLpdS = math_mul3333xx3333(dLpdFi,dFidS) + dLpdS !-------------------------------------------------------------------------------------------------- @@ -1210,34 +1209,20 @@ end function crystallite_stressTangent !-------------------------------------------------------------------------------------------------- !> @brief calculates orientations !-------------------------------------------------------------------------------------------------- -subroutine crystallite_orientations +subroutine crystallite_orientations(co,ip,el) - integer & + integer, intent(in) :: & co, & !< counter in integration point component loop ip, & !< counter in integration point loop el !< counter in element loop - !$OMP PARALLEL DO - do el = FEsolving_execElem(1),FEsolving_execElem(2) - do ip = FEsolving_execIP(1),FEsolving_execIP(2) - do co = 1,homogenization_Nconstituents(material_homogenizationAt(el)) - call crystallite_orientation(co,ip,el)%fromMatrix(transpose(math_rotationalPart(crystallite_Fe(1:3,1:3,co,ip,el)))) - enddo; enddo; enddo - !$OMP END PARALLEL DO + call crystallite_orientation(co,ip,el)%fromMatrix(transpose(math_rotationalPart(crystallite_Fe(1:3,1:3,co,ip,el)))) + + if (plasticState(material_phaseAt(1,el))%nonlocal) & + call plastic_nonlocal_updateCompatibility(crystallite_orientation, & + phase_plasticityInstance(material_phaseAt(1,el)),ip,el) - nonlocalPresent: if (any(plasticState%nonlocal)) then - !$OMP PARALLEL DO - do el = FEsolving_execElem(1),FEsolving_execElem(2) - if (plasticState(material_phaseAt(1,el))%nonlocal) then - do ip = FEsolving_execIP(1),FEsolving_execIP(2) - call plastic_nonlocal_updateCompatibility(crystallite_orientation, & - phase_plasticityInstance(material_phaseAt(1,el)),ip,el) - enddo - endif - enddo - !$OMP END PARALLEL DO - endif nonlocalPresent end subroutine crystallite_orientations diff --git a/src/grid/discretization_grid.f90 b/src/grid/discretization_grid.f90 index 1b3700c14..48ad5b7e1 100644 --- a/src/grid/discretization_grid.f90 +++ b/src/grid/discretization_grid.f90 @@ -19,7 +19,6 @@ module discretization_grid use results use discretization use geometry_plastic_nonlocal - use FEsolving implicit none private @@ -117,9 +116,6 @@ subroutine discretization_grid_init(restart) (grid(1)+1) * (grid(2)+1) * grid3,& ! ...unless not last process worldrank+1==worldsize)) - FEsolving_execElem = [1,product(myGrid)] ! parallel loop bounds set to comprise all elements - FEsolving_execIP = [1,1] ! parallel loop bounds set to comprise the only IP - !-------------------------------------------------------------------------------------------------- ! store geometry information for post processing if(.not. restart) then diff --git a/src/grid/grid_mech_FEM.f90 b/src/grid/grid_mech_FEM.f90 index cdf806b35..003f568c6 100644 --- a/src/grid/grid_mech_FEM.f90 +++ b/src/grid/grid_mech_FEM.f90 @@ -18,7 +18,6 @@ module grid_mech_FEM use math use rotations use spectral_utilities - use FEsolving use config use homogenization use discretization diff --git a/src/grid/grid_mech_spectral_basic.f90 b/src/grid/grid_mech_spectral_basic.f90 index ebaaf3b55..9bc36165f 100644 --- a/src/grid/grid_mech_spectral_basic.f90 +++ b/src/grid/grid_mech_spectral_basic.f90 @@ -18,7 +18,6 @@ module grid_mech_spectral_basic use math use rotations use spectral_utilities - use FEsolving use config use homogenization use discretization_grid diff --git a/src/grid/grid_mech_spectral_polarisation.f90 b/src/grid/grid_mech_spectral_polarisation.f90 index 9f2a17c97..7160c1adc 100644 --- a/src/grid/grid_mech_spectral_polarisation.f90 +++ b/src/grid/grid_mech_spectral_polarisation.f90 @@ -18,7 +18,6 @@ module grid_mech_spectral_polarisation use math use rotations use spectral_utilities - use FEsolving use config use homogenization use discretization_grid diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index c0c84233d..e8bae223a 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -810,9 +810,9 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,& print'(/,a)', ' ... evaluating constitutive response ......................................' flush(IO_STDOUT) - homogenization_F = reshape(F,[3,3,product(grid(1:2))*grid3]) ! set materialpoint target F to estimated field + homogenization_F = reshape(F,[3,3,product(grid(1:2))*grid3]) ! set materialpoint target F to estimated field - call materialpoint_stressAndItsTangent(timeinc) ! calculate P field + call materialpoint_stressAndItsTangent(timeinc,[1,1],[1,product(grid(1:2))*grid3]) ! calculate P field P = reshape(homogenization_P, [3,3,grid(1),grid(2),grid3]) P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt ! average of P diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 27fdb6064..ebf5fd50d 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -11,7 +11,6 @@ module homogenization use math use material use constitutive - use FEsolving use discretization use thermal_isothermal use thermal_conduction @@ -144,27 +143,29 @@ end subroutine homogenization_init !-------------------------------------------------------------------------------------------------- !> @brief parallelized calculation of stress and corresponding tangent at material points !-------------------------------------------------------------------------------------------------- -subroutine materialpoint_stressAndItsTangent(dt) +subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execElem) real(pReal), intent(in) :: dt !< time increment + integer, dimension(2), intent(in) :: FEsolving_execElem, FEsolving_execIP integer :: & NiterationHomog, & NiterationMPstate, & ip, & !< integration point number el, & !< element number - myNgrains, co, ce + myNgrains, co, ce, ho real(pReal) :: & subFrac, & subStep logical :: & - requested, & converged logical, dimension(2) :: & doneAndHappy -!$OMP PARALLEL DO PRIVATE(ce,myNgrains,NiterationMPstate,NiterationHomog,subFrac,converged,subStep,requested,doneAndHappy) +!$OMP PARALLEL DO PRIVATE(ce,ho,myNgrains,NiterationMPstate,NiterationHomog,subFrac,converged,subStep,doneAndHappy) do el = FEsolving_execElem(1),FEsolving_execElem(2) + ho = material_homogenizationAt(el) + myNgrains = homogenization_Nconstituents(ho) do ip = FEsolving_execIP(1),FEsolving_execIP(2) !-------------------------------------------------------------------------------------------------- @@ -174,21 +175,19 @@ subroutine materialpoint_stressAndItsTangent(dt) subFrac = 0.0_pReal converged = .false. ! pretend failed step ... subStep = 1.0_pReal/num%subStepSizeHomog ! ... larger then the requested calculation - requested = .true. ! everybody requires calculation - if (homogState(material_homogenizationAt(el))%sizeState > 0) & - homogState(material_homogenizationAt(el))%subState0(:,material_homogenizationMemberAt(ip,el)) = & - homogState(material_homogenizationAt(el))%State0( :,material_homogenizationMemberAt(ip,el)) + if (homogState(ho)%sizeState > 0) & + homogState(ho)%subState0(:,material_homogenizationMemberAt(ip,el)) = & + homogState(ho)%State0( :,material_homogenizationMemberAt(ip,el)) + + if (damageState(ho)%sizeState > 0) & + damageState(ho)%subState0(:,material_homogenizationMemberAt(ip,el)) = & + damageState(ho)%State0( :,material_homogenizationMemberAt(ip,el)) - if (damageState(material_homogenizationAt(el))%sizeState > 0) & - damageState(material_homogenizationAt(el))%subState0(:,material_homogenizationMemberAt(ip,el)) = & - damageState(material_homogenizationAt(el))%State0( :,material_homogenizationMemberAt(ip,el)) NiterationHomog = 0 cutBackLooping: do while (.not. terminallyIll .and. subStep > num%subStepMinHomog) - myNgrains = homogenization_Nconstituents(material_homogenizationAt(el)) - if (converged) then subFrac = subFrac + subStep subStep = min(1.0_pReal-subFrac,num%stepIncreaseHomog*subStep) ! introduce flexibility for step increase/acceleration @@ -198,22 +197,20 @@ subroutine materialpoint_stressAndItsTangent(dt) ! wind forward grain starting point call constitutive_windForward(ip,el) - if(homogState(material_homogenizationAt(el))%sizeState > 0) & - homogState(material_homogenizationAt(el))%subState0(:,material_homogenizationMemberAt(ip,el)) = & - homogState(material_homogenizationAt(el))%State (:,material_homogenizationMemberAt(ip,el)) - if(damageState(material_homogenizationAt(el))%sizeState > 0) & - damageState(material_homogenizationAt(el))%subState0(:,material_homogenizationMemberAt(ip,el)) = & - damageState(material_homogenizationAt(el))%State (:,material_homogenizationMemberAt(ip,el)) + if(homogState(ho)%sizeState > 0) & + homogState(ho)%subState0(:,material_homogenizationMemberAt(ip,el)) = & + homogState(ho)%State (:,material_homogenizationMemberAt(ip,el)) + if(damageState(ho)%sizeState > 0) & + damageState(ho)%subState0(:,material_homogenizationMemberAt(ip,el)) = & + damageState(ho)%State (:,material_homogenizationMemberAt(ip,el)) endif steppingNeeded - else if ( (myNgrains == 1 .and. subStep <= 1.0 ) .or. & ! single grain already tried internal subStepping in crystallite num%subStepSizeHomog * subStep <= num%subStepMinHomog ) then ! would require too small subStep ! cutback makes no sense - if (.not. terminallyIll) then ! so first signals terminally ill... + if (.not. terminallyIll) & ! so first signals terminally ill... print*, ' Integration point ', ip,' at element ', el, ' terminally ill' - endif terminallyIll = .true. ! ...and kills all others else ! cutback makes sense subStep = num%subStepSizeHomog * subStep ! crystallite had severe trouble, so do a significant cutback @@ -221,23 +218,19 @@ subroutine materialpoint_stressAndItsTangent(dt) call crystallite_restore(ip,el,subStep < 1.0_pReal) call constitutive_restore(ip,el) - if(homogState(material_homogenizationAt(el))%sizeState > 0) & - homogState(material_homogenizationAt(el))%State( :,material_homogenizationMemberAt(ip,el)) = & - homogState(material_homogenizationAt(el))%subState0(:,material_homogenizationMemberAt(ip,el)) - if(damageState(material_homogenizationAt(el))%sizeState > 0) & - damageState(material_homogenizationAt(el))%State( :,material_homogenizationMemberAt(ip,el)) = & - damageState(material_homogenizationAt(el))%subState0(:,material_homogenizationMemberAt(ip,el)) + if(homogState(ho)%sizeState > 0) & + homogState(ho)%State( :,material_homogenizationMemberAt(ip,el)) = & + homogState(ho)%subState0(:,material_homogenizationMemberAt(ip,el)) + if(damageState(ho)%sizeState > 0) & + damageState(ho)%State( :,material_homogenizationMemberAt(ip,el)) = & + damageState(ho)%subState0(:,material_homogenizationMemberAt(ip,el)) endif endif - if (subStep > num%subStepMinHomog) then - requested = .true. - doneAndHappy = [.false.,.true.] - endif - + if (subStep > num%subStepMinHomog) doneAndHappy = [.false.,.true.] NiterationMPstate = 0 - convergenceLooping: do while (.not. terminallyIll .and. requested & + convergenceLooping: do while (.not. terminallyIll & .and. .not. doneAndHappy(1) & .and. NiterationMPstate < num%nMPstate) NiterationMPstate = NiterationMPstate + 1 @@ -245,7 +238,7 @@ subroutine materialpoint_stressAndItsTangent(dt) !-------------------------------------------------------------------------------------------------- ! deformation partitioning - if(requested .and. .not. doneAndHappy(1)) then ! requested but not yet done + if (.not. doneAndHappy(1)) then ce = (el-1)*discretization_nIPs + ip call mech_partition(homogenization_F0(1:3,1:3,ce) & + (homogenization_F(1:3,1:3,ce)-homogenization_F0(1:3,1:3,ce))& @@ -255,10 +248,7 @@ subroutine materialpoint_stressAndItsTangent(dt) do co = 1, myNgrains converged = converged .and. crystallite_stress(dt*subStep,co,ip,el) enddo - endif - - if (requested .and. .not. doneAndHappy(1)) then if (.not. converged) then doneAndHappy = [.true.,.false.] else @@ -281,10 +271,14 @@ subroutine materialpoint_stressAndItsTangent(dt) !$OMP END PARALLEL DO if (.not. terminallyIll ) then - call crystallite_orientations() ! calculate crystal orientations - !$OMP PARALLEL DO + !$OMP PARALLEL DO PRIVATE(ho,myNgrains) elementLooping3: do el = FEsolving_execElem(1),FEsolving_execElem(2) + ho = material_homogenizationAt(el) + myNgrains = homogenization_Nconstituents(ho) IpLooping3: do ip = FEsolving_execIP(1),FEsolving_execIP(2) + do co = 1, myNgrains + call crystallite_orientations(co,ip,el) + enddo call mech_homogenize(ip,el) enddo IpLooping3 enddo elementLooping3 diff --git a/src/homogenization_mech.f90 b/src/homogenization_mech.f90 index 56f1e554f..e4499e9b7 100644 --- a/src/homogenization_mech.f90 +++ b/src/homogenization_mech.f90 @@ -128,35 +128,35 @@ module subroutine mech_homogenize(ip,el) integer, intent(in) :: & ip, & !< integration point el !< element number - integer :: c,m + integer :: co,ce real(pReal) :: dPdFs(3,3,3,3,homogenization_Nconstituents(material_homogenizationAt(el))) - m = (el-1)* discretization_nIPs + ip + ce = (el-1)* discretization_nIPs + ip chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el))) case (HOMOGENIZATION_NONE_ID) chosenHomogenization - homogenization_P(1:3,1:3,m) = crystallite_P(1:3,1:3,1,ip,el) - homogenization_dPdF(1:3,1:3,1:3,1:3,m) = crystallite_stressTangent(1,ip,el) + homogenization_P(1:3,1:3,ce) = crystallite_P(1:3,1:3,1,ip,el) + homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = crystallite_stressTangent(1,ip,el) case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization - do c = 1, homogenization_Nconstituents(material_homogenizationAt(el)) - dPdFs(:,:,:,:,c) = crystallite_stressTangent(c,ip,el) + do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) + dPdFs(:,:,:,:,co) = crystallite_stressTangent(co,ip,el) enddo call mech_isostrain_averageStressAndItsTangent(& - homogenization_P(1:3,1:3,m), & - homogenization_dPdF(1:3,1:3,1:3,1:3,m),& + homogenization_P(1:3,1:3,ce), & + homogenization_dPdF(1:3,1:3,1:3,1:3,ce),& crystallite_P(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), & dPdFs, & homogenization_typeInstance(material_homogenizationAt(el))) case (HOMOGENIZATION_RGC_ID) chosenHomogenization - do c = 1, homogenization_Nconstituents(material_homogenizationAt(el)) - dPdFs(:,:,:,:,c) = crystallite_stressTangent(c,ip,el) + do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) + dPdFs(:,:,:,:,co) = crystallite_stressTangent(co,ip,el) enddo call mech_RGC_averageStressAndItsTangent(& - homogenization_P(1:3,1:3,m), & - homogenization_dPdF(1:3,1:3,1:3,1:3,m),& + homogenization_P(1:3,1:3,ce), & + homogenization_dPdF(1:3,1:3,1:3,1:3,ce),& crystallite_P(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), & dPdFs, & homogenization_typeInstance(material_homogenizationAt(el))) diff --git a/src/marc/discretization_marc.f90 b/src/marc/discretization_marc.f90 index ca0b54b73..675e57bd3 100644 --- a/src/marc/discretization_marc.f90 +++ b/src/marc/discretization_marc.f90 @@ -12,7 +12,6 @@ module discretization_marc use DAMASK_interface use IO use config - use FEsolving use element use discretization use geometry_plastic_nonlocal @@ -89,9 +88,6 @@ subroutine discretization_marc_init if (debug_e < 1 .or. debug_e > nElems) call IO_error(602,ext_msg='element') if (debug_i < 1 .or. debug_i > elem%nIPs) call IO_error(602,ext_msg='IP') - FEsolving_execElem = [1,nElems] - FEsolving_execIP = [1,elem%nIPs] - allocate(cellNodeDefinition(elem%nNodes-1)) allocate(connectivity_cell(elem%NcellNodesPerCell,elem%nIPs,nElems)) call buildCells(connectivity_cell,cellNodeDefinition,& diff --git a/src/mesh/DAMASK_mesh.f90 b/src/mesh/DAMASK_mesh.f90 index 1e353892c..7369520c1 100644 --- a/src/mesh/DAMASK_mesh.f90 +++ b/src/mesh/DAMASK_mesh.f90 @@ -15,7 +15,6 @@ program DAMASK_mesh use IO use math use CPFEM2 - use FEsolving use config use discretization_mesh use FEM_Utilities diff --git a/src/mesh/FEM_utilities.f90 b/src/mesh/FEM_utilities.f90 index cb81f1f0c..2f3633e11 100644 --- a/src/mesh/FEM_utilities.f90 +++ b/src/mesh/FEM_utilities.f90 @@ -160,7 +160,7 @@ subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData) print'(/,a)', ' ... evaluating constitutive response ......................................' - call materialpoint_stressAndItsTangent(timeinc) ! calculate P field + call materialpoint_stressAndItsTangent(timeinc,[1,mesh_maxNips],[1,mesh_NcpElems]) ! calculate P field cutBack = .false. ! reset cutBack status diff --git a/src/mesh/discretization_mesh.f90 b/src/mesh/discretization_mesh.f90 index 7dbd05e46..21c5feace 100644 --- a/src/mesh/discretization_mesh.f90 +++ b/src/mesh/discretization_mesh.f90 @@ -18,7 +18,6 @@ module discretization_mesh use config use discretization use results - use FEsolving use FEM_quadrature use YAML_types use prec @@ -30,7 +29,7 @@ module discretization_mesh mesh_Nboundaries, & mesh_NcpElemsGlobal - integer :: & + integer, public, protected :: & mesh_NcpElems !< total number of CP elements in mesh !!!! BEGIN DEPRECATED !!!!! @@ -174,9 +173,6 @@ subroutine discretization_mesh_init(restart) if (debug_element < 1 .or. debug_element > mesh_NcpElems) call IO_error(602,ext_msg='element') if (debug_ip < 1 .or. debug_ip > mesh_maxNips) call IO_error(602,ext_msg='IP') - FEsolving_execElem = [1,mesh_NcpElems] ! parallel loop bounds set to comprise all DAMASK elements - FEsolving_execIP = [1,mesh_maxNips] - allocate(mesh_node0(3,mesh_Nnodes),source=0.0_pReal) call discretization_init(materialAt,& From 959c18c85e4ca397453dc3db15096ef7eb3e9108 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 6 Jan 2021 14:24:52 +0100 Subject: [PATCH 18/19] No crystallite _converged --- src/CPFEM.f90 | 5 +- src/CPFEM2.f90 | 1 - src/DAMASK_marc.f90 | 1 - src/FEsolving.f90 | 15 ---- src/commercialFEM_fileList.f90 | 1 - src/constitutive.f90 | 83 ++++++++------------ src/constitutive_mech.f90 | 73 +++++++---------- src/grid/discretization_grid.f90 | 4 - src/grid/grid_mech_FEM.f90 | 1 - src/grid/grid_mech_spectral_basic.f90 | 1 - src/grid/grid_mech_spectral_polarisation.f90 | 1 - src/grid/spectral_utilities.f90 | 4 +- src/homogenization.f90 | 76 +++++++++--------- src/homogenization_mech.f90 | 24 +++--- src/marc/discretization_marc.f90 | 4 - src/mesh/DAMASK_mesh.f90 | 1 - src/mesh/FEM_utilities.f90 | 2 +- src/mesh/discretization_mesh.f90 | 6 +- 18 files changed, 113 insertions(+), 190 deletions(-) delete mode 100644 src/FEsolving.f90 diff --git a/src/CPFEM.f90 b/src/CPFEM.f90 index abbcce04a..240688a8c 100644 --- a/src/CPFEM.f90 +++ b/src/CPFEM.f90 @@ -5,7 +5,6 @@ !-------------------------------------------------------------------------------------------------- module CPFEM use prec - use FEsolving use math use rotations use YAML_types @@ -197,11 +196,9 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS CPFEM_dcsde(1:6,1:6,ip,elCP) = ODD_JACOBIAN * math_eye(6) else validCalculation - FEsolving_execElem = elCP - FEsolving_execIP = ip if (debugCPFEM%extensive) & print'(a,i8,1x,i2)', '<< CPFEM >> calculation for elFE ip ',elFE,ip - call materialpoint_stressAndItsTangent(dt) + call materialpoint_stressAndItsTangent(dt,[ip,ip],[elCP,elCP]) terminalIllness: if (terminallyIll) then diff --git a/src/CPFEM2.f90 b/src/CPFEM2.f90 index 44b93d1cb..5a500875d 100644 --- a/src/CPFEM2.f90 +++ b/src/CPFEM2.f90 @@ -6,7 +6,6 @@ module CPFEM2 use prec use config - use FEsolving use math use rotations use YAML_types diff --git a/src/DAMASK_marc.f90 b/src/DAMASK_marc.f90 index ea7430c6b..0ad68445c 100644 --- a/src/DAMASK_marc.f90 +++ b/src/DAMASK_marc.f90 @@ -176,7 +176,6 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, & use DAMASK_interface use config use YAML_types - use FEsolving use discretization_marc use homogenization use CPFEM diff --git a/src/FEsolving.f90 b/src/FEsolving.f90 deleted file mode 100644 index 3fc1482d3..000000000 --- a/src/FEsolving.f90 +++ /dev/null @@ -1,15 +0,0 @@ -!-------------------------------------------------------------------------------------------------- -!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH -!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH -!> @brief global variables for flow control -!-------------------------------------------------------------------------------------------------- -module FEsolving - - implicit none - public - - integer, dimension(2) :: & - FEsolving_execElem, & !< for ping-pong scheme always whole range, otherwise one specific element - FEsolving_execIP !< for ping-pong scheme always range to max IP, otherwise one specific IP - -end module FEsolving diff --git a/src/commercialFEM_fileList.f90 b/src/commercialFEM_fileList.f90 index 08e7b9c1c..d8ab6390d 100644 --- a/src/commercialFEM_fileList.f90 +++ b/src/commercialFEM_fileList.f90 @@ -13,7 +13,6 @@ #include "math.f90" #include "quaternions.f90" #include "rotations.f90" -#include "FEsolving.f90" #include "element.f90" #include "HDF5_utilities.f90" #include "results.f90" diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 7e380f8cd..e65ce864d 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -16,7 +16,6 @@ module constitutive use parallelization use HDF5_utilities use DAMASK_interface - use FEsolving use results implicit none @@ -65,10 +64,6 @@ module constitutive real(pReal), dimension(:,:,:,:,:), allocatable, public :: & crystallite_partitionedF !< def grad to be reached at end of homog inc - logical, dimension(:,:,:), allocatable :: & - crystallite_converged !< convergence flag - - type :: tTensorContainer real(pReal), dimension(:,:,:), allocatable :: data end type @@ -186,10 +181,10 @@ module constitutive ! == cleaned:end =================================================================================== - module function crystallite_stress(dt,co,ip,el) + module function crystallite_stress(dt,co,ip,el) result(converged_) real(pReal), intent(in) :: dt integer, intent(in) :: co, ip, el - logical :: crystallite_stress + logical :: converged_ end function crystallite_stress module function constitutive_homogenizedC(co,ip,el) result(C) @@ -873,10 +868,8 @@ subroutine crystallite_init source = crystallite_partitionedF) allocate(crystallite_subdt(cMax,iMax,eMax),source=0.0_pReal) - allocate(crystallite_orientation(cMax,iMax,eMax)) - allocate(crystallite_converged(cMax,iMax,eMax), source=.true.) num_crystallite => config_numerics%get('crystallite',defaultVal=emptyDict) @@ -940,8 +933,8 @@ subroutine crystallite_init flush(IO_STDOUT) !$OMP PARALLEL DO PRIVATE(ph,me) - do el = FEsolving_execElem(1),FEsolving_execElem(2) - do ip = FEsolving_execIP(1), FEsolving_execIP(2); do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) + do el = 1, size(material_phaseMemberAt,3) + do ip = 1, size(material_phaseMemberAt,2); do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) ph = material_phaseAt(co,el) me = material_phaseMemberAt(co,ip,el) @@ -967,14 +960,14 @@ subroutine crystallite_init crystallite_partitionedF0 = crystallite_F0 crystallite_partitionedF = crystallite_F0 - call crystallite_orientations() !$OMP PARALLEL DO PRIVATE(ph,me) - do el = FEsolving_execElem(1),FEsolving_execElem(2) - do ip = FEsolving_execIP(1),FEsolving_execIP(2) + do el = 1, size(material_phaseMemberAt,3) + do ip = 1, size(material_phaseMemberAt,2) do co = 1,homogenization_Nconstituents(material_homogenizationAt(el)) ph = material_phaseAt(co,el) me = material_phaseMemberAt(co,ip,el) + call crystallite_orientations(co,ip,el) call constitutive_plastic_dependentState(crystallite_partitionedF0(1:3,1:3,co,ip,el),co,ip,el) ! update dependent state variables to be consistent with basic states enddo enddo @@ -1089,7 +1082,7 @@ function crystallite_stressTangent(co,ip,el) result(dPdF) el !< counter in element loop integer :: & o, & - p, pp, m + p, ph, me real(pReal), dimension(3,3) :: devNull, & invSubFp0,invSubFi0,invFp,invFi, & @@ -1109,19 +1102,19 @@ function crystallite_stressTangent(co,ip,el) result(dPdF) real(pReal), dimension(9,9):: temp_99 logical :: error - pp = material_phaseAt(co,el) - m = material_phaseMemberAt(co,ip,el) + ph = material_phaseAt(co,el) + me = material_phaseMemberAt(co,ip,el) call constitutive_hooke_SandItsTangents(devNull,dSdFe,dSdFi, & crystallite_Fe(1:3,1:3,co,ip,el), & - constitutive_mech_Fi(pp)%data(1:3,1:3,m),co,ip,el) + constitutive_mech_Fi(ph)%data(1:3,1:3,me),co,ip,el) call constitutive_LiAndItsTangents(devNull,dLidS,dLidFi, & crystallite_S (1:3,1:3,co,ip,el), & - constitutive_mech_Fi(pp)%data(1:3,1:3,m), & + constitutive_mech_Fi(ph)%data(1:3,1:3,me), & co,ip,el) - invFp = math_inv33(constitutive_mech_Fp(pp)%data(1:3,1:3,m)) - invFi = math_inv33(constitutive_mech_Fi(pp)%data(1:3,1:3,m)) + invFp = math_inv33(constitutive_mech_Fp(ph)%data(1:3,1:3,me)) + invFi = math_inv33(constitutive_mech_Fi(ph)%data(1:3,1:3,me)) invSubFp0 = math_inv33(crystallite_subFp0(1:3,1:3,co,ip,el)) invSubFi0 = math_inv33(crystallite_subFi0(1:3,1:3,co,ip,el)) @@ -1150,7 +1143,7 @@ function crystallite_stressTangent(co,ip,el) result(dPdF) call constitutive_plastic_LpAndItsTangents(devNull,dLpdS,dLpdFi, & crystallite_S (1:3,1:3,co,ip,el), & - constitutive_mech_Fi(pp)%data(1:3,1:3,m),co,ip,el) + constitutive_mech_Fi(ph)%data(1:3,1:3,me),co,ip,el) dLpdS = math_mul3333xx3333(dLpdFi,dFidS) + dLpdS !-------------------------------------------------------------------------------------------------- @@ -1210,34 +1203,20 @@ end function crystallite_stressTangent !-------------------------------------------------------------------------------------------------- !> @brief calculates orientations !-------------------------------------------------------------------------------------------------- -subroutine crystallite_orientations +subroutine crystallite_orientations(co,ip,el) - integer & + integer, intent(in) :: & co, & !< counter in integration point component loop ip, & !< counter in integration point loop el !< counter in element loop - !$OMP PARALLEL DO - do el = FEsolving_execElem(1),FEsolving_execElem(2) - do ip = FEsolving_execIP(1),FEsolving_execIP(2) - do co = 1,homogenization_Nconstituents(material_homogenizationAt(el)) - call crystallite_orientation(co,ip,el)%fromMatrix(transpose(math_rotationalPart(crystallite_Fe(1:3,1:3,co,ip,el)))) - enddo; enddo; enddo - !$OMP END PARALLEL DO + call crystallite_orientation(co,ip,el)%fromMatrix(transpose(math_rotationalPart(crystallite_Fe(1:3,1:3,co,ip,el)))) + + if (plasticState(material_phaseAt(1,el))%nonlocal) & + call plastic_nonlocal_updateCompatibility(crystallite_orientation, & + phase_plasticityInstance(material_phaseAt(1,el)),ip,el) - nonlocalPresent: if (any(plasticState%nonlocal)) then - !$OMP PARALLEL DO - do el = FEsolving_execElem(1),FEsolving_execElem(2) - if (plasticState(material_phaseAt(1,el))%nonlocal) then - do ip = FEsolving_execIP(1),FEsolving_execIP(2) - call plastic_nonlocal_updateCompatibility(crystallite_orientation, & - phase_plasticityInstance(material_phaseAt(1,el)),ip,el) - enddo - endif - enddo - !$OMP END PARALLEL DO - endif nonlocalPresent end subroutine crystallite_orientations @@ -1268,7 +1247,7 @@ end function crystallite_push33ToRef !> @brief integrate stress, state with adaptive 1st order explicit Euler method !> using Fixed Point Iteration to adapt the stepsize !-------------------------------------------------------------------------------------------------- -subroutine integrateSourceState(co,ip,el) +function integrateSourceState(co,ip,el) result(broken) integer, intent(in) :: & el, & !< element index in element loop @@ -1288,12 +1267,13 @@ subroutine integrateSourceState(co,ip,el) r ! state residuum real(pReal), dimension(constitutive_source_maxSizeDotState,2,maxval(phase_Nsources)) :: source_dotState logical :: & - broken + broken, converged_ ph = material_phaseAt(co,el) me = material_phaseMemberAt(co,ip,el) + converged_ = .true. broken = constitutive_thermal_collectDotState(ph,me) broken = broken .or. constitutive_damage_collectDotState(crystallite_S(1:3,1:3,co,ip,el), co,ip,el,ph,me) if(broken) return @@ -1328,19 +1308,20 @@ subroutine integrateSourceState(co,ip,el) - sourceState(ph)%p(so)%dotState (1:size_so(so),me) * crystallite_subdt(co,ip,el) sourceState(ph)%p(so)%state(1:size_so(so),me) = sourceState(ph)%p(so)%state(1:size_so(so),me) & - r(1:size_so(so)) - crystallite_converged(co,ip,el) = & - crystallite_converged(co,ip,el) .and. converged(r(1:size_so(so)), & - sourceState(ph)%p(so)%state(1:size_so(so),me), & - sourceState(ph)%p(so)%atol(1:size_so(so))) + converged_ = converged_ .and. converged(r(1:size_so(so)), & + sourceState(ph)%p(so)%state(1:size_so(so),me), & + sourceState(ph)%p(so)%atol(1:size_so(so))) enddo - if(crystallite_converged(co,ip,el)) then + if(converged_) then broken = constitutive_damage_deltaState(crystallite_Fe(1:3,1:3,co,ip,el),co,ip,el,ph,me) exit iteration endif enddo iteration + broken = broken .or. .not. converged_ + contains @@ -1364,7 +1345,7 @@ subroutine integrateSourceState(co,ip,el) end function damper -end subroutine integrateSourceState +end function integrateSourceState !-------------------------------------------------------------------------------------------------- diff --git a/src/constitutive_mech.f90 b/src/constitutive_mech.f90 index 7a2224ede..de6f2ae9f 100644 --- a/src/constitutive_mech.f90 +++ b/src/constitutive_mech.f90 @@ -951,7 +951,7 @@ end function integrateStress !> @brief integrate stress, state with adaptive 1st order explicit Euler method !> using Fixed Point Iteration to adapt the stepsize !-------------------------------------------------------------------------------------------------- -subroutine integrateStateFPI(F_0,F,Delta_t,co,ip,el) +function integrateStateFPI(F_0,F,Delta_t,co,ip,el) result(broken) real(pReal), intent(in),dimension(3,3) :: F_0,F real(pReal), intent(in) :: Delta_t @@ -1004,11 +1004,7 @@ subroutine integrateStateFPI(F_0,F,Delta_t,co,ip,el) - plasticState(ph)%dotState (1:size_pl,me) * Delta_t plasticState(ph)%state(1:size_pl,me) = plasticState(ph)%state(1:size_pl,me) & - r(1:size_pl) - crystallite_converged(co,ip,el) = converged(r(1:size_pl), & - plasticState(ph)%state(1:size_pl,me), & - plasticState(ph)%atol(1:size_pl)) - - if(crystallite_converged(co,ip,el)) then + if (converged(r(1:size_pl),plasticState(ph)%state(1:size_pl,me),plasticState(ph)%atol(1:size_pl))) then broken = constitutive_deltaState(crystallite_S(1:3,1:3,co,ip,el), & constitutive_mech_Fi(ph)%data(1:3,1:3,me),co,ip,el,ph,me) exit iteration @@ -1016,7 +1012,6 @@ subroutine integrateStateFPI(F_0,F,Delta_t,co,ip,el) enddo iteration - contains !-------------------------------------------------------------------------------------------------- @@ -1039,13 +1034,13 @@ subroutine integrateStateFPI(F_0,F,Delta_t,co,ip,el) end function damper -end subroutine integrateStateFPI +end function integrateStateFPI !-------------------------------------------------------------------------------------------------- !> @brief integrate state with 1st order explicit Euler method !-------------------------------------------------------------------------------------------------- -subroutine integrateStateEuler(F_0,F,Delta_t,co,ip,el) +function integrateStateEuler(F_0,F,Delta_t,co,ip,el) result(broken) real(pReal), intent(in),dimension(3,3) :: F_0,F real(pReal), intent(in) :: Delta_t @@ -1075,15 +1070,14 @@ subroutine integrateStateEuler(F_0,F,Delta_t,co,ip,el) if(broken) return broken = integrateStress(F,Delta_t,co,ip,el) - crystallite_converged(co,ip,el) = .not. broken -end subroutine integrateStateEuler +end function integrateStateEuler !-------------------------------------------------------------------------------------------------- !> @brief integrate stress, state with 1st order Euler method with adaptive step size !-------------------------------------------------------------------------------------------------- -subroutine integrateStateAdaptiveEuler(F_0,F,Delta_t,co,ip,el) +function integrateStateAdaptiveEuler(F_0,F,Delta_t,co,ip,el) result(broken) real(pReal), intent(in),dimension(3,3) :: F_0,F real(pReal), intent(in) :: Delta_t @@ -1123,24 +1117,22 @@ subroutine integrateStateAdaptiveEuler(F_0,F,Delta_t,co,ip,el) broken = mech_collectDotState(Delta_t, co,ip,el,ph,me) if(broken) return + broken = .not. converged(residuum_plastic(1:sizeDotState) + 0.5_pReal * plasticState(ph)%dotState(:,me) * Delta_t, & + plasticState(ph)%state(1:sizeDotState,me), & + plasticState(ph)%atol(1:sizeDotState)) - sizeDotState = plasticState(ph)%sizeDotState - crystallite_converged(co,ip,el) = converged(residuum_plastic(1:sizeDotState) & - + 0.5_pReal * plasticState(ph)%dotState(:,me) * Delta_t, & - plasticState(ph)%state(1:sizeDotState,me), & - plasticState(ph)%atol(1:sizeDotState)) - -end subroutine integrateStateAdaptiveEuler +end function integrateStateAdaptiveEuler !--------------------------------------------------------------------------------------------------- !> @brief Integrate state (including stress integration) with the classic Runge Kutta method !--------------------------------------------------------------------------------------------------- -subroutine integrateStateRK4(F_0,F,Delta_t,co,ip,el) +function integrateStateRK4(F_0,F,Delta_t,co,ip,el) result(broken) real(pReal), intent(in),dimension(3,3) :: F_0,F real(pReal), intent(in) :: Delta_t integer, intent(in) :: co,ip,el + logical :: broken real(pReal), dimension(3,3), parameter :: & A = reshape([& @@ -1153,19 +1145,20 @@ subroutine integrateStateRK4(F_0,F,Delta_t,co,ip,el) real(pReal), dimension(4), parameter :: & B = [1.0_pReal/6.0_pReal, 1.0_pReal/3.0_pReal, 1.0_pReal/3.0_pReal, 1.0_pReal/6.0_pReal] - call integrateStateRK(F_0,F,Delta_t,co,ip,el,A,B,C) + broken = integrateStateRK(F_0,F,Delta_t,co,ip,el,A,B,C) -end subroutine integrateStateRK4 +end function integrateStateRK4 !--------------------------------------------------------------------------------------------------- !> @brief Integrate state (including stress integration) with the Cash-Carp method !--------------------------------------------------------------------------------------------------- -subroutine integrateStateRKCK45(F_0,F,Delta_t,co,ip,el) +function integrateStateRKCK45(F_0,F,Delta_t,co,ip,el) result(broken) real(pReal), intent(in),dimension(3,3) :: F_0,F real(pReal), intent(in) :: Delta_t integer, intent(in) :: co,ip,el + logical :: broken real(pReal), dimension(5,5), parameter :: & A = reshape([& @@ -1185,16 +1178,16 @@ subroutine integrateStateRKCK45(F_0,F,Delta_t,co,ip,el) [2825.0_pReal/27648.0_pReal, .0_pReal, 18575.0_pReal/48384.0_pReal,& 13525.0_pReal/55296.0_pReal, 277.0_pReal/14336.0_pReal, 1._pReal/4._pReal] - call integrateStateRK(F_0,F,Delta_t,co,ip,el,A,B,C,DB) + broken = integrateStateRK(F_0,F,Delta_t,co,ip,el,A,B,C,DB) -end subroutine integrateStateRKCK45 +end function integrateStateRKCK45 !-------------------------------------------------------------------------------------------------- !> @brief Integrate state (including stress integration) with an explicit Runge-Kutta method or an !! embedded explicit Runge-Kutta method !-------------------------------------------------------------------------------------------------- -subroutine integrateStateRK(F_0,F,Delta_t,co,ip,el,A,B,C,DB) +function integrateStateRK(F_0,F,Delta_t,co,ip,el,A,B,C,DB) result(broken) real(pReal), intent(in),dimension(3,3) :: F_0,F real(pReal), intent(in) :: Delta_t @@ -1205,15 +1198,14 @@ subroutine integrateStateRK(F_0,F,Delta_t,co,ip,el,A,B,C,DB) el, & !< element index in element loop ip, & !< integration point index in ip loop co !< grain index in grain loop + logical :: broken - integer :: & + integer :: & stage, & ! stage index in integration stage loop n, & ph, & me, & sizeDotState - logical :: & - broken real(pReal), dimension(constitutive_plasticity_maxSizeDotState,size(B)) :: plastic_RKdotState @@ -1266,10 +1258,8 @@ subroutine integrateStateRK(F_0,F,Delta_t,co,ip,el,A,B,C,DB) if(broken) return broken = integrateStress(F,Delta_t,co,ip,el) - crystallite_converged(co,ip,el) = .not. broken - -end subroutine integrateStateRK +end function integrateStateRK !-------------------------------------------------------------------------------------------------- @@ -1479,15 +1469,14 @@ end function constitutive_homogenizedC !-------------------------------------------------------------------------------------------------- !> @brief calculate stress (P) !-------------------------------------------------------------------------------------------------- -module function crystallite_stress(dt,co,ip,el) +module function crystallite_stress(dt,co,ip,el) result(converged_) real(pReal), intent(in) :: dt integer, intent(in) :: & co, & ip, & el - - logical :: crystallite_stress + logical :: converged_ real(pReal) :: & formerSubStep @@ -1519,7 +1508,7 @@ module function crystallite_stress(dt,co,ip,el) subFrac = 0.0_pReal subStep = 1.0_pReal/num%subStepSizeCryst todo = .true. - crystallite_converged(co,ip,el) = .false. ! pretend failed step of 1/subStepSizeCryst + converged_ = .false. ! pretend failed step of 1/subStepSizeCryst todo = .true. NiterationCrystallite = 0 @@ -1528,7 +1517,7 @@ module function crystallite_stress(dt,co,ip,el) !-------------------------------------------------------------------------------------------------- ! wind forward - if (crystallite_converged(co,ip,el)) then + if (converged_) then formerSubStep = subStep subFrac = subFrac + subStep subStep = min(1.0_pReal - subFrac, num%stepIncreaseCryst * subStep) @@ -1579,17 +1568,13 @@ module function crystallite_stress(dt,co,ip,el) math_inv33(matmul(constitutive_mech_Fi(ph)%data(1:3,1:3,me), & constitutive_mech_Fp(ph)%data(1:3,1:3,me)))) crystallite_subdt(co,ip,el) = subStep * dt - crystallite_converged(co,ip,el) = .false. - call integrateState(subF0,crystallite_subF(1:3,1:3,co,ip,el),& - crystallite_subdt(co,ip,el),co,ip,el) - call integrateSourceState(co,ip,el) + converged_ = .not. integrateState(subF0,crystallite_subF(1:3,1:3,co,ip,el),& + crystallite_subdt(co,ip,el),co,ip,el) + converged_ = converged_ .and. .not. integrateSourceState(co,ip,el) endif enddo cutbackLooping -! return whether converged or not - crystallite_stress = crystallite_converged(co,ip,el) - end function crystallite_stress end submodule constitutive_mech diff --git a/src/grid/discretization_grid.f90 b/src/grid/discretization_grid.f90 index 1b3700c14..48ad5b7e1 100644 --- a/src/grid/discretization_grid.f90 +++ b/src/grid/discretization_grid.f90 @@ -19,7 +19,6 @@ module discretization_grid use results use discretization use geometry_plastic_nonlocal - use FEsolving implicit none private @@ -117,9 +116,6 @@ subroutine discretization_grid_init(restart) (grid(1)+1) * (grid(2)+1) * grid3,& ! ...unless not last process worldrank+1==worldsize)) - FEsolving_execElem = [1,product(myGrid)] ! parallel loop bounds set to comprise all elements - FEsolving_execIP = [1,1] ! parallel loop bounds set to comprise the only IP - !-------------------------------------------------------------------------------------------------- ! store geometry information for post processing if(.not. restart) then diff --git a/src/grid/grid_mech_FEM.f90 b/src/grid/grid_mech_FEM.f90 index cdf806b35..003f568c6 100644 --- a/src/grid/grid_mech_FEM.f90 +++ b/src/grid/grid_mech_FEM.f90 @@ -18,7 +18,6 @@ module grid_mech_FEM use math use rotations use spectral_utilities - use FEsolving use config use homogenization use discretization diff --git a/src/grid/grid_mech_spectral_basic.f90 b/src/grid/grid_mech_spectral_basic.f90 index ebaaf3b55..9bc36165f 100644 --- a/src/grid/grid_mech_spectral_basic.f90 +++ b/src/grid/grid_mech_spectral_basic.f90 @@ -18,7 +18,6 @@ module grid_mech_spectral_basic use math use rotations use spectral_utilities - use FEsolving use config use homogenization use discretization_grid diff --git a/src/grid/grid_mech_spectral_polarisation.f90 b/src/grid/grid_mech_spectral_polarisation.f90 index 9f2a17c97..7160c1adc 100644 --- a/src/grid/grid_mech_spectral_polarisation.f90 +++ b/src/grid/grid_mech_spectral_polarisation.f90 @@ -18,7 +18,6 @@ module grid_mech_spectral_polarisation use math use rotations use spectral_utilities - use FEsolving use config use homogenization use discretization_grid diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index c0c84233d..e8bae223a 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -810,9 +810,9 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,& print'(/,a)', ' ... evaluating constitutive response ......................................' flush(IO_STDOUT) - homogenization_F = reshape(F,[3,3,product(grid(1:2))*grid3]) ! set materialpoint target F to estimated field + homogenization_F = reshape(F,[3,3,product(grid(1:2))*grid3]) ! set materialpoint target F to estimated field - call materialpoint_stressAndItsTangent(timeinc) ! calculate P field + call materialpoint_stressAndItsTangent(timeinc,[1,1],[1,product(grid(1:2))*grid3]) ! calculate P field P = reshape(homogenization_P, [3,3,grid(1),grid(2),grid3]) P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt ! average of P diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 27fdb6064..ebf5fd50d 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -11,7 +11,6 @@ module homogenization use math use material use constitutive - use FEsolving use discretization use thermal_isothermal use thermal_conduction @@ -144,27 +143,29 @@ end subroutine homogenization_init !-------------------------------------------------------------------------------------------------- !> @brief parallelized calculation of stress and corresponding tangent at material points !-------------------------------------------------------------------------------------------------- -subroutine materialpoint_stressAndItsTangent(dt) +subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execElem) real(pReal), intent(in) :: dt !< time increment + integer, dimension(2), intent(in) :: FEsolving_execElem, FEsolving_execIP integer :: & NiterationHomog, & NiterationMPstate, & ip, & !< integration point number el, & !< element number - myNgrains, co, ce + myNgrains, co, ce, ho real(pReal) :: & subFrac, & subStep logical :: & - requested, & converged logical, dimension(2) :: & doneAndHappy -!$OMP PARALLEL DO PRIVATE(ce,myNgrains,NiterationMPstate,NiterationHomog,subFrac,converged,subStep,requested,doneAndHappy) +!$OMP PARALLEL DO PRIVATE(ce,ho,myNgrains,NiterationMPstate,NiterationHomog,subFrac,converged,subStep,doneAndHappy) do el = FEsolving_execElem(1),FEsolving_execElem(2) + ho = material_homogenizationAt(el) + myNgrains = homogenization_Nconstituents(ho) do ip = FEsolving_execIP(1),FEsolving_execIP(2) !-------------------------------------------------------------------------------------------------- @@ -174,21 +175,19 @@ subroutine materialpoint_stressAndItsTangent(dt) subFrac = 0.0_pReal converged = .false. ! pretend failed step ... subStep = 1.0_pReal/num%subStepSizeHomog ! ... larger then the requested calculation - requested = .true. ! everybody requires calculation - if (homogState(material_homogenizationAt(el))%sizeState > 0) & - homogState(material_homogenizationAt(el))%subState0(:,material_homogenizationMemberAt(ip,el)) = & - homogState(material_homogenizationAt(el))%State0( :,material_homogenizationMemberAt(ip,el)) + if (homogState(ho)%sizeState > 0) & + homogState(ho)%subState0(:,material_homogenizationMemberAt(ip,el)) = & + homogState(ho)%State0( :,material_homogenizationMemberAt(ip,el)) + + if (damageState(ho)%sizeState > 0) & + damageState(ho)%subState0(:,material_homogenizationMemberAt(ip,el)) = & + damageState(ho)%State0( :,material_homogenizationMemberAt(ip,el)) - if (damageState(material_homogenizationAt(el))%sizeState > 0) & - damageState(material_homogenizationAt(el))%subState0(:,material_homogenizationMemberAt(ip,el)) = & - damageState(material_homogenizationAt(el))%State0( :,material_homogenizationMemberAt(ip,el)) NiterationHomog = 0 cutBackLooping: do while (.not. terminallyIll .and. subStep > num%subStepMinHomog) - myNgrains = homogenization_Nconstituents(material_homogenizationAt(el)) - if (converged) then subFrac = subFrac + subStep subStep = min(1.0_pReal-subFrac,num%stepIncreaseHomog*subStep) ! introduce flexibility for step increase/acceleration @@ -198,22 +197,20 @@ subroutine materialpoint_stressAndItsTangent(dt) ! wind forward grain starting point call constitutive_windForward(ip,el) - if(homogState(material_homogenizationAt(el))%sizeState > 0) & - homogState(material_homogenizationAt(el))%subState0(:,material_homogenizationMemberAt(ip,el)) = & - homogState(material_homogenizationAt(el))%State (:,material_homogenizationMemberAt(ip,el)) - if(damageState(material_homogenizationAt(el))%sizeState > 0) & - damageState(material_homogenizationAt(el))%subState0(:,material_homogenizationMemberAt(ip,el)) = & - damageState(material_homogenizationAt(el))%State (:,material_homogenizationMemberAt(ip,el)) + if(homogState(ho)%sizeState > 0) & + homogState(ho)%subState0(:,material_homogenizationMemberAt(ip,el)) = & + homogState(ho)%State (:,material_homogenizationMemberAt(ip,el)) + if(damageState(ho)%sizeState > 0) & + damageState(ho)%subState0(:,material_homogenizationMemberAt(ip,el)) = & + damageState(ho)%State (:,material_homogenizationMemberAt(ip,el)) endif steppingNeeded - else if ( (myNgrains == 1 .and. subStep <= 1.0 ) .or. & ! single grain already tried internal subStepping in crystallite num%subStepSizeHomog * subStep <= num%subStepMinHomog ) then ! would require too small subStep ! cutback makes no sense - if (.not. terminallyIll) then ! so first signals terminally ill... + if (.not. terminallyIll) & ! so first signals terminally ill... print*, ' Integration point ', ip,' at element ', el, ' terminally ill' - endif terminallyIll = .true. ! ...and kills all others else ! cutback makes sense subStep = num%subStepSizeHomog * subStep ! crystallite had severe trouble, so do a significant cutback @@ -221,23 +218,19 @@ subroutine materialpoint_stressAndItsTangent(dt) call crystallite_restore(ip,el,subStep < 1.0_pReal) call constitutive_restore(ip,el) - if(homogState(material_homogenizationAt(el))%sizeState > 0) & - homogState(material_homogenizationAt(el))%State( :,material_homogenizationMemberAt(ip,el)) = & - homogState(material_homogenizationAt(el))%subState0(:,material_homogenizationMemberAt(ip,el)) - if(damageState(material_homogenizationAt(el))%sizeState > 0) & - damageState(material_homogenizationAt(el))%State( :,material_homogenizationMemberAt(ip,el)) = & - damageState(material_homogenizationAt(el))%subState0(:,material_homogenizationMemberAt(ip,el)) + if(homogState(ho)%sizeState > 0) & + homogState(ho)%State( :,material_homogenizationMemberAt(ip,el)) = & + homogState(ho)%subState0(:,material_homogenizationMemberAt(ip,el)) + if(damageState(ho)%sizeState > 0) & + damageState(ho)%State( :,material_homogenizationMemberAt(ip,el)) = & + damageState(ho)%subState0(:,material_homogenizationMemberAt(ip,el)) endif endif - if (subStep > num%subStepMinHomog) then - requested = .true. - doneAndHappy = [.false.,.true.] - endif - + if (subStep > num%subStepMinHomog) doneAndHappy = [.false.,.true.] NiterationMPstate = 0 - convergenceLooping: do while (.not. terminallyIll .and. requested & + convergenceLooping: do while (.not. terminallyIll & .and. .not. doneAndHappy(1) & .and. NiterationMPstate < num%nMPstate) NiterationMPstate = NiterationMPstate + 1 @@ -245,7 +238,7 @@ subroutine materialpoint_stressAndItsTangent(dt) !-------------------------------------------------------------------------------------------------- ! deformation partitioning - if(requested .and. .not. doneAndHappy(1)) then ! requested but not yet done + if (.not. doneAndHappy(1)) then ce = (el-1)*discretization_nIPs + ip call mech_partition(homogenization_F0(1:3,1:3,ce) & + (homogenization_F(1:3,1:3,ce)-homogenization_F0(1:3,1:3,ce))& @@ -255,10 +248,7 @@ subroutine materialpoint_stressAndItsTangent(dt) do co = 1, myNgrains converged = converged .and. crystallite_stress(dt*subStep,co,ip,el) enddo - endif - - if (requested .and. .not. doneAndHappy(1)) then if (.not. converged) then doneAndHappy = [.true.,.false.] else @@ -281,10 +271,14 @@ subroutine materialpoint_stressAndItsTangent(dt) !$OMP END PARALLEL DO if (.not. terminallyIll ) then - call crystallite_orientations() ! calculate crystal orientations - !$OMP PARALLEL DO + !$OMP PARALLEL DO PRIVATE(ho,myNgrains) elementLooping3: do el = FEsolving_execElem(1),FEsolving_execElem(2) + ho = material_homogenizationAt(el) + myNgrains = homogenization_Nconstituents(ho) IpLooping3: do ip = FEsolving_execIP(1),FEsolving_execIP(2) + do co = 1, myNgrains + call crystallite_orientations(co,ip,el) + enddo call mech_homogenize(ip,el) enddo IpLooping3 enddo elementLooping3 diff --git a/src/homogenization_mech.f90 b/src/homogenization_mech.f90 index 56f1e554f..e4499e9b7 100644 --- a/src/homogenization_mech.f90 +++ b/src/homogenization_mech.f90 @@ -128,35 +128,35 @@ module subroutine mech_homogenize(ip,el) integer, intent(in) :: & ip, & !< integration point el !< element number - integer :: c,m + integer :: co,ce real(pReal) :: dPdFs(3,3,3,3,homogenization_Nconstituents(material_homogenizationAt(el))) - m = (el-1)* discretization_nIPs + ip + ce = (el-1)* discretization_nIPs + ip chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el))) case (HOMOGENIZATION_NONE_ID) chosenHomogenization - homogenization_P(1:3,1:3,m) = crystallite_P(1:3,1:3,1,ip,el) - homogenization_dPdF(1:3,1:3,1:3,1:3,m) = crystallite_stressTangent(1,ip,el) + homogenization_P(1:3,1:3,ce) = crystallite_P(1:3,1:3,1,ip,el) + homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = crystallite_stressTangent(1,ip,el) case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization - do c = 1, homogenization_Nconstituents(material_homogenizationAt(el)) - dPdFs(:,:,:,:,c) = crystallite_stressTangent(c,ip,el) + do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) + dPdFs(:,:,:,:,co) = crystallite_stressTangent(co,ip,el) enddo call mech_isostrain_averageStressAndItsTangent(& - homogenization_P(1:3,1:3,m), & - homogenization_dPdF(1:3,1:3,1:3,1:3,m),& + homogenization_P(1:3,1:3,ce), & + homogenization_dPdF(1:3,1:3,1:3,1:3,ce),& crystallite_P(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), & dPdFs, & homogenization_typeInstance(material_homogenizationAt(el))) case (HOMOGENIZATION_RGC_ID) chosenHomogenization - do c = 1, homogenization_Nconstituents(material_homogenizationAt(el)) - dPdFs(:,:,:,:,c) = crystallite_stressTangent(c,ip,el) + do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) + dPdFs(:,:,:,:,co) = crystallite_stressTangent(co,ip,el) enddo call mech_RGC_averageStressAndItsTangent(& - homogenization_P(1:3,1:3,m), & - homogenization_dPdF(1:3,1:3,1:3,1:3,m),& + homogenization_P(1:3,1:3,ce), & + homogenization_dPdF(1:3,1:3,1:3,1:3,ce),& crystallite_P(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), & dPdFs, & homogenization_typeInstance(material_homogenizationAt(el))) diff --git a/src/marc/discretization_marc.f90 b/src/marc/discretization_marc.f90 index ca0b54b73..675e57bd3 100644 --- a/src/marc/discretization_marc.f90 +++ b/src/marc/discretization_marc.f90 @@ -12,7 +12,6 @@ module discretization_marc use DAMASK_interface use IO use config - use FEsolving use element use discretization use geometry_plastic_nonlocal @@ -89,9 +88,6 @@ subroutine discretization_marc_init if (debug_e < 1 .or. debug_e > nElems) call IO_error(602,ext_msg='element') if (debug_i < 1 .or. debug_i > elem%nIPs) call IO_error(602,ext_msg='IP') - FEsolving_execElem = [1,nElems] - FEsolving_execIP = [1,elem%nIPs] - allocate(cellNodeDefinition(elem%nNodes-1)) allocate(connectivity_cell(elem%NcellNodesPerCell,elem%nIPs,nElems)) call buildCells(connectivity_cell,cellNodeDefinition,& diff --git a/src/mesh/DAMASK_mesh.f90 b/src/mesh/DAMASK_mesh.f90 index 1e353892c..7369520c1 100644 --- a/src/mesh/DAMASK_mesh.f90 +++ b/src/mesh/DAMASK_mesh.f90 @@ -15,7 +15,6 @@ program DAMASK_mesh use IO use math use CPFEM2 - use FEsolving use config use discretization_mesh use FEM_Utilities diff --git a/src/mesh/FEM_utilities.f90 b/src/mesh/FEM_utilities.f90 index cb81f1f0c..2f3633e11 100644 --- a/src/mesh/FEM_utilities.f90 +++ b/src/mesh/FEM_utilities.f90 @@ -160,7 +160,7 @@ subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData) print'(/,a)', ' ... evaluating constitutive response ......................................' - call materialpoint_stressAndItsTangent(timeinc) ! calculate P field + call materialpoint_stressAndItsTangent(timeinc,[1,mesh_maxNips],[1,mesh_NcpElems]) ! calculate P field cutBack = .false. ! reset cutBack status diff --git a/src/mesh/discretization_mesh.f90 b/src/mesh/discretization_mesh.f90 index 7dbd05e46..21c5feace 100644 --- a/src/mesh/discretization_mesh.f90 +++ b/src/mesh/discretization_mesh.f90 @@ -18,7 +18,6 @@ module discretization_mesh use config use discretization use results - use FEsolving use FEM_quadrature use YAML_types use prec @@ -30,7 +29,7 @@ module discretization_mesh mesh_Nboundaries, & mesh_NcpElemsGlobal - integer :: & + integer, public, protected :: & mesh_NcpElems !< total number of CP elements in mesh !!!! BEGIN DEPRECATED !!!!! @@ -174,9 +173,6 @@ subroutine discretization_mesh_init(restart) if (debug_element < 1 .or. debug_element > mesh_NcpElems) call IO_error(602,ext_msg='element') if (debug_ip < 1 .or. debug_ip > mesh_maxNips) call IO_error(602,ext_msg='IP') - FEsolving_execElem = [1,mesh_NcpElems] ! parallel loop bounds set to comprise all DAMASK elements - FEsolving_execIP = [1,mesh_maxNips] - allocate(mesh_node0(3,mesh_Nnodes),source=0.0_pReal) call discretization_init(materialAt,& From b0e5936b7ad8892f8bad96ff6297c6dda2cbbf5f Mon Sep 17 00:00:00 2001 From: Test User Date: Wed, 6 Jan 2021 17:36:47 +0100 Subject: [PATCH 19/19] [skip ci] updated version information after successful test of v3.0.0-alpha2-153-gf8dd5df0c --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index cb10a9c4f..c6d828650 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v3.0.0-alpha2-101-gab2a4a987 +v3.0.0-alpha2-153-gf8dd5df0c