diff --git a/src/discretization.f90 b/src/discretization.f90 index 9e41d8411..c90ce16ca 100644 --- a/src/discretization.f90 +++ b/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_Nodes_disp + discretization_IPcoords0, & + 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) :: & - nElem, & - nIP, & - nNodes + integer, intent(in) :: & + nElem + real(pReal), dimension(:,:), intent(in) :: & + IPcoords0, & + NodeCoords0 write(6,'(/,a)') ' <<<+- discretization init -+>>>' - discretization_nElem = nElem - discretization_nIP = nIP - - allocate(discretization_Centers_disp(3,nIP),source = 0.0_pReal) - allocate(discretization_Nodes_disp( 3,nNodes),source = 0.0_pReal) + discretization_nElem = nElem + discretization_nIP = size(IPcoords0,2) + 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') - call results_writeDataset('current',discretization_Centers_disp,'u','IP displacements','m') + real(pReal), dimension(:,:), allocatable :: u + + 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 diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index f545eab4e..b11f6ebce 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -7,14 +7,19 @@ module spectral_utilities use, intrinsic :: iso_c_binding #include use PETScSys - use prec, only: & - pReal, & - pStringLen - use math, only: & - math_I3 - + use prec + use math + use IO + use mesh + 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 diff --git a/src/mesh_grid.f90 b/src/mesh_grid.f90 index b009fe066..6034686ff 100644 --- a/src/mesh_grid.f90 +++ b/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