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