Imposing a vector source - how to?

General
1

Dear all,

I am following up on building a case study for 3D magnetostatic problem involving a coil with an impressed current density. I am stuck on imposing the current density. it runs but it appears that I do not a good job in setting the current density right.

I have built a small 2D case to try out. I am trying to impose a “current density” J along one axis in this specific case scenario.

I am purposely using gmsh to build the mesh to make sure that the mesh is correctly read and handle on my side.

Here is the code:

from mpi4py import MPI
from dolfinx.io import gmshio
from dolfinx.io import XDMFFile
from dolfinx.fem import (Expression, Function, FunctionSpace, VectorFunctionSpace)
from dolfinx.fem import (dirichletbc, locate_dofs_topological)
from dolfinx import *

import ufl
from ufl import as_vector, dot, dx, grad, curl, inner

from petsc4py.PETSc import ScalarType
import numpy as np

import gmsh
import sys

rank = MPI.COMM_WORLD.rank

gmsh.initialize()

gmsh.model.add("exampleGmsh")

model_rank = 0
mesh_comm = MPI.COMM_WORLD
gdim = 2

di = 0.2
do = 2*di

lc = di / 10

if mesh_comm.rank == model_rank:
	factory = gmsh.model.geo

	factory.addPoint(-di, -di, 0, lc, 1)
	factory.addPoint(di, -di, 0, lc, 2)
	factory.addPoint(di, di, 0, lc, 3)
	factory.addPoint(-di, di, 0, lc, 4)
	factory.addLine(1, 2, 1)
	factory.addLine(2, 3, 2)
	factory.addLine(3, 4, 3)
	factory.addLine(4, 1, 4)
	factory.addCurveLoop([1, 2, 3, 4], 1)
	factory.addPlaneSurface([1], 1)

#	factory.mesh.setTransfiniteCurve(1, 10)
#	factory.mesh.setTransfiniteCurve(2, 10)
#	factory.mesh.setTransfiniteCurve(3, 10)
#	factory.mesh.setTransfiniteCurve(4, 10)
#	factory.mesh.setTransfiniteSurface(1, "Left", [1, 2, 3, 4])
#	factory.mesh.setRecombine(2, 1)

	factory.addPoint(-do, -do, 0, lc, 5)
	factory.addPoint(do, -do, 0, lc, 6)
	factory.addPoint(do, do, 0, lc, 7)
	factory.addPoint(-do, do, 0, lc, 8)
	factory.addLine(5, 6, 5)
	factory.addLine(6, 7, 6)
	factory.addLine(7, 8, 7)
	factory.addLine(8, 5, 8)
	factory.addCurveLoop([5, 6, 7, 8], 2)
	factory.addPlaneSurface([1, 2], 2)

	factory.removeAllDuplicates

	factory.synchronize()

	gmsh.model.addPhysicalGroup(gdim, [1], 1, name="Inner")
	gmsh.model.addPhysicalGroup(gdim, [2], 2, name="Outer")
	gmsh.model.addPhysicalGroup(gdim-1, [5, 6, 7, 8], 3, name="Boundary")

	gmsh.model.mesh.generate(gdim)

	gmsh.write("exampleGmsh.msh")

# Read in by dolfinx
mesh, cell_tags, facet_tags = gmshio.model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=gdim)
gmsh.finalize()

# Output DOLFINx meshes to file
with XDMFFile(MPI.COMM_WORLD, "meshCheck.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_meshtags(facet_tags)
    xdmf.write_meshtags(cell_tags)

# Function space: edge elements
V0 = FunctionSpace(mesh, ("N1curl", 1))

# Create facet to cell connectivity required to determine boundary facets
dofs_facets = locate_dofs_topological(V0, mesh.topology.dim-1, facet_tags.find(3))

# Working boundary condition
u0 = Function(V0)
u0.x.set(0.)
bc_working = dirichletbc(u0, dofs_facets)

# Trial and test functions
u = ufl.TrialFunction(V0)
v = ufl.TestFunction(V0)

##Define Current Density and Mu0
V1 = VectorFunctionSpace(mesh, ("CG", 1))
J = Function(V1, dtype=ScalarType, name='Current density')
cells1 = cell_tags.find(1)
J.interpolate(lambda x: np.vstack((0., 1.)), cells1)
#J.interpolate(lambda x: np.vstack((1e8*np.sin(np.arctan2(x[1],x[0])), 1e8*np.cos(np.arctan2(x[1],x[0])))), cells1)

#x = ufl.SpatialCoordinate(mesh)
#J = Expression(("J0*sin(atan2(x[1],x[0]))","J0*cos(atan2(x[1],x[0]))", "0.0"), J0 = 1e8)

#mu = 4*np.pi*10**-7
#J = 1e8

#a = (1 / mu) * inner(curl(u), curl(v)) * dx
#L = inner(J, v) * dx(subdomain_id=(1))

#A = Function(V0)
#problem = LinearProblem(a, L, u=A, bcs=[bc])
#problem.solve()

from dolfinx import io
with io.XDMFFile(mesh.comm, "J.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_meshtags(cell_tags) 
    xdmf.write_function(J)

Best,

Frederic

2

Consider:

J.interpolate(lambda x: np.vstack(
    [np.zeros(x.shape[1]), np.ones(x.shape[1])]), cells1)
J.x.scatter_forward()

The function you send in to interpolate takes in all coordinates x, in the shape of (3, num_points) at once, making it possible to call this function in a vectorized fashion. Thus, you need to return something of shape (2, num_points).

3

Thanks a lot for the help. I am moving forward to the 3D model and see how it goes.

Hopefully, the final product will be a contribution to the fenicsx community as an additional example.

Best,

Frederic

1 Like