Output values only from x-axis as array

Hello,

I have a question about the tutorial. Here is the link:

https://jorgensd.github.io/dolfinx-tutorial/chapter1/nitsche.html

Is it possible to output only the values from the x-axis as an array from the plot? So only the values of the lower edge? If so, how?

I mean from uh in the example.

Have you considered: Implementation — FEniCSx tutorial

The function bb_tree = geometry.BoundingBoxTree(domain, domain.topology.dim) does not work. There is no error message, but the programme does not run.

Is there an alternative to determine the values on the x-axis?

Edit: Sorry, the function does work. Thank you!

It doesn’t work after all, but I don’t understand what the error is. Here are the two pictures. Does anyone understand what is wrong?

membrane_rank0

How do I find out which points collides with the bb_tree?

So these points from the example:

tol = 0.001 # Avoid hitting the outside of the domain
y = np.linspace(-1 + tol, 1 - tol, 101)
points = np.zeros((3, 101))
points[1] = y

MWE:

#Mesh Erzeugung:(mit GMSH)

from mpi4py import MPI
from dolfinx import mesh
import warnings

warnings.filterwarnings("ignore")

import gmsh

gmsh.initialize()

from dolfinx import fem
from dolfinx import nls
from dolfinx import log
from dolfinx import io
from dolfinx import geometry
from dolfinx.cpp.mesh import to_type, cell_entity_type
from dolfinx.fem import Constant, Function, FunctionSpace, assemble_scalar, dirichletbc, form, locate_dofs_topological, set_bc
from dolfinx.fem.petsc import apply_lifting, assemble_matrix, assemble_vector, create_vector, set_bc
from dolfinx.graph import create_adjacencylist
from dolfinx.geometry import BoundingBoxTree, compute_collisions, compute_colliding_cells
from dolfinx.io import (XDMFFile, cell_perm_gmsh, distribute_entity_data, extract_gmsh_geometry,
                        extract_gmsh_topology_and_markers, ufl_mesh_from_gmsh)
from dolfinx.mesh import create_mesh, meshtags_from_entities
import ufl
from ufl import (FacetNormal, FiniteElement, Identity, Measure, TestFunction, TrialFunction, VectorElement,
                 as_vector, div, dot, ds, dx, inner, lhs, grad, nabla_grad, rhs, sym)
import numpy as np
import matplotlib.pyplot as plt
import tqdm.notebook
from mpi4py import MPI
from petsc4py import PETSc
import scipy.constants as c


#Domain-Definieren:
L = 1.
H = 1.
c_x =  0.50
c_y = 0.50
r = 0.01
gdim = 2
rank = MPI.COMM_WORLD.rank

if rank == 0:
    rectangle = gmsh.model.occ.addRectangle(0,0,0, L, H, tag=1)
    obstacle = gmsh.model.occ.addDisk(c_x, c_y, 0, r, r)

if rank == 0:
    Mesh_domain = gmsh.model.occ.cut([(gdim, rectangle)], [(gdim, obstacle)])
    gmsh.model.occ.synchronize()


#Verbinden aller Mesh-Punkte
Mesh_domain_marker = 1

if rank == 0:
    volumes = gmsh.model.getEntities(dim=gdim)
    gmsh.model.addPhysicalGroup(volumes[0][0], [volumes[0][1]], Mesh_domain_marker)
    gmsh.model.setPhysicalName(volumes[0][0], Mesh_domain_marker, "Mesh_domain")

inlet_marker, outlet_marker, wall_marker, obstacle_marker, wall_top = 2, 3, 4, 5, 6
inflow, outflow, walls, obstacle, top = [], [], [], [], []

if rank == 0: #Walls sind Ground, Obstical ist Leiter
    boundaries = gmsh.model.getBoundary(volumes, oriented=False)
    for boundary in boundaries:
        center_of_mass = gmsh.model.occ.getCenterOfMass(boundary[0], boundary[1])
        if np.allclose(center_of_mass, [0, H/2, 0]):
            inflow.append(boundary[1])
        elif np.allclose(center_of_mass, [L, H/2, 0]):
            outflow.append(boundary[1])
        elif  np.allclose(center_of_mass, [L/2, 0, 0]):
            walls.append(boundary[1])
        elif np.allclose(center_of_mass, [L/2, H, 0]):
            top.append(boundary[1])
        else:
            obstacle.append(boundary[1])
            
    gmsh.model.addPhysicalGroup(1, walls, wall_marker)
    gmsh.model.setPhysicalName(1, wall_marker, "Walls")
    gmsh.model.addPhysicalGroup(1, inflow, inlet_marker)
    gmsh.model.setPhysicalName(1, inlet_marker, "Inlet")
    gmsh.model.addPhysicalGroup(1, outflow, outlet_marker)
    gmsh.model.setPhysicalName(1, outlet_marker, "Outlet")
    gmsh.model.addPhysicalGroup(1, obstacle, obstacle_marker)
    gmsh.model.setPhysicalName(1, obstacle_marker, "Obstacle")
    gmsh.model.addPhysicalGroup(1, top, wall_top)
    gmsh.model.setPhysicalName(1, wall_top, "top")

if rank == 0:
    gmsh.option.setNumber("Mesh.RecombinationAlgorithm", 8)
    gmsh.option.setNumber("Mesh.RecombineAll", 2)
    gmsh.option.setNumber("Mesh.SubdivisionAlgorithm", 1)
    gmsh.model.mesh.generate(gdim)
    gmsh.model.mesh.optimize("Netgen") 

if MPI.COMM_WORLD.rank == 0:

    # Mesh Geometrie   
    x = extract_gmsh_geometry(gmsh.model)

    # Topologie fĂĽr jedes Element vom Mesh
    topologies = extract_gmsh_topology_and_markers(gmsh.model)
    num_cell_types = len(topologies.keys())
    cell_information = {}
    cell_dimensions = np.zeros(num_cell_types, dtype=np.int32)
    
    for i, element in enumerate(topologies.keys()):
        properties = gmsh.model.mesh.getElementProperties(element)
        name, dim, order, num_nodes, local_coords, _ = properties
        cell_information[i] = {"id": element, "dim": dim, "num_nodes": num_nodes}
        cell_dimensions[i] = dim

    # Sort elements by ascending dimension
    perm_sort = np.argsort(cell_dimensions)

    # Broadcast cell type data and geometric dimension
    cell_id = cell_information[perm_sort[-1]]["id"]
    tdim = cell_information[perm_sort[-1]]["dim"]
    num_nodes = cell_information[perm_sort[-1]]["num_nodes"]
    cell_id, num_nodes = MPI.COMM_WORLD.bcast([cell_id, num_nodes], root=0)

    if tdim - 1 in cell_dimensions:
        num_facet_nodes = MPI.COMM_WORLD.bcast( cell_information[perm_sort[-2]]["num_nodes"], root=0)
        gmsh_facet_id = cell_information[perm_sort[-2]]["id"]
        marked_facets = np.asarray(topologies[gmsh_facet_id]["topology"], dtype=np.int64)
        facet_values = np.asarray(topologies[gmsh_facet_id]["cell_data"], dtype=np.int32)

    cells = np.asarray(topologies[cell_id]["topology"], dtype=np.int64)
    cell_values = np.asarray(topologies[cell_id]["cell_data"], dtype=np.int32)

else:
    cell_id, num_nodes = MPI.COMM_WORLD.bcast([None, None], root=0)
    cells, x = np.empty([0, num_nodes], np.int64), np.empty([0, gdim])
    cell_values = np.empty((0,), dtype=np.int32)
    num_facet_nodes = MPI.COMM_WORLD.bcast(None, root=0)
    marked_facets = np.empty((0, num_facet_nodes), dtype=np.int64)
    facet_values = np.empty((0,), dtype=np.int32)

    
# Create distributed mesh
ufl_domain = ufl_mesh_from_gmsh(cell_id, gdim)
gmsh_cell_perm = cell_perm_gmsh(to_type(str(ufl_domain.ufl_cell())), num_nodes)
cells = cells[:, gmsh_cell_perm]
mesh = create_mesh(MPI.COMM_WORLD, cells, x[:, :gdim], ufl_domain)
tdim = mesh.topology.dim
fdim = tdim - 1

# Permute facets from MSH to DOLFINx ordering
facet_type = cell_entity_type(to_type(str(ufl_domain.ufl_cell())), fdim, 0)
gmsh_facet_perm = cell_perm_gmsh(facet_type, num_facet_nodes)
marked_facets = np.asarray(marked_facets[:, gmsh_facet_perm], dtype=np.int64)

local_entities, local_values = distribute_entity_data(mesh, fdim, marked_facets, facet_values)
mesh.topology.create_connectivity(fdim, tdim)
adj = create_adjacencylist(local_entities)
bb_tree = geometry.BoundingBoxTree(mesh, mesh.topology.dim)

# Create DOLFINx MeshTags
ft = meshtags_from_entities(mesh, fdim, adj, np.int32(local_values))
ft.name = "Facet tags"


#Funktionspace und Vektorfunktionspace defenieren
s_cg1 = FiniteElement("CG", mesh.ufl_cell(), 1)
V = FunctionSpace(mesh, s_cg1)
V_E= fem.VectorFunctionSpace(mesh, ("DG", 0))
fdim = mesh.topology.dim - 1

#Konstanten, Funktionen und Randbedingungen definieren:
phi_0 = 80000 #V
h = c_y

u1_0 = Constant(mesh, PETSc.ScalarType(phi_0))  

#Boundary conditions definieren
# Walls
u1_nonslip = np.array((0,) *mesh.geometry.dim, dtype=PETSc.ScalarType)
u11_nonslip = np.array((0,) *mesh.geometry.dim, dtype=PETSc.ScalarType)
wall_facets = ft.indices[ft.values == wall_marker]
bcu1_walls = dirichletbc(PETSc.ScalarType(0), locate_dofs_topological(V, fdim, wall_facets), V)

# Obstacle
obstacle_facets = ft.indices[ft.values == obstacle_marker]
bcu1_obstacle = dirichletbc(PETSc.ScalarType(8e4), locate_dofs_topological(V, fdim, obstacle_facets), V)

#Boundary conditions for equations: 
bcu1 = [bcu1_obstacle, bcu1_walls] # Spannung am Leiter und Ränder

#Laplace
u1 = ufl.TrialFunction(V) 
v1 = ufl.TestFunction(V)
k = Constant(mesh, PETSc.ScalarType(0))
f = k/c.epsilon_0
epsilon = c.epsilon_0
a = ufl.inner(ufl.grad(u1),ufl.grad(v1)) * ufl.dx 
L = ufl.inner(f , v1) * ufl.dx
u1 = fem.Function(V)
problem1 = fem.petsc.LinearProblem(a, L, bcu1, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
phi1 = problem1.solve()

with io.VTKFile(mesh.comm, "Laplace.pvd", "w") as vtk:
    vtk.write([phi1._cpp_object])

with io.XDMFFile(mesh.comm, "Laplace.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_function(phi1)

###https://jorgensd.github.io/dolfinx-tutorial/chapter1/membrane_code.html#making-curve-plots-throughout-the-domain
tol = 0.001
y = np.linspace(-1 + tol, 1 - tol, 101)
points = np.zeros((3, 101))
points[1] = y
u_values = []
p_values = []
cells = []
points_on_proc = []
# Find cells whose bounding-box collide with the the points
cell_candidates = geometry.compute_collisions(bb_tree, points.T)
# Choose one of the cells that contains the point
colliding_cells = geometry.compute_colliding_cells(mesh, cell_candidates, points.T)
for i, point in enumerate(points.T):
    if len(colliding_cells.links(i))>0:
        points_on_proc.append(point)
        cells.append(colliding_cells.links(i)[0])
points_on_proc = np.array(points_on_proc, dtype=np.float64)
u_values = phi1.eval(points_on_proc, cells)

import matplotlib.pyplot as plt
fig = plt.figure()
plt.plot(points_on_proc[:,1], u_values, "k", linewidth=2, label="Phi")
plt.grid(True)
plt.xlabel("x")
plt.legend()
# If run in parallel as a python file, we save a plot per processor
plt.savefig(f"membrane_rank{MPI.COMM_WORLD.rank:d}.png")

print("finish")

It works in version 0.4.1

As your code is far from minimal I would suggest the following:

  1. print points_on_proc
  2. print cells
  3. print u_values

You can also use the plot over line filter and paraview and compare the results

Here are the outputs:

  1. print points_on_proc
[[0.00000000e+00 1.11022302e-16 0.00000000e+00]
 [0.00000000e+00 1.99800000e-02 0.00000000e+00]
 [0.00000000e+00 3.99600000e-02 0.00000000e+00]
 [0.00000000e+00 5.99400000e-02 0.00000000e+00]
 [0.00000000e+00 7.99200000e-02 0.00000000e+00]
 [0.00000000e+00 9.99000000e-02 0.00000000e+00]
 [0.00000000e+00 1.19880000e-01 0.00000000e+00]
 [0.00000000e+00 1.39860000e-01 0.00000000e+00]
 [0.00000000e+00 1.59840000e-01 0.00000000e+00]
 [0.00000000e+00 1.79820000e-01 0.00000000e+00]
 [0.00000000e+00 1.99800000e-01 0.00000000e+00]
 [0.00000000e+00 2.19780000e-01 0.00000000e+00]
 [0.00000000e+00 2.39760000e-01 0.00000000e+00]
 [0.00000000e+00 2.59740000e-01 0.00000000e+00]
 [0.00000000e+00 2.79720000e-01 0.00000000e+00]
 [0.00000000e+00 2.99700000e-01 0.00000000e+00]
 [0.00000000e+00 3.19680000e-01 0.00000000e+00]
 [0.00000000e+00 3.39660000e-01 0.00000000e+00]
 [0.00000000e+00 3.59640000e-01 0.00000000e+00]
 [0.00000000e+00 3.79620000e-01 0.00000000e+00]
 [0.00000000e+00 3.99600000e-01 0.00000000e+00]
 [0.00000000e+00 4.19580000e-01 0.00000000e+00]
 [0.00000000e+00 4.39560000e-01 0.00000000e+00]
 [0.00000000e+00 4.59540000e-01 0.00000000e+00]
 [0.00000000e+00 4.79520000e-01 0.00000000e+00]
 [0.00000000e+00 4.99500000e-01 0.00000000e+00]
 [0.00000000e+00 5.19480000e-01 0.00000000e+00]
 [0.00000000e+00 5.39460000e-01 0.00000000e+00]
 [0.00000000e+00 5.59440000e-01 0.00000000e+00]
 [0.00000000e+00 5.79420000e-01 0.00000000e+00]
 [0.00000000e+00 5.99400000e-01 0.00000000e+00]
 [0.00000000e+00 6.19380000e-01 0.00000000e+00]
 [0.00000000e+00 6.39360000e-01 0.00000000e+00]
 [0.00000000e+00 6.59340000e-01 0.00000000e+00]
 [0.00000000e+00 6.79320000e-01 0.00000000e+00]
 [0.00000000e+00 6.99300000e-01 0.00000000e+00]
 [0.00000000e+00 7.19280000e-01 0.00000000e+00]
 [0.00000000e+00 7.39260000e-01 0.00000000e+00]
 [0.00000000e+00 7.59240000e-01 0.00000000e+00]
 [0.00000000e+00 7.79220000e-01 0.00000000e+00]
 [0.00000000e+00 7.99200000e-01 0.00000000e+00]
 [0.00000000e+00 8.19180000e-01 0.00000000e+00]
 [0.00000000e+00 8.39160000e-01 0.00000000e+00]
 [0.00000000e+00 8.59140000e-01 0.00000000e+00]
 [0.00000000e+00 8.79120000e-01 0.00000000e+00]
 [0.00000000e+00 8.99100000e-01 0.00000000e+00]
 [0.00000000e+00 9.19080000e-01 0.00000000e+00]
 [0.00000000e+00 9.39060000e-01 0.00000000e+00]
 [0.00000000e+00 9.59040000e-01 0.00000000e+00]
 [0.00000000e+00 9.79020000e-01 0.00000000e+00]
 [0.00000000e+00 9.99000000e-01 0.00000000e+00]]
  1. print cells
[16, 16, 16, 16, 29, 29, 29, 34, 34, 34, 30, 30, 30, 17, 17, 17, 18, 18, 18, 35, 35, 35, 36, 36, 36, 36, 37, 37, 37, 38, 38, 38, 19, 19, 19, 10, 10, 10, 6, 6, 6, 3, 3, 3, 1, 1, 1, 0, 0, 0, 0]
  1. print u_values
[[8.90387774e-12]
 [1.55501111e+03]
 [3.11002221e+03]
 [4.66503332e+03]
 [6.20445734e+03]
 [7.74159072e+03]
 [9.27872410e+03]
 [1.07864333e+04]
 [1.22840044e+04]
 [1.37815755e+04]
 [1.52410842e+04]
 [1.66768271e+04]
 [1.81125699e+04]
 [1.95023828e+04]
 [2.08439078e+04]
 [2.21854329e+04]
 [2.34855331e+04]
 [2.47117842e+04]
 [2.59380354e+04]
 [2.71352296e+04]
 [2.82358191e+04]
 [2.93364085e+04]
 [3.04174349e+04]
 [3.13282816e+04]
 [3.22391282e+04]
 [3.31499748e+04]
 [3.39269739e+04]
 [3.47005375e+04]
 [3.54741010e+04]
 [3.60847091e+04]
 [3.66658465e+04]
 [3.72469839e+04]
 [3.77378229e+04]
 [3.81933224e+04]
 [3.86488218e+04]
 [3.90313941e+04]
 [3.93634119e+04]
 [3.96954296e+04]
 [3.99797921e+04]
 [4.02087631e+04]
 [4.04377341e+04]
 [4.06398807e+04]
 [4.07886194e+04]
 [4.09373582e+04]
 [4.10723884e+04]
 [4.11546472e+04]
 [4.12369061e+04]
 [4.13147633e+04]
 [4.13406473e+04]
 [4.13665313e+04]
 [4.13924153e+04]]

I can’t do anything with this. How can I use plot-over-line filter and paraview to compare the results?

As you are already plotting along the line above, you can use the filter in paraview to do the same plot.

I’m not sure why you claim that the plot you showed above is wrong. The plot runs from 0 to 40 000, which seems sensible if you look at the plot of the contour.

However, I only wanted to display the values of the x-axis in the diagram and not the entire plot. As I understood the bb_tree, it should be responsible for this, or?

Isn’t that what you are doing in

The boundingbox tree is one of multiple steps used to plot over an axis.
The bounding box tree gives you a list of candidate cells for each 3D point whose bounding box collides with the point.
Then

Finds cells actually colliding with the points.

This filters out all points not found on the current process

This finds the value u for each point in the corresponding cell.

plots the data over they axis.

My plot is only in 2D.
What I don’t understand is that in the plot there are only values <20,000 on the x-axis. How can the value in the diagram increase to 40,000? And shouldn’t the plot, for example, look more like this?

A

You’ve applied a Dirichlet bc of u=0 on the x-axis, so I wouldn’t expect the solution to look like this.

With this code, you are plotting points along the y-axis (x=0). Using a ParaView “Plot Over Line” filter, I confirmed that the result you are getting with matplotlib is correct for the y-axis. Do you intend to plot over the line y=1? As you can see below in the ParaView plot, u has a similar shape on y=1 to what you state that you are expecting.

To plot over y=1, you should use

tol = 0.001
y = np.linspace(-1 + tol, 1 - tol, 101)
points = np.zeros((3, 101))
points[0] = y
points[1] = np.ones_like(y)
1 Like

Thank you, that’s exactly what I want. I only have the problem that my diagram looks different with the settings. It is only a straight line. I suspect that I only map the y value on the x axis in the diagram and not the x value, because the straight line is still at 1. Could someone else tell me why and fix this is?
And just for information, how would I have to change the code in order to get y=0?

membrane_rank0

tol = 0.001
y = np.linspace(-1 + tol, 1 - tol, 101)
points = np.zeros((3, 101))
points[0] = y
points[1] = np.ones_like(y)
u_values = []
cells = []
points_on_proc = []
cell_candidates = geometry.compute_collisions(bb_tree, points.T)
colliding_cells = geometry.compute_colliding_cells(mesh, cell_candidates, points.T)
for i, point in enumerate(points.T):
    if len(colliding_cells.links(i))>0:
        points_on_proc.append(point)
        cells.append(colliding_cells.links(i)[0])
points_on_proc = np.array(points_on_proc, dtype=np.float64)
u_values = phi1.eval(points_on_proc, cells)

import matplotlib.pyplot as plt
fig = plt.figure()
plt.plot(points_on_proc[:,1], u_values,"k", linewidth=2, label="Phi")
plt.grid(True)
plt.legend()

The x-component of points is points[0], and the y-component is points[1]. To extract the solution at y=0, simply remove the following line of code that sets the y-coordinate to 1:

The x-component of points_on_proc is points_on_proc[:, 0], and the y-component is points_on_proc[:, 1]. Therefore, the following code is plotting y on the abscissa.

To plot x on the abscissa, use

plt.plot(points_on_proc[:,0], u_values,"k", linewidth=2, label="Phi")
1 Like

Thank you very much that has helped me a lot

1 Like