underworldcode / UWGeodynamics

Underworld Geodynamics
Other
81 stars 32 forks source link

How could I extract an interface in a 3D model to plot a relative topography picture in UW2 ? #257

Closed CEShale closed 2 years ago

CEShale commented 2 years ago

..............................................

Model = GEO.Model(elementRes=(40, 20, 50), minCoord=(0. u.kilometer, 0. u.kilometer, 0 u.kilometer), maxCoord=(40. u.kilometer, 20. u.kilometer, 5.0 u.kilometer), gravity=(0., 0., -9.81 * u.meter / u.second ** 2))

..............................................

npoints = 1000

xx = np.linspace(GEO.nd(Model.minCoord[0]), GEO.nd(Model.maxCoord[0]), npoints) yy = np.linspace(GEO.nd(Model.minCoord[1]), GEO.nd(Model.maxCoord[1]), npoints) zz = GEO.nd(sand1.top)

xx,yy,zz = np.meshgrid(xx,yy,zz) coords = np.ndarray((xx.size, 3))

coords[:,0] = xx.ravel() coords[:,1] = yy.ravel() coords[:,2] = zz.ravel()

coords[:,2] =coords[:,2]-np.mean(coords[:, 2]))

interface1 = Model.add_passive_tracers(name="Interface1", vertices=coords)

np.savetxt("interface1.txt",interface1.data)

mesh = uw.mesh.FeMesh_Cartesian( elementType = ("Q1/dQ0"), elementRes = (npoints, npoints), minCoord = (0., 0.), maxCoord = (40.0, 20.0))

RelativeTopographyField = Model.mesh.add_variable( nodeDofCount=1 )

for index, coord in enumerate(mesh.data): for i in range(0,npoints-1): if coords[i,0]==(coord[0]-GEO.nd(Model.minCoord[0]))/(GEO.nd(Model.maxCoord[0])- GEO.nd(Model.minCoord[0]))npoints and coords[i,1]==(coord[1]- GEO.nd(Model.minCoord[1]))/(GEO.nd(Model.maxCoord[1])- GEO.nd(Model.minCoord[1]))npoints: RelativeTopographyField.data[index] =float(coords[i,2]-np.mean(coords[:, 2]))

fig = vis.Figure(figsize=(800,400),title="Relative Topography") fig.append( vis.objects.Mesh(mesh) ) fig.append( vis.objects.Surface( mesh, RelativeTopographyField, colours="blue white red" ) ) fig.show()

rbeucher commented 2 years ago

Hi @CEShale,

Can you clarify what you are trying to do?

Romain

rbeucher commented 2 years ago

I am closing this. Feel free to reopen if you need. Please clarify what you are trying to achieve.