How to mark boundaries to a circle and change solution from cartesian to polar?

1

Hello! I have a solution on the boundary of a circle in Cartesian coordinates. I want to change it to polar coordinates. But before that I need to solve one more issue to create a minimal working example with a built-in mesh. Since I always import my mesh from gmsh with physical tags attached, I am having trouble defining the markers to my domain which is simply a circle of radius 3 created through the Python script. Any help would be appreciated.

from mshr import *
from ufl import *
from dolfin import * 
import matplotlib.pyplot as plt

from dolfin import *
import meshio
from dolfin import Mesh, XDMFFile, File, MeshValueCollection, cpp

# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, \
               "eliminate_zeros": True, \
               "precompute_basis_const": True, \
               "precompute_ip_const": True}

mesh = generate_mesh(Circle(Point(0,0),3), 50) # Generating mesh for the domain
bmesh = BoundaryMesh(mesh,"exterior") # Generating mesh for the boundary

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

bndry=bmesh()
bndry.mark(boundary_markers, 1)
  
V = FunctionSpace(mesh, "Lagrange", 1)

# Define boundary condition
G , mu = 1, 0.1
u_D=Constant(0.0)

bc = DirichletBC(V, u_D, 1)


# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(-G/mu) # f=-G/mu
a = dot(grad(u), grad(v))*dx
L = -f*v*dx

# Compute solution
usol = Function(V)
solve(a == L, usol, bc)

The error reads:

bndry=bmesh()
TypeError: 'dolfin.cpp.mesh.BoundaryMesh' object is not callable
2

you can just import the mesh from gmsh and then in Fenics can create boundary markers. See one of my previous post (Occlude a cylinder with suction pressure), where I first imported a 3D cylinder mesh and then defined the boundary marker in Fenics.

3

I already have the solution to this problem in Cartesian coordinate form by importing the mesh and its physical tags from gmsh. But for people to reproduce the code, a built-in mesh is the preferred option. So I was trying to create a built-in mesh.

But the main question is that how to convert the solution (usol) to a polar reference frame?

4

In case anyone is interested, here is how to extract subdomains and meshes from structures created and tagged in GMSH with the tags intact:

Bulk = MeshFunction("size_t" , s5 , cf)  #saves the interior info of the mesh
Boundary = MeshFunction("size_t", s5 , mf)
Sub_Mesh=SubMesh(mesh,cf,5)
surface_marker = MeshFunction("size_t", Sub_Mesh, Sub_Mesh.topology().dim() - 1, 0)
boundary_marker = MeshFunction("size_t", Sub_Mesh, Sub_Mesh.topology().dim()-1, 0)
ncells = MeshFunction("size_t", Sub_Mesh, Sub_Mesh.topology().dim())
volume_marker = MeshFunction("size_t", Sub_Mesh, Sub_Mesh.topology().dim())
surface_marker = MeshFunction("size_t", Sub_Mesh, Sub_Mesh.topology().dim(), 0)
vmap = Sub_Mesh.data().array("parent_vertex_indices", 0)
cmap = Sub_Mesh.data().array("parent_cell_indices", Sub_Mesh.topology().dim())

for c in cells(Sub_Mesh):
  parent_cell = Cell(mesh, cmap[c.index()])
  volume_marker.array()[c.index()] = cf.array()[parent_cell.index()]
  surface_marker.array()[c.index()] = cf.array()[parent_cell.index()]
  for f in facets(parent_cell):
    for g in facets(c):
      g_vertices = vmap[g.entities(0)]
      if set(f.entities(0)) == set(g_vertices):
        surface_marker.array()[g.index()] = mf.array()[f.index()]
        boundary_marker.array()[g.index()] = mf.array()[f.index()]
    
ds = Measure('ds', subdomain_data=boundary_marker)
dx = Measure('dx', subdomain_data=surface_marker)
5

Hi, I have a question related to defining a hollow cylinder made of two different materials. To define a hollow cylinder I use an outer tube and subtract the inner one. Then, however, I want to assign a different material to part of the outline, say from -pi/4 to pi/4. This seems easiest to me in cylindrical coordinates, while the cylinder is initially defined in Cartesian coordinates. Is there a way to assign a different material to a part of the cylinder in this way?

from fenics import *
from mshr import *

R = 0.1
L = 1
curvature = -0.1

outer_part = Cylinder(Point(0, 0, L), Point(0, 0, 0), R, R)
inner_part = Cylinder(Point(0, 0, L), Point(0, 0, 0), R/2, R/2)
hole = Cylinder(Point(0, curvature*L, 3*L/4), Point(0, curvature*L, L/4), R, R)
mesh_hole = generate_mesh(outer_part - inner_part - hole, 20)
plot(mesh_hole)

I try to assign the material by defining two different subdomains. Is there a difference here between the 3D situation and the 2D situation?

hole2 = Cylinder(Point(0, curvature*L, 3*L/4), Point(0, curvature*L, L/4), R/2, R/2)
domain = (outer_part- inner_part)
domain2 = (hole - hole2)
domain.set_subdomain(1, domain)
domain.set_subdomain(2, domain2)

In that case I get the following error in the 3D situation, so I assume that there is not yet availability in 3D. Is there perhaps another way to do this in that case?;

RuntimeError                              Traceback (most recent call last)
Cell In[65], line 17
     15 plot(mesh_hole)
     16 mesh = generate_mesh(outer_part - inner_part, 20)
---> 17 domain.set_subdomain(1, domain)
     18 domain.set_subdomain(2, domain2)

RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to setting subdomain.
*** Reason:  Subdomains are currently supported only in 2D.
*** Where:   This error was encountered inside CSGGeometry.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.2.0.dev0
*** Git changeset:  ubuntu

Thank you in advance

6

You should not rely on mshr, as it has been deprecated for about 4 years. I would strongly suggest moving to external mesh generators with the capability of marking subdomains, such as gmsh.

1 Like