Apply Neumann boundary condition in a specific area, which is part of a surface

[Yanjun]

Dear all,

I try to add heat flux(Neumann condition) to a surface, while this suface is large and I only want to input heat flux to specific area.

Previous, I have successfully add Dirchlet boundary conditions, as shown in below

bc_list = []
for i, (x_c, y_c) in enumerate(zip(x_co, y_co), start=1):
    var_name = f"bc_facets{i}"
    
    tag= mesh.locate_entities_boundary(
       domain, fdim, lambda x: (x[0]-x_c)**2 + (x[1]-y_c)**2  <= 324 )
        
    globals()[var_name] = tag
    bc_list.append(globals()[var_name])
    # rubbing elements

bc_disc = mesh.locate_entities_boundary(
    domain, fdim , lambda x:  np.isclose(x[2], 20) )

bc_all = find_common_e(bc_disc, bc_list)

Through above codes, I can add dirichlet boundary condition. My variation form like below:

a = (rho*c)/dt * u * v *dx + k* dot(grad(u), grad(v)) *dx 
-  k * dot(grad(u)*v, n)*ds
# g is the heat flux, should applied to specific areas of the surface(21)
L = f *v *dx + (rho*c)/dt * u_n * v * dx   +  g* v *ds(21)

From tutorial about Combining Dirichlet and Neumann conditions, the neumann condition g can be expressed as

x = SpatialCoordinate(mesh)
g = -4 * x[1]
f = Constant(mesh, default_scalar_type(-6))
L = f * v * dx - g * v * ds

Then, how should I define my g? which should fuifill requirement like

g = 100 * (   (x[0]+214)**2 + (x[1]+27)**2 <= 324 && x[2] ==20   )

But above g, the python is not support.
Thanks for your help!

[francesco-ballarin]

Your g seems to be defined based on inequalities on the spatial coordinate x, which define multiple boundary parts. The simplest approach would be to assign multiple markers (rather than only the marker 21) so that each part is marked with a separate label.

Then define separate expressions for g, and use them with ds as shown eg in Setting multiple Dirichlet, Neumann, and Robin conditions — FEniCSx tutorial

[Yanjun]

Hi Francesco, thanks for your reply.

Exactly, my g is defined based on inequalities of x. This part is a circle on surface 21. This is a transient model and this circle will rotate rapidly, around 100 new positions per second. Using separate labels mean lots of work?
All of these circles are on surface 21, while each time step is a different coordinate. Is there any other method?

Thanks very much!

[dokken]

@Yanjun

g = 100 * ( (x[0]+214)**2 + (x[1]+27)**2 <= 324 && x[2] ==20 )

I would use ufl.conditional,
see for instance How to assign boundary conditions based on values of dolfinx.Function - #2 by dokken or Form compilation — FEniCS Tutorial @ Sorbonne

[Yanjun]

Thank you @dokken very much for help me! I haven’t try this method since I have to solve a new problem first: Apply two Neumann boundary conditons at same time.

I have read the Setting multiple Dirichlet, Neumann, and Robin conditions. I download this jupyternotebook and find for ds, ds(3) and ds(4) Mesh Tags objects are the same, as shown in below picture. I am a little confused since it should different?

.

I tried two neumann boundary condions, ds1(10) and ds2(20) in my project. When one one of them exist, it works. Like below picture1 adn picture2

a = (rho*c)/dt*inner(u,v)*dx+k*inner(grad(u),grad(v))*dx -k * dot(grad(u) *v, n) * ds1(10)
L = inner(f,v) *dx + (rho*c)/dt * inner(u_n, v) * dx  +  g  * v * ds1(10)

a = (rho*c)/dt*inner(u,v)*dx+k*inner(grad(u),grad(v))*dx -k * dot(grad(u) *v, n) * ds2(20)
L = inner(f,v) *dx + (rho*c)/dt * inner(u_n, v) * dx  +  g  * v * ds2(20)

ds1(10) and ds2(20) are different sources, but when I tried to add this two neumann conditions together,

a = (rho*c)/dt*inner(u,v)*dx+k*inner(grad(u),grad(v))*dx- k*dot(grad(u)*v,n)*ds1(10)- k*dot(grad(u)*v,n)*ds2(20)
    L = inner(f,v) *dx + (rho*c)/dt * inner(u_n, v) * dx  +g*v*ds1(10) +g*v*ds2(20)

The bilinear form (a) can not create, it seems ds1 and ds2 should combine together. Then I defined ds1(10) and ds2(20) togeher:

for i, (x_c, y_c) in enumerate(zip(x_co, y_co), start=1):
    vari_name = 10*i
    boundaries.append(  (vari_name, lambda x: (x[0]-x_c)**2 + (x[1]-y_c)**2  <= 324)  )
   ...
facet_tag = meshtags(domain, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
ds = Measure('ds', domain=domain, subdomain_data=facet_tag) 

print(boundaries)
[(10, <function <lambda> at 0x744a500971c0>), 
(20, <function <lambda> at 0x744a50cfdcf0>)]

print(ds(10))
print(ds(20))
ds(subdomain_id=10, domain=<Mesh #0>, subdomain_data=<dolfinx.mesh.MeshTags object at 0x744a50d11060>)
ds(subdomain_id=20, domain=<Mesh #0>, subdomain_data=<dolfinx.mesh.MeshTags object at 0x744a50d11060>)

I check the linera form (L) and (a), are all correct,

print(L)
{ c_0 * (conj((v_0))) } * dx(<Mesh #0>[everywhere], {})
  +  { 0.036266999999999994 * u_n * (conj((v_0))) } * dx(<Mesh #0>[everywhere], {})
  +  { 18.297971999999998 * (conj((v_0))) } * ds(<Mesh #0>[10], {})
  +  { 18.297971999999998 * (conj((v_0))) } * ds(<Mesh #0>[20], {})

But the result is still only one heat source, which is ds(10):

Thank you very much for help me with this issue! I appreciate all of your time.

[dokken]

@Yanjun

I have read the Setting multiple Dirichlet, Neumann, and Robin conditions. I download this jupyternotebook and find for ds, ds(3) and ds(4) Mesh Tags objects are the same, as shown in below picture. I am a little confused since it should different?

As you can see, the subdomain_id in each of them are different. this is what restricts the integration to those facets in facet_tag marked with either 3 or 4.

For the rest of your post, it is quite hard to give any guidance, as you haven’t provided a minimal working example.
Start by checking the output of dolfinx.fem.assemble_scalar(dolfinx.fem.form(1*ds(10)))
and
dolfinx.fem.assemble_scalar(dolfinx.fem.form(1*ds(20)))
to ensure that they are bot non-zero.

[Yanjun]

Thanks for swift reply!

regards
yanjun

[Yanjun]

both ds0 and ds1 are no zero. Now the issue is how to form bilinear form (a), I tried

bilinear_form = fem.form(a)
linear_form = fem.form(L)

which shows

File /usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/forms.py:183, in form.<locals>._create_form(form)
    180 """Recursively convert ufl.Forms to dolfinx.fem.Form, otherwise
    181 return form argument"""
    182 if isinstance(form, ufl.Form):
--> 183     return _form(form)

My variation formula is

a = (rho*c)/dt*inner(u,v)*dx+k*inner(grad(u),grad(v))*dx- k*dot(grad(u)*v,n)*(ds[0](10)) - k*dot(grad(u)*v,n)*(ds[1](20))

But only ds0 or ds1 works

a = (rho*c)/dt*inner(u,v)*dx+k*inner(grad(u),grad(v))*dx- k*dot(grad(u)*v,n)*(ds[0](10))
or
a = (rho*c)/dt*inner(u,v)*dx+k*inner(grad(u),grad(v))*dx- k*dot(grad(u)*v,n)*(ds[1](20))

Thanks for your help.

[dokken]

As you haven’t provided a minimal reproducible example it is really hard to give you any further guidance.

[Yanjun]

Hi dokken, thanks for reply. now I solved this issue by very stupid method, I use

boundary1= (10, lambda x: (x[0]-x_co[0])**2 + (x[1]-y_co[0])**2 <= 324)
boundary2= (20, lambda x: (x[0]-x_co[1])**2 + (x[1]-y_co[1])**2 <= 324) 
boundary3= (30, lambda x: (x[0]-x_co[2])**2 + (x[1]-y_co[2])**2 <= 324)
boundary4= (40, lambda x: (x[0]-x_co[3])**2 + (x[1]-y_co[3])**2 <= 324) 

boundaries = [boundary1,boundary2,boundary3,boundary4]

rather than

for i, (x_c, y_c) in enumerate(zip(x_co, y_co), start=1):
    vari_name = 10*i
    boundaries.append(  (vari_name, lambda x: (x[0]-x_c)**2 + (x[1]-y_c)**2  <= 324)  )

although when print these two ways, it shows the same, like

print(boundaries)
[(10, <function <lambda> at 0x7789984d51b0>), 
(20, <function <lambda> at 0x7789984d75b0>), 
(30, <function <lambda> at 0x7789984d6d40>), 
(40, <function <lambda> at 0x7789984d4c10>)]

But for the facets number, the second method is incorrect. As for the linear and bilinear form, I still use

a = (rho*c)/dt*inner(u,v)*dx+k*inner(grad(u),grad(v))*dx 
L = inner(f,v) *dx + (rho*c)/dt * inner(u_n, v) * dx  

for i in list(range(1,19)):
     a +=   -  k * dot( grad(u)* v, n) *ds(10*i)
     L +=   +  inner(g, v) *ds(10*i)

Now this issue has been solved, next I try ufl.condition. Hope I can finish model and submit an abstract for FEniCSx2024.

[Yanjun]

Find the solution to question
‘’
I try to add heat flux(Neumann condition) to a surface, while this suface is large and I only want to input heat flux to specific area.
‘’
The idea is to find commont facets indicies:

bc_disc = mesh.locate_entities_boundary(
    domain, fdim , lambda x:  np.isclose(x[2], 20) )
facet_indices = np.intersect1d(bc_disc, facet_indices)

Here np.intersect1d help to find the common facets indicies of two areas.

Thanks for all of your help!