Vector field within subdomains

Errors
1

Hello,

I am working on a mesh made up of two sub-domains.
I would like to define a vector field on the whole domain such as:
q = [0, 0] on subdomain 1 and q = [5, 0] on subdomain 2.

I have tried to implement this in the following way but I have an error, and obviously the way I am trying to do it is not appropriate.

Q = FunctionSpace(mesh, ("DG", 0))
q = Function(Q)
subdomain1 = (cell_tags.values == 5)
q.x.array[subdomain1, :] = np.full((subdomain1.sum(), 2), [0, 0])
subdomain2 = (cell_tags.values == 6)
q.x.array[subdomain2, :] = np.full((subdomain2.sum(), 2), [5, 0])

The error is:

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
Input In [17], in <cell line: 5>()
      3 q = Function(Q)
      4 subdomain1 = (cell_tags.values == 5)
----> 5 q.x.array[subdomain1, :] = np.full((subdomain1.sum(), 2), [0, 0])
      6 subdomain2 = (cell_tags.values == 6)
      7 q.x.array[subdomain2, :] = np.full((subdomain2.sum(), 2), [5, 0])

IndexError: too many indices for array: array is 1-dimensional, but 2 were indexed

Do you have any idea how to do this?

Thanks a lot!

2

You should give a look at Defining subdomains for different materials — FEniCSx tutorial
You have replace subdomain1.sum() with len(subdomain1)

1 Like
3

Hello,

Thank you for your reply @RaphaelC, but unfortunately I don’t think that’s where the problem lies. Indeed, even after modifying the code as suggested:

Q = FunctionSpace(mesh, ("DG", 0))
q = Function(Q)
subdomain1 = cell_tags.indices[cell_tags.values==5]
q.x.array[subdomain1] = np.full((len(subdomain1),2),  [0, 0])
subdomain2 = cell_tags.indices[cell_tags.values==6]
q.x.array[subdomain2]  = np.full((len(subdomain2),2), [5, 0])

I still have an error:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Input In [39], in <cell line: 5>()
      3 q = Function(Q)
      4 subdomain1 = cell_tags.indices[cell_tags.values==5]
----> 5 q.x.array[subdomain1] = np.full((len(subdomain1),2),  [0, 0])
      6 subdomain2 = cell_tags.indices[cell_tags.values==6]
      7 q.x.array[subdomain2]  = np.full((len(subdomain2),2), [5, 0])

ValueError: shape mismatch: value array of shape (5336,2) could not be broadcast to indexing result of shape (5336,)

The problem comes from the fact that q.x.array is a vector and not a two-column matrix, so it is not possible to define a vector field in this way.

Do you have any other ideas for solving this problem?

Thank you!

4

Oh, yes, you are matching a float against a vector .
Either Q is a VectorFunctionSpace not a FunctionSpace or tag should be float.

1 Like
5

Yes, that’s right, thank you!

For those who would like to know, I fixed my code as follows:

Q = VectorFunctionSpace(mesh, ("DG", 0))
q = Function(Q)
num_cells = mesh.topology.index_map(mesh.topology.dim).size_local
block_size = Q.dofmap.index_map_bs 
for i in range(num_cells):
    if cell_tags.values[i] == 5:
        q.x.array[[i*block_size, i*block_size+1]] = [0, 0]
    elif cell_tags.values[i] == 6:
        q.x.array[[i*block_size, i*block_size+1]] = [5, 0]