### Implementation of subdomains in the 2017.2.0 version changed (compared to 1.6.0) ?

379
views
-1
6 months ago by
Hello there,

I have thought about this quite a while and wasn't able to find a solution or a similar post here.
I have updated my Fenics/Dolfin version to the latest available one (2017.2.0) from 1.6.0 and found some different behaviour in my code concerning the loading and working with meshes and subdomains in 3D from gmsh and using dirichlet boundary conditions on them. But since the same problem occurs in 2D I use this as an example. I used the following code to generate the meshes in gmsh:

//Defining the circle points
Point(1) = {0, 0, 0, .01};
Point(2) = {-.1, 0, 0, .01};
Point(3) = {.1, 0, 0, .01};

//Defining the box points
Point(4) = {-1, -1, 0, .1};
Point(5) = {1, -1, 0, .1};
Point(6) = {1, 1, 0, .1};
Point(7) = {-1, 1, 0, .1};

//Defining the lines
Circle(1) = {2, 1, 3};
Circle(2) = {3, 1, 2};
Line(3) = {7, 4};
Line(4) = {4, 5};
Line(5) = {5, 6};
Line(6) = {6, 7};

// Defining the surfaces
Line Loop(7) = {3, 4, 5, 6};
Line Loop(8) = {2, 1};
Plane Surface(9) = {7, 8};
Plane Surface(10) = {8};

// Physical labels
Physical Line(11) = {2, 1};
Physical Line(12) = {3};
Physical Line(13) = {5};
Physical Line(14) = {4};
Physical Line(15) = {6};
Physical Surface(1) = {9};
Physical Surface(2) = {10};

After that I convert the .msh file with "dolfin-convert" into .xml files.
I would like to solve the Stokes-equation on a domain consisting of a outer box and a circle in this box. On the box the nonslip condition should hold (here only on one site) and on the circle just a random velocity.

from dolfin import *

mesh = Mesh("circle.xml")
subdomains = MeshFunction("size_t", mesh, "circle_physical_region.xml")
boundaries = MeshFunction("size_t", mesh, "circle_facet_region.xml")

# Defining the function space
P2 = VectorElement("Lagrange", mesh.ufl_cell(), 2)
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
TH = P2 * P1
W = FunctionSpace(mesh, TH)

# Defining boundary conditions
noslip = Constant((0.0, 0.0))
particle_boundary = Constant((2.0, 3.0))
bc0 = DirichletBC(W.sub(0), noslip, boundaries, 12)
bc1 = DirichletBC(W.sub(0), particle_boundary, boundaries, 11)
bcs = [bc0, bc1]

# Defining the variational problem
(u, p) = TrialFunctions(W)
(v, q) = TestFunctions(W)
f = Constant((0.0, 0.0))
L = inner(f, v)*dx

# Solving the problem
U = Function(W)
solve( a==L, U, bcs)

u,p = U.split()
ufile_pvd = File("velocity.pvd")
ufile_pvd << u
pfile_pvd = File("pressure.pvd")
pfile_pvd << p

In the following are some screenshots of the solution with the 2017.2.0 version and the 1.6.0 version

Doing this with the 1.6.0 version I just have to change the function space to

V = VectorFunctionSpace(mesh, "Lagrange", 2)
Q = FunctionSpace(mesh, "Lagrange", 1)
W = V * Q​
This yields the right solution.

I set Fenics up with Docker if this matters.
So is there something fundamental I'm missing or has the way fenics works with subdomains changed?

Community: FEniCS Project
Unless you provide the way to produce XML files, there's nothing we can do about it.
written 6 months ago by Jan Blechta
I edited the .geo code for generating the mesh. I hope nothing is missing now.
written 6 months ago by J.H.59
Did you put any effort into making this example minimal?
written 6 months ago by Jan Blechta
I did but it seems it wasn't enough, sorry for that. Since the problem is the same in 2D I just changed the code to the 2D one, which is shorter.
written 6 months ago by J.H.59