MixedElement vs MixedFunctionSpace

variational formulation
1

Hello everyone :smiley:

I’m trying to solve a system of two PDEs with two unknowns, called \mathbf{E} and \mathbf{J}. I have a 2D domain made up of two subdomains, that in the figure below are labeled with 1 and 2:

image

\mathbf{E} is defined in both subdomains, while \mathbf{J} only exists in the subdomain 1 .

My question is: should I formulate my problem by using MixedElement or MixedFunctionSpace?

I’m asking because I tried to formulate the problem with MixedElement but I’m getting the Solution failed to converge in 0 iterations error, and then I read that MixedFunctionSpace could be the way to go when you have unknowns only defined in a sub-mesh, but I’m not sure about that.

2

Anyone that knows how to proceed? :smiley:

3

Hi all,

I’m still struggling with this problem. Using MixedFunctionSpace does not give me a convergent solution, probably because I’m not using it correctly.

So I was wondering: is there any way to put a “volume” constraint on my unknown \mathbf{J}? In other words, is there an equivalent of DirichletBC but acting on cells instead of facets? It might solve my problem.

4

I wrote some examples that show how extending the domain causes the error (even if they do not describe my complete problem). Please note that the script can be run directly, the only requirement is to have the gmsh python API installed, nothing else (therefore, even if they might not be the shortest working examples possible, it is really easy to run them).

In the script below, I’m using only the domain 1, and everything works fine:

# DOMAIN 1 ---> WORKING FINE

from dolfin import *
import numpy as np
import gmsh
import sys
import meshio

um = 10**-6

def msh_to_mf():

    dim = 2
    cell_str = "triangle"
    facet_str = "line"

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

    cells = np.vstack(np.array([cells.data for cells in msh.cells
                                if cells.type == cell_str]))

    total_mesh = meshio.Mesh(points=msh.points,
                             cells=[(cell_str, cells)])
    total_mesh.prune_z_0()

    cells_data = msh.cell_data_dict["gmsh:physical"][cell_str]

    cells_mesh = meshio.Mesh(points=msh.points,
                                cells=[(cell_str, cells)],
                                cell_data={"name_to_read": [cells_data]})

    cells_mesh.prune_z_0()
    facet_cells = np.vstack(np.array([cells.data for cells in msh.cells
                                      if cells.type == facet_str]))

    facet_data = msh.cell_data_dict["gmsh:physical"][facet_str]

    facet_mesh = meshio.Mesh(points=msh.points,
                             cells=[(facet_str, facet_cells)],
                             cell_data={"name_to_read": [facet_data]})

    facet_mesh.prune_z_0()
    meshio.write("total.xdmf", total_mesh)
    meshio.write("facets.xdmf", facet_mesh)
    meshio.write("cells.xdmf", cells_mesh)

    parameters["ghost_mode"] = "shared_vertex"

    mesh = Mesh()

    with XDMFFile("cells.xdmf") as infile:
        infile.read(mesh)

    mvc_facets = MeshValueCollection("size_t", mesh, dim-1)

    with XDMFFile("facets.xdmf") as infile:
        infile.read(mvc_facets, "name_to_read")

    mf_facets = cpp.mesh.MeshFunctionSizet(mesh, mvc_facets)

    mvc_cells = MeshValueCollection("size_t", mesh, dim)

    with XDMFFile("cells.xdmf") as infile:
        infile.read(mvc_cells, "name_to_read")

    mf_cells = cpp.mesh.MeshFunctionSizet(mesh, mvc_cells)

    return mesh, mf_facets, mf_cells

gmsh.initialize(sys.argv)
gmsh.model.add("geometry")
gmsh.model.occ.addCircle(0, 0, 0, 0.005*um, tag=1)
gmsh.model.occ.addCurveLoop([1], tag=1)
gmsh.model.occ.addPlaneSurface([1], tag=1)
gmsh.model.occ.synchronize()
gmsh.model.addPhysicalGroup(2, [1], tag=1)
gmsh.model.addPhysicalGroup(1, [1], tag=2)
gmsh.model.mesh.setSize([(0, 1)], size=0.0001*um)
gmsh.model.mesh.generate(2)
gmsh.write("mesh.msh")

mesh, mf_facets, mf_cells = msh_to_mf()
mesh.init_cell_orientations(Expression(("x[0]", "x[1]", "x[2]"), degree = 2))
wl0 = 0.5*um

dx = Measure("dx", domain=mesh, subdomain_data=mf_cells)

div_el = FiniteElement('N1div', mesh.ufl_cell(), 2)

V = FunctionSpace(mesh, div_el*div_el)

(rP, iP) = TrialFunctions(V)
(rK, iK) = TestFunctions(V)

rPn = Expression(("0.0", "0.0"), domain=mesh, degree=2)
iPn = Expression(("0.0", "0.0"), domain=mesh, degree=2)

zero_normal_real = DirichletBC(V.sub(0), rPn, mf_facets, 2, method="geometric")
zero_normal_imag = DirichletBC(V.sub(1), iPn, mf_facets, 2, method="geometric")

bcs = [zero_normal_real, zero_normal_imag]

rEb = Expression(("0.0", "cos(-2*pi/wl*x[0])"),pi=np.pi, wl=wl0, domain=mesh, degree=2)
iEb = Expression(("0.0", "sin(-2*pi/wl*x[0])"),pi=np.pi, wl=wl0, domain=mesh, degree=2)

P_term = -(dot(rP, rK))*dx(1) \
           -(dot(iP, iK))*dx(1) \
           -(dot(rP, iK))*dx(1) \
           +(dot(iP, rK))*dx(1)

Eb_term = -(dot(rEb, rK))*dx(1) \
            -(dot(iEb, iK))*dx(1)

equation = P_term + Eb_term

F = equation

a, L = lhs(F), rhs(F)

U = Function(V)
solve(a == L, U, bcs, solver_parameters={'linear_solver':'mumps'})
(rPm, iPm) = U.split()

As soon as I add the domain 2, I’m getting the error, as running the following script show. I think this is because the weak form does not have any information on my trial functions in the domain 2.

# DOMAIN 1 + DOMAIN 2 ---> Solution failed to converge in 0 iterations

from dolfin import *
import numpy as np
import gmsh
import sys 
import meshio

um = 10**-6

def msh_to_mf():

    dim = 2 
    cell_str = "triangle"
    facet_str = "line"

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

    cells = np.vstack(np.array([cells.data for cells in msh.cells
                                if cells.type == cell_str]))

    total_mesh = meshio.Mesh(points=msh.points,
                             cells=[(cell_str, cells)])
    total_mesh.prune_z_0()

    cells_data = msh.cell_data_dict["gmsh:physical"][cell_str]

    cells_mesh = meshio.Mesh(points=msh.points,
                                cells=[(cell_str, cells)],
                                cell_data={"name_to_read": [cells_data]})

    cells_mesh.prune_z_0()
    facet_cells = np.vstack(np.array([cells.data for cells in msh.cells
                                      if cells.type == facet_str]))

    facet_data = msh.cell_data_dict["gmsh:physical"][facet_str]

    facet_mesh = meshio.Mesh(points=msh.points,
                             cells=[(facet_str, facet_cells)],
                             cell_data={"name_to_read": [facet_data]})

    facet_mesh.prune_z_0()
    meshio.write("total.xdmf", total_mesh)
    meshio.write("facets.xdmf", facet_mesh)
    meshio.write("cells.xdmf", cells_mesh)

    parameters["ghost_mode"] = "shared_vertex"

    mesh = Mesh()

    with XDMFFile("cells.xdmf") as infile:
        infile.read(mesh)

    mvc_facets = MeshValueCollection("size_t", mesh, dim-1)

    with XDMFFile("facets.xdmf") as infile:
        infile.read(mvc_facets, "name_to_read")

    mf_facets = cpp.mesh.MeshFunctionSizet(mesh, mvc_facets)

    mvc_cells = MeshValueCollection("size_t", mesh, dim)

    with XDMFFile("cells.xdmf") as infile:
        infile.read(mvc_cells, "name_to_read")

    mf_cells = cpp.mesh.MeshFunctionSizet(mesh, mvc_cells)

    return mesh, mf_facets, mf_cells

gmsh.initialize(sys.argv)
gmsh.model.add("geometry")
gmsh.model.occ.addCircle(0, 0, 0, 0.005*um, tag=1)
gmsh.model.occ.addCircle(0, 0, 0, 0.05*um, tag=2)   #added
gmsh.model.occ.addCurveLoop([1], tag=1)
gmsh.model.occ.addPlaneSurface([1], tag=1)
gmsh.model.occ.addCurveLoop([2], tag=2)             #added
gmsh.model.occ.addCurveLoop([1], tag=3)             #added
gmsh.model.occ.addPlaneSurface([2, 3], tag=2)       #added
gmsh.model.occ.synchronize()
gmsh.model.addPhysicalGroup(2, [1], tag=1)
gmsh.model.addPhysicalGroup(1, [1], tag=2)
gmsh.model.addPhysicalGroup(2, [2], tag=3)          #added
gmsh.model.mesh.setSize([(0, 1)], size=0.0001*um)
gmsh.model.mesh.setSize([(0, 2)], size=0.02*um)     #added
gmsh.model.mesh.generate(2)
gmsh.write("mesh.msh")

mesh, mf_facets, mf_cells = msh_to_mf()
mesh.init_cell_orientations(Expression(("x[0]", "x[1]", "x[2]"), degree = 2)) 
wl0 = 0.5*um

dx = Measure("dx", domain=mesh, subdomain_data=mf_cells)

div_el = FiniteElement('N1div', mesh.ufl_cell(), 2)

V = FunctionSpace(mesh, div_el*div_el)

(rP, iP) = TrialFunctions(V)
(rK, iK) = TestFunctions(V)

rPn = Expression(("0.0", "0.0"), domain=mesh, degree=2)
iPn = Expression(("0.0", "0.0"), domain=mesh, degree=2)

zero_normal_real = DirichletBC(V.sub(0), rPn, mf_facets, 2, method="geometric")
zero_normal_imag = DirichletBC(V.sub(1), iPn, mf_facets, 2, method="geometric")

bcs = [zero_normal_real, zero_normal_imag]

rEb = Expression(("0.0", "cos(-2*pi/wl*x[0])"),pi=np.pi, wl=wl0, domain=mesh, degree=2)
iEb = Expression(("0.0", "sin(-2*pi/wl*x[0])"),pi=np.pi, wl=wl0, domain=mesh, degree=2)

P_term = -(dot(rP, rK))*dx(1) \
           -(dot(iP, iK))*dx(1) \
           -(dot(rP, iK))*dx(1) \
           +(dot(iP, rK))*dx(1)

Eb_term = -(dot(rEb, rK))*dx(1) \
            -(dot(iEb, iK))*dx(1)

equation = P_term + Eb_term

F = equation

a, L = lhs(F), rhs(F)

U = Function(V)
solve(a == L, U, bcs, solver_parameters={'linear_solver':'mumps'})
(rPm, iPm) = U.split()

The error is this one:

Traceback (most recent call last):
  File "n1div_single_equation_not_working.py", line 125, in <module>
    solve(a == L, U, bcs, solver_parameters={'linear_solver':'mumps'})
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/solving.py", line 233, in solve
    _solve_varproblem(*args, **kwargs)
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/solving.py", line 273, in _solve_varproblem
    solver.solve()
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 solve linear system using PETSc Krylov solver.
*** Reason:  Solution failed to converge in 0 iterations (PETSc reason DIVERGED_PC_FAILED, residual norm ||r|| = 0.000000e+00).
*** Where:   This error was encountered inside PETScKrylovSolver.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.2.0.dev0
*** Git changeset:  unknown
*** -------------------------------------------------------------------------

Therefore I’m asking again: how can I specify to dolfin that my unknown is = 0 in the domain 2? Or is there any other way to make the simulation run? Thank you.

5

I suspect that you’re not posing the problem correctly in the mathematical sense. Unknowns don’t typically “exist” in one domain only. There needs to be appropriate continuity conditions at the interface. Ask a math friend to look over your problem statement.

1 Like
6

Thanks @peb for your answer! Now it’s been a while since I started working on this problem, and the mathematical framework and conditions are quite clear to me (I hope :smiley: !). As you said, at the boundary between the two domains, I need a proper continuity condition. In particular, in my case, the normal component of \mathbf{J} must be zero, i.e. \mathbf{J}\cdot\mathbf{n} = 0, therefore I’ve chosen the N1div elements for solving this problem.

For further details, you can have a look at this article, where there is a clear explanation of the mathematical framework for the problem. Having said that, I will look again at the math of the problem, to be sure that everything is well-posed.

In the meanwhile, I was able to run a simulation by using MixedFunctionSpace, but the solution that I get is not what I should obtain, therefore I think there is something missing or wrong in the way I used MixedFunctionSpace. I suspect it is in the DirichletBC domain, but I do not know what I should change. Here the plug-and-play script:

from dolfin import *
import numpy as np
import gmsh
import sys
import meshio
import scipy.constants as const
import os

#---------------------------------------------
um = 10**-6

def msh_to_mf():

    dim = 2
    cell_str = "triangle"
    facet_str = "line"

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

    cells = np.vstack(np.array([cells.data for cells in msh.cells
                                if cells.type == cell_str]))

    total_mesh = meshio.Mesh(points=msh.points,
                             cells=[(cell_str, cells)])
    total_mesh.prune_z_0()

    cells_data = msh.cell_data_dict["gmsh:physical"][cell_str]

    cells_mesh = meshio.Mesh(points=msh.points,
                                cells=[(cell_str, cells)],
                                cell_data={"name_to_read": [cells_data]})

    cells_mesh.prune_z_0()
    facet_cells = np.vstack(np.array([cells.data for cells in msh.cells
                                      if cells.type == facet_str]))

    facet_data = msh.cell_data_dict["gmsh:physical"][facet_str]

    facet_mesh = meshio.Mesh(points=msh.points,
                             cells=[(facet_str, facet_cells)],
                             cell_data={"name_to_read": [facet_data]})

    facet_mesh.prune_z_0()
    meshio.write("total.xdmf", total_mesh)
    meshio.write("facets.xdmf", facet_mesh)
    meshio.write("cells.xdmf", cells_mesh)

    parameters["ghost_mode"] = "shared_vertex"

    mesh = Mesh()

    with XDMFFile("cells.xdmf") as infile:
        infile.read(mesh)

    mvc_facets = MeshValueCollection("size_t", mesh, dim-1)

    with XDMFFile("facets.xdmf") as infile:
        infile.read(mvc_facets, "name_to_read")

    mf_facets = cpp.mesh.MeshFunctionSizet(mesh, mvc_facets)

    mvc_cells = MeshValueCollection("size_t", mesh, dim)

    with XDMFFile("cells.xdmf") as infile:
        infile.read(mvc_cells, "name_to_read")

    mf_cells = cpp.mesh.MeshFunctionSizet(mesh, mvc_cells)

    return mesh, mf_facets, mf_cells

gmsh.initialize(sys.argv)
gmsh.model.add("geometry")
gmsh.model.occ.addCircle(0, 0, 0, 0.005*um, tag=1)
gmsh.model.occ.addCircle(0, 0, 0, 0.05*um, tag=2)
gmsh.model.occ.addCurveLoop([1], tag=1)
gmsh.model.occ.addPlaneSurface([1], tag=1)
gmsh.model.occ.addCurveLoop([2], tag=2)
gmsh.model.occ.addCurveLoop([1], tag=3)
gmsh.model.occ.addPlaneSurface([2, 3], tag=2)
gmsh.model.occ.synchronize()
gmsh.model.addPhysicalGroup(2, [1], tag=1)
gmsh.model.addPhysicalGroup(1, [1], tag=2)
gmsh.model.addPhysicalGroup(2, [2], tag=3)
gmsh.model.mesh.setSize([(0, 1)], size=0.0004*um)
gmsh.model.mesh.setSize([(0, 2)], size=0.02*um)
gmsh.model.mesh.generate(2)
gmsh.write("mesh.msh")

#---------------------------------------------

mesh, mf_facets, mf_cells = msh_to_mf()
mesh.init_cell_orientations(Expression(("x[0]", "x[1]", "x[2]"), degree = 2))
wl0 = 0.5*um
k0 = 2*np.pi/wl0
omega0 = const.c*k0
mu0 = const.physical_constants["vacuum mag. permeability"][0]
eps0 = const.physical_constants["vacuum electric permittivity"][0]
omegaP = 9*10**15
gamma = 2.5*10**14
beta = 6*10**5

mu0_c = Constant(mu0)
eps0_c = Constant(eps0)
omegaP_c = Constant(omegaP)
gamma_c = Constant(gamma)
beta_c = Constant(beta)
k0_c = Constant(k0)
omega0_c = Constant(omega0)

submesh = MeshView.create(mf_cells, 1)

U1 = FunctionSpace(submesh, "N1div", 2)
U2 = U1
V = MixedFunctionSpace(U1, U2)
dx1 = Measure("dx", domain=V.sub_space(0).mesh())
dx2 = Measure("dx", domain=V.sub_space(1).mesh())

class OnBoundary(SubDomain):

    def inside(self, x, on_boundary):
        return bool(sqrt(x[0]**2+x[1]**2)<=(DOLFIN_EPS+0.005*um) or\
           sqrt(x[0]**2+x[1]**2)>=(0.005*um-DOLFIN_EPS))


onboundary = MeshFunction("size_t", submesh, 1)
onboundary.set_all(0)
myboundary = OnBoundary()
myboundary.mark(onboundary, 1)

(rP, iP) = TrialFunctions(V)
(rK, iK) = TestFunctions(V)

rPn = Expression(("0.0", "0.0"), domain=submesh, degree=2)
iPn = Expression(("0.0", "0.0"), domain=submesh, degree=2)

zero_normal_real = DirichletBC(V.sub_space(0), rPn, onboundary, 1, method='geometric')
zero_normal_imag = DirichletBC(V.sub_space(1), iPn, onboundary, 1, method='geometric')

bcs = [zero_normal_real, zero_normal_imag]

rEb = Expression(("0.0", "cos(-2*pi/wl*x[0])"),pi=np.pi, wl=wl0, domain=submesh, degree=2)
iEb = Expression(("0.0", "sin(-2*pi/wl*x[0])"),pi=np.pi, wl=wl0, domain=submesh, degree=2)

P_term_2 = -(omega0_c**2*dot(rP, rK))*dx1 \
           -(omega0_c**2*dot(iP, iK))*dx2 \
           +(gamma_c*omega0_c*dot(rP, iK))*dx1 \
           -(gamma_c*omega0_c*dot(iP, rK))*dx1

Eb_term_2 = -(eps0_c*omegaP_c**2*dot(rEb, rK))*dx1 \
            -(eps0_c*omegaP_c**2*dot(iEb, iK))*dx1

div_P_term_2 = +(beta_c**2*div(rP)*div(rK))*dx1 \
               +(beta_c**2*div(iP)*div(iK))*dx1

a = P_term_2 + div_P_term_2
L = Eb_term_2

U = Function(V)

solve(a == L, U, bcs)

(rPm, iPm) = U.split()
rPm = project(rPm, FunctionSpace(submesh, "N1div", 2))
iPm = project(iPm, FunctionSpace(submesh, "N1div", 2))
normP = project(sqrt(dot(rPm, rPm)+dot(iPm,iPm)), FunctionSpace(submesh, "CG", 2))

fout = XDMFFile("normP.xdmf")
fout.write_checkpoint(normP, "normP", 0, XDMFFile.Encoding.HDF5, False)
fout.close()

#os.system("paraview "+ "normP.xdmf")

Any suggestion?

7

Hi CastriMik

I wouldn’t mess around with meshio, stick to a simple mesh as proof of concept and mark your boundarys with simple fenics commands.

I haven’t looked in detail at your variational form but I can’t see any coupling there, for this to work there must be one or more instances where (for example) you multiply a test function from E with a trial function from J. Maybe I have missed it - how have you coupled E and J?

For the variational form you could do that by defining mixed elements for E and J over the whole circle, then mark the two subdomains and define your equations for J on domain 1 with dx(1) . Then mark the boundary with 3 and use dS(3) to define boundary conditions for J on the interface. You will probably need to use one of Nitsches, penalty, Lagrange multipliers methods to enforce this and some of (+), (-) jump or avg. I can’t help with that but I suspect Nitsches method would be best. You will also need to use A.ident_zero() to remedy the fact that your J elements are not defined outside domain 1.

1 Like
8

Hi @Thomas2,

Thank you for your answer! :smiley:

Honestly I’ve always struggled in finding the commands for generating meshes with fenics for geometries different from rectangles or cubes, therefore I’ve always used gmsh and meshio, also because gmsh gives me the possibility to quickly visualize my geometry. However, I totally understand that a MWE has much fewer lines by using inner fenics mesh commands (even if the core of the MWE to focus on actually remains unchanged). Therefore, if you know how I can define a circle of a certain radius in FEniCS, tell me, so that I can further shorten my example.

You’re right, in my example there is no coupling of variables, this is because the aim of my MWE wasn’t to show the complete problem, but I just wanted to show the error, that is given by the fact that the unknown \mathbf{J} is only defined in one subdomain. For this purpose, I wrote a simple MWE that is only dependent on \mathbf{J}. Showing the whole problem would be too much complicated to review by other users (even in light of the fact that it is apparently too much complicated to review now).

Thanks for the suggestions, I already heard about the Nitsche method. But to me the problem does not seem too much complicated. After all, what I need here are mainly two things:

  1. Imposing a boundary condition of the type \mathbf{J}\cdot\mathbf{n} = 0, and this can be easily done by using H(div) Nèdèlec elements (N1div ) together with common DirichletBC condition.
  2. Restricting the definition domain of \mathbf{J} to just one subdomain. This should be doable by using MixedFunctionSpace, even if I didn’t get yet why I’m not able to get a physical solution from it.

That being said, I admit that my arguments could be quite naive since I’m still lacking of a deep mathematical and computational knowledge about FEM and fenics and therefore things could be much more complicated that what I said, but, if that is the case, please point out the flaws in my method and show me the right way to go.

9

HI CasriMik

You’re right, in my example there is no coupling of variables, this is because the aim of my MWE wasn’t to show the complete problem,

OK I see, but if its deliberately not coupled then make sure you have an adequate set of BC for both E and J on their own.
I would first try to solve them separately before solving with mixed function spaces together (I guess you have already done that?)

Because you are only defining J on a subset of the full domain then you will probably need to use ident_zeros() - replace solve with something along the lines of

A = assemble(…
b = assemble( …
… …
A.ident_zeros()
…
sol = LUSolver()
sol.solve(A,…,b…

1 Like
10

The FEniCS tutorial includes an example of circular domains using built-in meshes. I found it helpful when I was first starting out with FEniCS (note: it was only a couple of months ago; I’m still new to FEniCS :grinning:). Note this requires the use of mshr.

With regards to your initial question, have you tried applying a DirichletBC to a SubDomain? Adapting Jørgen Dokken’s example on combining Dirichlet and Neumann conditions to your problem, you can create a SubDomain that represents subdomain 2, extract the sub-FunctionSpace V_J corresponding to the trial/test functions \textbf{J}, and create a DirichletBC which is applied only to that subdomain and subspace. Then use assemble_system to create the A matrix and b vector with your boundary conditions incorporated. The following simple example seems to work for me in FEniCS 2019.1.0 (I haven’t tried DOLFINx yet so I can’t say if this would work in that case):

from fenics import *
from mshr import Circle, generate_mesh
from numpy import sqrt

# Create mesh
R_domain    = 2
R_scatterer = 1
domain      = Circle(Point(0, 0), R_domain)
scatterer   = Circle(Point(0, 0), R_scatterer)
domain.set_subdomain(1, scatterer)
mesh        = generate_mesh(domain, 8)

# Define function space
elem       = FiniteElement('P', mesh.ufl_cell(), 1)
mixed_elem = MixedElement([elem, elem])
V          = FunctionSpace(mesh, mixed_elem)
V_J        = V.sub(0)
V_E        = V.sub(1)

# Define boundary condition
u_DJ = Constant(0)
u_DE = Expression('5 -3 * x[0]', degree=2)

def boundary(x, on_boundary):
    return on_boundary

bc0 = DirichletBC(V_J, u_DJ, boundary)
bc1 = DirichletBC(V_E, u_DE, boundary)

# Create a subdomain (R >= 1) in which the the first component
# of the solution is constrained to 0:
class ZeroRegion(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1e-1
        return sqrt(x[0]**2 + x[1]**2) >= R_scatterer-tol

zeroRegion = ZeroRegion()
uJ_zero = Constant(0)
# Create a DirichletBC that applies only to the defined SubDomain and
# to the V_J subspace only, i.e only to the "J" component of the solution
bc_zero = DirichletBC(V_J, uJ_zero, zeroRegion)

# Define variational problem
u_trial = TrialFunction(V)
u_test  = TestFunction(V)
f       = Constant((6,0))
a       = inner(grad(u_trial), grad(u_test))*dx
L       = dot(f,u_test)*dx
sol     = Function(V)
sol_vec = sol.vector()

# Build system matrix and load vector, removing all constrained dofs
A, b = assemble_system(a, L, bcs=[bc0, bc1, bc_zero]) 
solve(A, sol_vec, b)

# Plot the result
from matplotlib import pyplot as plt

u0, u1 = sol.split()
plt.subplot(1,2,1)
h0 = plot(u0)
plot(mesh)
plt.colorbar(h0)
plt.title("$u_0$")

plt.subplot(1,2,2)
h1 = plot(u1)
plot(mesh)
plt.colorbar(h1)
plt.title("$u_1$")
print(dir())
plt.show()

Hope this helps.
-Connor

2 Likes
11

Hi @Thomas2,

I don’t know how to thank you, I was FINALLY able to run the simulation thanks to your suggestions! For clarity, I needed to change your pseudo code a bit, this is the final version:

F = # my weak form
a, L = lhs(F), rhs(F)
A = assemble(a, keep_diagonal=True)
b = assemble(L)
for bc in bcs: bc.apply(A, b)
A.ident_zeros()
sol = PETScLUSolver()
x = Function(V)
sol.solve(A, x.vector(), b)

This solution does not need the use of MixedFunctionSpace but just of MixedElement, I also run the simulation for my system of PDEs and everything run smoothly :smiley:!

12

Hi @conpierce8 ,

Thank you so much for your answer! Finally I was able to find how to generate a circular domain with dolfin, it was hard to find it in the documentation of the libraries, but I completely forgot about this example in the FEniCS tutorial book (my fault! :grinning_face_with_smiling_eyes:).

Your solution looks like what I wanted, and it is wonderful! I still prefer working with gmsh and its labels, therefore the solution of @Thomas2 is the one I’m going to use. However, your approach has given me some useful insights about the function DirichletBC, I didn’t know it could be used to impose conditions over domains other than facets (and it is quite strange since I’m sure I’ve tried this solution in the past, and I’m sure I received some kind of errors that made me give up and try other solutions).

Thank you again, I really appreciate it :grin: