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

[Swani]

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.

(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

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

Declare the input variables

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

Setting preliminaries

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

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

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

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)

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

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)

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’)

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

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

[conpierce8]

1. Run the command

   python3 -c "import dolfinx as dfx;print(dfx.__version__);print(dfx.git_commit_hash)"
   

   which produces on my machine:

   0.5.0
   2aaf3b20dbaedcbd3925a9640c3859deec563e02
   

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):

   

   Keep in mind that the dof ordering in old FEniCS and DOLFINx may be different. Also, in your DOLFINx code you apply a traction to the entire right boundary (load_points = [(1, lambda x: np.logical_and(x[0]==nelx, x[1]<=2))]), while your schematic appears to show a point load on the bottom-right corner. Is your code correctly implementing the intended boundary condition?

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.:

   \int_\Omega {\psi \phi \, \mathrm{d}x} = \int_\Omega {\psi(u_h) \phi \, \mathrm{d}x}

   where \psi is the quantity you want to solve for, i.e. the DG-0 projection of the strain energy density, \phi is a test function, and \psi(u_h) is the strain energy density calculated from the (known) discrete solution u_h. This can be implemented as:

   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)
   

   producing

   [ 79.91107539 111.61468049 211.04645563  66.28836385 143.69740152
    122.77015115  28.56753585 308.35397476]
   

[Swani]

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:

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

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?

[Swani]

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?

[dokken]

@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:

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

@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.

[dparsons]

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

[Swani]

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

[Swani]

Looking forward to the next update! Thank you!

[conpierce8]

@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?

[Swani]

@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:

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

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,

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

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.

[conpierce8]

@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

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

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 and here). Unfortunately, the reimplementation in DOLFINx is not complete yet; see PointSource in dolfinx.

[Swani]

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.)