extra indices of fault_blocks

[Ali1990dashti]

Dear community, Thanks for making such usefull package. When I try to export the indices of my fault block, I see that the length of array is more than the length of my grid and also in high resolution cases of a model with only one fault, I see several values for my indices rather than only having 0. and 1.. I copied a simple example here. I expect to have a fault_block with length of 10*10*10 but for me it is not (it is 1008). And if I increase the resolution, np.unique(sol.fault_block) will not be [0., 1.]. This is my syntax (I am using gempy 2.2.9):

import gempy as gp
import numpy as np
geo_model=gp.create_model('model')
geo_model = gp.init_data(geo_model, extent=[0,100,0,30,0,100], resolution=[10,10,10])
geo_model.set_default_surfaces()
geo_model.add_surface_points(X=10, Y=15, Z=70, surface='surface1')
geo_model.add_surface_points(X=20, Y=15, Z=70, surface='surface1')
geo_model.add_surface_points(X=80, Y=15, Z=45, surface='surface1')
geo_model.add_surface_points(X=90, Y=15, Z=45, surface='surface1')
geo_model.add_orientations(X=10, Y=15, Z=70, surface='surface1', orientation=np.array([0,0,1]))
geo_model.add_features('Fault1')
geo_model.add_surfaces('fault')
geo_model.add_surface_points(X=63.36, Y=15, Z=80.43, surface='fault')
geo_model.add_surface_points(X=40, Y=15, Z=40, surface='fault')
geo_model.add_orientations(55, Y=15, Z=50, surface='fault', orientation=np.array([270,64.5,1]))
gp.map_stack_to_surfaces(geo_model, {'Fault1':'fault', 'Layers':('surface1','surface2')})
geo_model.set_is_fault('Fault1')
gp.set_interpolator(geo_model, output=['geology'], theano_optimizer='fast_compile')
sol=gp.compute_model(geo_model)
print (np.unique(sol.fault_block))
print (sol.fault_block.shape)
print (sol.grid.values.shape)

In advance, I do appreciate your consideration.

[javoha]

Hello @Ali1990dashti ,

I am actually not entirely sure where the additional 8 entries come from - but I checked it and you can just use sol.fault_block[0,:1000] to get the correct entries (the remaining 8 or not necessary). Alternatively you can directly get the right matrix using (only in your specific case): sol.block_matrix[0].

Hope this helps for now, I will try to see why these additional values are there.

[Ali1990dashti]

Dear @javoha , Thanks for your feedack. I am also using the same trick (sol.fault_block[0,:1000]) but I was not sure whether the remainings (after 1000.) are redundant or not. I do appreciate if you also check np.unique(sol.fault_block). Just increase the resolution then you will see non 0. and 1. values. Thanks again for being that much supportive.

[Japhiolite]

Dear @Ali1990dashti I'm currently looking through open issues. The extra points in the fault block are still there in the recent gempy version. They're connected to the number of input points given in a model. As for the values changing from 0 and 1 towards numbers in between (up to 1.5ish), I assume some drift. I'll check more into that.

However, best for the fault blocks would be to use sol.block_matrix[0] as Jan suggested.

Cheers the other Jan

[Ali1990dashti]

Dear @Japhiolite ,

Thanks for the update. I very appreciate your time.

Cheers Ali