Best way to express constant vectors and the position vector in a weak form

This looks fine. What is the problem?

--
Anders

On Wed, Jun 13, 2012 at 05:15:52PM -0000, Luis Linares wrote:
> Question #200331 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/200331
>
> Description changed to:
> Hello,
>
> I would like to know what is the best way to express constant 3D vectors
> for a Robin boundary condition in my problem. They are originally numpy
> arrays. Here's my approximation at the moment:
>
> import scipy as sp
> from dolfin import *
>
> k0 = 15.0
> k_dir = sp.r_[0.0, 0.0, -1.0]
> E0 = sp.r_[1.0, 0.0, 0.0]
>
> mesh = Mesh("some_mesh_file.xml")
>
> element_type = "Nedelec 1st kind H(curl)"
> element_order = 1
> V_c = FunctionSpace(mesh, element_type, element_order)
> V = MixedFunctionSpace([V_c, V_c])
> Ef = TrialFunctions(V)
> vf = TestFunctions(V)
> Er, Ei = Ef[0], Ef[1]
> vr, vi = vf[0], vf[1]
>
> # Mark boundaries
> boundaries = FacetFunction("uint", mesh)
> boundaries.set_all(0)
>
> class RadiationBoundary(SubDomain):
> def inside(self, x, on_boundary):
> return on_boundary
>
> radiation_bnd = RadiationBoundary()
> d_exc = 1; radiation_bnd.mark(boundaries, d_exc)
>
> # Integration over the boundaries
> ds = Measure('ds')[boundaries]
> n_v = FacetNormal(mesh)
>
> # Position or 'radius' vector
> r_v = V.cell().x
> # Constant vectors for weak form
> k_dir_c = Constant( k_dir.tolist() )
> E0_c = Constant( E0.tolist() )
>
> pde = dot(k0*cross(n, cross(n, Ei)) + k0*sin( k0*dot(k_dir_c, r_v)
> )*cross(n_v, cross(E0_c, k_dir_c-n_v)), vr)*ds(d_exc) + dot(-k0*cross(n,
> cross(n, Ei)) + k0*cos( k0*dot(k_dir_c, r_v) )*cross(n_v, cross(E0_c,
> k_dir_c-n_v)), vi)*ds(d_exc)
>
> # Thank you!
>

It seems not to work in the expected way. I have omitted unrelated code, but this is a 3D time-harmonic Maxwell simulation, with the boundary condition implying excitation with a plane wave.

In the solution generated with the code above, there are no wave solutions in response to the incident plane wave. But if I build the expressions below, it works:

r_v = Expression(['x[{0:d}]'.format(i) for i in range(3)], element=V_c.ufl_element()) #instead of V.cell().x
k_dir_c = Expression(['{0:g}'.format(k_dir[i]) for i in range(3)], element=V_c.ufl_element() ) #instead of Constant( k_dir.tolist() )
E0_c = Expression(['{0:g}'.format(E0[i]) for i in range(3)], element=V_c.ufl_element() ) #instead of Constant( E0.tolist() )

I came up with this after posting the question. What do you think?

On Wed, Jun 13, 2012 at 07:45:50PM -0000, Luis Linares wrote:
> Question #200331 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/200331
>
> Status: Answered => Solved
>
> Luis Linares confirmed that the question is solved:
> It seems not to work in the expected way. I have omitted unrelated code,
> but this is a 3D time-harmonic Maxwell simulation, with the boundary
> condition implying excitation with a plane wave.
>
> In the solution generated with the code above, there are no wave
> solutions in response to the incident plane wave. But if I build the
> expressions below, it works:
>
> r_v = Expression(['x[{0:d}]'.format(i) for i in range(3)], element=V_c.ufl_element()) #instead of V.cell().x
> k_dir_c = Expression(['{0:g}'.format(k_dir[i]) for i in range(3)], element=V_c.ufl_element() ) #instead of Constant( k_dir.tolist() )
> E0_c = Expression(['{0:g}'.format(E0[i]) for i in range(3)], element=V_c.ufl_element() ) #instead of Constant( E0.tolist() )
>
> I came up with this after posting the question. What do you think?

Does it work if you just enter in the values manually?

k_dir_c = Constant((0, 0, -1))

?

--
Anders

Yes, with works with the manual values. It also seems indifferent to either declaration of the position vector r_v:

r_v = Expression(['x[{0:d}]'.format(i) for i in range(3)], element=V_c.ufl_element())
# or
r_v = instead of V.cell().x

Seem to work.