diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index dda2ed1b6..ce645dcf5 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -75,7 +75,6 @@ variables: # ------------ Defaults ---------------------------------------------- MSC: "$MSC2019" IntelMarc: "$IntelCompiler17_8" - IntelAbaqus: "$IntelCompiler16_4" HDF5Marc: "HDF5/1.10.5/Intel-17.8" # ++++++++++++ Documentation ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ Doxygen1_8_15: "Documentation/Doxygen/1.8.15" @@ -476,15 +475,6 @@ createTar: - release ################################################################################################### -AbaqusStd: - stage: createDocumentation - script: - - module load $IntelCompiler $MPICH_Intel $PETSc_MPICH_Intel $Doxygen - - $DAMASKROOT/PRIVATE/documenting/runDoxygen.sh $DAMASKROOT abaqus - except: - - master - - release - Marc: stage: createDocumentation script: diff --git a/DAMASK_prerequisites.sh b/DAMASK_prerequisites.sh index 17ed70823..b72f19b7a 100755 --- a/DAMASK_prerequisites.sh +++ b/DAMASK_prerequisites.sh @@ -122,9 +122,3 @@ done firstLevel "CMake" getDetails cmake --version -firstLevel "Abaqus" -cd installation/mods_Abaqus # to have the right environment file -for executable in abaqus abq2017 abq2018 abq2019; do - getDetails $executable 'information=all' -done -cd ../.. diff --git a/PRIVATE b/PRIVATE index 64432754c..9dc7065be 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 64432754ce3c590c882cf4987695539cee524ee8 +Subproject commit 9dc7065beedf2a097aa60a656bbfa52e55d7147c diff --git a/VERSION b/VERSION index 7e4465e78..fc29e0602 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.3-1624-g47109b90 +v2.0.3-1639-g81ae6686 diff --git a/env/CONFIG b/env/CONFIG index 8da4d5b96..1b3d801bc 100644 --- a/env/CONFIG +++ b/env/CONFIG @@ -3,5 +3,3 @@ set DAMASK_NUM_THREADS = 4 set MSC_ROOT = /opt/msc set MARC_VERSION = 2019 - -set ABAQUS_VERSION = 2019 diff --git a/env/DAMASK.csh b/env/DAMASK.csh index fbfbfd1a7..a669a4ea0 100644 --- a/env/DAMASK.csh +++ b/env/DAMASK.csh @@ -51,5 +51,4 @@ else setenv PYTHONPATH $DAMASK_ROOT/python:$PYTHONPATH endif setenv MSC_ROOT -setenv ABAQUS_VERSION setenv MARC_VERSION diff --git a/examples/ConfigFiles/debug.config b/examples/ConfigFiles/debug.config index 57349580a..4736c87e1 100644 --- a/examples/ConfigFiles/debug.config +++ b/examples/ConfigFiles/debug.config @@ -17,7 +17,6 @@ homogenization # homogenization_*.f90 possible values: ba CPFEM # CPFEM.f90 possible values: basic, extensive, selective spectral # DAMASK_spectral.f90 possible values: basic, fft, restart, divergence, rotation, petsc marc # MSC.MARC FEM solver possible values: basic -abaqus # ABAQUS FEM solver possible values: basic # # Parameters for selective element 1 # selected element for debugging (synonymous: "el", "e") diff --git a/examples/ConfigFiles/numerics.config b/examples/ConfigFiles/numerics.config index 3a654513e..ced4ec6f1 100644 --- a/examples/ConfigFiles/numerics.config +++ b/examples/ConfigFiles/numerics.config @@ -9,7 +9,7 @@ pert_method 1 # perturbation method (1 = forward, 2 = b integrator 1 # integration method (1 = Fixed Point Iteration, 2 = Euler, 3 = Adaptive Euler, 4 = classical 4th order Runge-Kutta, 5 = 5th order Runge-Kutta Cash-Karp) integratorStiffness 1 # integration method used for stiffness (1 = Fixed Point Iteration, 2 = Euler, 3 = Adaptive Euler, 4 = classical 4th order Runge-Kutta, 5 = 5th order Runge-Kutta Cash-Karp) unitlength 1 # physical length of one computational length unit -usepingpong 1 # use the ping pong (collect <-> calc) scheme (always off for Abaqus exp, must be on for Spectral Solver) +usepingpong 1 # use the ping pong (collect <-> calc) scheme ## crystallite numerical parameters ## nCryst 20 # crystallite loop limit (only for debugging info, loop limit is determined by "subStepMinCryst") diff --git a/installation/mods_Abaqus/abaqus_v6.env b/installation/mods_Abaqus/abaqus_v6.env deleted file mode 100644 index 83cc2ed33..000000000 --- a/installation/mods_Abaqus/abaqus_v6.env +++ /dev/null @@ -1,56 +0,0 @@ -# -# DAMASK Abaqus Environment File -# -# ------------------------------------ -# originally taken from Abaqus ver. 6.11.1 -# -# -# Linux (Opteron/EM64T) Settings: -# -# Compile and Link command for user subroutines. -# Compile_cpp and link_exe for Abaqus make utility. -# -import os, re, glob, driverUtils - -if False: - # use hdf5 compiler wrapper in $PATH - fortCmd = os.popen('h5fc -shlib -show').read().replace('\n','') # complicated way needed to pass in DAMASKVERSION string - link_sl += fortCmd.split()[1:] - fortCmd +=" -DDAMASK_HDF5" -else: - # Use the version in $PATH - fortCmd = "ifort" - -# -free to use free-format FORTRAN 90 syntax -# -O <0-3> optimization level -# -fpp use FORTRAN preprocessor on source code -# -fopenmp build with openMP support -# -w90 -w95 suppress messages about use of non-standard Fortran (previous version of abaqus_v6.env only) -# -WB turn a compile-time bounds check into a warning (previous version of abaqus_v6.env only) -# -mP2OPT_hpo_vec_divbyzero=F inofficial compiler switch, proposed by abaqus but highly dubios (previous version of abaqus_v6.env only) -# -ftz flush underflow to zero -# -diag-disable 5268 disable warnings about line length > 132 (only comments there anyway) -# -implicitnone assume no implicit types (e.g. i for integer) -# -standard-semantics sets standard (Fortran 2008) and some other conventions -# -assume nostd_mod_proc_name avoid problems with libraries compiled without that option -# -real-size 64 assume size of real to be 8 bytes, matches our definition of pReal - -compile_fortran = (fortCmd + " -c -fPIC -auto -shared-intel " + - "-I%I -free -O3 -fpp -fopenmp " + - "-ftz -diag-disable 5268 " + - "-implicitnone -standard-semantics " + - "-assume nostd_mod_proc_name " + - "-real-size 64 " + - '-DDAMASKVERSION=\\\"n/a\\\"') - -# Abaqus/CAE will generate an input file without parts and assemblies. -cae_no_parts_input_file=ON -# Both the Abaqus/Explicit packager and analysis are run in double precision. -double_precision=BOTH -# The user will not be asked whether old job files of the same name should be deleted. -ask_delete=OFF - -# usub_lib_dir='your_prefered_location/abqlib' - -# Remove the temporary names from the namespace -del fortCmd diff --git a/installation/mods_Abaqus/abaqus_v6_debug.env b/installation/mods_Abaqus/abaqus_v6_debug.env deleted file mode 100644 index 943d0d10e..000000000 --- a/installation/mods_Abaqus/abaqus_v6_debug.env +++ /dev/null @@ -1,61 +0,0 @@ -# -# DAMASK Abaqus Environment File -# -# ------------------------------------ -# originally taken from Abaqus ver. 6.11.1 -# -# -# Linux (Opteron/EM64T) Settings: -# -# Compile and Link command for user subroutines. -# Compile_cpp and link_exe for Abaqus make utility. -# -import os, re, glob, driverUtils - -if False: - # use hdf5 compiler wrapper in $PATH - fortCmd = os.popen('h5fc -shlib -show').read().replace('\n','') # complicated way needed to pass in DAMASKVERSION string - link_sl += fortCmd.split()[1:] - fortCmd +=" -DDAMASK_HDF5" -else: - # Use the version in $PATH - fortCmd = "ifort" - -# -free to use free-format FORTRAN 90 syntax -# -O <0-3> optimization level -# -fpp use FORTRAN preprocessor on source code -# -fopenmp build with openMP support -# -w90 -w95 suppress messages about use of non-standard Fortran (previous version of abaqus_v6.env only) -# -WB turn a compile-time bounds check into a warning (previous version of abaqus_v6.env only) -# -mP2OPT_hpo_vec_divbyzero=F inofficial compiler switch, proposed by abaqus but highly dubios (previous version of abaqus_v6.env only) -# -ftz flush underflow to zero -# -diag-disable 5268 disable warnings about line length > 132 (only comments there anyway) -# -implicitnone assume no implicit types (e.g. i for integer) -# -standard-semantics sets standard (Fortran 2008) and some other conventions -# -assume nostd_mod_proc_name avoid problems with libraries compiled without that option -# -real-size 64 assume size of real to be 8 bytes, matches our definition of pReal - -# 'check pointers' does not work - -compile_fortran = (fortCmd + " -c -fPIC -auto -shared-intel " + - "-I%I -free -O0 -fpp " + - "-ftz -diag-disable 5268 " + - "-implicitnone -standard-semantics " + - "-assume nostd_mod_proc_name " + - "-real-size 64 " + - "-check bounds,format,output_conversion,uninit " + - "-ftrapuv -fpe-all0 " + - "-g -traceback -gen-interfaces -fp-stack-check -fp-model strict " + - '-DDAMASKVERSION=\\\"n/a\\\"') - -# Abaqus/CAE will generate an input file without parts and assemblies. -cae_no_parts_input_file=ON -# Both the Abaqus/Explicit packager and analysis are run in double precision. -double_precision=BOTH -# The user will not be asked whether old job files of the same name should be deleted. -ask_delete=OFF - -# usub_lib_dir='your_prefered_location/abqlib' - -# Remove the temporary names from the namespace -del fortCmd diff --git a/python/damask/environment.py b/python/damask/environment.py index 9964e9be7..2fde2c329 100644 --- a/python/damask/environment.py +++ b/python/damask/environment.py @@ -20,6 +20,5 @@ class Environment(): for item in ['DAMASK_NUM_THREADS', 'MSC_ROOT', 'MARC_VERSION', - 'ABAQUS_VERSION', ]: self.options[item] = os.environ[item] if item in os.environ else None diff --git a/python/damask/solver/__init__.py b/python/damask/solver/__init__.py index 70f48eb26..a1939eb55 100644 --- a/python/damask/solver/__init__.py +++ b/python/damask/solver/__init__.py @@ -2,4 +2,3 @@ from .solver import Solver # noqa from .marc import Marc # noqa -from .abaqus import Abaqus # noqa diff --git a/python/damask/solver/abaqus.py b/python/damask/solver/abaqus.py deleted file mode 100644 index 4c5f820d7..000000000 --- a/python/damask/solver/abaqus.py +++ /dev/null @@ -1,36 +0,0 @@ -import subprocess - -from .solver import Solver -import damask - - -class Abaqus(Solver): - """Wrapper to run DAMASK with Abaqus.""" - - def __init__(self,version=damask.Environment().options['ABAQUS_VERSION']): - """ - Create a Abaqus solver object. - - Parameters - ---------- - version : integer - Abaqus version - - """ - self.solver = 'Abaqus' - try: - self.version = int(version) - except TypeError: - self.version = -1 - - def return_run_command(self,model): - try: - cmd = 'abq{}'.format(self.version) - subprocess.check_output([cmd,'information=release']) - except OSError: # link to abqXXX not existing - cmd = 'abaqus' - process = subprocess.Popen([cmd,'information=release'],stdout = subprocess.PIPE,stderr = subprocess.PIPE) - detectedVersion = int(process.stdout.readlines()[1].split()[1].decode('utf-8')) - if self.version not in [detectedVersion,-1]: - raise Exception('found Abaqus version {}, but requested {}'.format(detectedVersion,self.version)) - return '{} -job {} -user {}/src/DAMASK_abaqus interactive'.format(cmd,model,damask.Environment().rootDir()) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 2e4336698..75b39e377 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -4,37 +4,36 @@ if (CMAKE_Fortran_COMPILER_ID STREQUAL "GNU") SET_SOURCE_FILES_PROPERTIES("lattice.f90" PROPERTIES COMPILE_FLAGS "-ffree-line-length-240") endif() -file (GLOB damask-sources *.f90 *.c) +file(GLOB damask-sources *.f90 *.c) -# probably we should have subfolders for abaqus and MSC.Marc -list (FILTER damask-sources EXCLUDE REGEX ".*CPFEM.f90") -list (FILTER damask-sources EXCLUDE REGEX ".*DAMASK_marc.*.f90") -list (FILTER damask-sources EXCLUDE REGEX ".*mesh_marc.*.f90") -list (FILTER damask-sources EXCLUDE REGEX ".*mesh_abaqus.*.f90") -list (FILTER damask-sources EXCLUDE REGEX ".*commercialFEM_fileList.*.f90") +# probably we should have a subfolder for MSC.Marc +list(FILTER damask-sources EXCLUDE REGEX ".*CPFEM.f90") +list(FILTER damask-sources EXCLUDE REGEX ".*DAMASK_marc.*.f90") +list(FILTER damask-sources EXCLUDE REGEX ".*mesh_marc.*.f90") +list(FILTER damask-sources EXCLUDE REGEX ".*commercialFEM_fileList.*.f90") if (PROJECT_NAME STREQUAL "damask-grid") - list (FILTER damask-sources EXCLUDE REGEX ".*mesh_FEM.*.f90") - file (GLOB grid-sources grid/*.f90) + list(FILTER damask-sources EXCLUDE REGEX ".*mesh_FEM.*.f90") + file(GLOB grid-sources grid/*.f90) - if (NOT CMAKE_BUILD_TYPE STREQUAL "SYNTAXONLY") - add_executable (DAMASK_spectral ${damask-sources} ${grid-sources}) + if(NOT CMAKE_BUILD_TYPE STREQUAL "SYNTAXONLY") + add_executable(DAMASK_spectral ${damask-sources} ${grid-sources}) install (TARGETS DAMASK_spectral RUNTIME DESTINATION bin) - else() - add_library (DAMASK_spectral OBJECT ${damask-sources} ${grid-sources}) + else() + add_library(DAMASK_spectral OBJECT ${damask-sources} ${grid-sources}) exec_program (mktemp OUTPUT_VARIABLE nothing) exec_program (mktemp ARGS -d OUTPUT_VARIABLE black_hole) install (PROGRAMS ${nothing} DESTINATION ${black_hole}) - endif() + endif() elseif (PROJECT_NAME STREQUAL "damask-mesh") - list (FILTER damask-sources EXCLUDE REGEX ".*mesh_grid.*.f90") - file (GLOB mesh-sources mesh/*.f90) + list(FILTER damask-sources EXCLUDE REGEX ".*mesh_grid.*.f90") + file(GLOB mesh-sources mesh/*.f90) - add_executable (DAMASK_FEM ${damask-sources} ${mesh-sources}) - install (TARGETS DAMASK_FEM RUNTIME DESTINATION bin) + add_executable(DAMASK_FEM ${damask-sources} ${mesh-sources}) + install (TARGETS DAMASK_FEM RUNTIME DESTINATION bin) endif() diff --git a/src/CPFEM.f90 b/src/CPFEM.f90 index fd406219f..1f76e6c25 100644 --- a/src/CPFEM.f90 +++ b/src/CPFEM.f90 @@ -39,8 +39,7 @@ module CPFEM integer(pInt), public :: & cycleCounter = 0_pInt, & !< needs description theInc = -1_pInt, & !< needs description - lastLovl = 0_pInt, & !< lovl in previous call to marc hypela2 - lastStep = 0_pInt !< kstep in previous call to abaqus umat + lastLovl = 0_pInt !< lovl in previous call to marc hypela2 real(pReal), public :: & theTime = 0.0_pReal, & !< needs description theDelta = 0.0_pReal diff --git a/src/DAMASK_abaqus.f b/src/DAMASK_abaqus.f deleted file mode 100644 index 0f663dde3..000000000 --- a/src/DAMASK_abaqus.f +++ /dev/null @@ -1,380 +0,0 @@ -!-------------------------------------------------------------------------------------------------- -!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH -!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH -!> @author Koen Janssens, Paul Scherrer Institut -!> @author Arun Prakash, Fraunhofer IWM -!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH -!> @brief interfaces DAMASK with Abaqus/Standard -!> @details put the included file abaqus_v6.env in either your home or model directory, -!> it is a minimum Abaqus environment file containing all changes necessary to use the -!> DAMASK subroutine (see Abaqus documentation for more information on the use of abaqus_v6.env) -!> @details Abaqus subroutines used: -!> @details - UMAT -!> @details - DFLUX -!-------------------------------------------------------------------------------------------------- -#define Abaqus - -#include "prec.f90" - -module DAMASK_interface - - implicit none - private - character(len=4), dimension(2), parameter, public :: INPUTFILEEXTENSION = ['.pes','.inp'] - - public :: & - DAMASK_interface_init, & - getSolverJobName - -contains - -!-------------------------------------------------------------------------------------------------- -!> @brief reports and sets working directory -!-------------------------------------------------------------------------------------------------- -subroutine DAMASK_interface_init -#if __INTEL_COMPILER >= 1800 - use, intrinsic :: iso_fortran_env, only: & - compiler_version, & - compiler_options -#endif - use ifport, only: & - CHDIR - - implicit none - integer, dimension(8) :: & - dateAndTime - integer :: lenOutDir,ierr - character(len=256) :: wd - - write(6,'(/,a)') ' <<<+- DAMASK_abaqus init -+>>>' - - write(6,'(/,a)') ' Roters et al., Computational Materials Science 158:420–478, 2019' - write(6,'(a)') ' https://doi.org/10.1016/j.commatsci.2018.04.030' - - write(6,'(/,a)') ' Version: '//DAMASKVERSION - -! https://github.com/jeffhammond/HPCInfo/blob/master/docs/Preprocessor-Macros.md -#if __INTEL_COMPILER >= 1800 - write(6,'(/,a)') ' Compiled with: '//compiler_version() - write(6,'(a)') ' Compiler options: '//compiler_options() -#else - write(6,'(/,a,i4.4,a,i8.8)') ' Compiled with Intel fortran version :', __INTEL_COMPILER,& - ', build date :', __INTEL_COMPILER_BUILD_DATE -#endif - - write(6,'(/,a)') ' Compiled on: '//__DATE__//' at '//__TIME__ - - call date_and_time(values = dateAndTime) - write(6,'(/,a,2(i2.2,a),i4.4)') ' Date: ',dateAndTime(3),'/',dateAndTime(2),'/', dateAndTime(1) - write(6,'(a,2(i2.2,a),i2.2)') ' Time: ',dateAndTime(5),':', dateAndTime(6),':', dateAndTime(7) - - call getoutdir(wd, lenOutDir) - ierr = CHDIR(wd) - if (ierr /= 0) then - write(6,'(a20,a,a16)') ' working directory "',trim(wd),'" does not exist' - call quit(1) - endif - -end subroutine DAMASK_interface_init - - -!-------------------------------------------------------------------------------------------------- -!> @brief using Abaqus/Standard function to get solver job name -!-------------------------------------------------------------------------------------------------- -character(1024) function getSolverJobName() - - implicit none - integer :: lenJobName - - getSolverJobName='' - call getJobName(getSolverJobName, lenJobName) - -end function getSolverJobName - - -end module DAMASK_interface - - - - -#include "commercialFEM_fileList.f90" - -!-------------------------------------------------------------------------------------------------- -!> @brief This is the Abaqus std user subroutine for defining material behavior -!-------------------------------------------------------------------------------------------------- -subroutine UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,& - RPL,DDSDDT,DRPLDE,DRPLDT,STRAN,DSTRAN,& - TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,NDI,NSHR,NTENS,& - NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,CELENT,& - DFGRD0,DFGRD1,NOEL,NPT,KSLAY,KSPT,KSTEP,KINC) - use prec, only: & - pReal, & - pInt - use numerics, only: & -!$ DAMASK_NumThreadsInt, & - usePingPong - use FEsolving, only: & - calcMode, & - terminallyIll, & - symmetricSolver - use debug, only: & - debug_info, & - debug_reset, & - debug_levelBasic, & - debug_level, & - debug_abaqus - use mesh, only: & - mesh_unitlength, & - mesh_FEasCP, & - mesh_ipCoordinates - use CPFEM, only: & - CPFEM_general, & - CPFEM_init_done, & - CPFEM_initAll, & - CPFEM_CALCRESULTS, & - CPFEM_AGERESULTS, & - CPFEM_COLLECT, & - CPFEM_BACKUPJACOBIAN, & - cycleCounter, & - theInc, & - theTime, & - theDelta, & - lastIncConverged, & - outdatedByNewInc, & - outdatedFFN1, & - lastStep - - implicit none - integer(pInt), intent(in) :: & - nDi, & !< Number of direct stress components at this point - nShr, & !< Number of engineering shear stress components at this point - nTens, & !< Size of the stress or strain component array (NDI + NSHR) - nStatV, & !< Number of solution-dependent state variables - nProps, & !< User-defined number of material constants - noEl, & !< element number - nPt,& !< integration point number - kSlay, & !< layer number (shell elements etc.) - kSpt, & !< section point within the current layer - kStep, & !< step number - kInc !< increment number - character(len=80), intent(in) :: & - cmname !< uses-specified material name, left justified - real(pReal), intent(in) :: & - DTIME, & - TEMP, & - DTEMP, & - CELENT - real(pReal), dimension(1), intent(in) :: & - PREDEF, & - DPRED - real(pReal), dimension(2), intent(in) :: & - TIME !< step time/total time at beginning of the current increment - real(pReal), dimension(3), intent(in) :: & - COORDS - real(pReal), dimension(nTens), intent(in) :: & - STRAN, & !< total strains at beginning of the increment - DSTRAN !< strain increments - real(pReal), dimension(nProps), intent(in) :: & - PROPS - real(pReal), dimension(3,3), intent(in) :: & - DROT, & !< rotation increment matrix - DFGRD0, & !< F at beginning of increment - DFGRD1 !< F at end of increment - real(pReal), intent(inout) :: & - PNEWDT, & !< ratio of suggested new time increment - SSE, & !< specific elastic strain engergy - SPD, & !< specific plastic dissipation - SCD, & !< specific creep dissipation - RPL, & !< volumetric heat generation per unit time at the end of the increment - DRPLDT !< varation of RPL with respect to the temperature - real(pReal), dimension(nTens), intent(inout) :: & - STRESS !< stress tensor at the beginning of the increment, needs to be updated - real(pReal), dimension(nStatV), intent(inout) :: & - STATEV !< solution-dependent state variables - real(pReal), dimension(nTens), intent(out) :: & - DDSDDT, & - DRPLDE - real(pReal), dimension(nTens,nTens), intent(out) :: & - DDSDDE !< Jacobian matrix of the constitutive model - - real(pReal) :: temperature ! temp by Abaqus is intent(in) - real(pReal), dimension(6) :: stress_h - real(pReal), dimension(6,6) :: ddsdde_h - integer(pInt) :: computationMode, i, cp_en - logical :: cutBack - -#ifdef _OPENMP - integer :: defaultNumThreadsInt !< default value set by Abaqus - include "omp_lib.h" - defaultNumThreadsInt = omp_get_num_threads() ! remember number of threads set by Marc - call omp_set_num_threads(DAMASK_NumThreadsInt) ! set number of threads for parallel execution set by DAMASK_NUM_THREADS -#endif - - temperature = temp ! temp is intent(in) - DDSDDT = 0.0_pReal - DRPLDE = 0.0_pReal - - if (iand(debug_level(debug_abaqus),debug_levelBasic) /= 0 .and. noel == 1 .and. npt == 1) then - write(6,*) 'el',noel,'ip',npt - write(6,*) 'got kInc as',kInc - write(6,*) 'got dStran',dStran - flush(6) - endif - - if (.not. CPFEM_init_done) call CPFEM_initAll(noel,npt) - - computationMode = 0 - cp_en = mesh_FEasCP('elem',noel) - if (time(2) > theTime .or. kInc /= theInc) then ! reached convergence - terminallyIll = .false. - cycleCounter = -1 ! first calc step increments this to cycle = 0 - if (kInc == 1) then ! >> start of analysis << - lastIncConverged = .false. ! no Jacobian backup - outdatedByNewInc = .false. ! no aging of state - calcMode = .false. ! pretend last step was collection - write (6,'(i8,1x,i2,1x,a)') noel,npt,'<< UMAT >> start of analysis..!';flush(6) - else if (kInc - theInc > 1) then ! >> restart of broken analysis << - lastIncConverged = .false. ! no Jacobian backup - outdatedByNewInc = .false. ! no aging of state - calcMode = .true. ! pretend last step was calculation - write (6,'(i8,1x,i2,1x,a)') noel,npt,'<< UMAT >> restart of analysis..!';flush(6) - else ! >> just the next inc << - lastIncConverged = .true. ! request Jacobian backup - outdatedByNewInc = .true. ! request aging of state - calcMode = .true. ! assure last step was calculation - write (6,'(i8,1x,i2,1x,a)') noel,npt,'<< UMAT >> new increment..!';flush(6) - endif - else if ( dtime < theDelta ) then ! >> cutBack << - lastIncConverged = .false. ! no Jacobian backup - outdatedByNewInc = .false. ! no aging of state - terminallyIll = .false. - cycleCounter = -1 ! first calc step increments this to cycle = 0 - calcMode = .true. ! pretend last step was calculation - write(6,'(i8,1x,i2,1x,a)') noel,npt,'<< UMAT >> cutback detected..!';flush(6) - endif ! convergence treatment end - - - if (usePingPong) then - calcMode(npt,cp_en) = .not. calcMode(npt,cp_en) ! ping pong (calc <--> collect) - if (calcMode(npt,cp_en)) then ! now --- CALC --- - computationMode = CPFEM_CALCRESULTS - if ( lastStep /= kStep ) then ! first after ping pong - call debug_reset() ! resets debugging - outdatedFFN1 = .false. - cycleCounter = cycleCounter + 1_pInt - endif - if(outdatedByNewInc) then - computationMode = ior(computationMode,CPFEM_AGERESULTS) ! calc and age results - outdatedByNewInc = .false. ! reset flag - endif - else ! now --- COLLECT --- - computationMode = CPFEM_COLLECT ! plain collect - if(lastStep /= kStep .and. .not. terminallyIll) & - call debug_info() ! first after ping pong reports (meaningful) debugging - if (lastIncConverged) then - computationMode = ior(computationMode,CPFEM_BACKUPJACOBIAN) ! collect and backup Jacobian after convergence - lastIncConverged = .false. ! reset flag - endif - mesh_ipCoordinates(1:3,npt,cp_en) = mesh_unitlength * COORDS - endif - else ! --- PLAIN MODE --- - computationMode = CPFEM_CALCRESULTS ! always calc - if (lastStep /= kStep) then - if (.not. terminallyIll) & - call debug_info() ! first reports (meaningful) debugging - call debug_reset() ! and resets debugging - outdatedFFN1 = .false. - cycleCounter = cycleCounter + 1_pInt - endif - if (outdatedByNewInc) then - computationMode = ior(computationMode,CPFEM_AGERESULTS) - outdatedByNewInc = .false. ! reset flag - endif - if (lastIncConverged) then - computationMode = ior(computationMode,CPFEM_BACKUPJACOBIAN) ! backup Jacobian after convergence - lastIncConverged = .false. ! reset flag - endif - endif - - - theTime = time(2) ! record current starting time - theDelta = dtime ! record current time increment - theInc = kInc ! record current increment number - lastStep = kStep ! record step number - - if (iand(debug_level(debug_abaqus),debug_levelBasic) /= 0) then - write(6,'(a16,1x,i2,1x,a,i8,a,i8,1x,i5,a)') 'computationMode',computationMode,'(',cp_en,':',noel,npt,')' - flush(6) - endif - - call CPFEM_general(computationMode,usePingPong,dfgrd0,dfgrd1,temperature,dtime,noel,npt,stress_h,ddsdde_h) - -! DAMASK: 11, 22, 33, 12, 23, 13 -! ABAQUS explicit: 11, 22, 33, 12, 23, 13 -! ABAQUS implicit: 11, 22, 33, 12, 13, 23 -! ABAQUS implicit: 11, 22, 33, 12 - - ddsdde = ddsdde_h(1:ntens,1:ntens) - stress = stress_h(1:ntens) - if(symmetricSolver) ddsdde = 0.5_pReal*(ddsdde + transpose(ddsdde)) - if(ntens == 6) then - stress_h = stress - stress(5) = stress_h(6) - stress(6) = stress_h(5) - ddsdde_h = ddsdde - ddsdde(:,5) = ddsdde_h(:,6) - ddsdde(:,6) = ddsdde_h(:,5) - ddsdde_h = ddsdde - ddsdde(5,:) = ddsdde_h(6,:) - ddsdde(6,:) = ddsdde_h(5,:) - end if - - statev = 0 - - if (terminallyIll) pnewdt = 0.5_pReal ! force cutback directly ? -!$ call omp_set_num_threads(defaultNumThreadsInt) ! reset number of threads to stored default value - -end subroutine UMAT - -!-------------------------------------------------------------------------------------------------- -!> @brief calculate internal heat generated due to inelastic energy dissipation -!-------------------------------------------------------------------------------------------------- -SUBROUTINE DFLUX(FLUX,SOL,KSTEP,KINC,TIME,NOEL,NPT,COORDS,& - JLTYP,TEMP,PRESS,SNAME) - use prec, only: & - pReal, & - pInt - use thermal_conduction, only: & - thermal_conduction_getSourceAndItsTangent - use mesh, only: & - mesh_FEasCP - - implicit none - character(len=1024) :: sname - real(pReal), dimension(2), intent(out) :: flux - real(pReal), intent(in) :: sol - Integer(pInt), intent(in) :: Kstep, Kinc, Noel, Npt - real(pReal), dimension(3), intent(in) :: coords - real(pReal), intent(in) :: time - real(pReal) :: Jltyp - real(pReal), intent(in) :: temp - real(pReal), intent(in) :: Press - - jltyp = 1 - call thermal_conduction_getSourceAndItsTangent(flux(1), flux(2), sol, npt,mesh_FEasCP('elem',noel)) -end subroutine DFLUX - -!-------------------------------------------------------------------------------------------------- -!> @brief calls the exit function of Abaqus/Standard -!-------------------------------------------------------------------------------------------------- -subroutine quit(DAMASK_error) - use prec, only: & - pInt - - implicit none - integer(pInt) :: DAMASK_error - - flush(6) - call xit - -end subroutine quit diff --git a/src/FEsolving.f90 b/src/FEsolving.f90 index fc04499fa..cce95a07e 100644 --- a/src/FEsolving.f90 +++ b/src/FEsolving.f90 @@ -15,7 +15,7 @@ module FEsolving 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 -#if defined(Marc4DAMASK) || defined(Abaqus) +#if defined(Marc4DAMASK) logical, dimension(:,:), allocatable :: & calcMode !< do calculation or simply collect when using ping pong scheme #endif diff --git a/src/IO.f90 b/src/IO.f90 index da3dcbf3c..380d907b2 100644 --- a/src/IO.f90 +++ b/src/IO.f90 @@ -33,16 +33,6 @@ module IO IO_lc, & IO_error, & IO_warning -#if defined(Marc4DAMASK) || defined(Abaqus) - public :: & - IO_open_inputFile -#if defined(Abaqus) - public :: & - IO_continuousIntValues, & - IO_extractValue, & - IO_countDataLines -#endif -#endif contains @@ -175,94 +165,6 @@ integer function IO_open_binary(fileName,mode) end function IO_open_binary -#if defined(Marc4DAMASK) || defined(Abaqus) -!-------------------------------------------------------------------------------------------------- -!> @brief opens FEM input file for reading located in current working directory to given unit -!-------------------------------------------------------------------------------------------------- -subroutine IO_open_inputFile(fileUnit) - - integer, intent(in) :: fileUnit !< file unit - - integer :: myStat - character(len=1024) :: path -#if defined(Abaqus) - integer :: fileType - - fileType = 1 ! assume .pes - path = trim(getSolverJobName())//inputFileExtension(fileType) ! attempt .pes, if it exists: it should be used - open(fileUnit+1,status='old',iostat=myStat,file=path,action='read',position='rewind') - if(myStat /= 0) then ! if .pes does not work / exist; use conventional extension, i.e.".inp" - fileType = 2 - path = trim(getSolverJobName())//inputFileExtension(fileType) - open(fileUnit+1,status='old',iostat=myStat,file=path,action='read',position='rewind') - endif - if (myStat /= 0) call IO_error(100,el=myStat,ext_msg=path) - - path = trim(getSolverJobName())//inputFileExtension(fileType)//'_assembly' - open(fileUnit,iostat=myStat,file=path) - if (myStat /= 0) call IO_error(100,el=myStat,ext_msg=path) - if (.not.abaqus_assembleInputFile(fileUnit,fileUnit+1)) call IO_error(103) ! strip comments and concatenate any "include"s - close(fileUnit+1) - - contains - - !-------------------------------------------------------------------------------------------------- - !> @brief create a new input file for abaqus simulations by removing all comment lines and - !> including "include"s - !-------------------------------------------------------------------------------------------------- - recursive function abaqus_assembleInputFile(unit1,unit2) result(createSuccess) - - integer, intent(in) :: unit1, & - unit2 - - - integer, allocatable, dimension(:) :: chunkPos - character(len=pStringLen :: line,fname - logical :: createSuccess,fexist - - - do - read(unit2,'(A)',END=220) line - chunkPos = IO_stringPos(line) - - if (IO_lc(IO_StringValue(line,chunkPos,1))=='*include') then - fname = trim(line(9+scan(line(9:),'='):)) - inquire(file=fname, exist=fexist) - if (.not.(fexist)) then - write(6,*)'ERROR: file does not exist error in abaqus_assembleInputFile' - write(6,*)'filename: ', trim(fname) - createSuccess = .false. - return - endif - open(unit2+1,err=200,status='old',file=fname) - if (abaqus_assembleInputFile(unit1,unit2+1)) then - createSuccess=.true. - close(unit2+1) - else - createSuccess=.false. - return - endif - else if (line(1:2) /= '**' .OR. line(1:8)=='**damask') then - write(unit1,'(A)') trim(line) - endif - enddo - -220 createSuccess = .true. - return - -200 createSuccess =.false. - - end function abaqus_assembleInputFile -#elif defined(Marc4DAMASK) - path = trim(getSolverJobName())//inputFileExtension - open(fileUnit,status='old',iostat=myStat,file=path) - if (myStat /= 0) call IO_error(100,el=myStat,ext_msg=path) -#endif - -end subroutine IO_open_inputFile -#endif - - !-------------------------------------------------------------------------------------------------- !> @brief identifies strings without content !-------------------------------------------------------------------------------------------------- @@ -419,7 +321,6 @@ end function IO_lc !-------------------------------------------------------------------------------------------------- !> @brief write error statements to standard out and terminate the Marc/spectral run with exit #9xxx -!> in ABAQUS either time step is reduced or execution terminated !-------------------------------------------------------------------------------------------------- subroutine IO_error(error_ID,el,ip,g,instance,ext_msg) @@ -608,29 +509,6 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg) case (894) msg = 'MPI error' -!------------------------------------------------------------------------------------------------- -! error messages related to parsing of Abaqus input file - case (900) - msg = 'improper definition of nodes in input file (Nnodes < 2)' - case (901) - msg = 'no elements defined in input file (Nelems = 0)' - case (902) - msg = 'no element sets defined in input file (No *Elset exists)' - case (903) - msg = 'no materials defined in input file (Look into section assigments)' - case (904) - msg = 'no elements could be assigned for Elset: ' - case (905) - msg = 'error in mesh_abaqus_map_materials' - case (906) - msg = 'error in mesh_abaqus_count_cpElements' - case (907) - msg = 'size of mesh_mapFEtoCPelem in mesh_abaqus_map_elements' - case (908) - msg = 'size of mesh_mapFEtoCPnode in mesh_abaqus_map_nodes' - case (909) - msg = 'size of mesh_node in mesh_abaqus_build_nodes not equal to mesh_Nnodes' - !------------------------------------------------------------------------------------------------- ! general error messages @@ -751,126 +629,6 @@ subroutine IO_warning(warning_ID,el,ip,g,ext_msg) end subroutine IO_warning -#ifdef Abaqus -!-------------------------------------------------------------------------------------------------- -!> @brief extracts string value from key=value pair and check whether key matches -!-------------------------------------------------------------------------------------------------- -character(len=300) pure function IO_extractValue(pair,key) - - character(len=*), intent(in) :: pair, & !< key=value pair - key !< key to be expected - - character(len=*), parameter :: SEP = achar(61) ! '=' - - integer :: myChunk !< position number of desired chunk - - IO_extractValue = '' - - myChunk = scan(pair,SEP) - if (myChunk > 0 .and. pair(:myChunk-1) == key) IO_extractValue = pair(myChunk+1:) ! extract value if key matches - -end function IO_extractValue - - -!-------------------------------------------------------------------------------------------------- -!> @brief count lines containig data up to next *keyword -!-------------------------------------------------------------------------------------------------- -integer function IO_countDataLines(fileUnit) - - integer, intent(in) :: fileUnit !< file handle - - - integer, allocatable, dimension(:) :: chunkPos - character(len=pStringLen) :: line, & - tmp - - IO_countDataLines = 0 - line = '' - - do while (trim(line) /= IO_EOF) - read(fileUnit,'(A)') line - chunkPos = IO_stringPos(line) - tmp = IO_lc(IO_stringValue(line,chunkPos,1)) - if (tmp(1:1) == '*' .and. tmp(2:2) /= '*') then ! found keyword - exit - else - if (tmp(2:2) /= '*') IO_countDataLines = IO_countDataLines + 1 - endif - enddo - backspace(fileUnit) - -end function IO_countDataLines - - -!-------------------------------------------------------------------------------------------------- -!> @brief return integer list corresponding to items in consecutive lines. -!! First integer in array is counter -!> @details Abaqus: triplet of start,stop,inc or named set -!-------------------------------------------------------------------------------------------------- -function IO_continuousIntValues(fileUnit,maxN,lookupName,lookupMap,lookupMaxN) - - integer, intent(in) :: maxN - integer, dimension(1+maxN) :: IO_continuousIntValues - - integer, intent(in) :: fileUnit, & - lookupMaxN - integer, dimension(:,:), intent(in) :: lookupMap - character(len=*), dimension(:), intent(in) :: lookupName - integer :: i,first,last - integer :: j,l,c - integer, allocatable, dimension(:) :: chunkPos - character(len=pStringLen) :: line - logical :: rangeGeneration - - IO_continuousIntValues = 0 - rangeGeneration = .false. - - c = IO_countDataLines(fileUnit) - do l = 1,c - backspace(fileUnit) - enddo - -!-------------------------------------------------------------------------------------------------- -! check if the element values in the elset are auto generated - backspace(fileUnit) - read(fileUnit,'(A)',end=100) line - chunkPos = IO_stringPos(line) - do i = 1,chunkPos(1) - if (IO_lc(IO_stringValue(line,chunkPos,i)) == 'generate') rangeGeneration = .true. - enddo - - do l = 1,c - read(fileUnit,'(A)',end=100) line - chunkPos = IO_stringPos(line) - if (verify(IO_stringValue(line,chunkPos,1),'0123456789') > 0) then ! a non-int, i.e. set names follow on this line - do i = 1,chunkPos(1) ! loop over set names in line - do j = 1,lookupMaxN ! look through known set names - if (IO_stringValue(line,chunkPos,i) == lookupName(j)) then ! found matching name - first = 2 + IO_continuousIntValues(1) ! where to start appending data - last = first + lookupMap(1,j) - 1 ! up to where to append data - IO_continuousIntValues(first:last) = lookupMap(2:1+lookupMap(1,j),j) ! add resp. entity list - IO_continuousIntValues(1) = IO_continuousIntValues(1) + lookupMap(1,j) ! count them - endif - enddo - enddo - else if (rangeGeneration) then ! range generation - do i = IO_intValue(line,chunkPos,1),& - IO_intValue(line,chunkPos,2),& - max(1,IO_intValue(line,chunkPos,3)) - IO_continuousIntValues(1) = IO_continuousIntValues(1) + 1 - IO_continuousIntValues(1+IO_continuousIntValues(1)) = i - enddo - else ! read individual elem nums - do i = 1,chunkPos(1) - IO_continuousIntValues(1) = IO_continuousIntValues(1) + 1 - IO_continuousIntValues(1+IO_continuousIntValues(1)) = IO_intValue(line,chunkPos,i) - enddo - endif - enddo - -100 end function IO_continuousIntValues -#endif - !-------------------------------------------------------------------------------------------------- ! internal helper functions diff --git a/src/commercialFEM_fileList.f90 b/src/commercialFEM_fileList.f90 index c22a37a0f..f3b681052 100644 --- a/src/commercialFEM_fileList.f90 +++ b/src/commercialFEM_fileList.f90 @@ -1,7 +1,7 @@ !-------------------------------------------------------------------------------------------------- !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @brief all DAMASK files without solver -!> @details List of files needed by MSC.Marc and Abaqus/Standard +!> @details List of files needed by MSC.Marc !-------------------------------------------------------------------------------------------------- #include "IO.f90" #include "numerics.f90" @@ -19,9 +19,6 @@ #include "results.f90" #include "geometry_plastic_nonlocal.f90" #include "discretization.f90" -#ifdef Abaqus -#include "mesh_abaqus.f90" -#endif #ifdef Marc4DAMASK #include "mesh_marc.f90" #endif diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 681691c9e..6f3c91fff 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -25,44 +25,44 @@ module constitutive use kinematics_cleavage_opening use kinematics_slipplane_opening use kinematics_thermal_expansion - + implicit none private - + integer, public, protected :: & constitutive_plasticity_maxSizeDotState, & constitutive_source_maxSizeDotState - + interface module subroutine plastic_none_init end subroutine plastic_none_init - + module subroutine plastic_isotropic_init end subroutine plastic_isotropic_init - + module subroutine plastic_phenopowerlaw_init end subroutine plastic_phenopowerlaw_init - + module subroutine plastic_kinehardening_init end subroutine plastic_kinehardening_init - + module subroutine plastic_dislotwin_init end subroutine plastic_dislotwin_init module subroutine plastic_disloUCLA_init end subroutine plastic_disloUCLA_init - + module subroutine plastic_nonlocal_init end subroutine plastic_nonlocal_init - - + + module subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) real(pReal), dimension(3,3), intent(out) :: & Lp !< plastic velocity gradient real(pReal), dimension(3,3,3,3), intent(out) :: & dLp_dMp !< derivative of Lp with respect to the Mandel stress - + real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & @@ -75,20 +75,20 @@ module constitutive Lp !< plastic velocity gradient real(pReal), dimension(3,3,3,3), intent(out) :: & dLp_dMp !< derivative of Lp with respect to the Mandel stress - + real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & instance, & of end subroutine plastic_phenopowerlaw_LpAndItsTangent - + pure module subroutine plastic_kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) real(pReal), dimension(3,3), intent(out) :: & Lp !< plastic velocity gradient real(pReal), dimension(3,3,3,3), intent(out) :: & dLp_dMp !< derivative of Lp with respect to the Mandel stress - + real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & @@ -101,7 +101,7 @@ module constitutive Lp !< plastic velocity gradient real(pReal), dimension(3,3,3,3), intent(out) :: & dLp_dMp !< derivative of Lp with respect to the Mandel stress - + real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress real(pReal), intent(in) :: & @@ -110,13 +110,13 @@ module constitutive instance, & of end subroutine plastic_dislotwin_LpAndItsTangent - + pure module subroutine plastic_disloUCLA_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of) real(pReal), dimension(3,3), intent(out) :: & Lp !< plastic velocity gradient real(pReal), dimension(3,3,3,3), intent(out) :: & dLp_dMp !< derivative of Lp with respect to the Mandel stress - + real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress real(pReal), intent(in) :: & @@ -125,14 +125,14 @@ module constitutive instance, & of end subroutine plastic_disloUCLA_LpAndItsTangent - + module subroutine plastic_nonlocal_LpAndItsTangent(Lp, dLp_dMp, & Mp, Temperature, volume, ip, el) real(pReal), dimension(3,3), intent(out) :: & Lp !< plastic velocity gradient real(pReal), dimension(3,3,3,3), intent(out) :: & dLp_dMp !< derivative of Lp with respect to the Mandel stress - + real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress real(pReal), intent(in) :: & @@ -149,15 +149,15 @@ module constitutive Li !< inleastic velocity gradient real(pReal), dimension(3,3,3,3), intent(out) :: & dLi_dMi !< derivative of Li with respect to Mandel stress - + real(pReal), dimension(3,3), intent(in) :: & - Mi !< Mandel stress + Mi !< Mandel stress integer, intent(in) :: & instance, & of end subroutine plastic_isotropic_LiAndItsTangent - - + + module subroutine plastic_isotropic_dotState(Mp,instance,of) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress @@ -165,7 +165,7 @@ module constitutive instance, & of end subroutine plastic_isotropic_dotState - + module subroutine plastic_phenopowerlaw_dotState(Mp,instance,of) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress @@ -173,7 +173,7 @@ module constitutive instance, & of end subroutine plastic_phenopowerlaw_dotState - + module subroutine plastic_kinehardening_dotState(Mp,instance,of) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress @@ -181,7 +181,7 @@ module constitutive instance, & of end subroutine plastic_kinehardening_dotState - + module subroutine plastic_dislotwin_dotState(Mp,T,instance,of) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress @@ -201,8 +201,8 @@ module constitutive instance, & of end subroutine plastic_disloUCLA_dotState - - module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & + + module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature, & timestep,ip,el) integer, intent(in) :: & ip, & !< current integration point @@ -213,11 +213,11 @@ module constitutive real(pReal), dimension(3,3), intent(in) ::& Mp !< MandelStress real(pReal), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem), intent(in) :: & - Fe, & !< elastic deformation gradient + F, & !< deformation gradient Fp !< plastic deformation gradient end subroutine plastic_nonlocal_dotState - + module subroutine plastic_dislotwin_dependentState(T,instance,of) integer, intent(in) :: & instance, & @@ -225,19 +225,19 @@ module constitutive real(pReal), intent(in) :: & T end subroutine plastic_dislotwin_dependentState - + module subroutine plastic_disloUCLA_dependentState(instance,of) integer, intent(in) :: & instance, & of end subroutine plastic_disloUCLA_dependentState - - module subroutine plastic_nonlocal_dependentState(Fe, Fp, ip, el) + + module subroutine plastic_nonlocal_dependentState(F, Fp, ip, el) integer, intent(in) :: & ip, & el real(pReal), dimension(3,3), intent(in) :: & - Fe, & + F, & Fp end subroutine plastic_nonlocal_dependentState @@ -249,7 +249,7 @@ module constitutive instance, & of end subroutine plastic_kinehardening_deltaState - + module subroutine plastic_nonlocal_deltaState(Mp,ip,el) integer, intent(in) :: & ip, & @@ -267,48 +267,48 @@ module constitutive ip, & !< integration point el !< element end function plastic_dislotwin_homogenizedC - - module subroutine plastic_nonlocal_updateCompatibility(orientation,i,e) + + module subroutine plastic_nonlocal_updateCompatibility(orientation,i,e) integer, intent(in) :: & i, & e type(rotation), dimension(1,discretization_nIP,discretization_nElem), intent(in) :: & orientation !< crystal orientation end subroutine plastic_nonlocal_updateCompatibility - - + + module subroutine plastic_isotropic_results(instance,group) integer, intent(in) :: instance character(len=*), intent(in) :: group end subroutine plastic_isotropic_results - + module subroutine plastic_phenopowerlaw_results(instance,group) integer, intent(in) :: instance character(len=*), intent(in) :: group end subroutine plastic_phenopowerlaw_results - + module subroutine plastic_kinehardening_results(instance,group) integer, intent(in) :: instance character(len=*), intent(in) :: group end subroutine plastic_kinehardening_results - + module subroutine plastic_dislotwin_results(instance,group) integer, intent(in) :: instance character(len=*), intent(in) :: group end subroutine plastic_dislotwin_results - + module subroutine plastic_disloUCLA_results(instance,group) integer, intent(in) :: instance character(len=*), intent(in) :: group end subroutine plastic_disloUCLA_results - + module subroutine plastic_nonlocal_results(instance,group) integer, intent(in) :: instance character(len=*), intent(in) :: group end subroutine plastic_nonlocal_results - + end interface - + public :: & plastic_nonlocal_updateCompatibility, & constitutive_init, & @@ -321,7 +321,7 @@ module constitutive constitutive_collectDotState, & constitutive_collectDeltaState, & constitutive_results - + contains @@ -355,18 +355,18 @@ subroutine constitutive_init if (any(phase_source == SOURCE_damage_isoDuctile_ID)) call source_damage_isoDuctile_init if (any(phase_source == SOURCE_damage_anisoBrittle_ID)) call source_damage_anisoBrittle_init if (any(phase_source == SOURCE_damage_anisoDuctile_ID)) call source_damage_anisoDuctile_init - + !-------------------------------------------------------------------------------------------------- ! initialize kinematic mechanisms if (any(phase_kinematics == KINEMATICS_cleavage_opening_ID)) call kinematics_cleavage_opening_init if (any(phase_kinematics == KINEMATICS_slipplane_opening_ID)) call kinematics_slipplane_opening_init if (any(phase_kinematics == KINEMATICS_thermal_expansion_ID)) call kinematics_thermal_expansion_init - + write(6,'(/,a)') ' <<<+- constitutive init -+>>>'; flush(6) - + constitutive_plasticity_maxSizeDotState = 0 constitutive_source_maxSizeDotState = 0 - + PhaseLoop2:do ph = 1,material_Nphase !-------------------------------------------------------------------------------------------------- ! partition and inititalize state @@ -398,7 +398,7 @@ function constitutive_homogenizedC(ipc,ip,el) ipc, & !< component-ID of integration point ip, & !< integration point el !< element - + plasticityType: select case (phase_plasticity(material_phaseAt(ipc,el))) case (PLASTICITY_DISLOTWIN_ID) plasticityType constitutive_homogenizedC = plastic_dislotwin_homogenizedC(ipc,ip,el) @@ -412,23 +412,23 @@ end function constitutive_homogenizedC !-------------------------------------------------------------------------------------------------- !> @brief calls microstructure function of the different constitutive models !-------------------------------------------------------------------------------------------------- -subroutine constitutive_microstructure(Fe, Fp, ipc, ip, el) +subroutine constitutive_microstructure(F, Fp, ipc, ip, el) integer, intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point el !< element real(pReal), intent(in), dimension(3,3) :: & - Fe, & !< elastic deformation gradient + F, & !< elastic deformation gradient Fp !< plastic deformation gradient integer :: & ho, & !< homogenization tme, & !< thermal member position instance, of - + ho = material_homogenizationAt(el) tme = thermalMapping(ho)%p(ip,el) - + plasticityType: select case (phase_plasticity(material_phaseAt(ipc,el))) case (PLASTICITY_DISLOTWIN_ID) plasticityType of = material_phasememberAt(ipc,ip,el) @@ -439,7 +439,7 @@ subroutine constitutive_microstructure(Fe, Fp, ipc, ip, el) instance = phase_plasticityInstance(material_phaseAt(ipc,el)) call plastic_disloUCLA_dependentState(instance,of) case (PLASTICITY_NONLOCAL_ID) plasticityType - call plastic_nonlocal_dependentState (Fe,Fp,ip,el) + call plastic_nonlocal_dependentState (F,Fp,ip,el) end select plasticityType end subroutine constitutive_microstructure @@ -451,7 +451,7 @@ end subroutine constitutive_microstructure ! Mp in, dLp_dMp out !-------------------------------------------------------------------------------------------------- subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & - S, Fi, ipc, ip, el) + S, Fi, ipc, ip, el) integer, intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point @@ -473,49 +473,49 @@ subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & tme !< thermal member position integer :: & i, j, instance, of - + ho = material_homogenizationAt(el) tme = thermalMapping(ho)%p(ip,el) - + Mp = matmul(matmul(transpose(Fi),Fi),S) - + plasticityType: select case (phase_plasticity(material_phaseAt(ipc,el))) - + case (PLASTICITY_NONE_ID) plasticityType Lp = 0.0_pReal dLp_dMp = 0.0_pReal - + case (PLASTICITY_ISOTROPIC_ID) plasticityType of = material_phasememberAt(ipc,ip,el) instance = phase_plasticityInstance(material_phaseAt(ipc,el)) call plastic_isotropic_LpAndItsTangent (Lp,dLp_dMp,Mp,instance,of) - + case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType of = material_phasememberAt(ipc,ip,el) instance = phase_plasticityInstance(material_phaseAt(ipc,el)) call plastic_phenopowerlaw_LpAndItsTangent (Lp,dLp_dMp,Mp,instance,of) - + case (PLASTICITY_KINEHARDENING_ID) plasticityType of = material_phasememberAt(ipc,ip,el) instance = phase_plasticityInstance(material_phaseAt(ipc,el)) call plastic_kinehardening_LpAndItsTangent (Lp,dLp_dMp, Mp,instance,of) - + case (PLASTICITY_NONLOCAL_ID) plasticityType call plastic_nonlocal_LpAndItsTangent (Lp,dLp_dMp,Mp, & temperature(ho)%p(tme),geometry_plastic_nonlocal_IPvolume0(ip,el),ip,el) - + case (PLASTICITY_DISLOTWIN_ID) plasticityType of = material_phasememberAt(ipc,ip,el) instance = phase_plasticityInstance(material_phaseAt(ipc,el)) call plastic_dislotwin_LpAndItsTangent (Lp,dLp_dMp,Mp,temperature(ho)%p(tme),instance,of) - + case (PLASTICITY_DISLOUCLA_ID) plasticityType of = material_phasememberAt(ipc,ip,el) instance = phase_plasticityInstance(material_phaseAt(ipc,el)) call plastic_disloucla_LpAndItsTangent (Lp,dLp_dMp,Mp,temperature(ho)%p(tme),instance,of) - + end select plasticityType - + do i=1,3; do j=1,3 dLp_dFi(i,j,1:3,1:3) = matmul(matmul(Fi,S),transpose(dLp_dMp(i,j,1:3,1:3))) + & matmul(matmul(Fi,dLp_dMp(i,j,1:3,1:3)),S) @@ -545,7 +545,7 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & real(pReal), intent(out), dimension(3,3,3,3) :: & dLi_dS, & !< derivative of Li with respect to S dLi_dFi - + real(pReal), dimension(3,3) :: & my_Li, & !< intermediate velocity gradient FiInv, & @@ -557,11 +557,11 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & integer :: & k, i, j, & instance, of - + Li = 0.0_pReal dLi_dS = 0.0_pReal dLi_dFi = 0.0_pReal - + plasticityType: select case (phase_plasticity(material_phaseAt(ipc,el))) case (PLASTICITY_isotropic_ID) plasticityType of = material_phasememberAt(ipc,ip,el) @@ -571,10 +571,10 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & my_Li = 0.0_pReal my_dLi_dS = 0.0_pReal end select plasticityType - + Li = Li + my_Li dLi_dS = dLi_dS + my_dLi_dS - + KinematicsLoop: do k = 1, phase_Nkinematics(material_phaseAt(ipc,el)) kinematicsType: select case (phase_kinematics(k,material_phaseAt(ipc,el))) case (KINEMATICS_cleavage_opening_ID) kinematicsType @@ -590,12 +590,12 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & Li = Li + my_Li dLi_dS = dLi_dS + my_dLi_dS enddo KinematicsLoop - + FiInv = math_inv33(Fi) detFi = math_det33(Fi) Li = matmul(matmul(Fi,Li),FiInv)*detFi !< push forward to intermediate configuration temp_33 = matmul(FiInv,Li) - + do i = 1,3; do j = 1,3 dLi_dS(1:3,1:3,i,j) = matmul(matmul(Fi,dLi_dS(1:3,1:3,i,j)),FiInv)*detFi dLi_dFi(1:3,1:3,i,j) = dLi_dFi(1:3,1:3,i,j) + Li*FiInv(j,i) @@ -621,10 +621,10 @@ pure function constitutive_initialFi(ipc, ip, el) integer :: & phase, & homog, offset - + constitutive_initialFi = math_I3 phase = material_phaseAt(ipc,el) - + KinematicsLoop: do k = 1, phase_Nkinematics(phase) !< Warning: small initial strain assumption kinematicsType: select case (phase_kinematics(k,phase)) case (KINEMATICS_thermal_expansion_ID) kinematicsType @@ -640,7 +640,7 @@ end function constitutive_initialFi !-------------------------------------------------------------------------------------------------- !> @brief returns the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to -!> the elastic/intermediate deformation gradients depending on the selected elastic law +!> the elastic/intermediate deformation gradients depending on the selected elastic law !! (so far no case switch because only Hooke is implemented) !-------------------------------------------------------------------------------------------------- subroutine constitutive_SandItsTangents(S, dS_dFe, dS_dFi, Fe, Fi, ipc, ip, el) @@ -690,17 +690,17 @@ subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, & d !< counter in degradation loop integer :: & i, j - + ho = material_homogenizationAt(el) C = math_66toSym3333(constitutive_homogenizedC(ipc,ip,el)) - + DegradationLoop: do d = 1, phase_NstiffnessDegradations(material_phaseAt(ipc,el)) degradationType: select case(phase_stiffnessDegradation(d,material_phaseAt(ipc,el))) case (STIFFNESS_DEGRADATION_damage_ID) degradationType C = C * damage(ho)%p(damageMapping(ho)%p(ip,el))**2 end select degradationType enddo DegradationLoop - + E = 0.5_pReal*(matmul(transpose(Fe),Fe)-math_I3) !< Green-Lagrange strain in unloaded configuration S = math_mul3333xx33(C,matmul(matmul(transpose(Fi),E),Fi)) !< 2PK stress in lattice configuration in work conjugate with GL strain pulled back to lattice configuration @@ -715,7 +715,7 @@ end subroutine constitutive_hooke_SandItsTangents !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the rate of change of microstructure !-------------------------------------------------------------------------------------------------- -subroutine constitutive_collectDotState(S, FeArray, Fi, FpArray, subdt, ipc, ip, el) +subroutine constitutive_collectDotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el) integer, intent(in) :: & ipc, & !< component-ID of integration point @@ -724,7 +724,7 @@ subroutine constitutive_collectDotState(S, FeArray, Fi, FpArray, subdt, ipc, ip, real(pReal), intent(in) :: & subdt !< timestep real(pReal), intent(in), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: & - FeArray, & !< elastic deformation gradient + FArray, & !< elastic deformation gradient FpArray !< plastic deformation gradient real(pReal), intent(in), dimension(3,3) :: & Fi !< intermediate deformation gradient @@ -771,7 +771,7 @@ subroutine constitutive_collectDotState(S, FeArray, Fi, FpArray, subdt, ipc, ip, call plastic_disloucla_dotState (Mp,temperature(ho)%p(tme),instance,of) case (PLASTICITY_NONLOCAL_ID) plasticityType - call plastic_nonlocal_dotState (Mp,FeArray,FpArray,temperature(ho)%p(tme), & + call plastic_nonlocal_dotState (Mp,FArray,FpArray,temperature(ho)%p(tme), & subdt,ip,el) end select plasticityType @@ -858,32 +858,32 @@ subroutine constitutive_results do p=1,size(config_name_phase) group = trim('current/constituent')//'/'//trim(config_name_phase(p)) call results_closeGroup(results_addGroup(group)) - + group = trim(group)//'/plastic' - - call results_closeGroup(results_addGroup(group)) + + call results_closeGroup(results_addGroup(group)) select case(phase_plasticity(p)) - + case(PLASTICITY_ISOTROPIC_ID) - call plastic_isotropic_results(phase_plasticityInstance(p),group) - + call plastic_isotropic_results(phase_plasticityInstance(p),group) + case(PLASTICITY_PHENOPOWERLAW_ID) - call plastic_phenopowerlaw_results(phase_plasticityInstance(p),group) - + call plastic_phenopowerlaw_results(phase_plasticityInstance(p),group) + case(PLASTICITY_KINEHARDENING_ID) - call plastic_kinehardening_results(phase_plasticityInstance(p),group) + call plastic_kinehardening_results(phase_plasticityInstance(p),group) case(PLASTICITY_DISLOTWIN_ID) - call plastic_dislotwin_results(phase_plasticityInstance(p),group) - + call plastic_dislotwin_results(phase_plasticityInstance(p),group) + case(PLASTICITY_DISLOUCLA_ID) - call plastic_disloUCLA_results(phase_plasticityInstance(p),group) - + call plastic_disloUCLA_results(phase_plasticityInstance(p),group) + case(PLASTICITY_NONLOCAL_ID) - call plastic_nonlocal_results(phase_plasticityInstance(p),group) + call plastic_nonlocal_results(phase_plasticityInstance(p),group) end select - - enddo + + enddo end subroutine constitutive_results diff --git a/src/constitutive_plastic_nonlocal.f90 b/src/constitutive_plastic_nonlocal.f90 index 42b8a34f5..13d629e2b 100644 --- a/src/constitutive_plastic_nonlocal.f90 +++ b/src/constitutive_plastic_nonlocal.f90 @@ -11,7 +11,7 @@ submodule(constitutive) plastic_nonlocal IPvolume => geometry_plastic_nonlocal_IPvolume0, & IParea => geometry_plastic_nonlocal_IParea0, & IPareaNormal => geometry_plastic_nonlocal_IPareaNormal0 - + real(pReal), parameter :: & KB = 1.38e-23_pReal !< Physical parameter, Boltzmann constant in J/Kelvin @@ -44,11 +44,11 @@ submodule(constitutive) plastic_nonlocal integer, dimension(:), allocatable :: & totalNslip !< total number of active slip systems for each instance !END DEPRECATED - + real(pReal), dimension(:,:,:,:,:,:), allocatable :: & compatibility !< slip system compatibility between me and my neighbors - - enum, bind(c) + + enum, bind(c) enumerator :: & undefined_ID, & rho_sgl_mob_edg_pos_ID, & @@ -73,7 +73,7 @@ submodule(constitutive) plastic_nonlocal v_scr_neg_ID, & gamma_ID end enum - + type :: tParameters !< container type for internal constitutive parameters real(pReal) :: & atomicVolume, & !< atomic volume @@ -135,27 +135,27 @@ submodule(constitutive) plastic_nonlocal integer, dimension(:) ,allocatable :: & Nslip,& colinearSystem !< colinear system to the active slip system (only valid for fcc!) - + logical :: & shortRangeStressCorrection, & !< flag indicating the use of the short range stress correction by a excess density gradient term probabilisticMultiplication - + integer(kind(undefined_ID)), dimension(:), allocatable :: & outputID !< ID of each post result output - + end type tParameters - + type :: tNonlocalMicrostructure real(pReal), allocatable, dimension(:,:) :: & - tau_pass, & - tau_Back + tau_pass, & + tau_Back end type tNonlocalMicrostructure - + type :: tNonlocalState real(pReal), pointer, dimension(:,:) :: & rho, & ! < all dislocations rhoSgl, & - rhoSglMobile, & ! iRhoU + rhoSglMobile, & ! iRhoU rho_sgl_mob_edg_pos, & rho_sgl_mob_edg_neg, & rho_sgl_mob_scr_pos, & @@ -176,14 +176,15 @@ submodule(constitutive) plastic_nonlocal v_scr_pos, & v_scr_neg end type tNonlocalState - + type(tNonlocalState), allocatable, dimension(:) :: & deltaState, & dotState, & - state - + state, & + state0 + type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance) - + type(tNonlocalMicrostructure), dimension(:), allocatable :: microstructure contains @@ -210,8 +211,8 @@ module subroutine plastic_nonlocal_init extmsg = '', & structure character(len=pStringLen), dimension(:), allocatable :: outputs - integer :: NofMyPhase - + integer :: NofMyPhase + write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_NONLOCAL_label//' init -+>>>' write(6,'(/,a)') ' Reuber et al., Acta Materialia 71:333–348, 2014' @@ -226,6 +227,7 @@ module subroutine plastic_nonlocal_init allocate(param(maxNinstances)) allocate(state(maxNinstances)) + allocate(state0(maxNinstances)) allocate(dotState(maxNinstances)) allocate(deltaState(maxNinstances)) allocate(microstructure(maxNinstances)) @@ -238,13 +240,14 @@ module subroutine plastic_nonlocal_init associate(prm => param(phase_plasticityInstance(p)), & dot => dotState(phase_plasticityInstance(p)), & stt => state(phase_plasticityInstance(p)), & + st0 => state0(phase_plasticityInstance(p)), & del => deltaState(phase_plasticityInstance(p)), & dst => microstructure(phase_plasticityInstance(p)), & config => config_phase(p)) prm%aTolRho = config%getFloat('atol_rho', defaultVal=0.0_pReal) prm%aTolShear = config%getFloat('atol_shear', defaultVal=0.0_pReal) - + structure = config%getString('lattice_structure') ! This data is read in already in lattice @@ -293,20 +296,20 @@ module subroutine plastic_nonlocal_init prm%colinearSystem(s1) = s2 enddo enddo - + prm%rhoSglEdgePos0 = config%getFloats('rhosgledgepos0', requiredSize=size(prm%Nslip)) prm%rhoSglEdgeNeg0 = config%getFloats('rhosgledgeneg0', requiredSize=size(prm%Nslip)) prm%rhoSglScrewPos0 = config%getFloats('rhosglscrewpos0', requiredSize=size(prm%Nslip)) prm%rhoSglScrewNeg0 = config%getFloats('rhosglscrewneg0', requiredSize=size(prm%Nslip)) prm%rhoDipEdge0 = config%getFloats('rhodipedge0', requiredSize=size(prm%Nslip)) prm%rhoDipScrew0 = config%getFloats('rhodipscrew0', requiredSize=size(prm%Nslip)) - + prm%lambda0 = config%getFloats('lambda0', requiredSize=size(prm%Nslip)) prm%burgers = config%getFloats('burgers', requiredSize=size(prm%Nslip)) - + prm%lambda0 = math_expand(prm%lambda0,prm%Nslip) prm%burgers = math_expand(prm%burgers,prm%Nslip) - + prm%minDipoleHeight_edge = config%getFloats('minimumdipoleheightedge', requiredSize=size(prm%Nslip)) prm%minDipoleHeight_screw = config%getFloats('minimumdipoleheightscrew', requiredSize=size(prm%Nslip)) prm%minDipoleHeight_edge = math_expand(prm%minDipoleHeight_edge,prm%Nslip) @@ -314,7 +317,7 @@ module subroutine plastic_nonlocal_init allocate(prm%minDipoleHeight(prm%totalNslip,2)) prm%minDipoleHeight(:,1) = prm%minDipoleHeight_edge prm%minDipoleHeight(:,2) = prm%minDipoleHeight_screw - + prm%peierlsstress_edge = config%getFloats('peierlsstressedge', requiredSize=size(prm%Nslip)) prm%peierlsstress_screw = config%getFloats('peierlsstressscrew', requiredSize=size(prm%Nslip)) prm%peierlsstress_edge = math_expand(prm%peierlsstress_edge,prm%Nslip) @@ -326,7 +329,7 @@ module subroutine plastic_nonlocal_init prm%significantRho = config%getFloat('significantrho') prm%significantN = config%getFloat('significantn', 0.0_pReal) prm%CFLfactor = config%getFloat('cflfactor',defaultVal=2.0_pReal) - + prm%atomicVolume = config%getFloat('atomicvolume') prm%Dsd0 = config%getFloat('selfdiffusionprefactor') !,'dsd0' prm%selfDiffusionEnergy = config%getFloat('selfdiffusionenergy') !,'qsd' @@ -336,7 +339,7 @@ module subroutine plastic_nonlocal_init prm%solidSolutionEnergy = config%getFloat('solidsolutionenergy') prm%solidSolutionSize = config%getFloat('solidsolutionsize') prm%solidSolutionConcentration = config%getFloat('solidsolutionconcentration') - + prm%p = config%getFloat('p') prm%q = config%getFloat('q') prm%viscosity = config%getFloat('viscosity') @@ -345,62 +348,62 @@ module subroutine plastic_nonlocal_init ! ToDo: discuss logic prm%rhoSglScatter = config%getFloat('rhosglscatter') prm%rhoSglRandom = config%getFloat('rhosglrandom',0.0_pReal) - if (config%keyExists('/rhosglrandom/')) & + if (config%keyExists('/rhosglrandom/')) & prm%rhoSglRandomBinning = config%getFloat('rhosglrandombinning',0.0_pReal) !ToDo: useful default? ! if (rhoSglRandom(instance) < 0.0_pReal) & ! if (rhoSglRandomBinning(instance) <= 0.0_pReal) & - + prm%surfaceTransmissivity = config%getFloat('surfacetransmissivity',defaultVal=1.0_pReal) prm%grainboundaryTransmissivity = config%getFloat('grainboundarytransmissivity',defaultVal=-1.0_pReal) prm%fEdgeMultiplication = config%getFloat('edgemultiplication') prm%shortRangeStressCorrection = config%keyExists('/shortrangestresscorrection/') - + !-------------------------------------------------------------------------------------------------- ! sanity checks if (any(prm%burgers < 0.0_pReal)) extmsg = trim(extmsg)//' burgers' if (any(prm%lambda0 <= 0.0_pReal)) extmsg = trim(extmsg)//' lambda0' - + if (any(prm%rhoSglEdgePos0 < 0.0_pReal)) extmsg = trim(extmsg)//' rhoSglEdgePos0' if (any(prm%rhoSglEdgeNeg0 < 0.0_pReal)) extmsg = trim(extmsg)//' rhoSglEdgeNeg0' if (any(prm%rhoSglScrewPos0 < 0.0_pReal)) extmsg = trim(extmsg)//' rhoSglScrewPos0' if (any(prm%rhoSglScrewNeg0 < 0.0_pReal)) extmsg = trim(extmsg)//' rhoSglScrewNeg0' if (any(prm%rhoDipEdge0 < 0.0_pReal)) extmsg = trim(extmsg)//' rhoDipEdge0' if (any(prm%rhoDipScrew0 < 0.0_pReal)) extmsg = trim(extmsg)//' rhoDipScrew0' - + if (any(prm%peierlsstress < 0.0_pReal)) extmsg = trim(extmsg)//' peierlsstress' if (any(prm%minDipoleHeight < 0.0_pReal)) extmsg = trim(extmsg)//' minDipoleHeight' - + if (prm%viscosity <= 0.0_pReal) extmsg = trim(extmsg)//' viscosity' if (prm%selfDiffusionEnergy <= 0.0_pReal) extmsg = trim(extmsg)//' selfDiffusionEnergy' if (prm%fattack <= 0.0_pReal) extmsg = trim(extmsg)//' fattack' if (prm%doublekinkwidth <= 0.0_pReal) extmsg = trim(extmsg)//' doublekinkwidth' if (prm%Dsd0 < 0.0_pReal) extmsg = trim(extmsg)//' Dsd0' if (prm%atomicVolume <= 0.0_pReal) extmsg = trim(extmsg)//' atomicVolume' ! ToDo: in disloUCLA/dislotwin, the atomic volume is given as a factor - + if (prm%significantN < 0.0_pReal) extmsg = trim(extmsg)//' significantN' if (prm%significantrho < 0.0_pReal) extmsg = trim(extmsg)//' significantrho' if (prm%atolshear <= 0.0_pReal) extmsg = trim(extmsg)//' atolshear' if (prm%atolrho <= 0.0_pReal) extmsg = trim(extmsg)//' atolrho' if (prm%CFLfactor < 0.0_pReal) extmsg = trim(extmsg)//' CFLfactor' - - if (prm%p <= 0.0_pReal .or. prm%p > 1.0_pReal) extmsg = trim(extmsg)//' p' - if (prm%q < 1.0_pReal .or. prm%q > 2.0_pReal) extmsg = trim(extmsg)//' q' - + + if (prm%p <= 0.0_pReal .or. prm%p > 1.0_pReal) extmsg = trim(extmsg)//' p' + if (prm%q < 1.0_pReal .or. prm%q > 2.0_pReal) extmsg = trim(extmsg)//' q' + if (prm%linetensionEffect < 0.0_pReal .or. prm%linetensionEffect > 1.0_pReal) & extmsg = trim(extmsg)//' linetensionEffect' if (prm%edgeJogFactor < 0.0_pReal .or. prm%edgeJogFactor > 1.0_pReal) & extmsg = trim(extmsg)//' edgeJogFactor' - + if (prm%solidSolutionEnergy <= 0.0_pReal) extmsg = trim(extmsg)//' solidSolutionEnergy' if (prm%solidSolutionSize <= 0.0_pReal) extmsg = trim(extmsg)//' solidSolutionSize' if (prm%solidSolutionConcentration <= 0.0_pReal) extmsg = trim(extmsg)//' solidSolutionConcentration' - if (prm%grainboundaryTransmissivity > 1.0_pReal) extmsg = trim(extmsg)//' grainboundaryTransmissivity' + if (prm%grainboundaryTransmissivity > 1.0_pReal) extmsg = trim(extmsg)//' grainboundaryTransmissivity' if (prm%surfaceTransmissivity < 0.0_pReal .or. prm%surfaceTransmissivity > 1.0_pReal) & - extmsg = trim(extmsg)//' surfaceTransmissivity' - + extmsg = trim(extmsg)//' surfaceTransmissivity' + if (prm%fEdgeMultiplication < 0.0_pReal .or. prm%fEdgeMultiplication > 1.0_pReal) & - extmsg = trim(extmsg)//' fEdgeMultiplication' + extmsg = trim(extmsg)//' fEdgeMultiplication' endif slipActive @@ -470,39 +473,40 @@ module subroutine plastic_nonlocal_init 'gamma ' ]) * prm%totalNslip !< "basic" microstructural state variables that are independent from other state variables sizeDependentState = size([ 'rhoForest ']) * prm%totalNslip !< microstructural state variables that depend on other state variables sizeState = sizeDotState + sizeDependentState & - + size([ 'velocityEdgePos ','velocityEdgeNeg ', & + + size([ 'velocityEdgePos ','velocityEdgeNeg ', & 'velocityScrewPos ','velocityScrewNeg ', & 'maxDipoleHeightEdge ','maxDipoleHeightScrew' ]) * prm%totalNslip !< other dependent state variables that are not updated by microstructure sizeDeltaState = sizeDotState - + call material_allocatePlasticState(p,NofMyPhase,sizeState,sizeDotState,sizeDeltaState) plasticState(p)%nonlocal = .true. plasticState(p)%offsetDeltaState = 0 ! ToDo: state structure does not follow convention - + totalNslip(phase_plasticityInstance(p)) = prm%totalNslip - + + st0%rho => plasticState(p)%state0 (0*prm%totalNslip+1:10*prm%totalNslip,:) stt%rho => plasticState(p)%state (0*prm%totalNslip+1:10*prm%totalNslip,:) dot%rho => plasticState(p)%dotState (0*prm%totalNslip+1:10*prm%totalNslip,:) del%rho => plasticState(p)%deltaState (0*prm%totalNslip+1:10*prm%totalNslip,:) plasticState(p)%aTolState(1:10*prm%totalNslip) = prm%aTolRho - + stt%rhoSgl => plasticState(p)%state (0*prm%totalNslip+1: 8*prm%totalNslip,:) dot%rhoSgl => plasticState(p)%dotState (0*prm%totalNslip+1: 8*prm%totalNslip,:) del%rhoSgl => plasticState(p)%deltaState (0*prm%totalNslip+1: 8*prm%totalNslip,:) - + stt%rhoSglMobile => plasticState(p)%state (0*prm%totalNslip+1: 4*prm%totalNslip,:) dot%rhoSglMobile => plasticState(p)%dotState (0*prm%totalNslip+1: 4*prm%totalNslip,:) del%rhoSglMobile => plasticState(p)%deltaState (0*prm%totalNslip+1: 4*prm%totalNslip,:) - + stt%rho_sgl_mob_edg_pos => plasticState(p)%state (0*prm%totalNslip+1: 1*prm%totalNslip,:) dot%rho_sgl_mob_edg_pos => plasticState(p)%dotState (0*prm%totalNslip+1: 1*prm%totalNslip,:) del%rho_sgl_mob_edg_pos => plasticState(p)%deltaState (0*prm%totalNslip+1: 1*prm%totalNslip,:) - + stt%rho_sgl_mob_edg_neg => plasticState(p)%state (1*prm%totalNslip+1: 2*prm%totalNslip,:) dot%rho_sgl_mob_edg_neg => plasticState(p)%dotState (1*prm%totalNslip+1: 2*prm%totalNslip,:) del%rho_sgl_mob_edg_neg => plasticState(p)%deltaState (1*prm%totalNslip+1: 2*prm%totalNslip,:) - + stt%rho_sgl_mob_scr_pos => plasticState(p)%state (2*prm%totalNslip+1: 3*prm%totalNslip,:) dot%rho_sgl_mob_scr_pos => plasticState(p)%dotState (2*prm%totalNslip+1: 3*prm%totalNslip,:) del%rho_sgl_mob_scr_pos => plasticState(p)%deltaState (2*prm%totalNslip+1: 3*prm%totalNslip,:) @@ -510,46 +514,46 @@ module subroutine plastic_nonlocal_init stt%rho_sgl_mob_scr_neg => plasticState(p)%state (3*prm%totalNslip+1: 4*prm%totalNslip,:) dot%rho_sgl_mob_scr_neg => plasticState(p)%dotState (3*prm%totalNslip+1: 4*prm%totalNslip,:) del%rho_sgl_mob_scr_neg => plasticState(p)%deltaState (3*prm%totalNslip+1: 4*prm%totalNslip,:) - + stt%rhoSglImmobile => plasticState(p)%state (4*prm%totalNslip+1: 8*prm%totalNslip,:) dot%rhoSglImmobile => plasticState(p)%dotState (4*prm%totalNslip+1: 8*prm%totalNslip,:) del%rhoSglImmobile => plasticState(p)%deltaState (4*prm%totalNslip+1: 8*prm%totalNslip,:) - + stt%rho_sgl_imm_edg_pos => plasticState(p)%state (4*prm%totalNslip+1: 5*prm%totalNslip,:) dot%rho_sgl_imm_edg_pos => plasticState(p)%dotState (4*prm%totalNslip+1: 5*prm%totalNslip,:) del%rho_sgl_imm_edg_pos => plasticState(p)%deltaState (4*prm%totalNslip+1: 5*prm%totalNslip,:) - + stt%rho_sgl_imm_edg_neg => plasticState(p)%state (5*prm%totalNslip+1: 6*prm%totalNslip,:) dot%rho_sgl_imm_edg_neg => plasticState(p)%dotState (5*prm%totalNslip+1: 6*prm%totalNslip,:) del%rho_sgl_imm_edg_neg => plasticState(p)%deltaState (5*prm%totalNslip+1: 6*prm%totalNslip,:) - + stt%rho_sgl_imm_scr_pos => plasticState(p)%state (6*prm%totalNslip+1: 7*prm%totalNslip,:) dot%rho_sgl_imm_scr_pos => plasticState(p)%dotState (6*prm%totalNslip+1: 7*prm%totalNslip,:) del%rho_sgl_imm_scr_pos => plasticState(p)%deltaState (6*prm%totalNslip+1: 7*prm%totalNslip,:) - + stt%rho_sgl_imm_scr_neg => plasticState(p)%state (7*prm%totalNslip+1: 8*prm%totalNslip,:) dot%rho_sgl_imm_scr_neg => plasticState(p)%dotState (7*prm%totalNslip+1: 8*prm%totalNslip,:) del%rho_sgl_imm_scr_neg => plasticState(p)%deltaState (7*prm%totalNslip+1: 8*prm%totalNslip,:) - + stt%rhoDip => plasticState(p)%state (8*prm%totalNslip+1:10*prm%totalNslip,:) dot%rhoDip => plasticState(p)%dotState (8*prm%totalNslip+1:10*prm%totalNslip,:) del%rhoDip => plasticState(p)%deltaState (8*prm%totalNslip+1:10*prm%totalNslip,:) - + stt%rho_dip_edg => plasticState(p)%state (8*prm%totalNslip+1: 9*prm%totalNslip,:) dot%rho_dip_edg => plasticState(p)%dotState (8*prm%totalNslip+1: 9*prm%totalNslip,:) del%rho_dip_edg => plasticState(p)%deltaState (8*prm%totalNslip+1: 9*prm%totalNslip,:) - + stt%rho_dip_scr => plasticState(p)%state (9*prm%totalNslip+1:10*prm%totalNslip,:) dot%rho_dip_scr => plasticState(p)%dotState (9*prm%totalNslip+1:10*prm%totalNslip,:) del%rho_dip_scr => plasticState(p)%deltaState (9*prm%totalNslip+1:10*prm%totalNslip,:) - + stt%gamma => plasticState(p)%state (10*prm%totalNslip + 1:11*prm%totalNslip ,1:NofMyPhase) dot%gamma => plasticState(p)%dotState (10*prm%totalNslip + 1:11*prm%totalNslip ,1:NofMyPhase) del%gamma => plasticState(p)%deltaState (10*prm%totalNslip + 1:11*prm%totalNslip ,1:NofMyPhase) plasticState(p)%aTolState(10*prm%totalNslip + 1:11*prm%totalNslip ) = prm%aTolShear plasticState(p)%slipRate => plasticState(p)%dotState (10*prm%totalNslip + 1:11*prm%totalNslip ,1:NofMyPhase) plasticState(p)%accumulatedSlip => plasticState(p)%state(10*prm%totalNslip + 1:11*prm%totalNslip ,1:NofMyPhase) - + stt%rho_forest => plasticState(p)%state (11*prm%totalNslip + 1:12*prm%totalNslip ,1:NofMyPhase) stt%v => plasticState(p)%state (12*prm%totalNslip + 1:16*prm%totalNslip ,1:NofMyPhase) stt%v_edg_pos => plasticState(p)%state (12*prm%totalNslip + 1:13*prm%totalNslip ,1:NofMyPhase) @@ -565,10 +569,10 @@ module subroutine plastic_nonlocal_init plasticState(p)%state0 = plasticState(p)%state enddo - + allocate(compatibility(2,maxval(totalNslip),maxval(totalNslip),nIPneighbors,& discretization_nIP,discretization_nElem), source=0.0_pReal) - + ! BEGIN DEPRECATED---------------------------------------------------------------------------------- allocate(iRhoU(maxval(totalNslip),4,maxNinstances), source=0) allocate(iRhoB(maxval(totalNslip),4,maxNinstances), source=0) @@ -601,7 +605,7 @@ module subroutine plastic_nonlocal_init iRhoD(s,c,phase_plasticityInstance(p)) = l enddo enddo - + l = l + param(phase_plasticityInstance(p))%totalNslip ! shear(rates) l = l + param(phase_plasticityInstance(p))%totalNslip ! rho_forest @@ -619,20 +623,20 @@ module subroutine plastic_nonlocal_init enddo if (iD(param(phase_plasticityInstance(p))%totalNslip,2,phase_plasticityInstance(p)) /= plasticState(p)%sizeState) & call IO_error(0, ext_msg = 'state indices not properly set ('//PLASTICITY_NONLOCAL_label//')') - - + + endif myPhase2 - + enddo initializeInstances ! END DEPRECATED------------------------------------------------------------------------------------ - + contains !-------------------------------------------------------------------------------------------------- !> @brief populates the initial dislocation density !-------------------------------------------------------------------------------------------------- subroutine stateInit(phase,NofMyPhase) - + integer,intent(in) ::& phase, & NofMyPhase @@ -647,7 +651,7 @@ module subroutine plastic_nonlocal_init phasemember real(pReal), dimension(2) :: & noise, & - rnd + rnd real(pReal) :: & meanDensity, & totalVolume, & @@ -655,14 +659,14 @@ module subroutine plastic_nonlocal_init minimumIpVolume real(pReal), dimension(NofMyPhase) :: & volume - - + + instance = phase_plasticityInstance(phase) associate(prm => param(instance), stt => state(instance)) - - ! randomly distribute dislocation segments on random slip system and of random type in the volume + + ! randomly distribute dislocation segments on random slip system and of random type in the volume if (prm%rhoSglRandom > 0.0_pReal) then - + ! get the total volume of the instance do e = 1,discretization_nElem do i = 1,discretization_nIP @@ -672,7 +676,7 @@ module subroutine plastic_nonlocal_init totalVolume = sum(volume) minimumIPVolume = minval(volume) densityBinning = prm%rhoSglRandomBinning / minimumIpVolume ** (2.0_pReal / 3.0_pReal) - + ! subsequently fill random ips with dislocation segments until we reach the desired overall density meanDensity = 0.0_pReal do while(meanDensity < prm%rhoSglRandom) @@ -701,9 +705,9 @@ module subroutine plastic_nonlocal_init enddo enddo endif - + end associate - + end subroutine stateInit end subroutine plastic_nonlocal_init @@ -712,15 +716,15 @@ end subroutine plastic_nonlocal_init !-------------------------------------------------------------------------------------------------- !> @brief calculates quantities characterizing the microstructure !-------------------------------------------------------------------------------------------------- -module subroutine plastic_nonlocal_dependentState(Fe, Fp, ip, el) - +module subroutine plastic_nonlocal_dependentState(F, Fp, ip, el) + integer, intent(in) :: & ip, & el real(pReal), dimension(3,3), intent(in) :: & - Fe, & + F, & Fp - + integer :: & ph, & !< phase of, & !< offset @@ -761,11 +765,12 @@ module subroutine plastic_nonlocal_dependentState(Fe, Fp, ip, el) rho_scr_delta real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el))),10) :: & rho, & - rho_neighbor + rho0, & + rho_neighbor0 real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el))), & totalNslip(phase_plasticityInstance(material_phaseAt(1,el)))) :: & myInteractionMatrix ! corrected slip interaction matrix - + real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el))),nIPneighbors) :: & rho_edg_delta_neighbor, & rho_scr_delta_neighbor @@ -774,30 +779,30 @@ module subroutine plastic_nonlocal_dependentState(Fe, Fp, ip, el) neighbor_rhoTotal ! total density at neighboring material point real(pReal), dimension(3,totalNslip(phase_plasticityInstance(material_phaseAt(1,el))),2) :: & m ! direction of dislocation motion - + ph = material_phaseAt(1,el) of = material_phasememberAt(1,ip,el) instance = phase_plasticityInstance(ph) associate(prm => param(instance),dst => microstructure(instance), stt => state(instance)) - + ns = prm%totalNslip - + rho = getRho(instance,of,ip,el) - + stt%rho_forest(:,of) = matmul(prm%forestProjection_Edge, sum(abs(rho(:,edg)),2)) & - + matmul(prm%forestProjection_Screw,sum(abs(rho(:,scr)),2)) - - - ! coefficients are corrected for the line tension effect + + matmul(prm%forestProjection_Screw,sum(abs(rho(:,scr)),2)) + + + ! coefficients are corrected for the line tension effect ! (see Kubin,Devincre,Hoc; 2008; Modeling dislocation storage rates and mean free paths in face-centered cubic crystals) - if (lattice_structure(ph) == LATTICE_bcc_ID .or. lattice_structure(ph) == LATTICE_fcc_ID) then ! only fcc and bcc + if (lattice_structure(ph) == LATTICE_bcc_ID .or. lattice_structure(ph) == LATTICE_fcc_ID) then do s = 1,ns correction = ( 1.0_pReal - prm%linetensionEffect & + prm%linetensionEffect & * log(0.35_pReal * prm%burgers(s) * sqrt(max(stt%rho_forest(s,of),prm%significantRho))) & / log(0.35_pReal * prm%burgers(s) * 1e6_pReal)) ** 2.0_pReal - myInteractionMatrix(1:ns,s) = correction * prm%interactionSlipSlip(1:ns,s) + myInteractionMatrix(1:ns,s) = correction * prm%interactionSlipSlip(1:ns,s) enddo else myInteractionMatrix = prm%interactionSlipSlip @@ -815,21 +820,21 @@ module subroutine plastic_nonlocal_dependentState(Fe, Fp, ip, el) ! ToDo: MD: this is most likely only correct for F_i = I !################################################################################################# - + rho0 = getRho0(instance,of,ip,el) if (.not. phase_localPlasticity(ph) .and. prm%shortRangeStressCorrection) then - invFe = math_inv33(Fe) invFp = math_inv33(Fp) - - rho_edg_delta = rho(:,mob_edg_pos) - rho(:,mob_edg_neg) - rho_scr_delta = rho(:,mob_scr_pos) - rho(:,mob_scr_neg) - + invFe = matmul(Fp,math_inv33(F)) + + rho_edg_delta = rho0(:,mob_edg_pos) - rho0(:,mob_edg_neg) + rho_scr_delta = rho0(:,mob_scr_pos) - rho0(:,mob_scr_neg) + rhoExcess(1,1:ns) = rho_edg_delta rhoExcess(2,1:ns) = rho_scr_delta - + FVsize = IPvolume(ip,el) ** (1.0_pReal/3.0_pReal) - + !* loop through my neighborhood and get the connection vectors (in lattice frame) and the excess densities - + nRealNeighbors = 0.0_pReal neighbor_rhoTotal = 0.0_pReal do n = 1,nIPneighbors @@ -839,16 +844,16 @@ module subroutine plastic_nonlocal_dependentState(Fe, Fp, ip, el) if (neighbor_el > 0 .and. neighbor_ip > 0) then neighbor_instance = phase_plasticityInstance(material_phaseAt(1,neighbor_el)) if (neighbor_instance == instance) then - + nRealNeighbors = nRealNeighbors + 1.0_pReal - rho_neighbor = getRho(instance,no,neighbor_ip,neighbor_el) - - rho_edg_delta_neighbor(:,n) = rho_neighbor(:,mob_edg_pos) - rho_neighbor(:,mob_edg_neg) - rho_scr_delta_neighbor(:,n) = rho_neighbor(:,mob_scr_pos) - rho_neighbor(:,mob_scr_neg) - - neighbor_rhoTotal(1,:,n) = sum(abs(rho_neighbor(:,edg)),2) - neighbor_rhoTotal(2,:,n) = sum(abs(rho_neighbor(:,scr)),2) - + rho_neighbor0 = getRho0(instance,no,neighbor_ip,neighbor_el) + + rho_edg_delta_neighbor(:,n) = rho_neighbor0(:,mob_edg_pos) - rho_neighbor0(:,mob_edg_neg) + rho_scr_delta_neighbor(:,n) = rho_neighbor0(:,mob_scr_pos) - rho_neighbor0(:,mob_scr_neg) + + neighbor_rhoTotal(1,:,n) = sum(abs(rho_neighbor0(:,edg)),2) + neighbor_rhoTotal(2,:,n) = sum(abs(rho_neighbor0(:,scr)),2) + connection_latticeConf(1:3,n) = matmul(invFe, discretization_IPcoords(1:3,neighbor_el+neighbor_ip-1) & - discretization_IPcoords(1:3,el+neighbor_ip-1)) normal_latticeConf = matmul(transpose(invFp), IPareaNormal(1:3,n,ip,el)) @@ -867,18 +872,18 @@ module subroutine plastic_nonlocal_dependentState(Fe, Fp, ip, el) rho_scr_delta_neighbor(:,n) = rho_scr_delta endif enddo - + neighbor_rhoExcess(1,:,:) = rho_edg_delta_neighbor neighbor_rhoExcess(2,:,:) = rho_scr_delta_neighbor - + !* loop through the slip systems and calculate the dislocation gradient by !* 1. interpolation of the excess density in the neighorhood !* 2. interpolation of the dead dislocation density in the central volume m(1:3,1:ns,1) = prm%slip_direction m(1:3,1:ns,2) = -prm%slip_transverse - + do s = 1,ns - + ! gradient from interpolation of neighboring excess density ... do c = 1,2 do dir = 1,3 @@ -891,27 +896,27 @@ module subroutine plastic_nonlocal_dependentState(Fe, Fp, ip, el) enddo invConnections = math_inv33(connections) if (all(dEq0(invConnections))) call IO_error(-1,ext_msg='back stress calculation: inversion error') - + rhoExcessGradient(c) = math_inner(m(1:3,s,c), matmul(invConnections,rhoExcessDifferences)) enddo - + ! ... plus gradient from deads ... rhoExcessGradient(1) = rhoExcessGradient(1) + sum(rho(s,imm_edg)) / FVsize rhoExcessGradient(2) = rhoExcessGradient(2) + sum(rho(s,imm_scr)) / FVsize - + ! ... normalized with the total density ... rhoTotal(1) = (sum(abs(rho(s,edg))) + sum(neighbor_rhoTotal(1,s,:))) / (1.0_pReal + nRealNeighbors) rhoTotal(2) = (sum(abs(rho(s,scr))) + sum(neighbor_rhoTotal(2,s,:))) / (1.0_pReal + nRealNeighbors) - + rhoExcessGradient_over_rho = 0.0_pReal where(rhoTotal > 0.0_pReal) & rhoExcessGradient_over_rho = rhoExcessGradient / rhoTotal - + ! ... gives the local stress correction when multiplied with a factor dst%tau_back(s,of) = - prm%mu * prm%burgers(s) / (2.0_pReal * pi) & * (rhoExcessGradient_over_rho(1) / (1.0_pReal - prm%nu) & + rhoExcessGradient_over_rho(2)) - + enddo endif @@ -932,7 +937,7 @@ end subroutine plastic_nonlocal_dependentState !-------------------------------------------------------------------------------------------------- -!> @brief calculates kinetics +!> @brief calculates kinetics !-------------------------------------------------------------------------------------------------- subroutine plastic_nonlocal_kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, & tauThreshold, c, Temperature, instance, of) @@ -945,17 +950,17 @@ subroutine plastic_nonlocal_kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, & tau, & !< resolved external shear stress (without non Schmid effects) tauNS, & !< resolved external shear stress (including non Schmid effects) tauThreshold !< threshold shear stress - + real(pReal), dimension(param(instance)%totalNslip), intent(out) :: & v, & !< velocity dv_dtau, & !< velocity derivative with respect to resolved shear stress (without non Schmid contributions) dv_dtauNS !< velocity derivative with respect to resolved shear stress (including non Schmid contributions) - + integer :: & ns, & !< short notation for the total number of active slip systems s !< index of my current slip system - real(pReal) :: & - tauRel_P, & + real(pReal) :: & + tauRel_P, & tauRel_S, & tauEff, & !< effective shear stress tPeierls, & !< waiting time in front of a peierls barriers @@ -976,22 +981,22 @@ subroutine plastic_nonlocal_kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, & criticalStress_P, & !< maximum obstacle strength criticalStress_S, & !< maximum obstacle strength mobility !< dislocation mobility - + associate(prm => param(instance)) ns = prm%totalNslip v = 0.0_pReal dv_dtau = 0.0_pReal dv_dtauNS = 0.0_pReal - - + + if (Temperature > 0.0_pReal) then do s = 1,ns if (abs(tau(s)) > tauThreshold(s)) then - + !* Peierls contribution !* Effective stress includes non Schmid constributions !* The derivative only gives absolute values; the correct sign is taken care of in the formula for the derivative of the velocity - + tauEff = max(0.0_pReal, abs(tauNS(s)) - tauThreshold(s)) ! ensure that the effective stress is positive meanfreepath_P = prm%burgers(s) jumpWidth_P = prm%burgers(s) @@ -1006,15 +1011,15 @@ subroutine plastic_nonlocal_kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, & if (tauEff < criticalStress_P) then dtPeierls_dtau = tPeierls * prm%p * prm%q * activationVolume_P / (KB * Temperature) & * (1.0_pReal - tauRel_P**prm%p)**(prm%q-1.0_pReal) & - * tauRel_P**(prm%p-1.0_pReal) + * tauRel_P**(prm%p-1.0_pReal) else dtPeierls_dtau = 0.0_pReal endif - - + + !* Contribution from solid solution strengthening !* The derivative only gives absolute values; the correct sign is taken care of in the formula for the derivative of the velocity - + tauEff = abs(tau(s)) - tauThreshold(s) meanfreepath_S = prm%burgers(s) / sqrt(prm%solidSolutionConcentration) jumpWidth_S = prm%solidSolutionSize * prm%burgers(s) @@ -1030,32 +1035,32 @@ subroutine plastic_nonlocal_kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, & dtSolidSolution_dtau = tSolidSolution * prm%p * prm%q & * activationVolume_S / (KB * Temperature) & * (1.0_pReal - tauRel_S**prm%p)**(prm%q-1.0_pReal) & - * tauRel_S**(prm%p-1.0_pReal) + * tauRel_S**(prm%p-1.0_pReal) else dtSolidSolution_dtau = 0.0_pReal endif - - + + !* viscous glide velocity - + tauEff = abs(tau(s)) - tauThreshold(s) mobility = prm%burgers(s) / prm%viscosity vViscous = mobility * tauEff - - - !* Mean velocity results from waiting time at peierls barriers and solid solution obstacles with respective meanfreepath of - !* free flight at glide velocity in between. + + + !* Mean velocity results from waiting time at peierls barriers and solid solution obstacles with respective meanfreepath of + !* free flight at glide velocity in between. !* adopt sign from resolved stress - + v(s) = sign(1.0_pReal,tau(s)) & / (tPeierls / meanfreepath_P + tSolidSolution / meanfreepath_S + 1.0_pReal / vViscous) dv_dtau(s) = v(s) * v(s) * (dtSolidSolution_dtau / meanfreepath_S & + mobility / (vViscous * vViscous)) - dv_dtauNS(s) = v(s) * v(s) * dtPeierls_dtau / meanfreepath_P + dv_dtauNS(s) = v(s) * v(s) * dtPeierls_dtau / meanfreepath_P endif enddo endif - + #ifdef DEBUGTODO write(6,'(a,/,12x,12(f12.5,1x))') '<< CONST >> tauThreshold / MPa', tauThreshold * 1e-6_pReal @@ -1089,8 +1094,8 @@ module subroutine plastic_nonlocal_LpAndItsTangent(Lp, dLp_dMp, & Lp !< plastic velocity gradient real(pReal), dimension(3,3,3,3), intent(out) :: & dLp_dMp !< derivative of Lp with respect to Tstar (9x9 matrix) - - + + integer :: & instance, & !< current instance of this plasticity ns, & !< short notation for the total number of active slip systems @@ -1114,20 +1119,20 @@ module subroutine plastic_nonlocal_LpAndItsTangent(Lp, dLp_dMp, & real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el)))) :: & tau, & !< resolved shear stress including backstress terms gdotTotal !< shear rate - + !*** shortcut for mapping ph = material_phaseAt(1,el) of = material_phasememberAt(1,ip,el) - + instance = phase_plasticityInstance(ph) associate(prm => param(instance),dst=>microstructure(instance),stt=>state(instance)) ns = prm%totalNslip - - !*** shortcut to state variables + + !*** shortcut to state variables rho = getRho(instance,of,ip,el) rhoSgl = rho(:,sgl) - - + + !*** get resolved shear stress !*** for screws possible non-schmid contributions are also taken into account do s = 1,ns @@ -1145,18 +1150,18 @@ module subroutine plastic_nonlocal_LpAndItsTangent(Lp, dLp_dMp, & forall (t = 1:4) & tauNS(1:ns,t) = tauNS(1:ns,t) + dst%tau_back(:,of) tau = tau + dst%tau_back(:,of) - - + + !*** get dislocation velocity and its tangent and store the velocity in the state array - - ! edges + + ! edges call plastic_nonlocal_kinetics(v(1:ns,1), dv_dtau(1:ns,1), dv_dtauNS(1:ns,1), & tau(1:ns), tauNS(1:ns,1), dst%tau_pass(1:ns,of), & 1, Temperature, instance, of) v(1:ns,2) = v(1:ns,1) dv_dtau(1:ns,2) = dv_dtau(1:ns,1) dv_dtauNS(1:ns,2) = dv_dtauNS(1:ns,1) - + !screws if (size(prm%nonSchmidCoeff) == 0) then forall(t = 3:4) @@ -1173,17 +1178,17 @@ module subroutine plastic_nonlocal_LpAndItsTangent(Lp, dLp_dMp, & endif stt%v(:,of) = pack(v,.true.) - + !*** Bauschinger effect forall (s = 1:ns, t = 5:8, rhoSgl(s,t) * v(s,t-4) < 0.0_pReal) & rhoSgl(s,t-4) = rhoSgl(s,t-4) + abs(rhoSgl(s,t)) - - + + gdotTotal = sum(rhoSgl(1:ns,1:4) * v, 2) * prm%burgers(1:ns) - + Lp = 0.0_pReal dLp_dMp = 0.0_pReal - + do s = 1,ns Lp = Lp + gdotTotal(s) * prm%Schmid(1:3,1:3,s) forall (i=1:3,j=1:3,k=1:3,l=1:3) & @@ -1191,13 +1196,13 @@ module subroutine plastic_nonlocal_LpAndItsTangent(Lp, dLp_dMp, & + prm%Schmid(i,j,s) * prm%Schmid(k,l,s) & * sum(rhoSgl(s,1:4) * dv_dtau(s,1:4)) * prm%burgers(s) & + prm%Schmid(i,j,s) & - * ( prm%nonSchmid_pos(k,l,s) * rhoSgl(s,3) * dv_dtauNS(s,3) & + * ( prm%nonSchmid_pos(k,l,s) * rhoSgl(s,3) * dv_dtauNS(s,3) & - prm%nonSchmid_neg(k,l,s) * rhoSgl(s,4) * dv_dtauNS(s,4)) * prm%burgers(s) enddo - - + + end associate - + end subroutine plastic_nonlocal_LpAndItsTangent @@ -1205,7 +1210,7 @@ end subroutine plastic_nonlocal_LpAndItsTangent !> @brief (instantaneous) incremental change of microstructure !-------------------------------------------------------------------------------------------------- module subroutine plastic_nonlocal_deltaState(Mp,ip,el) - + integer, intent(in) :: & ip, & el @@ -1234,23 +1239,23 @@ module subroutine plastic_nonlocal_deltaState(Mp,ip,el) dUpper, & ! current maximum stable dipole distance for edges and screws dUpperOld, & ! old maximum stable dipole distance for edges and screws deltaDUpper ! change in maximum stable dipole distance for edges and screws - + ph = material_phaseAt(1,el) of = material_phasememberAt(1,ip,el) instance = phase_plasticityInstance(ph) associate(prm => param(instance),dst => microstructure(instance),del => deltaState(instance)) ns = totalNslip(instance) - - !*** shortcut to state variables + + !*** shortcut to state variables forall (s = 1:ns, t = 1:4) & v(s,t) = plasticState(ph)%state(iV(s,t,instance),of) forall (s = 1:ns, c = 1:2) & dUpperOld(s,c) = plasticState(ph)%state(iD(s,c,instance),of) - + rho = getRho(instance,of,ip,el) rhoDip = rho(:,dip) - + !**************************************************************************** !*** dislocation remobilization (bauschinger effect) where(rho(:,imm) * v < 0.0_pReal) @@ -1261,48 +1266,48 @@ module subroutine plastic_nonlocal_deltaState(Mp,ip,el) elsewhere deltaRhoRemobilization(:,mob) = 0.0_pReal deltaRhoRemobilization(:,imm) = 0.0_pReal - endwhere + endwhere deltaRhoRemobilization(:,dip) = 0.0_pReal - + !**************************************************************************** !*** calculate dipole formation and dissociation by stress change - + !*** calculate limits for stable dipole height do s = 1,prm%totalNslip tau(s) = math_mul33xx33(Mp, prm%Schmid(1:3,1:3,s)) +dst%tau_back(s,of) if (abs(tau(s)) < 1.0e-15_pReal) tau(s) = 1.0e-15_pReal enddo - + dUpper(1:ns,1) = prm%mu * prm%burgers/(8.0_pReal * PI * (1.0_pReal - prm%nu) * abs(tau)) dUpper(1:ns,2) = prm%mu * prm%burgers/(4.0_pReal * PI * abs(tau)) - + where(dNeq0(sqrt(sum(abs(rho(:,edg)),2)))) & dUpper(1:ns,1) = min(1.0_pReal/sqrt(sum(abs(rho(:,edg)),2)),dUpper(1:ns,1)) - + where(dNeq0(sqrt(sum(abs(rho(:,scr)),2)))) & dUpper(1:ns,2) = min(1.0_pReal/sqrt(sum(abs(rho(:,scr)),2)),dUpper(1:ns,2)) - + dUpper = max(dUpper,prm%minDipoleHeight) deltaDUpper = dUpper - dUpperOld - - + + !*** dissociation by stress increase deltaRhoDipole2SingleStress = 0.0_pReal forall (c=1:2, s=1:ns, deltaDUpper(s,c) < 0.0_pReal .and. & dNeq0(dUpperOld(s,c) - prm%minDipoleHeight(s,c))) & deltaRhoDipole2SingleStress(s,8+c) = rhoDip(s,c) * deltaDUpper(s,c) & / (dUpperOld(s,c) - prm%minDipoleHeight(s,c)) - + forall (t=1:4) & deltaRhoDipole2SingleStress(1:ns,t) = -0.5_pReal & * deltaRhoDipole2SingleStress(1:ns,(t-1)/2+9) - + forall (s = 1:ns, c = 1:2) & - plasticState(ph)%state(iD(s,c,instance),of) = dUpper(s,c) - + plasticState(ph)%state(iD(s,c,instance),of) = dUpper(s,c) + plasticState(ph)%deltaState(:,of) = 0.0_pReal del%rho(:,of) = reshape(deltaRhoRemobilization + deltaRhoDipole2SingleStress, [10*ns]) - + #ifdef DEBUG if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0 & .and. ((debug_e == el .and. debug_i == ip)& @@ -1320,9 +1325,9 @@ end subroutine plastic_nonlocal_deltaState !--------------------------------------------------------------------------------------------------- !> @brief calculates the rate of change of microstructure !--------------------------------------------------------------------------------------------------- -module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & +module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature, & timestep,ip,el) - + integer, intent(in) :: & ip, & !< current integration point el !< current element number @@ -1332,11 +1337,11 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & real(pReal), dimension(3,3), intent(in) ::& Mp !< MandelStress real(pReal), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem), intent(in) :: & - Fe, & !< elastic deformation gradient + F, & !< elastic deformation gradient Fp !< plastic deformation gradient - + integer :: & - ph, & + ph, & instance, & !< current instance of this plasticity neighbor_instance, & !< instance of my neighbor's plasticity ns, & !< short notation for the total number of active slip systems @@ -1358,6 +1363,7 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & s !< index of my current slip system real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el))),10) :: & rho, & + rho0, & !< dislocation density at beginning of time step rhoDot, & !< density evolution rhoDotMultiplication, & !< density evolution by multiplication rhoDotFlux, & !< density evolution by flux @@ -1366,17 +1372,17 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & rhoDotThermalAnnihilation !< density evolution by thermal annihilation real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el))),8) :: & rhoSgl, & !< current single dislocation densities (positive/negative screw and edge without dipoles) - neighbor_rhoSgl, & !< current single dislocation densities of neighboring ip (positive/negative screw and edge without dipoles) - my_rhoSgl !< single dislocation densities of central ip (positive/negative screw and edge without dipoles) + neighbor_rhoSgl0, & !< current single dislocation densities of neighboring ip (positive/negative screw and edge without dipoles) + my_rhoSgl0 !< single dislocation densities of central ip (positive/negative screw and edge without dipoles) real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el))),4) :: & v, & !< current dislocation glide velocity - my_v, & !< dislocation glide velocity of central ip - neighbor_v, & !< dislocation glide velocity of enighboring ip + v0, & + neighbor_v0, & !< dislocation glide velocity of enighboring ip gdot !< shear rates real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el)))) :: & tau, & !< current resolved shear stress vClimb !< climb velocity of edge dipoles - + real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el))),2) :: & rhoDip, & !< current dipole dislocation densities (screw and edge dipoles) dLower, & !< minimum stable dipole distance for edges and screws @@ -1399,19 +1405,19 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & transmissivity, & !< overall transmissivity of dislocation flux to neighboring material point lineLength, & !< dislocation line length leaving the current interface selfDiffusion !< self diffusion - + logical :: & considerEnteringFlux, & considerLeavingFlux - + p = material_phaseAt(1,el) o = material_phasememberAt(1,ip,el) - + if (timestep <= 0.0_pReal) then plasticState(p)%dotState = 0.0_pReal return endif - + ph = material_phaseAt(1,el) instance = phase_plasticityInstance(ph) associate(prm => param(instance), & @@ -1419,25 +1425,27 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & dot => dotState(instance), & stt => state(instance)) ns = totalNslip(instance) - + tau = 0.0_pReal gdot = 0.0_pReal - - rho = getRho(instance,o,ip,el) + + rho = getRho(instance,o,ip,el) rhoSgl = rho(:,sgl) rhoDip = rho(:,dip) - + rho0 = getRho0(instance,o,ip,el) + my_rhoSgl0 = rho0(:,sgl) + forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(p)%state(iV (s,t,instance),o) endforall - - + + !**************************************************************************** !*** Calculate shear rate - + forall (t = 1:4) & gdot(1:ns,t) = rhoSgl(1:ns,t) * prm%burgers(1:ns) * v(1:ns,t) - + #ifdef DEBUG if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0 & .and. ((debug_e == el .and. debug_i == ip)& @@ -1446,29 +1454,29 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & write(6,'(a,/,4(12x,12(e12.5,1x),/))') '<< CONST >> gdot / 1/s',gdot endif #endif - - - + + + !**************************************************************************** !*** calculate limits for stable dipole height - + do s = 1,ns ! loop over slip systems tau(s) = math_mul33xx33(Mp, prm%Schmid(1:3,1:3,s)) + dst%tau_back(s,o) if (abs(tau(s)) < 1.0e-15_pReal) tau(s) = 1.0e-15_pReal enddo - + dLower = prm%minDipoleHeight dUpper(1:ns,1) = prm%mu * prm%burgers/(8.0_pReal * PI * (1.0_pReal - prm%nu) * abs(tau)) dUpper(1:ns,2) = prm%mu * prm%burgers/(4.0_pReal * PI * abs(tau)) - + where(dNeq0(sqrt(sum(abs(rho(:,edg)),2)))) & dUpper(1:ns,1) = min(1.0_pReal/sqrt(sum(abs(rho(:,edg)),2)),dUpper(1:ns,1)) - + where(dNeq0(sqrt(sum(abs(rho(:,scr)),2)))) & dUpper(1:ns,2) = min(1.0_pReal/sqrt(sum(abs(rho(:,scr)),2)),dUpper(1:ns,2)) dUpper = max(dUpper,dLower) - + !**************************************************************************** !*** calculate dislocation multiplication rhoDotMultiplication = 0.0_pReal @@ -1481,29 +1489,32 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & * sqrt(stt%rho_forest(s,o)) / prm%lambda0(s) ! & ! mean free path ! * 2.0_pReal * sum(abs(v(s,1:2))) / sum(abs(v(s,1:4))) ! ratio of edge to overall velocity determines screw generation endforall - + else isBCC rhoDotMultiplication(1:ns,1:4) = spread( & (sum(abs(gdot(1:ns,1:2)),2) * prm%fEdgeMultiplication + sum(abs(gdot(1:ns,3:4)),2)) & * sqrt(stt%rho_forest(:,o)) / prm%lambda0 / prm%burgers(1:ns), 2, 4) endif isBCC - - + + forall (s = 1:ns, t = 1:4) + v0(s,t) = plasticState(p)%state0(iV(s,t,instance),o) + endforall + !**************************************************************************** !*** calculate dislocation fluxes (only for nonlocal plasticity) rhoDotFlux = 0.0_pReal if (.not. phase_localPlasticity(material_phaseAt(1,el))) then - + !*** check CFL (Courant-Friedrichs-Lewy) condition for flux if (any( abs(gdot) > 0.0_pReal & ! any active slip system ... - .and. prm%CFLfactor * abs(v) * timestep & + .and. prm%CFLfactor * abs(v0) * timestep & > IPvolume(ip,el) / maxval(IParea(:,ip,el)))) then ! ...with velocity above critical value (we use the reference volume and area for simplicity here) #ifdef DEBUG - if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0) then + if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0) then write(6,'(a,i5,a,i2)') '<< CONST >> CFL condition not fullfilled at el ',el,' ip ',ip write(6,'(a,e10.3,a,e10.3)') '<< CONST >> velocity is at ', & - maxval(abs(v), abs(gdot) > 0.0_pReal & - .and. prm%CFLfactor * abs(v) * timestep & + maxval(abs(v0), abs(gdot) > 0.0_pReal & + .and. prm%CFLfactor * abs(v0) * timestep & > IPvolume(ip,el) / maxval(IParea(:,ip,el))), & ' at a timestep of ',timestep write(6,'(a)') '<< CONST >> enforcing cutback !!!' @@ -1512,70 +1523,69 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & plasticState(p)%dotState = IEEE_value(1.0_pReal,IEEE_quiet_NaN) ! -> return NaN and, hence, enforce cutback return endif - - + + !*** be aware of the definition of slip_transverse = slip_direction x slip_normal !!! !*** opposite sign to our p vector in the (s,p,n) triplet !!! - + m(1:3,1:ns,1) = prm%slip_direction m(1:3,1:ns,2) = -prm%slip_direction m(1:3,1:ns,3) = -prm%slip_transverse m(1:3,1:ns,4) = prm%slip_transverse - - my_Fe = Fe(1:3,1:3,1,ip,el) - my_F = matmul(my_Fe, Fp(1:3,1:3,1,ip,el)) - + + my_F = F(1:3,1:3,1,ip,el) + my_Fe = matmul(my_F, math_inv33(Fp(1:3,1:3,1,ip,el))) + neighbors: do n = 1,nIPneighbors - + neighbor_el = IPneighborhood(1,n,ip,el) neighbor_ip = IPneighborhood(2,n,ip,el) neighbor_n = IPneighborhood(3,n,ip,el) np = material_phaseAt(1,neighbor_el) no = material_phasememberAt(1,neighbor_ip,neighbor_el) - + opposite_neighbor = n + mod(n,2) - mod(n+1,2) opposite_el = IPneighborhood(1,opposite_neighbor,ip,el) opposite_ip = IPneighborhood(2,opposite_neighbor,ip,el) opposite_n = IPneighborhood(3,opposite_neighbor,ip,el) - + if (neighbor_n > 0) then ! if neighbor exists, average deformation gradient neighbor_instance = phase_plasticityInstance(material_phaseAt(1,neighbor_el)) - neighbor_Fe = Fe(1:3,1:3,1,neighbor_ip,neighbor_el) - neighbor_F = matmul(neighbor_Fe, Fp(1:3,1:3,1,neighbor_ip,neighbor_el)) + neighbor_F = F(1:3,1:3,1,neighbor_ip,neighbor_el) + neighbor_Fe = matmul(neighbor_F, math_inv33(Fp(1:3,1:3,1,neighbor_ip,neighbor_el))) Favg = 0.5_pReal * (my_F + neighbor_F) else ! if no neighbor, take my value as average Favg = my_F endif - - + + !* FLUX FROM MY NEIGHBOR TO ME - !* This is only considered, if I have a neighbor of nonlocal plasticity - !* (also nonlocal constitutive law with local properties) that is at least a little bit + !* This is only considered, if I have a neighbor of nonlocal plasticity + !* (also nonlocal constitutive law with local properties) that is at least a little bit !* compatible. !* If it's not at all compatible, no flux is arriving, because everything is dammed in front of !* my neighbor's interface. - !* The entering flux from my neighbor will be distributed on my slip systems according to the + !* The entering flux from my neighbor will be distributed on my slip systems according to the !* compatibility - + considerEnteringFlux = .false. - neighbor_v = 0.0_pReal ! needed for check of sign change in flux density below - neighbor_rhoSgl = 0.0_pReal + neighbor_v0 = 0.0_pReal ! needed for check of sign change in flux density below if (neighbor_n > 0) then if (phase_plasticity(material_phaseAt(1,neighbor_el)) == PLASTICITY_NONLOCAL_ID & .and. any(compatibility(:,:,:,n,ip,el) > 0.0_pReal)) & considerEnteringFlux = .true. endif - + enteringFlux: if (considerEnteringFlux) then forall (s = 1:ns, t = 1:4) - neighbor_v(s,t) = plasticState(np)%state(iV (s,t,neighbor_instance),no) - neighbor_rhoSgl(s,t) = max(plasticState(np)%state(iRhoU(s,t,neighbor_instance),no), & + neighbor_v0(s,t) = plasticState(np)%state0(iV (s,t,neighbor_instance),no) + neighbor_rhoSgl0(s,t) = max(plasticState(np)%state0(iRhoU(s,t,neighbor_instance),no), & 0.0_pReal) endforall - - where (neighbor_rhoSgl * IPvolume(neighbor_ip,neighbor_el) ** 0.667_pReal < prm%significantN & - .or. neighbor_rhoSgl < prm%significantRho) & - neighbor_rhoSgl = 0.0_pReal + + where (neighbor_rhoSgl0 * IPvolume(neighbor_ip,neighbor_el) ** 0.667_pReal < prm%significantN & + .or. neighbor_rhoSgl0 < prm%significantRho) & + neighbor_rhoSgl0 = 0.0_pReal normal_neighbor2me_defConf = math_det33(Favg) * matmul(math_inv33(transpose(Favg)), & IPareaNormal(1:3,neighbor_n,neighbor_ip,neighbor_el)) ! calculate the normal of the interface in (average) deformed configuration (now pointing from my neighbor to me!!!) normal_neighbor2me = matmul(transpose(neighbor_Fe), normal_neighbor2me_defConf) & @@ -1586,9 +1596,9 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & do t = 1,4 c = (t + 1) / 2 topp = t + mod(t,2) - mod(t+1,2) - if (neighbor_v(s,t) * math_inner(m(1:3,s,t), normal_neighbor2me) > 0.0_pReal & ! flux from my neighbor to me == entering flux for me - .and. v(s,t) * neighbor_v(s,t) >= 0.0_pReal ) then ! ... only if no sign change in flux density - lineLength = neighbor_rhoSgl(s,t) * neighbor_v(s,t) & + if (neighbor_v0(s,t) * math_inner(m(1:3,s,t), normal_neighbor2me) > 0.0_pReal & ! flux from my neighbor to me == entering flux for me + .and. v0(s,t) * neighbor_v0(s,t) >= 0.0_pReal ) then ! ... only if no sign change in flux density + lineLength = neighbor_rhoSgl0(s,t) * neighbor_v0(s,t) & * math_inner(m(1:3,s,t), normal_neighbor2me) * area ! positive line length that wants to enter through this interface where (compatibility(c,1:ns,s,n,ip,el) > 0.0_pReal) & ! positive compatibility... rhoDotFlux(1:ns,t) = rhoDotFlux(1:ns,t) & @@ -1602,107 +1612,104 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & enddo enddo endif enteringFlux - - + + !* FLUX FROM ME TO MY NEIGHBOR - !* This is not considered, if my opposite neighbor has a different constitutive law than nonlocal (still considered for nonlocal law with local properties). + !* This is not considered, if my opposite neighbor has a different constitutive law than nonlocal (still considered for nonlocal law with local properties). !* Then, we assume, that the opposite(!) neighbor sends an equal amount of dislocations to me. !* So the net flux in the direction of my neighbor is equal to zero: !* leaving flux to neighbor == entering flux from opposite neighbor !* In case of reduced transmissivity, part of the leaving flux is stored as dead dislocation density. !* That means for an interface of zero transmissivity the leaving flux is fully converted to dead dislocations. - + considerLeavingFlux = .true. if (opposite_n > 0) then if (phase_plasticity(material_phaseAt(1,opposite_el)) /= PLASTICITY_NONLOCAL_ID) & considerLeavingFlux = .false. endif - + leavingFlux: if (considerLeavingFlux) then - my_rhoSgl = rhoSgl - my_v = v - normal_me2neighbor_defConf = math_det33(Favg) & - * matmul(math_inv33(transpose(Favg)), & + * matmul(math_inv33(transpose(Favg)), & IPareaNormal(1:3,n,ip,el)) ! calculate the normal of the interface in (average) deformed configuration (pointing from me to my neighbor!!!) normal_me2neighbor = matmul(transpose(my_Fe), normal_me2neighbor_defConf) & / math_det33(my_Fe) ! interface normal in my lattice configuration area = IParea(n,ip,el) * norm2(normal_me2neighbor) - normal_me2neighbor = normal_me2neighbor / norm2(normal_me2neighbor) ! normalize the surface normal to unit length + normal_me2neighbor = normal_me2neighbor / norm2(normal_me2neighbor) ! normalize the surface normal to unit length do s = 1,ns do t = 1,4 c = (t + 1) / 2 - if (my_v(s,t) * math_inner(m(1:3,s,t), normal_me2neighbor) > 0.0_pReal ) then ! flux from me to my neighbor == leaving flux for me (might also be a pure flux from my mobile density to dead density if interface not at all transmissive) - if (my_v(s,t) * neighbor_v(s,t) >= 0.0_pReal) then ! no sign change in flux density + if (v0(s,t) * math_inner(m(1:3,s,t), normal_me2neighbor) > 0.0_pReal ) then ! flux from me to my neighbor == leaving flux for me (might also be a pure flux from my mobile density to dead density if interface not at all transmissive) + if (v0(s,t) * neighbor_v0(s,t) >= 0.0_pReal) then ! no sign change in flux density transmissivity = sum(compatibility(c,1:ns,s,n,ip,el)**2.0_pReal) ! overall transmissivity from this slip system to my neighbor else ! sign change in flux density means sign change in stress which does not allow for dislocations to arive at the neighbor transmissivity = 0.0_pReal endif - lineLength = my_rhoSgl(s,t) * my_v(s,t) & + lineLength = my_rhoSgl0(s,t) * v0(s,t) & * math_inner(m(1:3,s,t), normal_me2neighbor) * area ! positive line length of mobiles that wants to leave through this interface rhoDotFlux(s,t) = rhoDotFlux(s,t) - lineLength / IPvolume(ip,el) ! subtract dislocation flux from current type rhoDotFlux(s,t+4) = rhoDotFlux(s,t+4) & - + lineLength / IPvolume(ip,el) * (1.0_pReal - transmissivity) & - * sign(1.0_pReal, my_v(s,t)) ! dislocation flux that is not able to leave through interface (because of low transmissivity) will remain as immobile single density at the material point + + lineLength / IPvolume(ip,el) * (1.0_pReal - transmissivity) & + * sign(1.0_pReal, v0(s,t)) ! dislocation flux that is not able to leave through interface (because of low transmissivity) will remain as immobile single density at the material point endif enddo enddo endif leavingFlux - + enddo neighbors endif - - - + + + !**************************************************************************** !*** calculate dipole formation and annihilation - + !*** formation by glide - + do c = 1,2 rhoDotSingle2DipoleGlide(1:ns,2*c-1) = -2.0_pReal * dUpper(1:ns,c) / prm%burgers(1:ns) & * (rhoSgl(1:ns,2*c-1) * abs(gdot(1:ns,2*c)) & ! negative mobile --> positive mobile + rhoSgl(1:ns,2*c) * abs(gdot(1:ns,2*c-1)) & ! positive mobile --> negative mobile + abs(rhoSgl(1:ns,2*c+4)) * abs(gdot(1:ns,2*c-1))) ! positive mobile --> negative immobile - + rhoDotSingle2DipoleGlide(1:ns,2*c) = -2.0_pReal * dUpper(1:ns,c) / prm%burgers(1:ns) & * (rhoSgl(1:ns,2*c-1) * abs(gdot(1:ns,2*c)) & ! negative mobile --> positive mobile + rhoSgl(1:ns,2*c) * abs(gdot(1:ns,2*c-1)) & ! positive mobile --> negative mobile + abs(rhoSgl(1:ns,2*c+3)) * abs(gdot(1:ns,2*c))) ! negative mobile --> positive immobile - + rhoDotSingle2DipoleGlide(1:ns,2*c+3) = -2.0_pReal * dUpper(1:ns,c) / prm%burgers(1:ns) & * rhoSgl(1:ns,2*c+3) * abs(gdot(1:ns,2*c)) ! negative mobile --> positive immobile - + rhoDotSingle2DipoleGlide(1:ns,2*c+4) = -2.0_pReal * dUpper(1:ns,c) / prm%burgers(1:ns)& * rhoSgl(1:ns,2*c+4) * abs(gdot(1:ns,2*c-1)) ! positive mobile --> negative immobile - + rhoDotSingle2DipoleGlide(1:ns,c+8) = - rhoDotSingle2DipoleGlide(1:ns,2*c-1) & - rhoDotSingle2DipoleGlide(1:ns,2*c) & + abs(rhoDotSingle2DipoleGlide(1:ns,2*c+3)) & + abs(rhoDotSingle2DipoleGlide(1:ns,2*c+4)) enddo - - + + !*** athermal annihilation - + rhoDotAthermalAnnihilation = 0.0_pReal - - forall (c=1:2) & + + forall (c=1:2) & rhoDotAthermalAnnihilation(1:ns,c+8) = -2.0_pReal * dLower(1:ns,c) / prm%burgers(1:ns) & * ( 2.0_pReal * (rhoSgl(1:ns,2*c-1) * abs(gdot(1:ns,2*c)) + rhoSgl(1:ns,2*c) * abs(gdot(1:ns,2*c-1))) & ! was single hitting single + 2.0_pReal * (abs(rhoSgl(1:ns,2*c+3)) * abs(gdot(1:ns,2*c)) + abs(rhoSgl(1:ns,2*c+4)) * abs(gdot(1:ns,2*c-1))) & ! was single hitting immobile single or was immobile single hit by single + rhoDip(1:ns,c) * (abs(gdot(1:ns,2*c-1)) + abs(gdot(1:ns,2*c)))) ! single knocks dipole constituent ! annihilated screw dipoles leave edge jogs behind on the colinear system - + if (lattice_structure(ph) == LATTICE_fcc_ID) & forall (s = 1:ns, prm%colinearSystem(s) > 0) & rhoDotAthermalAnnihilation(prm%colinearSystem(s),1:2) = - rhoDotAthermalAnnihilation(s,10) & * 0.25_pReal * sqrt(stt%rho_forest(s,o)) * (dUpper(s,2) + dLower(s,2)) * prm%edgeJogFactor - - - + + + !*** thermally activated annihilation of edge dipoles by climb - + rhoDotThermalAnnihilation = 0.0_pReal selfDiffusion = prm%Dsd0 * exp(-prm%selfDiffusionEnergy / (KB * Temperature)) vClimb = prm%atomicVolume * selfDiffusion / ( KB * Temperature ) & @@ -1713,12 +1720,11 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & - rhoDip(s,1) / timestep - rhoDotAthermalAnnihilation(s,9) & - rhoDotSingle2DipoleGlide(s,9)) ! make sure that we do not annihilate more dipoles than we have - rhoDot = 0.0_pReal rhoDot = rhoDotFlux & + rhoDotMultiplication & + rhoDotSingle2DipoleGlide & + rhoDotAthermalAnnihilation & - + rhoDotThermalAnnihilation + + rhoDotThermalAnnihilation #ifdef DEBUG if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0 & @@ -1747,7 +1753,7 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & if ( any(rho(:,mob) + rhoDot(1:ns,1:4) * timestep < -prm%aTolRho) & .or. any(rho(:,dip) + rhoDot(1:ns,9:10) * timestep < -prm%aTolRho)) then #ifdef DEBUG - if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0) then + if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0) then write(6,'(a,i5,a,i2)') '<< CONST >> evolution rate leads to negative density at el ',el,' ip ',ip write(6,'(a)') '<< CONST >> enforcing cutback !!!' endif @@ -1765,21 +1771,21 @@ end subroutine plastic_nonlocal_dotState !-------------------------------------------------------------------------------------------------- !> @brief Compatibility update -!> @detail Compatibility is defined as normalized product of signed cosine of the angle between the slip +!> @detail Compatibility is defined as normalized product of signed cosine of the angle between the slip ! plane normals and signed cosine of the angle between the slip directions. Only the largest values ! that sum up to a total of 1 are considered, all others are set to zero. !-------------------------------------------------------------------------------------------------- -module subroutine plastic_nonlocal_updateCompatibility(orientation,i,e) - +module subroutine plastic_nonlocal_updateCompatibility(orientation,i,e) + integer, intent(in) :: & i, & e type(rotation), dimension(1,discretization_nIP,discretization_nElem), intent(in) :: & orientation ! crystal orientation - + integer :: & Nneighbors, & ! number of neighbors - n, & ! neighbor index + n, & ! neighbor index neighbor_e, & ! element index of my neighbor neighbor_i, & ! integration point index of my neighbor ph, & @@ -1792,13 +1798,13 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,i,e) s2 ! slip system index (my neighbor) real(pReal), dimension(2,totalNslip(phase_plasticityInstance(material_phaseAt(1,e))),& totalNslip(phase_plasticityInstance(material_phaseAt(1,e))),& - nIPneighbors) :: & - my_compatibility ! my_compatibility for current element and ip + nIPneighbors) :: & + my_compatibility ! my_compatibility for current element and ip real(pReal) :: & my_compatibilitySum, & thresholdValue, & nThresholdValues - logical, dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,e)))) :: & + logical, dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,e)))) :: & belowThreshold type(rotation) :: mis @@ -1808,32 +1814,32 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,i,e) instance = phase_plasticityInstance(ph) ns = totalNslip(instance) associate(prm => param(instance)) - + !*** start out fully compatible my_compatibility = 0.0_pReal - - forall(s1 = 1:ns) my_compatibility(1:2,s1,s1,1:Nneighbors) = 1.0_pReal - + + forall(s1 = 1:ns) my_compatibility(1:2,s1,s1,1:Nneighbors) = 1.0_pReal + !*** Loop thrugh neighbors and check whether there is any compatibility. - + neighbors: do n = 1,Nneighbors neighbor_e = IPneighborhood(1,n,i,e) neighbor_i = IPneighborhood(2,n,i,e) - - + + !* FREE SURFACE !* Set surface transmissivity to the value specified in the material.config - + if (neighbor_e <= 0 .or. neighbor_i <= 0) then forall(s1 = 1:ns) my_compatibility(1:2,s1,s1,n) = sqrt(prm%surfaceTransmissivity) cycle endif - - + + !* PHASE BOUNDARY - !* If we encounter a different nonlocal phase at the neighbor, + !* If we encounter a different nonlocal phase at the neighbor, !* we consider this to be a real "physical" phase boundary, so completely incompatible. - !* If one of the two phases has a local plasticity law, + !* If one of the two phases has a local plasticity law, !* we do not consider this to be a phase boundary, so completely compatible. neighbor_phase = material_phaseAt(1,neighbor_e) if (neighbor_phase /= ph) then @@ -1841,8 +1847,8 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,i,e) forall(s1 = 1:ns) my_compatibility(1:2,s1,s1,n) = 0.0_pReal cycle endif - - + + !* GRAIN BOUNDARY ! !* fixed transmissivity for adjacent ips with different texture (only if explicitly given in material.config) if (prm%grainboundaryTransmissivity >= 0.0_pReal) then @@ -1854,16 +1860,16 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,i,e) endif cycle endif - - + + !* GRAIN BOUNDARY ? !* Compatibility defined by relative orientation of slip systems: !* The my_compatibility value is defined as the product of the slip normal projection and the slip direction projection. - !* Its sign is always positive for screws, for edges it has the same sign as the slip normal projection. - !* Since the sum for each slip system can easily exceed one (which would result in a transmissivity larger than one), + !* Its sign is always positive for screws, for edges it has the same sign as the slip normal projection. + !* Since the sum for each slip system can easily exceed one (which would result in a transmissivity larger than one), !* only values above or equal to a certain threshold value are considered. This threshold value is chosen, such that - !* the number of compatible slip systems is minimized with the sum of the original compatibility values exceeding one. - !* Finally the smallest compatibility value is decreased until the sum is exactly equal to one. + !* the number of compatible slip systems is minimized with the sum of the original compatibility values exceeding one. + !* Finally the smallest compatibility value is decreased until the sum is exactly equal to one. !* All values below the threshold are set to zero. else mis = orientation(1,i,e)%misorientation(orientation(1,neighbor_i,neighbor_e)) @@ -1878,7 +1884,7 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,i,e) * abs(math_inner(prm%slip_direction(1:3,s1), & mis%rotate(prm%slip_direction(1:3,s2)))) enddo neighborSlipSystems - + my_compatibilitySum = 0.0_pReal belowThreshold = .true. do while (my_compatibilitySum < 1.0_pReal .and. any(belowThreshold(1:ns))) @@ -1896,33 +1902,33 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,i,e) where (belowThreshold(1:ns)) my_compatibility(2,1:ns,s1,n) = 0.0_pReal enddo mySlipSystems endif - + enddo neighbors - + compatibility(1:2,1:ns,1:ns,1:Nneighbors,i,e) = my_compatibility - + end associate - + end subroutine plastic_nonlocal_updateCompatibility - + !-------------------------------------------------------------------------------------------------- !> @brief returns copy of current dislocation densities from state !> @details raw values is rectified !-------------------------------------------------------------------------------------------------- function getRho(instance,of,ip,el) - + integer, intent(in) :: instance, of,ip,el real(pReal), dimension(param(instance)%totalNslip,10) :: getRho - + associate(prm => param(instance)) getRho = reshape(state(instance)%rho(:,of),[prm%totalNslip,10]) - + ! ensure positive densities (not for imm, they have a sign) getRho(:,mob) = max(getRho(:,mob),0.0_pReal) getRho(:,dip) = max(getRho(:,dip),0.0_pReal) - + where(abs(getRho) < max(prm%significantN/IPvolume(ip,el)**(2.0_pReal/3.0_pReal),prm%significantRho)) & getRho = 0.0_pReal @@ -1931,6 +1937,31 @@ function getRho(instance,of,ip,el) end function getRho +!-------------------------------------------------------------------------------------------------- +!> @brief returns copy of current dislocation densities from state +!> @details raw values is rectified +!-------------------------------------------------------------------------------------------------- +function getRho0(instance,of,ip,el) + + integer, intent(in) :: instance, of,ip,el + real(pReal), dimension(param(instance)%totalNslip,10) :: getRho0 + + associate(prm => param(instance)) + + getRho0 = reshape(state0(instance)%rho(:,of),[prm%totalNslip,10]) + + ! ensure positive densities (not for imm, they have a sign) + getRho0(:,mob) = max(getRho0(:,mob),0.0_pReal) + getRho0(:,dip) = max(getRho0(:,dip),0.0_pReal) + + where(abs(getRho0) < max(prm%significantN/IPvolume(ip,el)**(2.0_pReal/3.0_pReal),prm%significantRho)) & + getRho0 = 0.0_pReal + + end associate + +end function getRho0 + + !-------------------------------------------------------------------------------------------------- !> @brief writes results to HDF5 output file !-------------------------------------------------------------------------------------------------- diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 6e0a7f614..00d6f88b1 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -264,8 +264,8 @@ subroutine crystallite_init do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1),FEsolving_execIP(2) do c = 1,homogenization_Ngrains(material_homogenizationAt(e)) - call constitutive_microstructure(crystallite_Fe(1:3,1:3,c,i,e), & - crystallite_Fp(1:3,1:3,c,i,e), & + call constitutive_microstructure(crystallite_partionedF0(1:3,1:3,c,i,e), & + crystallite_partionedFp0(1:3,1:3,c,i,e), & c,i,e) ! update dependent state variables to be consistent with basic states enddo enddo @@ -1963,9 +1963,9 @@ subroutine update_dotState(timeFraction) !$OMP FLUSH(nonlocalStop) if ((crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) .and. .not. nonlocalStop) then call constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & - crystallite_Fe, & + crystallite_partionedF0, & crystallite_Fi(1:3,1:3,g,i,e), & - crystallite_Fp, & + crystallite_partionedFp0, & crystallite_subdt(g,i,e)*timeFraction, g,i,e) p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e) NaN = any(IEEE_is_NaN(plasticState(p)%dotState(:,c))) diff --git a/src/debug.f90 b/src/debug.f90 index 86bcfcd29..a0947c410 100644 --- a/src/debug.f90 +++ b/src/debug.f90 @@ -37,10 +37,9 @@ module debug debug_HOMOGENIZATION = 9, & debug_CPFEM = 10, & debug_SPECTRAL = 11, & - debug_MARC = 12, & - debug_ABAQUS = 13 + debug_MARC = 12 integer, parameter, private :: & - debug_MAXNTYPE = debug_ABAQUS !< must be set to the maximum defined debug type + debug_MAXNTYPE = debug_MARC !< must be set to the maximum defined debug type integer,protected, dimension(debug_maxNtype+2), public :: & ! specific ones, and 2 for "all" and "other" debug_level = 0 @@ -136,8 +135,6 @@ subroutine debug_init what = debug_SPECTRAL case ('marc') what = debug_MARC - case ('abaqus') - what = debug_ABAQUS case ('all') what = debug_MAXNTYPE + 1 case ('other') @@ -210,8 +207,6 @@ subroutine debug_init tag = ' Spectral solver' case (debug_MARC) tag = ' MSC.MARC FEM solver' - case (debug_ABAQUS) - tag = ' ABAQUS FEM solver' end select if(debug_level(i) /= 0) then diff --git a/src/mesh_abaqus.f90 b/src/mesh_abaqus.f90 deleted file mode 100644 index f4f5113e8..000000000 --- a/src/mesh_abaqus.f90 +++ /dev/null @@ -1,2822 +0,0 @@ -!-------------------------------------------------------------------------------------------------- -!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH -!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH -!> @author Christoph Koords, Max-Planck-Institut für Eisenforschung GmbH -!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH -!> @brief Sets up the mesh for the solvers MSC.Marc, Abaqus and the spectral solver -!-------------------------------------------------------------------------------------------------- -module mesh - use prec - use geometry_plastic_nonlocal - use discretization - use math - - implicit none - private - - integer, public, protected :: & - mesh_NcpElems, & !< total number of CP elements in local mesh - mesh_elemType, & !< Element type of the mesh (only support homogeneous meshes) - mesh_Nnodes, & !< total number of nodes in mesh - mesh_Ncellnodes, & !< total number of cell nodes in mesh (including duplicates) - mesh_Ncells, & !< total number of cells in mesh - mesh_maxNipNeighbors, & !< max number of IP neighbors in any CP element - mesh_maxNsharedElems !< max number of CP elements sharing a node -!!!! BEGIN DEPRECATED !!!!! - integer, public, protected :: & - mesh_maxNips, & !< max number of IPs in any CP element - mesh_maxNcellnodes !< max number of cell nodes in any CP element -!!!! BEGIN DEPRECATED !!!!! - - integer, dimension(:,:), allocatable, public, protected :: & - mesh_element, & !DEPRECATED - mesh_sharedElem, & !< entryCount and list of elements containing node - mesh_nodeTwins !< node twins are surface nodes that lie exactly on opposite sides of the mesh (surfaces nodes with equal coordinate values in two dimensions) - - integer, dimension(:,:,:,:), allocatable, public, protected :: & - mesh_ipNeighborhood !< 6 or less neighboring IPs as [element_num, IP_index, neighbor_index that points to me] - - real(pReal), public, protected :: & - mesh_unitlength !< physical length of one unit in mesh - - real(pReal), dimension(:,:), allocatable, public :: & - mesh_node, & !< node x,y,z coordinates (after deformation! ONLY FOR MARC!!!) - mesh_cellnode !< cell node x,y,z coordinates (after deformation! ONLY FOR MARC!!!) - - real(pReal), dimension(:,:), allocatable, public, protected :: & - mesh_ipVolume, & !< volume associated with IP (initially!) - mesh_node0 !< node x,y,z coordinates (initially!) - - real(pReal), dimension(:,:,:), allocatable, public, protected :: & - mesh_ipArea !< area of interface to neighboring IP (initially!) - - real(pReal), dimension(:,:,:), allocatable, public :: & - mesh_ipCoordinates !< IP x,y,z coordinates (after deformation!) - - real(pReal),dimension(:,:,:,:), allocatable, public, protected :: & - mesh_ipAreaNormal !< area normal of interface to neighboring IP (initially!) - - logical, dimension(3), public, protected :: mesh_periodicSurface !< flag indicating periodic outer surfaces (used for fluxes) - - integer, private :: & - mesh_maxNelemInSet, & - mesh_Nmaterials - - integer, dimension(2), private :: & - mesh_maxValStateVar = 0 - -integer, dimension(:,:), allocatable, private :: & - mesh_cellnodeParent !< cellnode's parent element ID, cellnode's intra-element ID - - integer,dimension(:,:,:), allocatable, private :: & - mesh_cell !< cell connectivity for each element,ip/cell - - integer, dimension(:,:,:), allocatable, private :: & - FE_nodesAtIP, & !< map IP index to node indices in a specific type of element - FE_ipNeighbor, & !< +x,-x,+y,-y,+z,-z list of intra-element IPs and(negative) neighbor faces per own IP in a specific type of element - FE_cell, & !< list of intra-element cell node IDs that constitute the cells in a specific type of element geometry - FE_cellface !< list of intra-cell cell node IDs that constitute the cell faces of a specific type of cell - - real(pReal), dimension(:,:,:), allocatable, private :: & - FE_cellnodeParentnodeWeights !< list of node weights for the generation of cell nodes - - integer, dimension(:,:,:,:), allocatable, private :: & - FE_subNodeOnIPFace - -! These definitions should actually reside in the FE-solver specific part (different for MARC/ABAQUS) -! Hence, I suggest to prefix with "FE_" - - integer, parameter, public :: & - FE_Nelemtypes = 13, & - FE_Ngeomtypes = 10, & - FE_Ncelltypes = 4, & - FE_maxNnodes = 20, & - FE_maxNips = 27, & - FE_maxNipNeighbors = 6, & - FE_maxmaxNnodesAtIP = 8, & !< max number of (equivalent) nodes attached to an IP - FE_maxNmatchingNodesPerFace = 4, & - FE_maxNfaces = 6, & - FE_maxNcellnodes = 64, & - FE_maxNcellnodesPerCell = 8, & - FE_maxNcellfaces = 6, & - FE_maxNcellnodesPerCellface = 4 - - integer, dimension(FE_Nelemtypes), parameter, public :: FE_geomtype = & !< geometry type of particular element type - int([ & - 1, & ! element 6 (2D 3node 1ip) - 2, & ! element 125 (2D 6node 3ip) - 3, & ! element 11 (2D 4node 4ip) - 4, & ! element 27 (2D 8node 9ip) - 3, & ! element 54 (2D 8node 4ip) - 5, & ! element 134 (3D 4node 1ip) - 6, & ! element 157 (3D 5node 4ip) - 6, & ! element 127 (3D 10node 4ip) - 7, & ! element 136 (3D 6node 6ip) - 8, & ! element 117 (3D 8node 1ip) - 9, & ! element 7 (3D 8node 8ip) - 9, & ! element 57 (3D 20node 8ip) - 10 & ! element 21 (3D 20node 27ip) - ],pInt) - - integer, dimension(FE_Ngeomtypes), parameter, public :: FE_celltype = & !< cell type that is used by each geometry type - int([ & - 1, & ! element 6 (2D 3node 1ip) - 2, & ! element 125 (2D 6node 3ip) - 2, & ! element 11 (2D 4node 4ip) - 2, & ! element 27 (2D 8node 9ip) - 3, & ! element 134 (3D 4node 1ip) - 4, & ! element 127 (3D 10node 4ip) - 4, & ! element 136 (3D 6node 6ip) - 4, & ! element 117 (3D 8node 1ip) - 4, & ! element 7 (3D 8node 8ip) - 4 & ! element 21 (3D 20node 27ip) - ],pInt) - - integer, dimension(FE_Ngeomtypes), parameter, public :: FE_dimension = & !< dimension of geometry type - int([ & - 2, & ! element 6 (2D 3node 1ip) - 2, & ! element 125 (2D 6node 3ip) - 2, & ! element 11 (2D 4node 4ip) - 2, & ! element 27 (2D 8node 9ip) - 3, & ! element 134 (3D 4node 1ip) - 3, & ! element 127 (3D 10node 4ip) - 3, & ! element 136 (3D 6node 6ip) - 3, & ! element 117 (3D 8node 1ip) - 3, & ! element 7 (3D 8node 8ip) - 3 & ! element 21 (3D 20node 27ip) - ],pInt) - - integer, dimension(FE_Nelemtypes), parameter, public :: FE_Nnodes = & !< number of nodes that constitute a specific type of element - int([ & - 3, & ! element 6 (2D 3node 1ip) - 6, & ! element 125 (2D 6node 3ip) - 4, & ! element 11 (2D 4node 4ip) - 8, & ! element 27 (2D 8node 9ip) - 8, & ! element 54 (2D 8node 4ip) - 4, & ! element 134 (3D 4node 1ip) - 5, & ! element 157 (3D 5node 4ip) - 10, & ! element 127 (3D 10node 4ip) - 6, & ! element 136 (3D 6node 6ip) - 8, & ! element 117 (3D 8node 1ip) - 8, & ! element 7 (3D 8node 8ip) - 20, & ! element 57 (3D 20node 8ip) - 20 & ! element 21 (3D 20node 27ip) - ],pInt) - - integer, dimension(FE_Ngeomtypes), parameter, public :: FE_Nfaces = & !< number of faces of a specific type of element geometry - int([ & - 3, & ! element 6 (2D 3node 1ip) - 3, & ! element 125 (2D 6node 3ip) - 4, & ! element 11 (2D 4node 4ip) - 4, & ! element 27 (2D 8node 9ip) - 4, & ! element 134 (3D 4node 1ip) - 4, & ! element 127 (3D 10node 4ip) - 5, & ! element 136 (3D 6node 6ip) - 6, & ! element 117 (3D 8node 1ip) - 6, & ! element 7 (3D 8node 8ip) - 6 & ! element 21 (3D 20node 27ip) - ],pInt) - - integer, dimension(FE_Ngeomtypes), parameter, private :: FE_NmatchingNodes = & !< number of nodes that are needed for face matching in a specific type of element geometry - int([ & - 3, & ! element 6 (2D 3node 1ip) - 3, & ! element 125 (2D 6node 3ip) - 4, & ! element 11 (2D 4node 4ip) - 4, & ! element 27 (2D 8node 9ip) - 4, & ! element 134 (3D 4node 1ip) - 4, & ! element 127 (3D 10node 4ip) - 6, & ! element 136 (3D 6node 6ip) - 8, & ! element 117 (3D 8node 1ip) - 8, & ! element 7 (3D 8node 8ip) - 8 & ! element 21 (3D 20node 27ip) - ],pInt) - - integer, dimension(FE_maxNfaces,FE_Ngeomtypes), parameter, private :: FE_NmatchingNodesPerFace = & !< number of matching nodes per face in a specific type of element geometry - reshape(int([ & - 2,2,2,0,0,0, & ! element 6 (2D 3node 1ip) - 2,2,2,0,0,0, & ! element 125 (2D 6node 3ip) - 2,2,2,2,0,0, & ! element 11 (2D 4node 4ip) - 2,2,2,2,0,0, & ! element 27 (2D 8node 9ip) - 3,3,3,3,0,0, & ! element 134 (3D 4node 1ip) - 3,3,3,3,0,0, & ! element 127 (3D 10node 4ip) - 3,4,4,4,3,0, & ! element 136 (3D 6node 6ip) - 4,4,4,4,4,4, & ! element 117 (3D 8node 1ip) - 4,4,4,4,4,4, & ! element 7 (3D 8node 8ip) - 4,4,4,4,4,4 & ! element 21 (3D 20node 27ip) - ],pInt),[FE_maxNipNeighbors,FE_Ngeomtypes]) - - integer, dimension(FE_maxNmatchingNodesPerFace,FE_maxNfaces,FE_Ngeomtypes), & - parameter, private :: FE_face = & !< List of node indices on each face of a specific type of element geometry - reshape(int([& - 1,2,0,0 , & ! element 6 (2D 3node 1ip) - 2,3,0,0 , & - 3,1,0,0 , & - 0,0,0,0 , & - 0,0,0,0 , & - 0,0,0,0 , & - 1,2,0,0 , & ! element 125 (2D 6node 3ip) - 2,3,0,0 , & - 3,1,0,0 , & - 0,0,0,0 , & - 0,0,0,0 , & - 0,0,0,0 , & - 1,2,0,0 , & ! element 11 (2D 4node 4ip) - 2,3,0,0 , & - 3,4,0,0 , & - 4,1,0,0 , & - 0,0,0,0 , & - 0,0,0,0 , & - 1,2,0,0 , & ! element 27 (2D 8node 9ip) - 2,3,0,0 , & - 3,4,0,0 , & - 4,1,0,0 , & - 0,0,0,0 , & - 0,0,0,0 , & - 1,2,3,0 , & ! element 134 (3D 4node 1ip) - 1,4,2,0 , & - 2,3,4,0 , & - 1,3,4,0 , & - 0,0,0,0 , & - 0,0,0,0 , & - 1,2,3,0 , & ! element 127 (3D 10node 4ip) - 1,4,2,0 , & - 2,4,3,0 , & - 1,3,4,0 , & - 0,0,0,0 , & - 0,0,0,0 , & - 1,2,3,0 , & ! element 136 (3D 6node 6ip) - 1,4,5,2 , & - 2,5,6,3 , & - 1,3,6,4 , & - 4,6,5,0 , & - 0,0,0,0 , & - 1,2,3,4 , & ! element 117 (3D 8node 1ip) - 2,1,5,6 , & - 3,2,6,7 , & - 4,3,7,8 , & - 4,1,5,8 , & - 8,7,6,5 , & - 1,2,3,4 , & ! element 7 (3D 8node 8ip) - 2,1,5,6 , & - 3,2,6,7 , & - 4,3,7,8 , & - 4,1,5,8 , & - 8,7,6,5 , & - 1,2,3,4 , & ! element 21 (3D 20node 27ip) - 2,1,5,6 , & - 3,2,6,7 , & - 4,3,7,8 , & - 4,1,5,8 , & - 8,7,6,5 & - ],pInt),[FE_maxNmatchingNodesPerFace,FE_maxNfaces,FE_Ngeomtypes]) - - integer, dimension(FE_Ngeomtypes), parameter, private :: FE_Ncellnodes = & !< number of cell nodes in a specific geometry type - int([ & - 3, & ! element 6 (2D 3node 1ip) - 7, & ! element 125 (2D 6node 3ip) - 9, & ! element 11 (2D 4node 4ip) - 16, & ! element 27 (2D 8node 9ip) - 4, & ! element 134 (3D 4node 1ip) - 15, & ! element 127 (3D 10node 4ip) - 21, & ! element 136 (3D 6node 6ip) - 8, & ! element 117 (3D 8node 1ip) - 27, & ! element 7 (3D 8node 8ip) - 64 & ! element 21 (3D 20node 27ip) - ],pInt) - - integer, dimension(FE_Ncelltypes), parameter, private :: FE_NcellnodesPerCell = & !< number of cell nodes in a specific cell type - int([ & - 3, & ! (2D 3node) - 4, & ! (2D 4node) - 4, & ! (3D 4node) - 8 & ! (3D 8node) - ],pInt) - - integer, dimension(FE_Ncelltypes), parameter, private :: FE_NcellnodesPerCellface = & !< number of cell nodes per cell face in a specific cell type - int([& - 2, & ! (2D 3node) - 2, & ! (2D 4node) - 3, & ! (3D 4node) - 4 & ! (3D 8node) - ],pInt) - - integer, dimension(FE_Ngeomtypes), parameter, public :: FE_Nips = & !< number of IPs in a specific type of element - int([ & - 1, & ! element 6 (2D 3node 1ip) - 3, & ! element 125 (2D 6node 3ip) - 4, & ! element 11 (2D 4node 4ip) - 9, & ! element 27 (2D 8node 9ip) - 1, & ! element 134 (3D 4node 1ip) - 4, & ! element 127 (3D 10node 4ip) - 6, & ! element 136 (3D 6node 6ip) - 1, & ! element 117 (3D 8node 1ip) - 8, & ! element 7 (3D 8node 8ip) - 27 & ! element 21 (3D 20node 27ip) - ],pInt) - - integer, dimension(FE_Ncelltypes), parameter, public :: FE_NipNeighbors = & !< number of ip neighbors / cell faces in a specific cell type - int([& - 3, & ! (2D 3node) - 4, & ! (2D 4node) - 4, & ! (3D 4node) - 6 & ! (3D 8node) - ],pInt) - - integer, dimension(FE_Ngeomtypes), parameter, private :: FE_maxNnodesAtIP = & !< maximum number of parent nodes that belong to an IP for a specific type of element - int([ & - 3, & ! element 6 (2D 3node 1ip) - 1, & ! element 125 (2D 6node 3ip) - 1, & ! element 11 (2D 4node 4ip) - 2, & ! element 27 (2D 8node 9ip) - 4, & ! element 134 (3D 4node 1ip) - 1, & ! element 127 (3D 10node 4ip) - 1, & ! element 136 (3D 6node 6ip) - 8, & ! element 117 (3D 8node 1ip) - 1, & ! element 7 (3D 8node 8ip) - 4 & ! element 21 (3D 20node 27ip) - ],pInt) - - integer, private :: & - mesh_Nelems, & !< total number of elements in mesh (including non-DAMASK elements) - mesh_maxNnodes, & !< max number of nodes in any CP element - mesh_NelemSets - character(len=64), dimension(:), allocatable, private :: & - mesh_nameElemSet, & !< names of elementSet - mesh_nameMaterial, & !< names of material in solid section - mesh_mapMaterial !< name of elementSet for material - integer, dimension(:,:), allocatable, private :: & - mesh_mapElemSet !< list of elements in elementSet - integer, dimension(:,:), allocatable, target, private :: & - mesh_mapFEtoCPelem, & !< [sorted FEid, corresponding CPid] - mesh_mapFEtoCPnode !< [sorted FEid, corresponding CPid] - logical, private :: noPart !< for cases where the ABAQUS input file does not use part/assembly information - - public :: & - mesh_init, & - mesh_build_cellnodes, & - mesh_build_ipVolumes, & - mesh_build_ipCoordinates, & - mesh_cellCenterCoordinates, & - mesh_FEasCP - - private :: & - mesh_get_damaskOptions, & - mesh_build_cellconnectivity, & - mesh_build_ipAreas, & - FE_mapElemtype, & - mesh_build_FEdata, & - mesh_build_nodeTwins, & - mesh_build_sharedElems, & - mesh_build_ipNeighborhood, & - mesh_abaqus_count_nodesAndElements, & - mesh_abaqus_count_elementSets, & - mesh_abaqus_count_materials, & - mesh_abaqus_map_elementSets, & - mesh_abaqus_map_materials, & - mesh_abaqus_count_cpElements, & - mesh_abaqus_map_elements, & - mesh_abaqus_map_nodes, & - mesh_abaqus_build_nodes, & - mesh_abaqus_count_cpSizes, & - mesh_abaqus_build_elements - - - type, public, extends(tMesh) :: tMesh_abaqus - - integer:: & - mesh_Nelems, & !< total number of elements in mesh (including non-DAMASK elements) - mesh_maxNnodes, & !< max number of nodes in any CP element - mesh_NelemSets, & - mesh_maxNelemInSet, & - mesh_Nmaterials - character(len=64), dimension(:), allocatable :: & - mesh_nameElemSet, & !< names of elementSet - mesh_nameMaterial, & !< names of material in solid section - mesh_mapMaterial !< name of elementSet for material - integer, dimension(:,:), allocatable :: & - mesh_mapElemSet !< list of elements in elementSet - logical:: noPart !< for cases where the ABAQUS input file does not use part/assembly information - - contains - procedure, pass(self) :: tMesh_abaqus_init - generic, public :: init => tMesh_abaqus_init - end type tMesh_abaqus - - type(tMesh_abaqus), public, protected :: theMesh - -contains - -subroutine tMesh_abaqus_init(self,elemType,nodes) - - - class(tMesh_abaqus) :: self - real(pReal), dimension(:,:), intent(in) :: nodes - integer, intent(in) :: elemType - - call self%tMesh%init('mesh',elemType,nodes) - -end subroutine tMesh_abaqus_init - -!-------------------------------------------------------------------------------------------------- -!> @brief initializes the mesh by calling all necessary private routines the mesh module -!! Order and routines strongly depend on type of solver -!-------------------------------------------------------------------------------------------------- -subroutine mesh_init(ip,el) - use DAMASK_interface - use IO, only: & - IO_open_InputFile, & - IO_error - use debug, only: & - debug_e, & - debug_i, & - debug_level, & - debug_mesh, & - debug_levelBasic - use numerics, only: & - usePingPong, & - numerics_unitlength, & - worldrank - use FEsolving, only: & - modelName, & - calcMode, & FEsolving_execElem, & - FEsolving_execIP - - - integer, parameter :: FILEUNIT = 222 - integer, intent(in), optional :: el, ip - integer :: j - logical :: myDebug - - write(6,'(/,a)') ' <<<+- mesh init -+>>>' - - mesh_unitlength = numerics_unitlength ! set physical extent of a length unit in mesh - - myDebug = (iand(debug_level(debug_mesh),debug_levelBasic) /= 0) - - call IO_open_inputFile(FILEUNIT) ! parse info from input file... - if (myDebug) write(6,'(a)') ' Opened input file'; flush(6) - noPart = hasNoPart(FILEUNIT) - call mesh_abaqus_count_nodesAndElements(FILEUNIT) - if (myDebug) write(6,'(a)') ' Counted nodes/elements'; flush(6) - call mesh_abaqus_count_elementSets(FILEUNIT) - if (myDebug) write(6,'(a)') ' Counted element sets'; flush(6) - call mesh_abaqus_count_materials(FILEUNIT) - if (myDebug) write(6,'(a)') ' Counted materials'; flush(6) - call mesh_abaqus_map_elementSets(FILEUNIT) - if (myDebug) write(6,'(a)') ' Mapped element sets'; flush(6) - call mesh_abaqus_map_materials(FILEUNIT) - if (myDebug) write(6,'(a)') ' Mapped materials'; flush(6) - call mesh_abaqus_count_cpElements(FILEUNIT) - if (myDebug) write(6,'(a)') ' Counted CP elements'; flush(6) - call mesh_abaqus_map_elements(FILEUNIT) - if (myDebug) write(6,'(a)') ' Mapped elements'; flush(6) - call mesh_abaqus_map_nodes(FILEUNIT) - if (myDebug) write(6,'(a)') ' Mapped nodes'; flush(6) - call mesh_abaqus_build_nodes(FILEUNIT) - if (myDebug) write(6,'(a)') ' Built nodes'; flush(6) - call mesh_abaqus_count_cpSizes(FILEUNIT) - if (myDebug) write(6,'(a)') ' Counted CP sizes'; flush(6) - call mesh_abaqus_build_elements(FILEUNIT) - if (myDebug) write(6,'(a)') ' Built elements'; flush(6) - call mesh_get_damaskOptions(mesh_periodicSurface,FILEUNIT) - if (myDebug) write(6,'(a)') ' Got DAMASK options'; flush(6) - close (FILEUNIT) - - call theMesh%init(mesh_element(2,1),mesh_node0) - call theMesh%setNelems(mesh_NcpElems) - call mesh_build_FEdata ! get properties of the different types of elements - - call mesh_build_cellconnectivity - if (myDebug) write(6,'(a)') ' Built cell connectivity'; flush(6) - mesh_cellnode = mesh_build_cellnodes(mesh_node,mesh_Ncellnodes) - if (myDebug) write(6,'(a)') ' Built cell nodes'; flush(6) - call mesh_build_ipCoordinates - if (myDebug) write(6,'(a)') ' Built IP coordinates'; flush(6) - call mesh_build_ipVolumes - if (myDebug) write(6,'(a)') ' Built IP volumes'; flush(6) - call mesh_build_ipAreas - if (myDebug) write(6,'(a)') ' Built IP areas'; flush(6) - call mesh_build_nodeTwins - if (myDebug) write(6,'(a)') ' Built node twins'; flush(6) - call mesh_build_sharedElems - if (myDebug) write(6,'(a)') ' Built shared elements'; flush(6) - call mesh_build_ipNeighborhood - if (myDebug) write(6,'(a)') ' Built IP neighborhood'; flush(6) - - if (usePingPong .and. (mesh_Nelems /= mesh_NcpElems)) & - call IO_error(600) ! ping-pong must be disabled when having non-DAMASK elements - if (debug_e < 1 .or. debug_e > mesh_NcpElems) & - call IO_error(602,ext_msg='element') ! selected element does not exist - if (debug_i < 1 .or. debug_i > FE_Nips(FE_geomtype(mesh_element(2,debug_e)))) & - call IO_error(602,ext_msg='IP') ! selected element does not have requested IP - FEsolving_execElem = [ 1,mesh_NcpElems ] ! parallel loop bounds set to comprise all DAMASK elements - allocate(FEsolving_execIP(2,mesh_NcpElems), source=1) ! parallel loop bounds set to comprise from first IP... - forall (j = 1:mesh_NcpElems) FEsolving_execIP(2,j) = FE_Nips(FE_geomtype(mesh_element(2,j))) ! ...up to own IP count for each element - allocate(calcMode(mesh_maxNips,mesh_NcpElems)) - calcMode = .false. ! pretend to have collected what first call is asking (F = I) - calcMode(ip,mesh_FEasCP('elem',el)) = .true. ! first ip,el needs to be already pingponged to "calc" - - -! better name - theMesh%homogenizationAt = mesh_element(3,:) - theMesh%microstructureAt = mesh_element(4,:) - - call discretization_init(mesh_element(3,:),mesh_element(4,:),& - reshape(mesh_ipCoordinates,[3,theMesh%elem%nIPs*theMesh%nElems]),& - mesh_node0) - call geometry_plastic_nonlocal_setIPvolume(mesh_ipVolume) - call geometry_plastic_nonlocal_setIPneighborhood(mesh_ipNeighborhood) - call geometry_plastic_nonlocal_setIParea(mesh_IParea) - call geometry_plastic_nonlocal_setIPareaNormal(mesh_IPareaNormal) - -contains - - -!-------------------------------------------------------------------------------------------------- -!> @brief check if the input file for Abaqus contains part info -!-------------------------------------------------------------------------------------------------- -logical function hasNoPart(fileUnit) - use IO, only: & - IO_stringPos, & - IO_stringValue, & - IO_lc - - - integer, intent(in) :: fileUnit - - integer, allocatable, dimension(:) :: chunkPos - character(len=65536) :: line - - hasNoPart = .true. - - rewind(fileUnit) - do - read(fileUnit,'(a65536)',END=620) line - chunkPos = IO_stringPos(line) - if (IO_lc(IO_stringValue(line,chunkPos,1)) == '*part' ) then - hasNoPart = .false. - exit - endif - enddo - -620 end function hasNoPart - -end subroutine mesh_init - - - - - - - - -!-------------------------------------------------------------------------------------------------- -!> @brief Count overall number of nodes and elements in mesh and stores them in -!! 'mesh_Nelems' and 'mesh_Nnodes' -!-------------------------------------------------------------------------------------------------- -subroutine mesh_abaqus_count_nodesAndElements(fileUnit) - - use IO, only: IO_lc, & - IO_stringValue, & - IO_stringPos, & - IO_countDataLines, & - IO_error - - - integer, intent(in) :: fileUnit - - integer, allocatable, dimension(:) :: chunkPos - character(len=300) :: line - integer :: myStat - logical :: inPart - - mesh_Nnodes = 0 - mesh_Nelems = 0 - - inPart = .false. - myStat = 0 - rewind(fileUnit) - do while(myStat == 0) - read (fileUnit,'(a300)',iostat=myStat) line - chunkPos = IO_stringPos(line) - if ( IO_lc(IO_stringValue(line,chunkPos,1)) == '*part' ) inPart = .true. - if ( IO_lc(IO_stringValue(line,chunkPos,1)) == '*end' .and. & - IO_lc(IO_stringValue(line,chunkPos,2)) == 'part' ) inPart = .false. - - if (inPart .or. noPart) then - select case ( IO_lc(IO_stringValue(line,chunkPos,1))) - case('*node') - if( & - IO_lc(IO_stringValue(line,chunkPos,2)) /= 'output' .and. & - IO_lc(IO_stringValue(line,chunkPos,2)) /= 'print' .and. & - IO_lc(IO_stringValue(line,chunkPos,2)) /= 'file' .and. & - IO_lc(IO_stringValue(line,chunkPos,2)) /= 'response' & - ) & - mesh_Nnodes = mesh_Nnodes + IO_countDataLines(fileUnit) - case('*element') - if( & - IO_lc(IO_stringValue(line,chunkPos,2)) /= 'output' .and. & - IO_lc(IO_stringValue(line,chunkPos,2)) /= 'matrix' .and. & - IO_lc(IO_stringValue(line,chunkPos,2)) /= 'response' & - ) then - mesh_Nelems = mesh_Nelems + IO_countDataLines(fileUnit) - endif - endselect - endif - enddo - - if (mesh_Nnodes < 2) call IO_error(error_ID=900) - if (mesh_Nelems == 0) call IO_error(error_ID=901) - -end subroutine mesh_abaqus_count_nodesAndElements - - -!-------------------------------------------------------------------------------------------------- -!> @brief count overall number of element sets in mesh and write 'mesh_NelemSets' and -!! 'mesh_maxNelemInSet' -!-------------------------------------------------------------------------------------------------- -subroutine mesh_abaqus_count_elementSets(fileUnit) - - use IO, only: IO_lc, & - IO_stringValue, & - IO_stringPos, & - IO_error - - - integer, intent(in) :: fileUnit - - integer, allocatable, dimension(:) :: chunkPos - character(len=300) :: line - integer :: myStat - logical :: inPart - - mesh_NelemSets = 0 - mesh_maxNelemInSet = mesh_Nelems ! have to be conservative, since Abaqus allows for recursive definitons - - inPart = .false. - myStat = 0 - rewind(fileUnit) - do while(myStat == 0) - read (fileUnit,'(a300)',iostat=myStat) line - chunkPos = IO_stringPos(line) - if ( IO_lc(IO_stringValue(line,chunkPos,1)) == '*part' ) inPart = .true. - if ( IO_lc(IO_stringValue(line,chunkPos,1)) == '*end' .and. & - IO_lc(IO_stringValue(line,chunkPos,2)) == 'part' ) inPart = .false. - - if ( (inPart .or. noPart) .and. IO_lc(IO_stringValue(line,chunkPos,1)) == '*elset' ) & - mesh_NelemSets = mesh_NelemSets + 1 - enddo - - if (mesh_NelemSets == 0) call IO_error(error_ID=902) - -end subroutine mesh_abaqus_count_elementSets - - -!-------------------------------------------------------------------------------------------------- -! count overall number of solid sections sets in mesh (Abaqus only) -! -! mesh_Nmaterials -!-------------------------------------------------------------------------------------------------- -subroutine mesh_abaqus_count_materials(fileUnit) - - use IO, only: IO_lc, & - IO_stringValue, & - IO_stringPos, & - IO_error - - - integer, intent(in) :: fileUnit - - integer, allocatable, dimension(:) :: chunkPos - character(len=300) :: line - integer :: myStat - logical :: inPart - - mesh_Nmaterials = 0 - - inPart = .false. - myStat = 0 - rewind(fileUnit) - do while(myStat == 0) - read (fileUnit,'(a300)',iostat=myStat) line - chunkPos = IO_stringPos(line) - if ( IO_lc(IO_stringValue(line,chunkPos,1)) == '*part' ) inPart = .true. - if ( IO_lc(IO_stringValue(line,chunkPos,1)) == '*end' .and. & - IO_lc(IO_stringValue(line,chunkPos,2)) == 'part' ) inPart = .false. - - if ( (inPart .or. noPart) .and. & - IO_lc(IO_StringValue(line,chunkPos,1)) == '*solid' .and. & - IO_lc(IO_StringValue(line,chunkPos,2)) == 'section' ) & - mesh_Nmaterials = mesh_Nmaterials + 1 - enddo - - if (mesh_Nmaterials == 0) call IO_error(error_ID=903) - -end subroutine mesh_abaqus_count_materials - - -!-------------------------------------------------------------------------------------------------- -! Build element set mapping -! -! allocate globals: mesh_nameElemSet, mesh_mapElemSet -!-------------------------------------------------------------------------------------------------- -subroutine mesh_abaqus_map_elementSets(fileUnit) - - use IO, only: IO_lc, & - IO_stringValue, & - IO_stringPos, & - IO_extractValue, & - IO_continuousIntValues, & - IO_error - - - integer, intent(in) :: fileUnit - - integer, allocatable, dimension(:) :: chunkPos - character(len=300) :: line - integer :: myStat - logical :: inPart - integer :: elemSet,i - - allocate (mesh_nameElemSet(mesh_NelemSets)); mesh_nameElemSet = '' - allocate (mesh_mapElemSet(1+mesh_maxNelemInSet,mesh_NelemSets),source=0) - - - elemSet = 0 - inPart = .false. - myStat = 0 - rewind(fileUnit) - do while(myStat == 0) - read (fileUnit,'(a300)',iostat=myStat) line - chunkPos = IO_stringPos(line) - if ( IO_lc(IO_stringValue(line,chunkPos,1)) == '*part' ) inPart = .true. - if ( IO_lc(IO_stringValue(line,chunkPos,1)) == '*end' .and. & - IO_lc(IO_stringValue(line,chunkPos,2)) == 'part' ) inPart = .false. - - if ( (inPart .or. noPart) .and. IO_lc(IO_stringValue(line,chunkPos,1)) == '*elset' ) then - elemSet = elemSet + 1 - mesh_nameElemSet(elemSet) = trim(IO_extractValue(IO_lc(IO_stringValue(line,chunkPos,2)),'elset')) - mesh_mapElemSet(:,elemSet) = IO_continuousIntValues(fileUnit,mesh_Nelems,mesh_nameElemSet,& - mesh_mapElemSet,elemSet-1) - endif - enddo - - do i = 1,elemSet - if (mesh_mapElemSet(1,i) == 0) call IO_error(error_ID=904,ext_msg=mesh_nameElemSet(i)) - enddo - -end subroutine mesh_abaqus_map_elementSets - - -!-------------------------------------------------------------------------------------------------- -! map solid section (Abaqus only) -! -! allocate globals: mesh_nameMaterial, mesh_mapMaterial -!-------------------------------------------------------------------------------------------------- -subroutine mesh_abaqus_map_materials(fileUnit) - - use IO, only: IO_lc, & - IO_stringValue, & - IO_stringPos, & - IO_extractValue, & - IO_error - - - integer, intent(in) :: fileUnit - - integer, allocatable, dimension(:) :: chunkPos - character(len=300) :: line - integer :: myStat - logical :: inPart - integer :: i,c - character(len=64) :: elemSetName,materialName - - allocate (mesh_nameMaterial(mesh_Nmaterials)); mesh_nameMaterial = '' - allocate (mesh_mapMaterial(mesh_Nmaterials)); mesh_mapMaterial = '' - - c = 0 - inPart = .false. - myStat = 0 - rewind(fileUnit) - do while(myStat == 0) - read (fileUnit,'(a300)',iostat=myStat) line - chunkPos = IO_stringPos(line) - if ( IO_lc(IO_stringValue(line,chunkPos,1)) == '*part' ) inPart = .true. - if ( IO_lc(IO_stringValue(line,chunkPos,1)) == '*end' .and. & - IO_lc(IO_stringValue(line,chunkPos,2)) == 'part' ) inPart = .false. - - if ( (inPart .or. noPart) .and. & - IO_lc(IO_StringValue(line,chunkPos,1)) == '*solid' .and. & - IO_lc(IO_StringValue(line,chunkPos,2)) == 'section' ) then - - elemSetName = '' - materialName = '' - - do i = 3,chunkPos(1) - if (IO_extractValue(IO_lc(IO_stringValue(line,chunkPos,i)),'elset') /= '') & - elemSetName = trim(IO_extractValue(IO_lc(IO_stringValue(line,chunkPos,i)),'elset')) - if (IO_extractValue(IO_lc(IO_stringValue(line,chunkPos,i)),'material') /= '') & - materialName = trim(IO_extractValue(IO_lc(IO_stringValue(line,chunkPos,i)),'material')) - enddo - - if (elemSetName /= '' .and. materialName /= '') then - c = c + 1 - mesh_nameMaterial(c) = materialName ! name of material used for this section - mesh_mapMaterial(c) = elemSetName ! mapped to respective element set - endif - endif - enddo - - if (c==0) call IO_error(error_ID=905) - do i=1,c - if (mesh_nameMaterial(i)=='' .or. mesh_mapMaterial(i)=='') call IO_error(error_ID=905) - enddo - - end subroutine mesh_abaqus_map_materials - - -!-------------------------------------------------------------------------------------------------- -!> @brief Count overall number of CP elements in mesh and stores them in 'mesh_NcpElems' -!-------------------------------------------------------------------------------------------------- -subroutine mesh_abaqus_count_cpElements(fileUnit) - - use IO, only: IO_lc, & - IO_stringValue, & - IO_stringPos, & - IO_error, & - IO_extractValue - - - integer, intent(in) :: fileUnit - - integer, allocatable, dimension(:) :: chunkPos - character(len=300) :: line - integer :: myStat - logical :: materialFound - integer :: i,k - character(len=64) ::materialName,elemSetName - - mesh_NcpElems = 0 - materialFound = .false. - myStat = 0 - rewind(fileUnit) - do while(myStat == 0) - read (fileUnit,'(a300)',iostat=myStat) line - chunkPos = IO_stringPos(line) - select case ( IO_lc(IO_stringValue(line,chunkPos,1)) ) - case('*material') - materialName = trim(IO_extractValue(IO_lc(IO_stringValue(line,chunkPos,2)),'name')) ! extract name=value - materialFound = materialName /= '' ! valid name? - case('*user') - if (IO_lc(IO_StringValue(line,chunkPos,2)) == 'material' .and. materialFound) then - do i = 1,mesh_Nmaterials ! look thru material names - if (materialName == mesh_nameMaterial(i)) then ! found one - elemSetName = mesh_mapMaterial(i) ! take corresponding elemSet - do k = 1,mesh_NelemSets ! look thru all elemSet definitions - if (elemSetName == mesh_nameElemSet(k)) & ! matched? - mesh_NcpElems = mesh_NcpElems + mesh_mapElemSet(1,k) ! add those elem count - enddo - endif - enddo - materialFound = .false. - endif - endselect - enddo - - if (mesh_NcpElems == 0) call IO_error(error_ID=906) - -end subroutine mesh_abaqus_count_cpElements - - -!-------------------------------------------------------------------------------------------------- -!> @brief Maps elements from FE ID to internal (consecutive) representation. -!! Allocates global array 'mesh_mapFEtoCPelem' -!-------------------------------------------------------------------------------------------------- -subroutine mesh_abaqus_map_elements(fileUnit) - - use math, only: math_sort - use IO, only: IO_lc, & - IO_stringValue, & - IO_stringPos, & - IO_extractValue, & - IO_error - - - integer, intent(in) :: fileUnit - - integer, allocatable, dimension(:) :: chunkPos - character(len=300) :: line - integer :: myStat - logical :: materialFound - integer ::i,j,k,cpElem - character (len=64) materialName,elemSetName ! why limited to 64? ABAQUS? - - allocate (mesh_mapFEtoCPelem(2,mesh_NcpElems), source = 0) - - cpElem = 0 - materialFound = .false. - myStat = 0 - rewind(fileUnit) - do while(myStat == 0) - read (fileUnit,'(a300)',iostat=myStat) line - chunkPos = IO_stringPos(line) - select case ( IO_lc(IO_stringValue(line,chunkPos,1)) ) - case('*material') - materialName = trim(IO_extractValue(IO_lc(IO_stringValue(line,chunkPos,2)),'name')) ! extract name=value - materialFound = materialName /= '' ! valid name? - case('*user') - if (IO_lc(IO_stringValue(line,chunkPos,2)) == 'material' .and. materialFound) then - do i = 1,mesh_Nmaterials ! look thru material names - if (materialName == mesh_nameMaterial(i)) then ! found one - elemSetName = mesh_mapMaterial(i) ! take corresponding elemSet - do k = 1,mesh_NelemSets ! look thru all elemSet definitions - if (elemSetName == mesh_nameElemSet(k)) then ! matched? - do j = 1,mesh_mapElemSet(1,k) - cpElem = cpElem + 1 - mesh_mapFEtoCPelem(1,cpElem) = mesh_mapElemSet(1+j,k) ! store FE id - mesh_mapFEtoCPelem(2,cpElem) = cpElem ! store our id - enddo - endif - enddo - endif - enddo - materialFound = .false. - endif - endselect - enddo - - call math_sort(mesh_mapFEtoCPelem,1,int(size(mesh_mapFEtoCPelem,2),pInt)) ! should be mesh_NcpElems - - if (int(size(mesh_mapFEtoCPelem),pInt) < 2) call IO_error(error_ID=907) - -end subroutine mesh_abaqus_map_elements - - -!-------------------------------------------------------------------------------------------------- -!> @brief Maps node from FE ID to internal (consecutive) representation. -!! Allocates global array 'mesh_mapFEtoCPnode' -!-------------------------------------------------------------------------------------------------- -subroutine mesh_abaqus_map_nodes(fileUnit) - - use math, only: math_sort - use IO, only: IO_lc, & - IO_stringValue, & - IO_stringPos, & - IO_countDataLines, & - IO_intValue, & - IO_error - - - integer, intent(in) :: fileUnit - - integer, allocatable, dimension(:) :: chunkPos - character(len=300) :: line - integer :: myStat - logical :: inPart - integer :: i,c,cpNode - - allocate (mesh_mapFEtoCPnode(2,mesh_Nnodes), source=0) - - cpNode = 0 - inPart = .false. - myStat = 0 - rewind(fileUnit) - do while(myStat == 0) - read (fileUnit,'(a300)',iostat=myStat) line - chunkPos = IO_stringPos(line) - if ( IO_lc(IO_stringValue(line,chunkPos,1)) == '*part' ) inPart = .true. - if ( IO_lc(IO_stringValue(line,chunkPos,1)) == '*end' .and. & - IO_lc(IO_stringValue(line,chunkPos,2)) == 'part' ) inPart = .false. - - if( (inPart .or. noPart) .and. & - IO_lc(IO_stringValue(line,chunkPos,1)) == '*node' .and. & - ( IO_lc(IO_stringValue(line,chunkPos,2)) /= 'output' .and. & - IO_lc(IO_stringValue(line,chunkPos,2)) /= 'print' .and. & - IO_lc(IO_stringValue(line,chunkPos,2)) /= 'file' .and. & - IO_lc(IO_stringValue(line,chunkPos,2)) /= 'response' ) & - ) then - c = IO_countDataLines(fileUnit) - do i = 1,c - backspace(fileUnit) - enddo - do i = 1,c - read (fileUnit,'(a300)') line - chunkPos = IO_stringPos(line) - cpNode = cpNode + 1 - mesh_mapFEtoCPnode(1,cpNode) = IO_intValue(line,chunkPos,1) - mesh_mapFEtoCPnode(2,cpNode) = cpNode - enddo - endif - enddo - - call math_sort(mesh_mapFEtoCPnode,1,int(size(mesh_mapFEtoCPnode,2),pInt)) - - if (int(size(mesh_mapFEtoCPnode),pInt) == 0) call IO_error(error_ID=908) - -end subroutine mesh_abaqus_map_nodes - - -!-------------------------------------------------------------------------------------------------- -!> @brief store x,y,z coordinates of all nodes in mesh. -!! Allocates global arrays 'mesh_node0' and 'mesh_node' -!-------------------------------------------------------------------------------------------------- -subroutine mesh_abaqus_build_nodes(fileUnit) - use IO, only: & - IO_lc, & - IO_stringValue, & - IO_floatValue, & - IO_stringPos, & - IO_error, & - IO_countDataLines, & - IO_intValue - - - integer, intent(in) :: fileUnit - - integer, allocatable, dimension(:) :: chunkPos - character(len=300) :: line - integer :: myStat - logical :: inPart - integer :: i,j,m,c - - allocate ( mesh_node0 (3,mesh_Nnodes), source=0.0_pReal) - allocate ( mesh_node (3,mesh_Nnodes), source=0.0_pReal) - - inPart = .false. - myStat = 0 - rewind(fileUnit) - do while(myStat == 0) - read (fileUnit,'(a300)',iostat=myStat) line - chunkPos = IO_stringPos(line) - if ( IO_lc(IO_stringValue(line,chunkPos,1)) == '*part' ) inPart = .true. - if ( IO_lc(IO_stringValue(line,chunkPos,1)) == '*end' .and. & - IO_lc(IO_stringValue(line,chunkPos,2)) == 'part' ) inPart = .false. - - if( (inPart .or. noPart) .and. & - IO_lc(IO_stringValue(line,chunkPos,1)) == '*node' .and. & - ( IO_lc(IO_stringValue(line,chunkPos,2)) /= 'output' .and. & - IO_lc(IO_stringValue(line,chunkPos,2)) /= 'print' .and. & - IO_lc(IO_stringValue(line,chunkPos,2)) /= 'file' .and. & - IO_lc(IO_stringValue(line,chunkPos,2)) /= 'response' ) & - ) then - c = IO_countDataLines(fileUnit) ! how many nodes are defined here? - do i = 1,c - backspace(fileUnit) ! rewind to first entry - enddo - do i = 1,c - read (fileUnit,'(a300)') line - chunkPos = IO_stringPos(line) - m = mesh_FEasCP('node',IO_intValue(line,chunkPos,1)) - do j=1, 3 - mesh_node0(j,m) = mesh_unitlength * IO_floatValue(line,chunkPos,j+1) - enddo - enddo - endif - enddo - - if (int(size(mesh_node0,2),pInt) /= mesh_Nnodes) call IO_error(error_ID=909) - mesh_node = mesh_node0 - -end subroutine mesh_abaqus_build_nodes - - -!-------------------------------------------------------------------------------------------------- -!> @brief Gets maximum count of nodes, IPs, IP neighbors, and subNodes among cpElements. -!! Sets global values 'mesh_maxNnodes', 'mesh_maxNips', 'mesh_maxNipNeighbors', -!! and 'mesh_maxNcellnodes' -!-------------------------------------------------------------------------------------------------- -subroutine mesh_abaqus_count_cpSizes(fileUnit) - - use IO, only: IO_lc, & - IO_stringValue, & - IO_stringPos, & - IO_extractValue ,& - IO_error, & - IO_countDataLines, & - IO_intValue - - - integer, intent(in) :: fileUnit - - integer, allocatable, dimension(:) :: chunkPos - character(len=300) :: line - integer :: myStat - logical :: inPart - integer :: i,c,t,g - - mesh_maxNnodes = 0 - mesh_maxNips = 0 - mesh_maxNipNeighbors = 0 - mesh_maxNcellnodes = 0 - - - inPart = .false. - myStat = 0 - rewind(fileUnit) - do while(myStat == 0) - read (fileUnit,'(a300)',iostat=myStat) line - chunkPos = IO_stringPos(line) - if ( IO_lc(IO_stringValue(line,chunkPos,1)) == '*part' ) inPart = .true. - if ( IO_lc(IO_stringValue(line,chunkPos,1)) == '*end' .and. & - IO_lc(IO_stringValue(line,chunkPos,2)) == 'part' ) inPart = .false. - - if( (inPart .or. noPart) .and. & - IO_lc(IO_stringValue(line,chunkPos,1)) == '*element' .and. & - ( IO_lc(IO_stringValue(line,chunkPos,2)) /= 'output' .and. & - IO_lc(IO_stringValue(line,chunkPos,2)) /= 'matrix' .and. & - IO_lc(IO_stringValue(line,chunkPos,2)) /= 'response' ) & - ) then - t = FE_mapElemtype(IO_extractValue(IO_lc(IO_stringValue(line,chunkPos,2)),'type')) ! remember elem type - g = FE_geomtype(t) - c = FE_celltype(g) - mesh_maxNnodes = max(mesh_maxNnodes,FE_Nnodes(t)) - mesh_maxNips = max(mesh_maxNips,FE_Nips(g)) - mesh_maxNipNeighbors = max(mesh_maxNipNeighbors,FE_NipNeighbors(c)) - mesh_maxNcellnodes = max(mesh_maxNcellnodes,FE_Ncellnodes(g)) - endif - enddo - -end subroutine mesh_abaqus_count_cpSizes - - -!-------------------------------------------------------------------------------------------------- -!> @brief Store FEid, type, mat, tex, and node list per elemen. -!! Allocates global array 'mesh_element' -!-------------------------------------------------------------------------------------------------- -subroutine mesh_abaqus_build_elements(fileUnit) - - use IO, only: IO_lc, & - IO_stringValue, & - IO_stringPos, & - IO_intValue, & - IO_extractValue, & - IO_floatValue, & - IO_countDataLines, & - IO_error - - - integer, intent(in) :: fileUnit - - integer, allocatable, dimension(:) :: chunkPos - character(len=300) :: line - integer :: myStat - logical :: inPart, materialFound - integer :: i,j,k,c,e,t,homog,micro, nNodesAlreadyRead - character (len=64) :: materialName,elemSetName - - allocate(mesh_element (4+mesh_maxNnodes,mesh_NcpElems), source=0) - mesh_elemType = -1 - - inPart = .false. - myStat = 0 - rewind(fileUnit) - do while(myStat == 0) - read (fileUnit,'(a300)',iostat=myStat) line - chunkPos = IO_stringPos(line) - if ( IO_lc(IO_stringValue(line,chunkPos,1)) == '*part' ) inPart = .true. - if ( IO_lc(IO_stringValue(line,chunkPos,1)) == '*end' .and. & - IO_lc(IO_stringValue(line,chunkPos,2)) == 'part' ) inPart = .false. - - if( (inPart .or. noPart) .and. & - IO_lc(IO_stringValue(line,chunkPos,1)) == '*element' .and. & - ( IO_lc(IO_stringValue(line,chunkPos,2)) /= 'output' .and. & - IO_lc(IO_stringValue(line,chunkPos,2)) /= 'matrix' .and. & - IO_lc(IO_stringValue(line,chunkPos,2)) /= 'response' ) & - ) then - t = FE_mapElemtype(IO_extractValue(IO_lc(IO_stringValue(line,chunkPos,2)),'type')) ! remember elem type - c = IO_countDataLines(fileUnit) - do i = 1,c - backspace(fileUnit) - enddo - do i = 1,c - read (fileUnit,'(a300)') line - chunkPos = IO_stringPos(line) ! limit to 64 nodes max - e = mesh_FEasCP('elem',IO_intValue(line,chunkPos,1)) - if (e /= 0) then ! disregard non CP elems - mesh_element(1,e) = -1 ! DEPRECATED - if (mesh_elemType /= t .and. mesh_elemType /= -1) & - call IO_error(191,el=t,ip=mesh_elemType) - mesh_elemType = t - mesh_element(2,e) = t ! elem type - nNodesAlreadyRead = 0 - do j = 1,chunkPos(1)-1 - mesh_element(4+j,e) = mesh_FEasCP('node',IO_intValue(line,chunkPos,1+j)) ! put CP ids of nodes to position 5: - enddo - nNodesAlreadyRead = chunkPos(1) - 1 - do while(nNodesAlreadyRead < FE_Nnodes(t)) ! read on if not all nodes in one line - read (fileUnit,'(a300)') line - chunkPos = IO_stringPos(line) - do j = 1,chunkPos(1) - mesh_element(4+nNodesAlreadyRead+j,e) & - = mesh_FEasCP('node',IO_IntValue(line,chunkPos,j)) ! CP ids of nodes - enddo - nNodesAlreadyRead = nNodesAlreadyRead + chunkPos(1) - enddo - endif - enddo - endif - enddo - - - rewind(fileUnit) ! just in case "*material" definitions apear before "*element" - - materialFound = .false. - myStat = 0 - rewind(fileUnit) - do while(myStat == 0) - read (fileUnit,'(a300)',iostat=myStat) line - chunkPos = IO_stringPos(line) - select case ( IO_lc(IO_StringValue(line,chunkPos,1))) - case('*material') - materialName = trim(IO_extractValue(IO_lc(IO_StringValue(line,chunkPos,2)),'name')) ! extract name=value - materialFound = materialName /= '' ! valid name? - case('*user') - if ( IO_lc(IO_StringValue(line,chunkPos,2)) == 'material' .and. & - materialFound ) then - read (fileUnit,'(a300)') line ! read homogenization and microstructure - chunkPos = IO_stringPos(line) - homog = nint(IO_floatValue(line,chunkPos,1),pInt) - micro = nint(IO_floatValue(line,chunkPos,2),pInt) - do i = 1,mesh_Nmaterials ! look thru material names - if (materialName == mesh_nameMaterial(i)) then ! found one - elemSetName = mesh_mapMaterial(i) ! take corresponding elemSet - do k = 1,mesh_NelemSets ! look thru all elemSet definitions - if (elemSetName == mesh_nameElemSet(k)) then ! matched? - do j = 1,mesh_mapElemSet(1,k) - e = mesh_FEasCP('elem',mesh_mapElemSet(1+j,k)) - mesh_element(3,e) = homog ! store homogenization - mesh_element(4,e) = micro ! store microstructure - mesh_maxValStateVar(1) = max(mesh_maxValStateVar(1),homog) - mesh_maxValStateVar(2) = max(mesh_maxValStateVar(2),micro) - enddo - endif - enddo - endif - enddo - materialFound = .false. - endif - endselect - enddo - -end subroutine mesh_abaqus_build_elements - - -!-------------------------------------------------------------------------------------------------- -!> @brief get any additional damask options from input file, sets mesh_periodicSurface -!-------------------------------------------------------------------------------------------------- -subroutine mesh_get_damaskOptions(periodic_surface,fileUnit) - -use IO, only: & - IO_lc, & - IO_stringValue, & - IO_stringPos - - - integer, intent(in) :: fileUnit - - integer, allocatable, dimension(:) :: chunkPos - character(len=300) :: line - integer :: myStat - integer :: chunk, Nchunks - character(len=300) :: v - logical, dimension(3) :: periodic_surface - - - periodic_surface = .false. - myStat = 0 - rewind(fileUnit) - do while(myStat == 0) - read (fileUnit,'(a300)',iostat=myStat) line - chunkPos = IO_stringPos(line) - Nchunks = chunkPos(1) - if (IO_lc(IO_stringValue(line,chunkPos,1)) == '**damask' .and. Nchunks > 1) then ! found keyword for damask option and there is at least one more chunk to read - select case(IO_lc(IO_stringValue(line,chunkPos,2))) - case('periodic') ! damask Option that allows to specify periodic fluxes - do chunk = 3,Nchunks ! loop through chunks (skipping the keyword) - v = IO_lc(IO_stringValue(line,chunkPos,chunk)) ! chunk matches keyvalues x,y, or z? - mesh_periodicSurface(1) = mesh_periodicSurface(1) .or. v == 'x' - mesh_periodicSurface(2) = mesh_periodicSurface(2) .or. v == 'y' - mesh_periodicSurface(3) = mesh_periodicSurface(3) .or. v == 'z' - enddo - endselect - endif - enddo - -end subroutine mesh_get_damaskOptions - - -!-------------------------------------------------------------------------------------------------- -!> @brief Split CP elements into cells. -!> @details Build a mapping between cells and the corresponding cell nodes ('mesh_cell'). -!> Cell nodes that are also matching nodes are unique in the list of cell nodes, -!> all others (currently) might be stored more than once. -!> Also allocates the 'mesh_node' array. -!-------------------------------------------------------------------------------------------------- -subroutine mesh_build_cellconnectivity - - - integer, dimension(:), allocatable :: & - matchingNode2cellnode - integer, dimension(:,:), allocatable :: & - cellnodeParent - integer, dimension(mesh_maxNcellnodes) :: & - localCellnode2globalCellnode - integer :: & - e,t,g,c,n,i, & - matchingNodeID, & - localCellnodeID - - allocate(mesh_cell(FE_maxNcellnodesPerCell,mesh_maxNips,mesh_NcpElems), source=0) - allocate(matchingNode2cellnode(mesh_Nnodes), source=0) - allocate(cellnodeParent(2,mesh_maxNcellnodes*mesh_NcpElems), source=0) - -!-------------------------------------------------------------------------------------------------- -! Count cell nodes (including duplicates) and generate cell connectivity list - mesh_Ncellnodes = 0 - mesh_Ncells = 0 - do e = 1,mesh_NcpElems ! loop over cpElems - t = mesh_element(2,e) ! get element type - g = FE_geomtype(t) ! get geometry type - c = FE_celltype(g) ! get cell type - localCellnode2globalCellnode = 0 - mesh_Ncells = mesh_Ncells + FE_Nips(g) - do i = 1,FE_Nips(g) ! loop over ips=cells in this element - do n = 1,FE_NcellnodesPerCell(c) ! loop over cell nodes in this cell - localCellnodeID = FE_cell(n,i,g) - if (localCellnodeID <= FE_NmatchingNodes(g)) then ! this cell node is a matching node - matchingNodeID = mesh_element(4+localCellnodeID,e) - if (matchingNode2cellnode(matchingNodeID) == 0) then ! if this matching node does not yet exist in the glbal cell node list ... - mesh_Ncellnodes = mesh_Ncellnodes + 1 ! ... count it as cell node ... - matchingNode2cellnode(matchingNodeID) = mesh_Ncellnodes ! ... and remember its global ID - cellnodeParent(1,mesh_Ncellnodes) = e ! ... and where it belongs to - cellnodeParent(2,mesh_Ncellnodes) = localCellnodeID - endif - mesh_cell(n,i,e) = matchingNode2cellnode(matchingNodeID) - else ! this cell node is no matching node - if (localCellnode2globalCellnode(localCellnodeID) == 0) then ! if this local cell node does not yet exist in the global cell node list ... - mesh_Ncellnodes = mesh_Ncellnodes + 1 ! ... count it as cell node ... - localCellnode2globalCellnode(localCellnodeID) = mesh_Ncellnodes ! ... and remember its global ID ... - cellnodeParent(1,mesh_Ncellnodes) = e ! ... and it belongs to - cellnodeParent(2,mesh_Ncellnodes) = localCellnodeID - endif - mesh_cell(n,i,e) = localCellnode2globalCellnode(localCellnodeID) - endif - enddo - enddo - enddo - - allocate(mesh_cellnodeParent(2,mesh_Ncellnodes)) - allocate(mesh_cellnode(3,mesh_Ncellnodes)) - forall(n = 1:mesh_Ncellnodes) - mesh_cellnodeParent(1,n) = cellnodeParent(1,n) - mesh_cellnodeParent(2,n) = cellnodeParent(2,n) - endforall - -end subroutine mesh_build_cellconnectivity - - -!-------------------------------------------------------------------------------------------------- -!> @brief Calculate position of cellnodes from the given position of nodes -!> Build list of cellnodes' coordinates. -!> Cellnode coordinates are calculated from a weighted sum of node coordinates. -!-------------------------------------------------------------------------------------------------- -function mesh_build_cellnodes(nodes,Ncellnodes) - - - integer, intent(in) :: Ncellnodes !< requested number of cellnodes - real(pReal), dimension(3,mesh_Nnodes), intent(in) :: nodes - real(pReal), dimension(3,Ncellnodes) :: mesh_build_cellnodes - - integer :: & - e,t,n,m, & - localCellnodeID - real(pReal), dimension(3) :: & - myCoords - - mesh_build_cellnodes = 0.0_pReal -!$OMP PARALLEL DO PRIVATE(e,localCellnodeID,t,myCoords) - do n = 1,Ncellnodes ! loop over cell nodes - e = mesh_cellnodeParent(1,n) - localCellnodeID = mesh_cellnodeParent(2,n) - t = mesh_element(2,e) ! get element type - myCoords = 0.0_pReal - do m = 1,FE_Nnodes(t) - myCoords = myCoords + nodes(1:3,mesh_element(4+m,e)) & - * FE_cellnodeParentnodeWeights(m,localCellnodeID,t) - enddo - mesh_build_cellnodes(1:3,n) = myCoords / sum(FE_cellnodeParentnodeWeights(:,localCellnodeID,t)) - enddo -!$OMP END PARALLEL DO - -end function mesh_build_cellnodes - - -!-------------------------------------------------------------------------------------------------- -!> @brief Calculates IP volume. Allocates global array 'mesh_ipVolume' -!> @details The IP volume is calculated differently depending on the cell type. -!> 2D cells assume an element depth of one in order to calculate the volume. -!> For the hexahedral cell we subdivide the cell into subvolumes of pyramidal -!> shape with a cell face as basis and the central ip at the tip. This subvolume is -!> calculated as an average of four tetrahedals with three corners on the cell face -!> and one corner at the central ip. -!-------------------------------------------------------------------------------------------------- -subroutine mesh_build_ipVolumes - use math, only: & - math_volTetrahedron, & - math_areaTriangle - - - integer :: e,t,g,c,i,m,f,n - real(pReal), dimension(FE_maxNcellnodesPerCellface,FE_maxNcellfaces) :: subvolume - - allocate(mesh_ipVolume(mesh_maxNips,mesh_NcpElems),source=0.0_pReal) - - !$OMP PARALLEL DO PRIVATE(t,g,c,m,subvolume) - do e = 1,mesh_NcpElems ! loop over cpElems - t = mesh_element(2,e) ! get element type - g = FE_geomtype(t) ! get geometry type - c = FE_celltype(g) ! get cell type - select case (c) - - case (1) ! 2D 3node - forall (i = 1:FE_Nips(g)) & ! loop over ips=cells in this element - mesh_ipVolume(i,e) = math_areaTriangle(mesh_cellnode(1:3,mesh_cell(1,i,e)), & - mesh_cellnode(1:3,mesh_cell(2,i,e)), & - mesh_cellnode(1:3,mesh_cell(3,i,e))) - - case (2) ! 2D 4node - forall (i = 1:FE_Nips(g)) & ! loop over ips=cells in this element - mesh_ipVolume(i,e) = math_areaTriangle(mesh_cellnode(1:3,mesh_cell(1,i,e)), & ! here we assume a planar shape, so division in two triangles suffices - mesh_cellnode(1:3,mesh_cell(2,i,e)), & - mesh_cellnode(1:3,mesh_cell(3,i,e))) & - + math_areaTriangle(mesh_cellnode(1:3,mesh_cell(3,i,e)), & - mesh_cellnode(1:3,mesh_cell(4,i,e)), & - mesh_cellnode(1:3,mesh_cell(1,i,e))) - - case (3) ! 3D 4node - forall (i = 1:FE_Nips(g)) & ! loop over ips=cells in this element - mesh_ipVolume(i,e) = math_volTetrahedron(mesh_cellnode(1:3,mesh_cell(1,i,e)), & - mesh_cellnode(1:3,mesh_cell(2,i,e)), & - mesh_cellnode(1:3,mesh_cell(3,i,e)), & - mesh_cellnode(1:3,mesh_cell(4,i,e))) - - case (4) ! 3D 8node - m = FE_NcellnodesPerCellface(c) - do i = 1,FE_Nips(g) ! loop over ips=cells in this element - subvolume = 0.0_pReal - forall(f = 1:FE_NipNeighbors(c), n = 1:FE_NcellnodesPerCellface(c)) & - subvolume(n,f) = math_volTetrahedron(& - mesh_cellnode(1:3,mesh_cell(FE_cellface( n ,f,c),i,e)), & - mesh_cellnode(1:3,mesh_cell(FE_cellface(1+mod(n ,m),f,c),i,e)), & - mesh_cellnode(1:3,mesh_cell(FE_cellface(1+mod(n+1,m),f,c),i,e)), & - mesh_ipCoordinates(1:3,i,e)) - mesh_ipVolume(i,e) = 0.5_pReal * sum(subvolume) ! each subvolume is based on four tetrahedrons, altough the face consists of only two triangles -> averaging factor two - enddo - - end select - enddo - !$OMP END PARALLEL DO - -end subroutine mesh_build_ipVolumes - - -!-------------------------------------------------------------------------------------------------- -!> @brief Calculates IP Coordinates. Allocates global array 'mesh_ipCoordinates' -! Called by all solvers in mesh_init in order to initialize the ip coordinates. -! Later on the current ip coordinates are directly prvided by the spectral solver and by Abaqus, -! so no need to use this subroutine anymore; Marc however only provides nodal displacements, -! so in this case the ip coordinates are always calculated on the basis of this subroutine. -! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -! FOR THE MOMENT THIS SUBROUTINE ACTUALLY CALCULATES THE CELL CENTER AND NOT THE IP COORDINATES, -! AS THE IP IS NOT (ALWAYS) LOCATED IN THE CENTER OF THE IP VOLUME. -! HAS TO BE CHANGED IN A LATER VERSION. -! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -!-------------------------------------------------------------------------------------------------- -subroutine mesh_build_ipCoordinates - - - integer :: e,t,g,c,i,n - real(pReal), dimension(3) :: myCoords - - if (.not. allocated(mesh_ipCoordinates)) & - allocate(mesh_ipCoordinates(3,mesh_maxNips,mesh_NcpElems),source=0.0_pReal) - - !$OMP PARALLEL DO PRIVATE(t,g,c,myCoords) - do e = 1,mesh_NcpElems ! loop over cpElems - t = mesh_element(2,e) ! get element type - g = FE_geomtype(t) ! get geometry type - c = FE_celltype(g) ! get cell type - do i = 1,FE_Nips(g) ! loop over ips=cells in this element - myCoords = 0.0_pReal - do n = 1,FE_NcellnodesPerCell(c) ! loop over cell nodes in this cell - myCoords = myCoords + mesh_cellnode(1:3,mesh_cell(n,i,e)) - enddo - mesh_ipCoordinates(1:3,i,e) = myCoords / real(FE_NcellnodesPerCell(c),pReal) - enddo - enddo - !$OMP END PARALLEL DO - -end subroutine mesh_build_ipCoordinates - - -!-------------------------------------------------------------------------------------------------- -!> @brief Calculates cell center coordinates. -!-------------------------------------------------------------------------------------------------- -pure function mesh_cellCenterCoordinates(ip,el) - - - integer, intent(in) :: el, & !< element number - ip !< integration point number - real(pReal), dimension(3) :: mesh_cellCenterCoordinates !< x,y,z coordinates of the cell center of the requested IP cell - integer :: t,g,c,n - - t = mesh_element(2,el) ! get element type - g = FE_geomtype(t) ! get geometry type - c = FE_celltype(g) ! get cell type - mesh_cellCenterCoordinates = 0.0_pReal - do n = 1,FE_NcellnodesPerCell(c) ! loop over cell nodes in this cell - mesh_cellCenterCoordinates = mesh_cellCenterCoordinates + mesh_cellnode(1:3,mesh_cell(n,ip,el)) - enddo - mesh_cellCenterCoordinates = mesh_cellCenterCoordinates / real(FE_NcellnodesPerCell(c),pReal) - - end function mesh_cellCenterCoordinates - - - - - - -!-------------------------------------------------------------------------------------------------- -!> @brief calculation of IP interface areas, allocate globals '_ipArea', and '_ipAreaNormal' -!-------------------------------------------------------------------------------------------------- -subroutine mesh_build_ipAreas - use math, only: & - math_cross - - - integer :: e,t,g,c,i,f,n,m - real(pReal), dimension (3,FE_maxNcellnodesPerCellface) :: nodePos, normals - real(pReal), dimension(3) :: normal - - allocate(mesh_ipArea(mesh_maxNipNeighbors,mesh_maxNips,mesh_NcpElems), source=0.0_pReal) - allocate(mesh_ipAreaNormal(3,mesh_maxNipNeighbors,mesh_maxNips,mesh_NcpElems), source=0.0_pReal) - - !$OMP PARALLEL DO PRIVATE(t,g,c,nodePos,normal,normals) - do e = 1,mesh_NcpElems ! loop over cpElems - t = mesh_element(2,e) ! get element type - g = FE_geomtype(t) ! get geometry type - c = FE_celltype(g) ! get cell type - select case (c) - - case (1,2) ! 2D 3 or 4 node - do i = 1,FE_Nips(g) ! loop over ips=cells in this element - do f = 1,FE_NipNeighbors(c) ! loop over cell faces - forall(n = 1:FE_NcellnodesPerCellface(c)) & - nodePos(1:3,n) = mesh_cellnode(1:3,mesh_cell(FE_cellface(n,f,c),i,e)) - normal(1) = nodePos(2,2) - nodePos(2,1) ! x_normal = y_connectingVector - normal(2) = -(nodePos(1,2) - nodePos(1,1)) ! y_normal = -x_connectingVector - normal(3) = 0.0_pReal - mesh_ipArea(f,i,e) = norm2(normal) - mesh_ipAreaNormal(1:3,f,i,e) = normal / norm2(normal) ! ensure unit length of area normal - enddo - enddo - - case (3) ! 3D 4node - do i = 1,FE_Nips(g) ! loop over ips=cells in this element - do f = 1,FE_NipNeighbors(c) ! loop over cell faces - forall(n = 1:FE_NcellnodesPerCellface(c)) & - nodePos(1:3,n) = mesh_cellnode(1:3,mesh_cell(FE_cellface(n,f,c),i,e)) - normal = math_cross(nodePos(1:3,2) - nodePos(1:3,1), & - nodePos(1:3,3) - nodePos(1:3,1)) - mesh_ipArea(f,i,e) = norm2(normal) - mesh_ipAreaNormal(1:3,f,i,e) = normal / norm2(normal) ! ensure unit length of area normal - enddo - enddo - - case (4) ! 3D 8node - ! for this cell type we get the normal of the quadrilateral face as an average of - ! four normals of triangular subfaces; since the face consists only of two triangles, - ! the sum has to be divided by two; this whole prcedure tries to compensate for - ! probable non-planar cell surfaces - m = FE_NcellnodesPerCellface(c) - do i = 1,FE_Nips(g) ! loop over ips=cells in this element - do f = 1,FE_NipNeighbors(c) ! loop over cell faces - forall(n = 1:FE_NcellnodesPerCellface(c)) & - nodePos(1:3,n) = mesh_cellnode(1:3,mesh_cell(FE_cellface(n,f,c),i,e)) - forall(n = 1:FE_NcellnodesPerCellface(c)) & - normals(1:3,n) = 0.5_pReal & - * math_cross(nodePos(1:3,1+mod(n ,m)) - nodePos(1:3,n), & - nodePos(1:3,1+mod(n+1,m)) - nodePos(1:3,n)) - normal = 0.5_pReal * sum(normals,2) - mesh_ipArea(f,i,e) = norm2(normal) - mesh_ipAreaNormal(1:3,f,i,e) = normal / norm2(normal) - enddo - enddo - - end select - enddo - !$OMP END PARALLEL DO - -end subroutine mesh_build_ipAreas - - -!-------------------------------------------------------------------------------------------------- -!> @brief assignment of twin nodes for each cp node, allocate globals '_nodeTwins' -!-------------------------------------------------------------------------------------------------- -subroutine mesh_build_nodeTwins - - - integer dir, & ! direction of periodicity - node, & - minimumNode, & - maximumNode, & - n1, & - n2 - integer, dimension(mesh_Nnodes+1) :: minimumNodes, maximumNodes ! list of surface nodes (minimum and maximum coordinate value) with first entry giving the number of nodes - real(pReal) minCoord, maxCoord, & ! extreme positions in one dimension - tolerance ! tolerance below which positions are assumed identical - real(pReal), dimension(3) :: distance ! distance between two nodes in all three coordinates - logical, dimension(mesh_Nnodes) :: unpaired - - allocate(mesh_nodeTwins(3,mesh_Nnodes)) - mesh_nodeTwins = 0 - - tolerance = 0.001_pReal * minval(mesh_ipVolume) ** 0.333_pReal - - do dir = 1,3 ! check periodicity in directions of x,y,z - if (mesh_periodicSurface(dir)) then ! only if periodicity is requested - - - !*** find out which nodes sit on the surface - !*** and have a minimum or maximum position in this dimension - - minimumNodes = 0 - maximumNodes = 0 - minCoord = minval(mesh_node0(dir,:)) - maxCoord = maxval(mesh_node0(dir,:)) - do node = 1,mesh_Nnodes ! loop through all nodes and find surface nodes - if (abs(mesh_node0(dir,node) - minCoord) <= tolerance) then - minimumNodes(1) = minimumNodes(1) + 1 - minimumNodes(minimumNodes(1)+1) = node - elseif (abs(mesh_node0(dir,node) - maxCoord) <= tolerance) then - maximumNodes(1) = maximumNodes(1) + 1 - maximumNodes(maximumNodes(1)+1) = node - endif - enddo - - - !*** find the corresponding node on the other side with the same position in this dimension - - unpaired = .true. - do n1 = 1,minimumNodes(1) - minimumNode = minimumNodes(n1+1) - if (unpaired(minimumNode)) then - do n2 = 1,maximumNodes(1) - maximumNode = maximumNodes(n2+1) - distance = abs(mesh_node0(:,minimumNode) - mesh_node0(:,maximumNode)) - if (sum(distance) - distance(dir) <= tolerance) then ! minimum possible distance (within tolerance) - mesh_nodeTwins(dir,minimumNode) = maximumNode - mesh_nodeTwins(dir,maximumNode) = minimumNode - unpaired(maximumNode) = .false. ! remember this node, we don't have to look for his partner again - exit - endif - enddo - endif - enddo - - endif - enddo - -end subroutine mesh_build_nodeTwins - - -!-------------------------------------------------------------------------------------------------- -!> @brief get maximum count of shared elements among cpElements and build list of elements shared -!! by each node in mesh. Allocate globals '_maxNsharedElems' and '_sharedElem' -!-------------------------------------------------------------------------------------------------- -subroutine mesh_build_sharedElems - - - integer(pint) e, & ! element index - g, & ! element type - node, & ! CP node index - n, & ! node index per element - myDim, & ! dimension index - nodeTwin ! node twin in the specified dimension - integer, dimension (mesh_Nnodes) :: node_count - integer, dimension(:), allocatable :: node_seen - - allocate(node_seen(maxval(FE_NmatchingNodes))) - - node_count = 0 - - do e = 1,mesh_NcpElems - g = FE_geomtype(mesh_element(2,e)) ! get elemGeomType - node_seen = 0 ! reset node duplicates - do n = 1,FE_NmatchingNodes(g) ! check each node of element - node = mesh_element(4+n,e) - if (all(node_seen /= node)) then - node_count(node) = node_count(node) + 1 ! if FE node not yet encountered -> count it - do myDim = 1,3 ! check in each dimension... - nodeTwin = mesh_nodeTwins(myDim,node) - if (nodeTwin > 0) & ! if I am a twin of some node... - node_count(nodeTwin) = node_count(nodeTwin) + 1 ! -> count me again for the twin node - enddo - endif - node_seen(n) = node ! remember this node to be counted already - enddo - enddo - - mesh_maxNsharedElems = int(maxval(node_count),pInt) ! most shared node - - allocate(mesh_sharedElem(1+mesh_maxNsharedElems,mesh_Nnodes),source=0) - - do e = 1,mesh_NcpElems - g = FE_geomtype(mesh_element(2,e)) ! get elemGeomType - node_seen = 0 - do n = 1,FE_NmatchingNodes(g) - node = mesh_element(4+n,e) - if (all(node_seen /= node)) then - mesh_sharedElem(1,node) = mesh_sharedElem(1,node) + 1 ! count for each node the connected elements - mesh_sharedElem(mesh_sharedElem(1,node)+1,node) = e ! store the respective element id - do myDim = 1,3 ! check in each dimension... - nodeTwin = mesh_nodeTwins(myDim,node) - if (nodeTwin > 0) then ! if i am a twin of some node... - mesh_sharedElem(1,nodeTwin) = mesh_sharedElem(1,nodeTwin) + 1 ! ...count me again for the twin - mesh_sharedElem(mesh_sharedElem(1,nodeTwin)+1,nodeTwin) = e ! store the respective element id - endif - enddo - endif - node_seen(n) = node - enddo - enddo - -end subroutine mesh_build_sharedElems - - -!-------------------------------------------------------------------------------------------------- -!> @brief build up of IP neighborhood, allocate globals '_ipNeighborhood' -!-------------------------------------------------------------------------------------------------- -subroutine mesh_build_ipNeighborhood - - integer :: myElem, & ! my CP element index - myIP, & - myType, & ! my element type - myFace, & - neighbor, & ! neighor index - neighboringIPkey, & ! positive integer indicating the neighboring IP (for intra-element) and negative integer indicating the face towards neighbor (for neighboring element) - candidateIP, & - neighboringType, & ! element type of neighbor - NlinkedNodes, & ! number of linked nodes - twin_of_linkedNode, & ! node twin of a specific linkedNode - NmatchingNodes, & ! number of matching nodes - dir, & ! direction of periodicity - matchingElem, & ! CP elem number of matching element - matchingFace, & ! face ID of matching element - a, anchor, & - neighboringIP, & - neighboringElem, & - pointingToMe - integer, dimension(FE_maxmaxNnodesAtIP) :: & - linkedNodes = 0, & - matchingNodes - logical checkTwins - - allocate(mesh_ipNeighborhood(3,mesh_maxNipNeighbors,mesh_maxNips,mesh_NcpElems)) - mesh_ipNeighborhood = 0 - - - do myElem = 1,mesh_NcpElems ! loop over cpElems - myType = FE_geomtype(mesh_element(2,myElem)) ! get elemGeomType - do myIP = 1,FE_Nips(myType) ! loop over IPs of elem - - do neighbor = 1,FE_NipNeighbors(FE_celltype(myType)) ! loop over neighbors of IP - neighboringIPkey = FE_ipNeighbor(neighbor,myIP,myType) - - !*** if the key is positive, the neighbor is inside the element - !*** that means, we have already found our neighboring IP - - if (neighboringIPkey > 0) then - mesh_ipNeighborhood(1,neighbor,myIP,myElem) = myElem - mesh_ipNeighborhood(2,neighbor,myIP,myElem) = neighboringIPkey - - - !*** if the key is negative, the neighbor resides in a neighboring element - !*** that means, we have to look through the face indicated by the key and see which element is behind that face - - elseif (neighboringIPkey < 0) then ! neighboring element's IP - myFace = -neighboringIPkey - call mesh_faceMatch(myElem, myFace, matchingElem, matchingFace) ! get face and CP elem id of face match - if (matchingElem > 0) then ! found match? - neighboringType = FE_geomtype(mesh_element(2,matchingElem)) - - !*** trivial solution if neighbor has only one IP - - if (FE_Nips(neighboringType) == 1) then - mesh_ipNeighborhood(1,neighbor,myIP,myElem) = matchingElem - mesh_ipNeighborhood(2,neighbor,myIP,myElem) = 1 - cycle - endif - - !*** find those nodes which build the link to the neighbor - - NlinkedNodes = 0 - linkedNodes = 0 - do a = 1,FE_maxNnodesAtIP(myType) ! figure my anchor nodes on connecting face - anchor = FE_nodesAtIP(a,myIP,myType) - if (anchor /= 0) then ! valid anchor node - if (any(FE_face(:,myFace,myType) == anchor)) then ! ip anchor sits on face? - NlinkedNodes = NlinkedNodes + 1 - linkedNodes(NlinkedNodes) = mesh_element(4+anchor,myElem) ! CP id of anchor node - else ! something went wrong with the linkage, since not all anchors sit on my face - NlinkedNodes = 0 - linkedNodes = 0 - exit - endif - endif - enddo - - !*** loop through the ips of my neighbor - !*** and try to find an ip with matching nodes - !*** also try to match with node twins - - checkCandidateIP: do candidateIP = 1,FE_Nips(neighboringType) - NmatchingNodes = 0 - matchingNodes = 0 - do a = 1,FE_maxNnodesAtIP(neighboringType) ! check each anchor node of that ip - anchor = FE_nodesAtIP(a,candidateIP,neighboringType) - if (anchor /= 0) then ! valid anchor node - if (any(FE_face(:,matchingFace,neighboringType) == anchor)) then ! sits on matching face? - NmatchingNodes = NmatchingNodes + 1 - matchingNodes(NmatchingNodes) = mesh_element(4+anchor,matchingElem) ! CP id of neighbor's anchor node - else ! no matching, because not all nodes sit on the matching face - NmatchingNodes = 0 - matchingNodes = 0 - exit - endif - endif - enddo - - if (NmatchingNodes /= NlinkedNodes) & ! this ip has wrong count of anchors on face - cycle checkCandidateIP - - !*** check "normal" nodes whether they match or not - - checkTwins = .false. - do a = 1,NlinkedNodes - if (all(matchingNodes /= linkedNodes(a))) then ! this linkedNode does not match any matchingNode - checkTwins = .true. - exit ! no need to search further - endif - enddo - - !*** if no match found, then also check node twins - - if(checkTwins) then - dir = int(maxloc(abs(mesh_ipAreaNormal(1:3,neighbor,myIP,myElem)),1),pInt) ! check for twins only in direction of the surface normal - do a = 1,NlinkedNodes - twin_of_linkedNode = mesh_nodeTwins(dir,linkedNodes(a)) - if (twin_of_linkedNode == 0 .or. & ! twin of linkedNode does not exist... - all(matchingNodes /= twin_of_linkedNode)) then ! ... or it does not match any matchingNode - cycle checkCandidateIP ! ... then check next candidateIP - endif - enddo - endif - - !*** we found a match !!! - - mesh_ipNeighborhood(1,neighbor,myIP,myElem) = matchingElem - mesh_ipNeighborhood(2,neighbor,myIP,myElem) = candidateIP - exit checkCandidateIP - enddo checkCandidateIP - endif ! end of valid external matching - endif ! end of internal/external matching - enddo - enddo - enddo - do myElem = 1,mesh_NcpElems ! loop over cpElems - myType = FE_geomtype(mesh_element(2,myElem)) ! get elemGeomType - do myIP = 1,FE_Nips(myType) ! loop over IPs of elem - do neighbor = 1,FE_NipNeighbors(FE_celltype(myType)) ! loop over neighbors of IP - neighboringElem = mesh_ipNeighborhood(1,neighbor,myIP,myElem) - neighboringIP = mesh_ipNeighborhood(2,neighbor,myIP,myElem) - if (neighboringElem > 0 .and. neighboringIP > 0) then ! if neighbor exists ... - neighboringType = FE_geomtype(mesh_element(2,neighboringElem)) - do pointingToMe = 1,FE_NipNeighbors(FE_celltype(neighboringType)) ! find neighboring index that points from my neighbor to myself - if ( myElem == mesh_ipNeighborhood(1,pointingToMe,neighboringIP,neighboringElem) & - .and. myIP == mesh_ipNeighborhood(2,pointingToMe,neighboringIP,neighboringElem)) then ! possible candidate - if (math_inner(mesh_ipAreaNormal(1:3,neighbor,myIP,myElem),& - mesh_ipAreaNormal(1:3,pointingToMe,neighboringIP,neighboringElem)) < 0.0_pReal) then ! area normals have opposite orientation (we have to check that because of special case for single element with two ips and periodicity. In this case the neighbor is identical in two different directions.) - mesh_ipNeighborhood(3,neighbor,myIP,myElem) = pointingToMe ! found match - exit ! so no need to search further - endif - endif - enddo - endif - enddo - enddo - enddo - - call geometry_plastic_nonlocal_set_IPneighborhood(mesh_ipNeighborhood) - - contains - !-------------------------------------------------------------------------------------------------- -!> @brief find face-matching element of same type -!-------------------------------------------------------------------------------------------------- -subroutine mesh_faceMatch(elem, face ,matchingElem, matchingFace) - - -integer, intent(out) :: matchingElem, & ! matching CP element ID - matchingFace ! matching face ID -integer, intent(in) :: face, & ! face ID - elem ! CP elem ID -integer, dimension(FE_NmatchingNodesPerFace(face,FE_geomtype(mesh_element(2,elem)))) :: & - myFaceNodes ! global node ids on my face -integer :: myType, & - candidateType, & - candidateElem, & - candidateFace, & - candidateFaceNode, & - minNsharedElems, & - NsharedElems, & - lonelyNode = 0, & - i, & - n, & - dir ! periodicity direction -integer, dimension(:), allocatable :: element_seen -logical checkTwins - -matchingElem = 0 -matchingFace = 0 -minNsharedElems = mesh_maxNsharedElems + 1 ! init to worst case -myType = FE_geomtype(mesh_element(2,elem)) ! figure elemGeomType - -do n = 1,FE_NmatchingNodesPerFace(face,myType) ! loop over nodes on face - myFaceNodes(n) = mesh_element(4+FE_face(n,face,myType),elem) ! CP id of face node - NsharedElems = mesh_sharedElem(1,myFaceNodes(n)) ! figure # shared elements for this node - if (NsharedElems < minNsharedElems) then - minNsharedElems = NsharedElems ! remember min # shared elems - lonelyNode = n ! remember most lonely node - endif -enddo - -allocate(element_seen(minNsharedElems)) -element_seen = 0 - -checkCandidate: do i = 1,minNsharedElems ! iterate over lonelyNode's shared elements - candidateElem = mesh_sharedElem(1+i,myFaceNodes(lonelyNode)) ! present candidate elem - if (all(element_seen /= candidateElem)) then ! element seen for the first time? - element_seen(i) = candidateElem - candidateType = FE_geomtype(mesh_element(2,candidateElem)) ! figure elemGeomType of candidate -checkCandidateFace: do candidateFace = 1,FE_maxNipNeighbors ! check each face of candidate - if (FE_NmatchingNodesPerFace(candidateFace,candidateType) & - /= FE_NmatchingNodesPerFace(face,myType) & ! incompatible face - .or. (candidateElem == elem .and. candidateFace == face)) then ! this is my face - cycle checkCandidateFace - endif - checkTwins = .false. - do n = 1,FE_NmatchingNodesPerFace(candidateFace,candidateType) ! loop through nodes on face - candidateFaceNode = mesh_element(4+FE_face(n,candidateFace,candidateType),candidateElem) - if (all(myFaceNodes /= candidateFaceNode)) then ! candidate node does not match any of my face nodes - checkTwins = .true. ! perhaps the twin nodes do match - exit - endif - enddo - if(checkTwins) then -checkCandidateFaceTwins: do dir = 1,3 - do n = 1,FE_NmatchingNodesPerFace(candidateFace,candidateType) ! loop through nodes on face - candidateFaceNode = mesh_element(4+FE_face(n,candidateFace,candidateType),candidateElem) - if (all(myFaceNodes /= mesh_nodeTwins(dir,candidateFaceNode))) then ! node twin does not match either - if (dir == 3) then - cycle checkCandidateFace - else - cycle checkCandidateFaceTwins ! try twins in next dimension - endif - endif - enddo - exit checkCandidateFaceTwins - enddo checkCandidateFaceTwins - endif - matchingFace = candidateFace - matchingElem = candidateElem - exit checkCandidate ! found my matching candidate - enddo checkCandidateFace - endif -enddo checkCandidate - -end subroutine mesh_faceMatch - -end subroutine mesh_build_ipNeighborhood - - -!-------------------------------------------------------------------------------------------------- -!> @brief mapping of FE element types to internal representation -!-------------------------------------------------------------------------------------------------- -integer function FE_mapElemtype(what) - use IO, only: IO_lc, IO_error - - - character(len=*), intent(in) :: what - - select case (IO_lc(what)) - case ( 'cpe4', & - 'cpe4t') - FE_mapElemtype = 3 ! Arbitrary Quadrilateral Plane-strain - case ( 'cpe8', & - 'cpe8t') - FE_mapElemtype = 4 ! Plane Strain, Eight-node Distorted Quadrilateral - case ( 'c3d4', & - 'c3d4t') - FE_mapElemtype = 6 ! Three-dimensional Four-node Tetrahedron - case ( 'c3d6', & - 'c3d6t') - FE_mapElemtype = 9 ! Three-dimensional Arbitrarily Distorted Pentahedral - case ( 'c3d8r', & - 'c3d8rt') - FE_mapElemtype = 10 ! Three-dimensional Arbitrarily Distorted linear hexahedral with reduced integration - case ( 'c3d8', & - 'c3d8t') - FE_mapElemtype = 11 ! Three-dimensional Arbitrarily Distorted Brick - case ( 'c3d20r', & - 'c3d20rt') - FE_mapElemtype = 12 ! Three-dimensional Arbitrarily Distorted quad hexahedral with reduced integration - case ( 'c3d20', & - 'c3d20t') - FE_mapElemtype = 13 ! Three-dimensional Arbitrarily Distorted quadratic hexahedral - case default - call IO_error(error_ID=190,ext_msg=IO_lc(what)) - end select - -end function FE_mapElemtype - - - - - -!-------------------------------------------------------------------------------------------------- -!> @brief get properties of different types of finite elements -!> @details assign globals: FE_nodesAtIP, FE_ipNeighbor, FE_cellnodeParentnodeWeights, FE_subNodeOnIPFace -!-------------------------------------------------------------------------------------------------- -subroutine mesh_build_FEdata - - - integer :: me - allocate(FE_nodesAtIP(FE_maxmaxNnodesAtIP,FE_maxNips,FE_Ngeomtypes), source=0) - allocate(FE_ipNeighbor(FE_maxNipNeighbors,FE_maxNips,FE_Ngeomtypes), source=0) - allocate(FE_cell(FE_maxNcellnodesPerCell,FE_maxNips,FE_Ngeomtypes), source=0) - allocate(FE_cellnodeParentnodeWeights(FE_maxNnodes,FE_maxNcellnodes,FE_Nelemtypes), source=0.0_pReal) - allocate(FE_cellface(FE_maxNcellnodesPerCellface,FE_maxNcellfaces,FE_Ncelltypes), source=0) - - - !*** fill FE_nodesAtIP with data *** - - me = 0 - - me = me + 1 - FE_nodesAtIP(1:FE_maxNnodesAtIP(me),1:FE_Nips(me),me) = & ! element 6 (2D 3node 1ip) - reshape(int([& - 1,2,3 & - ],pInt),[FE_maxNnodesAtIP(me),FE_Nips(me)]) - - me = me + 1 - FE_nodesAtIP(1:FE_maxNnodesAtIP(me),1:FE_Nips(me),me) = & ! element 125 (2D 6node 3ip) - reshape(int([& - 1, & - 2, & - 3 & - ],pInt),[FE_maxNnodesAtIP(me),FE_Nips(me)]) - - me = me + 1 - FE_nodesAtIP(1:FE_maxNnodesAtIP(me),1:FE_Nips(me),me) = & ! element 11 (2D 4node 4ip) - reshape(int([& - 1, & - 2, & - 4, & - 3 & - ],pInt),[FE_maxNnodesAtIP(me),FE_Nips(me)]) - - me = me + 1 - FE_nodesAtIP(1:FE_maxNnodesAtIP(me),1:FE_Nips(me),me) = & ! element 27 (2D 8node 9ip) - reshape(int([& - 1,0, & - 1,2, & - 2,0, & - 1,4, & - 0,0, & - 2,3, & - 4,0, & - 3,4, & - 3,0 & - ],pInt),[FE_maxNnodesAtIP(me),FE_Nips(me)]) - - me = me + 1 - FE_nodesAtIP(1:FE_maxNnodesAtIP(me),1:FE_Nips(me),me) = & ! element 134 (3D 4node 1ip) - reshape(int([& - 1,2,3,4 & - ],pInt),[FE_maxNnodesAtIP(me),FE_Nips(me)]) - - me = me + 1 - FE_nodesAtIP(1:FE_maxNnodesAtIP(me),1:FE_Nips(me),me) = & ! element 127 (3D 10node 4ip) - reshape(int([& - 1, & - 2, & - 3, & - 4 & - ],pInt),[FE_maxNnodesAtIP(me),FE_Nips(me)]) - - me = me + 1 - FE_nodesAtIP(1:FE_maxNnodesAtIP(me),1:FE_Nips(me),me) = & ! element 136 (3D 6node 6ip) - reshape(int([& - 1, & - 2, & - 3, & - 4, & - 5, & - 6 & - ],pInt),[FE_maxNnodesAtIP(me),FE_Nips(me)]) - - me = me + 1 - FE_nodesAtIP(1:FE_maxNnodesAtIP(me),1:FE_Nips(me),me) = & ! element 117 (3D 8node 1ip) - reshape(int([& - 1,2,3,4,5,6,7,8 & - ],pInt),[FE_maxNnodesAtIP(me),FE_Nips(me)]) - - me = me + 1 - FE_nodesAtIP(1:FE_maxNnodesAtIP(me),1:FE_Nips(me),me) = & ! element 7 (3D 8node 8ip) - reshape(int([& - 1, & - 2, & - 4, & - 3, & - 5, & - 6, & - 8, & - 7 & - ],pInt),[FE_maxNnodesAtIP(me),FE_Nips(me)]) - - me = me + 1 - FE_nodesAtIP(1:FE_maxNnodesAtIP(me),1:FE_Nips(me),me) = & ! element 21 (3D 20node 27ip) - reshape(int([& - 1,0, 0,0, & - 1,2, 0,0, & - 2,0, 0,0, & - 1,4, 0,0, & - 1,3, 2,4, & - 2,3, 0,0, & - 4,0, 0,0, & - 3,4, 0,0, & - 3,0, 0,0, & - 1,5, 0,0, & - 1,6, 2,5, & - 2,6, 0,0, & - 1,8, 4,5, & - 0,0, 0,0, & - 2,7, 3,6, & - 4,8, 0,0, & - 3,8, 4,7, & - 3,7, 0,0, & - 5,0, 0,0, & - 5,6, 0,0, & - 6,0, 0,0, & - 5,8, 0,0, & - 5,7, 6,8, & - 6,7, 0,0, & - 8,0, 0,0, & - 7,8, 0,0, & - 7,0, 0,0 & - ],pInt),[FE_maxNnodesAtIP(me),FE_Nips(me)]) - - - ! *** FE_ipNeighbor *** - ! is a list of the neighborhood of each IP. - ! It is sorted in (local) +x,-x, +y,-y, +z,-z direction. - ! Positive integers denote an intra-FE IP identifier. - ! Negative integers denote the interface behind which the neighboring (extra-FE) IP will be located. - me = 0 - - me = me + 1 - FE_ipNeighbor(1:FE_NipNeighbors(FE_celltype(me)),1:FE_Nips(me),me) = & ! element 6 (2D 3node 1ip) - reshape(int([& - -2,-3,-1 & - ],pInt),[FE_NipNeighbors(FE_celltype(me)),FE_Nips(me)]) - - me = me + 1 - FE_ipNeighbor(1:FE_NipNeighbors(FE_celltype(me)),1:FE_Nips(me),me) = & ! element 125 (2D 6node 3ip) - reshape(int([& - 2,-3, 3,-1, & - -2, 1, 3,-1, & - 2,-3,-2, 1 & - ],pInt),[FE_NipNeighbors(FE_celltype(me)),FE_Nips(me)]) - - me = me + 1 - FE_ipNeighbor(1:FE_NipNeighbors(FE_celltype(me)),1:FE_Nips(me),me) = & ! element 11 (2D 4node 4ip) - reshape(int([& - 2,-4, 3,-1, & - -2, 1, 4,-1, & - 4,-4,-3, 1, & - -2, 3,-3, 2 & - ],pInt),[FE_NipNeighbors(FE_celltype(me)),FE_Nips(me)]) - - me = me + 1 - FE_ipNeighbor(1:FE_NipNeighbors(FE_celltype(me)),1:FE_Nips(me),me) = & ! element 27 (2D 8node 9ip) - reshape(int([& - 2,-4, 4,-1, & - 3, 1, 5,-1, & - -2, 2, 6,-1, & - 5,-4, 7, 1, & - 6, 4, 8, 2, & - -2, 5, 9, 3, & - 8,-4,-3, 4, & - 9, 7,-3, 5, & - -2, 8,-3, 6 & - ],pInt),[FE_NipNeighbors(FE_celltype(me)),FE_Nips(me)]) - - me = me + 1 - FE_ipNeighbor(1:FE_NipNeighbors(FE_celltype(me)),1:FE_Nips(me),me) = & ! element 134 (3D 4node 1ip) - reshape(int([& - -1,-2,-3,-4 & - ],pInt),[FE_NipNeighbors(FE_celltype(me)),FE_Nips(me)]) - - me = me + 1 - FE_ipNeighbor(1:FE_NipNeighbors(FE_celltype(me)),1:FE_Nips(me),me) = & ! element 127 (3D 10node 4ip) - reshape(int([& - 2,-4, 3,-2, 4,-1, & - -2, 1, 3,-2, 4,-1, & - 2,-4,-3, 1, 4,-1, & - 2,-4, 3,-2,-3, 1 & - ],pInt),[FE_NipNeighbors(FE_celltype(me)),FE_Nips(me)]) - - me = me + 1 - FE_ipNeighbor(1:FE_NipNeighbors(FE_celltype(me)),1:FE_Nips(me),me) = & ! element 136 (3D 6node 6ip) - reshape(int([& - 2,-4, 3,-2, 4,-1, & - -3, 1, 3,-2, 5,-1, & - 2,-4,-3, 1, 6,-1, & - 5,-4, 6,-2,-5, 1, & - -3, 4, 6,-2,-5, 2, & - 5,-4,-3, 4,-5, 3 & - ],pInt),[FE_NipNeighbors(FE_celltype(me)),FE_Nips(me)]) - - me = me + 1 - FE_ipNeighbor(1:FE_NipNeighbors(FE_celltype(me)),1:FE_Nips(me),me) = & ! element 117 (3D 8node 1ip) - reshape(int([& - -3,-5,-4,-2,-6,-1 & - ],pInt),[FE_NipNeighbors(FE_celltype(me)),FE_Nips(me)]) - - me = me + 1 - FE_ipNeighbor(1:FE_NipNeighbors(FE_celltype(me)),1:FE_Nips(me),me) = & ! element 7 (3D 8node 8ip) - reshape(int([& - 2,-5, 3,-2, 5,-1, & - -3, 1, 4,-2, 6,-1, & - 4,-5,-4, 1, 7,-1, & - -3, 3,-4, 2, 8,-1, & - 6,-5, 7,-2,-6, 1, & - -3, 5, 8,-2,-6, 2, & - 8,-5,-4, 5,-6, 3, & - -3, 7,-4, 6,-6, 4 & - ],pInt),[FE_NipNeighbors(FE_celltype(me)),FE_Nips(me)]) - - me = me + 1 - FE_ipNeighbor(1:FE_NipNeighbors(FE_celltype(me)),1:FE_Nips(me),me) = & ! element 21 (3D 20node 27ip) - reshape(int([& - 2,-5, 4,-2,10,-1, & - 3, 1, 5,-2,11,-1, & - -3, 2, 6,-2,12,-1, & - 5,-5, 7, 1,13,-1, & - 6, 4, 8, 2,14,-1, & - -3, 5, 9, 3,15,-1, & - 8,-5,-4, 4,16,-1, & - 9, 7,-4, 5,17,-1, & - -3, 8,-4, 6,18,-1, & - 11,-5,13,-2,19, 1, & - 12,10,14,-2,20, 2, & - -3,11,15,-2,21, 3, & - 14,-5,16,10,22, 4, & - 15,13,17,11,23, 5, & - -3,14,18,12,24, 6, & - 17,-5,-4,13,25, 7, & - 18,16,-4,14,26, 8, & - -3,17,-4,15,27, 9, & - 20,-5,22,-2,-6,10, & - 21,19,23,-2,-6,11, & - -3,20,24,-2,-6,12, & - 23,-5,25,19,-6,13, & - 24,22,26,20,-6,14, & - -3,23,27,21,-6,15, & - 26,-5,-4,22,-6,16, & - 27,25,-4,23,-6,17, & - -3,26,-4,24,-6,18 & - ],pInt),[FE_NipNeighbors(FE_celltype(me)),FE_Nips(me)]) - - - ! *** FE_cell *** - me = 0 - - me = me + 1 - FE_cell(1:FE_NcellnodesPerCell(FE_celltype(me)),1:FE_Nips(me),me) = & ! element 6 (2D 3node 1ip) - reshape(int([& - 1,2,3 & - ],pInt),[FE_NcellnodesPerCell(FE_celltype(me)),FE_Nips(me)]) - - me = me + 1 - FE_cell(1:FE_NcellnodesPerCell(FE_celltype(me)),1:FE_Nips(me),me) = & ! element 125 (2D 6node 3ip) - reshape(int([& - 1, 4, 7, 6, & - 2, 5, 7, 4, & - 3, 6, 7, 5 & - ],pInt),[FE_NcellnodesPerCell(FE_celltype(me)),FE_Nips(me)]) - - me = me + 1 - FE_cell(1:FE_NcellnodesPerCell(FE_celltype(me)),1:FE_Nips(me),me) = & ! element 11 (2D 4node 4ip) - reshape(int([& - 1, 5, 9, 8, & - 5, 2, 6, 9, & - 8, 9, 7, 4, & - 9, 6, 3, 7 & - ],pInt),[FE_NcellnodesPerCell(FE_celltype(me)),FE_Nips(me)]) - - me = me + 1 - FE_cell(1:FE_NcellnodesPerCell(FE_celltype(me)),1:FE_Nips(me),me) = & ! element 27 (2D 8node 9ip) - reshape(int([& - 1, 5,13,12, & - 5, 6,14,13, & - 6, 2, 7,14, & - 12,13,16,11, & - 13,14,15,16, & - 14, 7, 8,15, & - 11,16,10, 4, & - 16,15, 9,10, & - 15, 8, 3, 9 & - ],pInt),[FE_NcellnodesPerCell(FE_celltype(me)),FE_Nips(me)]) - - me = me + 1 - FE_cell(1:FE_NcellnodesPerCell(FE_celltype(me)),1:FE_Nips(me),me) = & ! element 134 (3D 4node 1ip) - reshape(int([& - 1, 2, 3, 4 & - ],pInt),[FE_NcellnodesPerCell(FE_celltype(me)),FE_Nips(me)]) - - me = me + 1 - FE_cell(1:FE_NcellnodesPerCell(FE_celltype(me)),1:FE_Nips(me),me) = & ! element 127 (3D 10node 4ip) - reshape(int([& - 1, 5,11, 7, 8,12,15,14, & - 5, 2, 6,11,12, 9,13,15, & - 7,11, 6, 3,14,15,13,10, & - 8,12,15, 4, 4, 9,13,10 & - ],pInt),[FE_NcellnodesPerCell(FE_celltype(me)),FE_Nips(me)]) - - me = me + 1 - FE_cell(1:FE_NcellnodesPerCell(FE_celltype(me)),1:FE_Nips(me),me) = & ! element 136 (3D 6node 6ip) - reshape(int([& - 1, 7,16, 9,10,17,21,19, & - 7, 2, 8,16,17,11,18,21, & - 9,16, 8, 3,19,21,18,12, & - 10,17,21,19, 4,13,20,15, & - 17,11,18,21,13, 5,14,20, & - 19,21,18,12,15,20,14, 6 & - ],pInt),[FE_NcellnodesPerCell(FE_celltype(me)),FE_Nips(me)]) - - me = me + 1 - FE_cell(1:FE_NcellnodesPerCell(FE_celltype(me)),1:FE_Nips(me),me) = & ! element 117 (3D 8node 1ip) - reshape(int([& - 1, 2, 3, 4, 5, 6, 7, 8 & - ],pInt),[FE_NcellnodesPerCell(FE_celltype(me)),FE_Nips(me)]) - - me = me + 1 - FE_cell(1:FE_NcellnodesPerCell(FE_celltype(me)),1:FE_Nips(me),me) = & ! element 7 (3D 8node 8ip) - reshape(int([& - 1, 9,21,12,13,22,27,25, & - 9, 2,10,21,22,14,23,27, & - 12,21,11, 4,25,27,24,16, & - 21,10, 3,11,27,23,15,24, & - 13,22,27,25, 5,17,26,20, & - 22,14,23,27,17, 6,18,26, & - 25,27,24,16,20,26,19, 8, & - 27,23,15,24,26,18, 7,19 & - ],pInt),[FE_NcellnodesPerCell(FE_celltype(me)),FE_Nips(me)]) - - me = me + 1 - FE_cell(1:FE_NcellnodesPerCell(FE_celltype(me)),1:FE_Nips(me),me) = & ! element 21 (3D 20node 27ip) - reshape(int([& - 1, 9,33,16,17,37,57,44, & - 9,10,34,33,37,38,58,57, & - 10, 2,11,34,38,18,39,58, & - 16,33,36,15,44,57,60,43, & - 33,34,35,36,57,58,59,60, & - 34,11,12,35,58,39,40,59, & - 15,36,14, 4,43,60,42,20, & - 36,35,13,14,60,59,41,42, & - 35,12, 3,13,59,40,19,41, & - 17,37,57,44,21,45,61,52, & - 37,38,58,57,45,46,62,61, & - 38,18,39,58,46,22,47,62, & - 44,57,60,43,52,61,64,51, & - 57,58,59,60,61,62,63,64, & - 58,39,40,59,62,47,48,63, & - 43,60,42,20,51,64,50,24, & - 60,59,41,42,64,63,49,50, & - 59,40,19,41,63,48,23,49, & - 21,45,61,52, 5,25,53,32, & - 45,46,62,61,25,26,54,53, & - 46,22,47,62,26, 6,27,54, & - 52,61,64,51,32,53,56,31, & - 61,62,63,64,53,54,55,56, & - 62,47,48,63,54,27,28,55, & - 51,64,50,24,31,56,30, 8, & - 64,63,49,50,56,55,29,30, & - 63,48,23,49,55,28, 7,29 & - ],pInt),[FE_NcellnodesPerCell(FE_celltype(me)),FE_Nips(me)]) - - - ! *** FE_cellnodeParentnodeWeights *** - ! center of gravity of the weighted nodes gives the position of the cell node. - ! fill with 0. - ! example: face-centered cell node with face nodes 1,2,5,6 to be used in, - ! e.g., an 8 node element, would be encoded: - ! 1, 1, 0, 0, 1, 1, 0, 0 - me = 0 - - me = me + 1 - FE_cellnodeParentnodeWeights(1:FE_Nnodes(me),1:FE_Ncellnodes(FE_geomtype(me)),me) = & ! element 6 (2D 3node 1ip) - reshape(real([& - 1, 0, 0, & - 0, 1, 0, & - 0, 0, 1 & - ],pReal),[FE_Nnodes(me),FE_Ncellnodes(FE_geomtype(me))]) - - me = me + 1 - FE_cellnodeParentnodeWeights(1:FE_Nnodes(me),1:FE_Ncellnodes(FE_geomtype(me)),me) = & ! element 125 (2D 6node 3ip) - reshape(real([& - 1, 0, 0, 0, 0, 0, & - 0, 1, 0, 0, 0, 0, & - 0, 0, 1, 0, 0, 0, & - 0, 0, 0, 1, 0, 0, & - 0, 0, 0, 0, 1, 0, & - 0, 0, 0, 0, 0, 1, & - 1, 1, 1, 2, 2, 2 & - ],pReal),[FE_Nnodes(me),FE_Ncellnodes(FE_geomtype(me))]) - - me = me + 1 - FE_cellnodeParentnodeWeights(1:FE_Nnodes(me),1:FE_Ncellnodes(FE_geomtype(me)),me) = & ! element 11 (2D 4node 4ip) - reshape(real([& - 1, 0, 0, 0, & - 0, 1, 0, 0, & - 0, 0, 1, 0, & - 0, 0, 0, 1, & - 1, 1, 0, 0, & - 0, 1, 1, 0, & - 0, 0, 1, 1, & - 1, 0, 0, 1, & - 1, 1, 1, 1 & - ],pReal),[FE_Nnodes(me),FE_Ncellnodes(FE_geomtype(me))]) - - me = me + 1 - FE_cellnodeParentnodeWeights(1:FE_Nnodes(me),1:FE_Ncellnodes(FE_geomtype(me)),me) = & ! element 27 (2D 8node 9ip) - reshape(real([& - 1, 0, 0, 0, 0, 0, 0, 0, & - 0, 1, 0, 0, 0, 0, 0, 0, & - 0, 0, 1, 0, 0, 0, 0, 0, & - 0, 0, 0, 1, 0, 0, 0, 0, & - 1, 0, 0, 0, 2, 0, 0, 0, & - 0, 1, 0, 0, 2, 0, 0, 0, & - 0, 1, 0, 0, 0, 2, 0, 0, & - 0, 0, 1, 0, 0, 2, 0, 0, & - 0, 0, 1, 0, 0, 0, 2, 0, & - 0, 0, 0, 1, 0, 0, 2, 0, & - 0, 0, 0, 1, 0, 0, 0, 2, & - 1, 0, 0, 0, 0, 0, 0, 2, & - 4, 1, 1, 1, 8, 2, 2, 8, & - 1, 4, 1, 1, 8, 8, 2, 2, & - 1, 1, 4, 1, 2, 8, 8, 2, & - 1, 1, 1, 4, 2, 2, 8, 8 & - ],pReal),[FE_Nnodes(me),FE_Ncellnodes(FE_geomtype(me))]) - - me = me + 1 - FE_cellnodeParentnodeWeights(1:FE_Nnodes(me),1:FE_Ncellnodes(FE_geomtype(me)),me) = & ! element 54 (2D 8node 4ip) - reshape(real([& - 1, 0, 0, 0, 0, 0, 0, 0, & - 0, 1, 0, 0, 0, 0, 0, 0, & - 0, 0, 1, 0, 0, 0, 0, 0, & - 0, 0, 0, 1, 0, 0, 0, 0, & - 0, 0, 0, 0, 1, 0, 0, 0, & - 0, 0, 0, 0, 0, 1, 0, 0, & - 0, 0, 0, 0, 0, 0, 1, 0, & - 0, 0, 0, 0, 0, 0, 0, 1, & - 1, 1, 1, 1, 2, 2, 2, 2 & - ],pReal),[FE_Nnodes(me),FE_Ncellnodes(FE_geomtype(me))]) - - me = me + 1 - FE_cellnodeParentnodeWeights(1:FE_Nnodes(me),1:FE_Ncellnodes(FE_geomtype(me)),me) = & ! element 134 (3D 4node 1ip) - reshape(real([& - 1, 0, 0, 0, & - 0, 1, 0, 0, & - 0, 0, 1, 0, & - 0, 0, 0, 1 & - ],pReal),[FE_Nnodes(me),FE_Ncellnodes(FE_geomtype(me))]) - - me = me + 1 - FE_cellnodeParentnodeWeights(1:FE_Nnodes(me),1:FE_Ncellnodes(FE_geomtype(me)),me) = & ! element 157 (3D 5node 4ip) - reshape(real([& - 1, 0, 0, 0, 0, & - 0, 1, 0, 0, 0, & - 0, 0, 1, 0, 0, & - 0, 0, 0, 1, 0, & - 1, 1, 0, 0, 0, & - 0, 1, 1, 0, 0, & - 1, 0, 1, 0, 0, & - 1, 0, 0, 1, 0, & - 0, 1, 0, 1, 0, & - 0, 0, 1, 1, 0, & - 1, 1, 1, 0, 0, & - 1, 1, 0, 1, 0, & - 0, 1, 1, 1, 0, & - 1, 0, 1, 1, 0, & - 0, 0, 0, 0, 1 & - ],pReal),[FE_Nnodes(me),FE_Ncellnodes(FE_geomtype(me))]) - - me = me + 1 - FE_cellnodeParentnodeWeights(1:FE_Nnodes(me),1:FE_Ncellnodes(FE_geomtype(me)),me) = & ! element 127 (3D 10node 4ip) - reshape(real([& - 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, & - 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, & - 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, & - 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, & - 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, & - 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, & - 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, & - 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, & - 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, & - 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, & - 1, 1, 1, 0, 2, 2, 2, 0, 0, 0, & - 1, 1, 0, 1, 2, 0, 0, 2, 2, 0, & - 0, 1, 1, 1, 0, 2, 0, 0, 2, 2, & - 1, 0, 1, 1, 0, 0, 2, 2, 0, 2, & - 3, 3, 3, 3, 4, 4, 4, 4, 4, 4 & - ],pReal),[FE_Nnodes(me),FE_Ncellnodes(FE_geomtype(me))]) - - me = me + 1 - FE_cellnodeParentnodeWeights(1:FE_Nnodes(me),1:FE_Ncellnodes(FE_geomtype(me)),me) = & ! element 136 (3D 6node 6ip) - reshape(real([& - 1, 0, 0, 0, 0, 0, & - 0, 1, 0, 0, 0, 0, & - 0, 0, 1, 0, 0, 0, & - 0, 0, 0, 1, 0, 0, & - 0, 0, 0, 0, 1, 0, & - 0, 0, 0, 0, 0, 1, & - 1, 1, 0, 0, 0, 0, & - 0, 1, 1, 0, 0, 0, & - 1, 0, 1, 0, 0, 0, & - 1, 0, 0, 1, 0, 0, & - 0, 1, 0, 0, 1, 0, & - 0, 0, 1, 0, 0, 1, & - 0, 0, 0, 1, 1, 0, & - 0, 0, 0, 0, 1, 1, & - 0, 0, 0, 1, 0, 1, & - 1, 1, 1, 0, 0, 0, & - 1, 1, 0, 1, 1, 0, & - 0, 1, 1, 0, 1, 1, & - 1, 0, 1, 1, 0, 1, & - 0, 0, 0, 1, 1, 1, & - 1, 1, 1, 1, 1, 1 & - ],pReal),[FE_Nnodes(me),FE_Ncellnodes(FE_geomtype(me))]) - - me = me + 1 - FE_cellnodeParentnodeWeights(1:FE_Nnodes(me),1:FE_Ncellnodes(FE_geomtype(me)),me) = & ! element 117 (3D 8node 1ip) - reshape(real([& - 1, 0, 0, 0, 0, 0, 0, 0, & - 0, 1, 0, 0, 0, 0, 0, 0, & - 0, 0, 1, 0, 0, 0, 0, 0, & - 0, 0, 0, 1, 0, 0, 0, 0, & - 0, 0, 0, 0, 1, 0, 0, 0, & - 0, 0, 0, 0, 0, 1, 0, 0, & - 0, 0, 0, 0, 0, 0, 1, 0, & - 0, 0, 0, 0, 0, 0, 0, 1 & - ],pReal),[FE_Nnodes(me),FE_Ncellnodes(FE_geomtype(me))]) - - me = me + 1 - FE_cellnodeParentnodeWeights(1:FE_Nnodes(me),1:FE_Ncellnodes(FE_geomtype(me)),me) = & ! element 7 (3D 8node 8ip) - reshape(real([& - 1, 0, 0, 0, 0, 0, 0, 0, & ! - 0, 1, 0, 0, 0, 0, 0, 0, & ! - 0, 0, 1, 0, 0, 0, 0, 0, & ! - 0, 0, 0, 1, 0, 0, 0, 0, & ! - 0, 0, 0, 0, 1, 0, 0, 0, & ! 5 - 0, 0, 0, 0, 0, 1, 0, 0, & ! - 0, 0, 0, 0, 0, 0, 1, 0, & ! - 0, 0, 0, 0, 0, 0, 0, 1, & ! - 1, 1, 0, 0, 0, 0, 0, 0, & ! - 0, 1, 1, 0, 0, 0, 0, 0, & ! 10 - 0, 0, 1, 1, 0, 0, 0, 0, & ! - 1, 0, 0, 1, 0, 0, 0, 0, & ! - 1, 0, 0, 0, 1, 0, 0, 0, & ! - 0, 1, 0, 0, 0, 1, 0, 0, & ! - 0, 0, 1, 0, 0, 0, 1, 0, & ! 15 - 0, 0, 0, 1, 0, 0, 0, 1, & ! - 0, 0, 0, 0, 1, 1, 0, 0, & ! - 0, 0, 0, 0, 0, 1, 1, 0, & ! - 0, 0, 0, 0, 0, 0, 1, 1, & ! - 0, 0, 0, 0, 1, 0, 0, 1, & ! 20 - 1, 1, 1, 1, 0, 0, 0, 0, & ! - 1, 1, 0, 0, 1, 1, 0, 0, & ! - 0, 1, 1, 0, 0, 1, 1, 0, & ! - 0, 0, 1, 1, 0, 0, 1, 1, & ! - 1, 0, 0, 1, 1, 0, 0, 1, & ! 25 - 0, 0, 0, 0, 1, 1, 1, 1, & ! - 1, 1, 1, 1, 1, 1, 1, 1 & ! - ],pReal),[FE_Nnodes(me),FE_Ncellnodes(FE_geomtype(me))]) - - me = me + 1 - FE_cellnodeParentnodeWeights(1:FE_Nnodes(me),1:FE_Ncellnodes(FE_geomtype(me)),me) = & ! element 57 (3D 20node 8ip) - reshape(real([& - 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! - 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! - 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! - 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! - 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! 5 - 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! - 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! - 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! - 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! - 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! 10 - 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! - 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, & ! - 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, & ! - 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, & ! - 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, & ! 15 - 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, & ! - 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, & ! - 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, & ! - 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, & ! - 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, & ! 20 - 1, 1, 1, 1, 0, 0, 0, 0, 2, 2, 2, 2, 0, 0, 0, 0, 0, 0, 0, 0, & ! - 1, 1, 0, 0, 1, 1, 0, 0, 2, 0, 0, 0, 2, 0, 0, 0, 2, 2, 0, 0, & ! - 0, 1, 1, 0, 0, 1, 1, 0, 0, 2, 0, 0, 0, 2, 0, 0, 0, 2, 2, 0, & ! - 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 2, 0, 0, 0, 2, 0, 0, 0, 2, 2, & ! - 1, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 2, 0, 0, 0, 2, 2, 0, 0, 2, & ! 25 - 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 2, 2, 2, 2, 0, 0, 0, 0, & ! - 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4 & ! - ],pReal),[FE_Nnodes(me),FE_Ncellnodes(FE_geomtype(me))]) - - me = me + 1 - FE_cellnodeParentnodeWeights(1:FE_Nnodes(me),1:FE_Ncellnodes(FE_geomtype(me)),me) = & ! element 21 (3D 20node 27ip) - reshape(real([& - 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! - 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! - 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! - 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! - 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! 5 - 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! - 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! - 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! - 1, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! - 0, 1, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! 10 - 0, 1, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! - 0, 0, 1, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! - 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! - 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, & ! - 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, & ! 15 - 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, & ! - 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, & ! - 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, & ! - 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, & ! - 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, & ! 20 - 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, & ! - 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, & ! - 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, & ! - 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, & ! - 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, & ! 25 - 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, & ! - 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, & ! - 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, & ! - 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, & ! - 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, & ! 30 - 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, & ! - 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, & ! - 4, 1, 1, 1, 0, 0, 0, 0, 8, 2, 2, 8, 0, 0, 0, 0, 0, 0, 0, 0, & ! - 1, 4, 1, 1, 0, 0, 0, 0, 8, 8, 2, 2, 0, 0, 0, 0, 0, 0, 0, 0, & ! - 1, 1, 4, 1, 0, 0, 0, 0, 2, 8, 8, 2, 0, 0, 0, 0, 0, 0, 0, 0, & ! 35 - 1, 1, 1, 4, 0, 0, 0, 0, 2, 2, 8, 8, 0, 0, 0, 0, 0, 0, 0, 0, & ! - 4, 1, 0, 0, 1, 1, 0, 0, 8, 0, 0, 0, 2, 0, 0, 0, 8, 2, 0, 0, & ! - 1, 4, 0, 0, 1, 1, 0, 0, 8, 0, 0, 0, 2, 0, 0, 0, 2, 8, 0, 0, & ! - 0, 4, 1, 0, 0, 1, 1, 0, 0, 8, 0, 0, 0, 2, 0, 0, 0, 8, 2, 0, & ! - 0, 1, 4, 0, 0, 1, 1, 0, 0, 8, 0, 0, 0, 2, 0, 0, 0, 2, 8, 0, & ! 40 - 0, 0, 4, 1, 0, 0, 1, 1, 0, 0, 8, 0, 0, 0, 2, 0, 0, 0, 8, 2, & ! - 0, 0, 1, 4, 0, 0, 1, 1, 0, 0, 8, 0, 0, 0, 2, 0, 0, 0, 2, 8, & ! - 1, 0, 0, 4, 1, 0, 0, 1, 0, 0, 0, 8, 0, 0, 0, 2, 2, 0, 0, 8, & ! - 4, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 8, 0, 0, 0, 2, 8, 0, 0, 2, & ! - 1, 1, 0, 0, 4, 1, 0, 0, 2, 0, 0, 0, 8, 0, 0, 0, 8, 2, 0, 0, & ! 45 - 1, 1, 0, 0, 1, 4, 0, 0, 2, 0, 0, 0, 8, 0, 0, 0, 2, 8, 0, 0, & ! - 0, 1, 1, 0, 0, 4, 1, 0, 0, 2, 0, 0, 0, 8, 0, 0, 0, 8, 2, 0, & ! - 0, 1, 1, 0, 0, 1, 4, 0, 0, 2, 0, 0, 0, 8, 0, 0, 0, 2, 8, 0, & ! - 0, 0, 1, 1, 0, 0, 4, 1, 0, 0, 2, 0, 0, 0, 8, 0, 0, 0, 8, 2, & ! - 0, 0, 1, 1, 0, 0, 1, 4, 0, 0, 2, 0, 0, 0, 8, 0, 0, 0, 2, 8, & ! 50 - 1, 0, 0, 1, 1, 0, 0, 4, 0, 0, 0, 2, 0, 0, 0, 8, 2, 0, 0, 8, & ! - 1, 0, 0, 1, 4, 0, 0, 1, 0, 0, 0, 2, 0, 0, 0, 8, 8, 0, 0, 2, & ! - 0, 0, 0, 0, 4, 1, 1, 1, 0, 0, 0, 0, 8, 2, 2, 8, 0, 0, 0, 0, & ! - 0, 0, 0, 0, 1, 4, 1, 1, 0, 0, 0, 0, 8, 8, 2, 2, 0, 0, 0, 0, & ! - 0, 0, 0, 0, 1, 1, 4, 1, 0, 0, 0, 0, 2, 8, 8, 2, 0, 0, 0, 0, & ! 55 - 0, 0, 0, 0, 1, 1, 1, 4, 0, 0, 0, 0, 2, 2, 8, 8, 0, 0, 0, 0, & ! - 24, 8, 4, 8, 8, 4, 3, 4, 32,12,12,32, 12, 4, 4,12, 32,12, 4,12, & ! - 8,24, 8, 4, 4, 8, 4, 3, 32,32,12,12, 12,12, 4, 4, 12,32,12, 4, & ! - 4, 8,24, 8, 3, 4, 8, 4, 12,32,32,12, 4,12,12, 4, 4,12,32,12, & ! - 8, 4, 8,24, 4, 3, 4, 8, 12,12,32,32, 4, 4,12,12, 12, 4,12,32, & ! 60 - 8, 4, 3, 4, 24, 8, 4, 8, 12, 4, 4,12, 32,12,12,32, 32,12, 4,12, & ! - 4, 8, 4, 3, 8,24, 8, 4, 12,12, 4, 4, 32,32,12,12, 12,32,12, 4, & ! - 3, 4, 8, 4, 4, 8,24, 8, 4,12,12, 4, 12,32,32,12, 4,12,32,12, & ! - 4, 3, 4, 8, 8, 4, 8,24, 4, 4,12,12, 12,12,32,32, 12, 4,12,32 & ! - ],pReal),[FE_Nnodes(me),FE_Ncellnodes(FE_geomtype(me))]) - - - - ! *** FE_cellface *** - me = 0 - - me = me + 1 - FE_cellface(1:FE_NcellnodesPerCellface(me),1:FE_NipNeighbors(me),me) = & ! 2D 3node, VTK_TRIANGLE (5) - reshape(int([& - 2,3, & - 3,1, & - 1,2 & - ],pInt),[FE_NcellnodesPerCellface(me),FE_NipNeighbors(me)]) - - me = me + 1 - FE_cellface(1:FE_NcellnodesPerCellface(me),1:FE_NipNeighbors(me),me) = & ! 2D 4node, VTK_QUAD (9) - reshape(int([& - 2,3, & - 4,1, & - 3,4, & - 1,2 & - ],pInt),[FE_NcellnodesPerCellface(me),FE_NipNeighbors(me)]) - - me = me + 1 - FE_cellface(1:FE_NcellnodesPerCellface(me),1:FE_NipNeighbors(me),me) = & ! 3D 4node, VTK_TETRA (10) - reshape(int([& - 1,3,2, & - 1,2,4, & - 2,3,4, & - 1,4,3 & - ],pInt),[FE_NcellnodesPerCellface(me),FE_NipNeighbors(me)]) - - me = me + 1 - FE_cellface(1:FE_NcellnodesPerCellface(me),1:FE_NipNeighbors(me),me) = & ! 3D 8node, VTK_HEXAHEDRON (12) - reshape(int([& - 2,3,7,6, & - 4,1,5,8, & - 3,4,8,7, & - 1,2,6,5, & - 5,6,7,8, & - 1,4,3,2 & - ],pInt),[FE_NcellnodesPerCellface(me),FE_NipNeighbors(me)]) - - -end subroutine mesh_build_FEdata - - -!-------------------------------------------------------------------------------------------------- -!> @brief Gives the FE to CP ID mapping by binary search through lookup array -!! valid questions (what) are 'elem', 'node' -!-------------------------------------------------------------------------------------------------- -integer function mesh_FEasCP(what,myID) - use IO, only: & - IO_lc - - character(len=*), intent(in) :: what - integer, intent(in) :: myID - - integer, dimension(:,:), pointer :: lookupMap - integer :: lower,upper,center - - mesh_FEasCP = 0 - select case(IO_lc(what(1:4))) - case('elem') - lookupMap => mesh_mapFEtoCPelem - case('node') - lookupMap => mesh_mapFEtoCPnode - case default - return - endselect - - lower = 1 - upper = int(size(lookupMap,2),pInt) - - if (lookupMap(1,lower) == myID) then ! check at bounds QUESTION is it valid to extend bounds by 1 and just do binary search w/o init check at bounds? - mesh_FEasCP = lookupMap(2,lower) - return - elseif (lookupMap(1,upper) == myID) then - mesh_FEasCP = lookupMap(2,upper) - return - endif - binarySearch: do while (upper-lower > 1) - center = (lower+upper)/2 - if (lookupMap(1,center) < myID) then - lower = center - elseif (lookupMap(1,center) > myID) then - upper = center - else - mesh_FEasCP = lookupMap(2,center) - exit - endif - enddo binarySearch - -end function mesh_FEasCP - -end module mesh diff --git a/src/prec.f90 b/src/prec.f90 index 7ab63d8f9..8d8e63e09 100644 --- a/src/prec.f90 +++ b/src/prec.f90 @@ -12,11 +12,7 @@ module prec implicit none public ! https://software.intel.com/en-us/blogs/2017/03/27/doctor-fortran-in-it-takes-all-kinds -#ifdef Abaqus - integer, parameter, public :: pReal = selected_real_kind(15,307) !< number with 15 significant digits, up to 1e+-307 (typically 64 bit) -#else integer, parameter, public :: pReal = IEEE_selected_real_kind(15,307) !< number with 15 significant digits, up to 1e+-307 (typically 64 bit) -#endif #if(INT==8) integer, parameter, public :: pInt = selected_int_kind(18) !< number with at least up to +-1e18 (typically 64 bit) #else