Unable to get the values of a solution and plot a 2d plot in FEniCSx

dolfinx
1

Hi all,
I have solved a linear differential equation d2P/dx2 in fenicsx and got the correct answer too.

# Define physical parameters and boundary conditions
g1 = dolfinx.Constant(mesh, 0)

a = -ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx
f = dolfinx.Constant(mesh, 0)
L = ufl.inner(f, v) * ufl.dx + ufl.inner(g1, v) * ds(3) + ufl.inner(g1, v) * ds(4)

problem = dolfinx.fem.LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()

I am able to plot a contour plot of a solution using pyvista. Plesae find the attached fig.
But when I try to find the values of my solution at nodes using
uh.compute_point_values(). I am not getting the values in ordered manner.

array([[ 0.75],
       [ 0.75],
       [-0.  ],
       [-0.  ],
       [ 1.5 ],
       [ 1.5 ],
       [ 0.75],
       [-0.  ],
       [ 2.25],
       [ 2.25],
       [ 1.5 ],
       [ 0.75],
       [-0.  ],
       [ 3.  ],
       [ 3.  ],
       [ 2.25],
       [ 1.5 ],
       [ 0.75],
       [-0.  ],
       [ 3.  ],
       [ 2.25],
       [ 1.5 ],
       [ 3.  ],
       [ 2.25],
       [ 3.  ]])

If you see the solution in image, there are 25 nodes. The values at all left nodes are 3 and linearly decreasing to zero to the right nodes. I am not understanding why uh.compute_point_values() is not giving values in ordered manner. Why it is starting from 0.75?

Thanks

Pr
Pr1024×768 12.7 KB

2

Please attach a minimal example reproducing the behavior you are observing, as your code is not complete enough to reproduce what you are describing.

I think you are making assumptions on how the coordinates of your mesh is ordered, But its hard to tell from the current post

1 Like
3

I am solving a simple problem in FEniCsx. Here is the .inp mesh file for that.

*Heading
 F:\Intelimek\New folder\rect2.inp
*NODE
1, 0, 0, 0
2, 0.2, 0, 0
3, 0.2, 0.2, 0
4, 0, 0.2, 0
5, 0, 0.15, 0
6, 0, 0.1, 0
7, 0, 0.05, 0
8, 0.05, 0, 0
9, 0.1, 0, 0
10, 0.15, 0, 0
11, 0.2, 0.05, 0
12, 0.2, 0.1, 0
13, 0.2, 0.15, 0
14, 0.15, 0.2, 0
15, 0.1, 0.2, 0
16, 0.05, 0.2, 0
17, 0.05, 0.15, 0
18, 0.05, 0.1, 0
19, 0.05, 0.05, 0
20, 0.1, 0.15, 0
21, 0.1, 0.1, 0
22, 0.1, 0.05, 0
23, 0.15, 0.15, 0
24, 0.15, 0.1, 0
25, 0.15, 0.05, 0
******* E L E M E N T S *************
*ELEMENT, type=T3D2, ELSET=Line1
5, 4, 5
6, 5, 6
7, 6, 7
8, 7, 1
*ELEMENT, type=T3D2, ELSET=Line2
9, 1, 8
10, 8, 9
11, 9, 10
12, 10, 2
*ELEMENT, type=T3D2, ELSET=Line3
13, 2, 11
14, 11, 12
15, 12, 13
16, 13, 3
*ELEMENT, type=T3D2, ELSET=Line4
17, 3, 14
18, 14, 15
19, 15, 16
20, 16, 4
*ELEMENT, type=CPS4, ELSET=Surface1
21, 4, 5, 17, 16
22, 5, 6, 18, 17
23, 6, 7, 19, 18
24, 7, 1, 8, 19
25, 16, 17, 20, 15
26, 17, 18, 21, 20
27, 18, 19, 22, 21
28, 19, 8, 9, 22
29, 15, 20, 23, 14
30, 20, 21, 24, 23
31, 21, 22, 25, 24
32, 22, 9, 10, 25
33, 14, 23, 13, 3
34, 23, 24, 12, 13
35, 24, 25, 11, 12
36, 25, 10, 2, 11
*ELSET,ELSET=inlet
5, 6, 7, 8, 
*ELSET,ELSET=outlet
13, 14, 15, 16, 
*ELSET,ELSET=top
17, 18, 19, 20, 
*ELSET,ELSET=bottom
9, 10, 11, 12, 
*ELSET,ELSET=fluid
21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 
31, 32, 33, 34, 35, 36, 
*NSET,NSET=inlet
1, 4, 5, 6, 7, 
*NSET,NSET=outlet
2, 3, 11, 12, 13, 
*NSET,NSET=top
3, 4, 14, 15, 16, 
*NSET,NSET=bottom
1, 2, 8, 9, 10, 
*NSET,NSET=fluid
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 
11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 
21, 22, 23, 24, 25,

and here is the fenicsx code

# Problem parameters
u_L = 3
u_R = 0
import meshio
msh = meshio.abaqus.read("rect2.inp")
import numpy
def create_mesh2(mesh, cell_type, prune_z=False):
    cells = mesh.get_cells_type(cell_type)
    num_cells = cells.shape[0]
    boundary_tags = numpy.zeros(num_cells, int)
    cell_info = mesh.cell_sets_dict
    
    boundary_tags[cell_info["inlet"]["line"]]=1
    boundary_tags[cell_info["outlet"]["line"]]=2
    boundary_tags[cell_info["top"]["line"]]=3
    boundary_tags[cell_info["bottom"]["line"]]=4
    
    #cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
    #out_mesh = meshio.Mesh(points=mesh.points, cells={cell_type: cells}) #, cell_data={"name_to_read":[cell_data]})
    out_mesh = meshio.Mesh(points=mesh.points, cells={cell_type: cells}, cell_data={"boundaryMarks": [boundary_tags]})
    if prune_z:
        out_mesh.prune_z_0()
    return out_mesh

line_mesh = create_mesh2(msh, "line", prune_z=True)
meshio.write("facet_mesh.xdmf", line_mesh)

triangle_mesh = create_mesh2(msh, "quad", prune_z=True)
meshio.write("mesh.xdmf", triangle_mesh)
import dolfinx.io
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "mesh.xdmf", "r") as xdmf:
    mesh = xdmf.read_mesh(name="Grid")
    ct = xdmf.read_meshtags(mesh, name="Grid")
mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim-1)
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "facet_mesh.xdmf", "r") as xdmf:
    ft = xdmf.read_meshtags(mesh, name="Grid")

V = dolfinx.FunctionSpace(mesh, ("CG", 1))
u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
u_bc_left = dolfinx.Function(V)
left_facets = ft.indices[ft.values==1]
left_dofs = dolfinx.fem.locate_dofs_topological(V, mesh.topology.dim-1, left_facets)
u_bc_right = dolfinx.Function(V)
right_facets = ft.indices[ft.values==2]
right_dofs = dolfinx.fem.locate_dofs_topological(V, mesh.topology.dim-1, right_facets)

x = ufl.SpatialCoordinate(mesh)
with u_bc_left.vector.localForm() as loc:
    loc.set(u_L)
with u_bc_right.vector.localForm() as loc:
    loc.set(u_R)
left_bc = dolfinx.DirichletBC(u_bc_left, left_dofs)
right_bc = dolfinx.DirichletBC(u_bc_right, right_dofs)
bcs = [left_bc, right_bc]

ds = ufl.Measure("ds", domain=mesh, subdomain_data=ft)
g1 = dolfinx.Constant(mesh, 0)

a = -ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx
f = dolfinx.Constant(mesh, 0)
L = ufl.inner(f, v) * ufl.dx + ufl.inner(g1, v) * ds(3) + ufl.inner(g1, v) * ds(4)

# Visualize solution
import pyvista
# Start virtual framebuffer
pyvista.start_xvfb(wait=0.05)

import dolfinx.plot
pyvista_cells, cell_types = dolfinx.plot.create_vtk_topology(mesh, mesh.topology.dim)
grid = pyvista.UnstructuredGrid(pyvista_cells, cell_types, mesh.geometry.x)
point_values = uh.compute_point_values()
if np.iscomplexobj(point_values):
    point_values = point_values.real
grid.point_arrays["u"] = point_values
grid.set_active_scalars("u")

plotter = pyvista.Plotter()
plotter.add_text("uh", position="upper_edge", font_size=14, color="black")
plotter.add_mesh(grid, show_edges=False)
plotter.add_axes()
plotter.show_grid(grid=False)
plotter.view_xy()
plotter.screenshot("Pr.png")

if not pyvista.OFF_SCREEN:
    plotter.show()
else:
    figure = plotter.screenshot("robin_neumann_dirichlet.png")

I am getting the correct contour plot with this code, which I have shown in the previous post.
Now, I want to plot a 2d plot, say a variation of u with respect to the X-axis(basically a line graph). For which I tried to extract u values at nodes using uh.compute_point_values(). But this is giving me values in a random format as given in the above post.
What should I do, so that I will get the u values in ordered format.
I want in this manner,
x=0, u=3
x=0.05, u=2.25
x=0.1, u=1.5
x=0.15, u=0.75
x=0.2, u=0

Thanks

4

Please have a look at: Implementation — FEniCSx tutorial

1 Like
5

I found the solution in this tutorial.
Thank you very much.