Skip to content
Snippets Groups Projects
Liutex 4.07 KiB
Newer Older
Armand Duplex Guetse Tiogo's avatar
Armand Duplex Guetse Tiogo committed
# import the necessary libraries
import numpy as np
import pyvista as pv
from scipy.linalg import schur
from pyvista import examples

#Load a dataset
mesh = examples.download_carotid()
mesh

#Compute the gradients of the vectors (assuming you have a vector field in the point data of the mesh)
mesh_g = mesh.compute_derivative(scalars="vectors")



#Compute the gradients of the vectors (assuming you have a scalar field in the point data of the mesh)
mesh_g_s = mesh.compute_derivative(scalars="scalars")


#Access the gradient values for a vector field
gradients = mesh_g["gradient"]

"""
internal_mesh = gradients['internalMesh']
boundary_mesh = gradients['boundary']['movingWall']
"""
#Access the gradient values for a scalar field
gradients_s = mesh_g_s["gradient"]

#Organize the gradients into a dictionary for easier access
def gradients_to_dict(arr):
    """A helper method to label the gradients into a dictionary."""
    keys = np.array(
        ["du/dx", "du/dy", "du/dz", "dv/dx", "dv/dy", "dv/dz", "dw/dx", "dw/dy", "dw/dz"],
    )
    keys = keys.reshape((3, 3))[:, : arr.shape[1]].ravel()
    return dict(zip(keys, arr.T))

#for a vector field
gradients_dict = gradients_to_dict(gradients)

#for a scalar field
gradients_dict_s = gradients_to_dict(gradients_s)



# And we can add all of those components as individual arrays back to the mesh for a vector
# by:
mesh_g.point_data.update(gradients_dict)





# And we can add all of those components as individual arrays back to the mesh for a scalar
# by:
mesh_g_s.point_data.update(gradients_dict_s)

# Visualize the gradients
#keys for a vector
keys = np.array(list(gradients_dict.keys())).reshape(3, 3)

#keys for a scalar
keys_s = np.array(list(gradients_dict_s.keys())).reshape(1, 3)

# Transpose the matric keys
#A = np.transpose(keys, axes=None)
A = np.transpose(keys)


#Visualize the gradients of a vector
p = pv.Plotter(shape=keys.shape)
for (i, j), name in np.ndenumerate(keys):
    p.subplot(i, j)
    p.add_mesh(mesh_g.contour(scalars=name), scalars=name, opacity=0.75)
    p.add_mesh(mesh_g.outline(), color="k")
p.link_views()
p.view_isometric()
p.show()

points = pv.PointSet(mesh_g.points)
points

gradients_data = gradients.reshape(-1, 3)

#Visualize the transpose of the gradients of a vector
p = pv.Plotter(shape=A.shape)
for (i, j), name in np.ndenumerate(A):
    p.subplot(i, j)
    p.add_mesh(mesh_g.contour(scalars=name), scalars=name, opacity=0.75)
    p.add_mesh(mesh_g.outline(), color="k")
p.link_views()
p.view_isometric()
p.show()

#Print all the values of gradients_dict
print(gradients_dict.values()) 

#Print the items of gradients_dict
items = gradients_dict.items()

print(items) 
print(len(gradients_dict)) 

#get the value of each velocity  "du/dx", "du/dy", "du/dz" ....
value_vect = gradients_dict.get("du/dx")     
value_vect_1 = gradients_dict.get("du/dy")
value_vect_2 = gradients_dict.get("du/dz")
value_vect_3 = gradients_dict.get("dv/dx")
value_vect_4 = gradients_dict.get("dv/dy")
value_vect_5 = gradients_dict.get("dv/dz")
value_vect_6 = gradients_dict.get("dw/dx")
value_vect_7 = gradients_dict.get("dw/dy")
value_vect_8 = gradients_dict.get("dw/dz")

#print the value of "du/dx"
print(value_vect)

#change the gradient_dict to a matrix 
def create_matrix_from_parameters(params):
    if len(params)!=9:
        raise ValueError("Sie müssen genau 9 Parameter eingeben.")
        
    matrix = []
    for i in range(3):
        row=params[i*3:(i+1)*3]
        matrix.append(row)
    return matrix


  
#params = [value_vect, value_vect_1, value_vect_2, value_vect_3, value_vect_4, value_vect_5, value_vect_6, value_vect_7, value_vect_8]

params_trans = [value_vect, value_vect_3, value_vect_6, value_vect_1, value_vect_4, value_vect_7, value_vect_3, value_vect_5, value_vect_8]

matrix = create_matrix_from_parameters(params_trans)

print("Matrix 3x3", matrix )

# Compute the real Schur decomposition
Q_, S = schur(matrix, output='real')

print ("Matrix Q*", Q_)

print("Determinant of Q*", np.linalg.det(Q_))