0

I am solving this problem in python that consist in to obtain the total volume of a tetrahedral mesh of a knight.

This code that I wrote obtains the value of the volume but of the first tetrahedron of the mesh, the problem is that I do not know how to put it so that it reads all the tetrahedra and adds each of the volumes of this one. I know it is with a for loop but I don't know where to put it.

For a better understanding, here is attached a photo of the "tets" and "pts" matrices, the idea of the "tets" matrix is to store the row of "pts" where the coordinates of each of the 4 points are found that make up a tetrahedron.

Image 1

Image 2

And this is the code:

# -- coding: utf-8 --
"""
Editor de Spyder

Este es un archivo temporal.
"""
import numpy as np

import meshio

mesh = meshio.read("knight.msh")#read the mesh
pts = mesh.points #stores the mesh points in an array
tets = mesh.cells[0].data #"tets" stores the number of the row of "pts" in which the coordinates of the tetrahedron are found


#VERTICE A OF THE TETRAHEDRO
kA_i  = tets[0,0] #Take from the "test" matrix the value [n, 0] that corresponds to the place of "pts" #where the x, y, z coordinates of the point are stored


A = (pts[kA_i,0], pts[kA_i,1],pts[kA_i,2]) #Stores the coordinates of vertice A

#VERTICE B OF THE TETRAHEDRO
kB_i = tets[0,1]

B = (pts[kB_i,0],pts[kB_i,1],pts[kB_i,2])

#VERTICE C OF THE TETRAHEDRO
kC_i = tets[0,2]

C = (pts[kC_i,0],pts[kC_i,1],pts[kC_i,2])

#VERTICE D OF THE TETRAHEDRO
kD_i = tets[0,3]

D = (pts[kD_i,0],pts[kD_i,1],pts[kD_i,2])

#CALCULATION OF THE TETRAHEDRAL SEGMENTS

SAB = (B[0]-A[0],B[*1*]-A[*1*],B[*2*]-A[*2*]) #SEGMENT AB

U = tuple(SAB) #Directional vector U

SAC = (C[0]-A[0],C[*1*]-A[*1*],C[*2*]-A[*2*]) #SEGMENT AC

V = tuple(SAC) #Directional vector V

SAD = (D[0]-A[0],D[*1*]-A[*1*],D[*2*]-A[*2*]) #SEGMENT AD

W = tuple(SAD) #Directional vector W

#MATRIX GENERATION

M = np.array((U,V,W)) 

#CALCULATION OF THE VOLUME OF THE TETRAHEDRON: 

Vol = (1/6)*np.abs((M[0,0]*(M[1,1]*M[2,2] - M[2,1]*M[1,2]) - M[0,1]*(M[1,0]*M[2,2] - M[2,0]*M[1,2]) + M[0,2]*(M[1,0]*M[2,1] - M[2,0]*M[1,1])))
khelwood
  • 55,782
  • 14
  • 81
  • 108
Juan Alberto
  • 3
  • 1

1 Answers1

0

Make a procedure to compute the volume of a single tetrahedron

def tetra_vol(nodes):
    M = np.array( [ 1,*pts[node]] for node in nodes)
    return abs(np.linalg.det(M))/6

Feel free to replace this with the single steps that you used which compute exactly the same determinant value. The operations you employ are row operations that can be used to simplify the 4x4 determinant to a 3x3 determinant, and of that you used an explicit formula.

Then compute the sum over all tetrahedrons

Vol = sum(tetra_vol(tet) for tet in tets)
Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51