What's wrong for my program dimension?

dolfinx
1

I encountered the error message below. I am wondering is there a problem with the dimension of dx and ds not being the same. I would greatly appreciate any suggestions or guidance on how to resolve this issue.

Can't add expressions with different shapes.
ERROR:UFL:Can't add expressions with different shapes.
Traceback (most recent call last):
  File "/Users/victoriachan/Documents/CityU Doc/FYP/python_files/ProgramV2.py", line 61, in <module>
    + ufl.inner(w_h - ufl.grad(u_h), ufl.grad(kappa * ufl.nabla_div(p))) * ds \
  File "/Users/victoriachan/opt/anaconda3/envs/fenicsdolxconda/lib/python3.10/site-packages/ufl/exproperators.py", line 223, in _sub
    return Sum(self, -o)
  File "/Users/victoriachan/opt/anaconda3/envs/fenicsdolxconda/lib/python3.10/site-packages/ufl/algebra.py", line 41, in __new__
    error("Can't add expressions with different shapes.")
  File "/Users/victoriachan/opt/anaconda3/envs/fenicsdolxconda/lib/python3.10/site-packages/ufl/log.py", line 135, in error
    raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: Can't add expressions with different shapes.

The source code (without the import packages):

Screenshot 2023-03-23 at 2.42.07 AM
Screenshot 2023-03-23 at 2.42.07 AM1360×1306 162 KB

2

If you want to receive help, you need to supply a minimal code (in text form) that reproduces the error.
The code should be wrapped as

```python
import dolfinx
#add rest of code here

```

General comment to your code. If you want to solve your problem for both V and W, you need to use a MixedElement. See:
https://docs.fenicsproject.org/dolfinx/v0.6.0/python/demos/demo_stokes.html#non-blocked-direct-solver
or any of the approaches further up if you want a more technical (But better perfoming) code.

1 Like
3

Thanks. I am reading this example. Could I ask how can I change the boundary to a unit square boundary for the biharmonic equation? As I don’t need the velocity term now and I don’t know how to revise it. Thanks for your kind help.

The original boundary code from the link:

# Create mesh
msh = create_rectangle(MPI.COMM_WORLD,
                       [np.array([0, 0]), np.array([1, 1])],
                       [32, 32],
                       CellType.triangle, GhostMode.none)


# Function to mark x = 0, x = 1 and y = 0
def noslip_boundary(x):
    return np.logical_or(np.logical_or(np.isclose(x[0], 0.0),
                                       np.isclose(x[0], 1.0)),
                         np.isclose(x[1], 0.0))


# Function to mark the lid (y = 1)
def lid(x):
    return np.isclose(x[1], 1.0)


# Lid velocity
def lid_velocity_expression(x):
    return np.stack((np.ones(x.shape[1]), np.zeros(x.shape[1])))

# No-slip boundary condition for velocity field (`V`) on boundaries
# where x = 0, x = 1, and y = 0
noslip = np.zeros(msh.geometry.dim, dtype=PETSc.ScalarType)
facets = locate_entities_boundary(msh, 1, noslip_boundary)
bc0 = dirichletbc(noslip, locate_dofs_topological(V, 1, facets), V)

# Driving velocity condition u = (1, 0) on top boundary (y = 1)
lid_velocity = Function(V)
lid_velocity.interpolate(lid_velocity_expression)
facets = locate_entities_boundary(msh, 1, lid)
bc1 = dirichletbc(lid_velocity, locate_dofs_topological(V, 1, facets))

# Collect Dirichlet boundary conditions
bcs = [bc0, bc1]
4

I am not sure what you are asking for here.