Problem for computing the L matrix using SPM method

[BhaskarIlla]

Dear Prof. Giroux,

I am facing the problem for computing the L matrix using SPM. I tried couple of approaches but not succeed. Please find the attached code and generated errors information below.

1. "ValueError: Slowness vector has wrong size" occurred if I keep the cell_slowness as a 'False'

2."NotImplementedError: compute_L not implemented for mesh with slowness defined in cells" occurred if I keep the cell_slowness as a 'True'

code:
import gmsh
import numpy as np
from ttcrpy.tmesh import Mesh3d,Mesh2d
from ttcrpy import tmesh
import vtk
from vtk.util.numpy_support import vtk_to_numpy
import math
from scipy.sparse import csr_matrix, save_npz

gmsh.initialize()
mesh_file = "smesh_3d_new.msh"
gmsh.merge(mesh_file)
gmsh.model.mesh.renumberNodes()
gmsh.model.mesh.renumberElements()
num_elements = gmsh.model.mesh.getElements()
num_nodes = gmsh.model.mesh.getNodes()[1].shape[0]
physical_groups = gmsh.model.getPhysicalGroups()
for dim, tag in physical_groups:
name = gmsh.model.getPhysicalName(dim, tag)
print(f"Physical group (Dim {dim}, Tag {tag}): {name}")
entities = gmsh.model.getEntities()

slowness,mgrids=[],[]

for dim, tag in gmsh.model.getEntities():
elemTypes, elemTags, elemNodeTags = gmsh.model.mesh.getElements(dim, tag)
physicalTags = gmsh.model.getPhysicalGroupsForEntity(dim, tag)
if dim == 3:
elemTypes, elemTags, elemNodeTags = gmsh.model.mesh.getElements(dim, tag)
physicalTags = gmsh.model.getPhysicalGroupsForEntity(dim, tag)
name = gmsh.model.getPhysicalName(dim, physicalTags[0])
#print(name)
for n in range(len(elemTags[0])):
t = elemNodeTags[0][4n:(4n+4)]
mgrids.append(t)
if name == "volume7":
slowness.append(1/3.5)
elif name == "volume6":
slowness.append(1/3.65)
elif name == "volume5":
slowness.append(1/4.0)
elif name == "volume4":
slowness.append(1/4.63)
elif name == "volume3":
slowness.append(1/3.2)
elif name == "volume2":
slowness.append(1/3.9)
elif name == "volume1":
slowness.append(1/4.51)
print("shape\n")
slowness = np.array(slowness)
mgrids = np.array(mgrids)
uniqueTags = np.unique(mgrids)
equiv = np.empty((int(1+uniqueTags.max()),))
nodes = []
for n, tag in enumerate(uniqueTags):
equiv[tag] = n
node = gmsh.model.mesh.getNode(tag)
nodes.append(node[0])

print(mgrids.shape[0])
for n1 in range(mgrids.shape[0]):
for n2 in range(mgrids.shape[1]):
mgrids[n1, n2] = equiv[mgrids[n1, n2]]
nodes = np.array(nodes)
gmsh.finalize()

Source and receiver information

tx,ty,tz,rx,ry,rz,otime=np.genfromtxt("pick.dat",dtype=float,unpack=True,usecols=[1,3,2,4,6,5,7],delimiter=",") #z is lat,x is lon, and y is depth
#receivers
Rx=np.column_stack((rx,ry,rz))

source

Tx=np.column_stack((tx,ty,tz))

print("slowness shape ",np.shape(slowness))
print("nodes shape ",np.shape(nodes))
print("mgrids shape ",np.shape(mgrids))

#Error 1
mesh = Mesh3d(nodes, mgrids.astype(int), method='SPM',n_threads=8,cell_slowness=False,gradient_method=1,maxit=30,n_secondary=2,n_tertiary=2,translate_grid=False)
tt, L = mesh.raytrace(Tx,Rx,slowness=slowness_array,compute_L=True,return_rays=False)

#Error 2
mesh = Mesh3d(nodes, mgrids.astype(int), method='SPM',n_threads=8,cell_slowness=True,gradient_method=1,maxit=30,n_secondary=2,n_tertiary=2,translate_grid=False)
tt, L = mesh.raytrace(Tx,Rx,slowness=slowness_array,compute_L=True,return_rays=False)

[bernard-giroux]

Computing matrix L has not been implemented for tetrahedral meshes with slowness defined inside cells. I think that's something that could be implemented relatively easily, but not in the very short term.

[BhaskarIlla]

Dear Giroux,
Thank you very much for your response. Can you please suggest how to implement it?

[bernard-giroux]

Hi

I have started implementing the needed routines. I have to do some testing before releasing the code, but it should not be too long.

[BhaskarIlla]

Prof. Giroux, thank you very much for your concern and support.

[bernard-giroux]

Hi. It should work now (v1.3.3 uploaded on pipy...). Please let me know how it does.

Best regards

[BhaskarIlla]

Prof. Giroux, thank you very much for your concern. I will check it. Best regards Bhaskar

…

On Fri, Jan 26, 2024, 12:35 AM Bernard Giroux ***@***.***> wrote: Hi. It should work now. Please let me know how it does. Best regards — Reply to this email directly, view it on GitHub <#67 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/A22DKLVRJMH3AXRTNTP2AMTYQKUGLAVCNFSM6AAAAABA7TXQNGVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTSMJQHAYTGOJXHA> . You are receiving this because you authored the thread.Message ID: ***@***.***>