enabling regridding more than once by introducing deallocation of arrays

added J2 test (stub from Taymour)

code/DAMASK_run.py

@@ -46,7 +46,7 @@ start = 1
 exitCode=2
 print 'load case', options.loadcase
 print 'geometry', options.geometry
-res=numpy.array([24,24,24])
+res=numpy.array([32,32,32])
 while exitCode == 2:
   print 'restart at ', start
   proc=subprocess.Popen(executable='DAMASK_spectral',\
@@ -56,19 +56,21 @@ while exitCode == 2:
     myLine = proc.stdout.readline()
     if len(myLine)>1: print myLine[0:-1]   # print output without extra newline
   exitCode = proc.returncode
+  err = proc.stderr.readlines()
   if exitCode==2:
     os.system('rm -rf %i'%start)
     os.system('mkdir %i'%start)
     os.system('cp * %i/.'%start)
-    start =  int(string.split(re.search('restart at\s+\d+',proc.stderr.readlines()[-2]).group(0))[2]) #restart step is in second last line in stderr, search for restart at xxx
+    for i in xrange(len(err)):
+       if re.search('restart at\s+\d+',err[i]): start=int(string.split(err[i])[2])
 #------------regridding----------------------------------------------
 #--------------------------------------------------------------------
-    damask.core.prec.prec_init()
+    damask.core.prec.init()
-    damask.core.damask_interface.damask_interface_init(options.loadcase,options.geometry)
+    damask.core.DAMASK_interface.init(options.loadcase,options.geometry)
-    damask.core.io.io_init()
+    damask.core.IO.init()
-    damask.core.numerics.numerics_init()
+    damask.core.numerics.init()
-    damask.core.debug.debug_init()
+    damask.core.debug.init()
-    damask.core.math.math_init()
+    damask.core.math.init()
-    damask.core.fesolving.fe_init()
+    damask.core.FEsolving.init()
-    damask.core.mesh.mesh_init(1,1)
+    damask.core.mesh.init(1,1)
-    damask.core.mesh.mesh_regrid(resNewInput=res)
+    damask.core.mesh.regrid(adaptive=False,resNewInput=res)

code/IO.f90

@@ -1356,9 +1356,7 @@ subroutine IO_error(error_ID,e,i,g,ext_msg)
  case (880_pInt)
    msg = 'mismatch of microstructure count and a*b*c in geom file'
  case (890_pInt)
-   msg = 'invalid input for regridding a 2D model'
+   msg = 'invalid input for regridding'
- case (891_pInt)
-   msg = 'invalid input for regridding a 3D model'
 
 
  !* Error messages related to parsing of Abaqus input file

code/damask.core.pyf

@@ -20,15 +20,15 @@ python module core ! in
      end module prec
 
      module damask_interface ! in :damask_interface:DAMASK_spectral_interface.f90
-        subroutine damask_interface_init(loadcaseParameterIn,geometryParameterIn) ! in :damask_interface:DAMASK_spectral_interface.f90
+        subroutine DAMASK_interface_init(loadcaseParameterIn,geometryParameterIn) ! in :damask_interface:DAMASK_spectral_interface.f90
             character(len=1024), intent(in)   :: loadcaseParameterIn                                      
             character(len=1024), intent(in)   :: geometryParameterIn
-        end subroutine damask_interface_init            
+        end subroutine DAMASK_interface_init            
      end module damask_interface
 
      module io
-        subroutine io_init
+        subroutine IO_init
-        end subroutine io_init
+        end subroutine IO_init
      end module io
      
      module numerics
@@ -195,8 +195,8 @@ python module core ! in
      end module math
 
      module fesolving
-        subroutine fe_init
+        subroutine FE_init
-        end subroutine fe_init 
+        end subroutine FE_init 
      end module fesolving
 
      module mesh ! in :mesh:mesh.f90
@@ -205,7 +205,8 @@ python module core ! in
              integer, parameter             :: element = 1
          end subroutine mesh_init
 		
-        function mesh_regrid(resNewInput,minRes) ! in :mesh:mesh.f90
+         function mesh_regrid(adaptive,resNewInput,minRes) ! in :mesh:mesh.f90
+             logical, intent(in)                                  :: adaptive
              integer, dimension(3)                                :: mesh_regrid
              integer, dimension(3), intent(in), optional          :: resNewInput = -1
              integer, dimension(3), intent(in), optional          :: minRes = -1

code/mesh.f90

@@ -362,7 +362,25 @@ subroutine mesh_init(ip,element)
    write(6,*) '$Id$'
 #include "compilation_info.f90"
  !$OMP END CRITICAL (write2out)
-
+ if (allocated(mesh_mapFEtoCPelem)) deallocate(mesh_mapFEtoCPelem)
+ if (allocated(mesh_mapFEtoCPnode)) deallocate(mesh_mapFEtoCPnode)
+ if (allocated(mesh_node0)) deallocate(mesh_node0)
+ if (allocated(mesh_node)) deallocate(mesh_node)
+ if (allocated(mesh_element)) deallocate(mesh_element)
+ if (allocated(mesh_subNodeCoord)) deallocate(mesh_subNodeCoord)
+ if (allocated(mesh_ipCenterOfGravity)) deallocate(mesh_ipCenterOfGravity)
+ if (allocated(mesh_ipArea)) deallocate(mesh_ipArea)
+ if (allocated(mesh_ipAreaNormal)) deallocate(mesh_ipAreaNormal)
+ if (allocated(mesh_sharedElem)) deallocate(mesh_sharedElem)
+ if (allocated(mesh_ipNeighborhood)) deallocate(mesh_ipNeighborhood)
+ if (allocated(mesh_ipVolume)) deallocate(mesh_ipVolume)
+ if (allocated(mesh_nodeTwins)) deallocate(mesh_nodeTwins)
+ if (allocated(calcMode)) deallocate(calcMode)
+ if (allocated(FEsolving_execIP)) deallocate(FEsolving_execIP)
+ if (allocated(FE_nodesAtIP)) deallocate(FE_nodesAtIP)
+ if (allocated(FE_ipNeighbor))deallocate(FE_ipNeighbor)
+ if (allocated(FE_subNodeParent)) deallocate(FE_subNodeParent)
+ if (allocated(FE_subNodeOnIPFace)) deallocate(FE_subNodeOnIPFace)
  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...
@@ -418,9 +436,11 @@ subroutine mesh_init(ip,element)
  parallelExecution = (parallelExecution .and. (mesh_Nelems == mesh_NcpElems))                       ! plus potential killer from non-local constitutive
  
  FEsolving_execElem = [ 1_pInt,mesh_NcpElems]
+ if (allocated(FEsolving_execIP)) deallocate(FEsolving_execIP)
  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))
  
+ if (allocated(calcMode)) deallocate(calcMode)
  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"
@@ -821,7 +841,7 @@ end function mesh_spectral_getHomogenization
 !--------------------------------------------------------------------------------------------------
 !> @brief Performes a regridding from saved restart information
 !--------------------------------------------------------------------------------------------------
-function mesh_regrid(resNewInput,minRes)
+function mesh_regrid(adaptive,resNewInput,minRes)
  use prec, only: &
    pInt, &
    pReal
@@ -839,6 +859,7 @@ function mesh_regrid(resNewInput,minRes)
    deformed_FFT, &
    math_mul33x3
  character(len=1024):: formatString, N_Digits
+ logical, intent(in)                                    :: adaptive                                  ! if true, choose adaptive grid based on resNewInput, otherwise keep it constant
  integer(pInt), dimension(3), optional, intent(in)      :: resNewInput                               ! f2py cannot handle optional arguments correctly (they are always present)
  integer(pInt), dimension(3), optional, intent(in)      :: minRes
  integer(pInt), dimension(3)                            :: mesh_regrid, ratio
@@ -847,7 +868,7 @@ function mesh_regrid(resNewInput,minRes)
  integer(pInt), dimension(3)                            :: res, resNew
  integer(pInt), dimension(:),            allocatable    :: indices
  real(pReal),   dimension(3)                            :: geomdim,geomdimNew                                              
- real(pReal),   dimension(3,3)                          :: Favg 
+ real(pReal),   dimension(3,3)                          :: Favg, Favg_LastInc 
  real(pReal),   dimension(:,:),    allocatable :: & 
    coordinatesNew, &
    coordinatesLinear
@@ -857,6 +878,7 @@ function mesh_regrid(resNewInput,minRes)
    
  real(pReal),   dimension(:,:,:,:,:),    allocatable :: & 
    F,                  FNew, &
+   F_lastInc,  F_lastIncNew, &
    Fp,                FpNew, &
    Lp,                LpNew, &
    dcsdE,          dcsdENew
@@ -874,6 +896,15 @@ function mesh_regrid(resNewInput,minRes)
   res     = mesh_spectral_getResolution()
   geomdim = mesh_spectral_getDimension()
   Npoints = res(1)*res(2)*res(3)
+  
+  if (adaptive) then
+    if (present(resNewInput)) then
+      if (any (resNewInput<1)) call IO_error(890_pInt, ext_msg = 'resNewInput')
+    else
+      call IO_error(890_pInt, ext_msg = 'resNewInput')
+   endif
+  endif
+  
 !---------------------------------------------------------
   allocate(F(res(1),res(2),res(3),3,3))
   call IO_read_jobBinaryFile(777,'convergedSpectralDefgrad',trim(getSolverJobName()),size(F))
@@ -900,19 +931,28 @@ function mesh_regrid(resNewInput,minRes)
  ! ----For 2D /3D case----------------------------------                            
  if (res(3)== 1_pInt) then
    spatialDim = 2_pInt
-   if (minRes(3) > 1_pInt)  call IO_error(890_pInt)                      ! as f2py has problems with present, use pyf file for initialization to -1
+   if (present (minRes)) then
-   !check 1 and 2 for odd number, 3 should be even or 1
+     if (minRes(1) >1_pInt .and. minRes(2) > 1_pInt) then
+        if (minRes(3) /= 1_pInt .or. &
+           mod(minRes(1),2_pInt) /= 0_pInt .or. &
+           mod(minRes(2),2_pInt) /= 0_pInt)  call IO_error(890_pInt, ext_msg = '2D minRes')                      ! as f2py has problems with present, use pyf file for initialization to -1
+   endif; endif
  else
    spatialDim = 3_pInt
- !  if ( minRes(3) <= 1_pInt) call IO_error(891_pInt)
+   if (present (minRes)) then
-!check for odd numbers
+     if (all(minRes >1_pInt)) then
+        if (mod(minRes(1),2_pInt) /= 0_pInt.or. &
+            mod(minRes(2),2_pInt) /= 0_pInt .or. &
+            mod(minRes(3),2_pInt) /= 0_pInt)  call IO_error(890_pInt, ext_msg = '3D minRes')                      ! as f2py has problems with present, use pyf file for initialization to -1
+   endif; endif
  endif
    
  geomdimNew =  math_mul33x3(Favg,geomdim)
 !---- Automatic detection based on current geom -----------------
  
- if (any(resNewInput<0_pInt)) then
+ if (adaptive) then
-   ratio = floor(real(res,pReal) * (geomdimNew/geomdim), pInt)
+   ratio = floor(real(resNewInput,pReal) * (geomdimNew/geomdim), pInt)
+   
    possibleResNew = 1_pInt
  
    do i = 1_pInt, spatialDim
@@ -921,21 +961,26 @@ function mesh_regrid(resNewInput,minRes)
      else
        possibleResNew(i,1:2) = [ratio(i)-1_pInt, ratio(i) + 1_pInt]
      endif
-     if (.not.present(minRes)) then 
+     if (.not.present(minRes)) then                              ! calling from fortran, optional argument not given
        possibleResNew = possibleResNew
-     else
+     else                                                        ! optional argument is there
+       if (any(minRes<1_pInt)) then
+         possibleResNew = possibleResNew                     ! f2py calling, but without specification (or choosing invalid values), standard from pyf = -1
+       else                                                        ! given useful values
          do k = 1_pInt,3_pInt; do j = 1_pInt,3_pInt
             possibleResNew(j,k) = max(possibleResNew(j,k), minRes(j))
          enddo; enddo
        endif
+     endif
    enddo
    
   k = huge(1_pInt)
   do i = 0_pInt, 2_pInt**spatialDim - 1
-     j = possibleResNew(1,iand(i,1_pInt)/1_pInt + 1_pInt) &
+     j = abs( possibleResNew(1,iand(i,1_pInt)/1_pInt + 1_pInt) &
             * possibleResNew(2,iand(i,2_pInt)/2_pInt + 1_pInt) &
             * possibleResNew(3,iand(i,4_pInt)/4_pInt + 1_pInt) &
-       - Npoints
+            - resNewInput(1)*resNewInput(2)*resNewInput(3))
+       
     if (j < k) then 
       k = j
       resNew =[ possibleResNew(1,iand(i,1_pInt)/1_pInt + 1_pInt), &
@@ -944,11 +989,7 @@ function mesh_regrid(resNewInput,minRes)
     endif
   enddo 
  else
-   if(all(resNewInput==0_pInt)) then
    resNew = res
-   else
-     resNew = resNewInput
-   endif
  endif
  
  mesh_regrid = resNew
@@ -1012,7 +1053,6 @@ function mesh_regrid(resNewInput,minRes)
  close(777)
  
 !------Adjusting the point-to-grain association---------------------
- !write(N_Digits, '(I16.16)') 1_pInt + int(log10(real(mesh_element(4,1:Npoints),pReal)),pInt)
  write(N_Digits, '(I16.16)') 1_pInt+int(log10(real(maxval(mesh_element(4,1:Npoints)),pReal)),pInt)
  N_Digits = adjustl(N_Digits)
  formatString = '(I'//trim(N_Digits)//'.'//trim(N_Digits)//',a)'
@@ -1029,6 +1069,23 @@ function mesh_regrid(resNewInput,minRes)
  enddo
  close(777)
  
+!---------------------------------------------------------
+ allocate(F_lastInc(res(1),res(2),res(3),3,3))
+ allocate(F_lastIncNew(resNew(1),resNew(2),resNew(3),3,3))
+ call IO_read_jobBinaryFile(777,'convergedSpectralDefgrad_lastInc',trim(getSolverJobName()),size(F_lastInc))
+ read (777,rec=1) F_lastInc
+ close (777)
+ do i= 1_pInt,3_pInt; do j = 1_pInt,3_pInt
+  Favg_LastInc(i,j) = sum(F_lastInc(1:res(1),1:res(2),1:res(3),i,j)) / real(Npoints,pReal)
+ enddo; enddo
+ do k=1_pInt,resNew(3); do j=1_pInt, resNew(2); do i=1_pInt, resNew(1)
+   F_lastIncNew(i,j,k,1:3,1:3) = Favg
+ enddo; enddo; enddo
+ call IO_write_jobBinaryFile(777,'convergedSpectralDefgrad_lastInc',size(F_lastIncNew))
+ write (777,rec=1) F_lastIncNew
+ close (777)
+ deallocate(F_lastInc)
+ deallocate(F_lastIncNew)
 !-------------------------------------------------------------------
  allocate(material_phase    (1,1, Npoints))
  allocate(material_phaseNew (1,1, NpointsNew))
@@ -1060,7 +1117,6 @@ function mesh_regrid(resNewInput,minRes)
  close (777) 
  deallocate(F)
  deallocate(FNew)
- 
 !--------------------------------------------------------------------- 
  allocate(Fp       (3,3,1,1,Npoints))
  allocate(FpNew (3,3,1,1,NpointsNew))  
@@ -4126,6 +4182,13 @@ FE_ipNeighbor(1:FE_NipNeighbors(8),1:FE_Nips(8),8) = &  ! element 117
     ],pInt),[FE_NipFaceNodes,FE_NipNeighbors(10),FE_Nips(10)])
 
 end subroutine mesh_build_FEdata
-
-
 end module mesh
+
+   ! if (allocated(randInit)) deallocate(randInit)
+   ! if (allocated(FEsolving_execIP)) deallocate(FEsolving_execIP)
+   ! if (allocated(calcMode)) deallocate(calcMode)
+   ! if (allocated(FE_nodesAtIP)) deallocate(FE_nodesAtIP)
+   ! if (allocated(FE_ipNeighbor))deallocate(FE_ipNeighbor)
+   ! if (allocated(FE_subNodeParent)) deallocate(FE_subNodeParent)
+   ! if (allocated(FE_subNodeOnIPFace)) deallocate(FE_subNodeOnIPFace)
+

code/prec.f90

@@ -33,7 +33,7 @@ module prec
 !    *** Precision of real and integer variables ***
  integer,     parameter, public :: pReal = selected_real_kind(15,300)     ! 15 significant digits, up to 1e+-300
  integer,     parameter, public :: pInt  = selected_int_kind(9)           ! up to +- 1e9
- integer,     parameter, public :: pLongInt  = 8                          ! should be 64bit
+ integer,     parameter, public :: pLongInt  = selected_int_kind(12)      ! should be 64bit
  real(pReal), parameter, public :: tol_math_check = 1.0e-8_pReal
  real(pReal), parameter, public :: tol_gravityNodePos = 1.0e-100_pReal
  

lib/damask/__init__.py

@@ -13,6 +13,35 @@ from util import *
 
 try:
   from .          import core
+# cleaning up namespace
+###################################################################################################
+# capitalize according to convention
+  core.IO                            = core.io
+  core.FEsolving                     = core.fesolving
+  core.DAMASK_interface              = core.damask_interface
+# remove XXX_
+  core.prec.init                     = core.prec.prec_init
+  core.DAMASK_interface.init         = core.DAMASK_interface.DAMASK_interface_init
+  core.IO.init                       = core.IO.IO_init
+  core.numerics.init                 = core.numerics.numerics_init
+  core.debug.init                    = core.debug.debug_init
+  core.math.init                     = core.math.math_init
+  #core.math.volumeMismatch           = core.math.math_volumeMismatch
+  #core.math.shapeMismatch            = core.math.math_shapeMismatch
+  #core.math.deformedCoordsLin        = core.math.math_deformedCoordsLin
+  #core.math.deformedCoordsFFT        = core.math.math_deformedCoordsFFT
+  #core.math.curlFFT                  = core.math.math_curlFFT
+  #core.math.divergenceFFT            = core.math.math_divergenceFFT
+  #core.math.divergenceFDM            = core.math.math_divergenceFDM
+  #core.math.tensorAvg                = core.math.math_tensorAvg
+  #core.math.logStrainSpat            = core.math.math_logStrainSpat
+  #core.math.logStrainMat             = core.math.math_logStrainMat
+  #core.math.cauchyStress             = core.math.math_cauchyStress
+  #core.math.periodicNearestNeighbor  = core.math.math_periodicNearestNeighbor
+  core.FEsolving.init                = core.FEsolving.FE_init
+  core.mesh.init                     = core.mesh.mesh_init
+  core.mesh.regrid                   = core.mesh.mesh_regrid
+  #core.mesh.nodesAroundCentroids     = core.mesh.mesh_spectral_nodesAroundCentroids
 except ImportError:
   sys.stderr.write('\nWARNING: Core module (Fortran code) not available, try to run setup_processing.py\nError Message when importing core.so: \n\n')
   core = None # from http://www.python.org/dev/peps/pep-0008/