Error in defining boundaries and incorporating a scalar function into a matrix

1

In the MWE below, I am trying to solve a Poisson equation on a unit disc with a 0 Dirichlet BC. I have defined a Lagrange FunctionSpace on the domain. The solution is a scalar function defined at every point on the domain. Next, I am forming two other scalar functions, u1 and u2 from the solution u. Next, I want to form a matrix G that takes as input the two scalar functions u1 and u2.

Now, when I am trying to execute the code, it’s producing the error for each iteration.

*** Warning: Found no facets matching domain for boundary condition.
  *** Warning: Found no facets matching domain for boundary condition.

Python script:

C1 = Circle(Point(0, 0),1)

class b1(SubDomain):
  def inside(self, x, on_boundary):
      return near(x[0]**2 + x[1]**2, 1, eps=1e-3) and on_boundary

class s1(SubDomain):
    def inside(self, x, on_boundary):
        return x[0]**2 + x[1]**2 < 1 
        
        
mesh = generate_mesh(C1, 5)


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

b1().mark(boundary_markers, 1)

s1().mark(surface_markers, 8)

ds = Measure('dS', domain=mesh, subdomain_data=boundary_markers)
dx = Measure('dx', domain=mesh, subdomain_data=surface_markers)


W1 = FunctionSpace(mesh, "Lagrange", 1)


du = TrialFunction(W1)
v  = TestFunction(W1)             
u  = Function(W1)                 


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


bc = DirichletBC(W1, u_D, boundary_markers, 1)


# Define variational problem

f = Constant(-G/mu_flow) 
a = inner(grad(u), grad(v))*dx(domain=mesh) + f*v*dx(domain=mesh)


# Compute solution

solve(a == 0, u, bc)


u1 = Expression('(u**4)/((u**4)+(1.12*(10**-5))**4)')     
u2 = Expression('(0.25)*(u**4)/((u**4)+(1.12*(10**-5))**4)')


G1 = Expression((('1.0+u1','0.0','0.0'),
                     ('0.0','1.0+u2','0.0'),
                     ('0.0','0.0','1.0')))

My final aim is to define the matrix G at every vertex of the mesh. What am I doing wrong here? I know there is a problem with my defining the boundaries, but can’t quite figure out the problem here. Any help or suggestion would be highly appreciated.

2

Maybe adding <= does the trick or/and a padding with eps does the trick

class s1(SubDomain):
    def inside(self, x, on_boundary):  
        eps = 1e-3
        return x[0]**2 + x[1]**2 <= 1 + eps

the way I see it, b1 is just barely not in s1. Also your mesh wont be perfectly round, therefore a padding might be needed.

Maybe this helps :slight_smile:

Emanuel

3

Hi @Emanuel_Gebauer , thanks for the suggestion. However, the code you suggested also gives me the same warnings as before.

However, I created this MWE for a more complex problem, for which I import my mesh from gmsh. So the above problem doesn’t appear in that case. My main concern is how to calculate the determinant of an expression matrix and incorporate that into my strain-energy density function. Suppose I create the same mesh in gmsh and import it to FEniCS and introduce a FunctionSpace on the mesh and solve the variational form as

W1 = FunctionSpace(mesh, "Lagrange", 1)


du = TrialFunction(W1)
v  = TestFunction(W1)             
u  = Function(W1)                 


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


bc = DirichletBC(W1, u_D, boundary_markers, 1)


# Define variational problem

f = Constant(-G/mu_flow) 
a = inner(grad(u), grad(v))*dx(domain=mesh) + f*v*dx(domain=mesh)


# Compute solution

solve(a == 0, u, bc)

u1 = Expression(("pow(u_0,4)/(pow(u_0,4)+pow((1.12*pow(10,-5)),4))"), degree=1, u_0=u) 

u2 = Expression(("(0.25)*pow(u_0,4)/(pow(u_0,4)+pow(1.12*pow(10,-5),4))"),degree=1, u_0=u) 

g = Expression((('1.0+c1','0.0','0.0'),
                     ('0.0','1.0+c2','0.0'),
                     ('0.0','0.0','1.0')),degree=0,c1=u1,c2=u2)
                     
               
g1 = project(g,W1)

I need to find out the determinant of g at every node of the mesh. But it’s raising the error,

g1 = project(g,W1)
.....
ufl.log.UFLException: Shapes do not match:

So the main question is how to find the determinant, det(g1), for each of the nodes on the mesh?

4

Please make sure that your minimal example can reproduce your error.
As you have not supplied any mesh for this problem, the code above cannot be executed.
I would suggest you start with a UnitCubeMesh and work your way from there.

In general for having robust tagging of boundaries, I would suggest using an external mesh generator, such as GMSH, and import the mesh using meshio as explained in several other posts of the form. Then you can import your domain and boundary markers.

For your second issue, you are trying to project a tensor of shape (3,3) into a scalar space (W1 = FunctionSpace(mesh, "Lagrange", 1)). This cannot work, as you would need to project g into a TensorFunctionSpace of shape (3,3).

5

Thanks @dokken, that works.