Unable to apply BCs

1

Hello

I am trying to apply Dirichlet BCs on the upper-right and lower-left 10% corner edges of a square domain, see the following sketch and example code. There was no error in running the code, but by printing the max and min of the solution u, I expected min to be 0 and max be 1 but the code gave me both around 1. I might just not seeing my own bugs, could anyone help taking a look? Thank you!

Screen Shot 2021-03-24 at 10.49.57 PM

from dolfin import *

class ExternalSurf(SubDomain):
    def inside(self, x, on_boundary):
        top = ( (x[0] > 4.5) and (x[1] >5.-DOLFIN_EPS) ) 
        right = ( (x[1] > 4.5) and (x[0] > 5.-DOLFIN_EPS) )
        return   (top or right) and  on_boundary
class InternalSurf(SubDomain):
    def inside(self, x, on_boundary):
        bottom = ( (x[0] < 0.5) and (x[1] <DOLFIN_EPS) )
        left = ( (x[1] < 0.5) and (x[0] <DOLFIN_EPS) )
        return  (bottom or left) and on_boundary
mesh = RectangleMesh(Point(0., 0.), Point(5., 5.), 100, 100)
V = FunctionSpace(mesh, 'Lagrange', 2)
bc1 = DirichletBC(V, Constant(1.), InternalSurf())
bc2 = DirichletBC(V, Constant(0.), ExternalSurf())
bc = [ bc1, bc2 ]
u = Function(V)
v = TestFunction(V)
F = inner( grad(u), grad(v) )*dx
solve(F==0,u,bc)
print('min(u.vector()) ',min(u.vector()))
print('max(u.vector()) ',max(u.vector()))

the code does give the following warning though:

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

2

Replace your boundary checks with near function. Like this.

    from dolfin import *
    
    
    class ExternalSurf(SubDomain):
        def inside(self, x, on_boundary):
            top = (x[0] > 4.5) and near(x[1], 5.0)
            right = (x[1] > 4.5) and near(x[0], 5.0)
            return (top or right) and on_boundary
    
    
    class InternalSurf(SubDomain):
        def inside(self, x, on_boundary):
            bottom = (x[0] < 0.5) and near(x[1], 0.0)
            left = (x[1] < 0.5) and near(x[0], 0.0)
            return (bottom or left) and on_boundary
    
    
    mesh = RectangleMesh(Point(0.0, 0.0), Point(5.0, 5.0), 100, 100)
    V = FunctionSpace(mesh, "Lagrange", 2)
    bc1 = DirichletBC(V, Constant(1.0), InternalSurf())
    bc2 = DirichletBC(V, Constant(0.0), ExternalSurf())
    bc = [bc1, bc2]
    u = Function(V)
    v = TestFunction(V)
    F = inner(grad(u), grad(v)) * dx
    solve(F == 0, u, bc)
    print("min(u.vector()) ", min(u.vector()))
    print("max(u.vector()) ", max(u.vector()))

Output:

    No Jacobian form specified for nonlinear variational problem.
    Differentiating residual form F to obtain Jacobian J = F'.
    Solving nonlinear variational problem.
      Newton iteration 0: r (abs) = 1.596e+01 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
      Newton iteration 1: r (abs) = 1.210e-11 (tol = 1.000e-10) r (rel) = 7.581e-13 (tol = 1.000e-09)
      Newton solver finished in 1 iterations and 1 linear solver iterations.
    min(u.vector())  0.0
    max(u.vector())  1.0
1 Like