An error in assembling a variational problem

Question

Hi, I'm building a program to solve PDEs, I got stuck in assembling a variational problem.
here is a part of my code that causes the problem:

from dolfin import*
import numpy as np
import scipy as sc

mm = 1e-3
Dx = 100 * mm
Dy = 100 * mm
l = 2.1 * mm
e = 1.97 * mm

cond_x = 6.8e-8
cond_y = 9.9e-8

mesh = BoxMesh (-Dx/2.,-Dy/2.,-l,Dx/2,Dy/2,-l-e,10,10,2)

# Create mesh functions for c00, c01, c11, c22
c00 = MeshFunction("double", mesh, 3)
c01 = MeshFunction("double", mesh, 3)
c11 = MeshFunction("double", mesh, 3)
c22 = MeshFunction("double", mesh, 3)

for cell in cells(mesh):
    if (-Dx/2. < cell.midpoint().x() < Dx/2.) and \
       (-Dy/2. < cell.midpoint().y() < Dy/2.) and \
       (-l-e < cell.midpoint().z() < -l):
        c00[cell] = cond_x
        c01[cell] = 0.0
        c11[cell] = cond_y
        c22[cell] = cond_y

# Code for C++ evaluation of conductivity
conductivity_code = """
class Conductivity : public Expression
{
public:

  // Create expression with 4 components
  Conductivity() : Expression(4) {}

  // Function for evaluating expression on each cell
  void eval(Array<double>& values, const Array<double>& x, const ufc::cell& cell) const
  {
    const uint D = cell.topological_dimension;
    const uint cell_index = cell.index;
    values[0] = (*c00)[cell_index];
    values[1] = (*c01)[cell_index];
    values[2] = (*c11)[cell_index];
    values[3] = (*c22)[cell_index];
  }

  // The data stored in mesh functions
  std::shared_ptr<MeshFunction<double> > c00;
  std::shared_ptr<MeshFunction<double> > c01;
  std::shared_ptr<MeshFunction<double> > c11;
  std::shared_ptr<MeshFunction<double> > c22;

};
"""

c = Expression(cppcode=conductivity_code)
c.c00 = c00
c.c01 = c01
c.c11 = c11
c.c22 = c22
C = as_tensor(((c[0], c[1], c[1]), \
               (c[1], c[2], c[1]), \
               (c[1], c[1], c[3])))

# Define function spaces
V = VectorFunctionSpace(mesh, "CG", 1)
Q = FunctionSpace(mesh, "CG", 1)
W = MixedFunctionSpace([V, Q])

(u_r, p_r) = TrialFunctions(W)
(v_r, q_r) = TestFunctions(W)

a_r = inner(C*div(u_r) , div(v_r))*dx

A = assemble(a_r)

When I execute I get this following errors:

Shape mismatch.

Traceback (most recent call last):
  File "/home/chaouki/Documents/Exemples_FEniCS/exe01_ask_a_q.py", line 77, in <module>
    a_r = inner(C*div(u_r) , div(v_r))*dx
  File "/usr/lib/python2.7/dist-packages/ufl/operators.py", line 130, in inner
    return Inner(a, b)
  File "/usr/lib/python2.7/dist-packages/ufl/tensoralgebra.py", line 156, in __new__
    ufl_assert(ash == bsh, "Shape mismatch.")
  File "/usr/lib/python2.7/dist-packages/ufl/assertions.py", line 37, in ufl_assert
    if not condition: error(*message)
  File "/usr/lib/python2.7/dist-packages/ufl/log.py", line 151, in error
    raise self._exception_type(self._format_raw(*message))
UFLException: Shape mismatch.

I didn't understand what does the error massage mean, what should I do?

Answer 1

Well, the inner(.,.) is defined for objects of the same shape. C*div(u_r) is 3x3 object and div(v_r) is 1x1. What exactly do you want to accomplish with your computation?

Comments

The final result should be a matrix.
And if I compute this :

a_r = inner(div(u_r) , div(v_r))*dx

it yields this without any problem:

A.array()
array([[ 0.00016417,  0.        ,  0.00016417, ...,  0.        ,
         0.        ,  0.        ],
       [ 0.        ,  0.0004925 , -0.00016417, ...,  0.        ,
         0.        ,  0.        ],
       [ 0.00016417, -0.00016417,  0.00065667, ...,  0.        ,
         0.        ,  0.        ],
       ..., 
       [ 0.        ,  0.        ,  0.        , ...,  0.00032833,
         0.00166667,  0.        ],
       [ 0.        ,  0.        ,  0.        , ...,  0.00166667,
         0.01692047,  0.        ],
       [ 0.        ,  0.        ,  0.        , ...,  0.        ,
         0.        ,  0.        ]])

But with tensor C :

a_r = inner(C*div(u_r) , div(v_r))*dx

The problem happens.

---

The result you posted is of course a matrix. It contains values of integral of divergence of basis functions (that forms functions u_r and v_r) on each element. But that is something different.
First,

inner(div(u_r),div(v_r))*dx = div(u_r) *div(v_r) * dx`

since div(u_r), div(v_r) are scalar functions.
Second, once again, you cannot use inner product for differently shaped objects. You can write a_r = C * inner(div(u_r) , div(v_r))*dx which will not trigger the error you are talking about. But now you get an error that UFL can only integrate scalar expressions.

So, once again, what is your goal? Integrate a tensor form? Then I guess you have to do it component-wise, i.e. something like

a_r[i] = c[i] * div(u_r) * div(v_r) * dx
A[i] = assemble(a_r[i])