SubDomain

As said in my previous post, a missing feature in fem-fenics was the marking of subdomains. Indeed, I proposed an example that needed a file generated with a run of the corresponding Python code, which is not, honestly, the best approach. In order to address this issue, these days I have implemented a new class, subdomain, which can be used to mark mesh entities. In the following I will describe how to use this new functionality. Here is the code:

pkg load fem-fenics msh

ufl start Subdomains

ufl fe = FiniteElement "(""CG"", triangle, 2)"

ufl u = TrialFunction (fe)

ufl v = TestFunction (fe)

ufl

ufl a0 = Coefficient (fe)

ufl a1 = Coefficient (fe)

ufl g_L = Coefficient (fe)

ufl g_R = Coefficient (fe)

ufl f = Coefficient (fe)

ufl

ufl a = "inner(a0*grad(u), grad(v))*dx(0) + inner(a1*grad(u), grad(v))*dx(1)"

ufl L = g_L*v*ds(1) + g_R*v*ds(3) + f*v*dx(0) + f*v*dx(1)

ufl end

# Create mesh and define function space

x = y = linspace (0, 1, 65);

[msh, facets] = Mesh (msh2m_structured_mesh (x, y, 0, 4:-1:1));

V = FunctionSpace ("Subdomains", msh);

# Define boundary conditions

bc1 = DirichletBC (V, @(x, y) 5.0, facets, 2);

bc2 = DirichletBC (V, @(x, y) 0.0, facets, 4);

# Define problem coefficients

a0 = Constant ("a0", 1.0);

a1 = Constant ("a1", 0.01);

g_L = Expression ("g_L", @(x, y) - 10*exp(- (y - 0.5) ^ 2));

g_R = Constant ("g_R", 1.0);

f = Constant ("f", 1.0);

# Define subdomains - Here are the edits #

domains = MeshFunction ("dx", msh, 2, 0);

obstacle = SubDomain (@(x,y) (y >= 0.5) && (y <= 0.7) && ...

                      (x >= 0.2) && (x <= 1.0), false);

domains = mark (obstacle, domains, 1);

# Define variational form

a = BilinearForm ("Subdomains", V, V, a0, a1, domains);

L = LinearForm ("Subdomains", V, g_L, g_R, f, facets, domains);

# Assemble system

[A, b] = assemble_system (a, L, bc1, bc2);

sol = A \ b;

u = Function ("u", V, sol);

# Save solution in VTK format

save (u, "subdomains");

# Plot solution

[X, Y] = meshgrid (x, y);

U = u (X, Y);

surf (X, Y, U);

As you can see, it is basically the same as in the previous post, except the line used to import the meshfunction. I wrote in the corresponding comment where the edits are to be found. Now the workflow comprises these steps: first of all, a meshfunction needs to be created, then a subdomain, in the end we should mark cells. 

The call to MeshFunction is something new, since it is now possible to instantiate a meshfunction given a mesh, the required topological dimension and the value to initialise it with. Moreover, the optional label "dx" means that it can be used in calls to BilinearForm or LinearForm to supply markers for subsets of the integration domains. In the example, this instruction returns a meshfunction of dimension 2, which means it holds values associated with each triangle in the mesh, initialised to be 0 in every entry.

The subsequent instruction, instead, defines a subdomain, passing as arguments a function handle and a logical flag. The former will be the override of the dolfin::SubDomain::inside method, so it must return true for entities contained in the subset and false otherwise. In facts it checks whether the coordinates are inside the 2-interval defining the obstacle. The latter, instead, can be used to ask for a boundary subdomain, when set to true.

At last, mark is called to set the entries corresponding to cells inside the subdomain to 1, so that the returned meshfunction now represents the obstacle: after these lines, the variable named domains assumes value 1 on cells inside the obstacle region and 0 outside. Thus, it is now possible to solve a problem whose formulation entails subdomains entirely using fem-fenics.