Hi,
I am trying to define a source at a specific vertice of the mesh grid. How can I identify a specific point (x,y,z) to set a specific coefficient (e.g., 1.0+1j*0.0)?
bempp.api.GridFunction(div_space, coefficients=coefficients, dual_space=curl_space)
This is the code I have until now:
# Source driven - The DOFs of RWG and SNC spaces are at the midpoints of the edges of each cell.
# Point (vertice) to identify
pointV1 = [-0.25*OperationWavelengthPA, 0.0*OperationWavelengthPA, 0.0*OperationWavelengthPA]
#Get points
pointInd=[] # A list because maybe more than one point is found
for i in np.arange(0,grid.number_of_edges,1):
e=grid.edges[0,i]
V1axis=np.abs(grid.vertices[:, e]-pointV1)
DiffPoint1=np.sqrt(V1axis[0]**2+V1axis[1]**2+V1axis[2]**2)
e=grid.edges[1,i]
V2axis=np.abs(grid.vertices[:, e]-pointV1)
DiffPoint2=np.sqrt(V2axis[0]**2+V2axis[1]**2+V2axis[2]**2)
if (DiffPoint1<1.0*GridSizeh and DiffPoint2<1.0*GridSizeh):
pointInd.append(i)
print('pointInd: '+str(pointInd))
# Get DOF
dof=[]
"""
dof=pointInd
for iIter in range(0,len(pointInd),1):
for edges, dofs in zip(grid.edges, div_space.local2global):
for e, d in zip(edges, dofs):
if e == pointInd[iIter]:
dof[iIter] = d
"""
#Set-Up Coefficients
for iIter in range(0,len(pointInd),1):
for edges, dofs in zip(grid.element_edges.T, div_space.local2global):#, div_space.local_multipliers):
for i, e in enumerate(edges):
if (e == pointInd[iIter]):
dof.append(dofs[i])
print('dof: '+str(dof))
Improved code
# Study how to identify DOFs for coefficient excitation
pointV1 = [0.0*OperationWavelengthPA, 0.25*OperationWavelengthPA, 0.0*OperationWavelengthPA]
# DOF: Degree of freedom; are the points in the middle of the grid triangular-mesh-edges which the excitation coefficients are set.
# In particular, might be of interst to set coefficient for the three DOFs in the grid cells of a RWC/SNC space
# Point to identify the closest edge. The closest edge is the one that has the two vertices closest to the point of interest
pointV1aux=pointV1
gridAux=grid
div_spaceAux=div_space
GridSizehAux=GridSizeh
def BemppIdentifyDOFcellsFromPoint(pointV1aux,gridAux,div_spaceAux,GridSizehAux,CelldofAux=False):
## Seems better for coefficient excitation to excite a single edge, and correspondingly its DOF. Variable CelldofAux=False identifies the closests edge DOF. Variable CelldofAux=True identifies the cells DOFs
#Get edges
EdgeInd=[] # Number identifying the edge number. A list because maybe more than one point is found
v1EdgeInd=[] # List of vertice 1
v2EdgeInd=[] # List of vertice 2
Diff1EdgeInd=[] # List of distances of vertice 1
Diff2EdgeInd=[] # List of distances of vertice 2
for i in np.arange(0,gridAux.number_of_edges,1): # First identify a few candidate edges
v1=gridAux.edges[0,i] # vertice 1 of the edge
V1axis=np.abs(gridAux.vertices[:, v1]-pointV1aux)
DiffPoint1=np.sqrt(V1axis[0]**2+V1axis[1]**2+V1axis[2]**2)
v2=gridAux.edges[1,i] # vertice 2 of the edge
V2axis=np.abs(gridAux.vertices[:, v2]-pointV1aux)
DiffPoint2=np.sqrt(V2axis[0]**2+V2axis[1]**2+V2axis[2]**2)
#if (DiffPoint1<1.5*GridSizehAux and DiffPoint2<1.5*GridSizehAux):
EdgeInd.append(i)
Diff1EdgeInd.append(DiffPoint1)
Diff2EdgeInd.append(DiffPoint2)
v1EdgeInd.append(gridAux.vertices[:, v1])
v2EdgeInd.append(gridAux.vertices[:, v2])
# Order from minimum distance (the sum of the two distances between the two vertices of the edge) to maximum distance the identified edges
EdgeInd=np.array(EdgeInd)# Convert to numpy array
v1EdgeInd=np.array(v1EdgeInd)# Convert to numpy array
v2EdgeInd=np.array(v2EdgeInd)# Convert to numpy array
Diff1EdgeInd=np.array(Diff1EdgeInd)# Convert to numpy array
Diff2EdgeInd=np.array(Diff2EdgeInd)# Convert to numpy array
TotalDiffEdgeInd=Diff1EdgeInd+Diff2EdgeInd # Total distance of the edge vertices to the point
OrderIndicesEdgeInd=np.argsort(TotalDiffEdgeInd)
OrderedMinEdgeIndPot=EdgeInd[OrderIndicesEdgeInd]
OrderedMinv1EdgeIndPot=v1EdgeInd[OrderIndicesEdgeInd]
OrderedMinv2EdgeIndPot=v2EdgeInd[OrderIndicesEdgeInd]
OrderedMinDiff1EdgeIndPot=Diff1EdgeInd[OrderIndicesEdgeInd]
OrderedMinDiff2EdgeIndPot=Diff2EdgeInd[OrderIndicesEdgeInd]
# Consider the first two edges as the rulers. Now, order the other edges that share vertices
# Find the common vertices of the first two edges, and then find the opposite vertices that the third edge should have
OrderedMinEdgeInd=np.zeros_like(OrderedMinEdgeIndPot)
OrderedMinv1EdgeInd=np.zeros_like(OrderedMinv1EdgeIndPot)
OrderedMinv2EdgeInd=np.zeros_like(OrderedMinv2EdgeIndPot)
OrderedMinDiff1EdgeInd=np.zeros_like(OrderedMinDiff1EdgeIndPot)
OrderedMinDiff2EdgeInd=np.zeros_like(OrderedMinDiff2EdgeIndPot)
if (CelldofAux):
OrderedMinEdgeInd[0:2]=OrderedMinEdgeIndPot[0:2]
OrderedMinv1EdgeInd[:,0:2]=OrderedMinv1EdgeIndPot[:,0:2]
OrderedMinv2EdgeInd[:,0:2]=OrderedMinv2EdgeIndPot[:,0:2]
OrderedMinDiff1EdgeInd[0:2]=OrderedMinDiff1EdgeIndPot[0:2]
OrderedMinDiff2EdgeInd[0:2]=OrderedMinDiff2EdgeIndPot[0:2]
v1a=gridAux.edges[0,OrderedMinEdgeInd[0]] # vertice a of edge 1
v1b=gridAux.edges[1,OrderedMinEdgeInd[0]] # vertice b of edge 1
v2a=gridAux.edges[0,OrderedMinEdgeInd[1]] # vertice a of edge 2
v2b=gridAux.edges[1,OrderedMinEdgeInd[1]] # vertice b of edge 2
if (v1a==v2a):
v3a=v1b
v3b=v2b
elif(v1a==v2b):
v3a=v1b
v3b=v2a
elif(v1b==v2a):
v3a=v1a
v3b=v2b
else: # v1b==v2b
v3a=v1a
v3b=v2a
# Find the third edge
Edge3Ind=0
for iIter in range(2,len(OrderedMinEdgeInd),1):
if (((v3a==gridAux.edges[0,OrderedMinEdgeIndPot[iIter]]).all() and (v3b==gridAux.edges[1,OrderedMinEdgeIndPot[iIter]]).all() ) or ((v3b==gridAux.edges[0,OrderedMinEdgeIndPot[iIter]]).all() and (v3a==gridAux.edges[1,OrderedMinEdgeIndPot[iIter]]).all() )):
Edge3Ind=iIter
OrderedMinEdgeInd[2]=OrderedMinEdgeIndPot[Edge3Ind]
OrderedMinv1EdgeInd[2,:]=OrderedMinv1EdgeIndPot[Edge3Ind,:]
OrderedMinv2EdgeInd[2,:]=OrderedMinv2EdgeIndPot[Edge3Ind,:]
OrderedMinDiff1EdgeInd[2]=OrderedMinDiff1EdgeIndPot[Edge3Ind]
OrderedMinDiff2EdgeInd[2]=OrderedMinDiff2EdgeIndPot[Edge3Ind]
# Find a fourth edge, following the min distance criterion
if (Edge3Ind==2):
OrderedMinEdgeInd[3]=OrderedMinEdgeIndPot[3]
OrderedMinv1EdgeInd[3,:]=OrderedMinv1EdgeIndPot[3,:]
OrderedMinv2EdgeInd[3,:]=OrderedMinv2EdgeIndPot[3,:]
OrderedMinDiff1EdgeInd[3]=OrderedMinDiff1EdgeIndPot[3]
OrderedMinDiff2EdgeInd[3]=OrderedMinDiff2EdgeIndPot[3]
else:
OrderedMinEdgeInd[3]=OrderedMinEdgeIndPot[2]
OrderedMinv1EdgeInd[3,:]=OrderedMinv1EdgeIndPot[2,:]
OrderedMinv2EdgeInd[3,:]=OrderedMinv2EdgeIndPot[2,:]
OrderedMinDiff1EdgeInd[3]=OrderedMinDiff1EdgeIndPot[2]
OrderedMinDiff2EdgeInd[3]=OrderedMinDiff2EdgeIndPot[2]
#print('OrderedMinEdgeInd: '+str(OrderedMinEdgeInd))
# Should take the first three edges, that normally define a cell. In can happen that the point of inteest falls in and edge or vertice (then this situation should be handled)
# But take the first six edges, because usually at the grid border less dofs are identified
NpotentialEdges=4
EdgeIndCell=OrderedMinEdgeInd[0:NpotentialEdges]
v1IndCell=OrderedMinv1EdgeInd[0:NpotentialEdges]
v2IndCell=OrderedMinv2EdgeInd[0:NpotentialEdges]
EdgeDiff1IndCell=OrderedMinDiff1EdgeInd[0:NpotentialEdges]
EdgeDiff2IndCell=OrderedMinDiff2EdgeInd[0:NpotentialEdges]
# Checks
if (np.sum(EdgeDiff1IndCell==0.0)>0 or np.sum(EdgeDiff2IndCell==0.0)>0): # Hit vertice
print("Vertice hit. It should be dealt a side!")
DotProdArray1=np.zeros((3),dtype=np.float32)
DotProdArray2=np.zeros((3),dtype=np.float32)
AbsProdArray1=np.zeros((3),dtype=np.float32)
AbsProdArray2=np.zeros((3),dtype=np.float32)
for iIter in range(0,3,1):
DotProdArray1[iIter]=np.abs(np.dot(pointV1aux-v1IndCell[iIter],v2IndCell[iIter]-v1IndCell[iIter]))
DotProdArray2[iIter]=np.abs(np.dot(pointV1aux-v2IndCell[iIter],v2IndCell[iIter]-v1IndCell[iIter]))
AbsProdArray1[iIter]=np.sqrt((pointV1aux[0]-v1IndCell[iIter,0])**2+(pointV1aux[1]-v1IndCell[iIter,1])**2+(pointV1aux[2]-v1IndCell[iIter,2])**2)*np.sqrt((v2IndCell[iIter,0]-v1IndCell[iIter,0])**2+(v2IndCell[iIter,1]-v1IndCell[iIter,1])**2+(v2IndCell[iIter,2]-v1IndCell[iIter,2])**2)
AbsProdArray2[iIter]=np.sqrt((pointV1aux[0]-v2IndCell[iIter,0])**2+(pointV1aux[1]-v2IndCell[iIter,1])**2+(pointV1aux[2]-v2IndCell[iIter,2])**2)*np.sqrt((v2IndCell[iIter,0]-v1IndCell[iIter,0])**2+(v2IndCell[iIter,1]-v1IndCell[iIter,1])**2+(v2IndCell[iIter,2]-v1IndCell[iIter,2])**2)
# Probably now hitting an edge is well dealt
#if (np.sum(DotProdArray1==AbsProdArray1)>0 or np.sum(DotProdArray2==AbsProdArray2)>0): # Hit edge
# print("Edge hit. It should be dealt a side!")
# Get DOF
else:
EdgeIndCell=OrderedMinEdgeIndPot[0:1]
v1IndCell=OrderedMinv1EdgeIndPot[0:1]
v2IndCell=OrderedMinv2EdgeIndPot[0:1]
EdgeDiff1IndCell=OrderedMinDiff1EdgeIndPot[0:1]
EdgeDiff2IndCell=OrderedMinDiff2EdgeIndPot[0:1]
Potentialdof=[]
"""
dofAux=pointInd
for iIter in range(0,len(EdgeIndCell),1):
for edges, dofs in zip(gridAux.edges, div_spaceAux.local2global):
for e, d in zip(edges, dofs):
if e == EdgeIndCell[iIter]:
dof[iIter] = d
"""
#Set-Up Coefficients
for iIter in range(0,len(EdgeIndCell),1):
for edges, dofs in zip(gridAux.element_edges.T, div_spaceAux.local2global):#, div_spaceAux.local_multipliers):
for i, e in enumerate(edges):
if (e == EdgeIndCell[iIter]):
Potentialdof.append(dofs[i])
if (CelldofAux): # DOFs forming a cell
dofAux = ([*set(Potentialdof)])[0:3] # Remove duplicate DOFs from the list, and consider the first three dofs (to form a cell)
else: # DOF of the closest edge
dofAux = [([*set(Potentialdof)])[0]] # Remove duplicate DOFs from the list, and consider only the first dofs (of an edge)
print('dof: '+str(dofAux))
return dofAux
dof=BemppIdentifyDOFcellsFromPoint(pointV1,grid,div_space,GridSizeh,False)