extracting solution on a given plane

Main Forum Old-Forum - do not post here
1

Dear NGSolve Developers,

I am trying to extract the solution in a 3D problem on a given plane. So, once the solution is available from the solver, I set up a cloud of points defined over a 2D region [-R,R]^2:

comm.Barrier()
if (rank==0):
    Nx=100;
    Ny=100;
    xmin,xmax,ymin,ymax = [-R,R,-R,R]; # R is given
    plot_grid = np.mgrid[xmin:xmax:Nx*1j,ymin:ymax:Ny*1j];
    points = np.vstack((plot_grid[0].ravel(),plot_grid[1].ravel(),
    np.array([0]*plot_grid[0].size)));
    fem_xx = my_points[0,:]
    fem_xy = my_points[1,:]
    for ipt in range(0,len(fem_xx)):
        val = u(mesh(fem_xx[ipt],fem_xy[ipt],0.0))

While this works without mpirun, when i launch the script with more than 1 cores, i end up in SIGSEGV fault. Most likely because inside the for loop the partitioning information is lost. Is there a better way of extracting the solution on a given plane which works with mpirun? Other possibility which i could think of is to write the vtk file and then use vtkcutter to slice through the domain and save the solution.
Thanks in advance!

2

With MPI, every proc only has one part of the mesh, with rank 0 not keeping any part at all.

The code

    for ipt in range(0,len(fem_xx)):
        val = u(mesh(fem_xx[ipt],fem_xy[ipt],0.0))

will be called on every rank. Every rank will only have elements for SOME of the points in your plane.
Currently it segfaults for points that are not in the local part of the mesh.

I have attached a small example where I first check every point for whether it is in the mesh or not:

for ipt in range(0,len(fem_xx)):
    if mesh.Contains(fem_xx[ipt],fem_xy[ipt],0.0):
        mesh_pnt = mesh(fem_xx[ipt],fem_xy[ipt],0.0)
        val = gfu(mesh_pnt)
        pvs.append((mesh_pnt, val))

However, in that case you still have your points and values distributed among the MPI procs. If you want them all in one place, mpi4py might be what you are looking for.

Depending on what you want to do with the restriction of your solution, constructing your mesh such that it resolves the plane might be another solution.

Otherwise the VTK approach is also a way.

Attachment: mesh_eval.py