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][4*n:(4*n+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 https://github.com/groupeLIAMG/ttcr/issues/67#issuecomment-1910813978, or unsubscribe https://github.com/notifications/unsubscribe-auth/A22DKLVRJMH3AXRTNTP2AMTYQKUGLAVCNFSM6AAAAABA7TXQNGVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTSMJQHAYTGOJXHA . You are receiving this because you authored the thread.Message ID: @.***>