Time-dependent dirichlet condition issue : AttributeError

dolfinx
1

Dear community,

I recently switched from dolfin to dolfinx and I am trying to udpate the following “Expression” class.

V         = VectorFunctionSpace(mesh,'Lagrange', 1)
U_app      = Expression( ('t', '0.'), t = 0.0, degree = 1 ) # Constant((0.,0.))
bcu        = [DirichletBC( V , Constant((0.,0.)), "near(x[1],0) && on_boundary") ,
              DirichletBC( V , U_app, "near(x[1],1.) && on_boundary" )]

In the following script, I tried to follow what it already did in the discussion https://fenicsproject.discourse.group/t/equivalent-for-expression-in-dolfinx/2641/12, but I met trouble at the interpolation :

import ufl
import ffcx
import matplotlib.pyplot as plt
import numpy as np
from dolfinx import log, io, mesh, fem, cpp
from mpi4py import MPI
from petsc4py.PETSc import ScalarType

log.set_log_level(log.LogLevel.OFF)
comm  = MPI.COMM_WORLD

msh = mesh.create_unit_square(MPI.COMM_WORLD, 151, 151, mesh.CellType.quadrilateral)
ndim = msh.topology.dim

# Time discretization

T = 0.02                            # final time
num_steps = 200                     # number of time steps
dt = T / num_steps                  # time step size
time_step = np.linspace(0, T, num_steps)

# BC 

V_u         = fem.VectorFunctionSpace(msh,('CG', 1))
V_alpha   = fem.FunctionSpace(msh,('CG', 1))

class MyExpression:
    def __init__(self):
        self.t = 0.0
    def evaluate(self, x):
        return np.full(x.shape[1], self.t)

u_app = MyExpression()
u_app.t = 0.0
u_app_func = fem.Function(V_u)
# len(u_app_func.x.array) = 46208 = 2* len(msh.geometry.x)
#u_app_func = fem.Function(V_alpha)
u_app_func.interpolate(u_app.evaluate)


boundaries = [(1, lambda x: np.isclose(x[0], 0)),
              (2, lambda x: np.isclose(x[1], 1))]

facet_indices, facet_markers = [], []
bcu = []
fdim = ndim - 1
for (marker, locator) in boundaries:
    facets = mesh.locate_entities(msh, fdim, locator)
    if marker == 1:
        facets_dofs = fem.locate_dofs_geometrical(V_u,locator)
        u_zero = np.array((0. ,)*msh.geometry.dim, dtype=ScalarType) # print(u_zeros) = [0.0, 0.]
        bc = fem.dirichletbc( u_zero,facets_dofs, V_u)
    if marker == 2:
        facets_dofs = fem.locate_dofs_topological(V_u.sub(0), fdim, marker)
        bc = fem.dirichletbc( u_app_func, facets_dofs,V_u.sub(0))
    bcu.append(bc)
    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, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])

I get an error about the dimension of the array :

Traceback (most recent call last):
  File "/stck/lmersel/CODE_PROGRAMME/anaconda3/envs/fenics/lib/python3.10/site-packages/dolfinx/fem/function.py", line 318, in _interpolate
    self._cpp_object.interpolate(u._cpp_object, cells)
AttributeError: 'function' object has no attribute '_cpp_object'

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/scratchm/lmersel/FEniCSx/mesh_quad_newton_stagg_AT2_dlfx/test.py", line 47, in <module>
    u_app_func.interpolate(u_app.evaluate)
  File "/stck/lmersel/CODE_PROGRAMME/anaconda3/envs/fenics/lib/python3.10/site-packages/dolfinx/fem/function.py", line 336, in interpolate
    _interpolate(u, cells)
  File "/stck/lmersel/CODE_PROGRAMME/anaconda3/envs/fenics/lib/python3.10/functools.py", line 889, in wrapper
    return dispatch(args[0].__class__)(*args, **kw)
  File "/stck/lmersel/CODE_PROGRAMME/anaconda3/envs/fenics/lib/python3.10/site-packages/dolfinx/fem/function.py", line 320, in _interpolate
    self._cpp_object.interpolate(u, cells)
RuntimeError: Expected 2D array of data

I checked the length of u_app_func = fem.Function(V_u) and I noticed that it is twice the number of node even if I try u_app_func.sub(0) or u_app_func.sub(1)

In summary, I just want to apply ux(t)=t and uy = 0 on the top edge of the square.

Could you help me to understand the mistake ?

Thanks,

Best regards,

Lamia

2

I also tried an other way to apply a dirichlet condition evolving in time with the code :

msh = mesh.create_unit_square(MPI.COMM_WORLD, 151, 151, mesh.CellType.quadrilateral)
ndim = msh.topology.dim

# Time discretization

T         = 0.02                            # final time
num_steps = 200                     # number of time steps
time_disc = np.linspace(0, T, num_steps)
load_ref  = 1.
loads     = load_ref*np.linspace(0, T, num_steps)


V_u       = fem.VectorFunctionSpace(msh,('CG', 1))
V_alpha   = fem.FunctionSpace(msh,('CG', 1))

u_app_func = fem.Function(V_u.sub(0).collapse())

with u_app_func.vector.localForm() as bc_local:
    bc_local.set(0.0)

u_zero = np.array((0. ,)*msh.geometry.dim, dtype=ScalarType)
uy_zero = ScalarType(0.)

facets_dofs_1 = fem.locate_dofs_geometrical(V_u,lambda x: np.isclose(x[0], 0))
facets_dofs_2 = fem.locate_dofs_topological(V_u.sub(0), ndim - 1, lambda x: np.isclose(x[1], 1))
facets_dofs_3 = fem.locate_dofs_topological(V_u.sub(1), ndim - 1, lambda x: np.isclose(x[1], 1))
bc_1 = fem.dirichletbc( u_zero,facets_dofs_1, V_u)
bc_2 = fem.dirichletbc( u_app_func, facets_dofs_2,V_u.sub(0))
bc_3 = fem.dirichletbc( uy_zero, facets_dofs_3,V_u.sub(1))
bc = [bc_1, bc_2, bc_3]


for i_t, t in enumerate(loads):
    with u_app_func.vector.localForm() as bc_local:
        bc_local.set(t)

But I am stuck to the AttributeError :

Traceback (most recent call last):
  File "/scratchm/lmersel/FEniCSx/mesh_quad_newton_stagg_AT2_dlfx/test.py", line 43, in <module>
    u_app_func = fem.Function(V_u.sub(0).collapse())
  File "/stck/lmersel/CODE_PROGRAMME/anaconda3/envs/fenics/lib/python3.10/site-packages/dolfinx/fem/function.py", line 245, in __init__
    self._cpp_object = functiontype(dtype)(V._cpp_object)
AttributeError: 'tuple' object has no attribute '_cpp_object'

Thanks for the help,

Lamia

3

Consider the following code:

import numpy as np
from dolfinx import log, mesh, fem
from mpi4py import MPI
from petsc4py.PETSc import ScalarType

log.set_log_level(log.LogLevel.OFF)
comm = MPI.COMM_WORLD

msh = mesh.create_unit_square(comm, 2, 2, mesh.CellType.quadrilateral)
ndim = msh.topology.dim

V_u = fem.VectorFunctionSpace(msh, ('CG', 1))


t = 0.5
u_app = fem.Constant(msh, ScalarType((t, 0)))


boundaries = [(1, lambda x: np.isclose(x[0], 0)),
              (2, lambda x: np.isclose(x[1], 1))]

facet_indices, facet_markers = [], []
bcu = []
fdim = ndim - 1
for (marker, locator) in boundaries:
    facets = mesh.locate_entities(msh, fdim, locator)
    if marker == 1:
        facets_dofs = fem.locate_dofs_geometrical(V_u, locator)
        u_zero = np.array((0.,)*msh.geometry.dim, dtype=ScalarType)
        bc = fem.dirichletbc(u_zero, facets_dofs, V_u)
    if marker == 2:
        facets_dofs = fem.locate_dofs_topological(V_u, fdim, facets)
        bc = fem.dirichletbc(u_app, facets_dofs, V_u)
    bcu.append(bc)
    facet_indices.append(facets)
    facet_markers.append(np.full(len(facets), marker))

u_test = fem.Function(V_u)
u_test.x.array[:] = -1
fem.petsc.set_bc(u_test.vector, bcu)
print(u_test.x.array)
# Change value of t in the constant
u_app.value[0] = 0.1
fem.petsc.set_bc(u_test.vector, bcu)
print(u_test.x.array)

Topics such as these has also been covered in Test problem 2: Flow past a cylinder (DFG 2D-3 benchmark) — FEniCSx tutorial and various other problems in the tutorial.

4

Thanh you very much, this is exactly what I needed. But I wonder how to extract the x and y displacement. I thought that it was with sub(0) and sub(1). But u_test.x.array show exactly the same vector than u_test.sub(0).x.array

Besides at the second print :

# Change value of t in the constant
u_app.value[0] = 0.1
fem.petsc.set_bc(u_test.vector, bcu)
print(u_test.x.array)

I get

[ 0.   0.   0.   0.  -1.  -1.  -1.  -1.   0.1  0.   0.1  0.  -1.  -1.
 -1.  -1.   0.1  0. ]

Indeed, I have 3 nodes at the top where ux = 0.1 but I should have 9 nodes where uy = 0 and 6 nodes with ux = -1. So, I think I don’t get how the u_test vector is ordered.
If I indexe like u_test[0] and u_test[1], do I need to project to have the conponent?

Thank you for the explications,

Lamia

5

See: How to get the solution obtained by fenics as a vector? - #7 by dokken

6

Thank you for the reference. Indeed, this is what I have done with fenics 2019 such as

vector = np.zeros((num_sub_spaces, num_dofs_per_component))
for i in range(num_sub_spaces):
    vector[i] = u.sub(i, deepcopy=True).vector().get_local()

But with dolfinx I didn’t manage to extract with sub() function.

msh = mesh.create_unit_square(comm, 2, 2, mesh.CellType.quadrilateral)
ndim = msh.topology.dim

V_u = fem.VectorFunctionSpace(msh, ('CG', 1))
t = 0.5
u_app = fem.Constant(msh, ScalarType((t, 0)))
boundaries = [(1, lambda x: np.isclose(x[0], 0)),
              (2, lambda x: np.isclose(x[1], 1))]

facet_indices, facet_markers = [], []
bcu = []
fdim = ndim - 1
for (marker, locator) in boundaries:
    facets = mesh.locate_entities(msh, fdim, locator)
    if marker == 1:
        facets_dofs = fem.locate_dofs_geometrical(V_u, locator)
        u_zero = np.array((0.,)*msh.geometry.dim, dtype=ScalarType)
        bc = fem.dirichletbc(u_zero, facets_dofs, V_u)
    if marker == 2:
        facets_dofs = fem.locate_dofs_topological(V_u, fdim, facets)
        bc = fem.dirichletbc(u_app, facets_dofs, V_u)
    bcu.append(bc)
    facet_indices.append(facets)
    facet_markers.append(np.full(len(facets), marker))

u_test = fem.Function(V_u)
vector = np.zeros((V_u.dofmap.index_map_bs, V_u.dofmap.index_map.size_local)) # shape : (2,9)
for i in range(V_u.dofmap.index_map_bs):
      vector[i] = u_test.sub(i).vector.array

ValueError: could not broadcast input array from shape (18,) into shape (9,)

I am going to do the tutorial on th Deflection of the membrane where you plot the deflection (uy) thank to uh.eval().
But could you just specify what “BoundingBoxTree(domain, domain.topology.dim)” does when I apply the function to the mesh ‘msh’. It returns this :

In [37]: msh = mesh.create_unit_square(comm, 2, 2, mesh.CellType.quadrilateral)
In [38]: bbox_0 = geometry.BoundingBoxTree(msh, msh.topology.dim)

In [39]: print(bbox_0)
[0 0 0 ]->[1 1 0 ]
{[0 0 0 ]->[0.5 1 0 ]
{[0 0 0 ]->[0.5 0.5 0 ]
leaf containing entity (0),
[0 0.5 0 ]->[0.5 1 0 ]
leaf containing entity (1)}
,
[0.5 0 0 ]->[1 1 0 ]
{[0.5 0 0 ]->[1 0.5 0 ]
leaf containing entity (2),
[0.5 0.5 0 ]->[1 1 0 ]
leaf containing entity (3)}
}

I am going to do the tutorial properly, it would be better !

Thank you for the help.

Lamia

7

Consider:

assert(V_u.dofmap.index_map_bs == V_u.num_sub_spaces)
vector = np.zeros((V_u.dofmap.index_map_bs, V_u.dofmap.index_map.size_local))

for i in range(V_u.dofmap.index_map_bs):
    vector[i] = u_test.x.array[i::V_u.dofmap.index_map_bs]

coords = V_u.tabulate_dof_coordinates()
for x, val in zip(coords, vector.T):
    print(x, val)

The bounding box tree is a binary tree structure to quickly look up possible collisions between cells and other sets of points (as a single point or a bounding box tree).
Each cell can be bounding by a 3D box, which can be described by 6 floating values, i.e. bbox = (x_min, y_min, z_min), (x_max, y_max, z_max).
This is used to limit how often one uses an machine precision accurate collision algorithm to compute actual collisions.

8

Perfect, thank for the explanations. I understood that one out of every two components belongs to the same direction x or y. I wasn’t sure about this order.

Concerning the bounding box tree, I better understood the structure the function returns previous,

Thanks,

Lamia