DOF indices misaligned & boundaries do not transfer heat

[DanielNobbe]

Hi all,

I'm pretty new to using this library, and having a hard time adapting one of the examples to my use case.

I have adapted example 17 (heat transfer in a wire) to a case where there are four quadrants in a square. See an outcome from my test below. The left bottom square is heated.

I run into two problem:

• In basis0.get_dofs(elements=s), the indices returned by self.normalize_elements do not align with the indices in self.topo.t. The returned values from normalize_elements start at index 160, and the indices in the last quadrant go 160 over the max index in self.topo.t. I am not sure where to look for the cause of this, so any pointers would be great.
• By adding a subtraction of 160 inside the get_element_dofs call, I can work around this problem, but I am not sure what the result is. In that case, my analysis runs, but there is no heat propagated through the active boundary between the bottom-left and top-left quadrants. I'm sure I'm missing something here, since there is heat transfer through the boundary layer in the example. Any pointers would be great.

Thanks in advance!

See my code below:

r"""
Based on example 17: https://github.com/kinnala/scikit-fem/blob/7.0.1/docs/examples/ex17.py

"""
from pathlib import Path
from typing import Optional

import numpy as np

from skfem import *
from skfem.helpers import dot, grad
from skfem.models.poisson import mass, unit_load
from skfem.io.meshio import from_meshio
from skfem.io.json import from_file

import pygmsh

### Generate mesh with pygmsh
with pygmsh.geo.Geometry() as geom:

   # Four quadrants
    bl = geom.add_rectangle(0.0, 1.0, 0.0, 1.0, 0.0, 0.1)  # xmin, xmax, ymin, ymax, z, mesh_size
    br = geom.add_rectangle(1.0, 2.0, 0.0, 1.0, 0.0, 0.1)
    tl = geom.add_rectangle(0.0, 1.0, 1.0, 2.0, 0.0, 0.1)
    tr = geom.add_rectangle(1.0, 2.0, 1.0, 2.0, 0.0, 0.1)

    # Define subdomains
    geom.add_physical(bl, label='bottom-left')
    geom.add_physical(br, label='bottom-right')
    geom.add_physical(tl, label='top-left')
    geom.add_physical(tr, label='top-right')

    # Define each quadrant's edges as boundaries
    geom.add_physical(bl.lines[2], label='bl-top')
    geom.add_physical(br.lines[2], label='br-top')
    geom.add_physical(tl.lines[1], label='tl-right')
    geom.add_physical(tr.lines[3], label='tr-left')
    geom.add_physical(bl.lines[1], label='bl-right')
    geom.add_physical(br.lines[3], label='br-left')
    geom.add_physical(tl.lines[0], label='tl-bottom')
    geom.add_physical(tr.lines[0], label='tr-bottom')

    # Define each outer edge as boundary
    geom.add_physical([br.lines[1], tr.lines[1]], label='right')
    geom.add_physical([tr.lines[2], tl.lines[2]], label='top')
    geom.add_physical([tl.lines[3], bl.lines[3]], label='left')
    geom.add_physical([bl.lines[0], br.lines[0]], label='bottom')

    # Generate mesh with pygmsh
    mesh = geom.generate_mesh(dim=2)
mesh = from_meshio(mesh)

# adapted from example 17
joule_heating = 5.
heat_transfer_coefficient = 7.
# Define conductivity per quadrant
thermal_conductivity = {'bottom-right': 1.,  'bottom-left': 10., 'top-left': 5., 'top-right': 10000.}

@BilinearForm
def conduction(u, v, w):
    return dot(w['conductivity'] * grad(u), grad(v))

convection = mass

element = ElementTriP1()
basis = Basis(mesh, element)

basis0 = basis.with_element(ElementTriP1())

conductivity = basis0.zeros()

for s in mesh.subdomains:
    conductivity[basis0.get_dofs(elements=s)] = thermal_conductivity[s]

L = asm(conduction, basis, conductivity=basis0.interpolate(conductivity))

facet_basis = FacetBasis(mesh, element, facets=[mesh.boundaries['bl-top']])
H = heat_transfer_coefficient * asm(convection, facet_basis)

core_basis = Basis(mesh, basis.elem, elements=mesh.subdomains['bottom-left'])
f = joule_heating * asm(unit_load, core_basis)

temperature = solve(L + H, f)

if __name__ == '__main__':

    from os.path import splitext
    from sys import argv
    from skfem.visuals.matplotlib import draw, plot, savefig

    ax = draw(mesh)
    plot(basis, temperature, ax=ax, colorbar=True)
    savefig(splitext(argv[0])[0] + '_solution.png')

requirements.txt

[kinnala]

Perhaps you have duplicate nodes there on the interface, with no connectivity from one quadrant to another?

[kinnala]

mesh.draw(node_numbering=True).show() will generate quite a mess but perhaps with zooming you can see if that's the case.

[DanielNobbe]

Looks like there are indeed duplicate nodes on the interface, thank you for the suggestion. Do you have any ideas on how to link these duplicate nodes?

[kinnala]

You should preferably create a proper mesh right at the start. I don't use pygmsh and, therefore, I don't know how the code must be modified.

There is a hack to remove duplicate nodes via mesh = mesh + mesh but this will, as a side effect, remove any tags attached to the elements or the facets of your mesh. This means that those top-left etc. tags will vanish.

[kinnala]

It is of course possible to add any missing tags again using MeshTri.with_subdomains and MeshTri.with_boundaries commands but I think it is perhaps a suboptimal solution and I'd rather try fixing the original mesh.

It is also possible to create meshes manually by combining square blocks:

mesh1 = MeshTri().refined(5)
mesh2 = MeshTri().translated((1.0, 0.0)).refined(5)
mesh = mesh1 + mesh2
mesh = mesh.with_subdomains({
   'left': lambda x: x[0] < 1.,
   'right': lambda x: x[0] > 1.,
})

[kinnala]

If you are only interested in this simple geometry, it is also possible to skip all this mesh generation and tagging effort and simply hard code the following form with the heat conductivity defined inside:

@BilinearForm
def conduction(u, v, w):
    hc = ((w.x[0] < 1) * (w.x[1] < 1) * 10 +
          (w.x[0] > 1) * (w.x[1] > 1) * 10000 +
          (w.x[0] < 1) * (w.x[1] > 1) * 5 +
          (w.x[0] > 1) * (w.x[1] < 1) * 1)
    return dot(hc * grad(u), grad(v))