With legacy Fenics I used to calculate the heat flux from the solved temperature variable with
V = VectorFunctionSpace(mesh, ('CG', 1))
heat_flux = project(k*grad(T), V)
but I can’t figure out how to do that with Fenicsx, or did I miss the relevant documentation?
See Where to find 'project'-function in dolfinx? - #2 by nate.
See also interpolation of ufl forms in dolfinx, which I believe fits your requirements: https://github.com/FEniCS/dolfinx/blob/main/python/demo/demo_elasticity.py#L229-L232
Also bear in mind that depending on the underlying function spaces of k and T, projection or interpolation onto a C^0 basis is likely not appropriate. Consider a discontinuous space.
I had already tried Where to find 'project'-function in dolfinx? - #2 by nate, but it crashes with the error AttributeError: FunctionSpace object has no attribute function_space.
And interpolation complains that RuntimeError: Cannot interpolate Expression with Argument
So I’ll keep looking.
It would be easier to help you if you provide a minimal reproducible example, ie a stand alone code the produces the error messages you are posting above
My mistake was to use the legacy Fenics solve(a == L, T, bcs) instead of the new Fenicsx T = problem.solve() so that T post-processing was yielding these errors messages.
Maybe to try to upgrade my existing legacy Fenics script is not the best appraach and it is better to start from scratch.
So problem solved, thank you for your help!
After solving the temperature distribution in a linear heat transfer problem, I want to extract the heat flux density in the z direction.
In Fenics legacy the following method was yielding a correct result shown in the picture fenics_legacy_heat_flux.png below:
V = VectorFunctionSpace(mesh, 'CG', 1)
qw = project(k*grad(T), V, solver_type="cg")
qw_x, qw_y, qw_z = qw.split(deepcopy=True)
Heat_Flux_z = Function(V, name='Heat_Flux_z')
Heat_Flux_z.assign(qw_z)
Now with Fenicsx and the following method I get the incorrect results shown in fenicsx_heat_flux.png below:
V = VectorFunctionSpace(mesh, ("Discontinuous Lagrange", 0))
qw_expr = Expression(k*grad(T), V.element.interpolation_points())
qw = Function(V)
qw.interpolate(qw_expr)
qw_x, qw_y, qw_z = qw.split()
Heat_Flux_z = Function(V, name='Heat_Flux_z')
Heat_Flux_z.x.array[:] = qw_z.x.array
I also tried
V = VectorFunctionSpace(mesh, ("Discontinuous Lagrange", 0))
qw_expr = Expression(k*grad(T), V.element.interpolation_points())
qw = Function(V)
qw.interpolate(qw_expr)
geometry = qw.function_space.tabulate_dof_coordinates()
V = FunctionSpace(mesh, ("Discontinuous Lagrange", 0))
qw_z = Function(V, name='Heat_Flux_z')
qw_z.x.array[:] = qw.x.array.real.reshape((geometry.shape[0], len(qw)))[:, 2]
but it does not yield the correct results either. Any idea what I am doing wrong here?
I don’t follow the logic here.
You are creating a vector function space V for the gradient.
Then you say that the flux in the z direction is also part of this space.
Are you sure you don’t want
Heat_Flux_z= qw_z.collapse()
The solution was twofold:
qw_z.collapse() as you suggested('CG', 1) rather than ("Discontinuous Lagrange", 0) V = VectorFunctionSpace(mesh, ('CG', 1))
qw_expr = Expression(k*grad(T), V.element.interpolation_points())
qw = Function(V)
qw.interpolate(qw_expr)
qw_x, qw_y, qw_z = qw.split()
qw_z = qw_z.collapse()
V = FunctionSpace(mesh, ('CG', 1))
Heat_Flux_z = Function(V, name='Heat_Flux_z')
Heat_Flux_z.x.array[:] = qw_z.x.array
scipy.interpolate.griddata(), is a bit less smooth than with Fenics legacy, but absolutely usable.You should really be careful with
as you are interpolating into a function space that the expression is not in, yielding arbitrary behavior in parallel (and is also not well defined).
How to do it better? Using ("Discontinuous Lagrange", 0) yields completely wrong results.
Using
V = VectorFunctionSpace(mesh, ("Discontinuous Lagrange", 0))
qw_expr = Expression(k*grad(T), V.element.interpolation_points())
qw = Function(V)
qw.interpolate(qw_expr)
qw_x, qw_y, qw_z = qw.split()
qw_z = qw_z.collapse()
V = FunctionSpace(mesh, ("Discontinuous Lagrange", 0))
Heat_Flux_z = Function(V, name='Heat_Flux_z')
Heat_Flux_z.x.array[:] = qw_z.x.array
yields
You have not provided the code you use for plotting the results, so I cannot help you with that.
Secondly, you shouldn’t need all the lines I’ve quoted.
Simply using
Heat_Flux_z = qw_z.collapse()
is sufficient.
The reason why I am using
qw_z = qw_z.collapse()
V = FunctionSpace(mesh, ("Discontinuous Lagrange", 0))
Heat_Flux_z = Function(V, name='Heat_Flux_z')
Heat_Flux_z.x.array[:] = qw_z.x.array
is because I wanted the name Heat_Flux_z to appear in Paraview instead of f. But it’s simply cosmetic.
I am simulating heat spreading through a multilayered and multi-materials model imported from Gmsh. Regarding the plotting, I have first a function that extracts the temperature values_T and heat flux values_q density distribution and its maximum per layer in the domain (each layer has a different size and thickness):
def extract_values(T, qw_z, boundary, mesh, cell_tags, name_layers, PM_dict):
Tfs = T.function_space
coordinates = Tfs.tabulate_dof_coordinates()
eps = 1e-8
ln = name_layers.index(boundary)
H_list = np.array([PM_dict[l]['thickness'][0] for l in name_layers], dtype=float)
HL = H_list[:ln].sum()
BB = PM_dict[boundary]['BoundingBox']
X_list = np.array( [ np.array([bb.x1, bb.x2]).flatten() for bb in BB ] )
Y_list = np.array( [ np.array([bb.y1, bb.y2]).flatten() for bb in BB ] )
Xm_list = np.array( [ (bb.x1+bb.x2)/2. for bb in BB ] )
Ym_list = np.array( [ (bb.y1+bb.y2)/2. for bb in BB ] )
DX_list = np.array( [ (bb.x2-bb.x1) for bb in BB ] )
DY_list = np.array( [ (bb.y2-bb.y1) for bb in BB ] )
Xmin = X_list.flatten().min()
Xmax = X_list.flatten().max()
Ymin = Y_list.flatten().min()
Ymax = Y_list.flatten().max()
dx = 0.05 * Xmax
dy = 0.05 * Ymax
mXmin = coordinates[:,0] > Xmin + dx
mXmax = coordinates[:,0] < Xmax - dx
mYmin = coordinates[:,1] > Ymin + dy
mYmax = coordinates[:,1] < Ymax -dy
mZmax = abs(coordinates[:,2] - HL)<eps
if boundary == 'chip':
print("Chip level reached.")
HL = HL + H_list[-1]
if len(BB)==1:
print("BOUNDARY BOX USED")
d = np.where(mXmin & mXmax & mYmin & mYmax & mZmax)
else:
print("BOUNDARY BOX NOT USED")
d = np.where(mZmax)
xyz = coordinates[d]
values_q = qw_z.x.array[d]
values_T = T.x.array[d]
x = np.array([v[0] for v in xyz])
y = np.array([v[1] for v in xyz])
tdim = mesh.topology.dim
num_cells = mesh.topology.index_map(tdim).size_local
mesh = dolfinx.cpp.mesh.Mesh(mesh.comm, mesh.topology, mesh.geometry) # TODO remove with dolfinx 0.7.0
hmin = dolfinx.cpp.mesh.h(mesh, tdim, np.arange(num_cells, dtype=np.int32)).min()
Tmax = values_T.max()
imax = np.where(values_T==Tmax)
dm = np.where(values_T==Tmax)
qmax = values_q[dm].mean()
qavg = values_q.mean()
Tmax = values_T[dm].mean()
if 'PM_results' in PM_dict: # heater
dct = PM_dict['PM_results']
else: # Module
dct = PM_dict
if 'chip' in dct:
BB = dct['chip']['BoundingBox']
Xm_list = np.array( [ (bb.x1+bb.x2)/2. for bb in BB ] )
Ym_list = np.array( [ (bb.y1+bb.y2)/2. for bb in BB ] )
mXYc = (x!=x) # Initialize at False
for xc,yc in zip(Xm_list, Ym_list):
min_R = min(np.sqrt((x-xc)**2+(y-yc)**2))
mXYc = mXYc | (np.sqrt((x-xc)**2+(y-yc)**2)<=min_R)
ic = np.where(mXYc)
Tc_list = values_T[ic]
else:
Tc_list = Tmax
return x, y, values_q, values_T, qmax, qavg, Tmax, imax, Tc_list
Then the heat flux density map with its maximum is plotted for each layer with
def plot_map(x0, y0, f, imax, layer, name_layers ,PM_dict, saveDir, fn, unit=''):
print("Plot layer = ", layer)
# Map centering offsets
Dx = Dy = 0. #Chips and solder chip
if layer in PM_dict:
Dx = PM_dict[layer]['W'][0] / 2.
Dy = PM_dict[layer]['L'][0] / 2.
BB = PM_dict[layer]['BoundingBox']
ln = name_layers.index(layer)
# Center and convert to mm
x = (x0 + Dx) * 1e3
y = (y0 + Dy) * 1e3
xmax = max(x)
xmin = min(x)
ymax = max(y)
ymin = min(y)
xm = x[imax]
ym = y[imax]
fmax = f[imax]
fig = plt.figure()
for BBi in BB:
X_min = BBi.x1
X_max = BBi.x2
Y_min = BBi.y1
Y_max = BBi.y2
dx = 0
dy = 0
mX_min = (x0 >= (X_min+dx))
mX_max = (x0 <= (X_max-dx))
mY_min = (y0 >= (Y_min+dy))
mY_max = (y0 <= (Y_max-dy))
xbb = x[np.where(mX_min & mX_max & mY_min & mY_max)]
ybb = y[np.where(mX_min & mX_max & mY_min & mY_max)]
fbb = f[np.where(mX_min & mX_max & mY_min & mY_max)]
xbbmin = xbb.min()
xbbmax = xbb.max()
ybbmin = ybb.min()
ybbmax = ybb.max()
levels = np.linspace(min(fbb), max(fbb), 256)
points = (xbb, ybb)
grid_x, grid_y = np.mgrid[xbbmin:xbbmax:100j, ybbmin:ybbmax:100j]
grid = scipy.interpolate.griddata(points, fbb, (grid_x, grid_y), method='linear', fill_value=min(fbb))
plt.imshow(grid.T, extent=(xbbmin,xbbmax,ybbmin,ybbmax), origin='lower', cmap=plt.cm.jet)
plt.xticks([xbbmin,xbbmax])
plt.yticks([ybbmin,ybbmax])
cbar = plt.colorbar(pad=0.05)
cbar.set_label(unit)
plt.xlabel('x (mm)')
plt.ylabel('y (mm)')
plt.axis([xmin, xmax, ymin, ymax])
plt.axis('scaled')
plt.plot(xm, ym, 'w+', markersize=8)
plt.annotate(int(np.round(fmax,0)), (xm, ym), xycoords='data', xytext=(3,3), textcoords='offset points', size=7, color='w', **kwargs)
plt.savefig(saveDir+fn+'.pdf', format='pdf')
plt.close()
return 0
It’s far from being a minimal working example though.
As you state
This code is really hard to parse for anyone not familiar with the details of what you are doing.
Have you looked at the solution in Paraview? Then you should see that it is correct when using DG-0 (save Heat_Flux_z to XDMFFile)
Using DG-0, the Paraview solution look correct.
Furthermore if I test different ways of obtaining coordinates, I get
using temperature variable T
coordinates = T.function_space.tabulate_dof_coordinates()
coordinates = [[-2.25612395e-02 4.39310798e-02 3.74000000e-03]
[-2.06470588e-02 4.65000000e-02 3.74000000e-03]
[-2.04751368e-02 4.35168462e-02 3.74000000e-03]
...
[ 5.65000000e-02 -6.25000000e-02 2.50000000e-04]
[ 5.15869565e-02 -6.25000000e-02 1.00000000e-03]
[ 5.65000000e-02 -6.25000000e-02 1.38777878e-20]]
Using qw_z and DG0
coordinates = qw_z.function_space.tabulate_dof_coordinates()
coordinates = [[-0.02103964 0.04436619 0.00370875]
[-0.02000428 0.04439077 0.00370875]
[-0.02124815 0.04313583 0.00370875]
...
[ 0.05326228 -0.06173451 0.0008125 ]
[ 0.05449054 -0.06173451 0.000125 ]
[ 0.05449054 -0.06173451 0.000375 ]]
Using qw_z and CG1
coordinates = qw_z.function_space.tabulate_dof_coordinates()
coordinates = [[-2.25612395e-02 4.39310798e-02 3.74000000e-03]
[-2.06470588e-02 4.65000000e-02 3.74000000e-03]
[-2.04751368e-02 4.35168462e-02 3.74000000e-03]
...
[ 5.65000000e-02 -6.25000000e-02 2.50000000e-04]
[ 5.15869565e-02 -6.25000000e-02 1.00000000e-03]
[ 5.65000000e-02 -6.25000000e-02 1.38777878e-20]]
It would mean that the coordinates are ordered identically to the temperature variable with qw_z+CG1 but differently with qw_z+DG0. I guess that’s the issue.
Of course they are ordered differently. A first order Lagrange space has degrees of freedom at every vertex of the mesh. A DG 0 space consists of piecewise constants, they have dof coordinates at the midpoint of each cell.
Then it means that I cannot retrieve values at the interface between 2 layers with DG0 but only in the middle of the last cell before the interface?
Yes, because the value is not well defined on an interface, at it has different values on each side (as it is constant in each cell)
There is only one thing I still struggle to understand, is why I should not use CG1 in Expression(k*grad(T), V.element.interpolation_points())considering that the problem variable T was defined with
V = dolfinx.fem.FunctionSpace(mesh, ('CG', 1))
T = ufl.TrialFunction(V)
I used DG0 though to define a constant thermal conductivity in each layer.
I have finished the transition from Fenics legacy to Fenicsx and I must say that it is more pythonic, the interaction with numpy is lean and the mesh importation from gmsh straighforward now. It’s definitely worth the effort! Thanks a lot for your help!
When you differentiate a piecewise linear continuous function, you end you with cellwise constant basis functions. See for instance Introduction to the Finite Element Method (FEM) — Simula Summer School in Computational Physiology
I stumbled on this topic Defining subdomains from gmsh on CG elements - #2 by dokken where you say that " 1. Interpolate (not well defined) or project (well defined) the DG-0 data into a CG-1 space."
Does it mean that I can compute my heat flux qw safely in a CG1 space instead of DG0 using project?
Using https://github.com/michalhabera/dolfiny/blob/master/dolfiny/projection.py to make my own project function to perform qw = project(k*grad(T), V) I get the results below
interpolate: