This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

Mark segments between two different materials.

0 votes

Hi everyone

I have mesh and there are two different materials market in this mesh. It is 2D. I d like to mark the edges that share two elements with different marked domains or properties.

mesh = Mesh("Square.xml") 
subdomains = MeshFunction("size_t", mesh, "Square_physical_region.xml")

I have the subdomains marked 0 or 1.

I d like to do a looping between elements around each element in all mesh. And if there are one element around that has different properties (or subdomain marked) the segment between this two elements must be marked.

Or one looping in segments if this segment share 2 elements that are in different subdomains marked... this segment must be marked.

how can I mark segments between two different materials?

Kind Regards,
Leo

asked Feb 6, 2016 by LeoCosta FEniCS User (1,190 points)

1 Answer

+1 vote

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!

answered Feb 7, 2016 by LeoCosta FEniCS User (1,190 points)
edited Feb 7, 2016 by LeoCosta
...