diff --git a/processing/post/addDivergence b/processing/post/addDivergence index bff36b615..1151ded8b 100755 --- a/processing/post/addDivergence +++ b/processing/post/addDivergence @@ -1,6 +1,6 @@ #!/usr/bin/python -import os,re,sys,math,string +import os,re,sys,math,string,numpy,postprocessingMath from optparse import OptionParser, Option # ----------------------------- @@ -60,8 +60,12 @@ Deals with both vector- and tensor-valued fields. ) -parser.add_option('-a','--accuracy', dest='accuracy', type='int', \ - help='degree of central difference accuracy [%default]') +parser.add_option('--fdm', dest='accuracy', action='extend', type='string', \ + help='degree of central difference accuracy') +parser.add_option('--fft', dest='fft', action='store_true', \ + help='calculate divergence in Fourier space [%default]') +parser.add_option('-m','--memory', dest='memory', action='store_true', \ + help='memory efficient calculation (not possible for FFT based divergency [%default]') parser.add_option('-v','--vector', dest='vector', action='extend', type='string', \ help='heading of columns containing vector field values') parser.add_option('-t','--tensor', dest='tensor', action='extend', type='string', \ @@ -71,10 +75,12 @@ parser.add_option('-d','--dimension', dest='dim', type='float', nargs=3, \ parser.add_option('-r','--resolution', dest='res', type='int', nargs=3, \ help='resolution of data set in x (fast) y z (slow)') -parser.set_defaults(accuracy = 8) +parser.set_defaults(accuracy = []) +parser.set_defaults(memory = False) +parser.set_defaults(fft = False) parser.set_defaults(vector = []) parser.set_defaults(tensor = []) -parser.set_defaults(dim = [1.0,1.0,1.0]) +parser.set_defaults(dim = []) accuracyChoices = [2,4,6,8] (options,filenames) = parser.parse_args() @@ -85,9 +91,10 @@ if len(options.dim) < 3: parser.error('improper dimension specification...') if not options.res or len(options.res) < 3: parser.error('improper resolution specification...') -if options.accuracy not in accuracyChoices: - parser.error('accuracy must be chosen from %s...'%(', '.join(accuracyChoices))) -accuracy = options.accuracy//2-1 +for choice in options.accuracy: + if int(choice) not in accuracyChoices: + parser.error('accuracy must be chosen from %s...'%(', '.join(accuracyChoices))) +if options.fft: options.accuracy.append('fft') datainfo = { # list of requested labels per datatype 'vector': {'len':3, @@ -96,7 +103,6 @@ datainfo = { # lis 'label':[]}, } - if options.vector != None: datainfo['vector']['label'] += options.vector if options.tensor != None: datainfo['tensor']['label'] += options.tensor @@ -136,6 +142,7 @@ for file in files: active = {} column = {} values = {} + div_field ={} head = [] for datatype,info in datainfo.items(): @@ -148,12 +155,18 @@ for file in files: if datatype not in active: active[datatype] = [] if datatype not in column: column[datatype] = {} if datatype not in values: values[datatype] = {} + if datatype not in div_field: div_field[datatype] = {} active[datatype].append(label) column[datatype][label] = headers.index(key) - values[datatype][label] = [[0.0 for i in range(datainfo[datatype]['len'])] \ - for j in range(options.res[0]*options.res[1]*options.res[2])] - head += prefixMultiply('div%i(%s)'%(options.accuracy,label),datainfo[datatype]['len']/3) - + values[datatype][label] = numpy.array([0.0 for i in range(datainfo[datatype]['len']*\ + options.res[0]*options.res[1]*options.res[2])]).\ + reshape((options.res[0],options.res[1],options.res[2],\ + 3,datainfo[datatype]['len']//3)) + + for what in options.accuracy: # loop over all requested degrees of accuracy (plus potentially fft) + if not options.memory or what != 'fft': # FFT divergence excluded in memory saving mode + head += prefixMultiply('div%s(%s)'%(what,label),datainfo[datatype]['len']//3) + # ------------------------------------------ assemble header --------------------------------------- output = '%i\theader'%(headerlines+1) + '\n' + \ @@ -172,40 +185,87 @@ for file in files: for datatype,labels in active.items(): for label in labels: - values[datatype][label][idx] = map(float,items[column[datatype][label]: - column[datatype][label]+datainfo[datatype]['len']]) + values[datatype][label][location(idx,options.res)[0]][location(idx,options.res)[1]][location(idx,options.res)[2]]\ + = numpy.reshape(items[column[datatype][label]: + column[datatype][label]+datainfo[datatype]['len']],(3,datainfo[datatype]['len']//3)) idx += 1 - -# ------------------------------------------ read file --------------------------------------- - - idx = 0 - for line in data: - items = line.split()[:len(headers)] - if len(items) < len(headers): - continue - - output += '\t'.join(items) +# ------------------------------------------ read file --------------------------------------- + if options.memory: + idx = 0 + for line in data: + items = line.split()[:len(headers)] + if len(items) < len(headers): + continue - (x,y,z) = location(idx,options.res) + output += '\t'.join(items) + + (x,y,z) = location(idx,options.res) + + for datatype,labels in active.items(): + for label in labels: + for accuracy in options.accuracy: + if accuracy == 'fft': continue + for k in range(datainfo[datatype]['len']/3): # formulas from Mikhail Itskov: Tensor Algebra and Tensor Analysis for Engineers, Springer 2009, p 52 + theDiv = 0.0 + for a in range(int(accuracy)//2): + + theDiv += FDcoefficients[int(accuracy)//2-1][a] * \ + ( \ + (values[datatype][label][location(index([x+1+a,y,z],options.res),options.res)[0]] \ + [location(index([x+1+a,y,z],options.res),options.res)[1]] \ + [location(index([x+1+a,y,z],options.res),options.res)[2]][k][0] - \ + values[datatype][label][location(index([x-1-a,y,z],options.res),options.res)[0]] \ + [location(index([x-1-a,y,z],options.res),options.res)[1]] \ + [location(index([x-1-a,y,z],options.res),options.res)[2]][k][0]) * options.res[0] / options.dim[0] + \ + (values[datatype][label][location(index([x,y+1+a,z],options.res),options.res)[0]] \ + [location(index([x,y+1+a,z],options.res),options.res)[1]] \ + [location(index([x,y+1+a,z],options.res),options.res)[2]][k][1] - \ + values[datatype][label][location(index([x,y-1-a,z],options.res),options.res)[0]] \ + [location(index([x,y-1-a,z],options.res),options.res)[1]] \ + [location(index([x,y-1-a,z],options.res),options.res)[2]][k][1]) * options.res[1] / options.dim[1] + \ + (values[datatype][label][location(index([x,y,z+1+a],options.res),options.res)[0]] \ + [location(index([x,y,z+1+a],options.res),options.res)[1]] \ + [location(index([x,y,z+1+a],options.res),options.res)[2]][k][2]- \ + values[datatype][label][location(index([x,y,z-1-a],options.res),options.res)[0]] \ + [location(index([x,y,z-1-a],options.res),options.res)[1]] \ + [location(index([x,y,z-1-a],options.res),options.res)[2]][k][2]) * options.res[2] / options.dim[2] \ + ) + output += '\t%f'%theDiv + + output += '\n' + idx += 1 + else: for datatype,labels in active.items(): for label in labels: - for k in range(datainfo[datatype]['len']/3): # formulas from Mikhail Itskov: Tensor Algebra and Tensor Analysis for Engineers, Springer 2009, p 52 - theDiv = 0.0 - for a in range(1+accuracy): - theDiv += FDcoefficients[accuracy][a] * \ - ( \ - (values[datatype][label][index([x+1+a,y,z],options.res)][k*3+0] - \ - values[datatype][label][index([x-1-a,y,z],options.res)][k*3+0]) * options.res[0] / options.dim[0] + \ - (values[datatype][label][index([x,y+1+a,z],options.res)][k*3+1] - \ - values[datatype][label][index([x,y-1-a,z],options.res)][k*3+1]) * options.res[1] / options.dim[1] + \ - (values[datatype][label][index([x,y,z+1+a],options.res)][k*3+2] - \ - values[datatype][label][index([x,y,z-1-a],options.res)][k*3+2]) * options.res[2] / options.dim[2] \ - ) - output += '\t%f'%theDiv + if label not in div_field[datatype]: div_field[datatype][label] = {} + for accuracy in options.accuracy: + div_field[datatype][label][accuracy] = numpy.array([0.0 for i in range((datainfo[datatype]['len'])//3*\ + options.res[0]*options.res[1]*options.res[2])]).\ + reshape((options.res[0],options.res[1],options.res[2],\ + datainfo[datatype]['len']//3)) + if accuracy == 'fft': + div_field[datatype][label][accuracy] = postprocessingMath.divergence_fft(options.res[0],options.res[1],options.res[2],datainfo[datatype]['len']//3,options.dim,values[datatype][label]) + else: + div_field[datatype][label][accuracy] = postprocessingMath.divergence(options.res[0],options.res[1],options.res[2],datainfo[datatype]['len']//3,eval(accuracy)//2-1,options.dim,values[datatype][label]) + + idx = 0 + for line in data: + items = line.split()[:len(headers)] + if len(items) < len(headers): + continue - output += '\n' - idx += 1 + output += '\t'.join(items) + + + for datatype,labels in active.items(): + for label in labels: + for accuracy in options.accuracy: + for i in range(datainfo[datatype]['len']/3): + output += '\t%f'%div_field[datatype][label][accuracy][location(idx,options.res)[0]][location(idx,options.res)[1]][location(idx,options.res)[2]][i] + output += '\n' + idx += 1 + # ------------------------------------------ output result --------------------------------------- diff --git a/processing/post/make_postprocessingMath b/processing/post/make_postprocessingMath index 5fd56384e..9441a43fa 100755 --- a/processing/post/make_postprocessingMath +++ b/processing/post/make_postprocessingMath @@ -17,7 +17,7 @@ else fi cd $wd - +rm postprocessingMath.so f2py -c --f90flags="-heap-arrays 500000000" \ postprocessingMath.pyf \ postprocessingMath.f90 \ diff --git a/processing/post/postprocessingMath.f90 b/processing/post/postprocessingMath.f90 index 9e14412ae..7ec8f445f 100644 --- a/processing/post/postprocessingMath.f90 +++ b/processing/post/postprocessingMath.f90 @@ -266,13 +266,60 @@ end function math_mul33x33 END IF RETURN END SUBROUTINE math_spectral1 + +!************************************************************************** +! volume of tetrahedron given by four vertices +!************************************************************************** + pure function math_volTetrahedron(v1,v2,v3,v4) + + implicit none + + real*8 math_volTetrahedron + real*8, dimension (3), intent(in) :: v1,v2,v3,v4 + real*8, dimension (3,3) :: m + + m(:,1) = v1-v2 + m(:,2) = v2-v3 + m(:,3) = v3-v4 + + math_volTetrahedron = math_det3x3(m)/6.0 + + end function math_volTetrahedron + +!subroutines below are for postprocessing with python + +!two small helper functions for indexing +! CAREFULL, index and location runs from 0 to N-1 (python style) + + function mesh_location(idx,resolution) + integer, intent(in) :: idx + integer, intent(in) :: resolution(3) + integer :: mesh_location(3) + mesh_location = (/modulo(idx/ resolution(3) / resolution(2),resolution(1)), & + modulo(idx/ resolution(3), resolution(2)), & + modulo(idx, resolution(3))/) + end function mesh_location + + function mesh_index(location,resolution) + integer, intent(in) :: location(3) + integer, intent(in) :: resolution(3) + integer :: mesh_index + + mesh_index = modulo(location(3), resolution(3)) +& + (modulo(location(2), resolution(2)))*resolution(3) +& + (modulo(location(1), resolution(1)))*resolution(3)*resolution(2) + end function mesh_index end module math -!subroutines below are for postprocessing with python + + + !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ subroutine mesh(res_x,res_y,res_z,geomdim,defgrad_av,centroids,nodes) !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +! Routine to build a regular mesh of cubes for given coordinates (= center of the cubes) +! implicit none real*8 geomdim(3) @@ -335,6 +382,9 @@ end subroutine mesh !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ subroutine deformed(res_x,res_y,res_z,geomdim,defgrad,defgrad_av,coord_avgCorner) !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +! Routine to calculate coordinates in current configuration for given defgrad +! using linear interpolation (blurres out high frequency defomation) +! implicit none real*8 geomdim(3) integer res_x, res_y, res_z @@ -380,6 +430,9 @@ subroutine deformed(res_x,res_y,res_z,geomdim,defgrad,defgrad_av,coord_avgCorner integer rear(3), init(3), ones(3), oppo(3), me(3), res(3) integer i, j, k, s, o + print*, 'Restore geometry using linear integration' + print '(a,/,e12.5,e12.5,e12.5)', ' Dimension:', geomdim + print '(a,/,i5,i5,i5)', ' Resolution:', res_x,res_y,res_z ones = 1 fones = 1.0 coord_avgOrder=0.0 @@ -435,7 +488,9 @@ end subroutine deformed !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ subroutine deformed_fft(res_x,res_y,res_z,geomdim,defgrad,defgrad_av,scaling,coords) !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - +! Routine to calculate coordinates in current configuration for given defgrad +! using integration in Fourier space (more accurate than deformed(...)) +! implicit none integer res_x, res_y, res_z real*8 geomdim(3) @@ -452,7 +507,11 @@ subroutine deformed_fft(res_x,res_y,res_z,geomdim,defgrad,defgrad_av,scaling,coo real*8, parameter :: pi = 3.14159265358979323846264338327950288419716939937510 integer*8 :: plan_fft(2) include 'fftw3.f' !header file for fftw3 (declaring variables). Library files are also needed - + + print*, 'Restore geometry using FFT-based integration' + print '(a,/,e12.5,e12.5,e12.5)', ' Dimension:', geomdim + print '(a,/,i5,i5,i5)', ' Resolution:', res_x,res_y,res_z + call dfftw_plan_many_dft(plan_fft(1),3,(/res_x,res_y,res_z/),9,& defgrad_fft,(/res_x,res_y,res_z/),1,res_x*res_y*res_z,& defgrad_fft,(/res_x,res_y,res_z/),1,res_x*res_y*res_z,FFTW_FORWARD,FFTW_PATIENT) @@ -499,9 +558,99 @@ subroutine deformed_fft(res_x,res_y,res_z,geomdim,defgrad,defgrad_av,scaling,coo enddo; enddo; enddo end subroutine deformed_fft + +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +subroutine volume_compare(res_x,res_y,res_z,geomdim,nodes,defgrad,volume_mismatch) +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +! Routine to calculate the mismatch between volume of reconstructed (compatible) +! cube and determinant of defgrad at the FP + + use math + implicit none + + real*8 geomdim(3) + integer res_x, res_y, res_z + real*8 nodes(res_x+1,res_y+1,res_z+1,3) + real*8 defgrad(res_x ,res_y ,res_z ,3,3) + real*8 volume_mismatch(res_x ,res_y ,res_z ) + real*8 coords(8,3) + integer i,j,k + real*8, dimension(3,3) :: defgrad_av + real*8 vol_initial + + print*, 'Calculating volume mismatch' + vol_initial = geomdim(1)*geomdim(2)*geomdim(3)/real(res_x)/real(res_y)/real(res_z) + do k = 1,res_z + do j = 1,res_y + do i = 1,res_x + coords(1,:) = nodes(i ,j ,k ,:) + coords(2,:) = nodes(i+1,j ,k ,:) + coords(3,:) = nodes(i+1,j+1,k ,:) + coords(4,:) = nodes(i ,j+1,k ,:) + coords(5,:) = nodes(i ,j, k+1,:) + coords(6,:) = nodes(i+1,j ,k+1,:) + coords(7,:) = nodes(i+1,j+1,k+1,:) + coords(8,:) = nodes(i ,j+1,k+1,:) + volume_mismatch(i,j,k) = abs(math_volTetrahedron(coords(7,:),coords(1,:),coords(8,:),coords(4,:))) & + + abs(math_volTetrahedron(coords(7,:),coords(1,:),coords(8,:),coords(5,:))) & + + abs(math_volTetrahedron(coords(7,:),coords(1,:),coords(3,:),coords(4,:))) & + + abs(math_volTetrahedron(coords(7,:),coords(1,:),coords(3,:),coords(2,:))) & + + abs(math_volTetrahedron(coords(7,:),coords(5,:),coords(2,:),coords(6,:))) & + + abs(math_volTetrahedron(coords(7,:),coords(5,:),coords(2,:),coords(1,:))) + volume_mismatch(i,j,k) = volume_mismatch(i,j,k)/math_det3x3(defgrad(i,j,k,:,:)) + enddo; enddo; enddo + volume_mismatch = volume_mismatch/vol_initial +end subroutine volume_compare + +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +subroutine shape_compare(res_x,res_y,res_z,geomdim,nodes,centroids,defgrad,shape_mismatch) +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +! 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 + implicit none + + real*8 geomdim(3) + integer res_x, res_y, res_z + real*8 nodes(res_x+1,res_y+1,res_z+1,3) + real*8 centroids(res_x ,res_y ,res_z ,3) + real*8 defgrad(res_x ,res_y ,res_z ,3,3) + real*8 shape_mismatch(res_x ,res_y ,res_z) + real*8 coords(8,3) + real*8 coords_initial(8,3) + integer i,j,k,n + + print*, 'Calculating shape mismatch' + coords_initial(1,:) = (/-geomdim(1)/2.0/real(res_x),-geomdim(2)/2.0/real(res_y),-geomdim(3)/2.0/real(res_z)/) + coords_initial(2,:) = (/+geomdim(1)/2.0/real(res_x),-geomdim(2)/2.0/real(res_y),-geomdim(3)/2.0/real(res_z)/) + coords_initial(3,:) = (/+geomdim(1)/2.0/real(res_x),+geomdim(2)/2.0/real(res_y),-geomdim(3)/2.0/real(res_z)/) + coords_initial(4,:) = (/-geomdim(1)/2.0/real(res_x),+geomdim(2)/2.0/real(res_y),-geomdim(3)/2.0/real(res_z)/) + coords_initial(5,:) = (/-geomdim(1)/2.0/real(res_x),-geomdim(2)/2.0/real(res_y),+geomdim(3)/2.0/real(res_z)/) + coords_initial(6,:) = (/+geomdim(1)/2.0/real(res_x),-geomdim(2)/2.0/real(res_y),+geomdim(3)/2.0/real(res_z)/) + coords_initial(7,:) = (/+geomdim(1)/2.0/real(res_x),+geomdim(2)/2.0/real(res_y),+geomdim(3)/2.0/real(res_z)/) + coords_initial(8,:) = (/-geomdim(1)/2.0/real(res_x),+geomdim(2)/2.0/real(res_y),+geomdim(3)/2.0/real(res_z)/) + do i=1,8 + enddo + do k = 1,res_z + do j = 1,res_y + do i = 1,res_x + shape_mismatch(i,j,k) = sqrt(sum((nodes(i ,j ,k ,:) - centroids(i,j,k,:) - matmul(defgrad(i,j,k,:,:), coords_initial(1,:)))**2.0))& + + sqrt(sum((nodes(i+1,j ,k ,:) - centroids(i,j,k,:) - matmul(defgrad(i,j,k,:,:), coords_initial(2,:)))**2.0))& + + sqrt(sum((nodes(i+1,j+1,k ,:) - centroids(i,j,k,:) - matmul(defgrad(i,j,k,:,:), coords_initial(3,:)))**2.0))& + + sqrt(sum((nodes(i ,j+1,k ,:) - centroids(i,j,k,:) - matmul(defgrad(i,j,k,:,:), coords_initial(4,:)))**2.0))& + + sqrt(sum((nodes(i ,j, k+1,:) - centroids(i,j,k,:) - matmul(defgrad(i,j,k,:,:), coords_initial(5,:)))**2.0))& + + sqrt(sum((nodes(i+1,j ,k+1,:) - centroids(i,j,k,:) - matmul(defgrad(i,j,k,:,:), coords_initial(6,:)))**2.0))& + + sqrt(sum((nodes(i+1,j+1,k+1,:) - centroids(i,j,k,:) - matmul(defgrad(i,j,k,:,:), coords_initial(7,:)))**2.0))& + + sqrt(sum((nodes(i ,j+1,k+1,:) - centroids(i,j,k,:) - matmul(defgrad(i,j,k,:,:), coords_initial(8,:)))**2.0)) + enddo; enddo; enddo + end subroutine shape_compare + !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ subroutine inverse_reconstruction(res_x,res_y,res_z,reference_configuration,current_configuration,defgrad) !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +! Routine to calculate deformation gradient from reference and current configuration +! NOT WORKING BY NOW!!!!!!!!!!!!! +! use math implicit none integer res_x, res_y, res_z @@ -600,6 +749,8 @@ end subroutine inverse_reconstruction !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ subroutine tensor_avg(res_x,res_y,res_z,tensor,avg) !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +!calculate average of tensor field +! implicit none integer res_x, res_y, res_z @@ -618,6 +769,8 @@ end subroutine tensor_avg !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ subroutine logstrain_spat(res_x,res_y,res_z,defgrad,logstrain_field) !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +!calculate logarithmic strain in spatial configuration for given defgrad field +! use math implicit none integer res_x, res_y, res_z @@ -644,7 +797,9 @@ subroutine logstrain_spat(res_x,res_y,res_z,defgrad,logstrain_field) !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ subroutine logstrain_mat(res_x,res_y,res_z,defgrad,logstrain_field) -!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +!calculate logarithmic strain in material configuration for given defgrad field +! use math implicit none integer res_x, res_y, res_z @@ -669,7 +824,9 @@ subroutine logstrain_mat(res_x,res_y,res_z,defgrad,logstrain_field) !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ subroutine calculate_cauchy(res_x,res_y,res_z,defgrad,p_stress,c_stress) -!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +!calculate cauchy stress for given PK1 stress and defgrad field +! use math implicit none integer res_x, res_y, res_z @@ -687,7 +844,9 @@ end subroutine calculate_cauchy !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ subroutine calculate_mises(res_x,res_y,res_z,tensor,vm) -!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +!calculate von Mises equivalent of tensor field +! implicit none integer res_x, res_y, res_z integer i, j, k @@ -712,3 +871,134 @@ subroutine calculate_mises(res_x,res_y,res_z,tensor,vm) vm(i,j,k,:) = sqrt(3*J_2) enddo; enddo; enddo end subroutine calculate_mises + +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +subroutine divergence_fft(res_x,res_y,res_z,vec_tens,geomdim,field,divergence_field) +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +! calculates divergence field using integration in Fourier space +!use vec_tens to decide if tensor (3) or vector (1) + + implicit none + integer res_x, res_y, res_z, vec_tens + real*8 geomdim(3) + real*8 field(res_x,res_y,res_z,vec_tens,3) + real*8 field_copy(res_x,res_y,res_z,vec_tens,3) + real*8 xi(res_x,res_y,res_z,3) + real*8 divergence_field(res_x,res_y,res_z,vec_tens) + complex*16 divergence_field_fft(res_x/2+1,res_y,res_z,vec_tens) + complex*16 field_fft(res_x,res_y,res_z,vec_tens,3) + complex*16 img + integer i, j, k + integer k_s(3) + real*8, parameter :: pi = 3.14159265358979323846264338327950288419716939937510 + integer*8 :: plan_fft(2) + include 'fftw3.f' !header file for fftw3 (declaring variables). Library files are also needed + + img = cmplx(0.0,1.0) + + call dfftw_plan_many_dft_r2c(plan_fft(1),3,(/res_x,res_y,res_z/),vec_tens*3,& + field_copy,(/res_x,res_y,res_z/),1,res_x*res_y*res_z,& + field_fft,(/res_x/2+1,res_y,res_z/),1,(res_x/2+1)*res_y*res_z,FFTW_PATIENT) + + call dfftw_plan_many_dft_c2r(plan_fft(2),3,(/res_x,res_y,res_z/),vec_tens,& + divergence_field_fft,(/res_x/2+1,res_y,res_z/),1,(res_x/2+1)*res_y*res_z,& + divergence_field,(/res_x,res_y,res_z/),1,res_x*res_y*res_z,FFTW_PATIENT) + +! field_copy is destroyed during plan creation + field_copy = field + + call dfftw_execute_dft_r2c(plan_fft(1), field_copy, field_fft) + + xi = 0.0 +! Alternative calculation of discrete frequencies k_s, ordered as in FFTW (wrap around) +! do k = 0,res_z/2 -1 + ! do j = 0,res_y/2 -1 + ! do i = 0,res_x/2 -1 + ! xi(1+mod(res_x-i,res_x),1+mod(res_y-j,res_y),1+mod(res_z-k,res_z),:) = (/-i,-j,-k/)/geomdim + ! xi(1+i, 1+mod(res_y-j,res_y),1+mod(res_z-k,res_z),:) = (/ i,-j,-k/)/geomdim + ! xi(1+mod(res_x-i,res_x),1+j, 1+mod(res_z-k,res_z),:) = (/-i, j,-k/)/geomdim + ! xi(1+i, 1+j, 1+mod(res_z-k,res_z),:) = (/ i, j,-k/)/geomdim + ! xi(1+mod(res_x-i,res_x),1+mod(res_y-j,res_y),1+k, :) = (/-i,-j, k/)/geomdim + ! xi(1+i, 1+mod(res_y-j,res_y),1+k, :) = (/ i,-j, k/)/geomdim + ! xi(1+mod(res_x-i,res_x),1+j, 1+k, :) = (/-i, j, k/)/geomdim + ! xi(1+i, 1+j, 1+k, :) = (/ i, j, k/)/geomdim + ! xi(1+i, 1+j, 1+k, :) = (/ i, j, k/)/geomdim + ! xi(1+mod(res_x-i,res_x),1+j, 1+k, :) = (/-i, j, k/)/geomdim + ! xi(1+i, 1+mod(res_y-j,res_y),1+k, :) = (/ i,-j, k/)/geomdim + ! xi(1+mod(res_x-i,res_x),1+mod(res_y-j,res_y),1+k, :) = (/-i,-j, k/)/geomdim + ! xi(1+i, 1+j, 1+mod(res_z-k,res_z),:) = (/ i, j,-k/)/geomdim + ! xi(1+mod(res_x-i,res_x),1+j, 1+mod(res_z-k,res_z),:) = (/-i, j,-k/)/geomdim + ! xi(1+i, 1+mod(res_y-j,res_y),1+mod(res_z-k,res_z),:) = (/ i,-j,-k/)/geomdim + ! xi(1+mod(res_x-i,res_x),1+mod(res_y-j,res_y),1+mod(res_z-k,res_z),:) = (/-i,-j,-k/)/geomdim + ! enddo; enddo; enddo + + do k = 0, res_z-1 + do j = 0, res_y-1 + do i = 0, res_x/2 + xi(i+1,j+1,k+1,:) = (/real(i),real(j),real(k)/)/geomdim + if(k==res_z/2) xi(i+1,j+1,k+1,3)= 0.0 ! set highest frequencies to zero + if(j==res_y/2) xi(i+1,j+1,k+1,2)= 0.0 + if(i==res_x/2) xi(i+1,j+1,k+1,1)= 0.0 + enddo; enddo; enddo + + + do k = 1, res_z + do j = 1, res_y + do i = 1, res_x/2+1 + divergence_field_fft(i,j,k,1) = sum(field_fft(i,j,k,1,:)*xi(i,j,k,:)) + if(vec_tens == 3) then + divergence_field_fft(i,j,k,2) = sum(field_fft(i,j,k,2,:)*xi(i,j,k,:)) + divergence_field_fft(i,j,k,3) = sum(field_fft(i,j,k,3,:)*xi(i,j,k,:)) + endif + enddo; enddo; enddo + divergence_field_fft = divergence_field_fft*img*2.0*pi + + call dfftw_execute_dft_c2r(plan_fft(2), divergence_field_fft, divergence_field) + +end subroutine divergence_fft + +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +subroutine divergence(res_x,res_y,res_z,vec_tens,order,geomdim,field,divergence_field) +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +! calculates divergence field using FDM with variable accuracy +!use vec_tes to decide if tensor (3) or vector (1) + + use math + implicit none + integer res_x, res_y, res_z, vec_tens, order + integer coordinates(6,3) + real*8 geomdim(3) + real*8 field(res_x,res_y,res_z,vec_tens,3) + real*8 divergence_field(res_x,res_y,res_z,vec_tens) + integer i, j, k, m, l, x + integer k_s(3) + real*8, dimension(4,4) :: FDcoefficient = reshape((/ & !from http://en.wikipedia.org/wiki/Finite_difference_coefficients + 1.0/2.0, 0.0, 0.0, 0.0,& + 2.0/3.0,-1.0/12.0, 0.0, 0.0,& + 3.0/4.0,-3.0/20.0,1.0/ 60.0, 0.0,& + 4.0/5.0,-1.0/ 5.0,4.0/105.0,-1.0/280.0/),& + (/4,4/)) + divergence_field = 0.0 + order = order + 1 + do k = 0, res_z-1; do j = 0, res_y-1; do i = 0, res_x-1 + do m = 1, order + coordinates(1,:) = mesh_location(mesh_index((/i+m,j,k/),(/res_x,res_y,res_z/)),(/res_x,res_y,res_z/)) + (/1,1,1/) + coordinates(2,:) = mesh_location(mesh_index((/i-m,j,k/),(/res_x,res_y,res_z/)),(/res_x,res_y,res_z/)) + (/1,1,1/) + coordinates(3,:) = mesh_location(mesh_index((/i,j+m,k/),(/res_x,res_y,res_z/)),(/res_x,res_y,res_z/)) + (/1,1,1/) + coordinates(4,:) = mesh_location(mesh_index((/i,j-m,k/),(/res_x,res_y,res_z/)),(/res_x,res_y,res_z/)) + (/1,1,1/) + coordinates(5,:) = mesh_location(mesh_index((/i,j,k+m/),(/res_x,res_y,res_z/)),(/res_x,res_y,res_z/)) + (/1,1,1/) + coordinates(6,:) = mesh_location(mesh_index((/i,j,k-m/),(/res_x,res_y,res_z/)),(/res_x,res_y,res_z/)) + (/1,1,1/) + do l = 1, vec_tens + divergence_field(i+1,j+1,k+1,l) = divergence_field(i+1,j+1,k+1,l) + FDcoefficient(m,order) * & + ((field(coordinates(1,1),coordinates(1,2),coordinates(1,3),l,1)- & + field(coordinates(2,1),coordinates(2,2),coordinates(2,3),l,1))*real(res_x)/geomdim(1) +& + (field(coordinates(3,1),coordinates(3,2),coordinates(3,3),l,2)- & + field(coordinates(4,1),coordinates(4,2),coordinates(4,3),l,2))*real(res_y)/geomdim(2) +& + (field(coordinates(5,1),coordinates(5,2),coordinates(5,3),l,3)- & + field(coordinates(6,1),coordinates(6,2),coordinates(6,3),l,3))*real(res_z)/geomdim(3)) + enddo + enddo + enddo; enddo; enddo +end subroutine divergence + + diff --git a/processing/post/postprocessingMath.pyf b/processing/post/postprocessingMath.pyf index 695ff600c..b29a3fcfd 100644 --- a/processing/post/postprocessingMath.pyf +++ b/processing/post/postprocessingMath.pyf @@ -46,6 +46,22 @@ python module postprocessingMath ! in real*8 dimension(3,3),intent(out) :: eb2 real*8 dimension(3,3),intent(out) :: eb3 end subroutine math_spectral1 + function math_volTetrahedron(v1,v2,v3,v4) ! in :postprocessingMath:postprocessingMath.f90:math + real*8 dimension(3), intent(in) :: v1 + real*8 dimension(3), intent(in) :: v2 + real*8 dimension(3), intent(in) :: v3 + real*8 dimension(3), intent(in) :: v4 + end function math_volTetrahedron + function mesh_location(idx,resolution) ! in :postprocessingMath:postprocessingMath.f90:math + integer, intent(in) :: idx + integer dimension(3), intent(in) :: resolution + integer dimension(3) :: mesh_location + end function mesh_location + function mesh_index(location,resolution) ! in :postprocessingMath:postprocessingMath.f90:math + integer dimension(3), intent(in) :: resolution + integer dimension(3), intent(in) :: location + integer :: mesh_index + end function mesh_location end module math subroutine mesh(res_x,res_y,res_z,geomdim,defgrad_av,centroids,nodes) ! in :postprocessingMath:postprocessingMath.f90 integer intent(in) :: res_x @@ -75,12 +91,30 @@ python module postprocessingMath ! in real*8 intent(in) :: scaling real*8 dimension(3),intent(in) :: geomdim real*8 dimension(3,3),intent(in) :: defgrad_av - real*8 dimension(res_x,res_y,res_z,3),intent(out),depend(res_x,res_y,res_z) :: coords - complex*16 dimension(res_x/2+1,res_y,res_z,3),depend(res_x,res_y,res_z) :: coords_fft real*8 dimension(res_x,res_y,res_z,3,3),intent(in),depend(res_x,res_y,res_z) :: defgrad + real*8 dimension(res_x,res_y,res_z,3),intent(out),depend(res_x,res_y,res_z) :: coords complex*16 dimension(res_x,res_y,res_z,3,3),depend(res_x,res_y,res_z) :: defgrad_fft - complex*16 dimension(res_x/2+1,res_y,res_z,3),depend(res_x,res_y,res_z) :: integrator + complex*16 dimension(res_x/2+1,res_y,res_z,3),depend(res_x,res_y,res_z) :: coords_fft end subroutine deformed_fft + subroutine volume_compare(res_x,res_y,res_z,geomdim,nodes,defgrad,volume_mismatch) ! in :postprocessingMath:postprocessingMath.f90 + integer intent(in) :: res_x + integer intent(in) :: res_y + integer intent(in) :: res_z + real*8 dimension(3),intent(in) :: geomdim + real*8 dimension(res_x+1,res_y+1,res_z+1,3),intent(in),depend(res_x,res_y,res_z) :: nodes + real*8 dimension(res_x,res_y,res_z,3,3),intent(in),depend(res_x,res_y,res_z) :: defgrad + real*8 dimension(res_x,res_y,res_z),intent(out),depend(res_x,res_y,res_z) :: volume_mismatch + end subroutine volume_compare + subroutine shape_compare(res_x,res_y,res_z,geomdim,nodes,centroids,defgrad,shape_mismatch) ! in :postprocessingMath:postprocessingMath.f90 + integer intent(in) :: res_x + integer intent(in) :: res_y + integer intent(in) :: res_z + real*8 dimension(3),intent(in) :: geomdim + real*8 dimension(res_x,res_y,res_z,3,3),intent(in),depend(res_x,res_y,res_z) :: defgrad + real*8 dimension(res_x+1,res_y+1,res_z+1,3),intent(in),depend(res_x,res_y,res_z) :: nodes + real*8 dimension(res_x,res_y,res_z,3),intent(in),depend(res_x,res_y,res_z) :: centroids + real*8 dimension(res_x,res_y,res_z),intent(out),depend(res_x,res_y,res_z) :: shape_mismatch + end subroutine shape_compare subroutine inverse_reconstruction(res_x,res_y,res_z,reference_configuration,current_configuration,defgrad) ! in :postprocessingMath:postprocessingMath.f90 integer intent(in) :: res_x integer intent(in) :: res_y @@ -125,6 +159,27 @@ python module postprocessingMath ! in real*8 dimension(res_x,res_y,res_z,3,3),intent(in),depend(res_x,res_y,res_z) :: tensor real*8 dimension(res_x,res_y,res_z,1),intent(out),depend(res_x,res_y,res_z) :: vm end subroutine calculate_mises + subroutine divergence_fft(res_x,res_y,res_z,vec_tens,geomdim,field,divergence_field) ! in :postprocessingMath:postprocessingMath.f90 + integer intent(in) :: res_x + integer intent(in) :: res_y + integer intent(in) :: res_z + integer intent(in) :: vec_tens + real*8 dimension(3),intent(in) :: geomdim + real*8 dimension(res_x,res_y,res_z,3,vec_tens),intent(in),depend(res_x,res_y,res_z,vec_tens) :: field + real*8 dimension(res_x,res_y,res_z,vec_tens), intent(out),depend(res_x,res_y,res_z,vec_tens) :: divergence_field + complex*8 dimension(res_x,res_y,res_z,3,vec_tens),depend(res_x,res_y,res_z,vec_tens) :: field_fft + complex*8 dimension(res_x/2+1,res_y,res_z,vec_tens),depend(res_x,res_y,res_z,vec_tens) :: divergence_field_fft + end subroutine divergence_fft + subroutine divergence(res_x,res_y,res_z,vec_tens,order,geomdim,field,divergence_field) ! in :postprocessingMath:postprocessingMath.f90 + integer intent(in) :: res_x + integer intent(in) :: res_y + integer intent(in) :: res_z + integer intent(in) :: vec_tens + integer intent(in) :: order + real*8 dimension(3),intent(in) :: geomdim + real*8 dimension(res_x,res_y,res_z,vec_tens,3),intent(in),depend(res_x,res_y,res_z,vec_tens,3) :: field + real*8 dimension(res_x,res_y,res_z,vec_tens),intent(out),depend(res_x,res_y,res_z,vec_tens) :: divergence_field + end subroutine divergence end interface end python module postprocessingMath diff --git a/processing/setup/setup_processing.py b/processing/setup/setup_processing.py index b7fa01a47..901e2f4a4 100755 --- a/processing/setup/setup_processing.py +++ b/processing/setup/setup_processing.py @@ -27,6 +27,7 @@ bin_link = { \ 'addMises', 'addNorm', 'addStrainTensors', + 'addDebugInformation', 'averageDown', 'mentat_colorMap', 'postResults',