# Question about how to define a function

**URL:** <https://fenicsproject.discourse.group/t/question-about-how-to-define-a-function/14242>\
**Category:** General\
**Created:** [April 13, 2024, 3:43pm UTC](https://fenicsproject.discourse.group/t/question-about-how-to-define-a-function/14242 "2024-04-13T15:43:56Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![guangyi\_cai](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/guangyi_cai/32/6508_2.png) [@guangyi\_cai](https://fenicsproject.discourse.group/u/guangyi_cai)\
**Post date:** [April 13, 2024, 3:43pm UTC](https://fenicsproject.discourse.group/t/question-about-how-to-define-a-function/14242/1 "2024-04-13T15:43:56Z")

</div>

Im a rookie in Fenicsx, Im confused about how to define the function phi\_0 in picture 3 so that I can put phi\_0 into a ( u , v ) and do some advanced work.

 ![1](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/3/3ac5f5bbc91f4b20692a78054a57a112ffeaa4ee.png)  
 ![3](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/4/4d225ca9dcb0daf07d4ecf5d6b7145c3eab78f4a.png)  
 ![2](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/9/9d32a60966341c21d7f421d4e043aa8166a31af9.png)

my code is as followed:

```auto
# create a mesh
msh = create_rectangle(
    MPI.COMM_WORLD, [np.array([0, 0]), np.array([1, 1])], [32, 32], CellType.triangle
)

#function to mark x=0
def in_boundary(x):
    return np.isclose(x[0],0.0)

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

#function to mark x=1
def in_boundary(x):
    return np.isclose(x[0],1.0)

# in_velocity
def in_velocity_expression(x):
    return np.stack((np.ones(x.shape[1]), np.zeros(x.shape[1])))

#define V and Q
P2 = element("Lagrange", msh.basix_cell(), 2, shape=(msh.geometry.dim,))
P1 = element("Lagrange", msh.basix_cell(), 1)
V, Q = functionspace(msh, P2), functionspace(msh, P1)

# No-slip condition on boundaries where y=0,y=1
noslip = np.zeros(msh.geometry.dim, dtype=PETSc.ScalarType)
facets = locate_entities_boundary(msh, 1, wall_boundary)
bc0 = dirichletbc(noslip, locate_dofs_topological(V, 1, facets), V)

# Driving velocity condition on boundary where x=0
in_velocity = Function(V)
in_velocity.interpolate(in_velocity_expression)
facets = locate_entities_boundary(msh, 1, in_boundary)
bc1 = dirichletbc(in_velocity, locate_dofs_topological(V, 1, facets))

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

# alpha_0
alpha_0=Constant(mesh, PETSc.ScalarType(1))

# phi_0

# Define variational problem
(u, p) = ufl.TrialFunction(V), ufl.TrialFunction(Q)
(v, q) = ufl.TestFunction(V), ufl.TestFunction(Q)
f = Constant(msh, (PETSc.ScalarType(0), PETSc.ScalarType(0)))

a = form([[inner(grad(u), grad(v)) * dx + inner(alpha_0*(1-phi_0)*(1-phi_0)*u,v)*dx, -inner(p, div(v)) * dx], [inner(div(u), q) * dx, None]])
L = form([inner(f, v) * dx, inner(Constant(msh, PETSc.ScalarType(0)), q) * dx]) 

```

I would appreciate your help!

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [April 13, 2024, 4:14pm UTC](https://fenicsproject.discourse.group/t/question-about-how-to-define-a-function/14242/2 "2024-04-13T16:14:10Z")

</div>

It seems like you want phi\_0 to be a cellwise constant, i.e. Create a function in DG-0 that is 1 within the respective area and 0 elsewhere. See:  
[https://jsdokken.com/dolfinx-tutorial/chapter3/subdomains.html](https://jsdokken.com/dolfinx-tutorial/chapter3/subdomains.html)

---

<div class="post-metadata">

**Author:** ![guangyi\_cai](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/guangyi_cai/32/6508_2.png) [@guangyi\_cai](https://fenicsproject.discourse.group/u/guangyi_cai)\
**Post date:** [April 14, 2024, 10:38am UTC](https://fenicsproject.discourse.group/t/question-about-how-to-define-a-function/14242/4 "2024-04-14T10:38:21Z")

</div>

Appreciate your help！  
So the code would be:

```auto
#define phi_0
W = functionspace(msh, ("DG", 0))

def condition_0(x):
    return np.logical_or(x[0] <= 0.25, np.logical_and(x[1] >= 1/3, x[1] <= 2/3))

def condition_1(x):
    return np.logical_or(np.logical_and(x[0] > 0.25,x[1] < 1/3),np.logical_and(x[0] > 0.25,x[1] > 2/3))

phi_0 = Function(W)
cells_0 = locate_entities(msh, msh.topology.dim, condition_0)
cells_1 = locate_entities(msh, msh.topology.dim, condition_1)

phi_0.x.array[cells_0] = np.full_like(cells_0, 1, dtype=default_scalar_type)
phi_0.x.array[cells_1] = np.full_like(cells_1, 0, dtype=default_scalar_type)

# Define variational problem
(u, p) = ufl.TrialFunction(V), ufl.TrialFunction(Q)
(v, q) = ufl.TestFunction(V), ufl.TestFunction(Q)
f = Constant(msh, (PETSc.ScalarType(0), PETSc.ScalarType(0)))

a = form([[inner(grad(u), grad(v)) * dx + inner(alpha_0*(Constant(msh, PETSc.ScalarType(1))-phi_0)*(Constant(msh, PETSc.ScalarType(1))-phi_0)*u,v)*dx, -inner(p, div(v)) * dx], [inner(div(u), q) * dx, None]])
L = form([inner(f, v) * dx, inner(Constant(msh, PETSc.ScalarType(0)), q) * dx]) 

```

I’m not sure if my code is correct.

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [April 14, 2024, 10:52am UTC](https://fenicsproject.discourse.group/t/question-about-how-to-define-a-function/14242/5 "2024-04-14T10:52:24Z")

</div>

You can check for correctness by visualizing `phi_0` with either paraview or pyvista.

---

<div class="post-metadata">

**Author:** ![guangyi\_cai](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/guangyi_cai/32/6508_2.png) [@guangyi\_cai](https://fenicsproject.discourse.group/u/guangyi_cai)\
**Post date:** [April 14, 2024, 3:22pm UTC](https://fenicsproject.discourse.group/t/question-about-how-to-define-a-function/14242/7 "2024-04-14T15:22:39Z")

</div>

I have Defined variational problem:

```auto
# Define variational problem
(u, p) = ufl.TrialFunction(V), ufl.TrialFunction(Q)
(v, q) = ufl.TestFunction(V), ufl.TestFunction(Q)
f = Constant(msh, (PETSc.ScalarType(0), PETSc.ScalarType(0)))

a = form([[inner(grad(u), grad(v)) * dx + inner(alpha_0*(Constant(msh, PETSc.ScalarType(1.0))-phi_0)*(Constant(msh, PETSc.ScalarType(1.0))-phi_0)*u, v)*dx, -inner(p, div(v)) * dx], [inner(div(u), q) * dx, None]])
L = form([inner(f, v) * dx, inner(Constant(msh, PETSc.ScalarType(0)), q) * dx]) 

```

How can I solve the PDE in picture 1 with Taylor-Hood element,  
is there some references for me?

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [April 14, 2024, 5:49pm UTC](https://fenicsproject.discourse.group/t/question-about-how-to-define-a-function/14242/8 "2024-04-14T17:49:17Z")

</div>

It is unclear what your current issue is? What happens if you follow the demo  
[https://docs.fenicsproject.org/dolfinx/v0.7.3/python/demos/demo\_stokes.html](https://docs.fenicsproject.org/dolfinx/v0.7.3/python/demos/demo_stokes.html)

Do you get an error message.?

---

<div class="post-metadata">

**Author:** ![guangyi\_cai](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/guangyi_cai/32/6508_2.png) [@guangyi\_cai](https://fenicsproject.discourse.group/u/guangyi_cai)\
**Post date:** [April 15, 2024, 4:28am UTC](https://fenicsproject.discourse.group/t/question-about-how-to-define-a-function/14242/9 "2024-04-15T04:28:12Z")

</div>

I have a few questions about the demo：

1.we can see “-inner(p, div(v)) \* dx” in the weak formulation given in the document,  
and the code from the demo:

```auto
a = form([[inner(grad(u), grad(v)) * dx, inner(p, div(v)) * dx],
          [inner(div(u), q) * dx, None]])

```

why there is not a negative sign in front of " inner(p, div(v)) \* dx" ?

1. Under this boundary conditions, I get the right output:

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

# Function to mark the lid (x = 0)
def lid(x):
    return np.isclose(x[0], 0.0)

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

# No-slip condition on boundaries where y = 0, y = 1, x = 1
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 (lid) velocity condition on top boundary (x = 0)
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]

```

But once I want to remove the dirichlet boundary condition of x = 1 and I change the code into:

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

# Function to mark the lid (x = 0)
def lid(x):
    return np.isclose(x[0], 0.0)

```

(the rest parts of the code is as before)

I get an issue:

```auto
Traceback (most recent call last):
  File "/home/fenics/work/demo3.py", line 226, in <module>
    norm_u_0, norm_p_0 = nested_iterative_solver()
  File "/home/fenics/work/demo3.py", line 157, in nested_iterative_solver
    assert nsp.test(A)
AssertionError

```

I think I set the boundary conditions wrongly but I don’t know how to modify it.

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [April 15, 2024, 5:04am UTC](https://fenicsproject.discourse.group/t/question-about-how-to-define-a-function/14242/10 "2024-04-15T05:04:00Z")

</div>

> [@guangyi\_cai](#):
>
> why there is not a negative sign in front of " inner(p, div(v)) \* dx" ?

One usually solves for the negative pressure to ensure a symmetric system.

> [@guangyi\_cai](#):
>
> But once I want to remove the dirichlet boundary condition of x = 1 and I change the code into:
> 
> ```auto
> # Function to mark y = 0, y = 1
> def noslip_boundary(x):
> return np.logical_or(np.isclose(x[1], 0.0), np.isclose(x[1], 1.0))
> 
> # Function to mark the lid (x = 0)
> def lid(x):
> return np.isclose(x[0], 0.0)
> 
> ```

Just to make sure, here you are trying to enforce  
u(0,y)=(1,0),  
u(x,0)=u(x,1)=(0,0)  
-du/dn+pn=0 at (1,y), i.e. A square channel with (1,0) at the left inlet, outlet at right boundary and top and bottom being zero?

(Im currently not at a computer, so I cannot test your code).

Is the error messages stemming from you modifying the referenced demo or your code above ?

---

<div class="post-metadata">

**Author:** ![guangyi\_cai](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/guangyi_cai/32/6508_2.png) [@guangyi\_cai](https://fenicsproject.discourse.group/u/guangyi_cai)\
**Post date:** [April 16, 2024, 2:36am UTC](https://fenicsproject.discourse.group/t/question-about-how-to-define-a-function/14242/12 "2024-04-16T02:36:18Z")

</div>

The boundary conditions you mentioned is exactly what I want to set up!  
And I just changed the code snippet from the referenced demo **demo\_stokes.py** given in the document

```auto
# 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)

```

into:

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

# Function to mark the lid (x = 0)
def lid(x):
    return np.isclose(x[0], 0.0)

```

I tried to run the demo I modified and then I got the error messages I mentioned above.

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [April 16, 2024, 11:10am UTC](https://fenicsproject.discourse.group/t/question-about-how-to-define-a-function/14242/13 "2024-04-16T11:10:58Z")

</div>

The error message is for testing the null space.  
Simply removing the null space (as you no longer have a nullspace, as you have a natural boundary condition. I.e. remove

```python
    nsp = PETSc.NullSpace().create(vectors=[null_vec])
    assert nsp.test(A)
    A.setNullSpace(nsp)

```

this from the two places it is used.

---

<div class="post-metadata">

**Author:** ![guangyi\_cai](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/guangyi_cai/32/6508_2.png) [@guangyi\_cai](https://fenicsproject.discourse.group/u/guangyi_cai)\
**Post date:** [April 17, 2024, 5:11pm UTC](https://fenicsproject.discourse.group/t/question-about-how-to-define-a-function/14242/15 "2024-04-17T17:11:01Z")

</div>

Appreciate your help again！  
There are still some problem:  
I want to put the **u** that was solved out in the demo into my advanced work i.e.:

 ![d85636fe7bd1732458c814fa5eaa70b](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/4/48cda4cbe81b15f04dc248c4ecb4fdf1c6a07d74.png)

And there is also a constrain condition :  
The solutions are mismatched between the velocity field and the phase field due to employing different polynomials as the basis functions on the same mesh. Projection and interpolationare simply used to transport datas:

 ![3692683c4057a5cf9dff0fef5c5dbc7](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/7/71658c6c55c5a9cce185837b5d92438e5e6405cf.png)

I am confused about how to import **u** that was solved out in the demo into the weak formulation in my first picture.

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [April 18, 2024, 6:42am UTC](https://fenicsproject.discourse.group/t/question-about-how-to-define-a-function/14242/16 "2024-04-18T06:42:01Z")

</div>

> [@guangyi\_cai](#):
>
> I am confused about how to import **u** that was solved out in the demo into the weak formulation in my first picture.

You can use adios4dolfinx, see  
[http://jsdokken.com/adios4dolfinx/README.html](http://jsdokken.com/adios4dolfinx/README.html)

---

<div class="post-metadata">

**Author:** ![guangyi\_cai](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/guangyi_cai/32/6508_2.png) [@guangyi\_cai](https://fenicsproject.discourse.group/u/guangyi_cai)\
**Post date:** [April 27, 2024, 12:01pm UTC](https://fenicsproject.discourse.group/t/question-about-how-to-define-a-function/14242/17 "2024-04-27T12:01:52Z")

</div>

I just installed adios4dolfinx 0.7.2 and tried to run the domo **writing\_functions\_checkpoint.py** downloaded from [http://jsdokken.com/adios4dolfinx/docs/writing\_functions\_checkpoint.html](http://jsdokken.com/adios4dolfinx/docs/writing_functions_checkpoint.html)  
And then i got an error:

```auto
Traceback (most recent call last):
  File "/home/fenics/work/writing_functions_checkpoint.py", line 41, in <module>
    adios4dolfinx.write_mesh(filename, mesh)
  File "/usr/local/lib/python3.10/dist-packages/adios4dolfinx/checkpointing.py", line 86, in write_mesh
    num_xdofs_local = mesh.geometry.index_map().size_local
AttributeError: 'PosixPath' object has no attribute 'geometry'

```

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [April 27, 2024, 12:53pm UTC](https://fenicsproject.discourse.group/t/question-about-how-to-define-a-function/14242/18 "2024-04-27T12:53:01Z")

</div>

Switch the order of the arguments. The webpage is based on v0.8.x of adios4dolfinx, and we changed the API to be consistent across all functions

> <https://github.com/jorgensd/adios4dolfinx/pull/72/files>
>
> Resolves #21 #63 and #69
> 
> Remove \`typing\` library completely and use \`from \_\_f…uture\_\_ import annotations\` instead.

See api changes in [Release v0.8.0 · jorgensd/adios4dolfinx · GitHub](https://github.com/jorgensd/adios4dolfinx/releases/tag/v0.8.0)

---

<div class="post-metadata">

**Author:** ![guangyi\_cai](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/guangyi_cai/32/6508_2.png) [@guangyi\_cai](https://fenicsproject.discourse.group/u/guangyi_cai)\
**Post date:** [April 29, 2024, 5:59am UTC](https://fenicsproject.discourse.group/t/question-about-how-to-define-a-function/14242/21 "2024-04-29T05:59:15Z")

</div>

I just referred to [https://jsdokken.com/checkpointing-presentation/#/10](https://jsdokken.com/checkpointing-presentation/#/10)  
and just simply add

```auto
#save u with adios4dolfinx
    checkpoint_file = pathlib.Path("function_checkpoint.bp")
    adx.write_mesh(msh, checkpoint_file, engine="BP4")
    adx.write_function(u, checkpoint_file, engine="BP4")

```

behind

```auto
# Create finite element {py:class}`Function <dolfinx.fem.Function>`s
    # for the velocity (on the space `V`) and for the pressure (on the
    # space `Q`). The vectors for `u` and `p` are combined to form a
    # nested vector and the system is solved.
    u, p = Function(V), Function(Q)
    x = PETSc.Vec().createNest([la.create_petsc_vector_wrap(u.x),
                                la.create_petsc_vector_wrap(p.x)])
    ksp.solve(b, x)

```

in **demo\_stokes.py**  
downloaded from[https://docs.fenicsproject.org/dolfinx/v0.7.3/python/demos/demo\_stokes.html#](https://docs.fenicsproject.org/dolfinx/v0.7.3/python/demos/demo_stokes.html#):

```auto
# Create finite element {py:class}`Function <dolfinx.fem.Function>`s
    # for the velocity (on the space `V`) and for the pressure (on the
    # space `Q`). The vectors for `u` and `p` are combined to form a
    # nested vector and the system is solved.
    u, p = Function(V), Function(Q)
    x = PETSc.Vec().createNest([la.create_petsc_vector_wrap(u.x),
                                la.create_petsc_vector_wrap(p.x)])
    ksp.solve(b, x)

    #save u with adios4dolfinx
    checkpoint_file = pathlib.Path("function_checkpoint.bp")
    adx.write_mesh(msh, checkpoint_file, engine="BP4")
    adx.write_function(u, checkpoint_file, engine="BP4")

```

I’m not sure if my operation is correct (since i have seen many very complex codes in other demos)

Furthermore,if i want to import **u** into my advanced work ,can i just simply use the code here？

```auto
from mpi4py import MPI
import adios4dolfinx
from pathlib import Path
import dolfinx
checkpoint_file = Path("function_checkpoint.bp")
msh = adios4dolfinx.read_mesh(MPI.COMM_WORLD, checkpoint_file, engine='BP4', ghost_mode=dolfinx.mesh.GhostMode.shared_facet)

from basix.ufl import element
from dolfinx.fem import functionspace, Function
P2 = element("Lagrange", msh.basix_cell(), 2, shape=(msh.geometry.dim,))
V = functionspace(msh, P2)
u = Function(V)
adios4dolfinx.read_function(u, checkpoint_file, engine="BP4")

```

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [April 29, 2024, 2:59pm UTC](https://fenicsproject.discourse.group/t/question-about-how-to-define-a-function/14242/22 "2024-04-29T14:59:32Z")

</div>

Yes, that would be the general idea.

---

<div class="post-metadata">

**Author:** ![guangyi\_cai](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/guangyi_cai/32/6508_2.png) [@guangyi\_cai](https://fenicsproject.discourse.group/u/guangyi_cai)\
**Post date:** [April 29, 2024, 4:58pm UTC](https://fenicsproject.discourse.group/t/question-about-how-to-define-a-function/14242/23 "2024-04-29T16:58:41Z")

</div>

Appreciate for your patience！  
Now I want to use projection to transport data (Maybe from P2 to P1) ,so that I can complete the follow-up calculation of the weak formulation below：

 ![48cda4cbe81b15f04dc248c4ecb4fdf1c6a07d74_2_1035x505](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/1/106672c82eee4c1626007d968d8b715b74cf9c21.png)  
Im confused about how to do the projection work.Here is my code:

```auto
#read mesh
checkpoint_file = Path("function_checkpoint.bp")
msh = adios4dolfinx.read_mesh(MPI.COMM_WORLD, checkpoint_file, engine='BP4', ghost_mode=dolfinx.mesh.GhostMode.shared_facet)

#read u
P2 = element("Lagrange", msh.basix_cell(), 2, shape=(msh.geometry.dim,))
U = functionspace(msh, P2)
u = Function(U)
adios4dolfinx.read_function(u, checkpoint_file, engine="BP4")

P1 = element("Lagrange", msh.basix_cell(), 1)
V = functionspace(msh, P1)

#projection
u_projected = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
a_proj = inner(u_projected, v) * dx
L_proj = inner(u, v) * dx

A_proj = assemble_matrix(a_proj)
A_proj.assemble()
b_proj = assemble_vector(L_proj)
solver_proj = PETSc.KSP().create(msh.comm)
solver_proj.setOperators(A_proj)
solver_proj.solve(b_proj, u_projected.vector)
u_projected.x.scatter_forward()

```

And it failed with:

```auto
Traceback (most recent call last):
  File "/home/fenics/work/phase.py", line 33, in <module>
    L_proj = inner(u, v) * dx
  File "/usr/local/lib/python3.10/dist-packages/ufl/operators.py", line 167, in inner
    return Inner(a, b)
  File "/usr/local/lib/python3.10/dist-packages/ufl/tensoralgebra.py", line 162, in __new__
    raise ValueError(f"Shapes do not match: {ufl_err_str(a)} and {ufl_err_str(b)}")
ValueError: Shapes do not match: <Coefficient id=140108284108736> and <Argument id=140108318196832>

```

It seems that there is a shape mismatch between the operands in the `inner` function …? How should I modify the code ?  
And I alse referred to [L2 projection of fine mesh solution to a coarse mesh](https://fenicsproject.discourse.group/t/projection-between-meshes/457/5) ,maybe I should interpolate the function **u** first? If so ,how should I interpolate the function?

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [April 29, 2024, 5:00pm UTC](https://fenicsproject.discourse.group/t/question-about-how-to-define-a-function/14242/24 "2024-04-29T17:00:24Z")

</div>

> [@guangyi\_cai](#):
>
> ```auto
> P1 = element("Lagrange", msh.basix_cell(), 1)
> V = functionspace(msh, P1)
> 
> ```

You need to supply a shape to: i.e. make it a vector space like in

> [@guangyi\_cai](#):
>
> ```auto
> P2 = element("Lagrange", msh.basix_cell(), 2, shape=(msh.geometry.dim,))
> U = functionspace(msh, P2)
> 
> ```

---

<div class="post-metadata">

**Author:** ![guangyi\_cai](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/guangyi_cai/32/6508_2.png) [@guangyi\_cai](https://fenicsproject.discourse.group/u/guangyi_cai)\
**Post date:** [May 3, 2024, 4:43pm UTC](https://fenicsproject.discourse.group/t/question-about-how-to-define-a-function/14242/25 "2024-05-03T16:43:05Z")

</div>

Let’s return to the discussion of solving stokes equation,now I want to Update the velocity field **u^(n+1)** by the following state equation using the known phase-field function in the previous step.  
 ![06781a74682d5a9251431b5bb788e36](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/2/2f6d277e0deb1f609fcc690fbd4833fa063a2c1e.png)  
where α(Φ) = α\_0 (1 - Φ) (1 - Φ) and the boundary conditions are:

> [@dokken](#):
>
> u(0,y)=(1,0)u(0,y)=(1,0),  
> u(x,0)=u(x,1)=(0,0)u(x,0)=u(x,1)=(0,0)u(x,0)=u(x,1)=(0,0)  
> -du/dn+pn=0−du/dn+pn=0-du/dn+pn=0 at (1,y)

I just modified **demo\_stokes.py** downloaded from [https://docs.fenicsproject.org/dolfinx/v0.7.3/python/demos/demo\_stokes.html](https://docs.fenicsproject.org/dolfinx/v0.7.3/python/demos/demo_stokes.html)

The modifications I made were:  
1.changed boundary condituons in the demo  
2.defined the Φ function  
3.set

```auto
(u, p) = ufl.TrialFunction(V), ufl.TrialFunction(Q)
(v, q) = ufl.TestFunction(V), ufl.TestFunction(Q)
f = Constant(msh, (PETSc.ScalarType(0), PETSc.ScalarType(0))) # type: ignore
a = form([[inner(grad(u), grad(v)) * dx + inner(alpha_0*(Constant(msh, PETSc.ScalarType(1.0))-phi_0)*(Constant(msh, PETSc.ScalarType(1.0))-phi_0)*u, v)*dx, inner(p, div(v)) * dx], [inner(div(u), q) * dx, None]])
L = form([inner(f, v) * dx , inner(Constant(msh, PETSc.ScalarType(0)), q) * dx]) 

```

4.simply removed

> [@dokken](#):
>
> ```auto
> nsp = PETSc.NullSpace().create(vectors=[null_vec])
> assert nsp.test(A)
> A.setNullSpace(nsp)
> 
> ```

in the demo

Actually I got a solution by running the code I modified but there are some differences between the result of the output and the expected result…  
Maybe I set **a** wrongly or maybe there are more modifications needed to be done in the defined function **nested\_iterative\_solver():**?

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [May 3, 2024, 5:07pm UTC](https://fenicsproject.discourse.group/t/question-about-how-to-define-a-function/14242/26 "2024-05-03T17:07:39Z")

</div>

Please provide a minimal reproducible example that uses only one of the four implementation given in

> [@guangyi\_cai](#):
>
> I just modified **demo\_stokes.py** downloaded from [Stokes equations using Taylor-Hood elements — DOLFINx 0.7.3 documentation](https://docs.fenicsproject.org/dolfinx/v0.7.3/python/demos/demo_stokes.html)

[Next page](https://fenicsproject.discourse.group/t/question-about-how-to-define-a-function/14242.md?page=2)