# The problem in a simple example I made and how to 'project' the value assigned to 'FunctionSpace' into a numpy array in dolfinx environment?

**URL:** <https://fenicsproject.discourse.group/t/the-problem-in-a-simple-example-i-made-and-how-to-project-the-value-assigned-to-functionspace-into-a-numpy-array-in-dolfinx-environment/9198>\
**Category:** dolfinx\
**Created:** [September 3, 2022, 6:22am UTC](https://fenicsproject.discourse.group/t/the-problem-in-a-simple-example-i-made-and-how-to-project-the-value-assigned-to-functionspace-into-a-numpy-array-in-dolfinx-environment/9198 "2022-09-03T06:22:38Z")\
**Posts on this page:** 12\
**Page:** 1

<div class="post-metadata">

**Author:** ![Swani](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/swani/32/3940_2.png) [@Swani](https://fenicsproject.discourse.group/u/Swani)\
**Post date:** [September 3, 2022, 6:22am UTC](https://fenicsproject.discourse.group/t/the-problem-in-a-simple-example-i-made-and-how-to-project-the-value-assigned-to-functionspace-into-a-numpy-array-in-dolfinx-environment/9198/1 "2022-09-03T06:22:38Z")

</div>

hi 😀  
Currently, I am using the dolfinx library.  
I was studying the dolfinx library by myself, and many questions arose, so I wrote a post on the forum.  
Condensing the questions into three,

1] Minor question?: The dolfinx library does not use docker, it is directly installed on ubuntu 22.04, and the version is assumed to be 0.4.1. Is there a separate way to check the dolfinx version?

2] To study the dolfinx library, I made an example of applying an external force to the simple 2D cantilever (spanning [0,L]\times[0,H]) as shown below.  
Compared with the reliable results obtained from the dolfin library, there is a difference and I think that there is an error when grafting the given formula condition to the code.

![화면 캡처 2022-08-30 140142](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/8/861fe838f0a0a69471c4496edf11d4b3cb9501fd.jpeg)  
(Afterwards, in the example code, L(`nelx`) was applied as 2, and H(`nely`) was also applied as 2)

The formula can be rearranged as:  
governing equation: -\nabla\cdot\sigma(\mathbf{u})=0  
Dirichlet boundary condition(clamp): \mathbf{u}=0 at \partial\Omega\_{D} (left of cantilever beam)

External load: \textbf{F}\_{T}=-1 on \partial\Omega\_{N}(for y-direction)

Then, using Cauchy stress tensor and strain energy equation such that  
Cauchy stress tensor: \sigma(\mathbf{u})=\lambda\textrm{tr}(\epsilon(\mathbf{u}))\textit{I}\ +2\mu\epsilon(\mathbf{u})  
Strain energy equation: \psi(\mathbf{u})=\frac{\lambda}{2}\{\textrm{tr}(\epsilon(\mathbf{u}))\}^{2}\ + \mu\textrm{tr}(\epsilon(\mathbf{u})^{2})  
(\epsilon(\mathbf{u})=\frac{1}{2}(\nabla\mathbf{u}\ + (\nabla\mathbf{u})^{T})=sym(\nabla\mathbf{u}), \lambda=0.6, \mu=0.4)

The code using the dolfinx library is written as follows

```auto
import ufl
import numpy as np
from dolfinx import fem, mesh
from mpi4py import MPI
from petsc4py import PETSc

```

> Declare the input variables

```auto
nelx = 2
nely = 2
volfrac = 0.5
penal = 3

```

> Setting preliminaries

```auto
sigma = lambda _u: 2.0 * mu * ufl.sym(ufl.grad(_u)) + lmd * ufl.tr(ufl.sym(ufl.grad(_u))) * ufl.Identity(len(_u))
psi = lambda _u: lmd / 2 * (ufl.tr(ufl.sym(ufl.grad(_u))) ** 2) + mu * ufl.tr(ufl.sym(ufl.grad(_u)) * ufl.sym(ufl.grad(_u)))

mu, lmd = PETSc.ScalarType(0.4), PETSc.ScalarType(0.6)

```

> Prepare finite element analysis

```auto
msh = mesh.create_rectangle(MPI.COMM_WORLD, ((0.0, 0.0), (nelx, nely)), (nelx, nely), cell_type=mesh.CellType.triangle)
U = fem.VectorFunctionSpace(msh, ("CG", 1))
D = fem.FunctionSpace(msh, ("DG", 0))
u, v = ufl.TrialFunction(U), ufl.TestFunction(U)
u_sol, density = fem.Function(U), fem.Function(D)
density.x.array[:] = volfrac

```

> Define support

```auto
def left_clamp(x):
    return np.isclose(x[0], 0.0)
f_dim = msh.topology.dim - 1
bc_facets = mesh.locate_entities_boundary(msh, f_dim, left_clamp)
u_zero = np.array([0.0, 0.0], dtype=PETSc.ScalarType)
bc_l = fem.dirichletbc(u_zero, fem.locate_dofs_topological(U, f_dim, bc_facets), U)
bcs = [bc_l]

```

> Define external load

```auto
load_points = [(1, lambda x: np.logical_and(x[0]==nelx, x[1]<=2))]

facet_indices, facet_markers = [], []
for (marker, locator) in load_points:
    facets = mesh.locate_entities(msh, f_dim, locator)
    facet_indices.append(facets)
    facet_markers.append(np.full(len(facets), marker))
facet_indices = np.array(np.hstack(facet_indices), dtype=np.int32)
facet_markers = np.array(np.hstack(facet_markers), dtype=np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = mesh.meshtags(msh, f_dim, facet_indices[sorted_facets], facet_markers[sorted_facets])

ds = ufl.Measure("ds", domain=msh, subdomain_data=facet_tag)
f = ufl.dot(v, fem.Constant(msh, PETSc.ScalarType((0.0, -1.0)))) * ds(1)

```

> Set up the variational problem and solver (using SIMP)

```auto
k = ufl.inner(density ** penal * sigma(u), ufl.grad(v)) * ufl.dx
problem = fem.petsc.LinearProblem(k, f, bcs)
u_sol = problem.solve()

```

Then, the result for `u_sol` array is

```auto
array([-17.05014904, -31.16299225, -23.40409419, -70.10238634,
        -1.23442586, -66.45633149, 0.17561615, -26.80259435,
         0. , 0. , 20.43334505, -65.04736033,
         0. , 0. , 15.06233617, -27.34430426,
         0. , 0. ])

```

the result I want to get (reliable results from code written with the dolfin library)

```auto
array([ 2.19098143e+01, -6.77188329e+01, 4.26325641e-15, -6.44827586e+01,
        1.51458886e+01, -2.69893899e+01, -1.51458886e+01, -2.69893899e+01,
        0.00000000e+00, 0.00000000e+00, 1.08212540e-15, -2.81564987e+01,
       -2.19098143e+01, -6.77188329e+01, 0.00000000e+00, 0.00000000e+00,
        0.00000000e+00, 0.00000000e+00])

```

I think that different results were obtained because of wrong coding when applying boundary condition or external load. Which part of the code could be considered bad?

3] Declare ‘FunctionSpace’ in the dolfinx environment, and arosed question how to convert the calculated value stored in the ‘FunctionSpace’ location to a numpy array.

Extending the above code example,  
Assuming that the u\_sol value is well obtained, how to map this value to a numpy array?  
(In the dolfin library, mapping was done as ‘project’)

```auto
project(psi(u_sol), D0).vector()[:]

```

(but I do not know how to do it in the dolfinx library.)

---

<div class="post-metadata">

**Author:** ![conpierce8](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/conpierce8/32/2464_2.png) [@conpierce8](https://fenicsproject.discourse.group/u/conpierce8)\
**Post date:** [September 5, 2022, 2:08am UTC](https://fenicsproject.discourse.group/t/the-problem-in-a-simple-example-i-made-and-how-to-project-the-value-assigned-to-functionspace-into-a-numpy-array-in-dolfinx-environment/9198/2 "2022-09-05T02:08:31Z")

</div>

1. Run the command

2. I obtain identical results with the code you provided and a code that I developed in old FEniCS (old FEniCS on the left, DOLFINx on the right):

3. I assume that you are asking how to project the strain energy density into a `DG-0` function. To do that, you can set up and solve the linear variational problem that represents a projection, i.e.:

---

<div class="post-metadata">

**Author:** ![Swani](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/swani/32/3940_2.png) [@Swani](https://fenicsproject.discourse.group/u/Swani)\
**Post date:** [September 5, 2022, 7:27am UTC](https://fenicsproject.discourse.group/t/the-problem-in-a-simple-example-i-made-and-how-to-project-the-value-assigned-to-functionspace-into-a-numpy-array-in-dolfinx-environment/9198/3 "2022-09-05T07:27:31Z")

</div>

Hi, @conpierce8  
Thank you very much for the detailed reply. 🙂

1] I checked the version with the command you gave me!  
Using command produces on my PC:

```auto
0.4.1
unknown

```

→ As an additional question, my PC is using version 0.4.1, how do I update the version(ex. update to 0.5.0)?

2] The reason why I wrote the ‘Define external load’ part in the dolfinx code I wrote like the code above

```auto
load_marker = MeshFunction("size_t", msh, msh.topology().dim() - 1)
CompiledSubDomain("x[0]==l && x[1]<=2", l=nelx).mark(load_marker, 1)
ds = Measure("ds")(subdomain_data=load_marker)

```

the code that in the dolfin library as follows to apply the load to the bottom-right corner point. Knowing that the dolfinx library should be applied similarly, I wrote it by quoting the above code using ‘mark’. However, I didn’t know that the code applied loads to the entire right side boundary.

Then, how should I write the code to load the bottom-right corner point in dolfinx?

---

<div class="post-metadata">

**Author:** ![Swani](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/swani/32/3940_2.png) [@Swani](https://fenicsproject.discourse.group/u/Swani)\
**Post date:** [September 5, 2022, 3:01pm UTC](https://fenicsproject.discourse.group/t/the-problem-in-a-simple-example-i-made-and-how-to-project-the-value-assigned-to-functionspace-into-a-numpy-array-in-dolfinx-environment/9198/4 "2022-09-05T15:01:33Z")

</div>

In my opinion, the difference in `u_sol` value between the dolfin library code and the code made with dolfinx library is the dolfin library additionally allocates a very small subdomain when the load is applied to the bottom-right corner.

I think this part needs to be corrected in dolfinx code, is it not a big deal?

---

<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:** [September 5, 2022, 5:36pm UTC](https://fenicsproject.discourse.group/t/the-problem-in-a-simple-example-i-made-and-how-to-project-the-value-assigned-to-functionspace-into-a-numpy-array-in-dolfinx-environment/9198/5 "2022-09-05T17:36:54Z")

</div>

> [@conpierce8](#):
>
> I assume that you are asking how to project the strain energy density into a `DG-0` function. To do that, you can set up and solve the linear variational problem that represents a projection, i.e.:
> 
> \int\_\Omega {\psi \phi , \mathrm{d}x} = \int\_\Omega {\psi(u\_h) \phi , \mathrm{d}x}
> 
> ∫Ωψϕdx=∫Ωψ(uh)ϕdx \int\_\Omega {\psi \phi , \mathrm{d}x} = \int\_\Omega {\psi(u\_h) \phi , \mathrm{d}x}
> 
> where \psi ψ\psi is the quantity you want to solve for, i.e. the `DG-0` projection of the strain energy density, \phi ϕ\phi is a test function, and \psi(u\_h) ψ(uh)\psi(u\_h) is the strain energy density calculated from the (known) discrete solution u\_h uhu\_h . This can be implemented as:
> 
> ```auto
> psi_trial = ufl.TrialFunction(D)
> psi_test = ufl.TestFunction(D)
> a = psi_trial * psi_test * ufl.dx
> L = psi(u_sol) * psi_test * ufl.dx
> projection_problem = fem.petsc.LinearProblem(a, L, [])
> psi_sol = projection_problem.solve()
> print(psi_sol.x.array)
> 
> ```

Note that a more scalable approach is to use interpolation, as explained in [Implementation — FEniCSx tutorial](https://jorgensd.github.io/dolfinx-tutorial/chapter2/linearelasticity_code.html#stress-computation)

> [@Swani](#):
>
> As an additional question, my PC is using version 0.4.1, how do I update the version(ex. update to 0.5.0)?

I guess you used `apt-get` to install DOLFINx, and then you would have to use `apt-get update` and `apt-get upgrade` to get the latest release

> [@Swani](#):
>
> The dolfinx library does not use docker, it is directly installed on ubuntu 22.04, and the version is assumed to be 0.4.1.

I’m not quite sure what is the status of releases to apt, @dparsons can probably give a concrete answer to this.

---

<div class="post-metadata">

**Author:** ![dparsons](https://avatars.discourse-cdn.com/v4/letter/d/5e9695/32.png) [@dparsons](https://fenicsproject.discourse.group/u/dparsons)\
**Post date:** [September 5, 2022, 6:01pm UTC](https://fenicsproject.discourse.group/t/the-problem-in-a-simple-example-i-made-and-how-to-project-the-value-assigned-to-functionspace-into-a-numpy-array-in-dolfinx-environment/9198/6 "2022-09-05T18:01:04Z")

</div>

Debian has had some problems with the pmix (MPI) package. That’s now resolved so dolfinx 0.5 packages will be built soon.

---

<div class="post-metadata">

**Author:** ![Swani](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/swani/32/3940_2.png) [@Swani](https://fenicsproject.discourse.group/u/Swani)\
**Post date:** [September 5, 2022, 10:20pm UTC](https://fenicsproject.discourse.group/t/the-problem-in-a-simple-example-i-made-and-how-to-project-the-value-assigned-to-functionspace-into-a-numpy-array-in-dolfinx-environment/9198/7 "2022-09-05T22:20:53Z")

</div>

1] The answer to the question is now complete. Thank you! 😄

---

<div class="post-metadata">

**Author:** ![Swani](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/swani/32/3940_2.png) [@Swani](https://fenicsproject.discourse.group/u/Swani)\
**Post date:** [September 5, 2022, 10:23pm UTC](https://fenicsproject.discourse.group/t/the-problem-in-a-simple-example-i-made-and-how-to-project-the-value-assigned-to-functionspace-into-a-numpy-array-in-dolfinx-environment/9198/8 "2022-09-05T22:23:32Z")

</div>

Looking forward to the next update! Thank you! 🙂

---

<div class="post-metadata">

**Author:** ![conpierce8](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/conpierce8/32/2464_2.png) [@conpierce8](https://fenicsproject.discourse.group/u/conpierce8)\
**Post date:** [September 6, 2022, 1:47pm UTC](https://fenicsproject.discourse.group/t/the-problem-in-a-simple-example-i-made-and-how-to-project-the-value-assigned-to-functionspace-into-a-numpy-array-in-dolfinx-environment/9198/9 "2022-09-06T13:47:50Z")

</div>

> [@Swani](#):
>
> the dolfin library additionally allocates a very small subdomain when the load is applied to the bottom-right corner.

I would not expect this behavior when using a built-in mesh. Can you share your DOLFIN code?

---

<div class="post-metadata">

**Author:** ![Swani](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/swani/32/3940_2.png) [@Swani](https://fenicsproject.discourse.group/u/Swani)\
**Post date:** [September 7, 2022, 1:57am UTC](https://fenicsproject.discourse.group/t/the-problem-in-a-simple-example-i-made-and-how-to-project-the-value-assigned-to-functionspace-into-a-numpy-array-in-dolfinx-environment/9198/10 "2022-09-07T01:57:50Z")

</div>

> [@conpierce8](#):
>
> I would not expect this behavior when using a built-in mesh. Can you share your DOLFIN code?

I used a simple rectangular mesh in dolfin. Here’s the code:

```auto
mesh = RectangleMesh(Point(0, 0), Point(nelx, nely), nelx, nely, "right/left")

```

The constraint was applied only to the left part of the domain as follows

```auto
bc = CompiledSubDomain("near(x[0], 0.0, tol) && on_boundary", tol=1e-14)
bcs = [DirichletBC(U, Constant((0.0, 0.0)), bc)]

```

And as in the above reply, load condition was given,

```auto
load_marker = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
CompiledSubDomain("x[0]==l && x[1]<=2", l=nelx).mark(load_marker, 1)
ds = Measure("ds")(subdomain_data=load_marker)
f = dot(v, Constant((0.0, -1.0))) * ds(1)

```

and `u_sol` was obtained with

```auto
k = inner(density**penal*sigma(u), grad(v)) * dx
problem = LinearVariationalProblem(k, f, u_sol, bcs=bcs)
solver = LinearVariationalSolver(problem)

solver.solve()

```

I composed the dolfin code like this, and I thought that the same value (even if the dof indexing is different) should appear in the dolfinx code, but it did not come out.

I don’t know why the values are different. 😥

---

<div class="post-metadata">

**Author:** ![conpierce8](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/conpierce8/32/2464_2.png) [@conpierce8](https://fenicsproject.discourse.group/u/conpierce8)\
**Post date:** [September 7, 2022, 1:52pm UTC](https://fenicsproject.discourse.group/t/the-problem-in-a-simple-example-i-made-and-how-to-project-the-value-assigned-to-functionspace-into-a-numpy-array-in-dolfinx-environment/9198/11 "2022-09-07T13:52:21Z")

</div>

> [@Swani](#):
>
> `mesh = RectangleMesh(Point(0, 0), Point(nelx, nely), nelx, nely, "right/left")`

You are getting different results because your old DOLFIN code and DOLFINx code are using different meshes (`"right/left"` for old DOLFIN and `dolfinx.mesh.DiagonalType.right` for DOLFINx). See screenshot below

 ![screenshot2](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/e/e53bbe52d0333ca7cad7b3c651840004ea132173.png)

If you change your DOLFINx mesh generation code to the following, you should see that the solutions are identical

```python
msh = mesh.create_rectangle(
    MPI.COMM_WORLD,
    ((0.0, 0.0), (nelx, nely)),
    (nelx, nely),
    cell_type=mesh.CellType.triangle,
    diagonal=mesh.DiagonalType.right_left,
)

```

As a final note, your code in old DOLFIN is also applying the traction to the entire right boundary, not a point load at the bottom right corner. The following statement marks all facets on the right boundary because `x[1]<=2` is true for all points on the right boundary.

> [@Swani](#):
>
> `CompiledSubDomain("x[0]==l && x[1]<=2", l=nelx).mark(load_marker, 1)`

If you want to apply a point load at the bottom right corner, you should use `PointSource` in old DOLFIN (see [here](https://bitbucket.org/fenics-project/dolfin/src/74a6ac2095059ea6582fe4ea11686c1c723e14f5/python/dolfin/ __init__.py#lines-57) and [here](https://bitbucket.org/fenics-project/dolfin/src/master/dolfin/fem/PointSource.h)). Unfortunately, the reimplementation in DOLFINx is not complete yet; see [PointSource in dolfinx](https://fenicsproject.discourse.group/t/pointsource-in-dolfinx/8337).

---

<div class="post-metadata">

**Author:** ![Swani](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/swani/32/3940_2.png) [@Swani](https://fenicsproject.discourse.group/u/Swani)\
**Post date:** [September 8, 2022, 1:56am UTC](https://fenicsproject.discourse.group/t/the-problem-in-a-simple-example-i-made-and-how-to-project-the-value-assigned-to-functionspace-into-a-numpy-array-in-dolfinx-environment/9198/12 "2022-09-08T01:56:54Z")

</div>

The problem seems to have been resolved! I learned a lot about dolfin and dolfinx thanks to @conpierce8.  
Thank you so much!! 😄 😆

It was confirmed that almost same value was derived when the mesh was given in dolfinx. Awesome!  
I thought that dolfin was also giving an external load to the Point, but I found out that it wasn’t. (Sorry for the confusion.)