Mapping a 2D field distribution to a 3D domain

Question

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:

Does anyone know how to fix this issue? Is it maybe a bug of DirichletBC?

Kind regards and thank you for your valuable support
Erna

Comments

Dear FEniCS team,

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
Erna

Answer 1

Whoever is interested in this question, we solved this now (very very dirty) with the following code and it does work but is not a permanent solution:

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


global_element_order = 2

# Specify the grid
mesh_density = 50
points = 1000

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

p1 = Point(0,0,0)
p2 = Point(a,b,L)
box1 = Box(Point(0,0,0), Point(a,b,L))

mesh = RectangleMesh(p1,p2,mesh_density,mesh_density)


domain = box1
  # Resolution of mesh
mesh = generate_mesh(domain, mesh_density)

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

#define left boundary
class BoundaryLeft(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and abs(x[2]-L)<DOLFIN_EPS

boundaryLeft = BoundaryLeft()

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

u= interpolate(u, V_ned)

plot(u)

#plot

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

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