### Dot product requires non-scalar arguments, got arguments with ranks 0 and 2.

212
views
0
5 months ago by
Hello All,

I am trying to implement Vector Laplacian with fenics.

Governing Equation:

- Grad(div(u)) + curl(curl(u)) = f in OMEGA ,  $u.n=0$u.n=0 ,  $curlu\times n=0$curlu×n=0

Mixed Formulation:

$-div\left(u\right)=\sigma$div(u)=σ   ,   $grad\left(\sigma\right)+curl\left(curl\left(u\right)\right)=f$grad(σ)+curl(curl(u))=ƒ   in OMEGA

Weak Form:

$\int\sigma.\mu$σ.μ  dx  -   $\int u.grad\left(\mu\right)$u.grad(μ)  dx = 0   ,  $\mu$ μ  in    $H^1\left(OMEGA\right)$H1(OMEGA)
$\int grad\left(\sigma\right).v$grad(σ).v  dx +  $curl\left(u\right).curl\left(v\right)$curl(u).curl(v) dx  =  $\int f.v$ƒ .v   ,   $v$v  in    $H^c\left(OMEGA\right)$Hc(OMEGA)

from dolfin import *

# Create mesh
mesh = UnitSquareMesh(32, 32)

# Define function spaces

NED = FiniteElement("N1curl", mesh.ufl_cell(), 1)
CG  = FiniteElement("CG", mesh.ufl_cell(), 1)

W = NED * CG
V = FunctionSpace(mesh,W)

# Define trial and test functions

(sigma, u) = TrialFunctions(V)
(mio, v) = TestFunctions(V)

I have used the bilinear form as follows:

a = (dot(sigma, mio) - dot(u, grad(mio))+ dot(grad(sigma), v) + inner(cursl(u),curl(v)))*dx
L = f*v*dx​

But I have seen this error :

Dot product requires non-scalar arguments, got arguments with ranks 0 and 2.

I think the error roots in my bilinear form.

Please let me know how can I improve this line?

Thanks,
Ali

Community: FEniCS Project

7
5 months ago by
In the formulation you are using the primary variable u is a vector field, to be approximated by an element of the Nedelec space, and the auxiliary variable is   sigma = -div(u),  which is a scalar function to be approximated by an element of the Lagrange (CG) space.  Thus you should have defined W as CG * NED, not NED * CG.  Here is a working code.
from dolfin import *
# Create mesh
mesh = UnitSquareMesh(32, 32)
# Define function spaces
NED = FiniteElement("N1curl", mesh.ufl_cell(), 1)
CG = FiniteElement("CG", mesh.ufl_cell(), 1)
W = CG * NED
V = FunctionSpace(mesh,W)
# Define trial and test functions
(sigma, u) = TrialFunctions(V)
(tau, v) = TestFunctions(V)
# Define bilinear form and right hand side
f = Constant((1., 2.))
L = dot(f, v)*dx
# solve and plot
w = Function(V)
solve(a == L, w)
sigma, u = w.split()
plot(u)
import matplotlib.pyplot as plt
plt.show()
​
Dear Professor Arnold,

Thank you so much.

I am trying to see how the error will change by using polynomials with different degrees, (eg, 1, 2 3), so based on

f = Constant((1., 2.))​

I have defined an exact solution:

u_D = Expression(('0.5*pow(x[0],2)','pow(x[1],2)'),degree=2)
Thus our L_2 error will be:

error_L2 = errornorm(u_D, u, 'L2')​

Then by changing the degree of NED and CG elements and mesh size, the error does not change ('error_L2 =', 0.40631205910602214).

Thanks,
Ali

written 5 months ago by Ali
3
The problem is that your u_D is not the exact solution: it does not satisfy the boundary condition u.n = 0.  Change it to
u_D = Expression(('0.5*(x[0] - pow(x[0], 2))','x[1] - pow(x[1], 2)'),degree=2)
​

and everything will work.  With the same 32 x 32 mesh, the error is 0.0051.  Changing to degree 2 elements it drops to 0.000081.  Since the solution is a polynomial of degree 2, changing to degree 3 elements reduces the error to round-off (10^{-12}), since NED3 then contains the solution.

written 5 months ago by Douglas N Arnold