When set numerical function in variational form got error: 'shape mismatch' or 'invalid ranks'

Question

Hi, everyone,

I set term G and f in variational form as two functions, because in future, I will assign other value of G and f at each node based on calculation result u. I applied these two terms in variational form. For term f, when it product with v, I can write as inner(f, v). However, for term G, when I product G with sym(grad(u)), it always does not work and I tried different method then got errors mentioned in title.

I write a small sample only included my problem.

from dolfin import *

import numpy

parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["optimize"] = True

# Define mesh
mesh = UnitSquareMesh(3, 3)

# Define vector function space
V = VectorFunctionSpace(mesh, "CG", 1)

u = TrialFunction(V)
v = TestFunction(V)

G = interpolate(Expression(('1', '1')), V)
plot(G, title='G')

f = interpolate(Expression(('1', '1')), V)
# f = Expression(('1', '1'))

a = G*inner(sym(grad(u)), grad(v))*dx  # does not works

L = inner(f, v)*dx

bc = DirichletBC(V, Constant((0.0, 0.0)), DomainBoundary())

u = Function(V)

solve(a == L, u, bc)

plot(u, title='solution')

Comments

G is a vector (i.e. a tensor of rank 1), while sym(grad(u)), grad(v) are matrices (i.e. tensor of rank 2).

The two weak forms you wrote

a = inner(G*sym(grad(u)), grad(v))*dx  # does not works
a = inner(inner(G, sym(grad(u))), grad(v))*dx  # does not works

are not dimensionally correct.

For

a = inner(G*sym(grad(u)), grad(v))*dx

to make sense, G should either be a scalar or a matrix.

Could you explain the problem you want to solve? In particular please provide the weak form you want to discretize.

---

Another thing that I am curious about is why do you not want to use the non-linear solver considering it to be a non-linear problem. You say that it is a linear problem because you use the current strain to update the distribution. Has this method been used by anyone in past? Is there a reference or source for this?

---

Hi, Legolas,
See you again!
This problem can also be solved by nonlinear solver.
My first thought is I can solve my problem with the concept of iteration. Just like what I always do to solve a finite difference problem with explicit scheme. It is just my intuitive.
I will try both.
Thank you for the reminding!

Answer 1

I found a way to solve my problem. I got this solution from FEniCS course, lecture 3: static nonlinear PDEs.
If I want to set a parameter invariational form as a function of

1. Unknown
2. Result in previous step
3. Other parameters which is different at different nodes.

I should set this parameter as function.
In the example above, use G as an example.
I should not set G as:

G = interpolate(Expression(('1', '1')), V)

I should set G as a function, for example it is a function of unknown:

def G(u):
    return G0*inner(u, u)**0.5

Then it works.