easier to store positions than displacement

src/discretization.f90

@@ -11,43 +11,65 @@ module discretization
   private
   
   integer, public, protected :: &
-    discretization_nElem, &
+    discretization_nIP, &
-    discretization_nIP
+    discretization_nElem
+
 
   real(pReal), dimension(:,:), allocatable :: &
-    discretization_Centers_disp, &
+    discretization_IPcoords0, &
-    discretization_Nodes_disp
+    discretization_NodeCoords0, &
+    discretization_IPcoords, &
+    discretization_NodeCoords
 
   public :: &
     discretization_init, &
-    discretization_results
+    discretization_results, &
+    discretization_setIPcoords
 
 contains
   
 
-subroutine discretization_init(nElem,nIP,nNodes)
+subroutine discretization_init(nElem,IPcoords0,NodeCoords0)
 
-  integer, intent(in) :: &
+  integer,                     intent(in) :: &
-    nElem, &
+    nElem
-    nIP, &
+  real(pReal), dimension(:,:), intent(in) :: &
-    nNodes
+    IPcoords0, &
+    NodeCoords0
 
   write(6,'(/,a)')   ' <<<+-  discretization init  -+>>>'
    
-  discretization_nElem = nElem
+  discretization_nElem = nElem 
-  discretization_nIP   = nIP
+  discretization_nIP   = size(IPcoords0,2)
-
-  allocate(discretization_Centers_disp(3,nIP),source    = 0.0_pReal)
-  allocate(discretization_Nodes_disp(  3,nNodes),source = 0.0_pReal)
 
+  discretization_IPcoords0   = IPcoords0
+  discretization_IPcoords    = IPcoords0
+  discretization_NodeCoords0 = NodeCoords0
+  discretization_NodeCoords  = NodeCoords0
+  
 end subroutine discretization_init
 
 
 subroutine discretization_results
 
-  call results_writeDataset('current',discretization_Nodes_disp,'U',  'nodal displacements','m')
+  real(pReal), dimension(:,:), allocatable :: u
-  call results_writeDataset('current',discretization_Centers_disp,'u','IP displacements','m')
+  
+  u =  discretization_NodeCoords -discretization_NodeCoords0
+  call results_writeDataset('current',U,'U','nodal displacements','m')
+  
+  u = discretization_IPcoords -discretization_IPcoords0
+  call results_writeDataset('current',u,'u','IP displacements','m')
 
 end subroutine discretization_results
 
+
+subroutine discretization_setIPcoords(IPcoords)
+
+  real(pReal), dimension(:,:), intent(in) :: IPcoords
+  
+  discretization_IPcoords = IPcoords
+
+end subroutine discretization_setIPcoords
+
+
 end module discretization

src/grid/spectral_utilities.f90

@@ -7,14 +7,19 @@ module spectral_utilities
   use, intrinsic :: iso_c_binding
 #include <petsc/finclude/petscsys.h>
   use PETScSys
-  use prec, only: &
+  use prec
-    pReal, &
+  use math
-    pStringLen
+  use IO
-  use math, only: &
+  use mesh
-   math_I3
+  use numerics
-
+  use debug
+  use config
+  use discretization
+  use homogenization
+  
   implicit none
   private
+  
   include 'fftw3-mpi.f03'
 
   logical, public             :: cutBack = .false.                                                  !< cut back of BVP solver in case convergence is not achieved or a material point is terminally ill
@@ -170,32 +175,7 @@ contains
 !> Initializes FFTW.
 !--------------------------------------------------------------------------------------------------
 subroutine utilities_init
-  use IO, only: &
+
-    IO_error, &
-    IO_warning, &
-    IO_lc
-  use numerics, only: &
-    petsc_defaultOptions, &
-    petsc_options
-  use debug, only: &
-    debug_level, &
-    debug_SPECTRAL, &
-    debug_LEVELBASIC, &
-    debug_SPECTRALDIVERGENCE, &
-    debug_SPECTRALFFTW, &
-    debug_SPECTRALPETSC, &
-    debug_SPECTRALROTATION
-  use config, only: &
-    config_numerics
-  use debug, only: &
-    PETSCDEBUG
-  use math
-  use mesh, only: &
-    grid, &
-    grid3, &
-    grid3Offset, &
-    geomSize
- 
   PetscErrorCode :: ierr
   integer        :: i, j, k, &
     FFTW_planner_flag
@@ -412,17 +392,6 @@ end subroutine utilities_init
 !> Also writes out the current reference stiffness for restart.
 !---------------------------------------------------------------------------------------------------
 subroutine utilities_updateGamma(C,saveReference)
-  use IO, only: &
-    IO_open_jobFile_binary
-  use numerics, only: &
-    worldrank
-  use mesh, only: &
-    grid3Offset, &
-    grid3,&
-    grid
-  use math, only: &
-    math_det33, &
-    math_invert2
  
   real(pReal), intent(in), dimension(3,3,3,3) :: C                                                  !< input stiffness to store as reference stiffness
   logical    , intent(in)                     :: saveReference                                      !< save reference stiffness to file for restart
@@ -538,13 +507,6 @@ end subroutine utilities_FFTvectorBackward
 !> @brief doing convolution gamma_hat * field_real, ensuring that average value = fieldAim
 !--------------------------------------------------------------------------------------------------
 subroutine utilities_fourierGammaConvolution(fieldAim)
-  use math, only: &
-    math_det33, &
-    math_invert2
-  use mesh, only: &
-    grid3, &
-    grid, &
-    grid3Offset
  
   real(pReal), intent(in), dimension(3,3) :: fieldAim                                               !< desired average value of the field after convolution
   complex(pReal),          dimension(3,3) :: temp33_complex, xiDyad_cmplx
@@ -600,12 +562,7 @@ end subroutine utilities_fourierGammaConvolution
 !> @brief doing convolution DamageGreenOp_hat * field_real
 !--------------------------------------------------------------------------------------------------
 subroutine utilities_fourierGreenConvolution(D_ref, mobility_ref, deltaT)
-  use math, only: &
+
-    PI
-  use mesh, only: &
-    grid, &
-    grid3
- 
   real(pReal), dimension(3,3), intent(in) :: D_ref
   real(pReal),                 intent(in) :: mobility_ref, deltaT
   complex(pReal)                          :: GreenOp_hat
@@ -676,13 +633,7 @@ end function utilities_divergenceRMS
 !> @brief calculate max of curl of field_fourier
 !--------------------------------------------------------------------------------------------------
 real(pReal) function utilities_curlRMS()
-  use IO, only: &
+
-    IO_error
-  use mesh, only: &
-    geomSize, &
-    grid, &
-    grid3
- 
   integer  ::  i, j, k, l, ierr
   complex(pReal), dimension(3,3) :: curl_fourier
   complex(pReal), dimension(3)   :: rescaledGeom
@@ -743,17 +694,7 @@ end function utilities_curlRMS
 !> @brief calculates mask compliance tensor used to adjust F to fullfill stress BC
 !--------------------------------------------------------------------------------------------------
 function utilities_maskedCompliance(rot_BC,mask_stress,C)
-  use, intrinsic :: &
+
-    IEEE_arithmetic
-  use IO, only: &
-    IO_error
-  use math, only: &
-    math_3333to99, &
-    math_99to3333, &
-    math_rotate_forward3333, &
-    math_rotate_forward33, &
-    math_invert2
- 
   real(pReal),              dimension(3,3,3,3) :: utilities_maskedCompliance                        !< masked compliance
   real(pReal), intent(in) , dimension(3,3,3,3) :: C                                                 !< current average stiffness
   real(pReal), intent(in) , dimension(3,3)     :: rot_BC                                            !< rotation of load frame
@@ -844,9 +785,6 @@ end function utilities_maskedCompliance
 !> @brief calculate scalar gradient in fourier field
 !--------------------------------------------------------------------------------------------------
 subroutine utilities_fourierScalarGradient()
-  use mesh, only: &
-    grid3, &
-    grid
  
   integer :: i, j, k
  
@@ -861,9 +799,6 @@ end subroutine utilities_fourierScalarGradient
 !> @brief calculate vector divergence in fourier field
 !--------------------------------------------------------------------------------------------------
 subroutine utilities_fourierVectorDivergence()
-  use mesh, only: &
-    grid3, &
-    grid
  
   integer                :: i, j, k
  
@@ -879,9 +814,6 @@ end subroutine utilities_fourierVectorDivergence
 !> @brief calculate vector gradient in fourier field
 !--------------------------------------------------------------------------------------------------
 subroutine utilities_fourierVectorGradient()
-  use mesh, only: &
-    grid3, &
-    grid
  
   integer :: i, j, k, m, n
  
@@ -899,10 +831,7 @@ end subroutine utilities_fourierVectorGradient
 !> @brief calculate tensor divergence in fourier field
 !--------------------------------------------------------------------------------------------------
 subroutine utilities_fourierTensorDivergence()
-  use mesh, only: &
+
-    grid3, &
-    grid
- 
   integer :: i, j, k, m, n
  
   vectorField_fourier = cmplx(0.0_pReal,0.0_pReal,pReal)
@@ -921,21 +850,6 @@ end subroutine utilities_fourierTensorDivergence
 !--------------------------------------------------------------------------------------------------
 subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
                                           F,timeinc,rotation_BC)
-  use IO, only: &
-    IO_error
-  use numerics, only: &
-    worldrank
-  use math, only: &
-    math_rotate_forward33, &
-    math_det33
-  use mesh, only: &
-    grid,&
-    grid3
-  use homogenization, only: &
-    materialpoint_F, &
-    materialpoint_P, &
-    materialpoint_dPdF, &
-    materialpoint_stressAndItsTangent
  
   real(pReal),intent(out), dimension(3,3,3,3)                     :: C_volAvg, C_minmaxAvg          !< average stiffness
   real(pReal),intent(out), dimension(3,3)                         :: P_av                           !< average PK stress
@@ -1127,16 +1041,6 @@ end function utilities_getFreqDerivative
 ! convolution
 !--------------------------------------------------------------------------------------------------
 subroutine utilities_updateIPcoords(F)
-  use prec, only: &
-    cNeq
-  use IO, only: &
-    IO_error
-  use mesh, only: &
-    grid, &
-    grid3, &
-    grid3Offset, &
-    geomSize, &
-    mesh_ipCoordinates
  
   real(pReal),   dimension(3,3,grid(1),grid(2),grid3), intent(in) :: F
   integer :: i, j, k, m, ierr
@@ -1178,6 +1082,8 @@ subroutine utilities_updateIPcoords(F)
                                 + matmul(Favg,step*real([i,j,k+grid3Offset]-1,pReal))
     m = m+1
   enddo; enddo; enddo
+  
+  call discretization_setIPcoords(reshape(mesh_ipCoordinates,[3,grid(1)*grid(2)*grid3]))
 
 end subroutine utilities_updateIPcoords
 

src/mesh_grid.f90

@@ -197,7 +197,9 @@ subroutine mesh_init(ip,el)
  theMesh%homogenizationAt  = mesh_element(3,:)
  theMesh%microstructureAt  = mesh_element(4,:)
 !!!!!!!!!!!!!!!!!!!!!!!!
-  call discretization_init(product(grid),product(grid),mesh_Nnodes)
+  call discretization_init(product(grid), &
+                           reshape(mesh_ipCoordinates,[3,grid(1)*grid(2)*grid3]),&
+                           mesh_node0)
 
 end subroutine mesh_init