Mapping a 2D field distribution to a 3D domain

[Thomas Flisgen]

Dear FEniCS team,

I would like to map a 2D field distribution which is living on parts of the boundary mesh of a 3D domain to the respective 3D domain. For this purpose, I generated the following code:

from dolfin import *

global_element_order = 2

# Specify the grid
mesh_density = 10

a = 20e-3
b = 10e-3
L = 70e-3

p1 = Point(-a/2,-b/2,-L/2)
p2 = Point(a/2,b/2,L/2)

mesh = BoxMesh(p1,p2,mesh_density,mesh_density,mesh_density)
#mesh = Mesh("waveguide_dense.xml")

# Specify the ansatz function space
V = FunctionSpace(mesh, "Nedelec 1st kind H(curl)", global_element_order)
u = Function(V)

#define left boundary
class BoundaryLeft(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and \
            (near(x[2], +L/2))

boundaryLeft = BoundaryLeft()

bc = DirichletBC(V, Expression(("0.0","cos(pi*x[0]/width)","0.0"),width = a),boundaryLeft)
bc.apply(u.vector())
plot(u)
interactive()

Everything works nicely for the case of a structured grid, which is delivered by BoxMesh. Please have a look at the following plot:

When I use an unstructured mesh (please find the mesh here in order to reproduce my results), the resulting field distribution looks weird. Please have a look at the following plot:

I suppose that this is a bug of DirichletBC?

Kind regards and thank you for your valuable support Thomas

[Prof Garth Wells]

Use [e-mail] for support requests.

You're likely seeing a post-processing artefact; visualisation is vertex-based, which can lead spurious visual effects for H(div) and H(curl) elements.

[Thomas Flisgen]

Dear Garth,

thank you for the quick answer. I did some further investigations. Instead of observing the 3D field distribution vertex-based, I evaluated the x, y, and z-component of the field in the middle of the domain. My code looks as follows:

from dolfin import *
import pylab as pl
import numpy as np

global_element_order = 3

# Specify the grid
mesh_density = 10
points = 10001

a = 20e-3
b = 10e-3
L = 70e-3

p1 = Point(-a/2,-b/2,-L/2)
p2 = Point(a/2,b/2,L/2)

#mesh = BoxMesh(p1,p2,mesh_density,mesh_density,mesh_density)
mesh = Mesh("Test_Meshes/waveguide_dense.xml")

# Specify the ansatz function space
V = FunctionSpace(mesh, "Nedelec 1st kind H(curl)", global_element_order)
u = Function(V)

#define left boundary
class BoundaryLeft(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and \
            (near(x[2], +L/2))

boundaryLeft = BoundaryLeft()

bc = DirichletBC(V, Expression(("0.0","1.0","0.0")),boundaryLeft)
bc.apply(u.vector())
plot(u)

beamline = np.zeros((points,3))
location  = np.linspace(0.99*L/2,L/2,points)
for i in range(points):
            x = np.array([0, 0, location[i]]) 
            beamline[i,:] = u(x)

pl.plot(location,(beamline[:,0]),label='x component')            
pl.plot(location,(beamline[:,1]),label='y component')            
pl.plot(location,(beamline[:,2]),label='z component')
pl.legend(loc='lower left')
pl.show()

Despite the fact that I request the field on the left boundary to be [0,1,0], I obtain a z-component which is actually larger than the y-component (which is equal to 1.0 as required):

I'm wondering whether this behaviour is correct....

Kind regards Thomas

[Prof Garth Wells]

Please don't use the issue tracker for discussion. Use the mailing list or Q&A forum.

[Thomas Flisgen]

Dear Garth,

thank you. I wrote it to the mailing list.

Kind regards Thomas

[Anders Logg (Chalmers)]

Removing milestone: 1.7 (automated comment)