Define boundaries for intersection of two cylinders

variational formulation
1

Hello all,

I am trying to understand how can I define boundary conditions for intersections of two cylinders like below:

hwt
hwt1172×611 39.2 KB

fenics code to generate that:

N = 50                                #Number of divisions
c1 = 0.075                         #radius1 
c2 = 0.212                         #radius2
z1 = 0.06                           #thickness
# Create geometry and generate mesh
cylinder1 = Cylinder(Point(c1, c1, 0.0), Point(c1, c1, z1), 0.075, 0.075)
cylinder2 = Cylinder(Point(c2, c1, 0.0), Point(c2, c1, z1), 0.075, 0.075)
mold = cylinder1 + cylinder2 
mesh = generate_mesh(mold, N)

defining top and bottom looks straight forward as below:

def bottom_bc(x, on_boundary):
    tol = 1E-12   
    return on_boundary and abs(x[2]) < tol

def top_bc(x, on_boundary):
    tol = 1E-12   
    return on_boundary and abs(x[2] - z1) < tol

But I am a bit confused about how can I define sides?

P.S. I am going to assign DirichleteBC for side.

Thank you in advance.

1 Like
2

Please consider the following:

class cyl1(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and np.isclose(np.sqrt((x[0]-c1)**2+(x[1]-c1)**2), 0.075, rtol=3e-2)

class cyl2(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and np.isclose(np.sqrt((x[0]-c2)**2+(x[1]-c1)**2), 0.075, rtol=3e-2)


mf = MeshFunction("size_t", mesh, 2, 0)
cyl1().mark(mf, 2)
cyl2().mark(mf, 3)
File("mf.pvd") << mf

Then you can use the DirichletBC interface to enforce the boundary condition

3

On another note, i tend to recommend using external meshing software such as Gmsh (pygmsh) to do complex geometries. Then you can use meshio to convert the meshing files to XDMF and load them into dolfin. Then you Get a more robust domain marking strategy. There are several threads in the forum describing this procedure.

2 Likes
4

Thanks for for your helpful notes.
because my question is related to this problem I am writing that here:

I was wondering if you can give me a tip about marking a strip on top boundary to impose traction. as follows:

wtr
wtr6768×4477 2.36 MB

the amount of load is constant and I am defining by:

tr = Expression(("0.0", "0.0", "-rho*g"), rho = rho, g = g, degree = 2 + k)
#which is the orange wheel's weight

and I am trying to mark part of red strip while my dynamic load moving forward and backward on top_boundary as below:

A = 0.115                # the amplitude
omega = 2.618            # the frequency
g = Expression('A*sin(omega*t)', A = A,
                 omega = omega, t = 0.0, degree = 2)
class strip(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-12   # tolerance for coordinate comparisons
        return on_boundary and abs(x[2] - z1) < tol and abs(x[1] - 0.051505) > tol and abs(x[1] - 0.098495) < tol and \
        abs(x[0] - g) > tol and abs(x[0] - 1.1*g) < tol
Gamma_T = strip()
Gamma_T.mark(boundary_parts, 1)

So my marking should get update while time passes and compute integral of different parts of the red strip.

But I think the way I am doing this does not make sense, so could you please give me a hint how can I update the place of my loading (orange wheel) while time passing ? (load magnitude is constant)

I was thinking to instead of defining red strip define my loading tr as a user-defined JIT-compiled Expression and then do integration on top boundary instead but I am so confused that how can write ‘cppcode’ for this problem.

Thanks.

5

I have a similar problem, except my mesh is loaded from an external file. Therefor, I don’t know the geometric shape. I want to set a facet marker for the intersection of cut_mesh which overlaps mesh

    vtk_fullpath = "sph_coarse_hollow0001.xdmf"
    mesh = Mesh()
    with XDMFFile(vtk_fullpath) as infile:
            infile.read(mesh)

    vtk_fullpath = "sph_coarse_hollow0001_incube.xdmf"
    cut_mesh = Mesh()
    with XDMFFile(vtk_fullpath) as infile:
            infile.read(cut_mesh)

    boundary_subdomains = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
    boundary_subdomains.set_all(0)

    class SetBound(SubDomain):
        def inside(self, x, on_boundary):
            return on_boundary
    set_bound = SetBound()
    set_bound.mark(boundary_subdomains, 3)
6

ok, found a suitable workflow using trimesh which even uses an stl-file

pip install --user trimesh
sudo apt install python3-rtree

    hold_mesh = trimesh.load("./sph_coarse_hollow0001_incube.stl")
    class intsecHold(SubDomain):
        def inside(self, x, on_boundary):
            return hold_mesh.contains([x])[0] and on_boundary
    set_intsec_hold = intsecHold()
    set_intsec_hold.mark(boundary_subdomains, 1)

If there is a dolfin solution, I’m happy to hear it anyway.