Robin conditions

General
1

Hello everyone!
I am new to fenics and trying to improve my skills by trying to solve some basic problems. Currently I am trying to solve scattering problem on a cylinder. But I got stuck on a step with Robin boundary conditions. I dont understand how to implement it on outer boundary.

from fenics import *
import sympy as sym
import matplotlib.pyplot as plt
import numpy as np
from mshr import *
import math


def mu(u):
    "Return nonlinear coefficient"
    return 1 + u**2


x, y = sym.symbols('x[0], x[1]')
u = 1 + x + 2*y
f = - sym.diff(q(u)*sym.diff(u, x), x) - sym.diff(q(u)*sym.diff(u, y), y)
f = sym.simplify(f)
u_code = sym.printing.ccode(u)
f_code = sym.printing.ccode(f)
print('u =', u_code)
print('f =', f_code)




# Construct Mesh
r_in = 0.5
r_out = 1.5
r1 = Circle(Point(0.0, 0.0), r_in)
r2 = Circle(Point(0.0, 0.0), r_out)

# Define domain and resolution
domain = r2-r1
res = 64
# Create mesh and define function space
mesh = generate_mesh(domain, res)
plot(mesh)
plt.show()
print('Success mesh')

V = FunctionSpace(mesh, 'P', 1)


# Define boundary condition
u_D = Expression(u_code, degree=2)

def fixed_boundary(x,on_boundary):
    r = math.sqrt(x[0] * x[0] + x[1] * x[1])
    return on_boundary and (near(r,r_in,0.001))
def ABC(x,on_boundary):
    r = math.sqrt(x[0] * x[0] + x[1] * x[1])
    return on_boundary and near(r,r_out,0.001)

bc = DirichletBC(V, Constant(0), boundary)
print('Success DirichletBC')

The question is what should I put in func ABC, so that I apply that condition:
image

2

See for instance:
https://fenicsproject.org/pub/tutorial/html/._ftut1014.html
for Legacy DOLFIN and
https://jsdokken.com/dolfinx-tutorial/chapter3/robin_neumann_diriclet.html?highlight=robin
for DOLFINx

3

Thanks for useful links. Read everything for legacy, but applying new code involving weak formulation for my case I got error. I attach code and bug report below:

from fenics import *
import sympy as sym
import matplotlib.pyplot as plt
from mshr import *
import math
import numpy as np

def mu(u):
    "Return nonlinear coefficient"
    return 1+u**2

k_0 = 2*np.pi
eps = 1
beta = k_0*k_0*mu(0)*eps
E0 = 0.1

x, y = sym.symbols('x[0], x[1]')
u = 1 + x + 2*y
f = - sym.diff(mu(u)*sym.diff(u, x), x) - sym.diff(mu(u)*sym.diff(u, y), y)
f = sym.simplify(f)
u_code = sym.printing.ccode(u)
f_code = sym.printing.ccode(f)
print('u =', u_code)
print('f =', f_code)




# Construct Mesh
r_in = 0.5
r_out = 1.5
r1 = Circle(Point(0.0, 0.0), r_in)
r2 = Circle(Point(0.0, 0.0), r_out)

# Define domain and resolution
domain = r2-r1
res = 15
# Create mesh
mesh = generate_mesh(domain, res)
plot(mesh)
plt.show()
print('Success mesh')


class Boundary_in(SubDomain):
    def Dirichlet(x, on_boundary):
        r = math.sqrt(x[0] * x[0] + x[1] * x[1])
        return on_boundary and (near(r, r_in, 0.001))

class Boundary_out(SubDomain):
    def ABC(self,x, on_boundary):
        r = math.sqrt(x[0] * x[0] + x[1] * x[1])
        return on_boundary and near(r, r_out, 0.001)

# Mark boundaries
boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim() - 1, mesh.domains())
boundary_markers.set_all(3)
b_in = Boundary_in()
b_out = Boundary_out()
b_in.mark(boundary_markers, 0)
b_out.mark(boundary_markers, 1)

# Redefine boundary integration measure
ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)

# Define function space and trial and test functions
V = FunctionSpace(mesh, 'P', 1)

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


bc = DirichletBC(V, Constant(0), boundary_markers, 0)
print('Success DirichletBC')

#ABC
E_inc = E0*np.exp(-1j*k_0*x[0])
integrals_R_a = []
integrals_R_L = []

gamma = 1j*k_0 + 1/(2*r_out)
q = -1j*k_0*E_inc+gamma*E_inc

integrals_R_a.append(gamma*u*v*ds(1))
integrals_R_L.append(q*v*ds(1))

# Sum integrals to define variational problem
a = -1*dot(grad(u), grad(v))*dx + beta*dot(u,v)*dx - sum(integrals_R_a)
L = 0*f*v*dx - sum(integrals_R_L)

F = a-L

# Compute solution
solve(F == 0, u, bc)

# Plot solution
plot(u)
plt.show()

Bug:

  File "/home/apoh_regex/ptn/trial.py", line 69, in <module>
    bin.mark(boundary_markers, 0)
RuntimeError:

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to check whether point is inside subdomain.
*** Reason:  Function inside() not implemented by user.
*** Where:   This error was encountered inside SubDomain.cpp.
*** Process: 0
***
*** DOLFIN version: 2019.2.0.dev0
*** Git changeset:  unknown
*** -------------------------------------------------------------------------
4

This function should be called inside, and similarly for

5

Sorry, English is foreign language for me. What do you mean by inside?

6

Should be

class Boundary_in(SubDomain):
    def inside(x, on_boundary):

And the same for the Boundary_out class

7

Thanks. I have one more question, hope this one is the last. I looke through legacy documentation, and found out, that it doescnt work with complex functions. I understand that the only way to solve is to split up my equation in weak form into imaginary and real parts and then check absolute value and angle dependence. But as I am solving infinite problem I need to put Robin conditions on outgoing wave (scattered). So for total Ez component I will have E = Esc + Einc. And in Robin boundary condition I will have x variable. Is there any easy way to implement it (below I attach my by-hand derivation of boundary condition, which then inserts into weak form)? Also is there a way to project function on normal direction?

1
11080×268 28.1 KB

Thanks for your support!

8

I would strongly suggest using DOLFINx if you are new to FEniCS, especially as it has complex number support.

For a normal projection, simply define a facetnormal, and call

u = Function(V)
n = FacetNormal(mesh)
un = dot(u, n) * n

For tangential projections see Tangential gradient operator error: invalid ranks 1 and 1 in product - #2 by dokken

9

Thanks, but still, how could I deal with x-dependence in Robin BC. Putting it inside weak formulation symbolically causes error, cause LU solve doesnt know how to deal with it

10

with x-dependence i guess you mean that the Robin condition is dependent of the physical location in space? If so, you can use x=SpatialCoordinate(mesh) in your variational forms, where x[0] is the x coordinate, x[1] is y etc

11

Hello, I have followed your advice and now solving with dofinx version in docker. But now I have another problem. Since here we use not mshr, but pygmsh to create mesh, I come into some troubles. I followed one of your tutorials on mesh plotting (Mesh generation and conversion with GMSH and PYGMSH | Jørgen S. Dokken). I defined model and created mesh and wanted to create triangle mesh and write in in .xdmf format, so that I could use it further to define BCs and so on, but during running error appeared:

r_in = 0.5
r_out = 1.5
resolution = 0.04
c = [ 0, 0, 0 ]

# Initialize empty geometry using the build in kernel in GMSH
geometry = pygmsh.geo.Geometry()
# Fetch model we would like to add data to
model = geometry.__enter__()
# Add circle
circ_inner = model.add_circle(c, r_in, mesh_size=resolution)
circ_outer = model.add_circle(c, r_out, mesh_size=resolution)

# Create plane surface for meshing
plane_surface = model.add_plane_surface(
    circ_outer.curve_loop, holes=[circ_inner.curve_loop])

# Call gmsh kernel before add physical entities
model.synchronize()

volume_marker = 6
model.add_physical(circ_outer.curve_loop.curves, "Outer")
model.add_physical(circ_inner.curve_loop.curves, "Inner")

geometry.generate_mesh(dim=2)
gmsh.write("mesh.msh")
gmsh.clear()
geometry.__exit__()

mesh_from_file = meshio.read("mesh.msh")

def create_mesh(mesh, cell_type, prune_z=False):
    cells = mesh.get_cells_type(cell_type)
    cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
    points = mesh.points[:,:2] if prune_z else mesh.points
    out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
    return out_mesh

line_mesh = create_mesh(mesh_from_file, "line", prune_z=True)
meshio.write("facet_mesh.xdmf", line_mesh)

triangle_mesh = create_mesh(mesh_from_file, "triangle", prune_z=True)
meshio.write("mesh.xdmf", triangle_mesh)

and error txt:

Cell In [111], line 35, in create_mesh(mesh, cell_type, prune_z)
     33 def create_mesh(mesh, cell_type, prune_z=False):
     34     cells = mesh.get_cells_type(cell_type)
---> 35     cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
     36     points = mesh.points[:,:2] if prune_z else mesh.points
     37     out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})

File /usr/local/lib/python3.10/dist-packages/meshio/_mesh.py:249, in Mesh.get_cell_data(self, name, cell_type)
    248 def get_cell_data(self, name: str, cell_type: str):
--> 249     return np.concatenate(
    250         [d for c, d in zip(self.cells, self.cell_data[name]) if c.type == cell_type]
    251     )

File <__array_function__ internals>:180, in concatenate(*args, **kwargs)

ValueError: need at least one array to concatenate
12

You have not marked the physical surface (volume of the mesh, and thus it cannot load it), see for instance: Problem converting gmsh to xdmf for fenics - #2 by dokken

1 Like
13

Thanks for previous advices, now I created mesh for my problem successfully. But when I came up to solver, Incompatible mesh error appeared, what can be the reason why? Code:

import dolfinx
import gmsh
!pip install meshio
import meshio
!pip install pygmsh
import pygmsh
!pip install h5py
import h5py
from mpi4py import MPI
print(f"DOLFINx version: {dolfinx.__version__} based on GIT commit: {dolfinx.git_commit_hash} of https://github.com/FEniCS/dolfinx/")

# Radius of the inner and outer boubdary of the domain
r_in = 0.5
r_out = 1.5
c = [ 0, 0, 0 ]

model_rank = 0
mesh_comm = MPI.COMM_WORLD

# Initialize empty geometry using the build in kernel in GMSH
geometry = pygmsh.geo.Geometry()
# Fetch model we would like to add data to
model = geometry.__enter__()
# Add circle
circ_inner = model.add_circle(c, r_in, mesh_size=resolution)
circ_outer = model.add_circle(c, r_out, mesh_size=resolution)

# Create plane surface for meshing
plane_surface = model.add_plane_surface(
    circ_outer.curve_loop, holes=[circ_inner.curve_loop])

# Call gmsh kernel before add physical entities
model.synchronize()

volume_marker = 6
model.add_physical(circ_outer.curve_loop.curves, "3")
model.add_physical(circ_inner.curve_loop.curves, "Inner")
model.add_physical([plane_surface], "Domain")

geometry.generate_mesh(dim=2)
from dolfinx.io.gmshio import model_to_mesh
mesh, ct, ft = model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=2)
gmsh.finalize()

import pyvista
pyvista.set_jupyter_backend("pythreejs")
from dolfinx.plot import create_vtk_mesh

plotter = pyvista.Plotter()
grid = pyvista.UnstructuredGrid(*create_vtk_mesh(mesh, mesh.topology.dim))
num_local_cells = mesh.topology.index_map(mesh.topology.dim).size_local
grid.cell_data["Marker"] = ct.values[ct.indices<num_local_cells]
grid.set_active_scalars("Marker")
actor = plotter.add_mesh(grid, color='cyan',  show_edges=True)
plotter.view_xy()
if not pyvista.OFF_SCREEN:
    plotter.show()
else:
    pyvista.start_xvfb()
    cell_tag_fig = plotter.screenshot("cell_tags.png")

import numpy as np
!pip install sympy
import sympy
from petsc4py import PETSc
import ufl
from ufl import (FacetNormal, as_vector, conj, cross, curl, inner, lhs, rhs, sqrt)

if not np.issubdtype(PETSc.ScalarType, np.complexfloating):
    print("Demo should only be executed with DOLFINx complex mode")
    exit(0)

class BackgroundElectricField:

    def __init__(self, theta, n, k0):
        self.theta = theta
        self.k0 = k0
        self.n = n

    def eval(self, x):

        kx = self.n * self.k0 * np.cos(self.theta)
        ky = self.n  * self.k0 * np.sin(self.theta)
        phi = kx * x[0] + ky * x[1]

        ax = np.sin(self.theta)
        ay = np.cos(self.theta)

        return (-ax * np.exp(1j * phi), ay * np.exp(1j * phi))


def radial_distance(x):
    """Returns the radial distance from the origin"""
    return ufl.sqrt(x[0]**2 + x[1]**2)

degree = 3
curl_el = ufl.FiniteElement("N1curl", mesh.ufl_cell(), degree)
V = dolfinx.fem.FunctionSpace(mesh, curl_el)


f = BackgroundElectricField(theta=0, n=1, k0=2*np.pi)
Einc = dolfinx.fem.Function(V)
Einc.interpolate(f.eval)

x = ufl.SpatialCoordinate(mesh)
r = radial_distance(x)

# Definition of Trial and Test functions
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

# Definition of 3d fields
u_3d = as_vector((u[0], u[1], 0))
v_3d = as_vector((v[0], v[1], 0))

# Measures for subdomains
dx = ufl.Measure("dx", mesh, subdomain_data=ct)
ds = ufl.Measure("ds", mesh, subdomain_data=ft)

dsbc = ds(3)

# Normal to the boundary
n = FacetNormal(mesh)
n_3d = as_vector((n[0], n[1], 0))


# Weak form
k = 2*np.pi
eps = 1
mu = 1
beta = eps*mu*k**2
F = - inner(ufl.grad(u),ufl.grad(v))*dx - (1j * k + 1 / (2 * r_out))*inner(u,v)*ds + beta*inner(u,v)*dx


# Splitting in left-hand side and right-hand side
a, L = lhs(F), rhs(F)

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

Please actually state your error message (with the backtrace), and make sure that the code is executable (all variables defined.

There are several errors in your code, and it is also assuming that you are using colab or a notebook, something you should mention.
The following code, ran with the dolfinx docker image dolfinx/lab:nightly does not produce any errors:

import dolfinx

import gmsh
!pip install meshio
import meshio
!pip install pygmsh
import pygmsh
!pip install h5py
import h5py
from mpi4py import MPI
print(f"DOLFINx version: {dolfinx.__version__} based on GIT commit: {dolfinx.git_commit_hash} of https://github.com/FEniCS/dolfinx/")

# Radius of the inner and outer boubdary of the domain
resolution = 0.04
r_in = 0.5
r_out = 1.5
c = [ 0, 0, 0 ]

model_rank = 0
mesh_comm = MPI.COMM_WORLD

# Initialize empty geometry using the build in kernel in GMSH
geometry = pygmsh.geo.Geometry()
# Fetch model we would like to add data to
model = geometry.__enter__()
# Add circle
circ_inner = model.add_circle(c, r_in, mesh_size=resolution)
circ_outer = model.add_circle(c, r_out, mesh_size=resolution)

# Create plane surface for meshing
plane_surface = model.add_plane_surface(
    circ_outer.curve_loop, holes=[circ_inner.curve_loop])

# Call gmsh kernel before add physical entities
model.synchronize()

volume_marker = 6
model.add_physical(circ_outer.curve_loop.curves, "3")
model.add_physical(circ_inner.curve_loop.curves, "Inner")
model.add_physical([plane_surface], "Domain")

geometry.generate_mesh(dim=2)
from dolfinx.io.gmshio import model_to_mesh
mesh, ct, ft = model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=2)
gmsh.finalize()


import pyvista
pyvista.set_jupyter_backend("pythreejs")
from dolfinx.plot import create_vtk_mesh

plotter = pyvista.Plotter()
grid = pyvista.UnstructuredGrid(*create_vtk_mesh(mesh, mesh.topology.dim))
num_local_cells = mesh.topology.index_map(mesh.topology.dim).size_local
grid.cell_data["Marker"] = ct.values[ct.indices<num_local_cells]
grid.set_active_scalars("Marker")
actor = plotter.add_mesh(grid, color='cyan',  show_edges=True)
plotter.view_xy()
print(pyvista.OFF_SCREEN)
if not pyvista.OFF_SCREEN:
    plotter.show()
else:
    pyvista.start_xvfb()
    cell_tag_fig = plotter.screenshot("cell_tags.png")

import numpy as np
!pip install sympy
import sympy
from petsc4py import PETSc
import ufl
from ufl import (FacetNormal, as_vector, conj, cross, curl, inner, lhs, rhs, sqrt)

if not np.issubdtype(PETSc.ScalarType, np.complexfloating):
    print("Demo should only be executed with DOLFINx complex mode")
    exit(0)

class BackgroundElectricField:

    def __init__(self, theta, n, k0):
        self.theta = theta
        self.k0 = k0
        self.n = n

    def eval(self, x):

        kx = self.n * self.k0 * np.cos(self.theta)
        ky = self.n  * self.k0 * np.sin(self.theta)
        phi = kx * x[0] + ky * x[1]

        ax = np.sin(self.theta)
        ay = np.cos(self.theta)

        return (-ax * np.exp(1j * phi), ay * np.exp(1j * phi))


def radial_distance(x):
    """Returns the radial distance from the origin"""
    return ufl.sqrt(x[0]**2 + x[1]**2)

degree = 3
curl_el = ufl.FiniteElement("N1curl", mesh.ufl_cell(), degree)
V = dolfinx.fem.FunctionSpace(mesh, curl_el)


f = BackgroundElectricField(theta=0, n=1, k0=2*np.pi)
Einc = dolfinx.fem.Function(V)
Einc.interpolate(f.eval)

x = ufl.SpatialCoordinate(mesh)
r = radial_distance(x)

# Definition of Trial and Test functions
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

# Definition of 3d fields
u_3d = as_vector((u[0], u[1], 0))
v_3d = as_vector((v[0], v[1], 0))

# Measures for subdomains
dx = ufl.Measure("dx", mesh, subdomain_data=ct)
ds = ufl.Measure("ds", mesh, subdomain_data=ft)

dsbc = ds(3)

# Normal to the boundary
n = FacetNormal(mesh)
n_3d = as_vector((n[0], n[1], 0))


# Weak form
k = 2*np.pi
eps = 1
mu = 1
beta = eps*mu*k**2
F = - inner(ufl.grad(u),ufl.grad(v))*dx - (1j * k + 1 / (2 * r_out))*inner(u,v)*ds + beta*inner(u,v)*dx + inner(dolfinx.fem.Constant(mesh, PETSc.ScalarType((0,0))), v)*dx(9999)


# Splitting in left-hand side and right-hand side
a, L = lhs(F), rhs(F)

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

I am using dolfinx/lab:latest, 0.6.0 (GitHub - FEniCS/dolfinx: Next generation FEniCS problem solving environment) with jupyter notebook. Full error message:

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
Cell In [38], line 4
      1 # Splitting in left-hand side and right-hand side
      2 a, L = lhs(F), rhs(F)
----> 4 problem = dolfinx.fem.petsc.LinearProblem(a, L, bcs=[], petsc_options={
      5                                   "ksp_type": "preonly", "pc_type": "lu"})
      6 Eh = problem.solve()

File /usr/local/dolfinx-complex/lib/python3.10/dist-packages/dolfinx/fem/petsc.py:566, in LinearProblem.__init__(self, a, L, bcs, u, petsc_options, form_compiler_options, jit_options)
    538 def __init__(self, a: ufl.Form, L: ufl.Form, bcs: typing.List[DirichletBCMetaClass] = [],
    539              u: typing.Optional[_Function] = None, petsc_options={}, form_compiler_options={}, jit_options={}):
    540     """Initialize solver for a linear variational problem.
    541 
    542     Args:
   (...)
    564                                                "pc_factor_mat_solver_type": "mumps"})
    565     """
--> 566     self._a = _create_form(a, form_compiler_options=form_compiler_options, jit_options=jit_options)
    567     self._A = create_matrix(self._a)
    569     self._L = _create_form(L, form_compiler_options=form_compiler_options, jit_options=jit_options)

File /usr/local/dolfinx-complex/lib/python3.10/dist-packages/dolfinx/fem/forms.py:176, in form(form, dtype, form_compiler_options, jit_options)
    173         return list(map(lambda sub_form: _create_form(sub_form), form))
    174     return form
--> 176 return _create_form(form)

File /usr/local/dolfinx-complex/lib/python3.10/dist-packages/dolfinx/fem/forms.py:171, in form.<locals>._create_form(form)
    168 """Recursively convert ufl.Forms to dolfinx.fem.Form, otherwise
    169 return form argument"""
    170 if isinstance(form, ufl.Form):
--> 171     return _form(form)
    172 elif isinstance(form, collections.abc.Iterable):
    173     return list(map(lambda sub_form: _create_form(sub_form), form))

File /usr/local/dolfinx-complex/lib/python3.10/dist-packages/dolfinx/fem/forms.py:165, in form.<locals>._form(form)
    159 # Subdomain markers (possibly None for some dimensions)
    160 subdomains = {_cpp.fem.IntegralType.cell: subdomains.get("cell"),
    161               _cpp.fem.IntegralType.exterior_facet: subdomains.get("exterior_facet"),
    162               _cpp.fem.IntegralType.interior_facet: subdomains.get("interior_facet"),
    163               _cpp.fem.IntegralType.vertex: subdomains.get("vertex")}
--> 165 return formcls(ufcx_form, V, coeffs, constants, subdomains, mesh, module.ffi, code)

File /usr/local/dolfinx-complex/lib/python3.10/dist-packages/dolfinx/fem/forms.py:51, in FormMetaClass.__init__(self, form, V, coeffs, constants, subdomains, mesh, ffi, code)
     49 self._code = code
     50 self._ufcx_form = form
---> 51 super().__init__(ffi.cast("uintptr_t", ffi.addressof(self._ufcx_form)),
     52                  V, coeffs, constants, subdomains, mesh)

RuntimeError: Incompatible mesh
16

You should not use the latest tag, but nightly, as latest hasnt been updated for a month.
I would suggest calling docker pull dolfinx/lab:nightly prior to starting a new container.

17

Hello, now everything works without errors, but I get zero solution. Can it be because of mistake in stating boundary conditions or in weak formulation of my problem? Because for me everything looks fine:

import dolfinx
import gmsh
!pip install meshio
import meshio
!pip install pygmsh
import pygmsh
!pip install h5py
import h5py
from mpi4py import MPI
import numpy as np
print(f"DOLFINx version: {dolfinx.__version__} based on GIT commit: {dolfinx.git_commit_hash} of https://github.com/FEniCS/dolfinx/")

# Radius of the inner and outer boubdary of the domain
resolution = 0.04
r_in = 0.5
r_out = 1.5
c = [ 0, 0, 0 ]

gmsh.initialize()

model_rank = 0
mesh_comm = MPI.COMM_WORLD


# Define geometry for iron cylinder
circ_inner = gmsh.model.occ.addDisk(0, 0, 0, r_in,r_in)
circ_outer = gmsh.model.occ.addDisk(0, 0, 0, r_out,r_out)

gmsh.model.occ.synchronize()

# Resolve all boundaries of the different wires in the background domain
all_surfaces = [(2, circ_inner)]
whole_domain = gmsh.model.occ.fragment([(2, circ_outer)], all_surfaces)
gmsh.model.occ.synchronize()

# Create physical markers for the different wires.
# We use the following markers:
# - Inner: 0
# - Domain: 1
background_surfaces = []
other_surfaces = []
for domain in whole_domain[0]:
    com = gmsh.model.occ.getCenterOfMass(domain[0], domain[1])
    mass = gmsh.model.occ.getMass(domain[0], domain[1])
    # Identify iron circle by its mass
    if np.isclose(mass, np.pi*(r_in**2)):
        gmsh.model.addPhysicalGroup(domain[0], [domain[1]], tag=103)
        other_surfaces.append(domain)
    # Identify the domain circle by its center of mass
    elif np.allclose(com, [0, 0, 0]):
        background_surfaces.append(domain[1])

# Add marker for the vacuum
gmsh.model.addPhysicalGroup(2, background_surfaces, tag=104)

gmsh.model.mesh.field.add("Distance", 1)
edges = gmsh.model.getBoundary(other_surfaces, oriented=False)
gmsh.model.mesh.field.setNumbers(1, "EdgesList", [e[1] for e in edges])
gmsh.model.mesh.field.add("Threshold", 2)
gmsh.model.mesh.field.setNumber(2, "IField", r_out)
gmsh.model.mesh.field.setNumber(2, "LcMin", 0.04)
gmsh.model.mesh.field.setNumber(2, "LcMax", 0.04)
gmsh.model.mesh.field.setNumber(2, "DistMin", r_in)
gmsh.model.mesh.field.setNumber(2, "DistMax", r_out)
gmsh.model.mesh.field.setAsBackgroundMesh(2)


# Generate mesh
gmsh.model.mesh.generate(2)


from dolfinx.io.gmshio import model_to_mesh
mesh, ct, _ = model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=2)
gmsh.finalize()

import pyvista
pyvista.set_jupyter_backend("pythreejs")
from dolfinx.plot import create_vtk_mesh

plotter = pyvista.Plotter()
grid = pyvista.UnstructuredGrid(*create_vtk_mesh(mesh, mesh.topology.dim))
num_local_cells = mesh.topology.index_map(mesh.topology.dim).size_local
grid.cell_data["Marker"] = ct.values[ct.indices<num_local_cells]
grid.set_active_scalars("Marker")
actor = plotter.add_mesh(grid, show_edges=True)
plotter.view_xy()
if not pyvista.OFF_SCREEN:
    plotter.show()
else:
    pyvista.start_xvfb()
    cell_tag_fig = plotter.screenshot("cell_tags.png")

from dolfinx.fem import (dirichletbc, Expression, Function, FunctionSpace, 
                         VectorFunctionSpace, locate_dofs_topological)
from dolfinx.fem.petsc import LinearProblem
from dolfinx.mesh import locate_entities_boundary
from ufl import TestFunction, TrialFunction, as_vector, dot, dx, grad, inner, ds, lhs, rhs, SpatialCoordinate, sqrt, FacetNormal
from petsc4py.PETSc import ScalarType

if not np.issubdtype(ScalarType, np.complexfloating):
    print("Demo should only be executed with DOLFINx complex mode")
    exit(0)

class BackgroundElectricField:

    def __init__(self, n, k0):
        self.k0 = k0
        self.n = n

    def eval(self, x):
        return (np.exp(-1j * self.k0*self.n*x[0]))


def radial_distance(x):
    """Returns the radial distance from the origin"""
    return sqrt(x[0]**2 + x[1]**2)

n = FacetNormal(mesh)

V = FunctionSpace(mesh, ("CG", 2))
tdim = mesh.topology.dim
facets = locate_entities_boundary(mesh, tdim-1, lambda x: np.full(x.shape[1], True))
dofs = locate_dofs_topological(V, tdim-1, facets)
bc = dirichletbc(ScalarType(0), dofs, V)
 
x = SpatialCoordinate(mesh)
f = BackgroundElectricField(n=1, k0=2*np.pi)
Einc = dolfinx.fem.Function(V)
Einc.interpolate(f.eval)

k = 2*np.pi
eps = 1
mu = 1
beta = eps*mu*k**2

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

F = (beta*inner(u,v)-inner(grad(u),grad(v)))*dx + (1j * k + 1 / (2 * r_out))*(inner(Einc,v)-inner(u,v))*ds  + inner(inner(grad(Einc),n),v)*ds

# Splitting in left-hand side and right-hand side
a, L = lhs(F), rhs(F)

Esc = Function(V)
problem = LinearProblem(a, L, u=Esc, bcs=[bc])
problem.solve()

plotter = pyvista.Plotter()

Esc_grid = pyvista.UnstructuredGrid(*create_vtk_mesh(V))
Esc_grid.point_data["Esc"] = Esc.x.array
Esc_grid.set_active_scalars("Esc")
warp = Esc_grid.warp_by_scalar("Esc", factor=1e7)
actor = plotter.add_mesh(warp, show_edges=True)
if not pyvista.OFF_SCREEN:
    plotter.show()
else:
    pyvista.start_xvfb()
    Esc_fig = plotter.screenshot("Esc.png")

I expected getting something like this one:

image
image779×688 55.2 KB

and after that to plot angular dependence.