From b2fd3e1180155d0c9a4e7a3dd45394c15da8a90b Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 15 Jun 2012 16:10:21 +0000 Subject: [PATCH] introduced preprocessor identifiers Marc, Abaqus, and Spectral to enable conditional compilation. This allows deleted dummy functions that are used by one solver only. Mainly affected modules are IO and mesh. Most of the changes in mesh result from reordering the functions when grouping them depending on their solver. Further advantage is that FE solver do not need FFTW and kdtree2 anymore. The include files for these two libraries moved to DAMASKROO/lib now as I figured out how to use a include path in the Makefile. Put all the files I got when testing compilation with abaqus in a folder which to become the abaqus compilation test. --- code/CPFEM.f90 | 78 +- code/DAMASK_abaqus_exp.f | 49 +- code/DAMASK_abaqus_std.f | 39 +- code/DAMASK_marc.f90 | 26 +- code/DAMASK_spectral.f90 | 24 +- code/DAMASK_spectral_interface.f90 | 494 ++- code/FEsolving.f90 | 151 +- code/IO.f90 | 301 +- code/Makefile | 17 +- code/math.f90 | 12 +- code/mesh.f90 | 5540 ++++++++++++++-------------- code/numerics.f90 | 109 +- {code => lib}/fftw3.f03 | 0 {code => lib}/kdtree2.f90 | 0 lib/pathinfo | 1 + 15 files changed, 3459 insertions(+), 3382 deletions(-) rename {code => lib}/fftw3.f03 (100%) rename {code => lib}/kdtree2.f90 (100%) diff --git a/code/CPFEM.f90 b/code/CPFEM.f90 index 5baa98898..aebdd32cf 100644 --- a/code/CPFEM.f90 +++ b/code/CPFEM.f90 @@ -91,7 +91,9 @@ subroutine CPFEM_initAll(Temperature,element,IP) call crystallite_init(Temperature) ! (have to) use temperature of first IP for whole model call homogenization_init(Temperature) call CPFEM_init - if (trim(FEsolver)/='Spectral') call DAMASK_interface_init() ! Spectral solver is doing initialization earlier +#ifndef Spectral + call DAMASK_interface_init() ! Spectral solver is doing initialization earlier +#endif CPFEM_init_done = .true. CPFEM_init_inProgress = .false. else ! loser, loser... @@ -118,7 +120,7 @@ subroutine CPFEM_init use FEsolving, only: parallelExecution, & symmetricSolver, & restartRead, & - FEmodelGeometry + modelName use mesh, only: mesh_NcpElems, & mesh_maxNips use material, only: homogenization_maxNgrains, & @@ -149,31 +151,31 @@ subroutine CPFEM_init !$OMP END CRITICAL (write2out) endif - call IO_read_jobBinaryFile(777,'recordedPhase',FEmodelGeometry,size(material_phase)) + call IO_read_jobBinaryFile(777,'recordedPhase',modelName,size(material_phase)) read (777,rec=1) material_phase close (777) - call IO_read_jobBinaryFile(777,'convergedF',FEmodelGeometry,size(crystallite_F0)) + call IO_read_jobBinaryFile(777,'convergedF',modelName,size(crystallite_F0)) read (777,rec=1) crystallite_F0 close (777) - call IO_read_jobBinaryFile(777,'convergedFp',FEmodelGeometry,size(crystallite_Fp0)) + call IO_read_jobBinaryFile(777,'convergedFp',modelName,size(crystallite_Fp0)) read (777,rec=1) crystallite_Fp0 close (777) - call IO_read_jobBinaryFile(777,'convergedLp',FEmodelGeometry,size(crystallite_Lp0)) + call IO_read_jobBinaryFile(777,'convergedLp',modelName,size(crystallite_Lp0)) read (777,rec=1) crystallite_Lp0 close (777) - call IO_read_jobBinaryFile(777,'convergeddPdF',FEmodelGeometry,size(crystallite_dPdF0)) + call IO_read_jobBinaryFile(777,'convergeddPdF',modelName,size(crystallite_dPdF0)) read (777,rec=1) crystallite_dPdF0 close (777) - call IO_read_jobBinaryFile(777,'convergedTstar',FEmodelGeometry,size(crystallite_Tstar0_v)) + call IO_read_jobBinaryFile(777,'convergedTstar',modelName,size(crystallite_Tstar0_v)) read (777,rec=1) crystallite_Tstar0_v close (777) - call IO_read_jobBinaryFile(777,'convergedStateConst',FEmodelGeometry) + call IO_read_jobBinaryFile(777,'convergedStateConst',modelName) m = 0_pInt do i = 1,homogenization_maxNgrains; do j = 1,mesh_maxNips; do k = 1,mesh_NcpElems do l = 1,size(constitutive_state0(i,j,k)%p) @@ -183,7 +185,7 @@ subroutine CPFEM_init enddo; enddo; enddo close (777) - call IO_read_jobBinaryFile(777,'convergedStateHomog',FEmodelGeometry) + call IO_read_jobBinaryFile(777,'convergedStateHomog',modelName) m = 0_pInt do k = 1,mesh_NcpElems; do j = 1,mesh_maxNips do l = 1,homogenization_sizeState(j,k) @@ -193,7 +195,7 @@ subroutine CPFEM_init enddo; enddo close (777) - call IO_read_jobBinaryFile(777,'convergeddcsdE',FEmodelGeometry,size(CPFEM_dcsdE)) + call IO_read_jobBinaryFile(777,'convergeddcsdE',modelName,size(CPFEM_dcsdE)) read (777,rec=1) CPFEM_dcsdE close (777) restartRead = .false. @@ -346,15 +348,10 @@ subroutine CPFEM_general(mode, coords, ffn, ffn1, Temperature, dt, element, IP, H, & jacobian3333 ! jacobian in Matrix notation integer(pInt) cp_en, & ! crystal plasticity element number - i, & - j, & - k, & - l, & - m, & - n, & - e, & - node, & - FEnodeID + i, j, k, l, m, n, e +#ifdef Marc + integer(pInt):: node, FEnodeID +#endif logical updateJaco ! flag indicating if JAcobian has to be updated @@ -520,8 +517,8 @@ subroutine CPFEM_general(mode, coords, ffn, ffn1, Temperature, dt, element, IP, write(6,'(a,i8,1x,i2)') '<< CPFEM >> Calculation for element ip ',cp_en,IP !$OMP END CRITICAL (write2out) endif - call materialpoint_stressAndItsTangent(updateJaco, dt) ! calculate stress and its tangent - call materialpoint_postResults(dt) ! post results + call materialpoint_stressAndItsTangent(updateJaco, dt) ! calculate stress and its tangent + call materialpoint_postResults(dt) ! post results !* parallel computation and calulation not yet done @@ -531,21 +528,22 @@ subroutine CPFEM_general(mode, coords, ffn, ffn1, Temperature, dt, element, IP, write(6,'(a,i8,a,i8)') '<< CPFEM >> Calculation for elements ',FEsolving_execElem(1),' to ',FEsolving_execElem(2) !$OMP END CRITICAL (write2out) endif - if (FEsolver == 'Marc') then ! marc returns nodal coordinates, whereas Abaqus and spectral solver return ip coordinates. So for marc we have to calculate the ip coordinates from the nodal coordinates. - call mesh_build_subNodeCoords() ! update subnodal coordinates - call mesh_build_ipCoordinates() ! update ip coordinates - endif +#ifdef Marc +! marc returns nodal coordinates, whereas Abaqus and spectral solver return ip coordinates. So for marc we have to calculate the ip coordinates from the nodal coordinates. + call mesh_build_subNodeCoords() ! update subnodal coordinates + call mesh_build_ipCoordinates() ! update ip coordinates +#endif if (iand(debug_what(debug_CPFEM), debug_levelBasic) /= 0_pInt) then !$OMP CRITICAL (write2out) write(6,'(a,i8,a,i8)') '<< CPFEM >> Start stress and tangent ',FEsolving_execElem(1),' to ',FEsolving_execElem(2) !$OMP END CRITICAL (write2out) endif - call materialpoint_stressAndItsTangent(updateJaco, dt) ! calculate stress and its tangent (parallel execution inside) - call materialpoint_postResults(dt) ! post results + call materialpoint_stressAndItsTangent(updateJaco, dt) ! calculate stress and its tangent (parallel execution inside) + call materialpoint_postResults(dt) ! post results !$OMP PARALLEL DO - do e = FEsolving_execElem(1),FEsolving_execElem(2) ! loop over all parallely processed elements - if (microstructure_elemhomo(mesh_element(4,e))) then ! dealing with homogeneous element? - forall (i = 2:FE_Nips(mesh_element(2,e))) ! copy results of first IP to all others + do e = FEsolving_execElem(1),FEsolving_execElem(2) ! loop over all parallely processed elements + if (microstructure_elemhomo(mesh_element(4,e))) then ! dealing with homogeneous element? + forall (i = 2:FE_Nips(mesh_element(2,e))) ! copy results of first IP to all others materialpoint_P(1:3,1:3,i,e) = materialpoint_P(1:3,1:3,1,e) materialpoint_F(1:3,1:3,i,e) = materialpoint_F(1:3,1:3,1,e) materialpoint_dPdF(1:3,1:3,1:3,1:3,i,e) = materialpoint_dPdF(1:3,1:3,1:3,1:3,1,e) @@ -606,16 +604,14 @@ subroutine CPFEM_general(mode, coords, ffn, ffn1, Temperature, dt, element, IP, CPFEM_cs(1:6,IP,cp_en) = rnd * CPFEM_odd_stress CPFEM_dcsde(1:6,1:6,IP,cp_en) = CPFEM_odd_jacobian * math_identity2nd(6) CPFEM_calc_done = .false. - select case (FEsolver) - case ('Abaqus','Spectral') - mesh_ipCenterOfGravity(1:3,IP,cp_en) = coords(1:3,1) - case ('Marc') - do node = 1,FE_Nnodes(mesh_element(2,cp_en)) - FEnodeID = mesh_FEasCP('node',mesh_element(4+node,cp_en)) - mesh_node(1:3,FEnodeID) = mesh_node0(1:3,FEnodeID) + coords(1:3,node) - enddo - end select - +#ifndef Marc + mesh_ipCenterOfGravity(1:3,IP,cp_en) = coords(1:3,1) +#else + do node = 1,FE_Nnodes(mesh_element(2,cp_en)) + FEnodeID = mesh_FEasCP('node',mesh_element(4+node,cp_en)) + mesh_node(1:3,FEnodeID) = mesh_node0(1:3,FEnodeID) + coords(1:3,node) + enddo +#endif ! --+>> RECYCLING OF FORMER RESULTS (MARC SPECIALTY) <<+-- diff --git a/code/DAMASK_abaqus_exp.f b/code/DAMASK_abaqus_exp.f index 84172cd85..5d1f23634 100644 --- a/code/DAMASK_abaqus_exp.f +++ b/code/DAMASK_abaqus_exp.f @@ -38,15 +38,15 @@ !******************************************************************** #include "prec.f90" +#define Abaqus +module DAMASK_interface -MODULE DAMASK_interface - -character(len=64), parameter :: FEsolver = 'Abaqus' +implicit none character(len=4), parameter :: InputFileExtension = '.inp' character(len=4), parameter :: LogFileExtension = '.log' -CONTAINS +contains !-------------------- subroutine DAMASK_interface_init() @@ -55,13 +55,14 @@ subroutine DAMASK_interface_init() write(6,*) '<<<+- DAMASK_abaqus init -+>>>' write(6,*) '$Id$' write(6,*) - return -end subroutine + +end subroutine DAMASK_interface_init !-------------------- function getSolverWorkingDirectoryName() !-------------------- - use prec + use prec, only: pInt + implicit none character(1024) getSolverWorkingDirectoryName integer(pInt) LENOUTDIR @@ -69,32 +70,24 @@ function getSolverWorkingDirectoryName() getSolverWorkingDirectoryName='' CALL VGETOUTDIR( getSolverWorkingDirectoryName, LENOUTDIR ) getSolverWorkingDirectoryName=trim(getSolverWorkingDirectoryName)//'/' -end function - - -!-------------------- -function getModelName() -!-------------------- - character(1024) getModelName - - getModelName = getSolverJobName() -end function - + +end function getSolverWorkingDirectoryName !-------------------- function getSolverJobName() !-------------------- - use prec + use prec, only: pInt + implicit none - character(1024) getSolverJobName, JOBNAME integer(pInt) LENJOBNAME getSolverJobName='' CALL VGETJOBNAME(getSolverJobName , LENJOBNAME ) -end function + +end function getSolverJobName -END MODULE +end module DAMASK_interface #include "IO.f90" #include "numerics.f90" @@ -153,8 +146,7 @@ subroutine vumat (jblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal, & stressNew, stateNew, enerInternNew, enerInelasNew, & jblock(5), jblock(2)) - return - end subroutine + end subroutine vumat subroutine vumatXtrArg (nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal, & @@ -302,18 +294,17 @@ subroutine vumat (jblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal, & enddo - return - end subroutine + end subroutine vumatXtrArg !******************************************************************** ! This subroutine replaces the corresponding Marc subroutine !******************************************************************** subroutine quit(mpie_error) - use prec, only: pReal, & - pInt + use prec, only: pInt + implicit none integer(pInt) mpie_error call xit - end subroutine + end subroutine quit diff --git a/code/DAMASK_abaqus_std.f b/code/DAMASK_abaqus_std.f index 1337edc43..fcb26351e 100644 --- a/code/DAMASK_abaqus_std.f +++ b/code/DAMASK_abaqus_std.f @@ -38,15 +38,16 @@ !******************************************************************** #include "prec.f90" +#define Abaqus -MODULE DAMASK_interface +module DAMASK_interface -character(len=64), parameter :: FEsolver = 'Abaqus' +implicit none character(len=4), parameter :: InputFileExtension = '.inp' character(len=4), parameter :: LogFileExtension = '.log' -CONTAINS +contains !-------------------- subroutine DAMASK_interface_init() @@ -55,13 +56,14 @@ subroutine DAMASK_interface_init() write(6,*) '<<<+- DAMASK_abaqus init -+>>>' write(6,*) '$Id$' write(6,*) - return -end subroutine + +end subroutine DAMASK_interface_init !-------------------- function getSolverWorkingDirectoryName() !-------------------- use prec + implicit none character(1024) getSolverWorkingDirectoryName integer(pInt) LENOUTDIR @@ -70,33 +72,26 @@ function getSolverWorkingDirectoryName() CALL GETOUTDIR( getSolverWorkingDirectoryName, LENOUTDIR ) getSolverWorkingDirectoryName=trim(getSolverWorkingDirectoryName)//'/' ! write(6,*) 'getSolverWorkingDirectoryName', getSolverWorkingDirectoryName -end function - -!-------------------- -function getModelName() -!-------------------- - character(1024) getModelName - - getModelName = getSolverJobName() -end function +end function getSolverWorkingDirectoryName !-------------------- function getSolverJobName() !-------------------- use prec + implicit none - character(1024) getSolverJobName, JOBNAME integer(pInt) LENJOBNAME getSolverJobName='' CALL GETJOBNAME(getSolverJobName , LENJOBNAME ) ! write(6,*) 'getSolverJobName', getSolverJobName -end function -END MODULE +end function getSolverJobName + +end module DAMASK_interface #include "IO.f90" #include "numerics.f90" @@ -150,7 +145,6 @@ subroutine UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,& implicit none - CHARACTER*80 CMNAME integer(pInt) ndi, nshr, ntens, nstatv, nprops, noel, npt,& kslay, kspt, kstep, kinc @@ -290,18 +284,17 @@ subroutine UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,& if ( terminallyIll ) pnewdt = 0.5_pReal ! force cutback directly ? - return - end subroutine + end subroutine UMAT !******************************************************************** ! This subroutine replaces the corresponding Marc subroutine !******************************************************************** subroutine quit(mpie_error) - use prec, only: pReal, & - pInt + use prec, only: pInt + implicit none integer(pInt) mpie_error call xit - end subroutine + end subroutine quit diff --git a/code/DAMASK_marc.f90 b/code/DAMASK_marc.f90 index c70c5ecd5..ba0709bd5 100644 --- a/code/DAMASK_marc.f90 +++ b/code/DAMASK_marc.f90 @@ -56,13 +56,14 @@ !******************************************************************** ! #include "prec.f90" - +#define Marc module DAMASK_interface - - character(len=64), parameter :: FEsolver = 'Marc' - character(len=4), parameter :: InputFileExtension = '.dat' - character(len=4), parameter :: LogFileExtension = '.log' + use prec, only: pInt + + implicit none + character(len=4), parameter :: InputFileExtension = '.dat' + character(len=4), parameter :: LogFileExtension = '.log' contains @@ -95,19 +96,8 @@ function getSolverWorkingDirectoryName() end function getSolverWorkingDirectoryName - -function getModelName() - implicit none - character(1024) :: getModelName - - getModelName = getSolverJobName() - -end function getModelName - - function getSolverJobName() - use prec, only: pInt implicit none character(1024) :: getSolverJobName, inputName @@ -116,7 +106,7 @@ function getSolverJobName() getSolverJobName='' inputName='' - inquire(5, name=inputName) ! determine outputfile + inquire(5, name=inputName) ! determine inputfile extPos = len_trim(inputName)-4 getSolverJobName=inputName(scan(inputName,pathSep,back=.true.)+1:extPos) ! write(6,*) 'getSolverJobName', getSolverJobName @@ -241,7 +231,7 @@ subroutine hypela2(& !$ use numerics, only: DAMASK_NumThreadsInt ! number of threads set by DAMASK_NUM_THREADS implicit none - include "omp_lib.h" ! the openMP function library +!$ include "omp_lib.h" ! the openMP function library ! ** Start of generated type statements ** real(pReal) coord, d, de, disp, dispt, dt, e, eigvn, eigvn1, ffn, ffn1 real(pReal) frotn, frotn1, g diff --git a/code/DAMASK_spectral.f90 b/code/DAMASK_spectral.f90 index 4c72aa6c8..e0cd3a4dd 100644 --- a/code/DAMASK_spectral.f90 +++ b/code/DAMASK_spectral.f90 @@ -35,17 +35,19 @@ ! MPI fuer Eisenforschung, Duesseldorf #include "spectral_quit.f90" +!#ifdef PETSC +!#include "finclude/petscdef.h" +!#endif program DAMASK_spectral use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran >4.6 at the moment) use DAMASK_interface, only: & DAMASK_interface_init, & - getLoadcaseName, & + loadCaseFile, & + geometryFile, & getSolverWorkingDirectoryName, & - getSolverJobName, & - getModelName, & - inputFileExtension + getSolverJobName use prec, only: & pInt, & @@ -106,7 +108,7 @@ program DAMASK_spectral use homogenization, only: & materialpoint_sizeResults, & materialpoint_results - + implicit none !-------------------------------------------------------------------------------------------------- ! variables related to information from load case and geom file @@ -257,7 +259,7 @@ program DAMASK_spectral !-------------------------------------------------------------------------------------------------- ! reading the load case file and allocate data structure containing load cases - call IO_open_file(myUnit,trim(getLoadcaseName())) + call IO_open_file(myUnit,trim(loadCaseFile)) rewind(myUnit) do read(myUnit,'(a1024)',END = 100) line @@ -279,7 +281,7 @@ program DAMASK_spectral 100 N_Loadcases = N_n if ((N_l + N_Fdot /= N_n) .or. (N_n /= N_t)) & ! sanity check - call IO_error(error_ID=837_pInt,ext_msg = trim(getLoadcaseName())) ! error message for incomplete loadcase + call IO_error(error_ID=837_pInt,ext_msg = trim(loadCaseFile)) ! error message for incomplete loadcase allocate (bc(N_Loadcases)) !-------------------------------------------------------------------------------------------------- @@ -374,13 +376,13 @@ program DAMASK_spectral write(6,'(a)') 'The spectral method boundary value problem solver for' write(6,'(a)') 'the Duesseldorf Advanced Material Simulation Kit' write(6,'(a)') '#############################################################' - write(6,'(a)') 'geometry file: ',trim(getModelName())//InputFileExtension + write(6,'(a)') 'geometry file: ',trim(geometryFile) write(6,'(a)') '=============================================================' write(6,'(a,3(i12 ))') 'resolution a b c:', res write(6,'(a,3(f12.5))') 'dimension x y z:', geomdim write(6,'(a,i5)') 'homogenization: ',homog write(6,'(a)') '#############################################################' - write(6,'(a)') 'loadcase file: ',trim(getLoadcaseName()) + write(6,'(a)') 'loadcase file: ',trim(loadCaseFile) !-------------------------------------------------------------------------------------------------- ! consistency checks and output of load case @@ -580,9 +582,9 @@ C_ref = C * wgt ! write header of output file open(538,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())& //'.spectralOut',form='UNFORMATTED',status='REPLACE') - write(538) 'load', trim(getLoadcaseName()) + write(538) 'load', trim(loadCaseFile) write(538) 'workingdir', trim(getSolverWorkingDirectoryName()) - write(538) 'geometry', trim(getSolverJobName())//InputFileExtension + write(538) 'geometry', trim(geometryFile) write(538) 'resolution', res write(538) 'dimension', geomdim write(538) 'materialpoint_sizeResults', materialpoint_sizeResults diff --git a/code/DAMASK_spectral_interface.f90 b/code/DAMASK_spectral_interface.f90 index f0337447c..18e12c458 100644 --- a/code/DAMASK_spectral_interface.f90 +++ b/code/DAMASK_spectral_interface.f90 @@ -24,185 +24,191 @@ !! by DAMASK !-------------------------------------------------------------------------------------------------- module DAMASK_interface + use prec, only: & + pInt implicit none private - character(len=64), parameter, public :: FEsolver = 'Spectral' !< Keyword for spectral solver - character(len=5), parameter, public :: inputFileExtension = '.geom' !< File extension for geometry description - character(len=4), parameter, public :: logFileExtension = '.log' !< Dummy variable as the spectral solver has no log - character(len=1024), private :: geometryParameter, & !< Interpretated parameter given at command line - loadcaseParameter !< Interpretated parameter given at command line + logical, public :: & + appendToOutFile = .false. !< Append to existing spectralOut file (in case of restart, not in case of regridding) + integer(pInt), public :: & + spectralRestart = 1_pInt !< Increment at which calculation starts + character(len=1024), public :: & + geometryFile = '', & !< parameter given for geometry file + loadCaseFile = '' !< parameter given for load case file - public :: getSolverWorkingDirectoryName, & !< Interpretated parameter given at command line + public :: getSolverWorkingDirectoryName, & getSolverJobName, & - getLoadCase, & - getLoadCaseName, & - getModelName, & DAMASK_interface_init - private :: rectifyPath, & + private :: getGeometryFile, & + getLoadCaseFile, & + rectifyPath, & makeRelativePath, & - getPathSep + getPathSep, & + IO_stringValue, & + IO_intValue, & + IO_lc, & + IO_stringPos contains !-------------------------------------------------------------------------------------------------- !> @brief initializes the solver by interpreting the command line arguments. Also writes -!! information on computation on screen +!! information on computation to screen !-------------------------------------------------------------------------------------------------- -subroutine DAMASK_interface_init(loadcaseParameterIn,geometryParameterIn) +subroutine DAMASK_interface_init(loadCaseParameterIn,geometryParameterIn) use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) - use prec, only: pInt implicit none character(len=1024), optional, intent(in) :: & - loadcaseParameterIn, & + loadCaseParameterIn, & geometryParameterIn character(len=1024) :: & commandLine, & !< command line call as string + geometryParameter, & + loadCaseParameter, & hostName, & !< name of computer - userName !< name of user calling the executable + userName, & !< name of user calling the executable + tag integer :: & - i, & - start ,& - length + i + integer, parameter :: & + maxNchunks = 7 + integer, dimension(1+ 2* maxNchunks) :: & + positions integer, dimension(8) :: & dateAndTime ! type default integer + logical :: & + gotLoadCase = .false., & + gotGeometry = .false., & + gotRestart = .false. + write(6,'(a)') '' write(6,'(a)') '<<<+- DAMASK_spectral_interface init -+>>>' write(6,'(a)') '$Id$' #include "compilation_info.f90" + if ( present(loadcaseParameterIn) .and. present(geometryParameterIn)) then ! both mandatory parameters given in function call geometryParameter = geometryParameterIn loadcaseParameter = loadcaseParameterIn commandLine='n/a' - start = 3_pInt - else if ( .not.( present(loadcaseParameterIn) .and. present(geometryParameterIn))) then ! none parameters given in function call, trying to get them from comman line + gotLoadCase = .true. + gotGeometry = .true. + else if ( .not.( present(loadcaseParameterIn) .and. present(geometryParameterIn))) then ! none parameters given in function call, trying to get them from command line call get_command(commandLine) - call date_and_time(values = dateAndTime) - do i = 1,len(commandLine) ! remove capitals - if(64 0 .or. index(commandLine,' --help ',.true.) > 0) then ! search for ' -h ' or '--help' - write(6,'(a)') '#############################################################' - write(6,'(a)') 'DAMASK spectral:' - write(6,'(a)') 'The spectral method boundary value problem solver for' - write(6,'(a)') 'the Duesseldorf Advanced Material Simulation Kit' - write(6,'(a)') '#############################################################' - write(6,'(a)') 'Valid command line switches:' - write(6,'(a)') ' --geom (-g, --geometry)' - write(6,'(a)') ' --load (-l, --loadcase)' - write(6,'(a)') ' --restart (-r)' - write(6,'(a)') ' --help (-h)' - write(6,'(a)') ' ' - write(6,'(a)') 'Mandatory Arguments:' - write(6,'(a)') ' --load PathToLoadFile/NameOfLoadFile.load' - write(6,'(a)') ' "PathToGeomFile" will be the working directory.' - write(6,'(a)') ' Make sure the file "material.config" exists in the working' - write(6,'(a)') ' directory' - write(6,'(a)') ' For further configuration place "numerics.config"' - write(6,'(a)') ' and "numerics.config" in that directory.' - write(6,'(a)') ' ' - write(6,'(a)') ' --geom PathToGeomFile/NameOfGeom.geom' - write(6,'(a)') ' ' - write(6,'(a)') 'Optional Argument:' - write(6,'(a)') ' --restart XX' - write(6,'(a)') ' Reads in total increment No. XX-1 and continous to' - write(6,'(a)') ' calculate total increment No. XX.' - write(6,'(a)') ' Attention: Overwrites existing results file ' - write(6,'(a)') ' "NameOfGeom_NameOfLoadFile_spectralOut".' - write(6,'(a)') ' Works only if the restart information for total increment' - write(6,'(a)') ' No. XX-1 is available in the working directory.' - write(6,'(a)') 'Help:' - write(6,'(a)') ' --help' - write(6,'(a)') ' Prints this message and exits' - write(6,'(a)') ' ' - call quit(1_pInt) ! normal Termination - endif - if (.not.(command_argument_count()==4 .or. command_argument_count()==6)) then ! check for correct number of given arguments (no --help) - write(6,'(a)') 'Wrong Nr. of Arguments. Run DAMASK_spectral.exe --help' ! Could not find valid keyword (position 0 +3). Functions from IO.f90 are not available - call quit(100_pInt) ! abnormal termination - endif - start = index(commandLine,'-g',.true.) + 3 ! search for '-g' and jump to first char of geometry - if (index(commandLine,'--geom',.true.)>0) then ! if '--geom' is found, use that (contains '-g') - start = index(commandLine,'--geom',.true.) + 7 - endif - if (index(commandLine,'--geometry',.true.)>0) then ! again, now searching for --geometry' - start = index(commandLine,'--geometry',.true.) + 11 - endif - if(start==3_pInt) then ! Could not find valid keyword (position 0 +3). Functions from IO.f90 are not available - write(6,'(a)') 'No Geometry specified' - call quit(100_pInt) ! abnormal termination - endif - length = index(commandLine(start:len(commandLine)),' ',.false.) - - call get_command(commandLine) ! may contain capitals - geometryParameter = '' ! should be empty - geometryParameter(1:length)=commandLine(start:start+length) - - do i=1,len(commandLine) ! remove capitals - if(640) then ! if '--load' is found, use that (contains '-l') - start = index(commandLine,'--load',.true.) + 7 - endif - if (index(commandLine,'--loadcase',.true.)>0) then ! again, now searching for --loadcase' - start = index(commandLine,'--loadcase',.true.) + 11 - endif - if(start==3_pInt) then ! Could not find valid keyword (position 0 +3). Functions from IO.f90 are not available - write(6,'(a)') 'No Loadcase specified' - call quit(100_pInt) ! abnormal termination - endif - length = index(commandLine(start:len(commandLine)),' ',.false.) - - call get_command(commandLine) ! may contain capitals - loadcaseParameter = '' ! should be empty - loadcaseParameter(1:length)=commandLine(start:start+length) - - do i=1,len(commandLine) ! remove capitals - if(640) then ! if '--restart' is found, use that (contains '-l') - start = index(commandLine,'--restart',.true.) + 7 - endif - length = index(commandLine(start:len(commandLine)),' ',.false.) - - call get_command(commandLine) ! may contain capitals - else - write(6,'(a)') 'Wrong Nr. of Arguments!' ! Function call with wrong No of arguments - call quit(100_pInt) + endif + + if (.not. (gotLoadCase .and. gotGeometry)) then + write(6,'(a)') 'Please specify Geometry AND Load Case' + call quit(1_pInt) endif - call GET_ENVIRONMENT_VARIABLE('HOST',hostName) - call GET_ENVIRONMENT_VARIABLE('USER',userName) + geometryFile = getGeometryFile(geometryParameter) + loadCaseFile = getLoadCaseFile(loadCaseParameter) - 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) - write(6,'(a,a)') 'Host Name: ', trim(hostName) - write(6,'(a,a)') 'User Name: ', trim(userName) - write(6,'(a,a)') 'Path Separator: ', getPathSep() - write(6,'(a,a)') 'Command line call: ', trim(commandLine) - write(6,'(a,a)') 'Geometry Parameter: ', trim(geometryParameter) - write(6,'(a,a)') 'Loadcase Parameter: ', trim(loadcaseParameter) - if (start/=3_pInt) write(6,*) 'Restart Parameter: ', trim(commandLine(start:start+length)) + call get_environment_variable('HOST',hostName) + call get_environment_variable('USER',userName) + 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) + write(6,'(a,a)') 'Host name: ', trim(hostName) + write(6,'(a,a)') 'User name: ', trim(userName) + write(6,'(a,a)') 'Path separator: ', getPathSep() + write(6,'(a,a)') 'Command line call: ', trim(commandLine) + write(6,'(a,a)') 'Geometry parameter: ', trim(geometryParameter) + write(6,'(a,a)') 'Loadcase parameter: ', trim(loadcaseParameter) + if (SpectralRestart > 1) write(6,'(a,i6.6)') & + 'Restart at increment: ', spectralRestart + write(6,'(a,l1)') 'Append to result file: ', appendToOutFile end subroutine DAMASK_interface_init !-------------------------------------------------------------------------------------------------- !> @brief extract working directory from loadcase file possibly based on current working dir !-------------------------------------------------------------------------------------------------- - character(len=1024) function getSolverWorkingDirectoryName() +character(len=1024) function getSolverWorkingDirectoryName() implicit none character(len=1024) :: cwd @@ -210,11 +216,12 @@ end subroutine DAMASK_interface_init pathSep = getPathSep() - if (geometryParameter(1:1) == pathSep) then ! absolute path given as command line argument - getSolverWorkingDirectoryName = geometryParameter(1:scan(geometryParameter,pathSep,back=.true.)) + if (geometryFile(1:1) == pathSep) then ! absolute path given as command line argument + getSolverWorkingDirectoryName = geometryFile(1:scan(geometryFile,pathSep,back=.true.)) else - call getcwd(cwd) - getSolverWorkingDirectoryName = trim(cwd)//pathSep//geometryParameter(1:scan(geometryParameter,pathSep,back=.true.)) + call getcwd(cwd) ! relative path given as command line argument + getSolverWorkingDirectoryName = trim(cwd)//pathSep//& + geometryFile(1:scan(geometryFile,pathSep,back=.true.)) endif getSolverWorkingDirectoryName = rectifyPath(getSolverWorkingDirectoryName) @@ -223,93 +230,92 @@ end function getSolverWorkingDirectoryName !-------------------------------------------------------------------------------------------------- -!> @brief basename of geometry file from command line arguments +!> @brief solver job name (no extension) as combination of geometry and load case name !-------------------------------------------------------------------------------------------------- character(len=1024) function getSolverJobName() implicit none - getSolverJobName = trim(getModelName())//'_'//trim(getLoadCase()) + integer :: posExt,posSep + character :: pathSep + character(len=1024) :: tempString + + pathSep = getPathSep() + + tempString = geometryFile + posExt = scan(tempString,'.',back=.true.) + posSep = scan(tempString,pathSep,back=.true.) + + getSolverJobName = tempString(posSep+1:posExt-1) + + tempString = loadCaseFile + posExt = scan(tempString,'.',back=.true.) + posSep = scan(tempString,pathSep,back=.true.) + + getSolverJobName = trim(getSolverJobName)//'_'//tempString(posSep+1:posExt-1) end function getSolverJobName !-------------------------------------------------------------------------------------------------- -!> @brief basename of geometry file from command line arguments +!> @brief basename of geometry file with extension from command line arguments !-------------------------------------------------------------------------------------------------- -character(len=1024) function getModelName() - - use prec, only: pInt +character(len=1024) function getGeometryFile(geometryParameter) implicit none - character(len=1024) :: cwd - integer :: posExt,posSep + character(len=1024), intent(in) :: & + geometryParameter + character(len=1024) :: & + cwd + integer :: posExt, posSep character :: pathSep + getGeometryFile = geometryParameter pathSep = getPathSep() - posExt = scan(geometryParameter,'.',back=.true.) - posSep = scan(geometryParameter,pathSep,back=.true.) + posExt = scan(getGeometryFile,'.',back=.true.) + posSep = scan(getGeometryFile,pathSep,back=.true.) - if (posExt <= posSep) posExt = len_trim(geometryParameter)+1 ! no extension present - getModelName = geometryParameter(1:posExt-1_pInt) ! path to geometry file (excl. extension) - - if (scan(getModelName,pathSep) /= 1) then ! relative path given as command line argument + if (posExt <= posSep) getGeometryFile = trim(getGeometryFile)//('.geom') ! no extension present + if (scan(getGeometryFile,pathSep) /= 1) then ! relative path given as command line argument call getcwd(cwd) - getModelName = rectifyPath(trim(cwd)//'/'//getModelName) + getGeometryFile = rectifyPath(trim(cwd)//pathSep//getGeometryFile) else - getModelName = rectifyPath(getModelName) + getGeometryFile = rectifyPath(getGeometryFile) endif - getModelName = makeRelativePath(getSolverWorkingDirectoryName(),& - getModelName) - -end function getModelName + getGeometryFile = makeRelativePath(getSolverWorkingDirectoryName(), getGeometryFile) - -!-------------------------------------------------------------------------------------------------- -!> @brief name of load case file exluding extension -!-------------------------------------------------------------------------------------------------- -character(len=1024) function getLoadCase() - - implicit none - integer :: posExt,posSep - character :: pathSep - - pathSep = getPathSep() - posExt = scan(loadcaseParameter,'.',back=.true.) - posSep = scan(loadcaseParameter,pathSep,back=.true.) - - if (posExt <= posSep) posExt = len_trim(loadcaseParameter)+1 ! no extension present - getLoadCase = loadcaseParameter(posSep+1:posExt-1) ! name of load case file exluding extension - -end function getLoadCase +end function getGeometryFile !-------------------------------------------------------------------------------------------------- !> @brief relative path of loadcase from command line arguments !-------------------------------------------------------------------------------------------------- -character(len=1024) function getLoadcaseName() +character(len=1024) function getLoadCaseFile(loadCaseParameter) implicit none - character(len=1024) :: cwd - integer :: posExt = 0, posSep + character(len=1024), intent(in) :: & + loadCaseParameter + character(len=1024) :: & + cwd + integer :: posExt, posSep character :: pathSep - + + getLoadCaseFile = loadcaseParameter pathSep = getPathSep() - getLoadcaseName = loadcaseParameter - posExt = scan(getLoadcaseName,'.',back=.true.) - posSep = scan(getLoadcaseName,pathSep,back=.true.) + posExt = scan(getLoadCaseFile,'.',back=.true.) + posSep = scan(getLoadCaseFile,pathSep,back=.true.) - if (posExt <= posSep) getLoadcaseName = trim(getLoadcaseName)//('.load') ! no extension present - if (scan(getLoadcaseName,pathSep) /= 1) then ! relative path given as command line argument + if (posExt <= posSep) getLoadCaseFile = trim(getLoadCaseFile)//('.load') ! no extension present + if (scan(getLoadCaseFile,pathSep) /= 1) then ! relative path given as command line argument call getcwd(cwd) - getLoadcaseName = rectifyPath(trim(cwd)//pathSep//getLoadcaseName) + getLoadCaseFile = rectifyPath(trim(cwd)//pathSep//getLoadCaseFile) else - getLoadcaseName = rectifyPath(getLoadcaseName) + getLoadCaseFile = rectifyPath(getLoadCaseFile) endif - getLoadcaseName = makeRelativePath(getSolverWorkingDirectoryName(),& - getLoadcaseName) -end function getLoadcaseName + getLoadCaseFile = makeRelativePath(getSolverWorkingDirectoryName(), getLoadCaseFile) + +end function getLoadCaseFile !-------------------------------------------------------------------------------------------------- @@ -384,8 +390,6 @@ end function makeRelativePath !-------------------------------------------------------------------------------------------------- character function getPathSep() - use prec, only: pInt - implicit none character(len=2048) path integer(pInt) :: backslash = 0_pInt, slash = 0_pInt @@ -403,6 +407,106 @@ character function getPathSep() getPathSep = '/' endif -end function +end function getPathSep + +!******************************************************************** +! read string value at myPos from line +!******************************************************************** + pure function IO_stringValue(line,positions,myPos) + + implicit none + + integer(pInt), intent(in) :: positions(*), & + myPos + + character(len=1+positions(myPos*2+1)-positions(myPos*2)) :: IO_stringValue + + character(len=*), intent(in) :: line + + if (positions(1) < myPos) then + IO_stringValue = '' + else + IO_stringValue = line(positions(myPos*2):positions(myPos*2+1)) + endif + +end function IO_stringValue + +!******************************************************************** +! read int value at myPos from line +!******************************************************************** +integer(pInt) pure function IO_intValue(line,positions,myPos) + + implicit none + character(len=*), intent(in) :: line + integer(pInt), intent(in) :: positions(*), & + myPos + + if (positions(1) < myPos) then + IO_intValue = 0_pInt + else + read(UNIT=line(positions(myPos*2):positions(myPos*2+1)),ERR=100,FMT=*) IO_intValue + endif + return +100 IO_intValue = huge(1_pInt) + +end function IO_intValue + +!******************************************************************** +! change character in line to lower case +!******************************************************************** +pure function IO_lc(line) + + implicit none + character(26), parameter :: lower = 'abcdefghijklmnopqrstuvwxyz' + character(26), parameter :: upper = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ' + character(len=*), intent(in) :: line + character(len=len(line)) :: IO_lc + + integer :: i,n ! no pInt (len returns default integer) + + IO_lc = line + do i=1,len(line) + n = index(upper,IO_lc(i:i)) + if (n/=0) IO_lc(i:i) = lower(n:n) + enddo + +end function IO_lc + +!******************************************************************** +! locate at most N space-separated parts in line +! return array containing number of parts in line and +! the left/right positions of at most N to be used by IO_xxxVal +!******************************************************************** +pure function IO_stringPos(line,N) + + implicit none + integer(pInt), intent(in) :: N + integer(pInt) :: IO_stringPos(1_pInt+N*2_pInt) + + character(len=*), intent(in) :: line + + character(len=*), parameter :: sep=achar(44)//achar(32)//achar(9)//achar(10)//achar(13) ! comma and whitespaces + + integer :: left, right !no pInt (verify and scan return default integer) + + + IO_stringPos = -1_pInt + IO_stringPos(1) = 0_pInt + right = 0 + + do while (verify(line(right+1:),sep)>0) + left = right + verify(line(right+1:),sep) + right = left + scan(line(left:),sep) - 2 + if ( line(left:left) == '#' ) then + exit + endif + if ( IO_stringPos(1)0) & ! look for -r - start = index(commandLine,'-r ',.true.) + 3 ! set to position after trailing space - if (index(commandLine,'--restart ',.true.)>0) & ! look for --restart - start = index(commandLine,'--restart ',.true.) + 10 ! set to position after trailing space - if(start /= 0) then ! found something - length = verify(commandLine(start:len(commandLine)),'0123456789',.false.) ! where is first non number after argument? - read(commandLine(start:start+length),'(I12)') restartInc ! read argument - restartRead = restartInc > 0_pInt - if(restartInc <= 0_pInt) then - call IO_warning(warning_ID=34_pInt) - restartInc = 1_pInt - endif - endif - else +#ifdef Spectral + modelName = getSolverJobName() + restartInc = spectralRestart + if(restartInc <= 0_pInt) then + call IO_warning(warning_ID=34_pInt) + restartInc = 1_pInt + endif +#else + call IO_open_inputFile(fileunit,getSolverJobName()) + rewind(fileunit) + do + read (fileunit,'(a1024)',END=100) line + positions = IO_stringPos(line,maxNchunks) + tag = IO_lc(IO_stringValue(line,positions,1_pInt)) ! extract key + select case(tag) + case ('solver') + read (fileunit,'(a1024)',END=100) line ! next line + positions = IO_stringPos(line,maxNchunks) + symmetricSolver = (IO_intValue(line,positions,2_pInt) /= 1_pInt) + case ('restart') + read (fileunit,'(a1024)',END=100) line ! next line + positions = IO_stringPos(line,maxNchunks) + restartWrite = iand(IO_intValue(line,positions,1_pInt),1_pInt) > 0_pInt + restartRead = iand(IO_intValue(line,positions,1_pInt),2_pInt) > 0_pInt + case ('*restart') + do j=2_pInt,positions(1) + restartWrite = (IO_lc(IO_StringValue(line,positions,j)) == 'write') .or. restartWrite + restartRead = (IO_lc(IO_StringValue(line,positions,j)) == 'read') .or. restartRead + enddo + if(restartWrite) then + do j=2_pInt,positions(1) + restartWrite = (IO_lc(IO_StringValue(line,positions,j)) /= 'frequency=0') .and. restartWrite + enddo + endif + end select + enddo + 100 close(fileunit) + + if (restartRead) then +#ifdef Marc + call IO_open_logFile(fileunit) rewind(fileunit) do - read (fileunit,'(a1024)',END=100) line + read (fileunit,'(a1024)',END=200) line positions = IO_stringPos(line,maxNchunks) - tag = IO_lc(IO_stringValue(line,positions,1_pInt)) ! extract key - select case(tag) - case ('solver') - read (fileunit,'(a1024)',END=100) line ! next line - positions = IO_stringPos(line,maxNchunks) - symmetricSolver = (IO_intValue(line,positions,2_pInt) /= 1_pInt) - case ('restart') - read (fileunit,'(a1024)',END=100) line ! next line - positions = IO_stringPos(line,maxNchunks) - restartWrite = iand(IO_intValue(line,positions,1_pInt),1_pInt) > 0_pInt - restartRead = iand(IO_intValue(line,positions,1_pInt),2_pInt) > 0_pInt - case ('*restart') - do j=2_pInt,positions(1) - restartWrite = (IO_lc(IO_StringValue(line,positions,j)) == 'write') .or. restartWrite - restartRead = (IO_lc(IO_StringValue(line,positions,j)) == 'read') .or. restartRead - enddo - if(restartWrite) then - do j=2_pInt,positions(1) - restartWrite = (IO_lc(IO_StringValue(line,positions,j)) /= 'frequency=0') .and. restartWrite - enddo - endif - end select + if ( IO_lc(IO_stringValue(line,positions,1_pInt)) == 'restart' .and. & + IO_lc(IO_stringValue(line,positions,2_pInt)) == 'file' .and. & + IO_lc(IO_stringValue(line,positions,3_pInt)) == 'job' .and. & + IO_lc(IO_stringValue(line,positions,4_pInt)) == 'id' ) & + modelName = IO_StringValue(line,positions,6_pInt) enddo - endif - 100 close(fileunit) - - if (restartRead) then - if(FEsolver == 'Marc') then - call IO_open_logFile(fileunit) - rewind(fileunit) - do +#else + call IO_open_inputFile(fileunit,modelName) + rewind(fileunit) + do + read (fileunit,'(a1024)',END=200) line + positions = IO_stringPos(line,maxNchunks) + if ( IO_lc(IO_stringValue(line,positions,1_pInt))=='*heading') then read (fileunit,'(a1024)',END=200) line positions = IO_stringPos(line,maxNchunks) - if ( IO_lc(IO_stringValue(line,positions,1_pInt)) == 'restart' .and. & - IO_lc(IO_stringValue(line,positions,2_pInt)) == 'file' .and. & - IO_lc(IO_stringValue(line,positions,3_pInt)) == 'job' .and. & - IO_lc(IO_stringValue(line,positions,4_pInt)) == 'id' ) & - FEmodelGeometry = IO_StringValue(line,positions,6_pInt) - enddo - elseif (FEsolver == 'Abaqus') then - call IO_open_inputFile(fileunit,FEmodelGeometry) - rewind(fileunit) - do - read (fileunit,'(a1024)',END=200) line - positions = IO_stringPos(line,maxNchunks) - if ( IO_lc(IO_stringValue(line,positions,1_pInt))=='*heading') then - read (fileunit,'(a1024)',END=200) line - positions = IO_stringPos(line,maxNchunks) - FEmodelGeometry = IO_StringValue(line,positions,1_pInt) - endif - enddo - endif + modelName = IO_StringValue(line,positions,1_pInt) + endif + enddo +#endif + else + modelName = getSolverJobName() endif 200 close(fileunit) - +#endif !$OMP CRITICAL (write2out) write(6,*) write(6,*) '<<<+- FEsolving init -+>>>' @@ -190,7 +179,7 @@ subroutine FE_init if (iand(debug_what(debug_FEsolving),debug_levelBasic) /= 0_pInt) then write(6,*) 'restart writing: ', restartWrite write(6,*) 'restart reading: ', restartRead - if (restartRead) write(6,*) 'restart Job: ', trim(FEmodelGeometry) + if (restartRead) write(6,*) 'restart Job: ', trim(modelName) write(6,*) endif !$OMP END CRITICAL (write2out) diff --git a/code/IO.f90 b/code/IO.f90 index 9cbe14a7d..1fb3866ea 100644 --- a/code/IO.f90 +++ b/code/IO.f90 @@ -31,12 +31,9 @@ module IO IO_open_jobFile_stat, & IO_open_file, & IO_open_jobFile, & - IO_open_inputFile, & - IO_open_logFile, & IO_write_jobFile, & IO_write_jobBinaryFile, & IO_read_jobBinaryFile, & - IO_abaqus_hasNoPart, & IO_hybridIA, & IO_isBlank, & IO_getTag, & @@ -58,10 +55,20 @@ module IO IO_continuousIntValues, & IO_error, & IO_warning +#ifndef Spectral + public :: IO_open_inputFile, & + IO_open_logFile +#endif +#ifdef Abaqus + public :: IO_abaqus_hasNoPart +#endif + private :: IO_fixedFloatValue, & IO_lcInplace ,& - abaqus_assembleInputFile, & hybridIA_reps +#ifdef Abaqus + private :: abaqus_assembleInputFile +#endif contains @@ -187,7 +194,7 @@ subroutine IO_open_jobFile(myUnit,newExt) end subroutine IO_open_jobFile - +#ifndef Spectral !******************************************************************** ! open FEM inputfile to given myUnit ! AP: 12.07.10 @@ -196,10 +203,10 @@ end subroutine IO_open_jobFile !******************************************************************** subroutine IO_open_inputFile(myUnit,model) - use DAMASK_interface, only: getSolverWorkingDirectoryName,& - getSolverJobName, & - inputFileExtension, & - FEsolver + use DAMASK_interface, only: & + getSolverWorkingDirectoryName,& + getSolverJobName, & + inputFileExtension implicit none integer(pInt), intent(in) :: myUnit @@ -208,24 +215,22 @@ subroutine IO_open_inputFile(myUnit,model) integer(pInt) :: myStat character(len=1024) :: path - if (FEsolver == 'Abaqus') then - - path = trim(getSolverWorkingDirectoryName())//trim(model)//InputFileExtension - open(myUnit+1,status='old',iostat=myStat,file=path) - if (myStat /= 0_pInt) call IO_error(100_pInt,ext_msg=path) +#ifdef Abaqus + path = trim(getSolverWorkingDirectoryName())//trim(model)//InputFileExtension + open(myUnit+1,status='old',iostat=myStat,file=path) + if (myStat /= 0_pInt) call IO_error(100_pInt,ext_msg=path) - path = trim(getSolverWorkingDirectoryName())//trim(model)//InputFileExtension//'_assembly' - open(myUnit,iostat=myStat,file=path) - if (myStat /= 0_pInt) call IO_error(100_pInt,ext_msg=path) - if (.not.abaqus_assembleInputFile(myUnit,myUnit+1_pInt)) call IO_error(103_pInt) ! strip comments and concatenate any "include"s - close(myUnit+1_pInt) - - else + path = trim(getSolverWorkingDirectoryName())//trim(model)//InputFileExtension//'_assembly' + open(myUnit,iostat=myStat,file=path) + if (myStat /= 0_pInt) call IO_error(100_pInt,ext_msg=path) + if (.not.abaqus_assembleInputFile(myUnit,myUnit+1_pInt)) call IO_error(103_pInt) ! strip comments and concatenate any "include"s + close(myUnit+1_pInt) +#endif +#ifdef Marc path = trim(getSolverWorkingDirectoryName())//trim(model)//InputFileExtension open(myUnit,status='old',iostat=myStat,file=path) if (myStat /= 0_pInt) call IO_error(100_pInt,ext_msg=path) - - endif +#endif end subroutine IO_open_inputFile @@ -235,9 +240,10 @@ end subroutine IO_open_inputFile !******************************************************************** subroutine IO_open_logFile(myUnit) - use DAMASK_interface, only: getSolverWorkingDirectoryName, & - getSolverJobName, & - LogFileExtension + use DAMASK_interface, only: & + getSolverWorkingDirectoryName, & + getSolverJobName, & + LogFileExtension implicit none integer(pInt), intent(in) :: myUnit @@ -250,7 +256,7 @@ subroutine IO_open_logFile(myUnit) if (myStat /= 0) call IO_error(100_pInt,ext_msg=path) end subroutine IO_open_logFile - +#endif !******************************************************************** ! open (write) file related to current job @@ -334,7 +340,7 @@ subroutine IO_read_jobBinaryFile(myUnit,newExt,jobName,recMultiplier) end subroutine IO_read_jobBinaryFile - +#ifdef Abaqus !*********************************************************** ! check if the input file for Abaqus contains part info !*********************************************************** @@ -361,7 +367,7 @@ logical function IO_abaqus_hasNoPart(myUnit) enddo 620 end function IO_abaqus_hasNoPart - +#endif !******************************************************************** ! hybrid IA sampling of ODFfile @@ -995,57 +1001,51 @@ end function IO_countDataLines !******************************************************************** integer(pInt) function IO_countContinuousIntValues(myUnit) - use DAMASK_interface, only: FEsolver - implicit none integer(pInt), intent(in) :: myUnit integer(pInt), parameter :: maxNchunks = 8192_pInt - +#ifdef Abaqus integer(pInt) :: l,c +#endif integer(pInt), dimension(1+2*maxNchunks) :: myPos character(len=65536) :: line IO_countContinuousIntValues = 0_pInt - select case (FEsolver) - case ('Marc','Spectral') - - do - read(myUnit,'(A300)',end=100) line - myPos = IO_stringPos(line,maxNchunks) - if (IO_lc(IO_stringValue(line,myPos,2_pInt)) == 'to' ) then ! found range indicator - IO_countContinuousIntValues = 1_pInt + IO_intValue(line,myPos,3_pInt) - IO_intValue(line,myPos,1_pInt) - exit ! only one single range indicator allowed - else if (IO_lc(IO_stringValue(line,myPos,2_pInt)) == 'copies' .and. & - IO_lc(IO_stringValue(line,myPos,3_pInt)) == 'of' ) then ! found multiple entries indicator - IO_countContinuousIntValues = IO_intValue(line,myPos,1_pInt) - exit ! only one single multiplier allowed - else - IO_countContinuousIntValues = IO_countContinuousIntValues+myPos(1)-1_pInt ! add line's count when assuming 'c' - if ( IO_lc(IO_stringValue(line,myPos,myPos(1))) /= 'c' ) then ! line finished, read last value - IO_countContinuousIntValues = IO_countContinuousIntValues+1_pInt - exit ! data ended - endif - endif - enddo - - case('Abaqus') - - c = IO_countDataLines(myUnit) - do l = 1_pInt,c - backspace(myUnit) - enddo +#ifndef Abaqus + do + read(myUnit,'(A300)',end=100) line + myPos = IO_stringPos(line,maxNchunks) + if (IO_lc(IO_stringValue(line,myPos,2_pInt)) == 'to' ) then ! found range indicator + IO_countContinuousIntValues = 1_pInt + IO_intValue(line,myPos,3_pInt) - IO_intValue(line,myPos,1_pInt) + exit ! only one single range indicator allowed + else if (IO_lc(IO_stringValue(line,myPos,2_pInt)) == 'copies' .and. & + IO_lc(IO_stringValue(line,myPos,3_pInt)) == 'of' ) then ! found multiple entries indicator + IO_countContinuousIntValues = IO_intValue(line,myPos,1_pInt) + exit ! only one single multiplier allowed + else + IO_countContinuousIntValues = IO_countContinuousIntValues+myPos(1)-1_pInt ! add line's count when assuming 'c' + if ( IO_lc(IO_stringValue(line,myPos,myPos(1))) /= 'c' ) then ! line finished, read last value + IO_countContinuousIntValues = IO_countContinuousIntValues+1_pInt + exit ! data ended + endif + endif + enddo +#else + c = IO_countDataLines(myUnit) + do l = 1_pInt,c + backspace(myUnit) + enddo - do l = 1_pInt,c - read(myUnit,'(A300)',end=100) line - myPos = IO_stringPos(line,maxNchunks) - IO_countContinuousIntValues = IO_countContinuousIntValues + 1_pInt + & ! assuming range generation - (IO_intValue(line,myPos,2_pInt)-IO_intValue(line,myPos,1_pInt))/& - max(1_pInt,IO_intValue(line,myPos,3_pInt)) - enddo - - end select + do l = 1_pInt,c + read(myUnit,'(A300)',end=100) line + myPos = IO_stringPos(line,maxNchunks) + IO_countContinuousIntValues = IO_countContinuousIntValues + 1_pInt + & ! assuming range generation + (IO_intValue(line,myPos,2_pInt)-IO_intValue(line,myPos,1_pInt))/& + max(1_pInt,IO_intValue(line,myPos,3_pInt)) + enddo +#endif 100 end function IO_countContinuousIntValues @@ -1059,8 +1059,6 @@ integer(pInt) function IO_countContinuousIntValues(myUnit) !******************************************************************** function IO_continuousIntValues(myUnit,maxN,lookupName,lookupMap,lookupMaxN) - use DAMASK_interface, only: FEsolver - implicit none integer(pInt), intent(in) :: maxN integer(pInt), dimension(1+maxN) :: IO_continuousIntValues @@ -1069,10 +1067,12 @@ function IO_continuousIntValues(myUnit,maxN,lookupName,lookupMap,lookupMaxN) lookupMaxN integer(pInt), dimension(:,:), intent(in) :: lookupMap character(len=64), dimension(:), intent(in) :: lookupName - integer(pInt), parameter :: maxNchunks = 8192_pInt - - integer(pInt) :: i,j,l,c,first,last + integer(pInt) :: i +#ifdef Abaqus + integer(pInt) :: j,l,c,first,last +#endif + integer(pInt), dimension(1+2*maxNchunks) :: myPos character(len=65536) line logical rangeGeneration @@ -1080,88 +1080,83 @@ function IO_continuousIntValues(myUnit,maxN,lookupName,lookupMap,lookupMaxN) IO_continuousIntValues = 0_pInt rangeGeneration = .false. - select case (FEsolver) - case ('Marc','Spectral') - - do - read(myUnit,'(A65536)',end=100) line - myPos = IO_stringPos(line,maxNchunks) - if (verify(IO_stringValue(line,myPos,1_pInt),'0123456789') > 0) then ! a non-int, i.e. set name - do i = 1_pInt, lookupMaxN ! loop over known set names - if (IO_stringValue(line,myPos,1_pInt) == lookupName(i)) then ! found matching name - IO_continuousIntValues = lookupMap(:,i) ! return resp. entity list - exit - endif - enddo +#ifndef Abaqus + do + read(myUnit,'(A65536)',end=100) line + myPos = IO_stringPos(line,maxNchunks) + if (verify(IO_stringValue(line,myPos,1_pInt),'0123456789') > 0) then ! a non-int, i.e. set name + do i = 1_pInt, lookupMaxN ! loop over known set names + if (IO_stringValue(line,myPos,1_pInt) == lookupName(i)) then ! found matching name + IO_continuousIntValues = lookupMap(:,i) ! return resp. entity list exit - else if (myPos(1) > 2_pInt .and. IO_lc(IO_stringValue(line,myPos,2_pInt)) == 'to' ) then ! found range indicator - do i = IO_intValue(line,myPos,1_pInt),IO_intValue(line,myPos,3_pInt) - IO_continuousIntValues(1) = IO_continuousIntValues(1) + 1_pInt - IO_continuousIntValues(1+IO_continuousIntValues(1)) = i - enddo - exit - else if (myPos(1) > 3_pInt .and. IO_lc(IO_stringValue(line,myPos,2_pInt)) == 'copies' & - .and. IO_lc(IO_stringValue(line,myPos,3_pInt)) == 'of' ) then ! found multiple entries indicator - IO_continuousIntValues(1) = IO_intValue(line,myPos,1_pInt) - IO_continuousIntValues(2:IO_continuousIntValues(1)+1) = IO_intValue(line,myPos,4_pInt) - exit - else - do i = 1_pInt,myPos(1)-1_pInt ! interpret up to second to last value - IO_continuousIntValues(1) = IO_continuousIntValues(1) + 1_pInt - IO_continuousIntValues(1+IO_continuousIntValues(1)) = IO_intValue(line,myPos,i) - enddo - if ( IO_lc(IO_stringValue(line,myPos,myPos(1))) /= 'c' ) then ! line finished, read last value - IO_continuousIntValues(1) = IO_continuousIntValues(1) + 1_pInt - IO_continuousIntValues(1+IO_continuousIntValues(1)) = IO_intValue(line,myPos,myPos(1)) - exit - endif endif enddo - - case('Abaqus') - - c = IO_countDataLines(myUnit) - do l = 1_pInt,c - backspace(myUnit) + exit + else if (myPos(1) > 2_pInt .and. IO_lc(IO_stringValue(line,myPos,2_pInt)) == 'to' ) then ! found range indicator + do i = IO_intValue(line,myPos,1_pInt),IO_intValue(line,myPos,3_pInt) + IO_continuousIntValues(1) = IO_continuousIntValues(1) + 1_pInt + IO_continuousIntValues(1+IO_continuousIntValues(1)) = i enddo - -! check if the element values in the elset are auto generated - backspace(myUnit) - read(myUnit,'(A65536)',end=100) line - myPos = IO_stringPos(line,maxNchunks) - do i = 1_pInt,myPos(1) - if (IO_lc(IO_stringValue(line,myPos,i)) == 'generate') rangeGeneration = .true. - enddo - - do l = 1_pInt,c - read(myUnit,'(A65536)',end=100) line - myPos = IO_stringPos(line,maxNchunks) - if (verify(IO_stringValue(line,myPos,1_pInt),'0123456789') > 0) then ! a non-int, i.e. set names follow on this line - do i = 1_pInt,myPos(1) ! loop over set names in line - do j = 1_pInt,lookupMaxN ! look thru known set names - if (IO_stringValue(line,myPos,i) == lookupName(j)) then ! found matching name - first = 2_pInt + IO_continuousIntValues(1) ! where to start appending data - last = first + lookupMap(1,j) - 1_pInt ! 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,myPos,1_pInt),IO_intValue(line,myPos,2_pInt),max(1_pInt,IO_intValue(line,myPos,3_pInt)) - IO_continuousIntValues(1) = IO_continuousIntValues(1) + 1_pInt - IO_continuousIntValues(1+IO_continuousIntValues(1)) = i - enddo - else ! read individual elem nums - do i = 1_pInt,myPos(1) -! write(*,*)'IO_CIV-int',IO_intValue(line,myPos,i) - IO_continuousIntValues(1) = IO_continuousIntValues(1) + 1_pInt - IO_continuousIntValues(1+IO_continuousIntValues(1)) = IO_intValue(line,myPos,i) - enddo - endif + exit + else if (myPos(1) > 3_pInt .and. IO_lc(IO_stringValue(line,myPos,2_pInt)) == 'copies' & + .and. IO_lc(IO_stringValue(line,myPos,3_pInt)) == 'of' ) then ! found multiple entries indicator + IO_continuousIntValues(1) = IO_intValue(line,myPos,1_pInt) + IO_continuousIntValues(2:IO_continuousIntValues(1)+1) = IO_intValue(line,myPos,4_pInt) + exit + else + do i = 1_pInt,myPos(1)-1_pInt ! interpret up to second to last value + IO_continuousIntValues(1) = IO_continuousIntValues(1) + 1_pInt + IO_continuousIntValues(1+IO_continuousIntValues(1)) = IO_intValue(line,myPos,i) enddo + if ( IO_lc(IO_stringValue(line,myPos,myPos(1))) /= 'c' ) then ! line finished, read last value + IO_continuousIntValues(1) = IO_continuousIntValues(1) + 1_pInt + IO_continuousIntValues(1+IO_continuousIntValues(1)) = IO_intValue(line,myPos,myPos(1)) + exit + endif + endif + enddo +#else + c = IO_countDataLines(myUnit) + do l = 1_pInt,c + backspace(myUnit) + enddo - endselect + !heck if the element values in the elset are auto generated + backspace(myUnit) + read(myUnit,'(A65536)',end=100) line + myPos = IO_stringPos(line,maxNchunks) + do i = 1_pInt,myPos(1) + if (IO_lc(IO_stringValue(line,myPos,i)) == 'generate') rangeGeneration = .true. + enddo + + do l = 1_pInt,c + read(myUnit,'(A65536)',end=100) line + myPos = IO_stringPos(line,maxNchunks) + if (verify(IO_stringValue(line,myPos,1_pInt),'0123456789') > 0) then ! a non-int, i.e. set names follow on this line + do i = 1_pInt,myPos(1) ! loop over set names in line + do j = 1_pInt,lookupMaxN ! look thru known set names + if (IO_stringValue(line,myPos,i) == lookupName(j)) then ! found matching name + first = 2_pInt + IO_continuousIntValues(1) ! where to start appending data + last = first + lookupMap(1,j) - 1_pInt ! 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,myPos,1_pInt),IO_intValue(line,myPos,2_pInt),max(1_pInt,IO_intValue(line,myPos,3_pInt)) + IO_continuousIntValues(1) = IO_continuousIntValues(1) + 1_pInt + IO_continuousIntValues(1+IO_continuousIntValues(1)) = i + enddo + else ! read individual elem nums + do i = 1_pInt,myPos(1) + ! write(*,*)'IO_CIV-int',IO_intValue(line,myPos,i) + IO_continuousIntValues(1) = IO_continuousIntValues(1) + 1_pInt + IO_continuousIntValues(1+IO_continuousIntValues(1)) = IO_intValue(line,myPos,i) + enddo + endif + enddo +#endif 100 end function IO_continuousIntValues @@ -1426,6 +1421,8 @@ subroutine IO_warning(warning_ID,e,i,g,ext_msg) msg = 'invalid restart increment given' case (35_pInt) msg = 'could not get $DAMASK_NUM_THREADS' + case (40_pInt) + msg = 'Found Spectral solver parameter ' case (47_pInt) msg = 'No valid parameter for FFTW given, using FFTW_PATIENT' case (101_pInt) @@ -1467,8 +1464,10 @@ subroutine IO_warning(warning_ID,e,i,g,ext_msg) end subroutine IO_warning + ! INTERNAL (HELPER) FUNCTIONS: +#ifdef Abaqus !******************************************************************** ! AP: 12.07.10 ! create a new input file for abaqus simulations @@ -1525,7 +1524,7 @@ recursive function abaqus_assembleInputFile(unit1,unit2) result(createSuccess) 200 createSuccess =.false. end function abaqus_assembleInputFile - +#endif !******************************************************************** ! hybrid IA repetition counter diff --git a/code/Makefile b/code/Makefile index e65b8d580..9230348c1 100644 --- a/code/Makefile +++ b/code/Makefile @@ -44,14 +44,19 @@ LAPACKROOT :=/usr F90 ?=ifort COMPILERNAME ?= $(F90) INCLUDE_DIRS +=-I$(DAMASK_ROOT)/lib +ifdef PETSC_DIR +INCLUDE_DIRS +=-I$(PETSC_DIR)/include +INCLUDE_DIRS +=-I$(PETSC_DIR)/$(PETSC_ARCH)/include +INCLUDE_DIRS +=-DPETSC # just for the moment, as long as PETSC is non standard +endif ifeq "$(FASTBUILD)" "YES" OPENMP :=OFF OPTIMIZATION :=OFF -endif - +else OPENMP ?= ON OPTIMIZATION ?= DEFENSIVE +endif ifeq "$(OPTIMIZATION)" "OFF" OPTI := OFF @@ -107,13 +112,13 @@ endif endif ifdef STANDARD_CHECK -STANDARD_CHECK_ifort =$(STANDARD_CHECK) -DSTANDARD_CHECK -STANDARD_CHECK_gfortran =$(STANDARD_CHECK) -DSTANDARD_CHECK +STANDARD_CHECK_ifort =$(STANDARD_CHECK) +STANDARD_CHECK_gfortran =$(STANDARD_CHECK) endif ifneq "$(FASTBUILD)" "YES" -STANDARD_CHECK_ifort ?=-stand f08 -standard-semantics -warn stderrors -DSTANDARD_CHECK -STANDARD_CHECK_gfortran ?=-std=f2008 -fall-intrinsics -DSTANDARD_CHECK +STANDARD_CHECK_ifort ?=-stand f08 -standard-semantics -warn stderrors +STANDARD_CHECK_gfortran ?=-std=f2008 -fall-intrinsics endif #-fall-intrinsics: all intrinsic procedures (including the GNU-specific extensions) are accepted. This can be useful with -std=f95 to force standard-compliance # but get access to the full range of intrinsics available with gfortran. As a consequence, -Wintrinsics-std will be ignored and no user-defined diff --git a/code/math.f90 b/code/math.f90 index 23a6069de..89036429d 100644 --- a/code/math.f90 +++ b/code/math.f90 @@ -19,8 +19,9 @@ !############################################################## !* $Id$ !############################################################## - +#ifdef Spectral #include "kdtree2.f90" +#endif module math !############################################################## @@ -137,8 +138,10 @@ real(pReal), dimension(4,36), parameter, private :: & 0.0_pReal, 0.0_pReal, 0.0_pReal, 1.0_pReal & ],[4,36]) +#ifdef Spectral include 'fftw3.f03' - +#endif + public :: math_init, & qsort, & math_range, & @@ -3283,7 +3286,7 @@ subroutine deformed_linear(res,geomdim,defgrad_av,defgrad,coord_avgCorner) end subroutine deformed_linear - +#ifdef Spectral !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ subroutine deformed_fft(res,geomdim,defgrad_av,scaling,defgrad,coords) !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ @@ -3716,6 +3719,7 @@ subroutine divergence_fdm(res,geomdim,vec_tens,order,field,divergence) enddo; enddo; enddo end subroutine divergence_fdm +#endif !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ subroutine tensor_avg(res,tensor,avg) @@ -3826,6 +3830,7 @@ subroutine calculate_cauchy(res,defgrad,p_stress,c_stress) end subroutine calculate_cauchy +#ifdef Spectral !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ subroutine math_nearestNeighborSearch(spatialDim, Favg, geomdim, queryPoints, domainPoints, querySet, domainSet, indices) !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ @@ -3882,5 +3887,6 @@ subroutine math_nearestNeighborSearch(spatialDim, Favg, geomdim, queryPoints, do indices = indices -1_pInt ! let them run from 0 to domainPoints -1 end subroutine math_nearestNeighborSearch +#endif end module math diff --git a/code/mesh.f90 b/code/mesh.f90 index 4bfbfbddb..a7a0a3d53 100644 --- a/code/mesh.f90 +++ b/code/mesh.f90 @@ -16,77 +16,87 @@ ! You should have received a copy of the GNU General Public License ! along with DAMASK. If not, see . ! -!############################################################## +!-------------------------------------------------------------------------------------------------- !* $Id$ -!############################################################## - MODULE mesh -!############################################################## +!-------------------------------------------------------------------------------------------------- +!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH +!! Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH +!! Christoph Koords, Max-Planck-Institut für Eisenforschung GmbH +!! Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH +!! Krishna Komerla, 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, only: pReal, pInt + implicit none private integer(pInt), public :: & - mesh_NcpElems, & ! total number of CP elements in mesh + mesh_NcpElems, & !< total number of CP elements in mesh mesh_NelemSets, & mesh_maxNelemInSet, & mesh_Nmaterials, & - mesh_Nnodes, & ! total number of nodes in mesh - mesh_maxNnodes, & ! max number of nodes in any CP element - mesh_maxNips, & ! max number of IPs in any CP element - mesh_maxNipNeighbors, & ! max number of IP neighbors in any CP element - mesh_maxNsharedElems, & ! max number of CP elements sharing a node + mesh_Nnodes, & !< total number of nodes in mesh + mesh_maxNnodes, & !< max number of nodes in any CP element + mesh_maxNips, & !< max number of IPs in any CP element + mesh_maxNipNeighbors, & !< max number of IP neighbors in any CP element + mesh_maxNsharedElems, & !< max number of CP elements sharing a node mesh_maxNsubNodes integer(pInt), dimension(:,:), allocatable, public :: & - mesh_element, & ! FEid, type(internal representation), material, texture, node indices - 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) + mesh_element, & !< FEid, type(internal representation), material, texture, node indices + 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(pInt), dimension(:,:,:,:), allocatable, public :: & - mesh_ipNeighborhood ! 6 or less neighboring IPs as [element_num, IP_index] + mesh_ipNeighborhood !< 6 or less neighboring IPs as [element_num, IP_index] real(pReal), dimension(:,:), allocatable, public :: & - mesh_ipVolume, & ! volume associated with IP (initially!) - mesh_node0, & ! node x,y,z coordinates (initially!) - mesh_node ! node x,y,z coordinates (after deformation! ONLY FOR MARC!!!) + mesh_ipVolume, & !< volume associated with IP (initially!) + mesh_node0, & !< node x,y,z coordinates (initially!) + mesh_node !< node x,y,z coordinates (after deformation! ONLY FOR MARC!!!) real(pReal), dimension(:,:,:), allocatable, public :: & - mesh_ipCenterOfGravity, & ! center of gravity of IP (after deformation!) - mesh_ipArea ! area of interface to neighboring IP (initially!) + mesh_ipCenterOfGravity, & !< center of gravity of IP (after deformation!) + mesh_ipArea !< area of interface to neighboring IP (initially!) real(pReal),dimension(:,:,:,:), allocatable, public :: & - mesh_ipAreaNormal ! area normal of interface to neighboring IP (initially!) + mesh_ipAreaNormal !< area normal of interface to neighboring IP (initially!) - logical, dimension(3), public :: mesh_periodicSurface ! flag indicating periodic outer surfaces (used for fluxes) + logical, dimension(3), public :: mesh_periodicSurface !< flag indicating periodic outer surfaces (used for fluxes) integer(pInt), private :: & - mesh_Nelems, & ! total number of elements in mesh - hypoelasticTableStyle, & - initialcondTableStyle - + mesh_Nelems !< total number of elements in mesh +#ifdef Marc + integer(pInt), private :: & + hypoelasticTableStyle, & !< Table style (Marc only) + initialcondTableStyle !< Table style (Marc only) +#endif + integer(pInt), dimension(2), private :: & mesh_maxValStateVar = 0_pInt 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 + mesh_nameElemSet, & !< names of elementSet + mesh_nameMaterial, & !< names of material in solid section + mesh_mapMaterial !< name of elementSet for material integer(pInt), dimension(:,:), allocatable, private :: & - mesh_mapElemSet ! list of elements in elementSet + mesh_mapElemSet !< list of elements in elementSet integer(pInt), dimension(:,:), allocatable, target, private :: & - mesh_mapFEtoCPelem, & ! [sorted FEid, corresponding CPid] - mesh_mapFEtoCPnode ! [sorted FEid, corresponding CPid] + mesh_mapFEtoCPelem, & !< [sorted FEid, corresponding CPid] + mesh_mapFEtoCPnode !< [sorted FEid, corresponding CPid] real(pReal),dimension(:,:,:), allocatable, private :: & - mesh_subNodeCoord ! coordinates of subnodes per element + mesh_subNodeCoord !< coordinates of subnodes per element - logical, private :: noPart ! for cases where the ABAQUS input file does not use part/assembly information + logical, private :: noPart !< for cases where the ABAQUS input file does not use part/assembly information -! Thee definitions should actually reside in the FE-solver specific part (different for MARC/ABAQUS) +! These definitions should actually reside in the FE-solver specific part (different for MARC/ABAQUS) ! Hence, I suggest to prefix with "FE_" integer(pInt), parameter, public :: & @@ -95,10 +105,10 @@ FE_maxNsubNodes = 56_pInt, & FE_maxNips = 27_pInt, & FE_maxNipNeighbors = 6_pInt, & - FE_maxmaxNnodesAtIP = 8_pInt, & ! max number of (equivalent) nodes attached to an IP + FE_maxmaxNnodesAtIP = 8_pInt, & !< max number of (equivalent) nodes attached to an IP FE_NipFaceNodes = 4_pInt - integer(pInt), dimension(FE_Nelemtypes), parameter, public :: FE_Nnodes = & ! nodes in a specific type of element (how we use it) + integer(pInt), dimension(FE_Nelemtypes), parameter, public :: FE_Nnodes = & !< nodes in a specific type of element (how we use it) int([8, & ! element 7 4, & ! element 134 4, & ! element 11 @@ -110,7 +120,7 @@ 8, & ! element 57 (c3d20r == c3d8 --> copy of 7) 3 & ! element 155, 125, 128 ],pInt) - integer(pInt), dimension(FE_Nelemtypes), parameter, public :: FE_Nips = & ! IPs in a specific type of element + integer(pInt), dimension(FE_Nelemtypes), parameter, public :: FE_Nips = & !< IPs in a specific type of element int([8, & ! element 7 1, & ! element 134 4, & ! element 11 @@ -122,7 +132,7 @@ 8, & ! element 57 (c3d20r == c3d8 --> copy of 7) 3 & ! element 155, 125, 128 ],pInt) - integer(pInt), dimension(FE_Nelemtypes), parameter, public :: FE_NipNeighbors = & !IP neighbors in a specific type of element + integer(pInt), dimension(FE_Nelemtypes), parameter, public :: FE_NipNeighbors = & !< IP neighbors in a specific type of element int([6, & ! element 7 4, & ! element 134 4, & ! element 11 @@ -134,7 +144,7 @@ 6, & ! element 57 (c3d20r == c3d8 --> copy of 7) 4 & ! element 155, 125, 128 ],pInt) - integer(pInt), dimension(FE_Nelemtypes), parameter, private :: FE_NsubNodes = & ! subnodes required to fully define all IP volumes + integer(pInt), dimension(FE_Nelemtypes), parameter, private :: FE_NsubNodes = & !< subnodes required to fully define all IP volumes int([19,& ! element 7 ! order is +x,-x,+y,-y,+z,-z but meaning strongly depends on Elemtype 0, & ! element 134 5, & ! element 11 @@ -146,7 +156,7 @@ 19,& ! element 57 (c3d20r == c3d8 --> copy of 7) 4 & ! element 155, 125, 128 ],pInt) - integer(pInt), dimension(FE_Nelemtypes), parameter, private :: FE_NoriginalNodes = & ! nodes in a specific type of element (how it is originally defined by marc) + integer(pInt), dimension(FE_Nelemtypes), parameter, private :: FE_NoriginalNodes = & !< nodes in a specific type of element (how it is originally defined by marc) int([8, & ! element 7 4, & ! element 134 4, & ! element 11 @@ -158,7 +168,7 @@ 20,& ! element 57 (c3d20r == c3d8 --> copy of 7) 6 & ! element 155, 125, 128 ],pInt) - integer(pInt), dimension(FE_maxNipNeighbors,FE_Nelemtypes), parameter, private :: FE_NfaceNodes = &! nodes per face in a specific type of element + integer(pInt), dimension(FE_maxNipNeighbors,FE_Nelemtypes), parameter, private :: FE_NfaceNodes = &!< nodes per face in a specific type of element reshape(int([& 4,4,4,4,4,4, & ! element 7 3,3,3,3,0,0, & ! element 134 @@ -171,7 +181,7 @@ 4,4,4,4,4,4, & ! element 57 (c3d20r == c3d8 --> copy of 7) 2,2,2,0,0,0 & ! element 155, 125, 128 ],pInt),[FE_maxNipNeighbors,FE_Nelemtypes]) - integer(pInt), dimension(FE_Nelemtypes), parameter, private :: FE_maxNnodesAtIP = & ! map IP index to two node indices in a specific type of element + integer(pInt), dimension(FE_Nelemtypes), parameter, private :: FE_maxNnodesAtIP = & !< map IP index to two node indices in a specific type of element int([1, & ! element 7 4, & ! element 134 1, & ! element 11 @@ -184,7 +194,7 @@ 1 & ! element 155, 125, 128 ],pInt) integer(pInt), dimension(FE_NipFaceNodes,FE_maxNipNeighbors,FE_Nelemtypes), parameter, private :: & - FE_nodeOnFace = & ! List of node indices on each face of a specific type of element + FE_nodeOnFace = & !< List of node indices on each face of a specific type of element reshape(int([& 1,2,3,4 , & ! element 7 2,1,5,6 , & @@ -249,8 +259,8 @@ ],pInt),[FE_NipFaceNodes,FE_maxNipNeighbors,FE_Nelemtypes]) integer(pInt), dimension(:,:,:), allocatable, private :: & - FE_nodesAtIP, & ! map IP index to two 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_nodesAtIP, & !< map IP index to two 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_subNodeParent integer(pInt), dimension(:,:,:,:), allocatable, private :: & @@ -260,63 +270,87 @@ mesh_FEasCP, & mesh_build_subNodeCoords, & mesh_build_ipVolumes, & - mesh_build_ipCoordinates, & - mesh_regrid, & - mesh_spectral_getResolution, & + mesh_build_ipCoordinates +#ifdef Spectral + public :: mesh_spectral_getResolution, & mesh_spectral_getDimension, & - mesh_spectral_getHomogenization - private :: FE_mapElemtype, & - mesh_faceMatch, & - mesh_build_FEdata, & - mesh_marc_get_tableStyles, & - mesh_get_damaskOptions, & + mesh_spectral_getHomogenization, & + mesh_regrid +#endif + + private :: & +#ifdef Spectral mesh_spectral_count_nodesAndElements, & + mesh_spectral_count_cpElements, & + mesh_spectral_map_elements, & + mesh_spectral_map_nodes, & + mesh_spectral_count_cpSizes, & + mesh_spectral_build_nodes, & + mesh_spectral_build_elements, & +#endif +#ifdef Marc + mesh_marc_get_tableStyles, & mesh_marc_count_nodesAndElements, & + mesh_marc_count_elementSets, & + mesh_marc_map_elementSets, & + mesh_marc_count_cpElements, & + mesh_marc_map_Elements, & + mesh_marc_map_nodes, & + mesh_marc_build_nodes, & + mesh_marc_count_cpSizes, & + mesh_marc_build_elements, & +#endif +#ifdef Abaqus mesh_abaqus_count_nodesAndElements, & mesh_abaqus_count_elementSets, & mesh_abaqus_count_materials, & - mesh_spectral_count_cpElements, & - mesh_abaqus_count_cpElements, & - mesh_marc_map_elementSets, & mesh_abaqus_map_elementSets, & mesh_abaqus_map_materials, & - mesh_spectral_map_nodes, & - mesh_marc_map_nodes, & - mesh_abaqus_map_nodes, & - mesh_marc_map_elements, & + mesh_abaqus_count_cpElements, & mesh_abaqus_map_elements, & - mesh_spectral_count_cpSizes, & - mesh_marc_count_cpSizes, & - mesh_abaqus_count_cpSizes, & - mesh_spectral_build_nodes, & - mesh_marc_build_nodes, & + mesh_abaqus_map_nodes, & mesh_abaqus_build_nodes, & - mesh_spectral_build_elements, & - mesh_marc_build_elements, & + mesh_abaqus_count_cpSizes, & mesh_abaqus_build_elements, & - mesh_build_ipNeighborhood, & +#endif + mesh_get_damaskOptions, & mesh_build_ipAreas, & mesh_build_nodeTwins, & - mesh_tell_statistics + mesh_build_sharedElems, & + mesh_build_ipNeighborhood, & + mesh_tell_statistics, & + FE_mapElemtype, & + mesh_faceMatch, & + mesh_build_FEdata contains -!*********************************************************** -! initialization -!*********************************************************** +!-------------------------------------------------------------------------------------------------- +!> @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,element) use DAMASK_interface use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) - use IO, only: IO_error, & - IO_open_InputFile, & - IO_abaqus_hasNoPart - use FEsolving, only: parallelExecution, & - FEsolving_execElem, & - FEsolving_execIP, & - calcMode, & - lastMode, & - FEmodelGeometry + use IO, only: & + IO_error, & +#ifdef Abaqus + IO_abaqus_hasNoPart, & +#endif +#ifdef Spectral + IO_open_file +#else + IO_open_InputFile +#endif + + use FEsolving, only: & + parallelExecution, & + FEsolving_execElem, & + FEsolving_execIP, & + calcMode, & + lastMode, & + modelName implicit none integer(pInt), parameter :: fileUnit = 222_pInt @@ -330,44 +364,46 @@ subroutine mesh_init(ip,element) !$OMP END CRITICAL (write2out) call mesh_build_FEdata ! get properties of the different types of elements +#ifdef Spectral + call IO_open_file(fileUnit,geometryFile) ! parse info from geometry file... + + call mesh_spectral_count_nodesAndElements(fileUnit) + call mesh_spectral_count_cpElements + call mesh_spectral_map_elements + call mesh_spectral_map_nodes + call mesh_spectral_count_cpSizes + call mesh_spectral_build_nodes(fileUnit) + call mesh_spectral_build_elements(fileUnit) +#endif +#ifdef Marc + call IO_open_inputFile(fileUnit,modelName) ! parse info from input file... + call mesh_marc_get_tableStyles(fileUnit) + call mesh_marc_count_nodesAndElements(fileUnit) + call mesh_marc_count_elementSets(fileUnit) + call mesh_marc_map_elementSets(fileUnit) + call mesh_marc_count_cpElements(fileUnit) + call mesh_marc_map_elements(fileUnit) + call mesh_marc_map_nodes(fileUnit) + call mesh_marc_build_nodes(fileUnit) + call mesh_marc_count_cpSizes(fileunit) + call mesh_marc_build_elements(fileUnit) +#endif +#ifdef Abaqus + call IO_open_inputFile(fileUnit,modelName) ! parse info from input file... + noPart = IO_abaqus_hasNoPart(fileUnit) + call mesh_abaqus_count_nodesAndElements(fileUnit) + call mesh_abaqus_count_elementSets(fileUnit) + call mesh_abaqus_count_materials(fileUnit) + call mesh_abaqus_map_elementSets(fileUnit) + call mesh_abaqus_map_materials(fileUnit) + call mesh_abaqus_count_cpElements(fileUnit) + call mesh_abaqus_map_elements(fileUnit) + call mesh_abaqus_map_nodes(fileUnit) + call mesh_abaqus_build_nodes(fileUnit) + call mesh_abaqus_count_cpSizes(fileunit) + call mesh_abaqus_build_elements(fileUnit) +#endif - call IO_open_inputFile(fileUnit,FEmodelGeometry) ! parse info from input file... - - select case (FEsolver) - case ('Spectral') - call mesh_spectral_count_nodesAndElements(fileUnit) - call mesh_spectral_count_cpElements - call mesh_spectral_map_elements - call mesh_spectral_map_nodes - call mesh_spectral_count_cpSizes - call mesh_spectral_build_nodes(fileUnit) - call mesh_spectral_build_elements(fileUnit) - - case ('Marc') - call mesh_marc_get_tableStyles(fileUnit) - call mesh_marc_count_nodesAndElements(fileUnit) - call mesh_marc_count_elementSets(fileUnit) - call mesh_marc_map_elementSets(fileUnit) - call mesh_marc_count_cpElements(fileUnit) - call mesh_marc_map_elements(fileUnit) - call mesh_marc_map_nodes(fileUnit) - call mesh_marc_build_nodes(fileUnit) - call mesh_marc_count_cpSizes(fileunit) - call mesh_marc_build_elements(fileUnit) - case ('Abaqus') - noPart = IO_abaqus_hasNoPart(fileUnit) - call mesh_abaqus_count_nodesAndElements(fileUnit) - call mesh_abaqus_count_elementSets(fileUnit) - call mesh_abaqus_count_materials(fileUnit) - call mesh_abaqus_map_elementSets(fileUnit) - call mesh_abaqus_map_materials(fileUnit) - call mesh_abaqus_count_cpElements(fileUnit) - call mesh_abaqus_map_elements(fileUnit) - call mesh_abaqus_map_nodes(fileUnit) - call mesh_abaqus_build_nodes(fileUnit) - call mesh_abaqus_count_cpSizes(fileunit) - call mesh_abaqus_build_elements(fileUnit) - end select call mesh_get_damaskOptions(fileUnit) close (fileUnit) @@ -380,19 +416,2635 @@ subroutine mesh_init(ip,element) call mesh_build_ipNeighborhood call mesh_tell_statistics - parallelExecution = (parallelExecution .and. (mesh_Nelems == mesh_NcpElems)) ! plus potential killer from non-local constitutive + parallelExecution = (parallelExecution .and. (mesh_Nelems == mesh_NcpElems)) ! plus potential killer from non-local constitutive FEsolving_execElem = [ 1_pInt,mesh_NcpElems] allocate(FEsolving_execIP(2_pInt,mesh_NcpElems)); FEsolving_execIP = 1_pInt forall (e = 1_pInt:mesh_NcpElems) FEsolving_execIP(2,e) = FE_Nips(mesh_element(2,e)) allocate(calcMode(mesh_maxNips,mesh_NcpElems)) - calcMode = .false. ! pretend to have collected what first call is asking (F = I) - calcMode(ip,mesh_FEasCP('elem',element)) = .true. ! first ip,el needs to be already pingponged to "calc" - lastMode = .true. ! and its mode is already known... + calcMode = .false. ! pretend to have collected what first call is asking (F = I) + calcMode(ip,mesh_FEasCP('elem',element)) = .true. ! first ip,el needs to be already pingponged to "calc" + lastMode = .true. ! and its mode is already known... end subroutine mesh_init + +!-------------------------------------------------------------------------------------------------- +!> @brief Gives the FE to CP ID mapping by binary search through lookup array +!! valid questions (what) are 'elem', 'node' +!-------------------------------------------------------------------------------------------------- +integer(pInt) function mesh_FEasCP(what,myID) + use IO, only: IO_lc + + implicit none + character(len=*), intent(in) :: what + integer(pInt), intent(in) :: myID + + integer(pInt), dimension(:,:), pointer :: lookupMap + integer(pInt) :: lower,upper,center + + mesh_FEasCP = 0_pInt + select case(IO_lc(what(1:4))) + case('elem') + lookupMap => mesh_mapFEtoCPelem + case('node') + lookupMap => mesh_mapFEtoCPnode + case default + return + endselect + + lower = 1_pInt + upper = int(size(lookupMap,2_pInt),pInt) + + if (lookupMap(1_pInt,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_pInt,lower) + return + elseif (lookupMap(1_pInt,upper) == myID) then + mesh_FEasCP = lookupMap(2_pInt,upper) + return + endif + + do while (upper-lower > 1_pInt) ! binary search in between bounds + center = (lower+upper)/2_pInt + if (lookupMap(1_pInt,center) < myID) then + lower = center + elseif (lookupMap(1_pInt,center) > myID) then + upper = center + else + mesh_FEasCP = lookupMap(2_pInt,center) + exit + endif + enddo + +end function mesh_FEasCP + + +!-------------------------------------------------------------------------------------------------- +!> @brief Assigns coordinates for subnodes in each CP element. +!! Allocates global array 'mesh_subNodeCoord' +!-------------------------------------------------------------------------------------------------- +subroutine mesh_build_subNodeCoords + + implicit none + integer(pInt) e,t,n,p,Nparents + + if (.not. allocated(mesh_subNodeCoord)) then + allocate(mesh_subNodeCoord(3,mesh_maxNnodes+mesh_maxNsubNodes,mesh_NcpElems)) + endif + mesh_subNodeCoord = 0.0_pReal + + do e = 1_pInt,mesh_NcpElems ! loop over cpElems + t = mesh_element(2,e) ! get elemType + do n = 1_pInt,FE_Nnodes(t) + mesh_subNodeCoord(1:3,n,e) = mesh_node(1:3,mesh_FEasCP('node',mesh_element(4_pInt+n,e))) ! loop over nodes of this element type + enddo + do n = 1_pInt,FE_NsubNodes(t) ! now for the true subnodes + Nparents = count(FE_subNodeParent(1_pInt:FE_Nips(t),n,t) > 0_pInt) + do p = 1_pInt,Nparents ! loop through present parent nodes + mesh_subNodeCoord(1:3,FE_Nnodes(t)+n,e) & + = mesh_subNodeCoord(1:3,FE_Nnodes(t)+n,e) & + + mesh_node(1:3,mesh_FEasCP('node',mesh_element(4_pInt+FE_subNodeParent(p,n,t),e))) ! add up parents + enddo + mesh_subNodeCoord(1:3,n+FE_Nnodes(t),e) = mesh_subNodeCoord(1:3,n+FE_Nnodes(t),e)/real(Nparents,pReal) + enddo + enddo + +end subroutine mesh_build_subNodeCoords + + +!-------------------------------------------------------------------------------------------------- +!> @brief Calculates IP volume. Allocates global array 'mesh_ipVolume' +!-------------------------------------------------------------------------------------------------- +subroutine mesh_build_ipVolumes + + use math, only: math_volTetrahedron + implicit none + + integer(pInt) :: e,f,t,i,j,n + integer(pInt), parameter :: Ntriangles = FE_NipFaceNodes-2_pInt ! each interface is made up of this many triangles + real(pReal), dimension(3,FE_NipFaceNodes) :: nPos ! coordinates of nodes on IP face + real(pReal), dimension(Ntriangles,FE_NipFaceNodes) :: volume ! volumes of possible tetrahedra + + if (.not. allocated(mesh_ipVolume)) then + allocate(mesh_ipVolume(mesh_maxNips,mesh_NcpElems)) + endif + + mesh_ipVolume = 0.0_pReal + do e = 1_pInt,mesh_NcpElems ! loop over cpElems + t = mesh_element(2_pInt,e) ! get elemType + do i = 1_pInt,FE_Nips(t) ! loop over IPs of elem + do f = 1_pInt,FE_NipNeighbors(t) ! loop over interfaces of IP and add tetrahedra which connect to CoG + forall (n = 1_pInt:FE_NipFaceNodes) & + nPos(:,n) = mesh_subNodeCoord(:,FE_subNodeOnIPFace(n,f,i,t),e) + forall (n = 1_pInt:FE_NipFaceNodes, j = 1_pInt:Ntriangles) & ! start at each interface node and build valid triangles to cover interface + volume(j,n) = math_volTetrahedron(nPos(:,n), & ! calc volume of respective tetrahedron to CoG + nPos(:,1_pInt+mod(n-1_pInt +j ,FE_NipFaceNodes)),& ! start at offset j + nPos(:,1_pInt+mod(n-1_pInt +j+1_pInt,FE_NipFaceNodes)),& ! and take j's neighbor + mesh_ipCenterOfGravity(:,i,e)) + mesh_ipVolume(i,e) = mesh_ipVolume(i,e) + sum(volume) ! add contribution from this interface + enddo + mesh_ipVolume(i,e) = mesh_ipVolume(i,e) / FE_NipFaceNodes ! renormalize with interfaceNodeNum due to loop over them + enddo + enddo + +end subroutine mesh_build_ipVolumes + + +!-------------------------------------------------------------------------------------------------- +!> @brief Calculates IP Coordinates. Allocates global array 'mesh_ipCenterOfGravity' +!-------------------------------------------------------------------------------------------------- +subroutine mesh_build_ipCoordinates + + use prec, only: tol_gravityNodePos + + implicit none + integer(pInt) :: e,f,t,i,j,k,n + logical, dimension(mesh_maxNnodes+mesh_maxNsubNodes) :: gravityNode ! flagList to find subnodes determining center of grav + real(pReal), dimension(3,mesh_maxNnodes+mesh_maxNsubNodes) :: gravityNodePos ! coordinates of subnodes determining center of grav + real(pReal), dimension(3) :: centerOfGravity + + if (.not. allocated(mesh_ipCenterOfGravity)) allocate(mesh_ipCenterOfGravity(3,mesh_maxNips,mesh_NcpElems)) + + do e = 1_pInt,mesh_NcpElems ! loop over cpElems + t = mesh_element(2,e) ! get elemType + do i = 1_pInt,FE_Nips(t) ! loop over IPs of elem + gravityNode = .false. ! reset flagList + gravityNodePos = 0.0_pReal ! reset coordinates + do f = 1_pInt,FE_NipNeighbors(t) ! loop over interfaces of IP + do n = 1_pInt,FE_NipFaceNodes ! loop over nodes on interface + gravityNode(FE_subNodeOnIPFace(n,f,i,t)) = .true. + gravityNodePos(:,FE_subNodeOnIPFace(n,f,i,t)) = mesh_subNodeCoord(:,FE_subNodeOnIPFace(n,f,i,t),e) + enddo + enddo + + do j = 1_pInt,mesh_maxNnodes+mesh_maxNsubNodes-1_pInt ! walk through entire flagList except last + if (gravityNode(j)) then ! valid node index + do k = j+1_pInt,mesh_maxNnodes+mesh_maxNsubNodes ! walk through remainder of list + if (gravityNode(k) .and. all(abs(gravityNodePos(:,j) - gravityNodePos(:,k)) < tol_gravityNodePos)) then ! found duplicate + gravityNode(j) = .false. ! delete first instance + gravityNodePos(:,j) = 0.0_pReal + exit ! continue with next suspect + endif + enddo + endif + enddo + centerOfGravity = sum(gravityNodePos,2)/real(count(gravityNode),pReal) + mesh_ipCenterOfGravity(:,i,e) = centerOfGravity + enddo + enddo + +end subroutine mesh_build_ipCoordinates + + +#ifdef Spectral +!-------------------------------------------------------------------------------------------------- +!> @brief Reads resolution information from geometry file. If fileUnit is given, +!! assumes an opened file, otherwise tries to open the one specified in geometryFile +!-------------------------------------------------------------------------------------------------- +function mesh_spectral_getResolution(fileUnit) + use IO, only: & + IO_checkAndRewind, & + IO_open_file, & + IO_stringPos, & + IO_lc, & + IO_stringValue, & + IO_intValue, & + IO_floatValue, & + IO_error + use DAMASK_interface, only: & + geometryFile + + implicit none + integer(pInt), dimension(1_pInt + 7_pInt*2_pInt) :: positions ! for a,b c + 3 values + keyword + integer(pInt), intent(in), optional :: fileUnit + integer(pInt) :: headerLength = 0_pInt + integer(pInt), dimension(3) :: mesh_spectral_getResolution + character(len=1024) :: line, & + keyword + integer(pInt) :: i, j + logical :: gotResolution = .false. + integer(pInt) :: myUnit + + if(.not. present(fileUnit)) then + myUnit = 289_pInt + call IO_open_file(myUnit,trim(geometryFile)) + else + myUnit = fileUnit + endif + + call IO_checkAndRewind(myUnit) + + read(myUnit,'(a1024)') line + positions = IO_stringPos(line,2_pInt) + keyword = IO_lc(IO_StringValue(line,positions,2_pInt)) + if (keyword(1:4) == 'head') then + headerLength = IO_intValue(line,positions,1_pInt) + 1_pInt + else + call IO_error(error_ID=841_pInt, ext_msg='mesh_spectral_getResolution') + endif + rewind(myUnit) + do i = 1_pInt, headerLength + read(myUnit,'(a1024)') line + positions = IO_stringPos(line,7_pInt) + select case ( IO_lc(IO_StringValue(line,positions,1_pInt)) ) + case ('resolution') + gotResolution = .true. + do j = 2_pInt,6_pInt,2_pInt + select case (IO_lc(IO_stringValue(line,positions,j))) + case('a') + mesh_spectral_getResolution(1) = IO_intValue(line,positions,j+1_pInt) + case('b') + mesh_spectral_getResolution(2) = IO_intValue(line,positions,j+1_pInt) + case('c') + mesh_spectral_getResolution(3) = IO_intValue(line,positions,j+1_pInt) + end select + enddo + end select + enddo + + if(.not. present(fileUnit)) close(myUnit) + + if (.not. gotResolution) & + call IO_error(error_ID = 845_pInt, ext_msg='resolution') + if((mod(mesh_spectral_getResolution(1),2_pInt)/=0_pInt .or. & ! must be a even number + mesh_spectral_getResolution(1) < 2_pInt .or. & ! and larger than 1 + mod(mesh_spectral_getResolution(2),2_pInt)/=0_pInt .or. & ! -"- + mesh_spectral_getResolution(2) < 2_pInt .or. & ! -"- + (mod(mesh_spectral_getResolution(3),2_pInt)/=0_pInt .and. & + mesh_spectral_getResolution(3)/= 1_pInt)) .or. & ! third res might be 1 + mesh_spectral_getResolution(3) < 1_pInt) & + call IO_error(error_ID = 843_pInt, ext_msg='mesh_spectral_getResolution') + +end function mesh_spectral_getResolution + + +!-------------------------------------------------------------------------------------------------- +!> @brief Reads dimension information from geometry file. If fileUnit is given, +!! assumes an opened file, otherwise tries to open the one specified in geometryFile +!-------------------------------------------------------------------------------------------------- +function mesh_spectral_getDimension(fileUnit) + use IO, only: & + IO_checkAndRewind, & + IO_open_file, & + IO_stringPos, & + IO_lc, & + IO_stringValue, & + IO_intValue, & + IO_floatValue, & + IO_error + use DAMASK_interface, only: & + geometryFile + + implicit none + integer(pInt), dimension(1_pInt + 7_pInt*2_pInt) :: positions ! for a,b c + 3 values + keyword + integer(pInt), intent(in), optional :: fileUnit + integer(pInt) :: headerLength = 0_pInt + real(pReal), dimension(3) :: mesh_spectral_getDimension + character(len=1024) :: line, & + keyword + integer(pInt) :: i, j + logical :: gotDimension = .false. + integer(pInt) :: myUnit + + if(.not. present(fileUnit)) then + myUnit = 289_pInt + call IO_open_file(myUnit,trim(geometryFile)) + else + myUnit = fileUnit + endif + + call IO_checkAndRewind(myUnit) + + read(myUnit,'(a1024)') line + positions = IO_stringPos(line,2_pInt) + keyword = IO_lc(IO_StringValue(line,positions,2_pInt)) + if (keyword(1:4) == 'head') then + headerLength = IO_intValue(line,positions,1_pInt) + 1_pInt + else + call IO_error(error_ID=841_pInt, ext_msg='mesh_spectral_getDimension') + endif + rewind(myUnit) + do i = 1_pInt, headerLength + read(myUnit,'(a1024)') line + positions = IO_stringPos(line,7_pInt) + select case ( IO_lc(IO_StringValue(line,positions,1)) ) + case ('dimension') + gotDimension = .true. + do j = 2_pInt,6_pInt,2_pInt + select case (IO_lc(IO_stringValue(line,positions,j))) + case('x') + mesh_spectral_getDimension(1) = IO_floatValue(line,positions,j+1_pInt) + case('y') + mesh_spectral_getDimension(2) = IO_floatValue(line,positions,j+1_pInt) + case('z') + mesh_spectral_getDimension(3) = IO_floatValue(line,positions,j+1_pInt) + end select + enddo + end select + enddo + + if(.not. present(fileUnit)) close(myUnit) + + if (.not. gotDimension) & + call IO_error(error_ID = 845_pInt, ext_msg='dimension') + if (any(mesh_spectral_getDimension<=0.0_pReal)) & + call IO_error(error_ID = 844_pInt, ext_msg='mesh_spectral_getDimension') + +end function mesh_spectral_getDimension + + +!-------------------------------------------------------------------------------------------------- +!> @brief Reads homogenization information from geometry file. If fileUnit is given, +!! assumes an opened file, otherwise tries to open the one specified in geometryFile +!-------------------------------------------------------------------------------------------------- +function mesh_spectral_getHomogenization(fileUnit) + use IO, only: & + IO_checkAndRewind, & + IO_open_file, & + IO_stringPos, & + IO_lc, & + IO_stringValue, & + IO_intValue, & + IO_error + use DAMASK_interface, only: & + geometryFile + + implicit none + integer(pInt), dimension(1_pInt + 7_pInt*2_pInt) :: positions ! for a, b, c + 3 values + keyword + integer(pInt), intent(in), optional :: fileUnit + integer(pInt) :: headerLength = 0_pInt + integer(pInt) :: mesh_spectral_getHomogenization + character(len=1024) :: line, & + keyword + integer(pInt) :: i, j + logical :: gotHomogenization = .false. + integer(pInt) :: myUnit + + if(.not. present(fileUnit)) then + myUnit = 289_pInt + call IO_open_file(myUnit,trim(geometryFile)) + else + myUnit = fileUnit + endif + + call IO_checkAndRewind(myUnit) + + read(myUnit,'(a1024)') line + positions = IO_stringPos(line,2_pInt) + keyword = IO_lc(IO_StringValue(line,positions,2_pInt)) + if (keyword(1:4) == 'head') then + headerLength = IO_intValue(line,positions,1_pInt) + 1_pInt + else + call IO_error(error_ID=841_pInt, ext_msg='mesh_spectral_getHomogenization') + endif + rewind(myUnit) + do i = 1_pInt, headerLength + read(myUnit,'(a1024)') line + positions = IO_stringPos(line,7_pInt) + select case ( IO_lc(IO_StringValue(line,positions,1)) ) + case ('homogenization') + gotHomogenization = .true. + mesh_spectral_getHomogenization = IO_intValue(line,positions,2_pInt) + end select + enddo + + if(.not. present(fileUnit)) close(myUnit) + + if (.not. gotHomogenization ) & + call IO_error(error_ID = 845_pInt, ext_msg='homogenization') + if (mesh_spectral_getHomogenization<1_pInt) & + call IO_error(error_ID = 842_pInt, ext_msg='mesh_spectral_getHomogenization') + +end function mesh_spectral_getHomogenization + + +!-------------------------------------------------------------------------------------------------- +!> @brief Performes a regridding from saved restart information +!-------------------------------------------------------------------------------------------------- +function mesh_regrid(resNewInput) + use prec, only: & + pInt, & + pReal + use DAMASK_interface, only: & + getSolverJobName + use IO, only: & + IO_read_jobBinaryFile ,& + IO_write_jobBinaryFile + use math, only: & + math_nearestNeighborSearch, & + deformed_FFT, & + math_mul33x3 + + integer(pInt), dimension(3) :: mesh_regrid + integer(pInt), dimension(3), optional, intent(in) :: resNewInput + integer(pInt):: maxsize, i, j, k, ielem, Npoints, NpointsNew, spatialDim + integer(pInt), dimension(3) :: res, resNew + integer(pInt), dimension(:), allocatable :: indices + real(pReal), dimension(3) :: geomdim + real(pReal), dimension(3,3) :: Favg + real(pReal), dimension(:,:), allocatable :: & + coordinatesNew, & + coordinatesLinear + + real(pReal), dimension(:,:,:), allocatable :: & + material_phase, material_phaseNew + + real(pReal), dimension(:,:,:,:,:), allocatable :: & + F, FNew, & + Fp, FpNew, & + Lp, LpNew, & + dcsdE, dcsdENew + + real(pReal), dimension (:,:,:,:,:,:,:), allocatable :: & + dPdF, dPdFNew + real(pReal), dimension (:,:,:,:), allocatable :: & + coordinates, & + Tstar, TstarNew, & + stateHomog, & + stateConst + integer(pInt), dimension(:,:), allocatable :: & + sizeStateConst, sizeStateHomog + + res = mesh_spectral_getResolution() + geomdim = mesh_spectral_getDimension() + Npoints = res(1)*res(2)*res(3) + + + + allocate(F(res(1),res(2),res(3),3,3)) + call IO_read_jobBinaryFile(777,'convergedSpectralDefgrad',trim(getSolverJobName()),size(F)) + read (777,rec=1) F + close (777) + +! ----Calculate deformed configuration and average-------- + do i= 1_pInt,3_pInt; do j = 1_pInt,3_pInt + Favg(i,j) = sum(F(1:res(1),1:res(2),1:res(3),i,j)) / Npoints + enddo; enddo + allocate(coordinates(res(1),res(2),res(3),3)) + call deformed_fft(res,geomdim,Favg,1.0_pReal,F,coordinates) + deallocate(F) + +! ----Store coordinates into a linear list-------- + allocate(coordinatesLinear(3,Npoints)) + ielem = 0_pInt + do k=1_pInt,res(3); do j=1_pInt, res(2); do i=1_pInt, res(1) + ielem = ielem + 1_pInt + coordinatesLinear(1:3,ielem) = coordinates(i,j,k,1:3) + enddo; enddo; enddo + deallocate(coordinates) + + ! ----For 2D /3D case---------------------------------- + if (res(3)== 1_pInt) then + spatialDim = 2_pInt + else + spatialDim = 3_pInt + endif + +! ----determine/calculate new resolution-------------- + if (.not. present(resNewInput)) then + resNew = res + else + if(all(resNewInput==0.0_pReal)) then + mesh_regrid =[ 0,0,0] + return !no automatic determination right now + else + resNew = resNewInput + endif + endif + NpointsNew = resNew(1)*resNew(2)*resNew(3) + +! ----Calculate regular new coordinates------------ + allocate(coordinatesNew(3,NpointsNew)) + ielem = 0_pInt + do k=1_pInt,resNew(3); do j=1_pInt, resNew(2); do i=1_pInt, resNew(1) + ielem = ielem + 1_pInt + coordinatesNew(1:3,ielem) = math_mul33x3(Favg, geomdim/real(resNew,pReal)*real([i,j,k],pReal) & + - geomdim/real(2_pInt*resNew,pReal)) + enddo; enddo; enddo + +!----- Nearest neighbour search ----------------------- + allocate(indices(Npoints)) + call math_nearestNeighborSearch(spatialDim, Favg, geomdim, NpointsNew, Npoints, & + coordinatesNew, coordinatesLinear, indices) + deallocate(coordinatesNew) + indices = indices / 3_pInt**spatialDim +1_pInt + +!--------------------------------------------------------------------------- + allocate(material_phase (1,1, Npoints)) + allocate(material_phaseNew (1,1, NpointsNew)) + call IO_read_jobBinaryFile(777,'recordedPhase',trim(getSolverJobName()),size(material_phase)) + read (777,rec=1) material_phase + close (777) + do i = 1, NpointsNew + material_phaseNew(1,1,i) = material_phase(1,1,indices(i)) + enddo + + call IO_write_jobBinaryFile(777,'recordedPhase',size(material_phaseNew)) + write (777,rec=1) material_phaseNew + close (777) + deallocate(material_phase) + deallocate(material_phaseNew) + +!--------------------------------------------------------------------------- + allocate(F (3,3,1,1, Npoints)) + allocate(FNew (3,3,1,1, NpointsNew)) + call IO_read_jobBinaryFile(777,'convergedF',trim(getSolverJobName()),size(F)) + read (777,rec=1) F + close (777) + do i = 1, NpointsNew + FNew(1:3,1:3,1,1,i) = F(1:3,1:3,1,1,indices(i)) + enddo + + call IO_write_jobBinaryFile(777,'convergedF',size(FNew)) + write (777,rec=1) FNew + close (777) + deallocate(F) + deallocate(FNew) + +!--------------------------------------------------------------------- + allocate(Fp (3,3,1,1,Npoints)) + allocate(FpNew (3,3,1,1,NpointsNew)) + call IO_read_jobBinaryFile(777,'convergedFp',trim(getSolverJobName()),size(Fp)) + read (777,rec=1) Fp + close (777) + do i = 1, NpointsNew + FpNew(1:3,1:3,1,1,i) = Fp(1:3,1:3,1,1,indices(i)) + enddo + + call IO_write_jobBinaryFile(777,'convergedFp',size(FpNew)) + write (777,rec=1) FpNew + close (777) + deallocate(Fp) + deallocate(FpNew) + +!------------------------------------------------------------------------ + allocate(Lp (3,3,1,1,Npoints)) + allocate(LpNew (3,3,1,1,NpointsNew)) + call IO_read_jobBinaryFile(777,'convergedLp',trim(getSolverJobName()),size(Lp)) + read (777,rec=1) Lp + close (777) + do i = 1, NpointsNew + LpNew(1:3,1:3,1,1,i) = Lp(1:3,1:3,1,1,indices(i)) + enddo + call IO_write_jobBinaryFile(777,'convergedLp',size(LpNew)) + write (777,rec=1) LpNew + close (777) + deallocate(Lp) + deallocate(LpNew) + +!---------------------------------------------------------------------------- + allocate(dcsdE (6,6,1,1,Npoints)) + allocate(dcsdENew (6,6,1,1,NpointsNew)) + call IO_read_jobBinaryFile(777,'convergeddcsdE',trim(getSolverJobName()),size(dcsdE)) + read (777,rec=1) dcsdE + close (777) + do i = 1, NpointsNew + dcsdENew(1:6,1:6,1,1,i) = dcsdE(1:6,1:6,1,1,indices(i)) + enddo + call IO_write_jobBinaryFile(777,'convergeddcsdE',size(dcsdENew)) + write (777,rec=1) dcsdENew + close (777) + deallocate(dcsdE) + deallocate(dcsdENew) + +!--------------------------------------------------------------------------- + allocate(dPdF (3,3,3,3,1,1,Npoints)) + allocate(dPdFNew (3,3,3,3,1,1,NpointsNew)) + call IO_read_jobBinaryFile(777,'convergeddPdF',trim(getSolverJobName()),size(dPdF)) + read (777,rec=1) dPdF + close (777) + do i = 1, NpointsNew + dPdFNew(1:3,1:3,1:3,1:3,1,1,i) = dPdF(1:3,1:3,1:3,1:3,1,1,indices(i)) + enddo + call IO_write_jobBinaryFile(777,'convergeddPdF',size(dPdFNew)) + write (777,rec=1) dPdFNew + close (777) + deallocate(dPdF) + deallocate(dPdFNew) + +!--------------------------------------------------------------------------- + allocate(Tstar (6,1,1,Npoints)) + allocate(TstarNew (6,1,1,NpointsNew)) + call IO_read_jobBinaryFile(777,'convergedTstar',trim(getSolverJobName()),size(Tstar)) + read (777,rec=1) Tstar + close (777) + do i = 1, NpointsNew + TstarNew(1:6,1,1,i) = Tstar(1:6,1,1,indices(i)) + enddo + call IO_write_jobBinaryFile(777,'convergedTstar',size(TstarNew)) + write (777,rec=1) TstarNew + close (777) + deallocate(Tstar) + deallocate(TstarNew) + +!---------------------------------------------------------------------------- + allocate(sizeStateConst(1,Npoints)) + call IO_read_jobBinaryFile(777,'sizeStateConst',trim(getSolverJobName()),size(sizeStateConst)) + read (777,rec=1) sizeStateConst + close (777) + maxsize = maxval(sizeStateConst) + allocate(StateConst (1,1,Npoints,maxsize)) + + call IO_read_jobBinaryFile(777,'convergedStateConst',trim(getSolverJobName())) + k = 0_pInt + do i =1, Npoints + do j = 1,sizeStateConst(1,i) + k = k+1_pInt + read(777,rec=k) StateConst(1,1,i,j) + enddo + enddo + close(777) + + call IO_write_jobBinaryFile(777,'convergedStateConst') + k = 0_pInt + do i = 1,NpointsNew + do j = 1,sizeStateConst(1,i) + k=k+1_pInt + write(777,rec=k) StateConst(1,1,indices(i),j) + enddo + enddo + close (777) + deallocate(sizeStateConst) + deallocate(StateConst) + +!---------------------------------------------------------------------------- + allocate(sizeStateHomog(1,Npoints)) + call IO_read_jobBinaryFile(777,'sizeStateHomog',trim(getSolverJobName()),size(sizeStateHomog)) + read (777,rec=1) sizeStateHomog + close (777) + maxsize = maxval(sizeStateHomog) + allocate(stateHomog (1,1,Npoints,maxsize)) + + call IO_read_jobBinaryFile(777,'convergedStateHomog',trim(getSolverJobName())) + k = 0_pInt + do i =1, Npoints + do j = 1,sizeStateHomog(1,i) + k = k+1_pInt + read(777,rec=k) stateHomog(1,1,i,j) + enddo + enddo + close(777) + + call IO_write_jobBinaryFile(777,'convergedStateHomog') + k = 0_pInt + do i = 1,NpointsNew + do j = 1,sizeStateHomog(1,i) + k=k+1_pInt + write(777,rec=k) stateHomog(1,1,indices(i),j) + enddo + enddo + close (777) + deallocate(sizeStateHomog) + deallocate(stateHomog) + + deallocate(indices) + + mesh_regrid = resNew +end function mesh_regrid + + +!-------------------------------------------------------------------------------------------------- +!> @brief Count overall number of nodes and elements in mesh and stores them in +!! 'mesh_Nelems' and 'mesh_Nnodes' +!-------------------------------------------------------------------------------------------------- +subroutine mesh_spectral_count_nodesAndElements(myUnit) + + implicit none + integer(pInt), intent(in) :: myUnit + integer(pInt), dimension(3) :: res + + res = mesh_spectral_getResolution(myUnit) + mesh_Nelems = res(1)*res(2)*res(3) + mesh_Nnodes = (1_pInt + res(1))*(1_pInt + res(2))*(1_pInt + res(3)) + +end subroutine mesh_spectral_count_nodesAndElements + + +!-------------------------------------------------------------------------------------------------- +!> @brief Count overall number of CP elements in mesh and stores them in 'mesh_NcpElems' +!-------------------------------------------------------------------------------------------------- +subroutine mesh_spectral_count_cpElements + + implicit none + + mesh_NcpElems = mesh_Nelems + +end subroutine mesh_spectral_count_cpElements + + +!-------------------------------------------------------------------------------------------------- +!> @brief Maps elements from FE ID to internal (consecutive) representation. +!! Allocates global array 'mesh_mapFEtoCPelem' +!-------------------------------------------------------------------------------------------------- +subroutine mesh_spectral_map_elements + + implicit none + integer(pInt) :: i + + allocate (mesh_mapFEtoCPelem(2_pInt,mesh_NcpElems)) ; mesh_mapFEtoCPelem = 0_pInt + + forall (i = 1_pInt:mesh_NcpElems) & + mesh_mapFEtoCPelem(1:2,i) = i + +end subroutine mesh_spectral_map_elements + + +!-------------------------------------------------------------------------------------------------- +!> @brief Maps node from FE ID to internal (consecutive) representation. +!! Allocates global array 'mesh_mapFEtoCPnode' +!-------------------------------------------------------------------------------------------------- +subroutine mesh_spectral_map_nodes + + implicit none + integer(pInt) :: i + + allocate (mesh_mapFEtoCPnode(2_pInt,mesh_Nnodes)) ; mesh_mapFEtoCPnode = 0_pInt + + forall (i = 1_pInt:mesh_Nnodes) & + mesh_mapFEtoCPnode(1:2,i) = i + +end subroutine mesh_spectral_map_nodes + + +!-------------------------------------------------------------------------------------------------- +!> @brief Gets maximum count of nodes, IPs, IP neighbors, and subNodes among cpElements. +!! Allocates global arrays 'mesh_maxNnodes', 'mesh_maxNips', mesh_maxNipNeighbors', +!! and mesh_maxNsubNodes +!-------------------------------------------------------------------------------------------------- +subroutine mesh_spectral_count_cpSizes + + implicit none + integer(pInt) :: t + + t = FE_mapElemtype('C3D8R') ! fake 3D hexahedral 8 node 1 IP element + + mesh_maxNnodes = FE_Nnodes(t) + mesh_maxNips = FE_Nips(t) + mesh_maxNipNeighbors = FE_NipNeighbors(t) + mesh_maxNsubNodes = FE_NsubNodes(t) + +end subroutine mesh_spectral_count_cpSizes + + +!-------------------------------------------------------------------------------------------------- +!> @brief Store x,y,z coordinates of all nodes in mesh. +!! Allocates global arrays 'mesh_node0' and 'mesh_node' +!-------------------------------------------------------------------------------------------------- +subroutine mesh_spectral_build_nodes(myUnit) + + use IO, only: & + IO_error + + implicit none + integer(pInt), intent(in) :: myUnit + integer(pInt) :: n + integer(pInt), dimension(3) :: res = 1_pInt + real(pReal), dimension(3) :: geomdim = 1.0_pReal + + allocate ( mesh_node0 (3,mesh_Nnodes) ); mesh_node0 = 0.0_pReal + allocate ( mesh_node (3,mesh_Nnodes) ); mesh_node = 0.0_pReal + + res = mesh_spectral_getResolution(myUnit) + geomdim = mesh_spectral_getDimension(myUnit) + + forall (n = 0_pInt:mesh_Nnodes-1_pInt) + mesh_node0(1,n+1_pInt) = & + geomdim(1) * real(mod(n,(res(1)+1_pInt) ),pReal) & + / real(res(1),pReal) + mesh_node0(2,n+1_pInt) = & + geomdim(2) * real(mod(n/(res(1)+1_pInt),(res(2)+1_pInt)),pReal) & + / real(res(2),pReal) + mesh_node0(3,n+1_pInt) = & + geomdim(3) * real(mod(n/(res(1)+1_pInt)/(res(2)+1_pInt),(res(3)+1_pInt)),pReal) & + / real(res(3),pReal) + end forall + + mesh_node = mesh_node0 !why? + +end subroutine mesh_spectral_build_nodes + + +!-------------------------------------------------------------------------------------------------- +!> @brief Store FEid, type, material, texture, and node list per element. +!! Allocates global array 'mesh_element' +!-------------------------------------------------------------------------------------------------- +subroutine mesh_spectral_build_elements(myUnit) + + use IO, only: & + IO_checkAndRewind, & + IO_lc, & + IO_stringValue, & + IO_stringPos, & + IO_error, & + IO_continuousIntValues, & + IO_intValue, & + IO_countContinuousIntValues + + implicit none + integer(pInt), intent(in) :: myUnit + + integer(pInt), parameter :: maxNchunks = 7_pInt + integer(pInt), dimension (1_pInt+2_pInt*maxNchunks) :: myPos + integer(pInt), dimension(3) :: res + integer(pInt) :: e, i, j, homog = 0_pInt, headerLength = 0_pInt, maxIntCount + integer(pInt), dimension(:), allocatable :: microstructures + integer(pInt), dimension(1,1) :: dummySet = 0_pInt + character(len=65536) :: line,keyword + character(len=64), dimension(1) :: dummyName = '' + + res = mesh_spectral_getResolution(myUnit) + homog = mesh_spectral_getHomogenization(myUnit) + + call IO_checkAndRewind(myUnit) + + read(myUnit,'(a65536)') line + myPos = IO_stringPos(line,2_pInt) + keyword = IO_lc(IO_StringValue(line,myPos,2_pInt)) + if (keyword(1:4) == 'head') then + headerLength = IO_intValue(line,myPos,1_pInt) + 1_pInt + else + call IO_error(error_ID=841_pInt, ext_msg='mesh_spectral_build_elements') + endif + + rewind(myUnit) + do i = 1_pInt, headerLength + read(myUnit,'(a65536)') line + enddo + + maxIntCount = 0_pInt + i = 1_pInt + + do while (i > 0_pInt) + i = IO_countContinuousIntValues(myUnit) + maxIntCount = max(maxIntCount, i) + enddo + + rewind (myUnit) + do i=1_pInt,headerLength ! skip header + read(myUnit,'(a65536)') line + enddo + + allocate (mesh_element (4_pInt+mesh_maxNnodes,mesh_NcpElems)) ; mesh_element = 0_pInt + allocate (microstructures (1_pInt+maxIntCount)) ; microstructures = 2_pInt + + e = 0_pInt + do while (e < mesh_NcpElems .and. microstructures(1) > 0_pInt) ! fill expected number of elements, stop at end of data (or blank line!) + microstructures = IO_continuousIntValues(myUnit,maxIntCount,dummyName,dummySet,0_pInt) ! get affected elements + do i = 1_pInt,microstructures(1_pInt) + e = e+1_pInt ! valid element entry + mesh_element( 1,e) = e ! FE id + mesh_element( 2,e) = FE_mapElemtype('C3D8R') ! elem type + mesh_element( 3,e) = homog ! homogenization + mesh_element( 4,e) = microstructures(1_pInt+i) ! microstructure + mesh_element( 5,e) = e + (e-1_pInt)/res(1) + & + ((e-1_pInt)/(res(1)*res(2)))*(res(1)+1_pInt) ! base node + mesh_element( 6,e) = mesh_element(5,e) + 1_pInt + mesh_element( 7,e) = mesh_element(5,e) + res(1) + 2_pInt + mesh_element( 8,e) = mesh_element(5,e) + res(1) + 1_pInt + mesh_element( 9,e) = mesh_element(5,e) +(res(1) + 1_pInt) * (res(2) + 1_pInt) ! second floor base node + mesh_element(10,e) = mesh_element(9,e) + 1_pInt + mesh_element(11,e) = mesh_element(9,e) + res(1) + 2_pInt + mesh_element(12,e) = mesh_element(9,e) + res(1) + 1_pInt + mesh_maxValStateVar(1) = max(mesh_maxValStateVar(1),mesh_element(3,e)) !needed for statistics + mesh_maxValStateVar(2) = max(mesh_maxValStateVar(2),mesh_element(4,e)) + enddo + enddo + + deallocate(microstructures) + if (e /= mesh_NcpElems) call IO_error(880_pInt,e) + +end subroutine mesh_spectral_build_elements +#endif + + +#ifdef Marc +!-------------------------------------------------------------------------------------------------- +!> @brief Figures out table styles (Marc only) and stores to 'initialcondTableStyle' and +!! 'hypoelasticTableStyle' +!-------------------------------------------------------------------------------------------------- +subroutine mesh_marc_get_tableStyles(myUnit) + + use IO, only: & + IO_lc, & + IO_intValue, & + IO_stringValue, & + IO_stringPos + + implicit none + integer(pInt), intent(in) :: myUnit + + integer(pInt), parameter :: maxNchunks = 6_pInt + integer(pInt), dimension (1+2*maxNchunks) :: myPos + character(len=300) line + + initialcondTableStyle = 0_pInt + hypoelasticTableStyle = 0_pInt + +610 FORMAT(A300) + + rewind(myUnit) + do + read (myUnit,610,END=620) line + myPos = IO_stringPos(line,maxNchunks) + + if ( IO_lc(IO_stringValue(line,myPos,1_pInt)) == 'table' .and. myPos(1_pInt) .GT. 5) then + initialcondTableStyle = IO_intValue(line,myPos,4_pInt) + hypoelasticTableStyle = IO_intValue(line,myPos,5_pInt) + exit + endif + enddo + +620 end subroutine mesh_marc_get_tableStyles + + +!-------------------------------------------------------------------------------------------------- +!> @brief Count overall number of nodes and elements in mesh and stores them in +!! 'mesh_Nelems' and 'mesh_Nnodes' +!-------------------------------------------------------------------------------------------------- +subroutine mesh_marc_count_nodesAndElements(myUnit) + + use IO, only: IO_lc, & + IO_stringValue, & + IO_stringPos, & + IO_IntValue + + implicit none + integer(pInt), intent(in) :: myUnit + + integer(pInt), parameter :: maxNchunks = 4_pInt + integer(pInt), dimension (1+2*maxNchunks) :: myPos + character(len=300) line + + mesh_Nnodes = 0_pInt + mesh_Nelems = 0_pInt + +610 FORMAT(A300) + + rewind(myUnit) + do + read (myUnit,610,END=620) line + myPos = IO_stringPos(line,maxNchunks) + + if ( IO_lc(IO_StringValue(line,myPos,1_pInt)) == 'sizing') then + mesh_Nelems = IO_IntValue (line,myPos,3_pInt) + mesh_Nnodes = IO_IntValue (line,myPos,4_pInt) + exit + endif + enddo + +620 end subroutine mesh_marc_count_nodesAndElements + + +!-------------------------------------------------------------------------------------------------- +!> @brief Count overall number of element sets in mesh. Stores to 'mesh_NelemSets', and +!! 'mesh_maxNelemInSet' +!-------------------------------------------------------------------------------------------------- + subroutine mesh_marc_count_elementSets(myUnit) + + use IO, only: IO_lc, & + IO_stringValue, & + IO_stringPos, & + IO_countContinuousIntValues + + implicit none + integer(pInt), intent(in) :: myUnit + + integer(pInt), parameter :: maxNchunks = 2_pInt + integer(pInt), dimension (1+2*maxNchunks) :: myPos + character(len=300) line + + mesh_NelemSets = 0_pInt + mesh_maxNelemInSet = 0_pInt + +610 FORMAT(A300) + + rewind(myUnit) + do + read (myUnit,610,END=620) line + myPos = IO_stringPos(line,maxNchunks) + + if ( IO_lc(IO_StringValue(line,myPos,1_pInt)) == 'define' .and. & + IO_lc(IO_StringValue(line,myPos,2_pInt)) == 'element' ) then + mesh_NelemSets = mesh_NelemSets + 1_pInt + mesh_maxNelemInSet = max(mesh_maxNelemInSet, & + IO_countContinuousIntValues(myUnit)) + endif + enddo + +620 end subroutine mesh_marc_count_elementSets + + +!******************************************************************** +! map element sets +! +! allocate globals: mesh_nameElemSet, mesh_mapElemSet +!******************************************************************** +subroutine mesh_marc_map_elementSets(myUnit) + + use IO, only: IO_lc, & + IO_stringValue, & + IO_stringPos, & + IO_continuousIntValues + + implicit none + integer(pInt), intent(in) :: myUnit + + integer(pInt), parameter :: maxNchunks = 4_pInt + integer(pInt), dimension (1+2*maxNchunks) :: myPos + character(len=300) :: line + integer(pInt) :: elemSet = 0_pInt + + allocate (mesh_nameElemSet(mesh_NelemSets)) ; mesh_nameElemSet = '' + allocate (mesh_mapElemSet(1_pInt+mesh_maxNelemInSet,mesh_NelemSets)) ; mesh_mapElemSet = 0_pInt + +610 FORMAT(A300) + + rewind(myUnit) + do + read (myUnit,610,END=640) line + myPos = IO_stringPos(line,maxNchunks) + if( (IO_lc(IO_stringValue(line,myPos,1_pInt)) == 'define' ) .and. & + (IO_lc(IO_stringValue(line,myPos,2_pInt)) == 'element' ) ) then + elemSet = elemSet+1_pInt + mesh_nameElemSet(elemSet) = trim(IO_stringValue(line,myPos,4_pInt)) + mesh_mapElemSet(:,elemSet) = IO_continuousIntValues(myUnit,mesh_maxNelemInSet,mesh_nameElemSet,mesh_mapElemSet,mesh_NelemSets) + endif + enddo + +640 end subroutine mesh_marc_map_elementSets + + +!-------------------------------------------------------------------------------------------------- +!> @brief Count overall number of CP elements in mesh and stores them in 'mesh_NcpElems' +!-------------------------------------------------------------------------------------------------- +subroutine mesh_marc_count_cpElements(myUnit) + + use IO, only: IO_lc, & + IO_stringValue, & + IO_stringPos, & + IO_countContinuousIntValues + + implicit none + integer(pInt), intent(in) :: myUnit + + integer(pInt), parameter :: maxNchunks = 1_pInt + integer(pInt), dimension (1+2*maxNchunks) :: myPos + integer(pInt) :: i + character(len=300):: line + + mesh_NcpElems = 0_pInt + +610 FORMAT(A300) + + rewind(myUnit) + do + read (myUnit,610,END=620) line + myPos = IO_stringPos(line,maxNchunks) + + if ( IO_lc(IO_stringValue(line,myPos,1_pInt)) == 'hypoelastic') then + do i=1_pInt,3_pInt+hypoelasticTableStyle ! Skip 3 or 4 lines + read (myUnit,610,END=620) line + enddo + mesh_NcpElems = mesh_NcpElems + IO_countContinuousIntValues(myUnit) + exit + endif + enddo + +620 end subroutine mesh_marc_count_cpElements + + +!-------------------------------------------------------------------------------------------------- +!> @brief Maps elements from FE ID to internal (consecutive) representation. +!! Allocates global array 'mesh_mapFEtoCPelem' +!-------------------------------------------------------------------------------------------------- +subroutine mesh_marc_map_elements(myUnit) + + use math, only: qsort + use IO, only: IO_lc, & + IO_stringValue, & + IO_stringPos, & + IO_continuousIntValues + + implicit none + integer(pInt), intent(in) :: myUnit + + integer(pInt), parameter :: maxNchunks = 1_pInt + integer(pInt), dimension (1_pInt+2_pInt*maxNchunks) :: myPos + character(len=300) line + + integer(pInt), dimension (1_pInt+mesh_NcpElems) :: contInts + integer(pInt) :: i,cpElem = 0_pInt + + allocate (mesh_mapFEtoCPelem(2,mesh_NcpElems)) ; mesh_mapFEtoCPelem = 0_pInt + +610 FORMAT(A300) + + rewind(myUnit) + do + read (myUnit,610,END=660) line + myPos = IO_stringPos(line,maxNchunks) + if( IO_lc(IO_stringValue(line,myPos,1_pInt)) == 'hypoelastic' ) then + do i=1_pInt,3_pInt+hypoelasticTableStyle ! skip three (or four if new table style!) lines + read (myUnit,610,END=660) line + enddo + contInts = IO_continuousIntValues(myUnit,mesh_NcpElems,mesh_nameElemSet,& + mesh_mapElemSet,mesh_NelemSets) + do i = 1_pInt,contInts(1) + cpElem = cpElem+1_pInt + mesh_mapFEtoCPelem(1,cpElem) = contInts(1_pInt+i) + mesh_mapFEtoCPelem(2,cpElem) = cpElem + enddo + endif + enddo + +660 call qsort(mesh_mapFEtoCPelem,1_pInt,int(size(mesh_mapFEtoCPelem,2_pInt),pInt)) ! should be mesh_NcpElems + +end subroutine mesh_marc_map_elements + + +!-------------------------------------------------------------------------------------------------- +!> @brief Maps node from FE ID to internal (consecutive) representation. +!! Allocates global array 'mesh_mapFEtoCPnode' +!-------------------------------------------------------------------------------------------------- +subroutine mesh_marc_map_nodes(myUnit) + + use math, only: qsort + use IO, only: IO_lc, & + IO_stringValue, & + IO_stringPos, & + IO_fixedIntValue + + implicit none + integer(pInt), intent(in) :: myUnit + + integer(pInt), parameter :: maxNchunks = 1_pInt + integer(pInt), dimension (1_pInt+2_pInt*maxNchunks) :: myPos + character(len=300) line + + integer(pInt), dimension (mesh_Nnodes) :: node_count + integer(pInt) :: i + + allocate (mesh_mapFEtoCPnode(2_pInt,mesh_Nnodes)) ; mesh_mapFEtoCPnode = 0_pInt + +610 FORMAT(A300) + + node_count = 0_pInt + + rewind(myUnit) + do + read (myUnit,610,END=650) line + myPos = IO_stringPos(line,maxNchunks) + if( IO_lc(IO_stringValue(line,myPos,1_pInt)) == 'coordinates' ) then + read (myUnit,610,END=650) line ! skip crap line + do i = 1_pInt,mesh_Nnodes + read (myUnit,610,END=650) line + mesh_mapFEtoCPnode(1_pInt,i) = IO_fixedIntValue (line,[ 0_pInt,10_pInt],1_pInt) + mesh_mapFEtoCPnode(2_pInt,i) = i + enddo + exit + endif + enddo + +650 call qsort(mesh_mapFEtoCPnode,1_pInt,int(size(mesh_mapFEtoCPnode,2_pInt),pInt)) + +end subroutine mesh_marc_map_nodes + + +!-------------------------------------------------------------------------------------------------- +!> @brief store x,y,z coordinates of all nodes in mesh. +!! Allocates global arrays 'mesh_node0' and 'mesh_node' +!-------------------------------------------------------------------------------------------------- +subroutine mesh_marc_build_nodes(myUnit) + + use IO, only: IO_lc, & + IO_stringValue, & + IO_stringPos, & + IO_fixedIntValue, & + IO_fixedNoEFloatValue + + implicit none + integer(pInt), intent(in) :: myUnit + + integer(pInt), dimension(5), parameter :: node_ends = int([0,10,30,50,70],pInt) + integer(pInt), parameter :: maxNchunks = 1_pInt + integer(pInt), dimension (1_pInt+2_pInt*maxNchunks) :: myPos + character(len=300) :: line + integer(pInt) :: i,j,m + + allocate ( mesh_node0 (3,mesh_Nnodes) ); mesh_node0 = 0.0_pReal + allocate ( mesh_node (3,mesh_Nnodes) ); mesh_node = 0.0_pReal + +610 FORMAT(A300) + + rewind(myUnit) + do + read (myUnit,610,END=670) line + myPos = IO_stringPos(line,maxNchunks) + if( IO_lc(IO_stringValue(line,myPos,1_pInt)) == 'coordinates' ) then + read (myUnit,610,END=670) line ! skip crap line + do i=1_pInt,mesh_Nnodes + read (myUnit,610,END=670) line + m = mesh_FEasCP('node',IO_fixedIntValue(line,node_ends,1_pInt)) + forall (j = 1_pInt:3_pInt) mesh_node0(j,m) = IO_fixedNoEFloatValue(line,node_ends,j+1_pInt) + enddo + exit + endif + enddo + +670 mesh_node = mesh_node0 + +end subroutine mesh_marc_build_nodes + + +!-------------------------------------------------------------------------------------------------- +!> @brief Gets maximum count of nodes, IPs, IP neighbors, and subNodes among cpElements. +!! Allocates global arrays 'mesh_maxNnodes', 'mesh_maxNips', mesh_maxNipNeighbors', +!! and mesh_maxNsubNodes +!-------------------------------------------------------------------------------------------------- +subroutine mesh_marc_count_cpSizes(myUnit) + + use IO, only: IO_lc, & + IO_stringValue, & + IO_stringPos, & + IO_intValue, & + IO_skipChunks + + implicit none + integer(pInt), intent(in) :: myUnit + + integer(pInt), parameter :: maxNchunks = 2_pInt + integer(pInt), dimension (1_pInt+2_pInt*maxNchunks) :: myPos + character(len=300) :: line + integer(pInt) :: i,t,e + + mesh_maxNnodes = 0_pInt + mesh_maxNips = 0_pInt + mesh_maxNipNeighbors = 0_pInt + mesh_maxNsubNodes = 0_pInt + +610 FORMAT(A300) + rewind(myUnit) + do + read (myUnit,610,END=630) line + myPos = IO_stringPos(line,maxNchunks) + if( IO_lc(IO_stringValue(line,myPos,1_pInt)) == 'connectivity' ) then + read (myUnit,610,END=630) line ! Garbage line + do i=1_pInt,mesh_Nelems ! read all elements + read (myUnit,610,END=630) line + myPos = IO_stringPos(line,maxNchunks) ! limit to id and type + e = mesh_FEasCP('elem',IO_intValue(line,myPos,1_pInt)) + if (e /= 0_pInt) then + t = FE_mapElemtype(IO_stringValue(line,myPos,2_pInt)) + mesh_maxNnodes = max(mesh_maxNnodes,FE_Nnodes(t)) + mesh_maxNips = max(mesh_maxNips,FE_Nips(t)) + mesh_maxNipNeighbors = max(mesh_maxNipNeighbors,FE_NipNeighbors(t)) + mesh_maxNsubNodes = max(mesh_maxNsubNodes,FE_NsubNodes(t)) + call IO_skipChunks(myUnit,FE_NoriginalNodes(t)-(myPos(1_pInt)-2_pInt)) ! read on if FE_Nnodes exceeds node count present on current line + endif + enddo + exit + endif + enddo + +630 end subroutine mesh_marc_count_cpSizes + + +!-------------------------------------------------------------------------------------------------- +!> @brief Store FEid, type, mat, tex, and node list per elemen. +!! Allocates global array 'mesh_element' +!-------------------------------------------------------------------------------------------------- +subroutine mesh_marc_build_elements(myUnit) + + use IO, only: IO_lc, & + IO_stringValue, & + IO_fixedNoEFloatValue, & + IO_skipChunks, & + IO_stringPos, & + IO_intValue, & + IO_continuousIntValues + + implicit none + integer(pInt), intent(in) :: myUnit + + integer(pInt), parameter :: maxNchunks = 66_pInt + integer(pInt), dimension (1_pInt+2_pInt*maxNchunks) :: myPos + character(len=300) line + + integer(pInt), dimension(1_pInt+mesh_NcpElems) :: contInts + integer(pInt) :: i,j,sv,myVal,e + + allocate (mesh_element (4_pInt+mesh_maxNnodes,mesh_NcpElems)) ; mesh_element = 0_pInt + +610 FORMAT(A300) + + rewind(myUnit) + do + read (myUnit,610,END=620) line + myPos(1:1+2*1) = IO_stringPos(line,1_pInt) + if( IO_lc(IO_stringValue(line,myPos,1_pInt)) == 'connectivity' ) then + read (myUnit,610,END=620) line ! Garbage line + do i = 1_pInt,mesh_Nelems + read (myUnit,610,END=620) line + myPos = IO_stringPos(line,maxNchunks) ! limit to 64 nodes max (plus ID, type) + e = mesh_FEasCP('elem',IO_intValue(line,myPos,1_pInt)) + if (e /= 0_pInt) then ! disregard non CP elems + mesh_element(1,e) = IO_IntValue (line,myPos,1_pInt) ! FE id + mesh_element(2,e) = FE_mapElemtype(IO_StringValue(line,myPos,2_pInt)) ! elem type + forall (j = 1_pInt:FE_Nnodes(mesh_element(2,e))) & + mesh_element(j+4_pInt,e) = IO_IntValue(line,myPos,j+2_pInt) ! copy FE ids of nodes + call IO_skipChunks(myUnit,FE_NoriginalNodes(mesh_element(2_pInt,e))-(myPos(1_pInt)-2_pInt)) ! read on if FE_Nnodes exceeds node count present on current line + endif + enddo + exit + endif + enddo + +620 rewind(myUnit) ! just in case "initial state" apears before "connectivity" + read (myUnit,610,END=620) line + do + myPos(1:1+2*2) = IO_stringPos(line,2_pInt) + if( (IO_lc(IO_stringValue(line,myPos,1_pInt)) == 'initial') .and. & + (IO_lc(IO_stringValue(line,myPos,2_pInt)) == 'state') ) then + if (initialcondTableStyle == 2_pInt) read (myUnit,610,END=620) line ! read extra line for new style + read (myUnit,610,END=630) line ! read line with index of state var + myPos(1:1+2*1) = IO_stringPos(line,1_pInt) + sv = IO_IntValue(line,myPos,1_pInt) ! figure state variable index + if( (sv == 2_pInt).or.(sv == 3_pInt) ) then ! only state vars 2 and 3 of interest + read (myUnit,610,END=620) line ! read line with value of state var + myPos(1:1+2*1) = IO_stringPos(line,1_pInt) + do while (scan(IO_stringValue(line,myPos,1_pInt),'+-',back=.true.)>1) ! is noEfloat value? + myVal = nint(IO_fixedNoEFloatValue(line,[0_pInt,20_pInt],1_pInt),pInt) ! state var's value + mesh_maxValStateVar(sv-1_pInt) = max(myVal,mesh_maxValStateVar(sv-1_pInt)) ! remember max val of homogenization and microstructure index + if (initialcondTableStyle == 2_pInt) then + read (myUnit,610,END=630) line ! read extra line + read (myUnit,610,END=630) line ! read extra line + endif + contInts = IO_continuousIntValues& ! get affected elements + (myUnit,mesh_Nelems,mesh_nameElemSet,mesh_mapElemSet,mesh_NelemSets) + do i = 1_pInt,contInts(1) + e = mesh_FEasCP('elem',contInts(1_pInt+i)) + mesh_element(1_pInt+sv,e) = myVal + enddo + if (initialcondTableStyle == 0_pInt) read (myUnit,610,END=620) line ! ignore IP range for old table style + read (myUnit,610,END=630) line + myPos(1:1+2*1) = IO_stringPos(line,1_pInt) + enddo + endif + else + read (myUnit,610,END=630) line + endif + enddo + +630 end subroutine mesh_marc_build_elements +#endif + +#ifdef Abaqus +!-------------------------------------------------------------------------------------------------- +!> @brief Count overall number of nodes and elements in mesh and stores them in +!! 'mesh_Nelems' and 'mesh_Nnodes' +!-------------------------------------------------------------------------------------------------- +subroutine mesh_abaqus_count_nodesAndElements(myUnit) + + use IO, only: IO_lc, & + IO_stringValue, & + IO_stringPos, & + IO_countDataLines, & + IO_error + + implicit none + integer(pInt), intent(in) :: myUnit + + integer(pInt), parameter :: maxNchunks = 2_pInt + integer(pInt), dimension (1+2*maxNchunks) :: myPos + character(len=300) :: line + logical :: inPart + + mesh_Nnodes = 0_pInt + mesh_Nelems = 0_pInt + +610 FORMAT(A300) + + inPart = .false. + rewind(myUnit) + do + read (myUnit,610,END=620) line + myPos = IO_stringPos(line,maxNchunks) + if ( IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*part' ) inPart = .true. + if ( IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*end' .and. & + IO_lc(IO_stringValue(line,myPos,2_pInt)) == 'part' ) inPart = .false. + + if (inPart .or. noPart) then + select case ( IO_lc(IO_stringValue(line,myPos,1_pInt))) + case('*node') + if( & + IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'output' .and. & + IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'print' .and. & + IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'file' .and. & + IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'response' & + ) & + mesh_Nnodes = mesh_Nnodes + IO_countDataLines(myUnit) + case('*element') + if( & + IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'output' .and. & + IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'matrix' .and. & + IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'response' & + ) then + mesh_Nelems = mesh_Nelems + IO_countDataLines(myUnit) + endif + endselect + endif + enddo + +620 if (mesh_Nnodes < 2_pInt) call IO_error(error_ID=900_pInt) + if (mesh_Nelems == 0_pInt) call IO_error(error_ID=901_pInt) + +end subroutine mesh_abaqus_count_nodesAndElements + + +!******************************************************************** +! count overall number of element sets in mesh +! +! mesh_NelemSets, mesh_maxNelemInSet +!******************************************************************** +subroutine mesh_abaqus_count_elementSets(myUnit) + + use IO, only: IO_lc, & + IO_stringValue, & + IO_stringPos, & + IO_error + + implicit none + integer(pInt), intent(in) :: myUnit + + integer(pInt), parameter :: maxNchunks = 2_pInt + integer(pInt), dimension (1+2*maxNchunks) :: myPos + character(len=300) :: line + logical :: inPart + + mesh_NelemSets = 0_pInt + mesh_maxNelemInSet = mesh_Nelems ! have to be conservative, since Abaqus allows for recursive definitons + +610 FORMAT(A300) + + inPart = .false. + rewind(myUnit) + do + read (myUnit,610,END=620) line + myPos = IO_stringPos(line,maxNchunks) + if ( IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*part' ) inPart = .true. + if ( IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*end' .and. & + IO_lc(IO_stringValue(line,myPos,2_pInt)) == 'part' ) inPart = .false. + + if ( (inPart .or. noPart) .and. IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*elset' ) & + mesh_NelemSets = mesh_NelemSets + 1_pInt + enddo + +620 continue + if (mesh_NelemSets == 0) call IO_error(error_ID=902_pInt) + +end subroutine mesh_abaqus_count_elementSets + + +!******************************************************************** +! count overall number of solid sections sets in mesh (Abaqus only) +! +! mesh_Nmaterials +!******************************************************************** +subroutine mesh_abaqus_count_materials(myUnit) + + use IO, only: IO_lc, & + IO_stringValue, & + IO_stringPos, & + IO_error + + implicit none + integer(pInt), intent(in) :: myUnit + + integer(pInt), parameter :: maxNchunks = 2_pInt + integer(pInt), dimension (1_pInt+2_pInt*maxNchunks) :: myPos + character(len=300) :: line + logical inPart + + mesh_Nmaterials = 0_pInt + +610 FORMAT(A300) + + inPart = .false. + rewind(myUnit) + do + read (myUnit,610,END=620) line + myPos = IO_stringPos(line,maxNchunks) + if ( IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*part' ) inPart = .true. + if ( IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*end' .and. & + IO_lc(IO_stringValue(line,myPos,2_pInt)) == 'part' ) inPart = .false. + + if ( (inPart .or. noPart) .and. & + IO_lc(IO_StringValue(line,myPos,1_pInt)) == '*solid' .and. & + IO_lc(IO_StringValue(line,myPos,2_pInt)) == 'section' ) & + mesh_Nmaterials = mesh_Nmaterials + 1_pInt + enddo + +620 if (mesh_Nmaterials == 0_pInt) call IO_error(error_ID=903_pInt) + +end subroutine mesh_abaqus_count_materials + + +!******************************************************************** +! Build element set mapping +! +! allocate globals: mesh_nameElemSet, mesh_mapElemSet +!******************************************************************** +subroutine mesh_abaqus_map_elementSets(myUnit) + + use IO, only: IO_lc, & + IO_stringValue, & + IO_stringPos, & + IO_extractValue, & + IO_continuousIntValues, & + IO_error + + implicit none + integer(pInt), intent(in) :: myUnit + + integer(pInt), parameter :: maxNchunks = 4_pInt + integer(pInt), dimension (1_pInt+2_pInt*maxNchunks) :: myPos + character(len=300) :: line + integer(pInt) :: elemSet = 0_pInt,i + logical :: inPart = .false. + + allocate (mesh_nameElemSet(mesh_NelemSets)) ; mesh_nameElemSet = '' + allocate (mesh_mapElemSet(1_pInt+mesh_maxNelemInSet,mesh_NelemSets)) ; mesh_mapElemSet = 0_pInt + +610 FORMAT(A300) + + + rewind(myUnit) + do + read (myUnit,610,END=640) line + myPos = IO_stringPos(line,maxNchunks) + if ( IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*part' ) inPart = .true. + if ( IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*end' .and. & + IO_lc(IO_stringValue(line,myPos,2_pInt)) == 'part' ) inPart = .false. + + if ( (inPart .or. noPart) .and. IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*elset' ) then + elemSet = elemSet + 1_pInt + mesh_nameElemSet(elemSet) = trim(IO_extractValue(IO_lc(IO_stringValue(line,myPos,2_pInt)),'elset')) + mesh_mapElemSet(:,elemSet) = IO_continuousIntValues(myUnit,mesh_Nelems,mesh_nameElemSet,& + mesh_mapElemSet,elemSet-1_pInt) + endif + enddo + +640 do i = 1_pInt,elemSet + if (mesh_mapElemSet(1,i) == 0_pInt) call IO_error(error_ID=904_pInt,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(myUnit) + + use IO, only: IO_lc, & + IO_stringValue, & + IO_stringPos, & + IO_extractValue, & + IO_error + + implicit none + integer(pInt), intent(in) :: myUnit + + integer(pInt), parameter :: maxNchunks = 20_pInt + integer(pInt), dimension (1_pInt+2_pInt*maxNchunks) :: myPos + character(len=300) line + + integer(pInt) :: i,c = 0_pInt + logical :: inPart = .false. + character(len=64) :: elemSetName,materialName + + allocate (mesh_nameMaterial(mesh_Nmaterials)) ; mesh_nameMaterial = '' + allocate (mesh_mapMaterial(mesh_Nmaterials)) ; mesh_mapMaterial = '' + +610 FORMAT(A300) + + rewind(myUnit) + do + read (myUnit,610,END=620) line + myPos = IO_stringPos(line,maxNchunks) + if ( IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*part' ) inPart = .true. + if ( IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*end' .and. & + IO_lc(IO_stringValue(line,myPos,2_pInt)) == 'part' ) inPart = .false. + + if ( (inPart .or. noPart) .and. & + IO_lc(IO_StringValue(line,myPos,1_pInt)) == '*solid' .and. & + IO_lc(IO_StringValue(line,myPos,2_pInt)) == 'section' ) then + + elemSetName = '' + materialName = '' + + do i = 3_pInt,myPos(1_pInt) + if (IO_extractValue(IO_lc(IO_stringValue(line,myPos,i)),'elset') /= '') & + elemSetName = trim(IO_extractValue(IO_lc(IO_stringValue(line,myPos,i)),'elset')) + if (IO_extractValue(IO_lc(IO_stringValue(line,myPos,i)),'material') /= '') & + materialName = trim(IO_extractValue(IO_lc(IO_stringValue(line,myPos,i)),'material')) + enddo + + if (elemSetName /= '' .and. materialName /= '') then + c = c + 1_pInt + mesh_nameMaterial(c) = materialName ! name of material used for this section + mesh_mapMaterial(c) = elemSetName ! mapped to respective element set + endif + endif + enddo + +620 if (c==0_pInt) call IO_error(error_ID=905_pInt) + do i=1_pInt,c + if (mesh_nameMaterial(i)=='' .or. mesh_mapMaterial(i)=='') call IO_error(error_ID=905_pInt) + 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(myUnit) + + use IO, only: IO_lc, & + IO_stringValue, & + IO_stringPos, & + IO_error, & + IO_extractValue + + implicit none + integer(pInt), intent(in) :: myUnit + + integer(pInt), parameter :: maxNchunks = 2_pInt + integer(pInt), dimension (1+2*maxNchunks) :: myPos + character(len=300) line + integer(pInt) :: i,k + logical :: materialFound = .false. + character(len=64) ::materialName,elemSetName + + mesh_NcpElems = 0_pInt + +610 FORMAT(A300) + + rewind(myUnit) + do + read (myUnit,610,END=620) line + myPos = IO_stringPos(line,maxNchunks) + select case ( IO_lc(IO_stringValue(line,myPos,1_pInt)) ) + case('*material') + materialName = trim(IO_extractValue(IO_lc(IO_stringValue(line,myPos,2_pInt)),'name')) ! extract name=value + materialFound = materialName /= '' ! valid name? + case('*user') + if (IO_lc(IO_StringValue(line,myPos,2_pInt)) == 'material' .and. materialFound) then + do i = 1_pInt,mesh_Nmaterials ! look thru material names + if (materialName == mesh_nameMaterial(i)) then ! found one + elemSetName = mesh_mapMaterial(i) ! take corresponding elemSet + do k = 1_pInt,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 + +620 if (mesh_NcpElems == 0_pInt) call IO_error(error_ID=906_pInt) + +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(myUnit) + + use math, only: qsort + use IO, only: IO_lc, & + IO_stringValue, & + IO_stringPos, & + IO_extractValue, & + IO_error + + implicit none + integer(pInt), intent(in) :: myUnit + + integer(pInt), parameter :: maxNchunks = 2_pInt + integer(pInt), dimension (1_pInt+2_pInt*maxNchunks) :: myPos + character(len=300) :: line + integer(pInt) ::i,j,k,cpElem = 0_pInt + logical :: materialFound = .false. + character (len=64) materialName,elemSetName ! why limited to 64? ABAQUS? + + allocate (mesh_mapFEtoCPelem(2,mesh_NcpElems)) ; mesh_mapFEtoCPelem = 0_pInt + +610 FORMAT(A300) + + rewind(myUnit) + do + read (myUnit,610,END=660) line + myPos = IO_stringPos(line,maxNchunks) + select case ( IO_lc(IO_stringValue(line,myPos,1_pInt)) ) + case('*material') + materialName = trim(IO_extractValue(IO_lc(IO_stringValue(line,myPos,2_pInt)),'name')) ! extract name=value + materialFound = materialName /= '' ! valid name? + case('*user') + if (IO_lc(IO_stringValue(line,myPos,2_pInt)) == 'material' .and. materialFound) then + do i = 1_pInt,mesh_Nmaterials ! look thru material names + if (materialName == mesh_nameMaterial(i)) then ! found one + elemSetName = mesh_mapMaterial(i) ! take corresponding elemSet + do k = 1_pInt,mesh_NelemSets ! look thru all elemSet definitions + if (elemSetName == mesh_nameElemSet(k)) then ! matched? + do j = 1_pInt,mesh_mapElemSet(1,k) + cpElem = cpElem + 1_pInt + mesh_mapFEtoCPelem(1,cpElem) = mesh_mapElemSet(1_pInt+j,k) ! store FE id + mesh_mapFEtoCPelem(2,cpElem) = cpElem ! store our id + enddo + endif + enddo + endif + enddo + materialFound = .false. + endif + endselect + enddo + +660 call qsort(mesh_mapFEtoCPelem,1_pInt,int(size(mesh_mapFEtoCPelem,2_pInt),pInt)) ! should be mesh_NcpElems + + if (int(size(mesh_mapFEtoCPelem),pInt) < 2_pInt) call IO_error(error_ID=907_pInt) + +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(myUnit) + + use math, only: qsort + use IO, only: IO_lc, & + IO_stringValue, & + IO_stringPos, & + IO_countDataLines, & + IO_intValue, & + IO_error + + implicit none + integer(pInt), intent(in) :: myUnit + + integer(pInt), parameter :: maxNchunks = 2_pInt + integer(pInt), dimension (1_pInt+2_pInt*maxNchunks) :: myPos + character(len=300) line + + integer(pInt) :: i,c,cpNode = 0_pInt + logical :: inPart = .false. + + allocate (mesh_mapFEtoCPnode(2_pInt,mesh_Nnodes)) ; mesh_mapFEtoCPnode = 0_pInt + +610 FORMAT(A300) + + rewind(myUnit) + do + read (myUnit,610,END=650) line + myPos = IO_stringPos(line,maxNchunks) + if ( IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*part' ) inPart = .true. + if ( IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*end' .and. & + IO_lc(IO_stringValue(line,myPos,2_pInt)) == 'part' ) inPart = .false. + + if( (inPart .or. noPart) .and. & + IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*node' .and. & + ( IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'output' .and. & + IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'print' .and. & + IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'file' .and. & + IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'response' ) & + ) then + c = IO_countDataLines(myUnit) + do i = 1_pInt,c + backspace(myUnit) + enddo + do i = 1_pInt,c + read (myUnit,610,END=650) line + myPos = IO_stringPos(line,maxNchunks) + cpNode = cpNode + 1_pInt + mesh_mapFEtoCPnode(1_pInt,cpNode) = IO_intValue(line,myPos,1_pInt) + mesh_mapFEtoCPnode(2_pInt,cpNode) = cpNode + enddo + endif + enddo + +650 call qsort(mesh_mapFEtoCPnode,1_pInt,int(size(mesh_mapFEtoCPnode,2_pInt),pInt)) + + if (int(size(mesh_mapFEtoCPnode),pInt) == 0_pInt) call IO_error(error_ID=908_pInt) + +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(myUnit) + + use IO, only: IO_lc, & + IO_stringValue, & + IO_floatValue, & + IO_stringPos, & + IO_error, & + IO_countDataLines, & + IO_intValue + + implicit none + integer(pInt), intent(in) :: myUnit + + integer(pInt), parameter :: maxNchunks = 4_pInt + integer(pInt), dimension (1_pInt+2_pInt*maxNchunks) :: myPos + character(len=300) :: line + integer(pInt) :: i,j,m,c + logical :: inPart + + allocate ( mesh_node0 (3,mesh_Nnodes) ); mesh_node0 = 0.0_pReal + allocate ( mesh_node (3,mesh_Nnodes) ); mesh_node = 0.0_pReal + +610 FORMAT(A300) + + inPart = .false. + rewind(myUnit) + do + read (myUnit,610,END=670) line + myPos = IO_stringPos(line,maxNchunks) + if ( IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*part' ) inPart = .true. + if ( IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*end' .and. & + IO_lc(IO_stringValue(line,myPos,2_pInt)) == 'part' ) inPart = .false. + + if( (inPart .or. noPart) .and. & + IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*node' .and. & + ( IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'output' .and. & + IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'print' .and. & + IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'file' .and. & + IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'response' ) & + ) then + c = IO_countDataLines(myUnit) ! how many nodes are defined here? + do i = 1_pInt,c + backspace(myUnit) ! rewind to first entry + enddo + do i = 1_pInt,c + read (myUnit,610,END=670) line + myPos = IO_stringPos(line,maxNchunks) + m = mesh_FEasCP('node',IO_intValue(line,myPos,1_pInt)) + forall (j=1_pInt:3_pInt) mesh_node0(j,m) = IO_floatValue(line,myPos,j+1_pInt) + enddo + endif + enddo + +670 if (int(size(mesh_node0,2_pInt),pInt) /= mesh_Nnodes) call IO_error(error_ID=909_pInt) + mesh_node = mesh_node0 + +end subroutine mesh_abaqus_build_nodes + + +!-------------------------------------------------------------------------------------------------- +!> @brief Gets maximum count of nodes, IPs, IP neighbors, and subNodes among cpElements. +!! Allocates global arrays 'mesh_maxNnodes', 'mesh_maxNips', mesh_maxNipNeighbors', +!! and mesh_maxNsubNodes +!-------------------------------------------------------------------------------------------------- +subroutine mesh_abaqus_count_cpSizes(myUnit) + + use IO, only: IO_lc, & + IO_stringValue, & + IO_stringPos, & + IO_extractValue ,& + IO_error, & + IO_countDataLines, & + IO_intValue + + implicit none + integer(pInt), intent(in) :: myUnit + + integer(pInt), parameter :: maxNchunks = 2_pInt + integer(pInt), dimension (1_pInt+2_pInt*maxNchunks) :: myPos + character(len=300) :: line + integer(pInt) :: i,c,t + logical :: inPart + + mesh_maxNnodes = 0_pInt + mesh_maxNips = 0_pInt + mesh_maxNipNeighbors = 0_pInt + mesh_maxNsubNodes = 0_pInt + +610 FORMAT(A300) + + inPart = .false. + rewind(myUnit) + do + read (myUnit,610,END=620) line + myPos = IO_stringPos(line,maxNchunks) + if ( IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*part' ) inPart = .true. + if ( IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*end' .and. & + IO_lc(IO_stringValue(line,myPos,2_pInt)) == 'part' ) inPart = .false. + + if( (inPart .or. noPart) .and. & + IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*element' .and. & + ( IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'output' .and. & + IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'matrix' .and. & + IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'response' ) & + ) then + t = FE_mapElemtype(IO_extractValue(IO_lc(IO_stringValue(line,myPos,2_pInt)),'type')) ! remember elem type + if (t==0_pInt) call IO_error(error_ID=910_pInt,ext_msg='mesh_abaqus_count_cpSizes') + c = IO_countDataLines(myUnit) + do i = 1_pInt,c + backspace(myUnit) + enddo + do i = 1_pInt,c + read (myUnit,610,END=620) line + myPos = IO_stringPos(line,maxNchunks) ! limit to 64 nodes max + if (mesh_FEasCP('elem',IO_intValue(line,myPos,1_pInt)) /= 0_pInt) then ! disregard non CP elems + mesh_maxNnodes = max(mesh_maxNnodes,FE_Nnodes(t)) + mesh_maxNips = max(mesh_maxNips,FE_Nips(t)) + mesh_maxNipNeighbors = max(mesh_maxNipNeighbors,FE_NipNeighbors(t)) + mesh_maxNsubNodes = max(mesh_maxNsubNodes,FE_NsubNodes(t)) + endif + enddo + endif + enddo + +620 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(myUnit) + + use IO, only: IO_lc, & + IO_stringValue, & + IO_skipChunks, & + IO_stringPos, & + IO_intValue, & + IO_extractValue, & + IO_floatValue, & + IO_error, & + IO_countDataLines + + implicit none + integer(pInt), intent(in) :: myUnit + + integer(pInt), parameter :: maxNchunks = 65_pInt + integer(pInt), dimension (1_pInt+2_pInt*maxNchunks) :: myPos + + integer(pInt) :: i,j,k,c,e,t,homog,micro + logical inPart,materialFound + character (len=64) :: materialName,elemSetName + character(len=300) :: line + + allocate (mesh_element (4_pInt+mesh_maxNnodes,mesh_NcpElems)) ; mesh_element = 0_pInt + +610 FORMAT(A300) + + inPart = .false. + rewind(myUnit) + do + read (myUnit,610,END=620) line + myPos(1:1+2*2) = IO_stringPos(line,2_pInt) + if ( IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*part' ) inPart = .true. + if ( IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*end' .and. & + IO_lc(IO_stringValue(line,myPos,2_pInt)) == 'part' ) inPart = .false. + + if( (inPart .or. noPart) .and. & + IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*element' .and. & + ( IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'output' .and. & + IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'matrix' .and. & + IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'response' ) & + ) then + t = FE_mapElemtype(IO_extractValue(IO_lc(IO_stringValue(line,myPos,2_pInt)),'type')) ! remember elem type + if (t==0_pInt) call IO_error(error_ID=910_pInt,ext_msg='mesh_abaqus_build_elements') + c = IO_countDataLines(myUnit) + do i = 1_pInt,c + backspace(myUnit) + enddo + do i = 1_pInt,c + read (myUnit,610,END=620) line + myPos = IO_stringPos(line,maxNchunks) ! limit to 64 nodes max + e = mesh_FEasCP('elem',IO_intValue(line,myPos,1_pInt)) + if (e /= 0_pInt) then ! disregard non CP elems + mesh_element(1,e) = IO_intValue(line,myPos,1_pInt) ! FE id + mesh_element(2,e) = t ! elem type + forall (j=1_pInt:FE_Nnodes(t)) & + mesh_element(4_pInt+j,e) = IO_intValue(line,myPos,1_pInt+j) ! copy FE ids of nodes to position 5: + call IO_skipChunks(myUnit,FE_NoriginalNodes(t)-(myPos(1_pInt)-1_pInt)) ! read on (even multiple lines) if FE_NoriginalNodes exceeds required node count + endif + enddo + endif + enddo + + +620 rewind(myUnit) ! just in case "*material" definitions apear before "*element" + + materialFound = .false. + do + read (myUnit,610,END=630) line + myPos = IO_stringPos(line,maxNchunks) + select case ( IO_lc(IO_StringValue(line,myPos,1_pInt))) + case('*material') + materialName = trim(IO_extractValue(IO_lc(IO_StringValue(line,myPos,2_pInt)),'name')) ! extract name=value + materialFound = materialName /= '' ! valid name? + case('*user') + if ( IO_lc(IO_StringValue(line,myPos,2_pInt)) == 'material' .and. & + materialFound ) then + read (myUnit,610,END=630) line ! read homogenization and microstructure + myPos(1:1+2*2) = IO_stringPos(line,2_pInt) + homog = nint(IO_floatValue(line,myPos,1_pInt),pInt) + micro = nint(IO_floatValue(line,myPos,2_pInt),pInt) + do i = 1_pInt,mesh_Nmaterials ! look thru material names + if (materialName == mesh_nameMaterial(i)) then ! found one + elemSetName = mesh_mapMaterial(i) ! take corresponding elemSet + do k = 1_pInt,mesh_NelemSets ! look thru all elemSet definitions + if (elemSetName == mesh_nameElemSet(k)) then ! matched? + do j = 1_pInt,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 + +630 end subroutine mesh_abaqus_build_elements +#endif + + +!******************************************************************** +! get any additional damask options from input file +! +! mesh_periodicSurface +!******************************************************************** +subroutine mesh_get_damaskOptions(myUnit) + +use IO, only: & + IO_lc, & + IO_stringValue, & + IO_stringPos + +implicit none +integer(pInt), intent(in) :: myUnit + +integer(pInt), parameter :: maxNchunks = 5_pInt +integer(pInt), dimension (1+2*maxNchunks) :: myPos +integer(pInt) chunk, Nchunks +character(len=300) line, keyword, damaskOption, v + +mesh_periodicSurface = .false. + +610 FORMAT(A300) + +#ifdef Marc + keyword = '$damask' +#endif +#ifdef Abaqus + keyword = '**damask' +#endif + +rewind(myUnit) +do + read (myUnit,610,END=620) line + myPos = IO_stringPos(line,maxNchunks) + Nchunks = myPos(1) +#ifndef Spectral + if (IO_lc(IO_stringValue(line,myPos,1_pInt)) == keyword .and. Nchunks > 1_pInt) then ! found keyword for damask option and there is at least one more chunk to read + damaskOption = IO_lc(IO_stringValue(line,myPos,2_pInt)) + select case(damaskOption) + case('periodic') ! damask Option that allows to specify periodic fluxes + do chunk = 3_pInt,Nchunks ! loop through chunks (skipping the keyword) + v = IO_lc(IO_stringValue(line,myPos,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 +#else + damaskOption = IO_lc(IO_stringValue(line,myPos,1_pInt)) + select case(damaskOption) + case('periodic') ! damask Option that allows to specify periodic fluxes + do chunk = 2_pInt,Nchunks ! loop through chunks (skipping the keyword) + v = IO_lc(IO_stringValue(line,myPos,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 + +620 end subroutine mesh_get_damaskOptions + + +!*********************************************************** +! calculation of IP interface areas +! +! allocate globals +! _ipArea, _ipAreaNormal +!*********************************************************** +subroutine mesh_build_ipAreas + + use math, only: math_vectorproduct + + implicit none + integer(pInt) :: e,f,t,i,j,n + integer(pInt), parameter :: Ntriangles = FE_NipFaceNodes-2_pInt ! each interface is made up of this many triangles + real(pReal), dimension (3,FE_NipFaceNodes) :: nPos ! coordinates of nodes on IP face + real(pReal), dimension(3,Ntriangles,FE_NipFaceNodes) :: normal + real(pReal), dimension(Ntriangles,FE_NipFaceNodes) :: area + + allocate(mesh_ipArea(mesh_maxNipNeighbors,mesh_maxNips,mesh_NcpElems)) ; mesh_ipArea = 0.0_pReal + allocate(mesh_ipAreaNormal(3_pInt,mesh_maxNipNeighbors,mesh_maxNips,mesh_NcpElems)) ; mesh_ipAreaNormal = 0.0_pReal + do e = 1_pInt,mesh_NcpElems ! loop over cpElems + t = mesh_element(2,e) ! get elemType + do i = 1_pInt,FE_Nips(t) ! loop over IPs of elem + do f = 1_pInt,FE_NipNeighbors(t) ! loop over interfaces of IP + forall (n = 1_pInt:FE_NipFaceNodes) nPos(:,n) = mesh_subNodeCoord(:,FE_subNodeOnIPFace(n,f,i,t),e) + forall (n = 1_pInt:FE_NipFaceNodes, j = 1_pInt:Ntriangles) ! start at each interface node and build valid triangles to cover interface + normal(:,j,n) = math_vectorproduct(nPos(:,1_pInt+mod(n+j-1_pInt,FE_NipFaceNodes)) - nPos(:,n), & ! calc their normal vectors + nPos(:,1_pInt+mod(n+j-0_pInt,FE_NipFaceNodes)) - nPos(:,n)) + area(j,n) = sqrt(sum(normal(:,j,n)*normal(:,j,n))) ! and area + end forall + forall (n = 1_pInt:FE_NipFaceNodes, j = 1_pInt:Ntriangles, area(j,n) > 0.0_pReal) & + normal(1:3,j,n) = normal(1:3,j,n) / area(j,n) ! make myUnit normal + + mesh_ipArea(f,i,e) = sum(area) / (FE_NipFaceNodes*2.0_pReal) ! area of parallelograms instead of triangles + mesh_ipAreaNormal(:,f,i,e) = sum(sum(normal,3),2_pInt)/& ! average of all valid normals + real(count(area > 0.0_pReal),pReal) + enddo + enddo + enddo + + end subroutine mesh_build_ipAreas + + +!*********************************************************** +! assignment of twin nodes for each cp node +! +! allocate globals +! _nodeTwins +!*********************************************************** +subroutine mesh_build_nodeTwins + +implicit none +integer(pInt) dir, & ! direction of periodicity + node, & + minimumNode, & + maximumNode, & + n1, & + n2 +integer(pInt), 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_pInt + +tolerance = 0.001_pReal * minval(mesh_ipVolume) ** 0.333_pReal + +do dir = 1_pInt,3_pInt ! 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_pInt + maximumNodes = 0_pInt + minCoord = minval(mesh_node0(dir,:)) + maxCoord = maxval(mesh_node0(dir,:)) + do node = 1_pInt,mesh_Nnodes ! loop through all nodes and find surface nodes + if (abs(mesh_node0(dir,node) - minCoord) <= tolerance) then + minimumNodes(1) = minimumNodes(1) + 1_pInt + minimumNodes(minimumNodes(1)+1_pInt) = node + elseif (abs(mesh_node0(dir,node) - maxCoord) <= tolerance) then + maximumNodes(1) = maximumNodes(1) + 1_pInt + maximumNodes(maximumNodes(1)+1_pInt) = node + endif + enddo + + + !*** find the corresponding node on the other side with the same position in this dimension + + unpaired = .true. + do n1 = 1_pInt,minimumNodes(1) + minimumNode = minimumNodes(n1+1_pInt) + if (unpaired(minimumNode)) then + do n2 = 1_pInt,maximumNodes(1) + maximumNode = maximumNodes(n2+1_pInt) + 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 + + + + + + + + + + +!******************************************************************** +! get maximum count of shared elements among cpElements and +! build list of elements shared by each node in mesh +! +! _maxNsharedElems +! _sharedElem +!******************************************************************** +subroutine mesh_build_sharedElems + +implicit none +integer(pint) e, & ! element index + t, & ! element type + node, & ! CP node index + j, & ! node index per element + myDim, & ! dimension index + nodeTwin ! node twin in the specified dimension +integer(pInt), dimension (mesh_Nnodes) :: node_count +integer(pInt), dimension (:), allocatable :: node_seen + +allocate(node_seen(maxval(FE_Nnodes))) + + +node_count = 0_pInt + +do e = 1_pInt,mesh_NcpElems + t = mesh_element(2,e) ! get element type + + node_seen = 0_pInt ! reset node duplicates + do j = 1_pInt,FE_Nnodes(t) ! check each node of element + node = mesh_FEasCP('node',mesh_element(4+j,e)) ! translate to internal (consecutive) numbering + if (all(node_seen /= node)) then + node_count(node) = node_count(node) + 1_pInt ! if FE node not yet encountered -> count it + do myDim = 1_pInt,3_pInt ! check in each dimension... + nodeTwin = mesh_nodeTwins(myDim,node) + if (nodeTwin > 0_pInt) & ! if I am a twin of some node... + node_count(nodeTwin) = node_count(nodeTwin) + 1_pInt ! -> count me again for the twin node + enddo + endif + node_seen(j) = 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)) +mesh_sharedElem = 0_pInt + +do e = 1_pInt,mesh_NcpElems + t = mesh_element(2,e) + node_seen = 0_pInt + do j = 1_pInt,FE_Nnodes(t) + node = mesh_FEasCP('node',mesh_element(4_pInt+j,e)) + if (all(node_seen /= node)) then + mesh_sharedElem(1,node) = mesh_sharedElem(1,node) + 1_pInt ! count for each node the connected elements + mesh_sharedElem(mesh_sharedElem(1,node)+1_pInt,node) = e ! store the respective element id + do myDim = 1_pInt,3_pInt ! check in each dimension... + nodeTwin = mesh_nodeTwins(myDim,node) + if (nodeTwin > 0_pInt) then ! if i am a twin of some node... + mesh_sharedElem(1,nodeTwin) = mesh_sharedElem(1,nodeTwin) + 1_pInt ! ...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(j) = node + enddo +enddo + +deallocate(node_seen) + +end subroutine mesh_build_sharedElems + + +!*********************************************************** +! build up of IP neighborhood +! +! allocate globals +! _ipNeighborhood +!*********************************************************** +subroutine mesh_build_ipNeighborhood + +implicit none +integer(pInt) 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 +integer(pInt), dimension(FE_maxmaxNnodesAtIP) :: & + linkedNodes = 0_pInt, & + matchingNodes +logical checkTwins + +allocate(mesh_ipNeighborhood(2,mesh_maxNipNeighbors,mesh_maxNips,mesh_NcpElems)) +mesh_ipNeighborhood = 0_pInt + + +do myElem = 1_pInt,mesh_NcpElems ! loop over cpElems + myType = mesh_element(2,myElem) ! get elemType + do myIP = 1_pInt,FE_Nips(myType) ! loop over IPs of elem + + do neighbor = 1_pInt,FE_NipNeighbors(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_pInt) 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_pInt) 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_pInt) then ! found match? + neighboringType = mesh_element(2,matchingElem) + + !*** trivial solution if neighbor has only one IP + + if (FE_Nips(neighboringType) == 1_pInt) then + mesh_ipNeighborhood(1,neighbor,myIP,myElem) = matchingElem + mesh_ipNeighborhood(2,neighbor,myIP,myElem) = 1_pInt + cycle + endif + + !*** find those nodes which build the link to the neighbor + + NlinkedNodes = 0_pInt + linkedNodes = 0_pInt + do a = 1_pInt,FE_maxNnodesAtIP(myType) ! figure my anchor nodes on connecting face + anchor = FE_nodesAtIP(a,myIP,myType) + if (anchor /= 0_pInt) then ! valid anchor node + if (any(FE_nodeOnFace(:,myFace,myType) == anchor)) then ! ip anchor sits on face? + NlinkedNodes = NlinkedNodes + 1_pInt + linkedNodes(NlinkedNodes) = & + mesh_FEasCP('node',mesh_element(4_pInt+anchor,myElem)) ! CP id of anchor node + else ! something went wrong with the linkage, since not all anchors sit on my face + NlinkedNodes = 0_pInt + linkedNodes = 0_pInt + 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_pInt,FE_Nips(neighboringType) + NmatchingNodes = 0_pInt + matchingNodes = 0_pInt + do a = 1_pInt,FE_maxNnodesAtIP(neighboringType) ! check each anchor node of that ip + anchor = FE_nodesAtIP(a,candidateIP,neighboringType) + if (anchor /= 0_pInt) then ! valid anchor node + if (any(FE_nodeOnFace(:,matchingFace,neighboringType) == anchor)) then ! sits on matching face? + NmatchingNodes = NmatchingNodes + 1_pInt + matchingNodes(NmatchingNodes) = & + mesh_FEasCP('node',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_pInt + matchingNodes = 0_pInt + 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_pInt,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_pInt,NlinkedNodes + twin_of_linkedNode = mesh_nodeTwins(dir,linkedNodes(a)) + if (twin_of_linkedNode == 0_pInt .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 + +end subroutine mesh_build_ipNeighborhood + + +!*********************************************************** +! write statistics regarding input file parsing +! to the output file +! +!*********************************************************** +subroutine mesh_tell_statistics + + use math, only: math_range + use IO, only: IO_error + use debug, only: debug_what, & + debug_mesh, & + debug_levelBasic, & + debug_levelExtensive, & + debug_levelSelective, & + debug_e, & + debug_i + + implicit none + integer(pInt), dimension (:,:), allocatable :: mesh_HomogMicro + character(len=64) :: myFmt + integer(pInt) :: i,e,n,f,t, myDebug + + myDebug = debug_what(debug_mesh) + + if (mesh_maxValStateVar(1) < 1_pInt) call IO_error(error_ID=170_pInt) ! no homogenization specified + if (mesh_maxValStateVar(2) < 1_pInt) call IO_error(error_ID=180_pInt) ! no microstructure specified + + allocate (mesh_HomogMicro(mesh_maxValStateVar(1),mesh_maxValStateVar(2))); mesh_HomogMicro = 0_pInt +do e = 1_pInt,mesh_NcpElems + if (mesh_element(3,e) < 1_pInt) call IO_error(error_ID=170_pInt,e=e) ! no homogenization specified + if (mesh_element(4,e) < 1_pInt) call IO_error(error_ID=180_pInt,e=e) ! no microstructure specified + mesh_HomogMicro(mesh_element(3,e),mesh_element(4,e)) = & + mesh_HomogMicro(mesh_element(3,e),mesh_element(4,e)) + 1_pInt ! count combinations of homogenization and microstructure +enddo +!$OMP CRITICAL (write2out) + if (iand(myDebug,debug_levelBasic) /= 0_pInt) then + write (6,*) + write (6,*) 'Input Parser: STATISTICS' + write (6,*) + write (6,*) mesh_Nelems, ' : total number of elements in mesh' + write (6,*) mesh_NcpElems, ' : total number of CP elements in mesh' + write (6,*) mesh_Nnodes, ' : total number of nodes in mesh' + write (6,*) mesh_maxNnodes, ' : max number of nodes in any CP element' + write (6,*) mesh_maxNips, ' : max number of IPs in any CP element' + write (6,*) mesh_maxNipNeighbors, ' : max number of IP neighbors in any CP element' + write (6,*) mesh_maxNsubNodes, ' : max number of (additional) subnodes in any CP element' + write (6,*) mesh_maxNsharedElems, ' : max number of CP elements sharing a node' + write (6,*) + write (6,*) 'Input Parser: HOMOGENIZATION/MICROSTRUCTURE' + write (6,*) + write (6,*) mesh_maxValStateVar(1), ' : maximum homogenization index' + write (6,*) mesh_maxValStateVar(2), ' : maximum microstructure index' + write (6,*) + write (myFmt,'(a,i32.32,a)') '(9x,a2,1x,',mesh_maxValStateVar(2),'(i8))' + write (6,myFmt) '+-',math_range(mesh_maxValStateVar(2)) + write (myFmt,'(a,i32.32,a)') '(i8,1x,a2,1x,',mesh_maxValStateVar(2),'(i8))' + do i=1_pInt,mesh_maxValStateVar(1) ! loop over all (possibly assigned) homogenizations + write (6,myFmt) i,'| ',mesh_HomogMicro(i,:) ! loop over all (possibly assigned) microstructures + enddo + write(6,*) + write(6,*) 'Input Parser: ADDITIONAL MPIE OPTIONS' + write(6,*) + write(6,*) 'periodic surface : ', mesh_periodicSurface + write(6,*) + call flush(6) + endif + + if (iand(myDebug,debug_levelExtensive) /= 0_pInt) then + write (6,*) + write (6,*) 'Input Parser: SUBNODE COORDINATES' + write (6,*) + write(6,'(a8,1x,a5,1x,2(a15,1x),a20,3(1x,a12))')& + 'elem','IP','IP neighbor','IPFaceNodes','subNodeOnIPFace','x','y','z' + do e = 1_pInt,mesh_NcpElems ! loop over cpElems + if (iand(myDebug,debug_levelSelective) /= 0_pInt .and. debug_e /= e) cycle + t = mesh_element(2,e) ! get elemType + do i = 1_pInt,FE_Nips(t) ! loop over IPs of elem + if (iand(myDebug,debug_levelSelective) /= 0_pInt .and. debug_i /= i) cycle + do f = 1_pInt,FE_NipNeighbors(t) ! loop over interfaces of IP + do n = 1_pInt,FE_NipFaceNodes ! loop over nodes on interface + write(6,'(i8,1x,i5,2(1x,i15),1x,i20,3(1x,f12.8))') e,i,f,n,FE_subNodeOnIPFace(n,f,i,t),& + mesh_subNodeCoord(1,FE_subNodeOnIPFace(n,f,i,t),e),& + mesh_subNodeCoord(2,FE_subNodeOnIPFace(n,f,i,t),e),& + mesh_subNodeCoord(3,FE_subNodeOnIPFace(n,f,i,t),e) + enddo + enddo + enddo + enddo + write(6,*) + write(6,*) 'Input Parser: IP COORDINATES' + write(6,'(a8,1x,a5,3(1x,a12))') 'elem','IP','x','y','z' + do e = 1_pInt,mesh_NcpElems + if (iand(myDebug,debug_levelSelective) /= 0_pInt .and. debug_e /= e) cycle + do i = 1_pInt,FE_Nips(mesh_element(2,e)) + if (iand(myDebug,debug_levelSelective) /= 0_pInt .and. debug_i /= i) cycle + write (6,'(i8,1x,i5,3(1x,f12.8))') e, i, mesh_ipCenterOfGravity(:,i,e) + enddo + enddo + write (6,*) + write (6,*) 'Input Parser: ELEMENT VOLUME' + write (6,*) + write (6,'(a13,1x,e15.8)') 'total volume', sum(mesh_ipVolume) + write (6,*) + write (6,'(a8,1x,a5,1x,a15,1x,a5,1x,a15,1x,a16)') 'elem','IP','volume','face','area','-- normal --' + do e = 1_pInt,mesh_NcpElems + if (iand(myDebug,debug_levelSelective) /= 0_pInt .and. debug_e /= e) cycle + do i = 1_pInt,FE_Nips(mesh_element(2,e)) + if (iand(myDebug,debug_levelSelective) /= 0_pInt .and. debug_i /= i) cycle + write (6,'(i8,1x,i5,1x,e15.8)') e,i,mesh_IPvolume(i,e) + do f = 1_pInt,FE_NipNeighbors(mesh_element(2,e)) + write (6,'(i33,1x,e15.8,1x,3(f6.3,1x))') f,mesh_ipArea(f,i,e),mesh_ipAreaNormal(:,f,i,e) + enddo + enddo + enddo + write (6,*) + write (6,*) 'Input Parser: NODE TWINS' + write (6,*) + write(6,'(a6,3(3x,a6))') ' node','twin_x','twin_y','twin_z' + do n = 1_pInt,mesh_Nnodes ! loop over cpNodes + if (debug_e <= mesh_NcpElems) then + if (any(mesh_element(5:,debug_e) == n)) then + write(6,'(i6,3(3x,i6))') n, mesh_nodeTwins(1:3,n) + endif + endif + enddo + write(6,*) + write(6,*) 'Input Parser: IP NEIGHBORHOOD' + write(6,*) + write(6,'(a8,1x,a10,1x,a10,1x,a3,1x,a13,1x,a13)') 'elem','IP','neighbor','','elemNeighbor','ipNeighbor' + do e = 1_pInt,mesh_NcpElems ! loop over cpElems + if (iand(myDebug,debug_levelSelective) /= 0_pInt .and. debug_e /= e) cycle + t = mesh_element(2,e) ! get elemType + do i = 1_pInt,FE_Nips(t) ! loop over IPs of elem + if (iand(myDebug,debug_levelSelective) /= 0_pInt .and. debug_i /= i) cycle + do n = 1_pInt,FE_NipNeighbors(t) ! loop over neighbors of IP + write (6,'(i8,1x,i10,1x,i10,1x,a3,1x,i13,1x,i13)') e,i,n,'-->',mesh_ipNeighborhood(1,n,i,e),mesh_ipNeighborhood(2,n,i,e) + enddo + enddo + enddo + endif +!$OMP END CRITICAL (write2out) + + deallocate(mesh_HomogMicro) + +end subroutine mesh_tell_statistics + !*********************************************************** ! mapping of FE element types to internal representation @@ -440,62 +3092,7 @@ integer(pInt) function FE_mapElemtype(what) FE_mapElemtype = 0_pInt ! unknown element --> should raise an error upstream..! endselect - end function FE_mapElemtype - - - -!*********************************************************** -! FE to CP id mapping by binary search thru lookup array -! -! valid questions are 'elem', 'node' -!*********************************************************** -integer(pInt) function mesh_FEasCP(what,myID) - - use IO, only: IO_lc - - implicit none - character(len=*), intent(in) :: what - integer(pInt), intent(in) :: myID - - integer(pInt), dimension(:,:), pointer :: lookupMap - integer(pInt) :: lower,upper,center - - mesh_FEasCP = 0_pInt - select case(IO_lc(what(1:4))) - case('elem') - lookupMap => mesh_mapFEtoCPelem - case('node') - lookupMap => mesh_mapFEtoCPnode - case default - return - endselect - - lower = 1_pInt - upper = int(size(lookupMap,2_pInt),pInt) - - ! check at bounds QUESTION is it valid to extend bounds by 1 and just do binary search w/o init check at bounds? - if (lookupMap(1_pInt,lower) == myID) then - mesh_FEasCP = lookupMap(2_pInt,lower) - return - elseif (lookupMap(1_pInt,upper) == myID) then - mesh_FEasCP = lookupMap(2_pInt,upper) - return - endif - - ! binary search in between bounds - do while (upper-lower > 1_pInt) - center = (lower+upper)/2_pInt - if (lookupMap(1_pInt,center) < myID) then - lower = center - elseif (lookupMap(1_pInt,center) > myID) then - upper = center - else - mesh_FEasCP = lookupMap(2_pInt,center) - exit - endif - enddo - -end function mesh_FEasCP +end function FE_mapElemtype !*********************************************************** @@ -590,7 +3187,7 @@ deallocate(element_seen) end subroutine mesh_faceMatch - + !******************************************************************** ! get properties of different types of finite elements ! @@ -1432,2587 +4029,4 @@ FE_ipNeighbor(1:FE_NipNeighbors(8),1:FE_Nips(8),8) = & ! element 117 end subroutine mesh_build_FEdata -!******************************************************************** -! figure out table styles (Marc only) -! -! initialcondTableStyle, hypoelasticTableStyle -!******************************************************************** -subroutine mesh_marc_get_tableStyles(myUnit) - - use IO, only: & - IO_lc, & - IO_intValue, & - IO_stringValue, & - IO_stringPos - - implicit none - integer(pInt), intent(in) :: myUnit - - integer(pInt), parameter :: maxNchunks = 6_pInt - integer(pInt), dimension (1+2*maxNchunks) :: myPos - character(len=300) line - - initialcondTableStyle = 0_pInt - hypoelasticTableStyle = 0_pInt - -610 FORMAT(A300) - - rewind(myUnit) - do - read (myUnit,610,END=620) line - myPos = IO_stringPos(line,maxNchunks) - - if ( IO_lc(IO_stringValue(line,myPos,1_pInt)) == 'table' .and. myPos(1_pInt) .GT. 5) then - initialcondTableStyle = IO_intValue(line,myPos,4_pInt) - hypoelasticTableStyle = IO_intValue(line,myPos,5_pInt) - exit - endif - enddo - -620 end subroutine mesh_marc_get_tableStyles - - -!******************************************************************** -! get any additional damask options from input file -! -! mesh_periodicSurface -!******************************************************************** -subroutine mesh_get_damaskOptions(myUnit) - -use DAMASK_interface, only: FEsolver -use IO, only: IO_lc, & - IO_stringValue, & - IO_stringPos - -implicit none -integer(pInt), intent(in) :: myUnit - -integer(pInt), parameter :: maxNchunks = 5_pInt -integer(pInt), dimension (1+2*maxNchunks) :: myPos -integer(pInt) chunk, Nchunks -character(len=300) line, keyword, damaskOption, v - -mesh_periodicSurface = .false. - -610 FORMAT(A300) - -select case (FEsolver) - case ('Spectral') ! no special keyword needed, the damask option directly goes into the header - case ('Marc') - keyword = '$damask' - case ('Abaqus') - keyword = '**damask' -end select - -rewind(myUnit) -do - read (myUnit,610,END=620) line - myPos = IO_stringPos(line,maxNchunks) - Nchunks = myPos(1) - select case (FEsolver) - case ('Marc','Abaqus') - if (IO_lc(IO_stringValue(line,myPos,1_pInt)) == keyword .and. Nchunks > 1_pInt) then ! found keyword for damask option and there is at least one more chunk to read - damaskOption = IO_lc(IO_stringValue(line,myPos,2_pInt)) - select case(damaskOption) - case('periodic') ! damask Option that allows to specify periodic fluxes - do chunk = 3_pInt,Nchunks ! loop through chunks (skipping the keyword) - v = IO_lc(IO_stringValue(line,myPos,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 - case('Spectral') - damaskOption = IO_lc(IO_stringValue(line,myPos,1_pInt)) - select case(damaskOption) - case('periodic') ! damask Option that allows to specify periodic fluxes - do chunk = 2_pInt,Nchunks ! loop through chunks (skipping the keyword) - v = IO_lc(IO_stringValue(line,myPos,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 - endselect -enddo - -620 end subroutine mesh_get_damaskOptions - - -!******************************************************************** -! count overall number of nodes and elements in mesh -! -! mesh_Nelems, mesh_Nnodes -!******************************************************************** -subroutine mesh_spectral_count_nodesAndElements(myUnit) - - implicit none - integer(pInt), intent(in) :: myUnit - integer(pInt), dimension(3) :: res - res = mesh_spectral_getResolution(myUnit) - mesh_Nelems = res(1)*res(2)*res(3) - mesh_Nnodes = (1_pInt + res(1))*(1_pInt + res(2))*(1_pInt + res(3)) - -end subroutine mesh_spectral_count_nodesAndElements - -!******************************************************************** -! count overall number of nodes and elements in mesh -! -! mesh_Nelems, mesh_Nnodes -!******************************************************************** -subroutine mesh_marc_count_nodesAndElements(myUnit) - - use IO, only: IO_lc, & - IO_stringValue, & - IO_stringPos, & - IO_IntValue - - implicit none - integer(pInt), intent(in) :: myUnit - - integer(pInt), parameter :: maxNchunks = 4_pInt - integer(pInt), dimension (1+2*maxNchunks) :: myPos - character(len=300) line - - mesh_Nnodes = 0_pInt - mesh_Nelems = 0_pInt - -610 FORMAT(A300) - - rewind(myUnit) - do - read (myUnit,610,END=620) line - myPos = IO_stringPos(line,maxNchunks) - - if ( IO_lc(IO_StringValue(line,myPos,1_pInt)) == 'sizing') then - mesh_Nelems = IO_IntValue (line,myPos,3_pInt) - mesh_Nnodes = IO_IntValue (line,myPos,4_pInt) - exit - endif - enddo - -620 end subroutine mesh_marc_count_nodesAndElements - -!******************************************************************** -! count overall number of nodes and elements in mesh -! -! mesh_Nelems, mesh_Nnodes -!******************************************************************** -subroutine mesh_abaqus_count_nodesAndElements(myUnit) - - use IO, only: IO_lc, & - IO_stringValue, & - IO_stringPos, & - IO_countDataLines, & - IO_error - - implicit none - integer(pInt), intent(in) :: myUnit - - integer(pInt), parameter :: maxNchunks = 2_pInt - integer(pInt), dimension (1+2*maxNchunks) :: myPos - character(len=300) :: line - logical :: inPart - - mesh_Nnodes = 0_pInt - mesh_Nelems = 0_pInt - -610 FORMAT(A300) - - inPart = .false. - rewind(myUnit) - do - read (myUnit,610,END=620) line - myPos = IO_stringPos(line,maxNchunks) - if ( IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*part' ) inPart = .true. - if ( IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*end' .and. & - IO_lc(IO_stringValue(line,myPos,2_pInt)) == 'part' ) inPart = .false. - - if (inPart .or. noPart) then - select case ( IO_lc(IO_stringValue(line,myPos,1_pInt))) - case('*node') - if( & - IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'output' .and. & - IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'print' .and. & - IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'file' .and. & - IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'response' & - ) & - mesh_Nnodes = mesh_Nnodes + IO_countDataLines(myUnit) - case('*element') - if( & - IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'output' .and. & - IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'matrix' .and. & - IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'response' & - ) then - mesh_Nelems = mesh_Nelems + IO_countDataLines(myUnit) - endif - endselect - endif - enddo - -620 if (mesh_Nnodes < 2_pInt) call IO_error(error_ID=900_pInt) - if (mesh_Nelems == 0_pInt) call IO_error(error_ID=901_pInt) - -end subroutine mesh_abaqus_count_nodesAndElements - - -!******************************************************************** -! count overall number of element sets in mesh -! -! mesh_NelemSets, mesh_maxNelemInSet -!******************************************************************** - subroutine mesh_marc_count_elementSets(myUnit) - - use IO, only: IO_lc, & - IO_stringValue, & - IO_stringPos, & - IO_countContinuousIntValues - - implicit none - integer(pInt), intent(in) :: myUnit - - integer(pInt), parameter :: maxNchunks = 2_pInt - integer(pInt), dimension (1+2*maxNchunks) :: myPos - character(len=300) line - - mesh_NelemSets = 0_pInt - mesh_maxNelemInSet = 0_pInt - -610 FORMAT(A300) - - rewind(myUnit) - do - read (myUnit,610,END=620) line - myPos = IO_stringPos(line,maxNchunks) - - if ( IO_lc(IO_StringValue(line,myPos,1_pInt)) == 'define' .and. & - IO_lc(IO_StringValue(line,myPos,2_pInt)) == 'element' ) then - mesh_NelemSets = mesh_NelemSets + 1_pInt - mesh_maxNelemInSet = max(mesh_maxNelemInSet, & - IO_countContinuousIntValues(myUnit)) - endif - enddo - -620 end subroutine mesh_marc_count_elementSets - - -!******************************************************************** -! count overall number of element sets in mesh -! -! mesh_NelemSets, mesh_maxNelemInSet -!******************************************************************** -subroutine mesh_abaqus_count_elementSets(myUnit) - - use IO, only: IO_lc, & - IO_stringValue, & - IO_stringPos, & - IO_error - - implicit none - integer(pInt), intent(in) :: myUnit - - integer(pInt), parameter :: maxNchunks = 2_pInt - integer(pInt), dimension (1+2*maxNchunks) :: myPos - character(len=300) :: line - logical :: inPart - - mesh_NelemSets = 0_pInt - mesh_maxNelemInSet = mesh_Nelems ! have to be conservative, since Abaqus allows for recursive definitons - -610 FORMAT(A300) - - inPart = .false. - rewind(myUnit) - do - read (myUnit,610,END=620) line - myPos = IO_stringPos(line,maxNchunks) - if ( IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*part' ) inPart = .true. - if ( IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*end' .and. & - IO_lc(IO_stringValue(line,myPos,2_pInt)) == 'part' ) inPart = .false. - - if ( (inPart .or. noPart) .and. IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*elset' ) & - mesh_NelemSets = mesh_NelemSets + 1_pInt - enddo - -620 continue - if (mesh_NelemSets == 0) call IO_error(error_ID=902_pInt) - -end subroutine mesh_abaqus_count_elementSets - - -!******************************************************************** -! count overall number of solid sections sets in mesh (Abaqus only) -! -! mesh_Nmaterials -!******************************************************************** -subroutine mesh_abaqus_count_materials(myUnit) - - use IO, only: IO_lc, & - IO_stringValue, & - IO_stringPos, & - IO_error - - implicit none - integer(pInt), intent(in) :: myUnit - - integer(pInt), parameter :: maxNchunks = 2_pInt - integer(pInt), dimension (1_pInt+2_pInt*maxNchunks) :: myPos - character(len=300) :: line - logical inPart - - mesh_Nmaterials = 0_pInt - -610 FORMAT(A300) - - inPart = .false. - rewind(myUnit) - do - read (myUnit,610,END=620) line - myPos = IO_stringPos(line,maxNchunks) - if ( IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*part' ) inPart = .true. - if ( IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*end' .and. & - IO_lc(IO_stringValue(line,myPos,2_pInt)) == 'part' ) inPart = .false. - - if ( (inPart .or. noPart) .and. & - IO_lc(IO_StringValue(line,myPos,1_pInt)) == '*solid' .and. & - IO_lc(IO_StringValue(line,myPos,2_pInt)) == 'section' ) & - mesh_Nmaterials = mesh_Nmaterials + 1_pInt - enddo - -620 if (mesh_Nmaterials == 0_pInt) call IO_error(error_ID=903_pInt) - -end subroutine mesh_abaqus_count_materials - - -!******************************************************************** -! count overall number of cpElements in mesh -! -! mesh_NcpElems -!******************************************************************** -subroutine mesh_spectral_count_cpElements - - implicit none - - mesh_NcpElems = mesh_Nelems - -end subroutine mesh_spectral_count_cpElements - - -!******************************************************************** -! count overall number of cpElements in mesh -! -! mesh_NcpElems -!******************************************************************** -subroutine mesh_marc_count_cpElements(myUnit) - - use IO, only: IO_lc, & - IO_stringValue, & - IO_stringPos, & - IO_countContinuousIntValues - - implicit none - integer(pInt), intent(in) :: myUnit - - integer(pInt), parameter :: maxNchunks = 1_pInt - integer(pInt), dimension (1+2*maxNchunks) :: myPos - integer(pInt) :: i - character(len=300):: line - - mesh_NcpElems = 0_pInt - -610 FORMAT(A300) - - rewind(myUnit) - do - read (myUnit,610,END=620) line - myPos = IO_stringPos(line,maxNchunks) - - if ( IO_lc(IO_stringValue(line,myPos,1_pInt)) == 'hypoelastic') then - do i=1_pInt,3_pInt+hypoelasticTableStyle ! Skip 3 or 4 lines - read (myUnit,610,END=620) line - enddo - mesh_NcpElems = mesh_NcpElems + IO_countContinuousIntValues(myUnit) - exit - endif - enddo - -620 end subroutine mesh_marc_count_cpElements - - -!******************************************************************** -! count overall number of cpElements in mesh -! -! mesh_NcpElems -!******************************************************************** -subroutine mesh_abaqus_count_cpElements(myUnit) - - use IO, only: IO_lc, & - IO_stringValue, & - IO_stringPos, & - IO_error, & - IO_extractValue - - implicit none - integer(pInt), intent(in) :: myUnit - - integer(pInt), parameter :: maxNchunks = 2_pInt - integer(pInt), dimension (1+2*maxNchunks) :: myPos - character(len=300) line - integer(pInt) :: i,k - logical :: materialFound = .false. - character(len=64) ::materialName,elemSetName - - mesh_NcpElems = 0_pInt - -610 FORMAT(A300) - - rewind(myUnit) - do - read (myUnit,610,END=620) line - myPos = IO_stringPos(line,maxNchunks) - select case ( IO_lc(IO_stringValue(line,myPos,1_pInt)) ) - case('*material') - materialName = trim(IO_extractValue(IO_lc(IO_stringValue(line,myPos,2_pInt)),'name')) ! extract name=value - materialFound = materialName /= '' ! valid name? - case('*user') - if (IO_lc(IO_StringValue(line,myPos,2_pInt)) == 'material' .and. materialFound) then - do i = 1_pInt,mesh_Nmaterials ! look thru material names - if (materialName == mesh_nameMaterial(i)) then ! found one - elemSetName = mesh_mapMaterial(i) ! take corresponding elemSet - do k = 1_pInt,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 - -620 if (mesh_NcpElems == 0_pInt) call IO_error(error_ID=906_pInt) - -end subroutine mesh_abaqus_count_cpElements - - -!******************************************************************** -! map element sets -! -! allocate globals: mesh_nameElemSet, mesh_mapElemSet -!******************************************************************** -subroutine mesh_marc_map_elementSets(myUnit) - - use IO, only: IO_lc, & - IO_stringValue, & - IO_stringPos, & - IO_continuousIntValues - - implicit none - integer(pInt), intent(in) :: myUnit - - integer(pInt), parameter :: maxNchunks = 4_pInt - integer(pInt), dimension (1+2*maxNchunks) :: myPos - character(len=300) :: line - integer(pInt) :: elemSet = 0_pInt - - allocate (mesh_nameElemSet(mesh_NelemSets)) ; mesh_nameElemSet = '' - allocate (mesh_mapElemSet(1_pInt+mesh_maxNelemInSet,mesh_NelemSets)) ; mesh_mapElemSet = 0_pInt - -610 FORMAT(A300) - - rewind(myUnit) - do - read (myUnit,610,END=640) line - myPos = IO_stringPos(line,maxNchunks) - if( (IO_lc(IO_stringValue(line,myPos,1_pInt)) == 'define' ) .and. & - (IO_lc(IO_stringValue(line,myPos,2_pInt)) == 'element' ) ) then - elemSet = elemSet+1_pInt - mesh_nameElemSet(elemSet) = trim(IO_stringValue(line,myPos,4_pInt)) - mesh_mapElemSet(:,elemSet) = IO_continuousIntValues(myUnit,mesh_maxNelemInSet,mesh_nameElemSet,mesh_mapElemSet,mesh_NelemSets) - endif - enddo - -640 end subroutine mesh_marc_map_elementSets - - -!******************************************************************** -! Build element set mapping -! -! allocate globals: mesh_nameElemSet, mesh_mapElemSet -!******************************************************************** -subroutine mesh_abaqus_map_elementSets(myUnit) - - use IO, only: IO_lc, & - IO_stringValue, & - IO_stringPos, & - IO_extractValue, & - IO_continuousIntValues, & - IO_error - - implicit none - integer(pInt), intent(in) :: myUnit - - integer(pInt), parameter :: maxNchunks = 4_pInt - integer(pInt), dimension (1_pInt+2_pInt*maxNchunks) :: myPos - character(len=300) :: line - integer(pInt) :: elemSet = 0_pInt,i - logical :: inPart = .false. - - allocate (mesh_nameElemSet(mesh_NelemSets)) ; mesh_nameElemSet = '' - allocate (mesh_mapElemSet(1_pInt+mesh_maxNelemInSet,mesh_NelemSets)) ; mesh_mapElemSet = 0_pInt - -610 FORMAT(A300) - - - rewind(myUnit) - do - read (myUnit,610,END=640) line - myPos = IO_stringPos(line,maxNchunks) - if ( IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*part' ) inPart = .true. - if ( IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*end' .and. & - IO_lc(IO_stringValue(line,myPos,2_pInt)) == 'part' ) inPart = .false. - - if ( (inPart .or. noPart) .and. IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*elset' ) then - elemSet = elemSet + 1_pInt - mesh_nameElemSet(elemSet) = trim(IO_extractValue(IO_lc(IO_stringValue(line,myPos,2_pInt)),'elset')) - mesh_mapElemSet(:,elemSet) = IO_continuousIntValues(myUnit,mesh_Nelems,mesh_nameElemSet,& - mesh_mapElemSet,elemSet-1_pInt) - endif - enddo - -640 do i = 1_pInt,elemSet -! write(6,*)'elemSetName: ',mesh_nameElemSet(i) -! write(6,*)'elems in Elset',mesh_mapElemSet(:,i) - if (mesh_mapElemSet(1,i) == 0_pInt) call IO_error(error_ID=904_pInt,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(myUnit) - - use IO, only: IO_lc, & - IO_stringValue, & - IO_stringPos, & - IO_extractValue, & - IO_error - - implicit none - integer(pInt), intent(in) :: myUnit - - integer(pInt), parameter :: maxNchunks = 20_pInt - integer(pInt), dimension (1_pInt+2_pInt*maxNchunks) :: myPos - character(len=300) line - - integer(pInt) :: i,c = 0_pInt - logical :: inPart = .false. - character(len=64) :: elemSetName,materialName - - allocate (mesh_nameMaterial(mesh_Nmaterials)) ; mesh_nameMaterial = '' - allocate (mesh_mapMaterial(mesh_Nmaterials)) ; mesh_mapMaterial = '' - -610 FORMAT(A300) - - rewind(myUnit) - do - read (myUnit,610,END=620) line - myPos = IO_stringPos(line,maxNchunks) - if ( IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*part' ) inPart = .true. - if ( IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*end' .and. & - IO_lc(IO_stringValue(line,myPos,2_pInt)) == 'part' ) inPart = .false. - - if ( (inPart .or. noPart) .and. & - IO_lc(IO_StringValue(line,myPos,1_pInt)) == '*solid' .and. & - IO_lc(IO_StringValue(line,myPos,2_pInt)) == 'section' ) then - - elemSetName = '' - materialName = '' - - do i = 3_pInt,myPos(1_pInt) - if (IO_extractValue(IO_lc(IO_stringValue(line,myPos,i)),'elset') /= '') & - elemSetName = trim(IO_extractValue(IO_lc(IO_stringValue(line,myPos,i)),'elset')) - if (IO_extractValue(IO_lc(IO_stringValue(line,myPos,i)),'material') /= '') & - materialName = trim(IO_extractValue(IO_lc(IO_stringValue(line,myPos,i)),'material')) - enddo - - if (elemSetName /= '' .and. materialName /= '') then - c = c + 1_pInt - mesh_nameMaterial(c) = materialName ! name of material used for this section - mesh_mapMaterial(c) = elemSetName ! mapped to respective element set - endif - endif - enddo - -620 if (c==0_pInt) call IO_error(error_ID=905_pInt) - do i=1_pInt,c -! write(6,*)'name of materials: ',i,mesh_nameMaterial(i) -! write(6,*)'name of elemSets: ',i,mesh_mapMaterial(i) - if (mesh_nameMaterial(i)=='' .or. mesh_mapMaterial(i)=='') call IO_error(error_ID=905_pInt) - enddo - - end subroutine mesh_abaqus_map_materials - - - -!******************************************************************** -! map nodes from FE id to internal (consecutive) representation -! -! allocate globals: mesh_mapFEtoCPnode -!******************************************************************** -subroutine mesh_spectral_map_nodes - - implicit none - integer(pInt) :: i - - allocate (mesh_mapFEtoCPnode(2_pInt,mesh_Nnodes)) ; mesh_mapFEtoCPnode = 0_pInt - - forall (i = 1_pInt:mesh_Nnodes) & - mesh_mapFEtoCPnode(1:2,i) = i - -end subroutine mesh_spectral_map_nodes - - - -!******************************************************************** -! map nodes from FE id to internal (consecutive) representation -! -! allocate globals: mesh_mapFEtoCPnode -!******************************************************************** -subroutine mesh_marc_map_nodes(myUnit) - - use math, only: qsort - use IO, only: IO_lc, & - IO_stringValue, & - IO_stringPos, & - IO_fixedIntValue - - implicit none - integer(pInt), intent(in) :: myUnit - - integer(pInt), parameter :: maxNchunks = 1_pInt - integer(pInt), dimension (1_pInt+2_pInt*maxNchunks) :: myPos - character(len=300) line - - integer(pInt), dimension (mesh_Nnodes) :: node_count - integer(pInt) :: i - - allocate (mesh_mapFEtoCPnode(2_pInt,mesh_Nnodes)) ; mesh_mapFEtoCPnode = 0_pInt - -610 FORMAT(A300) - - node_count = 0_pInt - - rewind(myUnit) - do - read (myUnit,610,END=650) line - myPos = IO_stringPos(line,maxNchunks) - if( IO_lc(IO_stringValue(line,myPos,1_pInt)) == 'coordinates' ) then - read (myUnit,610,END=650) line ! skip crap line - do i = 1_pInt,mesh_Nnodes - read (myUnit,610,END=650) line - mesh_mapFEtoCPnode(1_pInt,i) = IO_fixedIntValue (line,[ 0_pInt,10_pInt],1_pInt) - mesh_mapFEtoCPnode(2_pInt,i) = i - enddo - exit - endif - enddo - -650 call qsort(mesh_mapFEtoCPnode,1_pInt,int(size(mesh_mapFEtoCPnode,2_pInt),pInt)) - -end subroutine mesh_marc_map_nodes - - - -!******************************************************************** -! map nodes from FE id to internal (consecutive) representation -! -! allocate globals: mesh_mapFEtoCPnode -!******************************************************************** -subroutine mesh_abaqus_map_nodes(myUnit) - - use math, only: qsort - use IO, only: IO_lc, & - IO_stringValue, & - IO_stringPos, & - IO_countDataLines, & - IO_intValue, & - IO_error - - implicit none - integer(pInt), intent(in) :: myUnit - - integer(pInt), parameter :: maxNchunks = 2_pInt - integer(pInt), dimension (1_pInt+2_pInt*maxNchunks) :: myPos - character(len=300) line - - integer(pInt) :: i,c,cpNode = 0_pInt - logical :: inPart = .false. - - allocate (mesh_mapFEtoCPnode(2_pInt,mesh_Nnodes)) ; mesh_mapFEtoCPnode = 0_pInt - -610 FORMAT(A300) - - rewind(myUnit) - do - read (myUnit,610,END=650) line - myPos = IO_stringPos(line,maxNchunks) - if ( IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*part' ) inPart = .true. - if ( IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*end' .and. & - IO_lc(IO_stringValue(line,myPos,2_pInt)) == 'part' ) inPart = .false. - - if( (inPart .or. noPart) .and. & - IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*node' .and. & - ( IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'output' .and. & - IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'print' .and. & - IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'file' .and. & - IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'response' ) & - ) then - c = IO_countDataLines(myUnit) - do i = 1_pInt,c - backspace(myUnit) - enddo - do i = 1_pInt,c - read (myUnit,610,END=650) line - myPos = IO_stringPos(line,maxNchunks) - cpNode = cpNode + 1_pInt - mesh_mapFEtoCPnode(1_pInt,cpNode) = IO_intValue(line,myPos,1_pInt) - mesh_mapFEtoCPnode(2_pInt,cpNode) = cpNode - enddo - endif - enddo - -650 call qsort(mesh_mapFEtoCPnode,1_pInt,int(size(mesh_mapFEtoCPnode,2_pInt),pInt)) - - if (int(size(mesh_mapFEtoCPnode),pInt) == 0_pInt) call IO_error(error_ID=908_pInt) - -end subroutine mesh_abaqus_map_nodes - - -!******************************************************************** -! map elements from FE id to internal (consecutive) representation -! -! allocate globals: mesh_mapFEtoCPelem -!******************************************************************** -subroutine mesh_spectral_map_elements - - implicit none - integer(pInt) :: i - - allocate (mesh_mapFEtoCPelem(2_pInt,mesh_NcpElems)) ; mesh_mapFEtoCPelem = 0_pInt - - forall (i = 1_pInt:mesh_NcpElems) & - mesh_mapFEtoCPelem(1:2,i) = i - -end subroutine mesh_spectral_map_elements - - - -!******************************************************************** -! map elements from FE id to internal (consecutive) representation -! -! allocate globals: mesh_mapFEtoCPelem -!******************************************************************** -subroutine mesh_marc_map_elements(myUnit) - - use math, only: qsort - use IO, only: IO_lc, & - IO_stringValue, & - IO_stringPos, & - IO_continuousIntValues - - implicit none - integer(pInt), intent(in) :: myUnit - - integer(pInt), parameter :: maxNchunks = 1_pInt - integer(pInt), dimension (1_pInt+2_pInt*maxNchunks) :: myPos - character(len=300) line - - integer(pInt), dimension (1_pInt+mesh_NcpElems) :: contInts - integer(pInt) :: i,cpElem = 0_pInt - - allocate (mesh_mapFEtoCPelem(2,mesh_NcpElems)) ; mesh_mapFEtoCPelem = 0_pInt - -610 FORMAT(A300) - - rewind(myUnit) - do - read (myUnit,610,END=660) line - myPos = IO_stringPos(line,maxNchunks) - if( IO_lc(IO_stringValue(line,myPos,1_pInt)) == 'hypoelastic' ) then - do i=1_pInt,3_pInt+hypoelasticTableStyle ! skip three (or four if new table style!) lines - read (myUnit,610,END=660) line - enddo - contInts = IO_continuousIntValues(myUnit,mesh_NcpElems,mesh_nameElemSet,mesh_mapElemSet,mesh_NelemSets) - do i = 1_pInt,contInts(1) - cpElem = cpElem+1_pInt - mesh_mapFEtoCPelem(1,cpElem) = contInts(1_pInt+i) - mesh_mapFEtoCPelem(2,cpElem) = cpElem - enddo - endif - enddo - -660 call qsort(mesh_mapFEtoCPelem,1_pInt,int(size(mesh_mapFEtoCPelem,2_pInt),pInt)) ! should be mesh_NcpElems - -end subroutine mesh_marc_map_elements - - -!******************************************************************** -! map elements from FE id to internal (consecutive) representation -! -! allocate globals: mesh_mapFEtoCPelem -!******************************************************************** -subroutine mesh_abaqus_map_elements(myUnit) - - use math, only: qsort - use IO, only: IO_lc, & - IO_stringValue, & - IO_stringPos, & - IO_extractValue, & - IO_error - - implicit none - integer(pInt), intent(in) :: myUnit - - integer(pInt), parameter :: maxNchunks = 2_pInt - integer(pInt), dimension (1_pInt+2_pInt*maxNchunks) :: myPos - character(len=300) :: line - integer(pInt) ::i,j,k,cpElem = 0_pInt - logical :: materialFound = .false. - character (len=64) materialName,elemSetName ! why limited to 64? ABAQUS? - - allocate (mesh_mapFEtoCPelem(2,mesh_NcpElems)) ; mesh_mapFEtoCPelem = 0_pInt - -610 FORMAT(A300) - - rewind(myUnit) - do - read (myUnit,610,END=660) line - myPos = IO_stringPos(line,maxNchunks) - select case ( IO_lc(IO_stringValue(line,myPos,1_pInt)) ) - case('*material') - materialName = trim(IO_extractValue(IO_lc(IO_stringValue(line,myPos,2_pInt)),'name')) ! extract name=value - materialFound = materialName /= '' ! valid name? - case('*user') - if (IO_lc(IO_stringValue(line,myPos,2_pInt)) == 'material' .and. materialFound) then - do i = 1_pInt,mesh_Nmaterials ! look thru material names - if (materialName == mesh_nameMaterial(i)) then ! found one - elemSetName = mesh_mapMaterial(i) ! take corresponding elemSet - do k = 1_pInt,mesh_NelemSets ! look thru all elemSet definitions - if (elemSetName == mesh_nameElemSet(k)) then ! matched? - do j = 1_pInt,mesh_mapElemSet(1,k) - cpElem = cpElem + 1_pInt - mesh_mapFEtoCPelem(1,cpElem) = mesh_mapElemSet(1_pInt+j,k) ! store FE id - mesh_mapFEtoCPelem(2,cpElem) = cpElem ! store our id - enddo - endif - enddo - endif - enddo - materialFound = .false. - endif - endselect - enddo - -660 call qsort(mesh_mapFEtoCPelem,1_pInt,int(size(mesh_mapFEtoCPelem,2_pInt),pInt)) ! should be mesh_NcpElems - - if (int(size(mesh_mapFEtoCPelem),pInt) < 2_pInt) call IO_error(error_ID=907_pInt) - -end subroutine mesh_abaqus_map_elements - - -!******************************************************************** -! get maximum count of nodes, IPs, IP neighbors, and subNodes -! among cpElements -! -! _maxNnodes, _maxNips, _maxNipNeighbors, _maxNsubNodes -!******************************************************************** -subroutine mesh_spectral_count_cpSizes - - implicit none - integer(pInt) :: t - - t = FE_mapElemtype('C3D8R') ! fake 3D hexahedral 8 node 1 IP element - - mesh_maxNnodes = FE_Nnodes(t) - mesh_maxNips = FE_Nips(t) - mesh_maxNipNeighbors = FE_NipNeighbors(t) - mesh_maxNsubNodes = FE_NsubNodes(t) - -end subroutine mesh_spectral_count_cpSizes - - -!******************************************************************** -! get maximum count of nodes, IPs, IP neighbors, and subNodes -! among cpElements -! -! _maxNnodes, _maxNips, _maxNipNeighbors, _maxNsubNodes -!******************************************************************** -subroutine mesh_marc_count_cpSizes(myUnit) - - use IO, only: IO_lc, & - IO_stringValue, & - IO_stringPos, & - IO_intValue, & - IO_skipChunks - - implicit none - integer(pInt), intent(in) :: myUnit - - integer(pInt), parameter :: maxNchunks = 2_pInt - integer(pInt), dimension (1_pInt+2_pInt*maxNchunks) :: myPos - character(len=300) :: line - integer(pInt) :: i,t,e - - mesh_maxNnodes = 0_pInt - mesh_maxNips = 0_pInt - mesh_maxNipNeighbors = 0_pInt - mesh_maxNsubNodes = 0_pInt - -610 FORMAT(A300) - rewind(myUnit) - do - read (myUnit,610,END=630) line - myPos = IO_stringPos(line,maxNchunks) - if( IO_lc(IO_stringValue(line,myPos,1_pInt)) == 'connectivity' ) then - read (myUnit,610,END=630) line ! Garbage line - do i=1_pInt,mesh_Nelems ! read all elements - read (myUnit,610,END=630) line - myPos = IO_stringPos(line,maxNchunks) ! limit to id and type - e = mesh_FEasCP('elem',IO_intValue(line,myPos,1_pInt)) - if (e /= 0_pInt) then - t = FE_mapElemtype(IO_stringValue(line,myPos,2_pInt)) - mesh_maxNnodes = max(mesh_maxNnodes,FE_Nnodes(t)) - mesh_maxNips = max(mesh_maxNips,FE_Nips(t)) - mesh_maxNipNeighbors = max(mesh_maxNipNeighbors,FE_NipNeighbors(t)) - mesh_maxNsubNodes = max(mesh_maxNsubNodes,FE_NsubNodes(t)) - call IO_skipChunks(myUnit,FE_NoriginalNodes(t)-(myPos(1_pInt)-2_pInt)) ! read on if FE_Nnodes exceeds node count present on current line - endif - enddo - exit - endif - enddo - -630 end subroutine mesh_marc_count_cpSizes - - -!******************************************************************** -! get maximum count of nodes, IPs, IP neighbors, and subNodes -! among cpElements -! -! _maxNnodes, _maxNips, _maxNipNeighbors, _maxNsubNodes -!******************************************************************** -subroutine mesh_abaqus_count_cpSizes(myUnit) - - use IO, only: IO_lc, & - IO_stringValue, & - IO_stringPos, & - IO_extractValue ,& - IO_error, & - IO_countDataLines, & - IO_intValue - - implicit none - integer(pInt), intent(in) :: myUnit - - integer(pInt), parameter :: maxNchunks = 2_pInt - integer(pInt), dimension (1_pInt+2_pInt*maxNchunks) :: myPos - character(len=300) :: line - integer(pInt) :: i,c,t - logical :: inPart - - mesh_maxNnodes = 0_pInt - mesh_maxNips = 0_pInt - mesh_maxNipNeighbors = 0_pInt - mesh_maxNsubNodes = 0_pInt - -610 FORMAT(A300) - - inPart = .false. - rewind(myUnit) - do - read (myUnit,610,END=620) line - myPos = IO_stringPos(line,maxNchunks) - if ( IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*part' ) inPart = .true. - if ( IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*end' .and. & - IO_lc(IO_stringValue(line,myPos,2_pInt)) == 'part' ) inPart = .false. - - if( (inPart .or. noPart) .and. & - IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*element' .and. & - ( IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'output' .and. & - IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'matrix' .and. & - IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'response' ) & - ) then - t = FE_mapElemtype(IO_extractValue(IO_lc(IO_stringValue(line,myPos,2_pInt)),'type')) ! remember elem type - if (t==0_pInt) call IO_error(error_ID=910_pInt,ext_msg='mesh_abaqus_count_cpSizes') - c = IO_countDataLines(myUnit) - do i = 1_pInt,c - backspace(myUnit) - enddo - do i = 1_pInt,c - read (myUnit,610,END=620) line - myPos = IO_stringPos(line,maxNchunks) ! limit to 64 nodes max - if (mesh_FEasCP('elem',IO_intValue(line,myPos,1_pInt)) /= 0_pInt) then ! disregard non CP elems - mesh_maxNnodes = max(mesh_maxNnodes,FE_Nnodes(t)) - mesh_maxNips = max(mesh_maxNips,FE_Nips(t)) - mesh_maxNipNeighbors = max(mesh_maxNipNeighbors,FE_NipNeighbors(t)) - mesh_maxNsubNodes = max(mesh_maxNsubNodes,FE_NsubNodes(t)) - endif - enddo - endif - enddo - -620 end subroutine mesh_abaqus_count_cpSizes - - -!******************************************************************** -! store x,y,z coordinates of all nodes in mesh -! -! allocate globals: -! _node -!******************************************************************** -subroutine mesh_spectral_build_nodes(myUnit) - - use IO, only: & - IO_error - - implicit none - integer(pInt), intent(in) :: myUnit - integer(pInt) :: n - integer(pInt), dimension(3) :: res = 1_pInt - real(pReal), dimension(3) :: geomdim = 1.0_pReal - - allocate ( mesh_node0 (3,mesh_Nnodes) ); mesh_node0 = 0.0_pReal - allocate ( mesh_node (3,mesh_Nnodes) ); mesh_node = 0.0_pReal - - res = mesh_spectral_getResolution(myUnit) - geomdim = mesh_spectral_getDimension(myUnit) - - forall (n = 0_pInt:mesh_Nnodes-1_pInt) - mesh_node0(1,n+1_pInt) = geomdim(1) * real(mod(n,(res(1)+1_pInt) ),pReal) & - / real(res(1),pReal) - mesh_node0(2,n+1_pInt) = geomdim(2) * real(mod(n/(res(1)+1_pInt),(res(2)+1_pInt)),pReal) & - / real(res(2),pReal) - mesh_node0(3,n+1_pInt) = geomdim(3) * real(mod(n/(res(1)+1_pInt)/(res(2)+1_pInt),(res(3)+1_pInt)),pReal) & - / real(res(3),pReal) - end forall - - mesh_node = mesh_node0 !why? - -end subroutine mesh_spectral_build_nodes - - -!******************************************************************** -! store x,y,z coordinates of all nodes in mesh -! -! allocate globals: -! _node -!******************************************************************** -subroutine mesh_marc_build_nodes(myUnit) - - use IO, only: IO_lc, & - IO_stringValue, & - IO_stringPos, & - IO_fixedIntValue, & - IO_fixedNoEFloatValue - - implicit none - integer(pInt), intent(in) :: myUnit - - integer(pInt), dimension(5), parameter :: node_ends = int([0,10,30,50,70],pInt) - integer(pInt), parameter :: maxNchunks = 1_pInt - integer(pInt), dimension (1_pInt+2_pInt*maxNchunks) :: myPos - character(len=300) :: line - integer(pInt) :: i,j,m - - allocate ( mesh_node0 (3,mesh_Nnodes) ); mesh_node0 = 0.0_pReal - allocate ( mesh_node (3,mesh_Nnodes) ); mesh_node = 0.0_pReal - -610 FORMAT(A300) - - rewind(myUnit) - do - read (myUnit,610,END=670) line - myPos = IO_stringPos(line,maxNchunks) - if( IO_lc(IO_stringValue(line,myPos,1_pInt)) == 'coordinates' ) then - read (myUnit,610,END=670) line ! skip crap line - do i=1_pInt,mesh_Nnodes - read (myUnit,610,END=670) line - m = mesh_FEasCP('node',IO_fixedIntValue(line,node_ends,1_pInt)) - forall (j = 1_pInt:3_pInt) mesh_node0(j,m) = IO_fixedNoEFloatValue(line,node_ends,j+1_pInt) - enddo - exit - endif - enddo - -670 mesh_node = mesh_node0 - -end subroutine mesh_marc_build_nodes - - -!******************************************************************** -! store x,y,z coordinates of all nodes in mesh -! -! allocate globals: -! _node -!******************************************************************** -subroutine mesh_abaqus_build_nodes(myUnit) - - use IO, only: IO_lc, & - IO_stringValue, & - IO_floatValue, & - IO_stringPos, & - IO_error, & - IO_countDataLines, & - IO_intValue - - implicit none - integer(pInt), intent(in) :: myUnit - - integer(pInt), parameter :: maxNchunks = 4_pInt - integer(pInt), dimension (1_pInt+2_pInt*maxNchunks) :: myPos - character(len=300) :: line - integer(pInt) :: i,j,m,c - logical :: inPart - - allocate ( mesh_node0 (3,mesh_Nnodes) ); mesh_node0 = 0.0_pReal - allocate ( mesh_node (3,mesh_Nnodes) ); mesh_node = 0.0_pReal - -610 FORMAT(A300) - - inPart = .false. - rewind(myUnit) - do - read (myUnit,610,END=670) line - myPos = IO_stringPos(line,maxNchunks) - if ( IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*part' ) inPart = .true. - if ( IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*end' .and. & - IO_lc(IO_stringValue(line,myPos,2_pInt)) == 'part' ) inPart = .false. - - if( (inPart .or. noPart) .and. & - IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*node' .and. & - ( IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'output' .and. & - IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'print' .and. & - IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'file' .and. & - IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'response' ) & - ) then - c = IO_countDataLines(myUnit) ! how many nodes are defined here? - do i = 1_pInt,c - backspace(myUnit) ! rewind to first entry - enddo - do i = 1_pInt,c - read (myUnit,610,END=670) line - myPos = IO_stringPos(line,maxNchunks) - m = mesh_FEasCP('node',IO_intValue(line,myPos,1_pInt)) - forall (j=1_pInt:3_pInt) mesh_node0(j,m) = IO_floatValue(line,myPos,j+1_pInt) - enddo - endif - enddo - -670 if (int(size(mesh_node0,2_pInt),pInt) /= mesh_Nnodes) call IO_error(error_ID=909_pInt) - mesh_node = mesh_node0 - -end subroutine mesh_abaqus_build_nodes - - -!******************************************************************** -! store FEid, type, mat, tex, and node list per element -! -! allocate globals: -! _element -!******************************************************************** -subroutine mesh_spectral_build_elements(myUnit) - - use IO, only: & - IO_checkAndRewind, & - IO_lc, & - IO_stringValue, & - IO_stringPos, & - IO_error, & - IO_continuousIntValues, & - IO_intValue, & - IO_countContinuousIntValues - - implicit none - integer(pInt), intent(in) :: myUnit - - integer(pInt), parameter :: maxNchunks = 7_pInt - integer(pInt), dimension (1_pInt+2_pInt*maxNchunks) :: myPos - integer(pInt), dimension(3) :: res - integer(pInt) :: e, i, j, homog = 0_pInt, headerLength = 0_pInt, maxIntCount - integer(pInt), dimension(:), allocatable :: microstructures - integer(pInt), dimension(1,1) :: dummySet = 0_pInt - character(len=65536) :: line,keyword - character(len=64), dimension(1) :: dummyName = '' - - res = mesh_spectral_getResolution(myUnit) - homog = mesh_spectral_getHomogenization(myUnit) - - call IO_checkAndRewind(myUnit) - - read(myUnit,'(a65536)') line - myPos = IO_stringPos(line,2_pInt) - keyword = IO_lc(IO_StringValue(line,myPos,2_pInt)) - if (keyword(1:4) == 'head') then - headerLength = IO_intValue(line,myPos,1_pInt) + 1_pInt - else - call IO_error(error_ID=841_pInt, ext_msg='mesh_spectral_build_elements') - endif - - rewind(myUnit) - do i = 1_pInt, headerLength - read(myUnit,'(a65536)') line - enddo - - maxIntCount = 0_pInt - i = 1_pInt - - do while (i > 0_pInt) - i = IO_countContinuousIntValues(myUnit) - maxIntCount = max(maxIntCount, i) - enddo - - rewind (myUnit) - do i=1_pInt,headerLength ! skip header - read(myUnit,'(a65536)') line - enddo - - allocate (mesh_element (4_pInt+mesh_maxNnodes,mesh_NcpElems)) ; mesh_element = 0_pInt - allocate (microstructures (1_pInt+maxIntCount)) ; microstructures = 2_pInt - - e = 0_pInt - do while (e < mesh_NcpElems .and. microstructures(1) > 0_pInt) ! fill expected number of elements, stop at end of data (or blank line!) - microstructures = IO_continuousIntValues(myUnit,maxIntCount,dummyName,dummySet,0_pInt) ! get affected elements - do i = 1_pInt,microstructures(1_pInt) - e = e+1_pInt ! valid element entry - mesh_element( 1,e) = e ! FE id - mesh_element( 2,e) = FE_mapElemtype('C3D8R') ! elem type - mesh_element( 3,e) = homog ! homogenization - mesh_element( 4,e) = microstructures(1_pInt+i) ! microstructure - mesh_element( 5,e) = e + (e-1_pInt)/res(1) + & - ((e-1_pInt)/(res(1)*res(2)))*(res(1)+1_pInt) ! base node - mesh_element( 6,e) = mesh_element(5,e) + 1_pInt - mesh_element( 7,e) = mesh_element(5,e) + res(1) + 2_pInt - mesh_element( 8,e) = mesh_element(5,e) + res(1) + 1_pInt - mesh_element( 9,e) = mesh_element(5,e) +(res(1) + 1_pInt) * (res(2) + 1_pInt) ! second floor base node - mesh_element(10,e) = mesh_element(9,e) + 1_pInt - mesh_element(11,e) = mesh_element(9,e) + res(1) + 2_pInt - mesh_element(12,e) = mesh_element(9,e) + res(1) + 1_pInt - mesh_maxValStateVar(1) = max(mesh_maxValStateVar(1),mesh_element(3,e)) !needed for statistics - mesh_maxValStateVar(2) = max(mesh_maxValStateVar(2),mesh_element(4,e)) - enddo - enddo - - deallocate(microstructures) - if (e /= mesh_NcpElems) call IO_error(880_pInt,e) - -end subroutine mesh_spectral_build_elements - - -!******************************************************************** -! store FEid, type, mat, tex, and node list per element -! -! allocate globals: -! _element -!******************************************************************** -subroutine mesh_marc_build_elements(myUnit) - - use IO, only: IO_lc, & - IO_stringValue, & - IO_fixedNoEFloatValue, & - IO_skipChunks, & - IO_stringPos, & - IO_intValue, & - IO_continuousIntValues - - implicit none - integer(pInt), intent(in) :: myUnit - - integer(pInt), parameter :: maxNchunks = 66_pInt - integer(pInt), dimension (1_pInt+2_pInt*maxNchunks) :: myPos - character(len=300) line - - integer(pInt), dimension(1_pInt+mesh_NcpElems) :: contInts - integer(pInt) :: i,j,sv,myVal,e - - allocate (mesh_element (4_pInt+mesh_maxNnodes,mesh_NcpElems)) ; mesh_element = 0_pInt - -610 FORMAT(A300) - - rewind(myUnit) - do - read (myUnit,610,END=620) line - myPos(1:1+2*1) = IO_stringPos(line,1_pInt) - if( IO_lc(IO_stringValue(line,myPos,1_pInt)) == 'connectivity' ) then - read (myUnit,610,END=620) line ! Garbage line - do i = 1_pInt,mesh_Nelems - read (myUnit,610,END=620) line - myPos = IO_stringPos(line,maxNchunks) ! limit to 64 nodes max (plus ID, type) - e = mesh_FEasCP('elem',IO_intValue(line,myPos,1_pInt)) - if (e /= 0_pInt) then ! disregard non CP elems - mesh_element(1,e) = IO_IntValue (line,myPos,1_pInt) ! FE id - mesh_element(2,e) = FE_mapElemtype(IO_StringValue(line,myPos,2_pInt)) ! elem type - forall (j = 1_pInt:FE_Nnodes(mesh_element(2,e))) & - mesh_element(j+4_pInt,e) = IO_IntValue(line,myPos,j+2_pInt) ! copy FE ids of nodes - call IO_skipChunks(myUnit,FE_NoriginalNodes(mesh_element(2_pInt,e))-(myPos(1_pInt)-2_pInt)) ! read on if FE_Nnodes exceeds node count present on current line - endif - enddo - exit - endif - enddo - -620 rewind(myUnit) ! just in case "initial state" apears before "connectivity" - read (myUnit,610,END=620) line - do - myPos(1:1+2*2) = IO_stringPos(line,2_pInt) - if( (IO_lc(IO_stringValue(line,myPos,1_pInt)) == 'initial') .and. & - (IO_lc(IO_stringValue(line,myPos,2_pInt)) == 'state') ) then - if (initialcondTableStyle == 2_pInt) read (myUnit,610,END=620) line ! read extra line for new style - read (myUnit,610,END=630) line ! read line with index of state var - myPos(1:1+2*1) = IO_stringPos(line,1_pInt) - sv = IO_IntValue(line,myPos,1_pInt) ! figure state variable index - if( (sv == 2_pInt).or.(sv == 3_pInt) ) then ! only state vars 2 and 3 of interest - read (myUnit,610,END=620) line ! read line with value of state var - myPos(1:1+2*1) = IO_stringPos(line,1_pInt) - do while (scan(IO_stringValue(line,myPos,1_pInt),'+-',back=.true.)>1) ! is noEfloat value? - myVal = nint(IO_fixedNoEFloatValue(line,[0_pInt,20_pInt],1_pInt),pInt) ! state var's value - mesh_maxValStateVar(sv-1_pInt) = max(myVal,mesh_maxValStateVar(sv-1_pInt)) ! remember max val of homogenization and microstructure index - if (initialcondTableStyle == 2_pInt) then - read (myUnit,610,END=630) line ! read extra line - read (myUnit,610,END=630) line ! read extra line - endif - contInts = IO_continuousIntValues& ! get affected elements - (myUnit,mesh_Nelems,mesh_nameElemSet,mesh_mapElemSet,mesh_NelemSets) - do i = 1_pInt,contInts(1) - e = mesh_FEasCP('elem',contInts(1_pInt+i)) - mesh_element(1_pInt+sv,e) = myVal - enddo - if (initialcondTableStyle == 0_pInt) read (myUnit,610,END=620) line ! ignore IP range for old table style - read (myUnit,610,END=630) line - myPos(1:1+2*1) = IO_stringPos(line,1_pInt) - enddo - endif - else - read (myUnit,610,END=630) line - endif - enddo - -630 end subroutine mesh_marc_build_elements - -!******************************************************************** -! store FEid, type, mat, tex, and node list per element -! -! allocate globals: -! _element -!******************************************************************** -subroutine mesh_abaqus_build_elements(myUnit) - - use IO, only: IO_lc, & - IO_stringValue, & - IO_skipChunks, & - IO_stringPos, & - IO_intValue, & - IO_extractValue, & - IO_floatValue, & - IO_error, & - IO_countDataLines - - implicit none - integer(pInt), intent(in) :: myUnit - - integer(pInt), parameter :: maxNchunks = 65_pInt - integer(pInt), dimension (1_pInt+2_pInt*maxNchunks) :: myPos - - integer(pInt) :: i,j,k,c,e,t,homog,micro - logical inPart,materialFound - character (len=64) :: materialName,elemSetName - character(len=300) :: line - - allocate (mesh_element (4_pInt+mesh_maxNnodes,mesh_NcpElems)) ; mesh_element = 0_pInt - -610 FORMAT(A300) - - inPart = .false. - rewind(myUnit) - do - read (myUnit,610,END=620) line - myPos(1:1+2*2) = IO_stringPos(line,2_pInt) - if ( IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*part' ) inPart = .true. - if ( IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*end' .and. & - IO_lc(IO_stringValue(line,myPos,2_pInt)) == 'part' ) inPart = .false. - - if( (inPart .or. noPart) .and. & - IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*element' .and. & - ( IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'output' .and. & - IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'matrix' .and. & - IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'response' ) & - ) then - t = FE_mapElemtype(IO_extractValue(IO_lc(IO_stringValue(line,myPos,2_pInt)),'type')) ! remember elem type - if (t==0_pInt) call IO_error(error_ID=910_pInt,ext_msg='mesh_abaqus_build_elements') - c = IO_countDataLines(myUnit) - do i = 1_pInt,c - backspace(myUnit) - enddo - do i = 1_pInt,c - read (myUnit,610,END=620) line - myPos = IO_stringPos(line,maxNchunks) ! limit to 64 nodes max - e = mesh_FEasCP('elem',IO_intValue(line,myPos,1_pInt)) - if (e /= 0_pInt) then ! disregard non CP elems - mesh_element(1,e) = IO_intValue(line,myPos,1_pInt) ! FE id - mesh_element(2,e) = t ! elem type - forall (j=1_pInt:FE_Nnodes(t)) & - mesh_element(4_pInt+j,e) = IO_intValue(line,myPos,1_pInt+j) ! copy FE ids of nodes to position 5: - call IO_skipChunks(myUnit,FE_NoriginalNodes(t)-(myPos(1_pInt)-1_pInt)) ! read on (even multiple lines) if FE_NoriginalNodes exceeds required node count - endif - enddo - endif - enddo - - -620 rewind(myUnit) ! just in case "*material" definitions apear before "*element" - - materialFound = .false. - do - read (myUnit,610,END=630) line - myPos = IO_stringPos(line,maxNchunks) - select case ( IO_lc(IO_StringValue(line,myPos,1_pInt))) - case('*material') - materialName = trim(IO_extractValue(IO_lc(IO_StringValue(line,myPos,2_pInt)),'name')) ! extract name=value - materialFound = materialName /= '' ! valid name? - case('*user') - if ( IO_lc(IO_StringValue(line,myPos,2_pInt)) == 'material' .and. & - materialFound ) then - read (myUnit,610,END=630) line ! read homogenization and microstructure - myPos(1:1+2*2) = IO_stringPos(line,2_pInt) - homog = nint(IO_floatValue(line,myPos,1_pInt),pInt) - micro = nint(IO_floatValue(line,myPos,2_pInt),pInt) - do i = 1_pInt,mesh_Nmaterials ! look thru material names - if (materialName == mesh_nameMaterial(i)) then ! found one - elemSetName = mesh_mapMaterial(i) ! take corresponding elemSet - do k = 1_pInt,mesh_NelemSets ! look thru all elemSet definitions - if (elemSetName == mesh_nameElemSet(k)) then ! matched? - do j = 1_pInt,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 - -630 end subroutine mesh_abaqus_build_elements - - -!******************************************************************** -! get maximum count of shared elements among cpElements and -! build list of elements shared by each node in mesh -! -! _maxNsharedElems -! _sharedElem -!******************************************************************** -subroutine mesh_build_sharedElems - -implicit none -integer(pint) e, & ! element index - t, & ! element type - node, & ! CP node index - j, & ! node index per element - myDim, & ! dimension index - nodeTwin ! node twin in the specified dimension -integer(pInt), dimension (mesh_Nnodes) :: node_count -integer(pInt), dimension (:), allocatable :: node_seen - -allocate(node_seen(maxval(FE_Nnodes))) - - -node_count = 0_pInt - -do e = 1_pInt,mesh_NcpElems - t = mesh_element(2,e) ! get element type - - node_seen = 0_pInt ! reset node duplicates - do j = 1_pInt,FE_Nnodes(t) ! check each node of element - node = mesh_FEasCP('node',mesh_element(4+j,e)) ! translate to internal (consecutive) numbering - if (all(node_seen /= node)) then - node_count(node) = node_count(node) + 1_pInt ! if FE node not yet encountered -> count it - do myDim = 1_pInt,3_pInt ! check in each dimension... - nodeTwin = mesh_nodeTwins(myDim,node) - if (nodeTwin > 0_pInt) & ! if I am a twin of some node... - node_count(nodeTwin) = node_count(nodeTwin) + 1_pInt ! -> count me again for the twin node - enddo - endif - node_seen(j) = 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)) -mesh_sharedElem = 0_pInt - -do e = 1_pInt,mesh_NcpElems - t = mesh_element(2,e) - node_seen = 0_pInt - do j = 1_pInt,FE_Nnodes(t) - node = mesh_FEasCP('node',mesh_element(4_pInt+j,e)) - if (all(node_seen /= node)) then - mesh_sharedElem(1,node) = mesh_sharedElem(1,node) + 1_pInt ! count for each node the connected elements - mesh_sharedElem(mesh_sharedElem(1,node)+1_pInt,node) = e ! store the respective element id - do myDim = 1_pInt,3_pInt ! check in each dimension... - nodeTwin = mesh_nodeTwins(myDim,node) - if (nodeTwin > 0_pInt) then ! if i am a twin of some node... - mesh_sharedElem(1,nodeTwin) = mesh_sharedElem(1,nodeTwin) + 1_pInt ! ...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(j) = node - enddo -enddo - -deallocate(node_seen) - -end subroutine mesh_build_sharedElems - - -!*********************************************************** -! build up of IP neighborhood -! -! allocate globals -! _ipNeighborhood -!*********************************************************** -subroutine mesh_build_ipNeighborhood - -implicit none -integer(pInt) 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 -integer(pInt), dimension(FE_maxmaxNnodesAtIP) :: & - linkedNodes = 0_pInt, & - matchingNodes -logical checkTwins - -allocate(mesh_ipNeighborhood(2,mesh_maxNipNeighbors,mesh_maxNips,mesh_NcpElems)) -mesh_ipNeighborhood = 0_pInt - - -do myElem = 1_pInt,mesh_NcpElems ! loop over cpElems - myType = mesh_element(2,myElem) ! get elemType - do myIP = 1_pInt,FE_Nips(myType) ! loop over IPs of elem - - do neighbor = 1_pInt,FE_NipNeighbors(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_pInt) 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_pInt) 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_pInt) then ! found match? - neighboringType = mesh_element(2,matchingElem) - - !*** trivial solution if neighbor has only one IP - - if (FE_Nips(neighboringType) == 1_pInt) then - mesh_ipNeighborhood(1,neighbor,myIP,myElem) = matchingElem - mesh_ipNeighborhood(2,neighbor,myIP,myElem) = 1_pInt - cycle - endif - - !*** find those nodes which build the link to the neighbor - - NlinkedNodes = 0_pInt - linkedNodes = 0_pInt - do a = 1_pInt,FE_maxNnodesAtIP(myType) ! figure my anchor nodes on connecting face - anchor = FE_nodesAtIP(a,myIP,myType) - if (anchor /= 0_pInt) then ! valid anchor node - if (any(FE_nodeOnFace(:,myFace,myType) == anchor)) then ! ip anchor sits on face? - NlinkedNodes = NlinkedNodes + 1_pInt - linkedNodes(NlinkedNodes) = & - mesh_FEasCP('node',mesh_element(4_pInt+anchor,myElem)) ! CP id of anchor node - else ! something went wrong with the linkage, since not all anchors sit on my face - NlinkedNodes = 0_pInt - linkedNodes = 0_pInt - 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_pInt,FE_Nips(neighboringType) - NmatchingNodes = 0_pInt - matchingNodes = 0_pInt - do a = 1_pInt,FE_maxNnodesAtIP(neighboringType) ! check each anchor node of that ip - anchor = FE_nodesAtIP(a,candidateIP,neighboringType) - if (anchor /= 0_pInt) then ! valid anchor node - if (any(FE_nodeOnFace(:,matchingFace,neighboringType) == anchor)) then ! sits on matching face? - NmatchingNodes = NmatchingNodes + 1_pInt - matchingNodes(NmatchingNodes) = & - mesh_FEasCP('node',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_pInt - matchingNodes = 0_pInt - 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_pInt,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_pInt,NlinkedNodes - twin_of_linkedNode = mesh_nodeTwins(dir,linkedNodes(a)) - if (twin_of_linkedNode == 0_pInt .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 - -end subroutine mesh_build_ipNeighborhood - - - -!*********************************************************** -! assignment of coordinates for subnodes in each cp element -! -! allocate globals -! _subNodeCoord -!*********************************************************** -subroutine mesh_build_subNodeCoords - - implicit none - integer(pInt) e,t,n,p,Nparents - - if (.not. allocated(mesh_subNodeCoord)) then - allocate(mesh_subNodeCoord(3,mesh_maxNnodes+mesh_maxNsubNodes,mesh_NcpElems)) - endif - mesh_subNodeCoord = 0.0_pReal - - do e = 1_pInt,mesh_NcpElems ! loop over cpElems - t = mesh_element(2,e) ! get elemType - do n = 1_pInt,FE_Nnodes(t) - mesh_subNodeCoord(1:3,n,e) = mesh_node(1:3,mesh_FEasCP('node',mesh_element(4_pInt+n,e))) ! loop over nodes of this element type - enddo - do n = 1_pInt,FE_NsubNodes(t) ! now for the true subnodes - Nparents = count(FE_subNodeParent(1_pInt:FE_Nips(t),n,t) > 0_pInt) - do p = 1_pInt,Nparents ! loop through present parent nodes - mesh_subNodeCoord(1:3,FE_Nnodes(t)+n,e) = mesh_subNodeCoord(1:3,FE_Nnodes(t)+n,e) & - + mesh_node(1:3,mesh_FEasCP('node',mesh_element(4_pInt+FE_subNodeParent(p,n,t),e))) ! add up parents - enddo - mesh_subNodeCoord(1:3,n+FE_Nnodes(t),e) = mesh_subNodeCoord(1:3,n+FE_Nnodes(t),e)/real(Nparents,pReal) - enddo - enddo - -end subroutine mesh_build_subNodeCoords - - -!*********************************************************** -! calculation of IP coordinates -! -! allocate globals -! _ipCenterOfGravity -!*********************************************************** -subroutine mesh_build_ipCoordinates - - use prec, only: tol_gravityNodePos - - implicit none - integer(pInt) :: e,f,t,i,j,k,n - logical, dimension(mesh_maxNnodes+mesh_maxNsubNodes) :: gravityNode ! flagList to find subnodes determining center of grav - real(pReal), dimension(3,mesh_maxNnodes+mesh_maxNsubNodes) :: gravityNodePos ! coordinates of subnodes determining center of grav - real(pReal), dimension(3) :: centerOfGravity - - if (.not. allocated(mesh_ipCenterOfGravity)) allocate(mesh_ipCenterOfGravity(3,mesh_maxNips,mesh_NcpElems)) - - do e = 1_pInt,mesh_NcpElems ! loop over cpElems - t = mesh_element(2,e) ! get elemType - do i = 1_pInt,FE_Nips(t) ! loop over IPs of elem - gravityNode = .false. ! reset flagList - gravityNodePos = 0.0_pReal ! reset coordinates - do f = 1_pInt,FE_NipNeighbors(t) ! loop over interfaces of IP - do n = 1_pInt,FE_NipFaceNodes ! loop over nodes on interface - gravityNode(FE_subNodeOnIPFace(n,f,i,t)) = .true. - gravityNodePos(:,FE_subNodeOnIPFace(n,f,i,t)) = mesh_subNodeCoord(:,FE_subNodeOnIPFace(n,f,i,t),e) - enddo - enddo - - do j = 1_pInt,mesh_maxNnodes+mesh_maxNsubNodes-1_pInt ! walk through entire flagList except last - if (gravityNode(j)) then ! valid node index - do k = j+1_pInt,mesh_maxNnodes+mesh_maxNsubNodes ! walk through remainder of list - if (gravityNode(k) .and. all(abs(gravityNodePos(:,j) - gravityNodePos(:,k)) < tol_gravityNodePos)) then ! found duplicate - gravityNode(j) = .false. ! delete first instance - gravityNodePos(:,j) = 0.0_pReal - exit ! continue with next suspect - endif - enddo - endif - enddo - centerOfGravity = sum(gravityNodePos,2)/real(count(gravityNode),pReal) - mesh_ipCenterOfGravity(:,i,e) = centerOfGravity - enddo - enddo - -end subroutine mesh_build_ipCoordinates - - -!*********************************************************** -! calculation of IP volume -! -! allocate globals -! _ipVolume -!*********************************************************** -subroutine mesh_build_ipVolumes - - use math, only: math_volTetrahedron - implicit none - - integer(pInt) :: e,f,t,i,j,n - integer(pInt), parameter :: Ntriangles = FE_NipFaceNodes-2_pInt ! each interface is made up of this many triangles - real(pReal), dimension(3,FE_NipFaceNodes) :: nPos ! coordinates of nodes on IP face - real(pReal), dimension(Ntriangles,FE_NipFaceNodes) :: volume ! volumes of possible tetrahedra - - if (.not. allocated(mesh_ipVolume)) then - allocate(mesh_ipVolume(mesh_maxNips,mesh_NcpElems)) - endif - - mesh_ipVolume = 0.0_pReal - do e = 1_pInt,mesh_NcpElems ! loop over cpElems - t = mesh_element(2_pInt,e) ! get elemType - do i = 1_pInt,FE_Nips(t) ! loop over IPs of elem - do f = 1_pInt,FE_NipNeighbors(t) ! loop over interfaces of IP and add tetrahedra which connect to CoG - forall (n = 1_pInt:FE_NipFaceNodes) & - nPos(:,n) = mesh_subNodeCoord(:,FE_subNodeOnIPFace(n,f,i,t),e) - forall (n = 1_pInt:FE_NipFaceNodes, j = 1_pInt:Ntriangles) & ! start at each interface node and build valid triangles to cover interface - volume(j,n) = math_volTetrahedron(nPos(:,n), & ! calc volume of respective tetrahedron to CoG - nPos(:,1_pInt+mod(n-1_pInt +j ,FE_NipFaceNodes)), & ! start at offset j - nPos(:,1_pInt+mod(n-1_pInt +j+1_pInt,FE_NipFaceNodes)), & ! and take j's neighbor - mesh_ipCenterOfGravity(:,i,e)) - mesh_ipVolume(i,e) = mesh_ipVolume(i,e) + sum(volume) ! add contribution from this interface - enddo - mesh_ipVolume(i,e) = mesh_ipVolume(i,e) / FE_NipFaceNodes ! renormalize with interfaceNodeNum due to loop over them - enddo - enddo - -end subroutine mesh_build_ipVolumes - - -!*********************************************************** -! calculation of IP interface areas -! -! allocate globals -! _ipArea, _ipAreaNormal -!*********************************************************** -subroutine mesh_build_ipAreas - - use math, only: math_vectorproduct - - implicit none - integer(pInt) :: e,f,t,i,j,n - integer(pInt), parameter :: Ntriangles = FE_NipFaceNodes-2_pInt ! each interface is made up of this many triangles - real(pReal), dimension (3,FE_NipFaceNodes) :: nPos ! coordinates of nodes on IP face - real(pReal), dimension(3,Ntriangles,FE_NipFaceNodes) :: normal - real(pReal), dimension(Ntriangles,FE_NipFaceNodes) :: area - - allocate(mesh_ipArea(mesh_maxNipNeighbors,mesh_maxNips,mesh_NcpElems)) ; mesh_ipArea = 0.0_pReal - allocate(mesh_ipAreaNormal(3_pInt,mesh_maxNipNeighbors,mesh_maxNips,mesh_NcpElems)) ; mesh_ipAreaNormal = 0.0_pReal - do e = 1_pInt,mesh_NcpElems ! loop over cpElems - t = mesh_element(2,e) ! get elemType - do i = 1_pInt,FE_Nips(t) ! loop over IPs of elem - do f = 1_pInt,FE_NipNeighbors(t) ! loop over interfaces of IP - forall (n = 1_pInt:FE_NipFaceNodes) nPos(:,n) = mesh_subNodeCoord(:,FE_subNodeOnIPFace(n,f,i,t),e) - forall (n = 1_pInt:FE_NipFaceNodes, j = 1_pInt:Ntriangles) ! start at each interface node and build valid triangles to cover interface - normal(:,j,n) = math_vectorproduct(nPos(:,1_pInt+mod(n+j-1_pInt,FE_NipFaceNodes)) - nPos(:,n), & ! calc their normal vectors - nPos(:,1_pInt+mod(n+j-0_pInt,FE_NipFaceNodes)) - nPos(:,n)) - area(j,n) = sqrt(sum(normal(:,j,n)*normal(:,j,n))) ! and area - end forall - forall (n = 1_pInt:FE_NipFaceNodes, j = 1_pInt:Ntriangles, area(j,n) > 0.0_pReal) & - normal(1:3,j,n) = normal(1:3,j,n) / area(j,n) ! make myUnit normal - - mesh_ipArea(f,i,e) = sum(area) / (FE_NipFaceNodes*2.0_pReal) ! area of parallelograms instead of triangles - mesh_ipAreaNormal(:,f,i,e) = sum(sum(normal,3),2_pInt)/& ! average of all valid normals - real(count(area > 0.0_pReal),pReal) - enddo - enddo - enddo - - end subroutine mesh_build_ipAreas - - -!*********************************************************** -! assignment of twin nodes for each cp node -! -! allocate globals -! _nodeTwins -!*********************************************************** -subroutine mesh_build_nodeTwins - -implicit none -integer(pInt) dir, & ! direction of periodicity - node, & - minimumNode, & - maximumNode, & - n1, & - n2 -integer(pInt), 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_pInt - -tolerance = 0.001_pReal * minval(mesh_ipVolume) ** 0.333_pReal - -do dir = 1_pInt,3_pInt ! 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_pInt - maximumNodes = 0_pInt - minCoord = minval(mesh_node0(dir,:)) - maxCoord = maxval(mesh_node0(dir,:)) - do node = 1_pInt,mesh_Nnodes ! loop through all nodes and find surface nodes - if (abs(mesh_node0(dir,node) - minCoord) <= tolerance) then - minimumNodes(1) = minimumNodes(1) + 1_pInt - minimumNodes(minimumNodes(1)+1_pInt) = node - elseif (abs(mesh_node0(dir,node) - maxCoord) <= tolerance) then - maximumNodes(1) = maximumNodes(1) + 1_pInt - maximumNodes(maximumNodes(1)+1_pInt) = node - endif - enddo - - - !*** find the corresponding node on the other side with the same position in this dimension - - unpaired = .true. - do n1 = 1_pInt,minimumNodes(1) - minimumNode = minimumNodes(n1+1_pInt) - if (unpaired(minimumNode)) then - do n2 = 1_pInt,maximumNodes(1) - maximumNode = maximumNodes(n2+1_pInt) - 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 - - - -!*********************************************************** -! write statistics regarding input file parsing -! to the output file -! -!*********************************************************** -subroutine mesh_tell_statistics - - use math, only: math_range - use IO, only: IO_error - use debug, only: debug_what, & - debug_mesh, & - debug_levelBasic, & - debug_levelExtensive, & - debug_levelSelective, & - debug_e, & - debug_i - - implicit none - integer(pInt), dimension (:,:), allocatable :: mesh_HomogMicro - character(len=64) :: myFmt - integer(pInt) :: i,e,n,f,t, myDebug - - myDebug = debug_what(debug_mesh) - - if (mesh_maxValStateVar(1) < 1_pInt) call IO_error(error_ID=170_pInt) ! no homogenization specified - if (mesh_maxValStateVar(2) < 1_pInt) call IO_error(error_ID=180_pInt) ! no microstructure specified - - allocate (mesh_HomogMicro(mesh_maxValStateVar(1),mesh_maxValStateVar(2))); mesh_HomogMicro = 0_pInt -do e = 1_pInt,mesh_NcpElems - if (mesh_element(3,e) < 1_pInt) call IO_error(error_ID=170_pInt,e=e) ! no homogenization specified - if (mesh_element(4,e) < 1_pInt) call IO_error(error_ID=180_pInt,e=e) ! no microstructure specified - mesh_HomogMicro(mesh_element(3,e),mesh_element(4,e)) = & - mesh_HomogMicro(mesh_element(3,e),mesh_element(4,e)) + 1_pInt ! count combinations of homogenization and microstructure -enddo -!$OMP CRITICAL (write2out) - if (iand(myDebug,debug_levelBasic) /= 0_pInt) then - write (6,*) - write (6,*) 'Input Parser: STATISTICS' - write (6,*) - write (6,*) mesh_Nelems, ' : total number of elements in mesh' - write (6,*) mesh_NcpElems, ' : total number of CP elements in mesh' - write (6,*) mesh_Nnodes, ' : total number of nodes in mesh' - write (6,*) mesh_maxNnodes, ' : max number of nodes in any CP element' - write (6,*) mesh_maxNips, ' : max number of IPs in any CP element' - write (6,*) mesh_maxNipNeighbors, ' : max number of IP neighbors in any CP element' - write (6,*) mesh_maxNsubNodes, ' : max number of (additional) subnodes in any CP element' - write (6,*) mesh_maxNsharedElems, ' : max number of CP elements sharing a node' - write (6,*) - write (6,*) 'Input Parser: HOMOGENIZATION/MICROSTRUCTURE' - write (6,*) - write (6,*) mesh_maxValStateVar(1), ' : maximum homogenization index' - write (6,*) mesh_maxValStateVar(2), ' : maximum microstructure index' - write (6,*) - write (myFmt,'(a,i32.32,a)') '(9x,a2,1x,',mesh_maxValStateVar(2),'(i8))' - write (6,myFmt) '+-',math_range(mesh_maxValStateVar(2)) - write (myFmt,'(a,i32.32,a)') '(i8,1x,a2,1x,',mesh_maxValStateVar(2),'(i8))' - do i=1_pInt,mesh_maxValStateVar(1) ! loop over all (possibly assigned) homogenizations - write (6,myFmt) i,'| ',mesh_HomogMicro(i,:) ! loop over all (possibly assigned) microstructures - enddo - write(6,*) - write(6,*) 'Input Parser: ADDITIONAL MPIE OPTIONS' - write(6,*) - write(6,*) 'periodic surface : ', mesh_periodicSurface - write(6,*) - call flush(6) - endif - - if (iand(myDebug,debug_levelExtensive) /= 0_pInt) then - write (6,*) - write (6,*) 'Input Parser: SUBNODE COORDINATES' - write (6,*) - write(6,'(a8,1x,a5,1x,2(a15,1x),a20,3(1x,a12))')& - 'elem','IP','IP neighbor','IPFaceNodes','subNodeOnIPFace','x','y','z' - do e = 1_pInt,mesh_NcpElems ! loop over cpElems - if (iand(myDebug,debug_levelSelective) /= 0_pInt .and. debug_e /= e) cycle - t = mesh_element(2,e) ! get elemType - do i = 1_pInt,FE_Nips(t) ! loop over IPs of elem - if (iand(myDebug,debug_levelSelective) /= 0_pInt .and. debug_i /= i) cycle - do f = 1_pInt,FE_NipNeighbors(t) ! loop over interfaces of IP - do n = 1_pInt,FE_NipFaceNodes ! loop over nodes on interface - write(6,'(i8,1x,i5,2(1x,i15),1x,i20,3(1x,f12.8))') e,i,f,n,FE_subNodeOnIPFace(n,f,i,t),& - mesh_subNodeCoord(1,FE_subNodeOnIPFace(n,f,i,t),e),& - mesh_subNodeCoord(2,FE_subNodeOnIPFace(n,f,i,t),e),& - mesh_subNodeCoord(3,FE_subNodeOnIPFace(n,f,i,t),e) - enddo - enddo - enddo - enddo - write(6,*) - write(6,*) 'Input Parser: IP COORDINATES' - write(6,'(a8,1x,a5,3(1x,a12))') 'elem','IP','x','y','z' - do e = 1_pInt,mesh_NcpElems - if (iand(myDebug,debug_levelSelective) /= 0_pInt .and. debug_e /= e) cycle - do i = 1_pInt,FE_Nips(mesh_element(2,e)) - if (iand(myDebug,debug_levelSelective) /= 0_pInt .and. debug_i /= i) cycle - write (6,'(i8,1x,i5,3(1x,f12.8))') e, i, mesh_ipCenterOfGravity(:,i,e) - enddo - enddo - write (6,*) - write (6,*) 'Input Parser: ELEMENT VOLUME' - write (6,*) - write (6,'(a13,1x,e15.8)') 'total volume', sum(mesh_ipVolume) - write (6,*) - write (6,'(a8,1x,a5,1x,a15,1x,a5,1x,a15,1x,a16)') 'elem','IP','volume','face','area','-- normal --' - do e = 1_pInt,mesh_NcpElems - if (iand(myDebug,debug_levelSelective) /= 0_pInt .and. debug_e /= e) cycle - do i = 1_pInt,FE_Nips(mesh_element(2,e)) - if (iand(myDebug,debug_levelSelective) /= 0_pInt .and. debug_i /= i) cycle - write (6,'(i8,1x,i5,1x,e15.8)') e,i,mesh_IPvolume(i,e) - do f = 1_pInt,FE_NipNeighbors(mesh_element(2,e)) - write (6,'(i33,1x,e15.8,1x,3(f6.3,1x))') f,mesh_ipArea(f,i,e),mesh_ipAreaNormal(:,f,i,e) - enddo - enddo - enddo - write (6,*) - write (6,*) 'Input Parser: NODE TWINS' - write (6,*) - write(6,'(a6,3(3x,a6))') ' node','twin_x','twin_y','twin_z' - do n = 1_pInt,mesh_Nnodes ! loop over cpNodes - if (debug_e <= mesh_NcpElems) then - if (any(mesh_element(5:,debug_e) == n)) then - write(6,'(i6,3(3x,i6))') n, mesh_nodeTwins(1:3,n) - endif - endif - enddo - write(6,*) - write(6,*) 'Input Parser: IP NEIGHBORHOOD' - write(6,*) - write(6,'(a8,1x,a10,1x,a10,1x,a3,1x,a13,1x,a13)') 'elem','IP','neighbor','','elemNeighbor','ipNeighbor' - do e = 1_pInt,mesh_NcpElems ! loop over cpElems - if (iand(myDebug,debug_levelSelective) /= 0_pInt .and. debug_e /= e) cycle - t = mesh_element(2,e) ! get elemType - do i = 1_pInt,FE_Nips(t) ! loop over IPs of elem - if (iand(myDebug,debug_levelSelective) /= 0_pInt .and. debug_i /= i) cycle - do n = 1_pInt,FE_NipNeighbors(t) ! loop over neighbors of IP - write (6,'(i8,1x,i10,1x,i10,1x,a3,1x,i13,1x,i13)') e,i,n,'-->',mesh_ipNeighborhood(1,n,i,e),mesh_ipNeighborhood(2,n,i,e) - enddo - enddo - enddo - endif -!$OMP END CRITICAL (write2out) - -deallocate(mesh_HomogMicro) - -end subroutine mesh_tell_statistics - - - -function mesh_regrid(resNewInput) !use new_res=[0.0,0.0,0.0] for automatic determination of new grid - use prec, only: & - pInt, & - pReal - use DAMASK_interface, only: & - getSolverJobName - use IO, only: & - IO_read_jobBinaryFile ,& - IO_write_jobBinaryFile - use math, only: & - math_nearestNeighborSearch, & - deformed_FFT, & - math_mul33x3 - - integer(pInt), dimension(3) :: mesh_regrid - integer(pInt), dimension(3), optional, intent(in) :: resNewInput - integer(pInt):: maxsize, i, j, k, ielem, Npoints, NpointsNew, spatialDim - integer(pInt), dimension(3) :: res, resNew - integer(pInt), dimension(:), allocatable :: indices - real(pReal), dimension(3) :: geomdim - real(pReal), dimension(3,3) :: Favg - real(pReal), dimension(:,:), allocatable :: & - coordinatesNew, & - coordinatesLinear - - real(pReal), dimension(:,:,:), allocatable :: & - material_phase, material_phaseNew - - real(pReal), dimension(:,:,:,:,:), allocatable :: & - F, FNew, & - Fp, FpNew, & - Lp, LpNew, & - dcsdE, dcsdENew - - real(pReal), dimension (:,:,:,:,:,:,:), allocatable :: & - dPdF, dPdFNew - real(pReal), dimension (:,:,:,:), allocatable :: & - coordinates, & - Tstar, TstarNew, & - stateHomog, & - stateConst - integer(pInt), dimension(:,:), allocatable :: & - sizeStateConst, sizeStateHomog - - res = mesh_spectral_getResolution() - geomdim = mesh_spectral_getDimension() - Npoints = res(1)*res(2)*res(3) - - - - allocate(F(res(1),res(2),res(3),3,3)) - call IO_read_jobBinaryFile(777,'convergedSpectralDefgrad',trim(getSolverJobName()),size(F)) - read (777,rec=1) F - close (777) - -! ----Calculate deformed configuration and average-------- - do i= 1_pInt,3_pInt; do j = 1_pInt,3_pInt - Favg(i,j) = sum(F(1:res(1),1:res(2),1:res(3),i,j)) / Npoints - enddo; enddo - allocate(coordinates(res(1),res(2),res(3),3)) - call deformed_fft(res,geomdim,Favg,1.0_pReal,F,coordinates) - deallocate(F) - -! ----Store coordinates into a linear list-------- - allocate(coordinatesLinear(3,Npoints)) - ielem = 0_pInt - do k=1_pInt,res(3); do j=1_pInt, res(2); do i=1_pInt, res(1) - ielem = ielem + 1_pInt - coordinatesLinear(1:3,ielem) = coordinates(i,j,k,1:3) - enddo; enddo; enddo - deallocate(coordinates) - - ! ----For 2D /3D case---------------------------------- - if (res(3)== 1_pInt) then - spatialDim = 2_pInt - else - spatialDim = 3_pInt - endif - -! ----determine/calculate new resolution-------------- - if (.not. present(resNewInput)) then - resNew = res - else - if(all(resNewInput==0.0_pReal)) then - mesh_regrid =[ 0,0,0] - return !no automatic determination right now - else - resNew = resNewInput - endif - endif - NpointsNew = resNew(1)*resNew(2)*resNew(3) - -! ----Calculate regular new coordinates------------ - allocate(coordinatesNew(3,NpointsNew)) - ielem = 0_pInt - do k=1_pInt,resNew(3); do j=1_pInt, resNew(2); do i=1_pInt, resNew(1) - ielem = ielem + 1_pInt - coordinatesNew(1:3,ielem) = math_mul33x3(Favg, geomdim/real(resNew,pReal)*real([i,j,k],pReal) & - - geomdim/real(2_pInt*resNew,pReal)) - enddo; enddo; enddo - -!----- Nearest neighbour search ----------------------- - allocate(indices(Npoints)) - call math_nearestNeighborSearch(spatialDim, Favg, geomdim, NpointsNew, Npoints, & - coordinatesNew, coordinatesLinear, indices) - deallocate(coordinatesNew) - indices = indices / 3_pInt**spatialDim +1_pInt - -!--------------------------------------------------------------------------- - allocate(material_phase (1,1, Npoints)) - allocate(material_phaseNew (1,1, NpointsNew)) - call IO_read_jobBinaryFile(777,'recordedPhase',trim(getSolverJobName()),size(material_phase)) - read (777,rec=1) material_phase - close (777) - do i = 1, NpointsNew - material_phaseNew(1,1,i) = material_phase(1,1,indices(i)) - enddo - - call IO_write_jobBinaryFile(777,'recordedPhase',size(material_phaseNew)) - write (777,rec=1) material_phaseNew - close (777) - deallocate(material_phase) - deallocate(material_phaseNew) - -!--------------------------------------------------------------------------- - allocate(F (3,3,1,1, Npoints)) - allocate(FNew (3,3,1,1, NpointsNew)) - call IO_read_jobBinaryFile(777,'convergedF',trim(getSolverJobName()),size(F)) - read (777,rec=1) F - close (777) - do i = 1, NpointsNew - FNew(1:3,1:3,1,1,i) = F(1:3,1:3,1,1,indices(i)) - enddo - - call IO_write_jobBinaryFile(777,'convergedF',size(FNew)) - write (777,rec=1) FNew - close (777) - deallocate(F) - deallocate(FNew) - -!--------------------------------------------------------------------- - allocate(Fp (3,3,1,1,Npoints)) - allocate(FpNew (3,3,1,1,NpointsNew)) - call IO_read_jobBinaryFile(777,'convergedFp',trim(getSolverJobName()),size(Fp)) - read (777,rec=1) Fp - close (777) - do i = 1, NpointsNew - FpNew(1:3,1:3,1,1,i) = Fp(1:3,1:3,1,1,indices(i)) - enddo - - call IO_write_jobBinaryFile(777,'convergedFp',size(FpNew)) - write (777,rec=1) FpNew - close (777) - deallocate(Fp) - deallocate(FpNew) - -!------------------------------------------------------------------------ - allocate(Lp (3,3,1,1,Npoints)) - allocate(LpNew (3,3,1,1,NpointsNew)) - call IO_read_jobBinaryFile(777,'convergedLp',trim(getSolverJobName()),size(Lp)) - read (777,rec=1) Lp - close (777) - do i = 1, NpointsNew - LpNew(1:3,1:3,1,1,i) = Lp(1:3,1:3,1,1,indices(i)) - enddo - call IO_write_jobBinaryFile(777,'convergedLp',size(LpNew)) - write (777,rec=1) LpNew - close (777) - deallocate(Lp) - deallocate(LpNew) - -!---------------------------------------------------------------------------- - allocate(dcsdE (6,6,1,1,Npoints)) - allocate(dcsdENew (6,6,1,1,NpointsNew)) - call IO_read_jobBinaryFile(777,'convergeddcsdE',trim(getSolverJobName()),size(dcsdE)) - read (777,rec=1) dcsdE - close (777) - do i = 1, NpointsNew - dcsdENew(1:6,1:6,1,1,i) = dcsdE(1:6,1:6,1,1,indices(i)) - enddo - call IO_write_jobBinaryFile(777,'convergeddcsdE',size(dcsdENew)) - write (777,rec=1) dcsdENew - close (777) - deallocate(dcsdE) - deallocate(dcsdENew) - -!--------------------------------------------------------------------------- - allocate(dPdF (3,3,3,3,1,1,Npoints)) - allocate(dPdFNew (3,3,3,3,1,1,NpointsNew)) - call IO_read_jobBinaryFile(777,'convergeddPdF',trim(getSolverJobName()),size(dPdF)) - read (777,rec=1) dPdF - close (777) - do i = 1, NpointsNew - dPdFNew(1:3,1:3,1:3,1:3,1,1,i) = dPdF(1:3,1:3,1:3,1:3,1,1,indices(i)) - enddo - call IO_write_jobBinaryFile(777,'convergeddPdF',size(dPdFNew)) - write (777,rec=1) dPdFNew - close (777) - deallocate(dPdF) - deallocate(dPdFNew) - -!--------------------------------------------------------------------------- - allocate(Tstar (6,1,1,Npoints)) - allocate(TstarNew (6,1,1,NpointsNew)) - call IO_read_jobBinaryFile(777,'convergedTstar',trim(getSolverJobName()),size(Tstar)) - read (777,rec=1) Tstar - close (777) - do i = 1, NpointsNew - TstarNew(1:6,1,1,i) = Tstar(1:6,1,1,indices(i)) - enddo - call IO_write_jobBinaryFile(777,'convergedTstar',size(TstarNew)) - write (777,rec=1) TstarNew - close (777) - deallocate(Tstar) - deallocate(TstarNew) - -!---------------------------------------------------------------------------- - allocate(sizeStateConst(1,Npoints)) - call IO_read_jobBinaryFile(777,'sizeStateConst',trim(getSolverJobName()),size(sizeStateConst)) - read (777,rec=1) sizeStateConst - close (777) - maxsize = maxval(sizeStateConst) - allocate(StateConst (1,1,Npoints,maxsize)) - - call IO_read_jobBinaryFile(777,'convergedStateConst',trim(getSolverJobName())) - k = 0_pInt - do i =1, Npoints - do j = 1,sizeStateConst(1,i) - k = k+1_pInt - read(777,rec=k) StateConst(1,1,i,j) - enddo - enddo - close(777) - - call IO_write_jobBinaryFile(777,'convergedStateConst') - k = 0_pInt - do i = 1,NpointsNew - do j = 1,sizeStateConst(1,i) - k=k+1_pInt - write(777,rec=k) StateConst(1,1,indices(i),j) - enddo - enddo - close (777) - deallocate(sizeStateConst) - deallocate(StateConst) - -!---------------------------------------------------------------------------- - allocate(sizeStateHomog(1,Npoints)) - call IO_read_jobBinaryFile(777,'sizeStateHomog',trim(getSolverJobName()),size(sizeStateHomog)) - read (777,rec=1) sizeStateHomog - close (777) - maxsize = maxval(sizeStateHomog) - allocate(stateHomog (1,1,Npoints,maxsize)) - - call IO_read_jobBinaryFile(777,'convergedStateHomog',trim(getSolverJobName())) - k = 0_pInt - do i =1, Npoints - do j = 1,sizeStateHomog(1,i) - k = k+1_pInt - read(777,rec=k) stateHomog(1,1,i,j) - enddo - enddo - close(777) - - call IO_write_jobBinaryFile(777,'convergedStateHomog') - k = 0_pInt - do i = 1,NpointsNew - do j = 1,sizeStateHomog(1,i) - k=k+1_pInt - write(777,rec=k) stateHomog(1,1,indices(i),j) - enddo - enddo - close (777) - deallocate(sizeStateHomog) - deallocate(stateHomog) - - deallocate(indices) - - mesh_regrid = resNew -end function mesh_regrid - - -function mesh_spectral_getDimension(fileUnit) - use IO, only: & - IO_checkAndRewind, & - IO_open_file, & - IO_stringPos, & - IO_lc, & - IO_stringValue, & - IO_intValue, & - IO_floatValue, & - IO_error - use DAMASK_interface, only: & - getModelName, & - InputFileExtension - - implicit none - integer(pInt), dimension(1_pInt + 7_pInt*2_pInt) :: positions ! for a,b c + 3 values + keyword - integer(pInt), intent(in), optional :: fileUnit - integer(pInt) :: headerLength = 0_pInt - real(pReal), dimension(3) :: mesh_spectral_getDimension - character(len=1024) :: line, & - keyword - integer(pInt) :: i, j - logical :: gotDimension = .false. - integer(pInt) :: myUnit - - if(.not. present(fileUnit)) then - myUnit = 289_pInt - call IO_open_file(myUnit,trim(getModelName())//InputFileExtension) - else - myUnit = fileUnit - endif - - call IO_checkAndRewind(myUnit) - - read(myUnit,'(a1024)') line - positions = IO_stringPos(line,2_pInt) - keyword = IO_lc(IO_StringValue(line,positions,2_pInt)) - if (keyword(1:4) == 'head') then - headerLength = IO_intValue(line,positions,1_pInt) + 1_pInt - else - call IO_error(error_ID=841_pInt, ext_msg='mesh_spectral_getDimension') - endif - rewind(myUnit) - do i = 1_pInt, headerLength - read(myUnit,'(a1024)') line - positions = IO_stringPos(line,7_pInt) - select case ( IO_lc(IO_StringValue(line,positions,1)) ) - case ('dimension') - gotDimension = .true. - do j = 2_pInt,6_pInt,2_pInt - select case (IO_lc(IO_stringValue(line,positions,j))) - case('x') - mesh_spectral_getDimension(1) = IO_floatValue(line,positions,j+1_pInt) - case('y') - mesh_spectral_getDimension(2) = IO_floatValue(line,positions,j+1_pInt) - case('z') - mesh_spectral_getDimension(3) = IO_floatValue(line,positions,j+1_pInt) - end select - enddo - end select - enddo - - if(.not. present(fileUnit)) close(myUnit) - - if (.not. gotDimension) & - call IO_error(error_ID = 845_pInt, ext_msg='dimension') - if (any(mesh_spectral_getDimension<=0.0_pReal)) & - call IO_error(error_ID = 844_pInt, ext_msg='mesh_spectral_getDimension') - -end function mesh_spectral_getDimension - - -function mesh_spectral_getResolution(fileUnit) - use IO, only: & - IO_checkAndRewind, & - IO_open_file, & - IO_stringPos, & - IO_lc, & - IO_stringValue, & - IO_intValue, & - IO_floatValue, & - IO_error - use DAMASK_interface, only: & - getModelName, & - InputFileExtension - - implicit none - integer(pInt), dimension(1_pInt + 7_pInt*2_pInt) :: positions ! for a,b c + 3 values + keyword - integer(pInt), intent(in), optional :: fileUnit - integer(pInt) :: headerLength = 0_pInt - integer(pInt), dimension(3) :: mesh_spectral_getResolution - character(len=1024) :: line, & - keyword - integer(pInt) :: i, j - logical :: gotResolution = .false. - integer(pInt) :: myUnit - - if(.not. present(fileUnit)) then - myUnit = 289_pInt - call IO_open_file(myUnit,trim(getModelName())//InputFileExtension) - else - myUnit = fileUnit - endif - - call IO_checkAndRewind(myUnit) - - read(myUnit,'(a1024)') line - positions = IO_stringPos(line,2_pInt) - keyword = IO_lc(IO_StringValue(line,positions,2_pInt)) - if (keyword(1:4) == 'head') then - headerLength = IO_intValue(line,positions,1_pInt) + 1_pInt - else - call IO_error(error_ID=841_pInt, ext_msg='mesh_spectral_getResolution') - endif - rewind(myUnit) - do i = 1_pInt, headerLength - read(myUnit,'(a1024)') line - positions = IO_stringPos(line,7_pInt) - select case ( IO_lc(IO_StringValue(line,positions,1_pInt)) ) - case ('resolution') - gotResolution = .true. - do j = 2_pInt,6_pInt,2_pInt - select case (IO_lc(IO_stringValue(line,positions,j))) - case('a') - mesh_spectral_getResolution(1) = IO_intValue(line,positions,j+1_pInt) - case('b') - mesh_spectral_getResolution(2) = IO_intValue(line,positions,j+1_pInt) - case('c') - mesh_spectral_getResolution(3) = IO_intValue(line,positions,j+1_pInt) - end select - enddo - end select - enddo - - if(.not. present(fileUnit)) close(myUnit) - - if (.not. gotResolution) & - call IO_error(error_ID = 845_pInt, ext_msg='resolution') - if((mod(mesh_spectral_getResolution(1),2_pInt)/=0_pInt .or. & ! must be a even number - mesh_spectral_getResolution(1) < 2_pInt .or. & ! and larger than 1 - mod(mesh_spectral_getResolution(2),2_pInt)/=0_pInt .or. & ! -"- - mesh_spectral_getResolution(2) < 2_pInt .or. & ! -"- - (mod(mesh_spectral_getResolution(3),2_pInt)/=0_pInt .and. & - mesh_spectral_getResolution(3)/= 1_pInt)) .or. & ! third res might be 1 - mesh_spectral_getResolution(3) < 1_pInt) & - call IO_error(error_ID = 843_pInt, ext_msg='mesh_spectral_getResolution') - -end function mesh_spectral_getResolution - -function mesh_spectral_getHomogenization(fileUnit) - use IO, only: & - IO_checkAndRewind, & - IO_open_file, & - IO_stringPos, & - IO_lc, & - IO_stringValue, & - IO_intValue, & - IO_error - use DAMASK_interface, only: & - getModelName, & - InputFileExtension - - implicit none - integer(pInt), dimension(1_pInt + 7_pInt*2_pInt) :: positions ! for a, b, c + 3 values + keyword - integer(pInt), intent(in), optional :: fileUnit - integer(pInt) :: headerLength = 0_pInt - integer(pInt) :: mesh_spectral_getHomogenization - character(len=1024) :: line, & - keyword - integer(pInt) :: i, j - logical :: gotHomogenization = .false. - integer(pInt) :: myUnit - - if(.not. present(fileUnit)) then - myUnit = 289_pInt - call IO_open_file(myUnit,trim(getModelName())//InputFileExtension) - else - myUnit = fileUnit - endif - - call IO_checkAndRewind(myUnit) - - read(myUnit,'(a1024)') line - positions = IO_stringPos(line,2_pInt) - keyword = IO_lc(IO_StringValue(line,positions,2_pInt)) - if (keyword(1:4) == 'head') then - headerLength = IO_intValue(line,positions,1_pInt) + 1_pInt - else - call IO_error(error_ID=841_pInt, ext_msg='mesh_spectral_getHomogenization') - endif - rewind(myUnit) - do i = 1_pInt, headerLength - read(myUnit,'(a1024)') line - positions = IO_stringPos(line,7_pInt) - select case ( IO_lc(IO_StringValue(line,positions,1)) ) - case ('homogenization') - gotHomogenization = .true. - mesh_spectral_getHomogenization = IO_intValue(line,positions,2_pInt) - end select - enddo - - if(.not. present(fileUnit)) close(myUnit) - - if (.not. gotHomogenization ) & - call IO_error(error_ID = 845_pInt, ext_msg='homogenization') - if (mesh_spectral_getHomogenization<1_pInt) & - call IO_error(error_ID = 842_pInt, ext_msg='mesh_spectral_getHomogenization') - -end function mesh_spectral_getHomogenization - end module mesh diff --git a/code/numerics.f90 b/code/numerics.f90 index dfa94a217..444f60d7f 100644 --- a/code/numerics.f90 +++ b/code/numerics.f90 @@ -64,9 +64,17 @@ real(pReal) :: relevantStrain = 1.0e-7_pReal, & maxdRelax_RGC = 1.0e+0_pReal, & ! threshold of maximum relaxation vector increment (if exceed this then cutback) maxVolDiscr_RGC = 1.0e-5_pReal, & ! threshold of maximum volume discrepancy allowed volDiscrMod_RGC = 1.0e+12_pReal, & ! stiffness of RGC volume discrepancy (zero = without volume discrepancy constraint) - volDiscrPow_RGC = 5.0_pReal, & ! powerlaw penalty for volume discrepancy + volDiscrPow_RGC = 5.0_pReal ! powerlaw penalty for volume discrepancy +logical :: analyticJaco = .false. ! use analytic Jacobian or perturbation, Default .false.: calculate Jacobian using perturbations +!* Random seeding parameters +integer(pInt) :: fixedSeed = 0_pInt ! fixed seeding for pseudo-random number generator, Default 0: use random seed +!* OpenMP variable +integer(pInt) :: DAMASK_NumThreadsInt = 0_pInt ! value stored in environment variable DAMASK_NUM_THREADS, set to zero if no OpenMP directive + + !* spectral parameters: - err_div_tol = 0.1_pReal, & ! Div(P)/avg(P)*meter +#ifdef Spectral +real(pReal) :: err_div_tol = 0.1_pReal, & ! Div(P)/avg(P)*meter err_stress_tolrel = 0.01_pReal, & ! relative tolerance for fullfillment of stress BC, Default: 0.01 allowing deviation of 1% of maximum stress err_stress_tolabs = huge(1.0_pReal), & ! absolute tolerance for fullfillment of stress BC, Default: 0.01 allowing deviation of 1% of maximum stress fftw_timelimit = -1.0_pReal, & ! sets the timelimit of plan creation for FFTW, see manual on www.fftw.org, Default -1.0: disable timelimit @@ -77,18 +85,11 @@ integer(pInt) :: fftw_planner_flag = 32_pInt, & itmin = 2_pInt ! minimum number of iterations logical :: memory_efficient = .true., & ! for fast execution (pre calculation of gamma_hat), Default .true.: do not precalculate divergence_correction = .false., & ! correct divergence calculation in fourier space, Default .false.: no correction - update_gamma = .false., & ! update gamma operator with current stiffness, Default .false.: use initial stiffness -!* end of spectral parameters: - analyticJaco = .false. ! use analytic Jacobian or perturbation, Default .false.: calculate Jacobian using perturbations + update_gamma = .false. ! update gamma operator with current stiffness, Default .false.: use initial stiffness +#endif -!* Random seeding parameters -integer(pInt) :: fixedSeed = 0_pInt ! fixed seeding for pseudo-random number generator, Default 0: use random seed -!* OpenMP variable -integer(pInt) :: DAMASK_NumThreadsInt = 0_pInt ! value stored in environment variable DAMASK_NUM_THREADS, set to zero if no OpenMP directive - - CONTAINS !******************************************* @@ -106,12 +107,12 @@ subroutine numerics_init IO_floatValue, & IO_intValue, & IO_warning -#ifdef STANDARD_CHECK ! If STANDARD_CHECK is defined (as in the makefile for the spectral solver by setting -!$ use OMP_LIB ! -DSTANDARD_CHECK use the module file for the openMP function library -#endif ! REASON: module file crashes with Marc but omp_lib.h is not standard conform - implicit none ! and ifort will does not compile it (gfortran seems to have an exeption) -#ifndef STANDARD_CHECK ! if STANDARD_CHECK is not defined (e.g. when compiling with Marc or Abaqus) -!$ include "omp_lib.h" ! use this file for the openMP function library +#ifndef Marc ! Use the standard conforming module file for omp if not using Marc +!$ use OMP_LIB, only: omp_set_num_threads +#endif + implicit none +#ifdef Marc ! use the non F90 standard include file because some versions of Marc crash when using the module +!$ include "omp_lib.h" #endif integer(pInt), parameter :: fileunit = 300_pInt ,& maxNchunks = 2_pInt @@ -121,12 +122,11 @@ subroutine numerics_init character(len=1024) :: line !$ character(len=6) DAMASK_NumThreadsString ! environment variable DAMASK_NUM_THREADS -!$OMP CRITICAL (write2out) write(6,*) write(6,*) '<<<+- numerics init -+>>>' write(6,*) '$Id$' #include "compilation_info.f90" -!$OMP END CRITICAL (write2out) + !$ call GET_ENVIRONMENT_VARIABLE(NAME='DAMASK_NUM_THREADS',VALUE=DAMASK_NumThreadsString,STATUS=gotDAMASK_NUM_THREADS) ! get environment variable DAMASK_NUM_THREADS... !$ if(gotDAMASK_NUM_THREADS /= 0) call IO_warning(47_pInt,ext_msg=DAMASK_NumThreadsString) @@ -136,12 +136,9 @@ subroutine numerics_init ! try to open the config file if(IO_open_file_stat(fileunit,numerics_configFile)) then - - !$OMP CRITICAL (write2out) - write(6,*) ' ... using values from config file' - write(6,*) - !$OMP END CRITICAL (write2out) - + + write(6,*) ' ... using values from config file' + write(6,*) !* read variables from config file and overwrite parameters @@ -229,9 +226,11 @@ subroutine numerics_init volDiscrMod_RGC = IO_floatValue(line,positions,2_pInt) case ('discrepancypower_rgc') volDiscrPow_RGC = IO_floatValue(line,positions,2_pInt) - + !* Random seeding parameters + case ('fixed_seed') + fixedSeed = IO_intValue(line,positions,2_pInt) !* spectral parameters - +#ifdef Spectral case ('err_div_tol') err_div_tol = IO_floatValue(line,positions,2_pInt) case ('err_stress_tolrel') @@ -254,12 +253,13 @@ subroutine numerics_init divergence_correction = IO_intValue(line,positions,2_pInt) > 0_pInt case ('update_gamma') update_gamma = IO_intValue(line,positions,2_pInt) > 0_pInt - - !* Random seeding parameters - - case ('fixed_seed') - fixedSeed = IO_intValue(line,positions,2_pInt) - +#endif +#ifndef Spectral + case ('err_div_tol','err_stress_tolrel','err_stress_tolabs',& + 'itmax', 'itmin','memory_efficient','fftw_timelimit','fftw_plan_mode', & + 'rotation_tol','divergence_correction','update_gamma') + call IO_warning(40_pInt,ext_msg=tag) +#endif case default call IO_error(300_pInt,ext_msg=tag) endselect @@ -268,14 +268,10 @@ subroutine numerics_init ! no config file, so we use standard values else - - !$OMP CRITICAL (write2out) - write(6,*) ' ... using standard values' - write(6,*) - !$OMP END CRITICAL (write2out) - + write(6,*) ' ... using standard values' + write(6,*) endif - +#ifdef Spectral select case(IO_lc(fftw_plan_mode)) ! setting parameters for the plan creation of FFTW. Basically a translation from fftw3.f case('estimate','fftw_estimate') ! ordered from slow execution (but fast plan creation) to fast execution fftw_planner_flag = 64_pInt @@ -289,11 +285,9 @@ subroutine numerics_init call IO_warning(warning_ID=47_pInt,ext_msg=trim(IO_lc(fftw_plan_mode))) fftw_planner_flag = 32_pInt end select - +#endif !* writing parameters to output file - - !$OMP CRITICAL (write2out) write(6,'(a24,1x,es8.1)') ' relevantStrain: ',relevantStrain write(6,'(a24,1x,es8.1)') ' defgradTolerance: ',defgradTolerance @@ -334,9 +328,14 @@ subroutine numerics_init write(6,'(a24,1x,es8.1)') ' maxVolDiscrepancy_RGC: ',maxVolDiscr_RGC write(6,'(a24,1x,es8.1)') ' volDiscrepancyMod_RGC: ',volDiscrMod_RGC write(6,'(a24,1x,es8.1,/)') ' discrepancyPower_RGC: ',volDiscrPow_RGC + !* Random seeding parameters + + write(6,'(a24,1x,i16,/)') ' fixed_seed: ',fixedSeed + !* openMP parameter + !$ write(6,'(a24,1x,i8,/)') ' number of threads: ',DAMASK_NumThreadsInt - !* spectral parameters - + !* spectral parameters +#ifdef Spectral write(6,'(a24,1x,es8.1)') ' err_div_tol: ',err_div_tol write(6,'(a24,1x,es8.1)') ' err_stress_tolrel: ',err_stress_tolrel write(6,'(a24,1x,es8.1)') ' err_stress_tolabs: ',err_stress_tolabs @@ -353,18 +352,8 @@ subroutine numerics_init write(6,'(a24,1x,es8.1)') ' rotation_tol: ',rotation_tol write(6,'(a24,1x,L8,/)') ' divergence_correction: ',divergence_correction write(6,'(a24,1x,L8,/)') ' update_gamma: ',update_gamma - - !* Random seeding parameters - - write(6,'(a24,1x,i16,/)') ' fixed_seed: ',fixedSeed - - !$OMP END CRITICAL (write2out) +#endif - - !* openMP parameter -!$ write(6,'(a24,1x,i8,/)') ' number of threads: ',DAMASK_NumThreadsInt - - !* sanity check if (relevantStrain <= 0.0_pReal) call IO_error(301_pInt,ext_msg='relevantStrain') @@ -393,7 +382,6 @@ subroutine numerics_init call IO_error(301_pInt,ext_msg='integrator') !* RGC parameters - if (absTol_RGC <= 0.0_pReal) call IO_error(301_pInt,ext_msg='absTol_RGC') if (relTol_RGC <= 0.0_pReal) call IO_error(301_pInt,ext_msg='relTol_RGC') if (absMax_RGC <= 0.0_pReal) call IO_error(301_pInt,ext_msg='absMax_RGC') @@ -409,7 +397,7 @@ subroutine numerics_init if (volDiscrPow_RGC <= 0.0_pReal) call IO_error(301_pInt,ext_msg='volDiscrPw_RGC') !* spectral parameters - +#ifdef Spectral if (err_div_tol <= 0.0_pReal) call IO_error(301_pInt,ext_msg='err_div_tol') if (err_stress_tolrel <= 0.0_pReal) call IO_error(301_pInt,ext_msg='err_stress_tolrel') if (err_stress_tolabs <= 0.0_pReal) call IO_error(301_pInt,ext_msg='err_stress_tolabs') @@ -417,10 +405,9 @@ subroutine numerics_init if (itmin > itmax) call IO_error(301_pInt,ext_msg='itmin') if (update_gamma .and. & .not. memory_efficient) call IO_error(error_ID = 847_pInt) +#endif if (fixedSeed <= 0_pInt) then - !$OMP CRITICAL (write2out) - write(6,'(a)') ' Random is random!' - !$OMP END CRITICAL (write2out) + write(6,'(a)') ' Random is random!' endif end subroutine numerics_init diff --git a/code/fftw3.f03 b/lib/fftw3.f03 similarity index 100% rename from code/fftw3.f03 rename to lib/fftw3.f03 diff --git a/code/kdtree2.f90 b/lib/kdtree2.f90 similarity index 100% rename from code/kdtree2.f90 rename to lib/kdtree2.f90 diff --git a/lib/pathinfo b/lib/pathinfo index adcf3115a..3374e938a 100644 --- a/lib/pathinfo +++ b/lib/pathinfo @@ -1,5 +1,6 @@ # possible options are MSC, FFTW, IKML, ACML, LAPACK +#IKLM /opt/intel/composerxe/mkl #ACML /opt/acml4.4.0 LAPACK /usr FFTW ./fftw