Is there a method for integrating a field along a line? For example, is it possible to compute the integral for Ampere’s Law in 3D?
I = \oint \mathbf{H}\cdot d\mathbf{l}
The path can lie on a boundary or be just inside the volumetric region, the idea being the calculation of current in a conductor immersed in an electromagnetic field.
Cheers!
Hello!
You can integrate along arbitrary lines using ufl.Measure.
The main idea is to mark your line in Gmsh using Physical Line and then use the line ID in a custom measure.
For example, in 2D I want to integrate constant function u(x,y)=1 along the red line.
For this I create a custom measure dS where I specify MeshTags for it to use
dS = ufl.Measure('dS', domain=mesh, subdomain_data=boundaries)
After that I could select appropriate line by supplying its ID as a measure argument. In this example red line’s ID is 2.
print(dolfinx.fem.assemble_scalar(u * dS(2)))
This particular line computes \int_L\operatorname{u}dS. You could use this measure in other forms as well.
Gmsh script for this example (Conversion between .msh and .xdmf is described here.):
Geometry.OCCTargetUnit = "M";
r1 = 0.5;
r2 = 1;
dx1 = 0;
dy1 = 0;
dx2 = 0;
dy2 = 0;
mesh1 = 0.04;
mesh2 = 0.04;
p0 = newp; Point(p0) = {dx1,dy1,0,mesh1};
p1 = newp; Point(p1) = {r1+dx1,dy1,0,mesh1};
p2 = newp; Point(p2) = {dx1,r1+dy1,0,mesh1};
p3 = newp; Point(p3) = {-r1+dx1,dy1,0,mesh1};
p4 = newp; Point(p4) = {dx1,-r1+dy1,0,mesh1};
l1 = newl; Circle(l1) = {p1,p0,p2};
l2 = newl; Circle(l2) = {p2,p0,p3};
l3 = newl; Circle(l3) = {p3,p0,p4};
l4 = newl; Circle(l4) = {p4,p0,p1};
ll1 = newll; Line loop(ll1) = {l1,l2,l3,l4};
s1 = news; Surface(s1) = {ll1};
p00 = newp; Point(p00) = {dx2,dy2,0,mesh2};
p10 = newp; Point(p10) = {r2+dx2,dy2,0,mesh2};
p20 = newp; Point(p20) = {dx2,r2+dy2,0,mesh2};
p30 = newp; Point(p30) = {-r2+dx2,dy2,0,mesh2};
p40 = newp; Point(p40) = {dx2,-r2+dy2,0,mesh2};
l10 = newl; Circle(l10) = {p10,p00,p20};
l20 = newl; Circle(l20) = {p20,p00,p30};
l30 = newl; Circle(l30) = {p30,p00,p40};
l40 = newl; Circle(l40) = {p40,p00,p10};
ll10 = newll; Line loop(ll10) = {l10,l20,l30,l40};
s10 = news; Surface(s10) = {ll10, ll1};
// Assign IDs to the mesh regions
Physical Surface(1) = {s10,s1};
Physical Line(2) = {l1,l2};
DolfinX code:
import dolfinx
import dolfinx.io
import ufl
from mpi4py import MPI
# Read mesh and meshtags
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "mesh.xdmf", "r") as xdmf:
mesh = xdmf.read_mesh(name="Grid")
mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "mesh_boundary.xdmf", "r") as xdmf:
boundaries = xdmf.read_meshtags(mesh, name="Grid")
V = dolfinx.FunctionSpace(mesh, ("Lagrange", 1))
# Set constant function
u = dolfinx.Function(V)
with u.vector.localForm() as loc:
loc.set(1)
# Create measure
dS = ufl.Measure('dS', domain=mesh, subdomain_data=boundaries)
# Integrate
print(dolfinx.fem.assemble_scalar(u * dS(2)))
Hi Community,
I have a similar problem: I am trying to calculate the length of an 1D line of a 3D geometry. Is there a way to do this? To the best of my knowledge, the integral types of ufl.Measure are restricted either to dx (3D) or ds/dS (2D), but there is no integral type for lines.
Here is the (non-functioning) MWE:
import numpy as np
from mpi4py import MPI
from dolfinx.fem import assemble_scalar, form
from ufl import Measure
import dolfinx
comm = MPI.COMM_WORLD
domain = dolfinx.mesh.create_unit_cube(comm, 10,10,10, cell_type=dolfinx.mesh.CellType.hexahedron)
def line_1D(x):
return np.logical_and(np.isclose(x[0], 0), np.isclose(x[1], 0))
lines = dolfinx.mesh.locate_entities_boundary(domain, 1, line_1D)
line_tag = dolfinx.mesh.meshtags(domain, 1, lines, 9)
dLine = Measure("ds", domain, subdomain_data=line_tag)
Length_line = domain.comm.allreduce(assemble_scalar(form(1 * dLine(9))), op=MPI.SUM)
I would create a sub-mesh of the facets you want to measure, so that you can use the dx measures on the sub mesh. Similar to: https://github.com/FEniCS/dolfinx/blob/88bf353382240b5c1a8826c34be79c142e6ab28a/python/test/unit/fem/test_assemble_submesh.py#L85-L104 just with the dim being 1 instead of 2
This works! Thank you.
I guess the right way (regarding the original post) would be to (i) create the submesh, (ii) interpolate the quantities to be integrated from the original mesh to the submesh and (iii) integrate the quantity over the submesh.