I did !!!
from dolfin import*
import numpy as np
#mesh = Mesh("My_mesh.xml")
#subdomains = MeshFunction("size_t", mesh, "My_region_region.xml")
Lc=1e-3
p1 = Point(0.0,0.0)
p2 = Point(Lc,Lc)
#my simple mesh for test
mesh_density = 100
mesh = RectangleMesh(p1,p2, mesh_density, mesh_density)
subdomains = MeshFunction("size_t", mesh, 2)
subdomains.set_all(0)
face=BoundaryMesh(mesh, "exterior", order=True)
edge1 = EdgeFunction("size_t", mesh)
edge1.set_all(0)
#we just mark here bc we didnt bring it from xlm
class Alfa(SubDomain):
def inside(self, x, on_boundary):
return (between(x[0], (0., 0.5*Lc)) and between(x[1], (0.5*Lc, Lc)))
alfa = Alfa()
alfa.mark(subdomains,1)
mesh.init()
for v1 in edges(mesh):
if len(v1.entities(2))==2: #non-boundary elements
j=v1.index()
print "The segment",v1.index(), " share with this elements", v1.entities(2)
u1=int(v1.entities(2)[0])
u2=int(v1.entities(2)[1])
if subdomains[u1] != subdomains[u2]:
edge1[j]=1
plot(edge1)
I post it bc if anyone had the same question anyone would find it.
Thanks for all community!