### A problem with the weak form

121

views

0

Hi guys,

So I have this chemotaxis term in my integral given by

$\nabla.\left(M\nabla P\right)$∇.(

`M`∇

`P`)

with the Neumann boundary condition:

$\frac{\partial P}{\partial n}=0$∂

`P`∂

`n`=0

I tried the weak form

$\int\left(\nabla M.\nabla P\right)v\text{ dx}-\int\nabla P.\nabla\left(Mv\right)\text{ dx}$∫(∇

`M`.∇

`P`)

`v`dx−∫∇

`P`.∇(

`M`

`v`) dx

with the formula

`dot(grad(M),grad(P))*v*dx-dot(grad(P),grad(M*v))*dx`

but I receive the error "unsupported operand type(s) for *: 'Measure' and 'Indexed'Can someone tell me what is wrong?

Here is the complete code:

```
from dolfin import *
import numpy as np
Domain=Mesh('general_boundary_wenrui+5.xml')
Bulk = MeshFunction('size_t' , Domain , 'general_boundary_wenrui+5_physical_region.xml' ) Boundary = MeshFunction('size_t', Domain , 'general_boundary_wenrui+5_facet_region.xml')
ds = Measure('ds', subdomain_data=Boundary)
dx = Measure('dx',subdomain_data= Bulk)
Dh = 3.93
Dl = 29.89
Dr = 0.205
Dlox = 29.89
Dm = 8.64e-7
Dp = 17.28
kl = 2.35e-4
kh = 5.29e-6
L0 = 9
r0 = 0.26
H0 = 7
alphaH = 1
alphaL = 1
alphaM = 0.2
Lambda = 10
LambdaPE = 8.65e-10
dM = 0.015
dP = 1.73
M0 = 5e-5
kaic = 1
KLox = 0.5
P1 = FiniteElement('P', tetrahedron, 1)
element = MixedElement([P1, P1, P1, P1, P1, P1])
Mixed_Space = FunctionSpace(Domain, element)
v1 , v2, v3, v4, v5, v6=TestFunction(Mixed_Space)
U = Function(Mixed_Space)
L, H, r, Lox, M, PP = split(U)
F=kl*r*L*v1*dx-Dl*alphaL*(L0-L)*v1*ds(1)+Dl*dot(grad(L),grad(v1))*dx\
+kh*r*H*v2*dx-Dh*alphaH*(H0-H)*v2*ds(1)+Dh*dot(grad(H),grad(v2))*dx\
+r*(kl*L+kh*H)*v3*dx-r0*v3*dx+Dr*dot(grad(r),grad(v3))*dx\
+Lambda*M0*Lox*v4*dx-kl*r*L*v4*dx+Dlox*dot(grad(Lox),grad(v4))*dx\
+ dM*M*v5*dx-Dm*alphaM*M0*v5*ds(1)+Dm*alphaM*M*v5*ds(1)+Dm*dot(grad(M),grad(v5))*dx\
+kaic*M*dot(grad(PP),grad(v5))*dx\
+ dP*PP*v6*dx-LambdaPE*(Lox/(KLox+Lox))*v6*dx+Dp*dot(grad(PP),grad(v6))*dx
solve(F == 0, U, solver_parameters={"newton_solver": {"relative_tolerance": 1e-12,
"absolute_tolerance": 1e-12,"maximum_iterations": 500}})
```

and the mesh files are:File attached: general_boundary_wenrui+5.xml (1.36 MB)

File attached: general_boundary_wenrui+5_facet_region.xml (1.1 MB)

File attached: general_boundary_wenrui+5_physical_region.xml (509.85 KB)

Community: FEniCS Project

Please login to add an answer/comment or follow this question.

M∇Pv.n|_{surƒ ace}−∫M∇P.∇vdx.Also try using nabla_grad(v) instead of grad(v) for the volume integral term.

If you are still not able to solve the issue, try inner() instead of dot(). Essentially both do the same thing, but I use inner(). Maybe it will work for you too.

Why would this not be your form (Poisson's eqn., M scalar or tensor):

`a = inner(grad(v), M*grad(p))*dx`