Problems with setting vector parameter

dolfinx
1

Hi! I just learned FEniCS for a short time. I am trying to calculate the magnetic field produced by a permanent magnet in vacuum. Since no current is involved, I chose to calculate the magnetic scalar potential to solve this problem.The partial differential equations satisfy Poisson’s equation or Laplace’s equation,
\nabla\cdot\nabla\Phi=\nabla\cdot \vec{M}
I followed the Poisson Problem tutorial to convert it to weak form,
\int_\Omega {\nabla \Phi \cdot \nabla vdx}=\int_\Omega {\nabla\cdot \vec{M}dx}+\int_{\partial\Omega}{{\partial\Phi} \over {\partial n}}vds
What confuses me is how to set vector parameter in FEniCS, which in my case refers to the magnetization vector \vec{M} .The following code is still not working. And, is it possible to select specific nodes to set the vector parameter？

import math
import numpy as np
from dolfinx import *
from dolfinx.mesh import create_unit_square, locate_entities, locate_entities_boundary, meshtags
from mpi4py import MPI
from petsc4py.PETSc import ScalarType
from ufl import (FacetNormal, Measure, SpatialCoordinate, TestFunction, TrialFunction, 
                 div, dot, dx, grad, inner, lhs, rhs)
#creat mesh
mesh = mesh.create_unit_square(MPI.COMM_WORLD, 20, 20)

V = fem.FunctionSpace(mesh,("DG",1))
#Q = VectorFunctionSpace(mesh, ("DG", 2))

u = TrialFunction(V)
v = TestFunction(V)

u_0 = 4*math.pi*1e-7
Ms = 1e4/0.012566

Mv = fem.Function(V)
def Omega_m(x):
    return np.logical_and(np.logical_and(x[0]>=0.25, x[0]<=0.75),
                          np.logical_and(x[1]>=0.25, x[1]<=0.75))
def Omega_v(x):
    return np.logical_or(np.logical_or(x[0]<=0.25, x[0]>=0.75),
                         np.logical_or(x[1]<=0.25, x[1]>=0.75))
def D_bc(x):
    return np.logical_or(np.logical_or(np.isclose(x[0], 0), np.isclose(x[0], 1)),
                         np.logical_or(np.isclose(x[1], 0), np.isclose(x[1], 1)))

cells_m = locate_entities(mesh, mesh.topology.dim, Omega_m)
cells_v = locate_entities(mesh, mesh.topology.dim, Omega_v)

#magnetization vector for magnet and vacuum，I know it is not a vector yet.
Mv.x.array[cells_m] = np.full_like(cells_m, Ms, dtype=ScalarType)
Mv.x.array[cells_v] = np.full_like(cells_v, 0, dtype=ScalarType)

dofs = fem.locate_dofs_geometrical(V, D_bc)
bcs = [fem.dirichletbc(ScalarType(0), dofs, V)]

dx = Measure('dx', mesh)
a = inner(grad(u), grad(v)) * dx
f = div(Mv)
L = f * v * dx

problem = fem.petsc.LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()

Looking forward to everyone’s help, many thanks！

2

There are a few ways to go about this. The simplest way is to use interpolation over subsets of cells:

import numpy as np
from mpi4py import MPI
import dolfinx.io
import dolfinx.mesh
import dolfinx.fem
mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 20, 20)

V = dolfinx.fem.FunctionSpace(mesh, ("DG", 1))


Ms = 1e4/0.012566


def Omega_m(x):
    return np.logical_and(np.logical_and(x[0] >= 0.25, x[0] <= 0.75),
                          np.logical_and(x[1] >= 0.25, x[1] <= 0.75))


def Omega_v(x):
    return np.logical_or(np.logical_or(x[0] <= 0.25, x[0] >= 0.75),
                         np.logical_or(x[1] <= 0.25, x[1] >= 0.75))


def D_bc(x):
    return np.logical_or(np.logical_or(np.isclose(x[0], 0), np.isclose(x[0], 1)),
                         np.logical_or(np.isclose(x[1], 0), np.isclose(x[1], 1)))


cells_m = dolfinx.mesh.locate_entities(mesh, mesh.topology.dim, Omega_m)
cells_v = dolfinx.mesh.locate_entities(mesh, mesh.topology.dim, Omega_v)


def f_m(x):
    return np.full((2, x.shape[1]), Ms)


def f_v(x):
    return np.full((2, x.shape[1]), 0)


W = dolfinx.fem.VectorFunctionSpace(mesh, ("DG", 0))
M = dolfinx.fem.Function(W)
M.interpolate(f_m, cells=cells_m)
M.interpolate(f_v, cells=cells_v)
M.x.scatter_forward()
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "M.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_function(M)
3

Thanks a lot for your reply!
However, there is a problem when I try to solve it. Is this due to the multiple function spaces? How can I make the program work as expected variational form？ \int_\Omega{\nabla\Phi\cdot\nabla v}dx=\int_{\Omega}{\nabla\cdot\vec{M}}dx 

import numpy as np
from mpi4py import MPI
import dolfinx.io
import dolfinx.mesh
import dolfinx.fem
from petsc4py.PETSc import ScalarType
from ufl import (FacetNormal, Measure, SpatialCoordinate, TestFunction, TrialFunction, 
                 div, dot, dx, grad, inner, lhs, rhs, nabla_div)

mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 20, 20)

V = dolfinx.fem.FunctionSpace(mesh, ("DG", 1))
u = TrialFunction(V)
v = TestFunction(V)

Ms = 1e4/0.012566


def Omega_m(x):
    return np.logical_and(np.logical_and(x[0] >= 0.25, x[0] <= 0.75),
                          np.logical_and(x[1] >= 0.25, x[1] <= 0.75))


def Omega_v(x):
    return np.logical_or(np.logical_or(x[0] <= 0.25, x[0] >= 0.75),
                         np.logical_or(x[1] <= 0.25, x[1] >= 0.75))


def D_bc(x):
    return np.logical_or(np.logical_or(np.isclose(x[0], 0), np.isclose(x[0], 1)),
                         np.logical_or(np.isclose(x[1], 0), np.isclose(x[1], 1)))


cells_m = dolfinx.mesh.locate_entities(mesh, mesh.topology.dim, Omega_m)
cells_v = dolfinx.mesh.locate_entities(mesh, mesh.topology.dim, Omega_v)


def f_m(x):
    return (np.ones((x.shape[1],2))*np.array([0, Ms])).T


def f_v(x):
    return np.full((2, x.shape[1]), 0)


W = dolfinx.fem.VectorFunctionSpace(mesh, ("DG", 0))
M = dolfinx.fem.Function(W)
M.interpolate(f_m, cells=cells_m)
M.interpolate(f_v, cells=cells_v)
M.x.scatter_forward()

dofs = dolfinx.fem.locate_dofs_geometrical(V, D_bc)
bcs = [dolfinx.fem.dirichletbc(ScalarType(0), dofs, V)]

dx = Measure('dx', mesh)
a = inner(grad(u), grad(v)) * dx
f = div(M)
L = f*v*dx

problem = dolfinx.fem.petsc.LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
Input In [1], in <cell line: 60>()
     57 f = div(M)
     58 L = f*v*dx
---> 60 problem = dolfinx.fem.petsc.LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
     61 uh = problem.solve()

File /usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/petsc.py:565, in LinearProblem.__init__(self, a, L, bcs, u, petsc_options, form_compiler_params, jit_params)
    562 self._A = create_matrix(self._a)
    564 self._L = _create_form(L, form_compiler_params=form_compiler_params, jit_params=jit_params)
--> 565 self._b = create_vector(self._L)
    567 if u is None:
    568     # Extract function space from TrialFunction (which is at the
    569     # end of the argument list as it is numbered as 1, while the
    570     # Test function is numbered as 0)
    571     self.u = _Function(a.arguments()[-1].ufl_function_space())

File /usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/petsc.py:79, in create_vector(L)
     69 def create_vector(L: form_types) -> PETSc.Vec:
     70     """Create a PETSc vector that is compaible with a linear form.
     71 
     72     Args:
   (...)
     77 
     78     """
---> 79     dofmap = L.function_spaces[0].dofmap
     80     return la.create_petsc_vector(dofmap.index_map, dofmap.index_map_bs)

IndexError: list index out of range
4

Forgive me for one more question. When I consider the Neumann boundary conditions at the interface between the material and vaccum, as given in the variational form earlier, it seems to have no effect and the results are all 0.
\int_\Omega{\nabla\Phi\cdot\nabla v}dx=\int_{\partial\Omega}{\nabla\cdot\vec{M}}dx+\int_{\Gamma m}{\vec{n}\cdot\vec{M}}ds 

import numpy as np
import pyvista

from dolfinx.fem import (Constant, Function, FunctionSpace, 
                         assemble_scalar, dirichletbc, form, locate_dofs_geometrical, locate_dofs_topological,Expression)
from dolfinx.fem.petsc import LinearProblem
from dolfinx.mesh import create_unit_square, locate_entities, locate_entities_boundary, meshtags
from dolfinx.io import XDMFFile
from mpi4py import MPI
from petsc4py.PETSc import ScalarType
from ufl import (FacetNormal, Measure, SpatialCoordinate, TestFunction, TrialFunction, 
                 div, dot, dx, grad, inner, lhs, rhs)
#crear mesh
mesh = create_unit_square(MPI.COMM_WORLD, 20, 20)

l1 = 1  #length of domain
l2 = 0.5   #length of magnet
l3 = (l1-l2)/2

#mark boundaries
x = SpatialCoordinate(mesh)
boundaries = [(1, lambda x: np.logical_and(np.isclose(x[1], l3),    np.logical_and(l3 <= x[0], x[0] <= l3+l2))),
              (2, lambda x: np.logical_and(np.isclose(x[1], l3+l2), np.logical_and(l3 <= x[0], x[0] <= l3+l2))),
              (3, lambda x: np.logical_and(np.isclose(x[0], l3),    np.logical_and(l3 <= x[1], x[1] <= l3+l2))),
              (4, lambda x: np.logical_and(np.isclose(x[0], l3+l2), np.logical_and(l3 <= x[1], x[1] <= l3+l2))),
              (5, lambda x: np.logical_or(np.logical_or(np.isclose(x[0], 0), np.isclose(x[0], l1)),
                                          np.logical_or(np.isclose(x[1], 0), np.isclose(x[1], l1))))]

facet_indices, facet_markers = [], []
fdim = mesh.topology.dim - 1
for (marker, locator) in boundaries:
    facets = locate_entities(mesh, fdim, locator)
    facet_indices.append(facets)
    facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = meshtags(mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])

V = FunctionSpace(mesh, ("CG", 1))
u = TrialFunction(V)
v = TestFunction(V)
ds = Measure("ds", domain=mesh, subdomain_data=facet_tag)

#set boundary condition
class BoundaryCondition():
    def __init__(self, type, marker, values):
        self._type = type
        if type == "Dirichlet":
            u_D = Function(V)
            u_D.interpolate(values)
            facets = facet_tag.find(marker)
            dofs = locate_dofs_topological(V, fdim, facets)
            self._bc = dirichletbc(u_D, dofs)
        elif type == "Neumann":
                self._bc = inner(values, v) * ds(marker)
        elif type == "Robin":
            self._bc = values[0] * inner(u-values[1], v)* ds(marker)
        else:
            raise TypeError("Unknown boundary condition: {0:s}".format(type))
    @property
    def bc(self):
        return self._bc

    @property
    def type(self):
        return self._type

Ms = 1
u_ex =  lambda x: Ms + x[0]*0
u_ex2 = lambda x: -Ms + x[0]*0
u_ex3 = lambda x: 0+x[0]*0
g1 = u_ex(x)
g2 = u_ex2(x)
g3 = u_ex3(x)

boundary_conditions = [BoundaryCondition("Neumann", 1, g1),
                       BoundaryCondition("Neumann", 2, g2),
                       BoundaryCondition("Neumann", 3, g3),
                       BoundaryCondition("Neumann", 4, g3),
                       BoundaryCondition("Dirichlet", 5, u_ex3)]

dx = Measure('dx', mesh)
a = inner(grad(u), grad(v)) * dx
f = Constant(mesh, ScalarType(0))
L = f*v*dx

#add boundary condition
bcs = []
for condition in boundary_conditions:
    if condition.type == "Dirichlet":
        bcs.append(condition.bc)
    else:
        L += condition.bc

#solve problem
problem = LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()