How to represent non-uniform first and second type boundary conditions?

[French_fries]

Hello everyone, considering the simplest heat conduction equation, suppose the geometric model is a rectangle; for a boundary of the rectangle, if the temperature on that boundary is T0, it can be represented in Fenics as:

aT = dot(grad(dT), grad(T_))*dx
LT = Constant(0)*T_*dx
bc = DirichletBC(V, Constant(T0.), bottom)
solve(aT == LT, T, bc)

If the heat flux density on that boundary is q0, it can be represented in the variational formulation：

L=Constant(0)*T_*dx - q0*T_*ds

I want to know how to modify the above two parts of code for non-uniform boundary conditions(such as (T0, T1, T2…) or (q0, q1, q2…)), Not considering it can be written in the form of an expression.

I have looked at this method:https://fenicsproject.discourse.group/t/how-to-impose-a-pressure-distribution-from-a-data-file-as-a-boundary-condition-loading/9232/5?u=french_fries , but it does not involve the use of boundary condition statements such as bc = DirichletBC( ) or modifying the variational formulation for the heat flux density condition.

[dokken]

@French_fries

I have looked at this method:https://fenicsproject.discourse.group/t/how-to-impose-a-pressure-distribution-from-a-data-file-as-a-boundary-condition-loading/9232/5?u=french_fries , but it does not involve the use of boundary condition statements such as bc = DirichletBC( ) or modifying the variational formulation for the heat flux density condition.

As the post states, you can set the data to a dolfin.Function, that you in turn can use in your variational formulation (as q) or as a function in your DirichletBC.

[French_fries]

Thank you, I will now attempt to solve this simple problem based on the information provided. I will ask for your help with reviewing it afterwards.

[French_fries]

I made some attempts, the examples you provided were for applying data to the entire domain, but I encountered many issues when applying the data to the bottom boundary. Here is my code:

from dolfin import *
import numpy as np
mesh = UnitSquareMesh(4, 4)

VT = FunctionSpace(mesh, "CG", 1)
T_data = np.array([0, 50, 100, 150, 0])
T_coordinates = np.array([[0., 0.25.], [0, 0.5], [0., 0.75.], [0., 1]])
ind2 = np.lexsort((T_coordinates[:, 1], T_coordinates[:, 0]))
print(ind2)

class BottomBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 0)
bottomBoundary = BottomBoundary()
boundarymesh = BoundaryMesh(mesh, 'exterior')
bottom_mesh = SubMesh(boundarymesh, bottomBoundary)

WT = FunctionSpace(bottom_mesh, "CG", 1)
T_load = Function(WT)
dof_coords = WT.tabulate_dof_coordinates()
ind = np.lexsort((dof_coords[:, 1], dof_coords[:, 0]))
for (in_index, T_load_index) in zip(ind2, ind):
    T_load.vector()[T_load_index] = T_data[in_index]

def bottom(x, on_boundary):
    return near(x[1], 0) and on_boundary
def lateral_sides(x, on_boundary):
    return (near(x[0], 0) or near(x[0], 1)) and on_boundary
def top(x, on_boundary):
    return near(x[1], 1) and on_boundary

T_ = TestFunction(VT)
dT = TrialFunction(VT)
Delta_T = Function(VT, name="Temperature increase")
aT = dot(grad(dT), grad(T_))*dx
LT=Constant(0)*T_*dx

bcT = [DirichletBC(VT, T_load, bottom),
       DirichletBC(VT, Constant(0), lateral_sides),
       DirichletBC(VT, Constant(0), top)]
solve(aT == LT, Delta_T, bcT)

Firstly, T_coordinates cannot take the value [0, 0.25], even though I have divided the boundary into four equal parts; Then, I created a sub-mesh representing the bottom boundary and tried to calculate the degrees of freedom, but encountered some issues: What I need here are the degrees of freedom belonging to y=0 on the bottom boundary mesh; Perhaps you could modify my code to apply non-uniform data to the bottom boundary for the solution? I hope you can provide assistance, thank you!

[dokken]

First of all, your coordinates make no sense, as

@French_fries

T_coordinates = np.array([[0., 0.25.], [0, 0.5], [0., 0.75.], [0., 1]])

will never match with

@French_fries

class BottomBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 0)

Please pay special attention to details before posting code.
Your code is also way to complicated, introducing boundary meshes, sub meshes of boundary meshes etc.
Here is a minimal working example of what you want to achieve:

from dolfin import *
import numpy as np
mesh = UnitSquareMesh(4, 4)

VT = FunctionSpace(mesh, "CG", 1)
T_data = np.array([0, 50, 100, 150, 0])
T_coordinates = np.array([[0, 0], [0., 0.25], [0., 0.5], [0., 0.75], [0., 1.]])
assert len(T_data) == len(T_coordinates)


class LeftBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0) and on_boundary


leftBoundary = LeftBoundary()

vt = MeshFunction("size_t", mesh, 0, 0)
leftBoundary.mark(vt, 1)

vertices = vt.where_equal(1)
vertex_to_dof_map = vertex_to_dof_map(VT)

T_load = Function(VT)

# Get all boundary coordinates
dof_coords = VT.tabulate_dof_coordinates()
boundary_dofs = [vertex_to_dof_map[vertex] for vertex in vertices]
boundary_coords = np.array(
    [dof_coords[boundary_dof] for boundary_dof in boundary_dofs])

# Compute distances between boundary coordinates and T coordinates
points_A = np.expand_dims(boundary_coords, 1)
points_B = np.expand_dims(T_coordinates, 0)
distances = np.sum(np.square(points_A-points_B), axis=2)
is_close = distances < 1e2 * DOLFIN_EPS
positions = np.nonzero(is_close)
for row, col in zip(*positions):
    T_load.vector()[boundary_dofs[row]] = T_data[col]
    print(row, col, boundary_coords[row], T_coordinates[col], T_data[col])

bcT = [DirichletBC(VT, T_load, leftBoundary),
       ]

u = Function(VT)
[bc.apply(u.vector()) for bc in bcT]

File("mwe.pvd") << u

giving

0 0 [0. 0.] [0. 0.] 0
1 1 [0.   0.25] [0.   0.25] 50
2 2 [0.  0.5] [0.  0.5] 100
3 3 [0.   0.75] [0.   0.75] 150
4 4 [0. 1.] [0. 1.] 0

and function

[French_fries]

I’m sorry, my carelessness. I intended to apply the condition at the bottom but used the coordinates on the left side. I have a question regarding the two regions as shown in the figure (continuous conditions need to be satisfied at the dashed boundary). Suppose I want to transfer the computed results (temperature or heat flux density) at the dashed boundary of region 1 to the dashed boundary of region 2 in order to solve for region 2; should I transfer the discrete results to region 2 as we discussed above, or do I need to interpolate all the results on the dashed boundary of region 1 to the dashed boundary of region 2? (My understanding is that when solving the variational equation, discrete data points are used, and it is not necessary to transfer all the continuous results to the boundary of region 2, but doing so cannot guarantee that the data on the two boundaries are completely consistent). Can you explain this issue? If only discrete data points need to be transferred, I will use the modification results of your code for temperature, and I will handle the heat flux density in the same way in the future.

[dokken]

Since you transfer the data into another function space, you are just transferring the coefficients c_i of the expression uh=\sum_{i=1}^Nc_i\phi_i(x) where \phi_i is the $i$th global basis function. Thus by transferring the discrete values, you set them so that uh can be evaluated at any point.

[French_fries]

Thank you for your advice. I will use the modified code you provided and then apply the same method to the heat flux density. Thank you very much for your time!

[French_fries]

Hi，dokken.I looked at an example of partitioned coupled heat transfer.
https://github.com/precice/tutorials/tree/master/partitioned-heat-conduction
where the coupling boundary is fully matched (same length and number of nodes), but why do still need to use RBF interpolation when transferring temperature and heat flux? Here are snippets of the mesh and coupling boundary expressions.

def get_geometry(domain_part):
    nx = ny = 9

    if domain_part is DomainPart.LEFT:
        p0 = Point(x_left, y_bottom)
        p1 = Point(x_coupling, y_top)
    elif domain_part is DomainPart.RIGHT:
        p0 = Point(x_coupling, y_bottom)
        p1 = Point(x_right, y_top)
    else:
        raise Exception("invalid domain_part: {}".format(domain_part))

    mesh = RectangleMesh(p0, p1, nx, ny, diagonal="left")
    coupling_boundary = StraightBoundary()
    remaining_boundary = ExcludeStraightBoundary()

    return mesh, coupling_boundary, remaining_boundary

class SegregatedRBFInterpolationExpression(CouplingExpression):
    """
    Uses polynomial quadratic fit + RBF interpolation for implementation of CustomExpression.interpolate.
Allows for arbitrary coupling interfaces.

If possible, sorry for taking up some of your time.Based on our discussion yesterday and the current example, I want to mimic this example to solve the partitioned heat transfer problem I am working on. However, I am not sure whether interpolation is needed for the boundary with matched mesh nodes.

[French_fries]

Hi，dokken，I made some modifications based on this, used to solve the problem of non-uniform heat flux density boundary; the printed results show that the data has been successfully matched with the coordinates, and then I calculated the heat flux density for verification, but it was not successful (I modified the data and tried again, but the heat flux density results at the boundary still do not match the predetermined values).

from dolfin import *
import numpy as np
mesh = UnitSquareMesh(4, 4)

VT = FunctionSpace(mesh, "CG", 1)
T_ = TestFunction(VT)
dT = TrialFunction(VT)
q_data = np.array([0, 0.5, 1, 1.5, 0])
q_coordinates = np.array([[0, 0], [0., 0.25], [0., 0.5], [0., 0.75], [0., 1.]])
assert len(q_data) == len(q_coordinates)

class LeftBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0) and on_boundary
leftBoundary = LeftBoundary()
vt = MeshFunction("size_t", mesh, 0, 0)
leftBoundary.mark(vt, 1)

def other_sides(x, on_boundary):
    return near(x[0], 1) or near(x[1], 0) or near(x[1], 1) and on_boundary

vertices = vt.where_equal(1)
vertex_to_dof_map = vertex_to_dof_map(VT)
q_load = Function(VT)

# Get all boundary coordinates
dof_coords = VT.tabulate_dof_coordinates()
boundary_dofs = [vertex_to_dof_map[vertex] for vertex in vertices]
boundary_coords = np.array(
    [dof_coords[boundary_dof] for boundary_dof in boundary_dofs])

# Compute distances between boundary coordinates and T coordinates
points_A = np.expand_dims(boundary_coords, 1)
points_B = np.expand_dims(q_coordinates, 0)
distances = np.sum(np.square(points_A-points_B), axis=2)
is_close = distances < 1e2 * DOLFIN_EPS
positions = np.nonzero(is_close)
for row, col in zip(*positions):
    q_load.vector()[boundary_dofs[row]] = q_data[col]
    print(row, col, boundary_coords[row], q_coordinates[col], q_data[col])

ds = Measure('ds', domain=mesh, subdomain_data=vt)

lambda_ = 1   # thermal conductivity
Delta_T = Function(VT)
aT = dot(grad(dT), grad(T_))*dx
LT = Constant(1)*T_*dx - q_load*T_*ds(1)
bcT = [DirichletBC(VT, Constant(100), other_sides)]
solve(aT == LT, Delta_T, bcT)

#  Compute heat flux
Q = VectorFunctionSpace(mesh, 'CG', 1)
w = TrialFunction(Q)
v = TestFunction(Q)
aQ = inner(w, v)*dx
LQ = inner(-lambda_*grad(Delta_T), v)*dx
flux = Function(Q)
solve(aQ == LQ, flux)
File("q.pvd") << flux

And this is the result of the heat flux density at the boundary：

[French_fries]

I feel that the non-uniform heat flux density boundary problem cannot be achieved with just simple modifications

[dokken]

I dont think you can enforce both u=u_0 and du/dn=t_0 on the same boundary using a strong Dirichlet condition.

You might be able to enforce it weakly, but to me it Seems like your problem is over-constrained.

[French_fries]

Dear dokken.Sorry, except for the boundary where I need to impose a non-uniform heat flux density condition, the other conditions are arbitrary to me. I don’t care if this simple example has a solution or not. I just want to know how to apply the non-uniform heat flux density condition and solve it in FEniCS. I am interested in the process represented, just like the example you provided to me about the temperature condition. Do you have an example or any ideas?

[dokken]

Simply start by solving a manufactured problem (on a single domain), where a non-uniform flux is created to satisfy the boundary condition, see for instance Combining Dirichlet and Neumann conditions — FEniCSx tutorial

[French_fries]

I’m really sorry for the misunderstanding. The reason I asked you this question is that I need to calculate heat transfer in different regions, which involves using the heat flux density results from one region’s boundary as the boundary condition for solving another region (these heat flux density data are irregular points and cannot be written as an expression). Your previous example on temperature was very helpful! I have already learned from you how to extract data from the boundary in the previous post, so I posted this to ask how to apply these discrete data (temperature or heat flux density) to the boundary of another region and solve it.

[French_fries]

For the non-uniform heat flux density boundary condition issue, it would be great if you could provide me with an example similar to the one you provided for the non-uniform temperature problem. Thank you.

[dokken]

I simply dont have then bandwidth to write all the examples for you.

For fluxes I would:

1. Interpolate/project the computed flux on mesh 1 into a suitable function space.
2. transfer fluxes from mesh 1 to mesh 2 in the same way as for temperature
3. add these fluxes to the variational form

[xiaohe]

Thanks for your reply！
I’d like to continue asking questions based on the discussion above and another post.[coupling problem](https://fenicsproject.discourse.group/t/some-problems-about-solving-coupling-problem/13686)
mesh 1 and mesh 2 have a boundary in the same space location (or very very close ,consider machine errors). This boundary is called interface, but the coordinates of the two mesh on the interface do not match. In this case ,how to transfer a fem function u1(defined mesh1) to a fem function u2 (defined mesh2) on the interface?
A simple code as follows:

from fenics import *
from mshr import *
domain01 = Rectangle(Point(0.0, 0.0), Point(1.0, 1.0))
domain02 = Rectangle(Point(1.0, 0.0), Point(2.0, 1.0))
mesh01 = generate_mesh(domain01, 0.05)
mesh02 = generate_mesh(domain02, 0.01)

# mark boundary
class Interface(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], 1.0)
# note: the coordinates of mesh1 and mesh1 on the interface do not match
tb01 = Interface()
bm01 = MeshFunction("size_t", mesh01, mesh01.topology().dim()-1, 0)
bm01.set_all(0)
tb01.mark(bm01, 1)
ds01 = Measure("ds", domain= mesh01, subdomain_data= bm01)

bm02 = MeshFunction("size_t", mesh02, mesh02.topology().dim()-1, 0)
bm02.set_all(0)
tb01.mark(bm02, 2)
ds02 = Measure("ds", domain= mesh02, subdomain_data= bm02)

V01 = FunctionSpace(mesh01, "P", 2)
u01 = Function(V01)
u01.set_allow_extrapolation(True)
# note :Suppose u1 is solved from a pde system
V02 = FunctionSpace(mesh02, "P", 2)
u02 = Function(V02)
# u02.interpolate(u01)
u02 = project(u01, V02)

[xiaohe]

Another question is about the fenics Function MeshFunction().
I’d like to know what the differences is between the following two code:

# mesh01 is a 2d mesh 
bm01 = MeshFunction("size_t", mesh01, mesh01.topology().dim()-1, 0)
and 
bm01 = MeshFunction("size_t", mesh01, 0, 0)

They all seem to be used to mark boundarys?

[dokken]

@xiaohe

bm01 = MeshFunction("size_t", mesh01, mesh01.topology().dim()-1, 0)

This marks all facets of a mesh with the value 0.

@xiaohe

bm01 = MeshFunction("size_t", mesh01, 0, 0)

This marks all vertices with the value 0.

@xiaohe

mesh 1 and mesh 2 have a boundary in the same space location (or very very close ,consider machine errors). This boundary is called interface, but the coordinates of the two mesh on the interface do not match. In this case ,how to transfer a fem function u1(defined mesh1) to a fem function u2 (defined mesh2) on the interface?

Again, this is covered with How to represent non-uniform first and second type boundary conditions? - #6 by dokken