From 60a3ac5b0462dc6533c5cc27253ea5c9fd5335c1 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 24 Mar 2016 12:35:33 +0100 Subject: [PATCH] copied fortran code --- processing/post/addCompatibilityMismatch.py | 145 +++++++++++++++++++- 1 file changed, 144 insertions(+), 1 deletion(-) diff --git a/processing/post/addCompatibilityMismatch.py b/processing/post/addCompatibilityMismatch.py index f237f7d19..20cdd284c 100755 --- a/processing/post/addCompatibilityMismatch.py +++ b/processing/post/addCompatibilityMismatch.py @@ -9,6 +9,150 @@ import damask scriptName = os.path.splitext(os.path.basename(__file__))[0] scriptID = ' '.join([scriptName,damask.version]) + +def volTetrahedron(vertices=None, sides=None): + """ + Return the volume of the tetrahedron with given vertices or sides. If + vertices are given they must be in a NumPy array with shape (4,3): the + position vectors of the 4 vertices in 3 dimensions; if the six sides are + given, they must be an array of length 6. If both are given, the sides + will be used in the calculation. + + This method implements + Tartaglia's formula using the Cayley-Menger determinant: + |0 1 1 1 1 | + |1 0 s1^2 s2^2 s3^2| + 288 V^2 = |1 s1^2 0 s4^2 s5^2| + |1 s2^2 s4^2 0 s6^2| + |1 s3^2 s5^2 s6^2 0 | + where s1, s2, ..., s6 are the tetrahedron side lengths. + + from http://codereview.stackexchange.com/questions/77593/calculating-the-volume-of-a-tetrahedron + + """ + + # The indexes of rows in the vertices array corresponding to all + # possible pairs of vertices + vertex_pair_indexes = np.array(((0, 1), (0, 2), (0, 3), + (1, 2), (1, 3), (2, 3))) + + # Get all the squares of all side lengths from the differences between + # the 6 different pairs of vertex positions + vertex1, vertex2 = vertex_pair_indexes[:,0], vertex_pair_indexes[:,1] + sides_squared = np.sum((vertices[vertex1] - vertices[vertex2])**2,axis=-1) + + + # Set up the Cayley-Menger determinant + M = np.zeros((5,5)) + # Fill in the upper triangle of the matrix + M[0,1:] = 1 + # The squared-side length elements can be indexed using the vertex + # pair indices (compare with the determinant illustrated above) + M[tuple(zip(*(vertex_pair_indexes + 1)))] = sides_squared + + # The matrix is symmetric, so we can fill in the lower triangle by + # adding the transpose + M = M + M.T + return np.sqrt(det / 288) + + +def mesh_volumeMismatch(size,F,nodes): + """ + calculates the mismatch between volume of reconstructed (compatible) cube and + determinant of defgrad at the FP + """ + 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,8) :: coords + + volInitial = size.prod()/grid.prod() + +#-------------------------------------------------------------------------------------------------- +# calculate actual volume and volume resulting from deformation gradient + for k in xrange(grid[2]): + for j in xrange(grid[1]): + for i in xrange(grid[0]): + coords(0:3,0) = nodes[0:3,i, j, k ] + coords(0:3,1) = nodes[0:3,i+1,j, k ] + coords(0:3,2) = nodes[0:3,i+1,j+1,k ] + coords(0:3,3) = nodes[0:3,i, j+1,k ] + coords(0:3,4) = nodes[0:3,i, j, k+1] + coords(0:3,5) = nodes[0:3,i+1,j, k+1] + coords(0:3,6) = nodes[0:3,i+1,j+1,k+1] + coords(0:3,7) = nodes[0:3,i, j+1,k+1] + vMismatch[i,j,k] = & + abs(volTetrahedron(coords[0:3,6],coords[0:3,0],coords[0:3,7],coords[0:3,3])) & + + abs(volTetrahedron(coords[0:3,6],coords[0:3,0],coords[0:3,7],coords[0:3,4])) & + + abs(volTetrahedron(coords[0:3,6],coords[0:3,0],coords[0:3,2],coords[0:3,3])) & + + abs(volTetrahedron(coords[0:3,6],coords[0:3,0],coords[0:3,2],coords[0:3,1])) & + + abs(volTetrahedron(coords[0:3,6],coords[0:3,4],coords[0:3,1],coords[0:3,5])) & + + abs(volTetrahedron(coords[0:3,6],coords[0:3,4],coords[0:3,1],coords[0:3,0])) + vMismatch[i,j,k] = vMismatch[i,j,k]/math_det33(F(1:3,1:3,i,j,k)) + enddo; enddo; enddo + + return vMismatch/volInitial + + + +def mesh_shapeMismatch(gDim,F,nodes,centres): + """ + 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(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,8) :: coordsInitial + integer(pInt) i,j,k + +!-------------------------------------------------------------------------------------------------- +! 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 + return sMismatch + # -------------------------------------------------------------------- # MAIN # -------------------------------------------------------------------- @@ -117,7 +261,6 @@ for name in filenames: idx += 1 F[0:3,0:3,x,y,z] = np.array(map(float,table.data[column:column+9]),'d').reshape(3,3) - Favg = damask.core.math.tensorAvg(F) centres = damask.core.mesh.deformedCoordsFFT(size,F,Favg,[1.0,1.0,1.0]) nodes = damask.core.mesh.nodesAroundCentres(size,Favg,centres)