Hello,
Suppose I have a list of fluxes, one flux q_i per face A_i on the surface Γ of my mesh. My goal is to prescribe a Neumann boundary condition g on Γ in my Continuous Galerkin FEM solver for the heat equation, where g must satisfy
\int_{A_i} g ds = |A_i|*q_i.
So g is essentially the projection of a piece-wise constant function q in L2(Γ) onto whatever space we choose to define g in. (I suppose g in H1(Γ) would be most natural for a CG solver.)
Given a list of pairs (i, q_i) corresponding to the value of q on the face i, can you recommend an efficient way to obtain the boundary condition g?
I have looked into manually assembling q as a GridFunction, but I got stuck trying to match the indices in q.vec with the ids of the mesh vertices. For GridFunctions in H1, there is a nice one-to-one correspondence, but not so with L2. And if I got past this step, I would still be unsure about what is the best approach to project q from L2 into H1.
For context: I am asking because I have a solver for non-local radiation problems that outputs pairs (i, q_i), and I need to convert these into a suitable boundary condition for my heat equation solver.
Thanks in advance for any guidance!
Best regards,
Sindre
In case my description was not entirely clear, I have attempted to construct a kind of conceptual MWE.
Let us assume that my domain is the equilateral triangle with vertices v_1=(0,0), v_2=(1,0), and v_3=(0.5, \sqrt{3}/2), and that I mesh it using only a single element. The mesh of the surface Γ then consists of three faces: A_1=v_1 \rightarrow v_2, A_2=v_2 \rightarrow v_3, and A_3=v_3 \rightarrow v_1.
Suppose that my radiation solver outputs q_1 = 1, q_2 = 2 and q_3 = 3. Then, I want to obtain a function g that averages 1 over A_1, 2 over A_2 and 3 over A_3. Moreover, with
g = ng.GridFunction(ng.H1(order=1)),
g should, in this example, be the unique piecewise linear function satisying g(v_1) = 2, g(v_2) = 0, g(v_3) = 4.
The purpose of obtaining g is to assign it as a boundary condition in my solver for the heat equation.
You need to formulate a variational problem for the L^2 projection on the boundary. That way you can solve for the g function as if you are solving an elliptic problem.
Since you want piecewise constants on the boundary, you will need the FacetFESpace (instead of L2). You need your q function in that space.
For the g function, I suppose you want it in a V=H1 space on the boundary (to make it continuous). You will then define a bilinear form to be an L^2 inner product on V.
Sorry this is not clear. I don’t know a specific example or demo for this. Also FacetFESpace is not that well documented, but you can always do help(FacetFESpace) in python.
Perhaps another user can give some pointers.
Thanks for your input, walker_tr9!
Your suggested approach sounds like it should work, so I had a go at implementing it. Unfortunately, my implementation (given in the first code block below) does not produce a GridFunction g with the correct face-wise averages. Why I do not know.
from netgen.geom2d import unit_square
import ngsolve as ng
import numpy as np
mesh = ng.Mesh(unit_square.GenerateMesh(maxh=0.1))
# Define q
L2_bnd = ng.FacetFESpace(mesh, order=0, definedon=mesh.Boundaries("top"), dgjumps=True)
q = ng.GridFunction(L2_bnd)
q_vals = [i+1 for i in range(len(q.vec))]
for i in range(len(q.vec)):
q.vec[i] = q_vals[i]
print(q.vec)
# Calculate g
V = ng.H1(mesh, order=1, definedon=mesh.Boundaries("top"))
g = ng.GridFunction(V)
u,v = V.TnT()
a = ng.BilinearForm(u*v*ng.ds(definedon=mesh.Boundaries("top"))).Assemble()
f = ng.LinearForm(q*v*ng.ds(definedon=mesh.Boundaries("top"))).Assemble()
g.vec.data = a.mat.Inverse(freedofs=V.FreeDofs()) * f.vec
print(g.vec)
# Validate (fails)
for i, face in enumerate(mesh.Boundaries("top").Elements()):
node_vals = [g.vec[v.nr] for v in face.vertices]
print("node_vals", node_vals)
avg = sum(node_vals)/len(node_vals)
np.testing.assert_approx_equal(avg, q_vals[i])
I have also tried a different approach based on np.linalg.lstsq. This approach produces a piece-wise linear g with correct face-wise averages for a surface mesh of triangles, but cannot take into account higher-order surface approximations of curved surfaces. In the implementation below, I have used standard numpy arrays. For large problems, I imagine it would be better to use some sparse matrix structure instead.
import numpy as np
from netgen.geom2d import unit_square
import ngsolve as ng
from ngsolve.webgui import Draw
mesh = ng.Mesh(unit_square.GenerateMesh(maxh=0.1))
# Extract faces and nodes on the relevant boundary segment.
top_faces = []
top_nodes = {}
for face in mesh.Boundaries("top").Elements():
top_faces.append(face.nr)
top_nodes[face.nr] = [v.nr for v in face.vertices]
# Construct connectivity matrix.
N = len(mesh.nodes(ng.VERTEX))
A = np.zeros((len(top_faces), N))
c = 1./mesh.dim
for i, f_nr in enumerate(top_faces):
for v_nr in top_nodes[f_nr]:
A[i][v_nr] = c
# Define desired mean of flux across each face.
#qs = np.random.uniform(1, 12, size=len(top_faces))
qs = [i+1 for i in range(len(top_faces))]
# Construct GridFunction from solution of least squares problem.
g = ng.GridFunction(ng.H1(mesh, order=1))
g_data = np.linalg.lstsq(A, qs)[0]
for i in range(g_data.shape[0]):
g.vec[i] = g_data[i]
Draw(g, mesh)
# Check correctness.
for i, face in enumerate(mesh.Boundaries("top").Elements()):
node_vals = [g.vec[v.nr] for v in face.vertices]
avg = sum(node_vals)/len(node_vals)
np.testing.assert_approx_equal(avg, qs[i])
I’m not surprised the first approach did not keep the face averages. You are projecting a piecewise constant function to a continuous piecewise linear. At best, you would preserve the average flux over groups of elements that share a vertex.
I’m not sure about the other approach, but if it does what you need, then fine. Are you doing a least squares fit over groups of elements that share a vertex? I guess you are just averaging the nodal values over a triangle.