From 21cfe0941230605fae556673a4d1e6478ad0897a Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 27 Jun 2016 18:15:56 +0200 Subject: [PATCH] removed a bunch of unneded functions --- code/mesh.f90 | 332 +------------------------------------------------- 1 file changed, 3 insertions(+), 329 deletions(-) diff --git a/code/mesh.f90 b/code/mesh.f90 index ada304d9e..5b999cdb4 100644 --- a/code/mesh.f90 +++ b/code/mesh.f90 @@ -640,7 +640,7 @@ subroutine mesh_init(ip,el) if (myDebug) write(6,'(a)') ' Built shared elements'; flush(6) call mesh_build_ipNeighborhood #else - call mesh_spectral_build_ipNeighborhood(FILEUNIT) + call mesh_spectral_build_ipNeighborhood #endif if (myDebug) write(6,'(a)') ' Built IP neighborhood'; flush(6) @@ -1388,11 +1388,9 @@ end subroutine mesh_spectral_build_elements !> @brief build neighborhood relations for spectral !> @details assign globals: mesh_ipNeighborhood !-------------------------------------------------------------------------------------------------- -subroutine mesh_spectral_build_ipNeighborhood(fileUnit) +subroutine mesh_spectral_build_ipNeighborhood implicit none - integer(pInt), intent(in) :: & - fileUnit integer(pInt) :: & x,y,z, & e @@ -1529,332 +1527,8 @@ function mesh_nodesAroundCentres(gDim,Favg,centres) result(nodes) nodes = nodes/8.0_pReal end function mesh_nodesAroundCentres - - -!-------------------------------------------------------------------------------------------------- -!> @brief calculate coordinates in current configuration for given defgrad -! using integration in Fourier space -!-------------------------------------------------------------------------------------------------- -function mesh_deformedCoordsFFT(gDim,F,FavgIn,scalingIn) result(coords) - use IO, only: & - IO_error - use numerics, only: & - fftw_timelimit, & - fftw_planner_flag - use debug, only: & - debug_mesh, & - debug_level, & - debug_levelBasic - use math, only: & - PI, & - math_mul33x3 - - implicit none - real(pReal), intent(in), dimension(:,:,:,:,:) :: F - real(pReal), dimension(3,size(F,3),size(F,4),size(F,5)) :: coords - real(pReal), intent(in), dimension(3) :: gDim - real(pReal), intent(in), dimension(3,3), optional :: FavgIn - real(pReal), intent(in), dimension(3), optional :: scalingIn - -! allocatable arrays for fftw c routines - type(C_PTR) :: planForth, planBack - type(C_PTR) :: coords_fftw, defgrad_fftw - real(pReal), dimension(:,:,:,:,:), pointer :: F_real - complex(pReal), dimension(:,:,:,:,:), pointer :: F_fourier - real(pReal), dimension(:,:,:,:), pointer :: coords_real - complex(pReal), dimension(:,:,:,:), pointer :: coords_fourier - ! other variables - integer(pInt) :: i, j, k, m, res1Red - integer(pInt), dimension(3) :: k_s, iRes - real(pReal), dimension(3) :: scaling, step, offset_coords, integrator - real(pReal), dimension(3,3) :: Favg - integer(pInt), dimension(2:3,2) :: Nyquist ! highest frequencies to be removed (1 if even, 2 if odd) - - if (present(scalingIn)) then - where (scalingIn < 0.0_pReal) ! invalid values. in case of f2py -1 if not present - scaling = [1.0_pReal,1.0_pReal,1.0_pReal] - elsewhere - scaling = scalingIn - end where - else - scaling = 1.0_pReal - endif - - iRes = [size(F,3),size(F,4),size(F,5)] - integrator = gDim / 2.0_pReal / PI ! see notes where it is used - res1Red = iRes(1)/2_pInt + 1_pInt ! size of complex array in first dimension (c2r, r2c) - step = gDim/real(iRes, pReal) - Nyquist(2,1:2) = [iRes(2)/2_pInt + 1_pInt, iRes(2)/2_pInt + 1_pInt + mod(iRes(2),2_pInt)] - Nyquist(3,1:2) = [iRes(3)/2_pInt + 1_pInt, iRes(3)/2_pInt + 1_pInt + mod(iRes(3),2_pInt)] - -!-------------------------------------------------------------------------------------------------- -! report - if (iand(debug_level(debug_mesh),debug_levelBasic) /= 0_pInt) then - write(6,'(a)') ' Restore geometry using FFT-based integration' - write(6,'(a,3(i12 ))') ' grid a b c: ', iRes - write(6,'(a,3(es12.5))') ' size x y z: ', gDim - endif - -!-------------------------------------------------------------------------------------------------- -! sanity check - if (pReal /= C_DOUBLE .or. pInt /= C_INT) & - call IO_error(0_pInt,ext_msg='Fortran to C in mesh_deformedCoordsFFT') - -!-------------------------------------------------------------------------------------------------- -! allocation and FFTW initialization - defgrad_fftw = fftw_alloc_complex(int(res1Red *iRes(2)*iRes(3)*9_pInt,C_SIZE_T)) ! C_SIZE_T is of type integer(8) - coords_fftw = fftw_alloc_complex(int(res1Red *iRes(2)*iRes(3)*3_pInt,C_SIZE_T)) ! C_SIZE_T is of type integer(8) - call c_f_pointer(defgrad_fftw, F_real, & - [iRes(1)+2_pInt-mod(iRes(1),2_pInt),iRes(2),iRes(3),3_pInt,3_pInt]) - call c_f_pointer(defgrad_fftw, F_fourier, & - [res1Red, iRes(2),iRes(3),3_pInt,3_pInt]) - call c_f_pointer(coords_fftw, coords_real, & - [iRes(1)+2_pInt-mod(iRes(1),2_pInt),iRes(2),iRes(3),3_pInt]) - call c_f_pointer(coords_fftw, coords_fourier, & - [res1Red, iRes(2),iRes(3),3_pInt]) - - call fftw_set_timelimit(fftw_timelimit) - planForth = fftw_plan_many_dft_r2c(3_pInt,[iRes(3),iRes(2) ,iRes(1)],9_pInt,& ! dimensions , length in each dimension in reversed order - F_real,[iRes(3),iRes(2) ,iRes(1)+2_pInt-mod(iRes(1),2_pInt)],& ! input data , physical length in each dimension in reversed order - 1_pInt, iRes(3)*iRes(2)*(iRes(1)+2_pInt-mod(iRes(1),2_pInt)),& ! striding , product of physical lenght in the 3 dimensions - F_fourier,[iRes(3),iRes(2) ,res1Red],& - 1_pInt, iRes(3)*iRes(2)* res1Red,fftw_planner_flag) - - planBack = fftw_plan_many_dft_c2r(3_pInt,[iRes(3),iRes(2) ,iRes(1)],3_pInt,& - coords_fourier,[iRes(3),iRes(2) ,res1Red],& - 1_pInt, iRes(3)*iRes(2)* res1Red,& - coords_real,[iRes(3),iRes(2) ,iRes(1)+2_pInt-mod(iRes(1),2_pInt)],& - 1_pInt, iRes(3)*iRes(2)*(iRes(1)+2_pInt-mod(iRes(1),2_pInt)),& - fftw_planner_flag) - F_real(1:iRes(1),1:iRes(2),1:iRes(3),1:3,1:3) = & - reshape(F,[iRes(1),iRes(2),iRes(3),3,3], order = [4,5,1,2,3]) ! F_real is overwritten during plan creatio, is larger (padding) and has different order - -!-------------------------------------------------------------------------------------------------- -! FFT - call fftw_execute_dft_r2c(planForth, F_real, F_fourier) - -!-------------------------------------------------------------------------------------------------- -! if no average F is given, compute it in Fourier space - if (present(FavgIn)) then - if (all(FavgIn < 0.0_pReal)) then - Favg = real(F_fourier(1,1,1,1:3,1:3),pReal)/real(product(iRes),pReal) !the f2py way to tell it is not present - else - Favg = FavgIn - endif - else - Favg = real(F_fourier(1,1,1,1:3,1:3),pReal)/real(product(iRes),pReal) - endif - -!-------------------------------------------------------------------------------------------------- -! remove highest frequency in each direction, in third direction only if not 2D - - if(iRes(1)/=1_pInt) & ! do not delete the whole slice in case of 2D calculation - F_fourier (res1Red, 1:iRes(2), 1:iRes(3), 1:3,1:3) & - = cmplx(0.0_pReal,0.0_pReal,pReal) - if(iRes(2)/=1_pInt) & ! do not delete the whole slice in case of 2D calculation - F_fourier (1:res1Red,Nyquist(2,1):Nyquist(2,2),1:iRes(3), 1:3,1:3) & - = cmplx(0.0_pReal,0.0_pReal,pReal) - if(iRes(3)/=1_pInt) & ! do not delete the whole slice in case of 2D calculation - F_fourier (1:res1Red,1:iRes(2), Nyquist(3,1):Nyquist(3,2),1:3,1:3) & - = cmplx(0.0_pReal,0.0_pReal,pReal) - -!-------------------------------------------------------------------------------------------------- -! integration in Fourier space - coords_fourier = cmplx(0.0_pReal,0.0_pReal,pReal) - do k = 1_pInt, iRes(3) - k_s(3) = k-1_pInt - if(k > iRes(3)/2_pInt+1_pInt) k_s(3) = k_s(3)-iRes(3) - do j = 1_pInt, iRes(2) - k_s(2) = j-1_pInt - if(j > iRes(2)/2_pInt+1_pInt) k_s(2) = k_s(2)-iRes(2) - do i = 1_pInt, res1Red - k_s(1) = i-1_pInt - do m = 1_pInt,3_pInt - coords_fourier(i,j,k,m) = sum(F_fourier(i,j,k,m,1:3)*& - cmplx(0.0_pReal,real(k_s,pReal)*integrator,pReal)) - enddo - if (any(k_s /= 0_pInt)) coords_fourier(i,j,k,1:3) = & - coords_fourier(i,j,k,1:3) / cmplx(-sum(k_s*k_s),0.0_pReal,pReal) - enddo; enddo; enddo - -!-------------------------------------------------------------------------------------------------- -! iFFT and freeing memory - call fftw_execute_dft_c2r(planBack,coords_fourier,coords_real) - coords = reshape(coords_real(1:iRes(1),1:iRes(2),1:iRes(3),1:3), [3,iRes(1),iRes(2),iRes(3)], & - order = [2,3,4,1])/real(product(iRes),pReal) ! weight and change order - call fftw_destroy_plan(planForth) - call fftw_destroy_plan(planBack) - call fftw_free(defgrad_fftw) - call fftw_free(coords_fftw) - -!-------------------------------------------------------------------------------------------------- -! add average to scaled fluctuation and put (0,0,0) on (0,0,0) - offset_coords = math_mul33x3(F(1:3,1:3,1,1,1),step/2.0_pReal) - scaling*coords(1:3,1,1,1) - forall(k = 1_pInt:iRes(3), j = 1_pInt:iRes(2), i = 1_pInt:iRes(1)) & - coords(1:3,i,j,k) = scaling(1:3)*coords(1:3,i,j,k) & - + offset_coords & - + math_mul33x3(Favg,step*real([i,j,k]-1_pInt,pReal)) - -end function mesh_deformedCoordsFFT - - -!-------------------------------------------------------------------------------------------------- -!> @brief calculates the mismatch between volume of reconstructed (compatible) cube and -! determinant of defgrad at the FP -!-------------------------------------------------------------------------------------------------- -function mesh_volumeMismatch(gDim,F,nodes) result(vMismatch) - use IO, only: & - IO_error - use debug, only: & - debug_mesh, & - debug_level, & - debug_levelBasic - use math, only: & - math_det33, & - math_volTetrahedron - - implicit none - real(pReal), intent(in), dimension(:,:,:,:,:) :: & - F - real(pReal), dimension(size(F,3),size(F,4),size(F,5)) :: & - vMismatch - real(pReal), intent(in), dimension(:,:,:,:) :: & - nodes - real(pReal), dimension(3) :: & - gDim - integer(pInt), dimension(3) :: & - iRes - real(pReal), dimension(3,8) :: coords - integer(pInt) :: i,j,k - real(pReal) :: volInitial - - iRes = [size(F,3),size(F,4),size(F,5)] - volInitial = product(gDim)/real(product(iRes), pReal) - -!-------------------------------------------------------------------------------------------------- -! report and check - if (iand(debug_level(debug_mesh),debug_levelBasic) /= 0_pInt) then - write(6,'(a)') ' Calculating volume mismatch' - write(6,'(a,3(i12 ))') ' grid a b c: ', iRes - write(6,'(a,3(es12.5))') ' size x y z: ', gDim - endif - - if (any([iRes/=size(nodes,2)-1_pInt,iRes/=size(nodes,3)-1_pInt,iRes/=size(nodes,4)-1_pInt]))& - call IO_error(0_pInt,ext_msg='Arrays F and nodes in mesh_volumeMismatch') - -!-------------------------------------------------------------------------------------------------- -! calculate actual volume and volume resulting from deformation gradient - do k = 1_pInt,iRes(3) - do j = 1_pInt,iRes(2) - do i = 1_pInt,iRes(1) - coords(1:3,1) = nodes(1:3,i, j, k ) - coords(1:3,2) = nodes(1:3,i+1_pInt,j, k ) - coords(1:3,3) = nodes(1:3,i+1_pInt,j+1_pInt,k ) - coords(1:3,4) = nodes(1:3,i, j+1_pInt,k ) - coords(1:3,5) = nodes(1:3,i, j, k+1_pInt) - coords(1:3,6) = nodes(1:3,i+1_pInt,j, k+1_pInt) - coords(1:3,7) = nodes(1:3,i+1_pInt,j+1_pInt,k+1_pInt) - coords(1:3,8) = nodes(1:3,i, j+1_pInt,k+1_pInt) - vMismatch(i,j,k) = & - abs(math_volTetrahedron(coords(1:3,7),coords(1:3,1),coords(1:3,8),coords(1:3,4))) & - + abs(math_volTetrahedron(coords(1:3,7),coords(1:3,1),coords(1:3,8),coords(1:3,5))) & - + abs(math_volTetrahedron(coords(1:3,7),coords(1:3,1),coords(1:3,3),coords(1:3,4))) & - + abs(math_volTetrahedron(coords(1:3,7),coords(1:3,1),coords(1:3,3),coords(1:3,2))) & - + abs(math_volTetrahedron(coords(1:3,7),coords(1:3,5),coords(1:3,2),coords(1:3,6))) & - + abs(math_volTetrahedron(coords(1:3,7),coords(1:3,5),coords(1:3,2),coords(1:3,1))) - vMismatch(i,j,k) = vMismatch(i,j,k)/math_det33(F(1:3,1:3,i,j,k)) - enddo; enddo; enddo - vMismatch = vMismatch/volInitial - -end function mesh_volumeMismatch - - -!-------------------------------------------------------------------------------------------------- -!> @brief Routine to calculate the mismatch between the vectors from the central point to -! the corners of reconstructed (combatible) volume element and the vectors calculated by deforming -! the initial volume element with the current deformation gradient -!-------------------------------------------------------------------------------------------------- -function mesh_shapeMismatch(gDim,F,nodes,centres) result(sMismatch) - use IO, only: & - IO_error - use debug, only: & - debug_mesh, & - debug_level, & - debug_levelBasic - use math, only: & - math_mul33x3 - - implicit none - real(pReal), intent(in), dimension(:,:,:,:,:) :: & - F - real(pReal), dimension(size(F,3),size(F,4),size(F,5)) :: & - sMismatch - real(pReal), intent(in), dimension(:,:,:,:) :: & - nodes, & - centres - real(pReal), dimension(3) :: & - gDim, & - fRes - integer(pInt), dimension(3) :: & - iRes - real(pReal), dimension(3,8) :: coordsInitial - integer(pInt) i,j,k - - iRes = [size(F,3),size(F,4),size(F,5)] - fRes = real(iRes,pReal) - -!-------------------------------------------------------------------------------------------------- -! report and check - if (iand(debug_level(debug_mesh),debug_levelBasic) /= 0_pInt) then - write(6,'(a)') ' Calculating shape mismatch' - write(6,'(a,3(i12 ))') ' grid a b c: ', iRes - write(6,'(a,3(es12.5))') ' size x y z: ', gDim - endif - - if(any([iRes/=size(nodes,2)-1_pInt,iRes/=size(nodes,3)-1_pInt,iRes/=size(nodes,4)-1_pInt]) .or.& - any([iRes/=size(centres,2), iRes/=size(centres,3), iRes/=size(centres,4)]))& - call IO_error(0_pInt,ext_msg='Arrays F and nodes/centres in mesh_shapeMismatch') - -!-------------------------------------------------------------------------------------------------- -! initial positions - coordsInitial(1:3,1) = [-gDim(1)/fRes(1),-gDim(2)/fRes(2),-gDim(3)/fRes(3)] - coordsInitial(1:3,2) = [+gDim(1)/fRes(1),-gDim(2)/fRes(2),-gDim(3)/fRes(3)] - coordsInitial(1:3,3) = [+gDim(1)/fRes(1),+gDim(2)/fRes(2),-gDim(3)/fRes(3)] - coordsInitial(1:3,4) = [-gDim(1)/fRes(1),+gDim(2)/fRes(2),-gDim(3)/fRes(3)] - coordsInitial(1:3,5) = [-gDim(1)/fRes(1),-gDim(2)/fRes(2),+gDim(3)/fRes(3)] - coordsInitial(1:3,6) = [+gDim(1)/fRes(1),-gDim(2)/fRes(2),+gDim(3)/fRes(3)] - coordsInitial(1:3,7) = [+gDim(1)/fRes(1),+gDim(2)/fRes(2),+gDim(3)/fRes(3)] - coordsInitial(1:3,8) = [-gDim(1)/fRes(1),+gDim(2)/fRes(2),+gDim(3)/fRes(3)] - coordsInitial = coordsInitial/2.0_pReal - -!-------------------------------------------------------------------------------------------------- -! compare deformed original and deformed positions to actual positions - do k = 1_pInt,iRes(3) - do j = 1_pInt,iRes(2) - do i = 1_pInt,iRes(1) - sMismatch(i,j,k) = & - sqrt(sum((nodes(1:3,i, j, k ) - centres(1:3,i,j,k)& - - math_mul33x3(F(1:3,1:3,i,j,k), coordsInitial(1:3,1)))**2.0_pReal))& - + sqrt(sum((nodes(1:3,i+1_pInt,j, k ) - centres(1:3,i,j,k)& - - math_mul33x3(F(1:3,1:3,i,j,k), coordsInitial(1:3,2)))**2.0_pReal))& - + sqrt(sum((nodes(1:3,i+1_pInt,j+1_pInt,k ) - centres(1:3,i,j,k)& - - math_mul33x3(F(1:3,1:3,i,j,k), coordsInitial(1:3,3)))**2.0_pReal))& - + sqrt(sum((nodes(1:3,i, j+1_pInt,k ) - centres(1:3,i,j,k)& - - math_mul33x3(F(1:3,1:3,i,j,k), coordsInitial(1:3,4)))**2.0_pReal))& - + sqrt(sum((nodes(1:3,i, j, k+1_pInt) - centres(1:3,i,j,k)& - - math_mul33x3(F(1:3,1:3,i,j,k), coordsInitial(1:3,5)))**2.0_pReal))& - + sqrt(sum((nodes(1:3,i+1_pInt,j, k+1_pInt) - centres(1:3,i,j,k)& - - math_mul33x3(F(1:3,1:3,i,j,k), coordsInitial(1:3,6)))**2.0_pReal))& - + sqrt(sum((nodes(1:3,i+1_pInt,j+1_pInt,k+1_pInt) - centres(1:3,i,j,k)& - - math_mul33x3(F(1:3,1:3,i,j,k), coordsInitial(1:3,7)))**2.0_pReal))& - + sqrt(sum((nodes(1:3,i, j+1_pInt,k+1_pInt) - centres(1:3,i,j,k)& - - math_mul33x3(F(1:3,1:3,i,j,k), coordsInitial(1:3,8)))**2.0_pReal)) - enddo; enddo; enddo - -end function mesh_shapeMismatch #endif - - + #ifdef Marc4DAMASK !-------------------------------------------------------------------------------------------------- !> @brief Figures out table styles (Marc only) and stores to 'initialcondTableStyle' and