In the below MWE, I have an annular mesh (gmsh script provided below). Here is a screenshot of the mesh:
I need to find the dofs of the interior nodes only. Here is a summary of what I tried:
(a) Defined a Lagrange vector function space V of order 1 on the whole mesh.
(b) Defined vertex_to_dof_map(V) and dof_to_vertex_map(V).
(c) Reshaped the vertex_to_dof_map(V) to a matrix such that each row corresponds to the dofs of a coordinate of a node on the mesh.
(d) Found the coordinates of the whole mesh.
(e) Found the coordinates of the boundaries.
(f) Deleted the boundary coordinates from the coordinates of the whole mesh to get only the interior coordinates.
Now what I want to do is scan through the rows of the vertex_to_dof_map matrix and check if a row of dofs corresponds to the coordinates of a node on the boundary. If so, I need to delete that row so that I end up with a dof matrix that corresponds to the interior nodes only. Now I am not sure how to execute this or if there is a smarter way to do this.
MWE:
import numpy
def create_mesh(mesh, cell_type, prune_z=False):
cells = mesh.get_cells_type(cell_type)
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]})
return out_mesh
#***************************************************************************************************************************************************************************************
#**************************************************************IMPORTING GEOMETRY FROM GMSH*********************************************************************************************
msh = meshio.read("MWE.msh")
triangle_mesh = create_mesh(msh, "triangle", True)
line_mesh = create_mesh(msh, "line", True)
meshio.write("mesh.xdmf", triangle_mesh)
meshio.write("mf.xdmf", line_mesh)
from dolfin import *
mesh = Mesh()
xdmf = XDMFFile(mesh.mpi_comm(),"mesh.xdmf")
xdmf.read(mesh)
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim())
with XDMFFile("mesh.xdmf") as infile:
infile.read(mvc, "name_to_read")
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
xdmf.close()
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim()-1)
with XDMFFile("mf.xdmf") as infile:
infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
ds_custom = Measure("dS", domain=mesh, subdomain_data=mf) # Defining global boundary measure, for global surface measure, plug in "dS" instead of "ds", which is used for subdomain boundary measures
dx_custom = Measure("dx", domain=mesh, subdomain_data=cf)
#################################################################
V = VectorFunctionSpace(mesh, 'P', 1) # Defining vector function space on the whole mesh
W = FunctionSpace(mesh, 'P', 3)
all_coords = mesh.coordinates() # Finding coordinates of all the nodes including the boundary nodes
bndy = BoundaryMesh(mesh,"exterior") # Finding boundary mesh
bndy_coords = bndy.coordinates() # Finding coordinates of the boundary nodes
v2d = np.array(vertex_to_dof_map(V), dtype = int) # Defining vertex to dof map on the whole mesh
d2v = np.array(dof_to_vertex_map(V), dtype = int) # Defining dof to vertex map on the whole mesh
#*********** interior media coordinates ******************
indices = [] # Creating an empty list to store indices of the coordinates from the whole mesh that matches the boundary coordinates
indexlist = np.array(indices) # Converting the list to an array
for i in range(0,len(all_coords)): # Looping through all the nodes on the mesh
for j in range(0,len(bndy_coords)): # Looping through all the nodes on the boundary
if mesh.coordinates()[i][0]==bndy.coordinates()[j][0] and mesh.coordinates()[i][1]==bndy.coordinates()[j][1]: # Checking if the coordinates on the whole mesh matches the coordinates on the boundary
indexlist = np.append(indexlist,i) # If there is a match, storing the index of the coordinates the match
intcoords=np.delete(all_coords, [indexlist[:]], 0) # From all the coordinates of the whole mesh, deleting the coordinates of the boundary through the index list
print(bndy_coords,len(bndy_coords)) # Making sure if all add up to the total number of nodes
print(intcoords,len(intcoords))
print(all_coords,len(all_coords))
###########################################################
########## FINDING DOFS OF ALL NODES ############
v2d_matrix_all = v2d.reshape(len(all_coords), 3) # Creating a matrix of vertex_to_dof map such that each row of the matrix corresponds to the dof of a coordinate on the whole mesh
print(v2d_matrix_all,v2d_matrix_all.shape)
####################################################
# Now I want to scan each row of v2d_matrix_all and check if that row of dof corresponds to a node on the boundary. If so, I will delete that row. Doing so will give me the dofs of the interior nodes only. I am not sure how to do this.
GMSH scipt:
// Gmsh project created on Thu May 19 12:23:52 2022
SetFactory("OpenCASCADE");
//+
Circle(1) = {0, 0, 0, 1, 0, 2*Pi};
//+
Circle(2) = {0, 0, 0, 2, 0, 2*Pi};
//+
Curve Loop(1) = {2};
//+
Curve Loop(2) = {1};
//+
Plane Surface(1) = {1, 2};
//+
Physical Curve("C1", 1) = {1};
//+
Physical Curve("C1", 1) += {1};
//+
Physical Curve("C2", 2) = {2};
//+
Physical Surface("S1", 3) = {1};
Any help or suggestion would be greatly appreciated!