Hey,
I am new to BEMPP, and I am trying to solve a driven source problem. Specifically, I am trying to define a source at a specific edge of the mesh of my geometry.
How would I go about and identifying the edge to flag it and then defining a source there?
Thank you,
Ricardo
Hello,
Once you have defined your grid, you can use the following to print the coordinates of the points at the end of each edge.
for i, e in enumerate(grid.edges.T):
print(i, grid.vertices[:, e[0]], grid.vertices[:, e[1]])
By putting an if statement in the loop, you could use something like this to identify the edges you are interested in.
You can then see which DOF number in your space corresponds to your chosen edge by looking at grid.element_edges and space.local2global. These tell you the three edge numbers for each triangular cell and the three DOF numbers for each triangular cell. (For RWG and SNC spaces, the DOF numbers on each edge are currently always the same as the number of the edge, so these two will contain the same data. I guess there’s no guarantee that this won’t change in future, but I can’t see any reason why it’ll change in the near future.)
You can then initialise a GridFunction using a vector of coefficients with the coefficients set at your chosen index.
Hope this helps. If this isn’t what you’re looking for, do you have a small/simplified example of the kind of thing you’re trying to do?
All the best,
Matthew
Thanks Matthew,
I was able to identify and obtain the edge of interest, but I am having an issue defining the coefficients for the GridFunction.
For reference, I am trying to excite the middle edge of a rectangle sheet of length L, and width W. Once I identify the edge, I tried initializing a GridFunction with a vector of coefficients, where it is a vector of zeros and a specific value at the edge index. I get an error when attempting to solve the LU decomposition:
ValueError: dimension mismatch
Here is the code:
#Constants
frequency = 500E6 #MHZ
vacuum_permittivity = 8.854187817E-12
vacuum_permeability = 4 * np.pi * 1E-7
w = 2 * np.pi * frequency
k = w * np.sqrt(vacuum_permittivity * vacuum_permeability)
wavelength = speed_of_light / frequency
#Rectangle Micro-strip
W = 0.003
L = 1
plate1 = np.array([[-L/2, W/2, 0], [0, W/2, 0], [0, -W/2, 0], [-L/2, -W/2, 0]])
plate2 = np.array([[0, W/2, 0], [0, -W/2, 0], [L/2, -W/2, 0], [L/2, W/2, 0]])
#Edges to Identify
edgeV1 = [0, -W/2, 0]
edgeV2 = [0, W/2, 0]
grid1 = bempp.api.shapes.shapes.screen(plate1, h=(1/10)*wavelength)
grid2 = bempp.api.shapes.shapes.screen(plate2, h=(1/10)*wavelength)
grid = bempp.api.grid.union([grid1, grid2])
#Get Edges
for i, e in enumerate(grid.edges.T):
if ((grid.vertices[:, e[0]] == edgeV1).all() or (grid.vertices[:, e[0]] == edgeV2).all()):
if ((grid.vertices[:, e[1]] == edgeV1).all() or (grid.vertices[:, e[1]] == edgeV2).all()):
edgeInd = i
#Basis Functions to represent problem
div_space = bempp.api.function_space(grid, "RWG", 0)
curl_space = bempp.api.function_space(grid, "SNC", 0)
#Set-up Appropriate Operators - EFIE
elec = bempp.api.operators.boundary.maxwell.electric_field(div_space, div_space, curl_space, k)
#Set-Up Coefficients
coefficients=np.zeros((grid.edges.shape[1],1))
coefficients[edgeInd, 0] = 1
trace_fun = bempp.api.GridFunction(div_space, coefficients= coefficients, dual_space=curl_space)
#Solve System
from bempp.api.linalg import lu
solved_system = lu(elec, trace_fun)
It’s a simple problem, and I thought including the code would help. I suspect that BEMPP is removing the perimeter edges of the rectangle, for continuity purposes, and thus reducing the size of the elec matrix.
In this case, it may be re-numbering the edges. Is there a function that can tell me which edges will be removed, or which will remain? Such that I could re-locate the edge of interest, and initialize a GridFunction with the appropriate sized vector of coefficients.
Best,
-Ricardo
Hello,
It looks like grid.edges.shape[1] isn’t the same as the number of DOFs in your space, I guess due to the non-inclusion of DOFs on boundary edges or something like that.
If you change the setting up of coefficients to the following, it should work. (I also made coefficients a 1D array to slightly simplify things.)
coefficients=np.zeros(div_space.global_dof_count)
coefficients[edgeInd] = 1
trace_fun = bempp.api.GridFunction(div_space, coefficients= coefficients, dual_space=curl_space)
Hi,
I see, this method tells me the DOF, which is 34, of the space to create a proper coefficient vector. But the problem is how would I relate this DOF value in my space to the coefficient index that represents the edge of interest, edgeInd. The previous method, grid.edges.T gives me the global (I think) of the edges, which is from 0-73, my edge of interest in this scenario is 38. I tried relating this using the space.local2global and grid.element_edges, but this gives me a different DOF, 36.
Instead, I went manually and looked at each edge and their neighbors, using grid.edge_neighbors and identified the 34 edges with neighbors and the one that corresponds to the edge of interest. Although this works, this assumes that the edge numeration stays the same, i.e the order of edges is maintained after the boundary ones are removed and they are enumerated as is.
Here is a bit of the code that I used:
count = 0
for ind, n in enumerate(grid.edge_neighbors):
if len(n) > 1:
if (n == localEdgeVals).all():
place = count #Edge of interest
count = count +1 #Coefficient length
coefficients=np.zeros(count, dtype = 'complex')
coefficients[place] = 1
Basically, how do I relate the space.global_dof_count or the DOF in the space to a specific edge? To make sure that the coefficient I am ‘exciting’ pertains to my edge of interest. Perhaps I am using the other functions, space.local2global and grid.element_edges, incorrectly.
Thank you,
-Ricardo
space.local2global tells you the global DOF numbers on each cell of the mesh. So if you know the number of the edge, you can get the DOF on that edge doing something like:
for edges, dofs in zip(grid.element_edges.T, space.local2global):
for e, d in zip(edges, dofs):
if e == edgeInd:
dof = d
break
else:
continue
break
(Sorry, this is a bit complicated as you have to work out which DOF is on each edge by going via the cells.)
In this code, edges and dofs will be arrays of size 3 containing the edge numbers and DOF numbers (respectively) for a cell. e and d then iterate through these arrays, and if the edge number is edgeInd it storted the dof number d and breaks.
So this code tells us what DOF corresponds to a specific edge. But I do not understand why different edges can have the same DOF. Take the following output for example.
The following is some of the output from printing the ‘grid.element_edges.T’ and ‘space.local2global’ :
Edge: [0 1 2] DOF: [2 2 0]
Edge: [3 4 5] DOF: [1 1 1]
Edge: [4 6 7] DOF: [1 0 2]
Edge: [8 2 9] DOF: [3 0 0]
Edge: [3 4 5] all have DOF = 1. Why is this? I am trying to understand how to interpret this information.
Thank you! 
This doesn’t look right. Could you send your mesh so I can look into whether this is an error in your mesh or a bug in Bempp?
Okay. Of course.
Here is a link to the mesh file:
Does this work for you?
The mesh was created using the same code as in the sample I posted earlier.
(Sorry, I could’ve made that mesh myself using your example, but confused this with another example I’ve been looking at that imports a mesh.)
This looks like a bug in the space setup when the grid has edges (eg a screen not a sphere), and I get the same error if I just use a single screen. I’ve made an issue on GitHub (https://github.com/bempp/bempp-cl/issues/72) and I’m currently looking into fixing it
Update: This is a feature not a bug.
The repeated values are for edges that are on the edge of the screen, so no DOF should be on those edges. But instead of putting -1 in the local2global information, 0 is put in the local_multipliers so that the values from these edges are removed when an evaluation is done. This gives a faster implementation than checking that the DOF is not -1.
So, to print the actual DOF numbers for each edge, we should do something like:
for edges, dofs, mults in zip(grid.element_edges.T, space.local2global, space.local_multipliers):
print("Edge:", edges)
print("DOF:", [None if j == 0 else i for i, j in zip(dofs, mults)]
I’ve opened a PR to bempp-cl that adds a convenience function cell_dofs to the space, which would allow the following simpler way to do this:
for i, edges enumerate(grid.element_edges.T):
print("Edge:", edges, "DOF:", space.cell_dofs(i))
I’ll update this reply if/when this functionality is merged
I see, this makes more sense. Thank you!
This now makes clear a different issue I am encountering. Regarding my original problem, I wanted to excite the center edge of my rectangle microstrip. I figured I could force a center edge to exist if I made two separate sheets (left-side and right-side) and united them using grid.union. I noticed that when I looked for my edge and its DOF in the space, it is not assigned a value. Similarly, when I looked at that specific edge in grid.edge_neighbors, the edge does not have two neighbors. This is odd considering it is positioned at the center, with other edges adjacent to it.
Original problem:
Essentially, there seems to be a discontinuity at the center when I use this function grid.union to unite the two rectangle halves. The idea was that the two rectangles would share a common edge, the center edge, and it would be appropriately included when I united the halves. Thus, my question, is there a proper way to unite the meshes to not have this issue? Or is there a way to mesh this microstrip to ensure that I will have an edge in a certain place?
Assuming the position chosen follows the boundary conditions of the RWG functions.
Hi,
the grid.union functions was not meant to join grids that share vertices. All it does is take the two grids, add all vertices from the second grid and modify the element numbering to take into account the new vertex numbers. This works well for grids that are physically separated, but not for grids that are meant to share vertices.
Hello,
So is there a way to properly unite two meshes when they are part of the same object?
Or would I need a different mesher, and then import the mesh?
You’ll need to create the mesh elsewhere and import.
Or you could use the gmsh functions that bempp uses internally to make a mesh, something like this:
from bempp.api.shapes.shapes import __generate_grid_from_geo_string
h = 0.1
W = 0.003
L = 1
stub = f"""
cl = {h};
Point(1) = {{ {-L/2}, {W/2}, 0, cl}};
Point(2) = {{ 0, {W/2}, 0, cl}};
Point(3) = {{ 0, {-W/2}, 0, cl}};
Point(4) = {{ {-L/2}, {-W/2}, 0, cl}};
Point(5) = {{ {L/2}, {-W/2}, 0, cl}};
Point(6) = {{ {L/2}, {W/2}, 0, cl}};
Line(1) = {{1,2}};
Line(2) = {{2,3}};
Line(3) = {{3,4}};
Line(4) = {{4,1}};
Line(5) = {{3,5}};
Line(6) = {{5,6}};
Line(7) = {{6,2}};
Line Loop(1) = {{1, 2, 3, 4}};
Line Loop(2) = {{2,5,6,7}};
Plane Surface(1) = {{1}};
Plane Surface(2) = {{2}};
Physical Surface(1) = {{1,2}};
Mesh.Algorithm = 6;
"""
grid = __generate_grid_from_geo_string(stub)