What is the best approach to mesh refinement on boundary?

Question

I'm solving Laplace equation with discontinuous Dirichlet conditions. So it would be great if there were more cells near the points of breaks than in other areas. What is the best approach to do it (both to UnitSquareMesh and mesh built using generate_mesh)?

from dolfin import *
from mshr import *

nx, ny = 30, 30
plx, prx = 0.2, 0.8
u0, u1 = 0.0, 0.05

#mesh = UnitSquareMesh(nx, ny)
mesh = generate_mesh(Rectangle(Point(0, 0), Point(1, 1)), 50)

class Plate(SubDomain):
    def inside(self, x, on_boundary):
        return (plx < x[0] < prx) and near(x[1], 1.0)

class Boundary(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0.0) or near(x[0], 1.0) or near(x[1], 0.0) or (near(x[1], 1.0) and not (plx < x[0] < prx))

# Initialize sub-domain instances
plate = Plate()
boudary = Boundary()


# Initialize mesh function for boundary domains
boundaries = FacetFunction("size_t", mesh)
boundaries.set_all(0)
plate.mark(boundaries, 1)
boudary.mark(boundaries, 2)

# Define function space and basis functions
V = FunctionSpace(mesh, "Lagrange", 1)
u = TrialFunction(V)
v = TestFunction(V)

# Define Dirichlet boundary conditions
bcs = [DirichletBC(V, Constant(u0), boundaries, 2),
       DirichletBC(V, Constant(u1), boundaries, 1)]

a = inner(grad(u), grad(v)) * dx
L = Constant(0.0) * v * dx

u = Function(V)

problem = LinearVariationalProblem(a, L, u, bcs)
solver = LinearVariationalSolver(problem)
solver.solve()

plot(u)

E = -project(grad(u), VectorFunctionSpace(mesh, "Lagrange", 1))

plot(E)

plot(mesh)

interactive()

Comments

This might not be an answer to your question (I never used the built-in mesh capabilities) but I know it is easy to do in GMSH.

---

Thank for the tip!
I would be very grateful if you sent me on joraillusion@yandex.by some example of usage gmsh with fenics if you have.

Answer 1

Hi,
also you can try something like this:

from dolfin import *

mesh = UnitSquareMesh(20, 20, "crossed")

# Break point
p   = Point(0.0, 0.5)
tol = 0.05

# Selecting edges to refine
class Border(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], p.x(), tol) and near(x[1], p.y(), tol) and on_boundary

Border = Border()

# Number of refinements
nor = 3

for i in range(nor):
    edge_markers = EdgeFunction("bool", mesh)
    Border.mark(edge_markers, True)

    adapt(mesh, edge_markers)
    mesh = mesh.child()

plot(mesh, interactive=True)

Comments

Perfect! Thanks!

Answer 2

I made a small example that might help you and hopefully also others.

• Create a folder and add two files: unitsquare.geo and gmsh_example.py. You can find their contents below.

• Provided GMSH is installed, run the following in the terminal to create the mesh:

gmsh -2 unitsquare.geo

• A new file unitsquare.msh appears. Convert this to a FEniCS friendly mesh format running:

dolfin-convert unitsquare.msh unitsquare.xml

• Run the Python file gmsh_example.py. I used your example and adjusted it a little. You can also call the terminal command withing your Python script, if you like to have everything automatized.

That should get you started. Let me know if it works and good luck!

unitsquare.geo

// GMSH allows you to define variables
l = 0.1;
lsmall = 0.01;
// Points are defined by 3 coordinates
// The 4-th entry describes the required
// element size around this point
Point(1) = {0, 0, 0, l};
Point(2) = {0, 1, 0, l};
Point(3) = {1, 1, 0, l};
Point(4) = {1, 0, 0, l};
Point(5) = {0.2, 1, 0, lsmall};
Point(6) = {0.8, 1, 0, lsmall};
Point(7) = {0.1, 1, 0, l};
Point(8) = {0.9, 1, 0, l};
Point(9) = {-0, 0.7, 0, l};
Point(10) = {1, 0.7, 0, l};
Line(1) = {2, 7};
Line(2) = {7, 5};
Line(3) = {5, 6};
Line(4) = {6, 8};
Line(5) = {8, 3};
Line(6) = {3, 10};
Line(7) = {10, 4};
Line(8) = {4, 1};
Line(9) = {1, 9};
Line(10) = {9, 2};
Line(11) = {9, 10};
Line Loop(12) = {1, 2, 3, 4, 5, 6, -11, 10};
Plane Surface(13) = {12};
Line Loop(14) = {11, 7, 8, 9};
Plane Surface(15) = {14};
// These Physical Lines will be used to define
// the boundaries; if you start counting from 1
// then all internal boundaries will be labeled 0
Physical Line(1) = {10, 1, 9, 8, 7, 6, 5};
Physical Line(2) = {3};
Physical Line(3) = {2, 4};
// Physical Surfaces define regions
// I usually start counting from 0
// You need at least 2 or dolfin-convert ignores these
Physical Surface(0) = {13};
Physical Surface(1) = {15};

gmsh_example.py

from dolfin import *

# Defining mesh, boundaries and regions
# This assumes some preprocessing (see explanation)
mesh = Mesh('unitsquare.xml')
boundaries = MeshFunction('size_t', mesh, 'unitsquare_facet_region.xml')
regions = MeshFunction('size_t', mesh, 'unitsquare_physical_region.xml')
dx = Measure('dx')[regions]
ds = Measure('ds')[boundaries]

# Define function space and basis functions
V = FunctionSpace(mesh, "Lagrange", 1)
u = TrialFunction(V)
v = TestFunction(V)

# Define Dirichlet boundary conditions
# The boundary labeled 3 can be uncommented to
# demand a discontinuous jump on the boundary
u0, u1 = 0.0, 0.05
bcs = [DirichletBC(V, Constant(u0), boundaries, 1),
       # DirichletBC(V, Constant(u0), boundaries, 3),
       DirichletBC(V, Constant(u1), boundaries, 2)]

# Define weak form and solve
# Although not used here, dx(i) and ds(i) could be used
# to specify an integral over a specific region/boundary
a = inner(grad(u), grad(v)) * dx
L = Constant(0.0) * v * dx
u = Function(V)
problem = LinearVariationalProblem(a, L, u, bcs)
solver = LinearVariationalSolver(problem)
solver.solve()
E = -project(grad(u), VectorFunctionSpace(mesh, "Lagrange", 1))

# Plotting
plot(mesh)
plot(boundaries)
plot(regions)
plot(u)
plot(E)
interactive()

Comments

For clarity, I created the unitsquare.geo file using GMSH's graphical interface. You only need to remember to change the numbering of the physical lines and physical regions if you do so.

---

Thank you for a detailed answer!