Problem when assigning values to layers in a mesh type mesh

[GabrielMoreno09]

Hi guys, I have doubts about how to assign values to the layers when I have already built my model in a tensor mesh.
I did it in a very usual way, here you have the example of how I did it.
If you have any other ideas on how to do that procedure I would be very grateful.
Thanks

Pos: If you need the notebook just just ask for it.

dx = 35
ncx = 40
dy = 37
ncy = 15
dz = 200
ncz = 40
dpadh = 3.5
exp_h = 1.5

hx = [(dx, dpadh, -exp_h), (dx, ncx), (dx, dpadh, exp_h)]
hy = [(dy, dpadh, -exp_h), (dy, ncy), (dy, dpadh, exp_h)]
hz = [(dz, 3, -exp_h), (dz, ncz)]
meshes = TensorMesh([hx, hy, hz], x0=[853614,987543,-6000])
cell_centers = meshes.cell_centers
xy_ravel=cell_centers

meshes.plot_grid(color="midnightblue", linewidth=0.1)
plt.plot(x_topo, y_topo, z_topo, '.', color='tomato')

array_sol = geo_model.solutions.raw_arrays.lith_block.reshape((50,50,50), order='F')

lith_block = gp.compute_model_at(geo_model, cell_centers)

Definir los valores de contraste de densidad para cada unidad en g/cm³

background_density = 0.5
densities = {
0: background_density, # Densidad del complejo metamórfico (g/cm³) #2.9
1: -0.2, # Densidad de Complejo_Dacitico #2.2
2: 0, # Densidad de Cuerpo_Intrusivo #2.4
3: 0.1, # Densidad de Flujo_Piroclastico #2.5
4: 0.2, # Densidad de Sucesion_Piroclastica #2.6
5: 0.3 # Densidad de Sucesion_Volcanica #2.7
}

Encontrar los índices para las celdas de malla activas (por ejemplo, celdas debajo de la superficie)

ind_active = active_from_xyz(meshes, topo_xyz)

Definir el mapeo del modelo a las celdas activas

nC = int(ind_active.sum())
model_map = maps.IdentityMap(nP=nC)

Definir el modelo. Los modelos en SimPEG son arrays vectoriales.

model = np.zeros(nC)

Mapear los valores de contraste de densidad del modelo GemPy al modelo SimPEG

for i, density in densities.items():
ind_lith = lith_block[ind_active] == i
model[ind_lith] = density - background_density

[javoha]

Hello @GabrielMoreno09,
thanks for your question. Just to make sure I understand correctly: You basically are trying to use gempy with a custom grid (defined by your variable cell_centers and then assign density values per lithology?
I think the notebook would help so that we can take a closer look.
Cheers
Jan

[GabrielMoreno09]

You are right, I am trying to use gempy with a custom grid from simPEG, to assign density values to each of the lithologies.
Here are the link to the repository to you can check it.
https://github.com/GabrielMoreno09/Modelamiento_Inversion_VCM_Project/blob/master/Modelamiento_directo/Modelamiento_SimPEG%20(4).ipynb

I would like to know if there is a better way to do this.

[javoha]

Hi,
to be honest I think this looks like a smart way to do it. There may be a way to do it without the loop but I think your approach is perfectly reasonable (and should not be very expensive computationally). Is there anything specific about it that you dislike?