Defining an expression on mesh edges

mesh
1

I am comparing simulation results for open cell foams obtained using the cell method to results achieved with FEM. For the FEM simulations, I use triagulations of the cell complexes used in the cell method simulations where the vertices of the cell complexes correspond to the metal in the foam (i.e. the metal parts of the foam are found on certain edges of my FEM mesh).

I would like to define the material parameters around the edges which correspond to the metal in my mesh (i.e. have an expression where if the degree of freedom is on an edge where there is metal, the parameter takes one value, and if it is not on an edge where there is metal, it takes another value).

I am starting with a 2D example with heat conduction only for proof of concept. I begin with a .geo file that defines my mesh and identifies my boundaries and the edges that correspond to the metal. I use gmsh and meshio to generate a .msh file and then two .xdmf files, one for the mesh and one for the boundaries and the metal edges.

Here is the .geo file for the dummy mesh I am using to get the implementation right

Point(1) = {0, 0, 0, 0.5};
Point(2) = {0.25, 0, 0, 0.5};
Point(3) = {0.5, 0, 0, 0.5};
Point(4) = {0., 0.25, 0, 0.5};
Point(5) = {0.25, 0.25, 0, 0.5};
Point(6) = {0.5, 0.25, 0, 0.5};
Point(7) = {0., 0.5, 0, 0.5};
Point(8) = {0.25, 0.5, 0, 0.5};
Point(9) = {0.5, 0.5, 0, 0.5};

Line(1) = {1, 4};
Line(2) = {4, 5};
Line(3) = {5, 1};
Line(4) = {5, 2};
Line(5) = {2, 1};
Line(6) = {5, 6};
Line(7) = {6, 2};
Line(8) = {6, 3};
Line(9) = {3, 2};
Line(10) = {4, 7};
Line(11) = {7, 8};
Line(12) = {8, 4};
Line(13) = {8, 5};
Line(14) = {8, 9};
Line(15) = {9, 5};
Line(16) = {9, 6};

Curve Loop(1) = {11, 12, 10};
Curve Loop(2) = {1, 2, 3};
Surface(1) = {2};

Curve Loop(3) = {3, -5, -4};
Surface(2) = {3};
Curve Loop(4) = {4, -7, -6};
Surface(3) = {4};
Curve Loop(5) = {7, -9, -8};
Surface(4) = {5};
Surface(5) = {1};
Curve Loop(6) = {12, 2, -13};
Surface(6) = {6};
Curve Loop(7) = {13, -15, -14};
Surface(7) = {7};
Curve Loop(8) = {15, 6, -16};
Surface(8) = {8};

Physical Surface("Interior", 1) = {1, 2, 3, 4, 5, 6, 7, 8};
Physical Curve("Cold", 2) = {1, 10};
Physical Curve("Hot", 3) = {8, 16};
Physical Curve("Adiabatic", 4) = {5, 9, 11, 14};
Physical Curve("Metal", 5) = {11, 14, 12, 13, 15, 6, 2, 16, 3, 4, 5, 1, 7, 9};

Here is my minimal code for the heat conduction after I have generated the .xdmf files.

# import
import fenics as fn
import matplotlib.pyplot as plt

# Define model dimensions
fullDim = 2 #2D model

# Import mesh
mesh = fn.Mesh()
mvc = fn.MeshValueCollection("size_t", mesh, fullDim)
with fn.XDMFFile("mesh.xdmf") as infile:
   infile.read(mesh)
   infile.read(mvc, "name_to_read")
mf = fn.cpp.mesh.MeshFunctionSizet(mesh, mvc)

# Import boundaries
mvc2 = fn.MeshValueCollection("size_t", mesh, fullDim-1)
with fn.XDMFFile("facet_mesh.xdmf") as infile:
   infile.read(mvc2, "name_to_read")
cf = fn.cpp.mesh.MeshFunctionSizet(mesh, mvc2)

# Define function space
V = fn.FunctionSpace(mesh, 'Lagrange', 1)
dx = fn.Measure('dx', domain=mesh, subdomain_data=mf)

# Trial function
u = fn.TrialFunction(V)

# Test function
v = fn.TestFunction(V)

# Define boundary condition
bc = fn.DirichletBC(V, fn.Constant(500), cf, 2)
bcs = [bc]
bc = fn.DirichletBC(V, fn.Constant(600), cf, 3)
bcs.append(bc)

# Define constants
a_c = fn.Constant(1)

# some plot setup
plt.figure()

# equation solving
F = u*v*dx + a_c*fn.dot(fn.grad(u), fn.grad(v))*dx
a, L = fn.lhs(F), fn.rhs(F)
u = fn.Function(V)
fn.solve(a == L, u, bcs)

c = fn.plot(u, title='Temperature')
plt.colorbar(c)
plt.show()

How can I redefine a_c in terms of an Expression so that the degrees of freedom on the edges defined by Physical Curve("Metal", 5) in my .geo file have a value which is different than the degrees of freedom not on these edges?

Any help you can offer would be greatly appreciated.

1 Like
2

Consider the following minimal example using a built in mesh, where we assign a value to the line x=0.5 that goes through the mesh:

from dolfin import (Constant, DirichletBC, File, Function, FunctionSpace,
                    MeshFunction, SubDomain, UnitSquareMesh, near)

mesh = UnitSquareMesh(10, 10)

mf = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)


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


interface = InternalInterface()
interface.mark(mf, 1)


V = FunctionSpace(mesh, "CG", 1)
u = Function(V)
u.vector()[:] = 2

a_c = Constant(5)
bc = DirichletBC(V, a_c, mf, 1)
bc.apply(u.vector())

File("u.pvd") << u
1 Like
3

Thanks so much @dokken, this works perfectly!