Defining local boundary conditions is not working in the current version of FEniCS

1

In previous versions of FEniCS, I was able to define specific boundary conditions in the domain, like this.

import numpy as np
from dolfin import *
from mshr import *
import matplotlib.pyplot as plt

pcoorx_obs = [1800.0]
pcoory_obs = [1600.0]
hp_obs = [20.9135966045882817]
pcoorx_v = [2100.0]
pcoory_v = [2100.0]
flow = [-2000.0]  # m3/s
valork = 20.0  
b = 10.0  
u_w = Constant(11.0)
u_e = Constant(11.0)

x0 = 0.0
y0 = 0.0
x1 = 4200.0
y1 = 4200.0
resolution = 30


rect = Rectangle(Point(x0, y0), Point(x1, y1))
mesh = generate_mesh(rect, resolution)

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


def uw_boundary(x):
    return x[0] < DOLFIN_EPS

def ue_boundary(x):
    return near(x[0], x1, DOLFIN_EPS)

def condition_p0(x): 
    return near(x[0], pcoorx_obs[0], DOLFIN_EPS) and near(x[1], pcoory_obs[0], DOLFIN_EPS)
     
     
bcw = DirichletBC(V, u_w, uw_boundary)
bce = DirichletBC(V, u_e, ue_boundary)
bcp_0 = DirichletBC(V, hp_obs[0], condition_p0, method='pointwise')  # This condition is not working!!!!!

bc_obs = [bcw, bce, bcp_0]


class K(UserExpression):
    def __init__(self, mesh, **kwargs):
        self.mesh = mesh
        super().__init__(**kwargs)
        
    def eval(self, values, x):
        T = valork * b
        values[0] = T


T = K(mesh, degree=1)


u_obs = TrialFunction(V)
v_obs = TestFunction(V)
a_obs = inner(T*nabla_grad(u_obs), nabla_grad(v_obs))*dx
L_obs = Constant(0.0)*v_obs*dx
A_obs = assemble(a_obs)
rhs1_obs = None
rhs_obs = assemble(L_obs, tensor=rhs1_obs)
u_obs = Function(V)
for bc in bc_obs:
    bc.apply(A_obs, rhs_obs)
for i in range(len(flow)):
    delta_obs = PointSource(V, Point(pcoorx_v[i], pcoory_v[i]), flow[i])
    delta_obs.apply(rhs_obs)
solve(A_obs, u_obs.vector(), rhs_obs)


plt.figure()
plot(mesh)

plt.figure()
plot(u_obs, title='Head')

plt.show()

In older versions this condition worked bcp_0 = DirichletBC(V, hp_obs[0], condition_p0, method='pointwise'). but now it has no effect.

2

I figured out one of the reasons why it was not working. In the old code (dolfin 1.6) I used this function to change the coordinates of some points of the mesh to the desired location, then created an xml file of the mesh. The problem is that this function is not working in this version (2019.1).

def MeshPassingToWellAsVertex(mesh,coorx,coory):
    ''' author: jpmdossantos@yahoo.com.br
    Based on code from Fenics website and Q&A'''
    #Example From: https://fenicsproject.org/qa/9763/method-compute_closest_point
    coords=mesh.coordinates()
    print(np.shape(coords))
    
    tree = mesh.bounding_box_tree()
    
    point_cloud = [dolfin.Point(point) for point in mesh.coordinates()]
    tree.build(point_cloud, 2)
    p=Point(coorx[0],coory[0])
    
    p_i, distance = tree.compute_closest_point(p)
    print (p_i)
    
    lengthCoor=np.shape(coords)[0]
    print('lenght coordinates matrix Xdirection', lengthCoor)
    for i in range(len(coorx)):
        print (i)
        p=Point(coorx[i],coory[i])
        #tree = mesh.bounding_box_tree()
        #point_cloud = [dolfin.Point(point) for point in mesh.coordinates()]
        #tree.build(point_cloud, 2)
        p_i, distance = tree.compute_closest_point(p)   
        print ("p:", p.str())
        print ("closest coordinate point by reference to index",coords[p_i])
        print ("closest point:",p_i,point_cloud[p_i].str())
        
        
        coords[p_i][0]=coorx[i];coords[p_i][1]=coory[i];
          
        #XMesh= coords[vertex.index()][0];YMesh= coords[vertex.index()][1]
    return mesh


mesh = generate_mesh(rect, resolution)
mesh = MeshPassingToWellAsVertex(mesh, pcoorx_obs, pcoory_obs)


File("meshToWellAsVertex.xml")<<mesh
mesh=Mesh("meshToWellAsVertex.xml")

the result I expected to achieve was something like this.

head
head1058×782 86.1 KB