Applying Multipoint Constraint on subspace using dolfinx_mpc

General
1

Dear. All
I would like to use dolfinx_mpc extension to apply multipoint constraint on the right edge of 2D rectangular domain. The purpose is to ensure that the right edge remain straight throughout simulation. I implies that the displacement in the x direction of all nodes on this boundary is the same. To achieve this goal, I first specified the degree of freedom associated to the right edge and then selected one degree of freedom as the master and the others as the salves DOFs. After running the code I received the following error:

Traceback (most recent call last):
File “/home/mohhammad/Documents/files/simulations/090/fenics-reaction-diffusion/reaction-diffusion-test2D2-scaled-fenicsx-stress-mpc.py”, line 172, in
MPC_right.add_constraint(ME.sub(2).sub(0), np.array([slave], dtype= np.int32), master_dof, np.array([1.0]),0,0)
~~~~~~~~~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File “/home/mohhammad/anaconda3/envs/fenicsx/lib/python3.13/site-packages/dolfinx_mpc/multipointconstraint.py”, line 145, in add_constraint
assert V == self.V
^^^^^^^^^^^
AssertionError
any assistance will be appreciated. I enclosed the relevant part of code:


from mpi4py import MPI

import dolfinx

from dolfinx.mesh import create_rectangle

import ufl

import basix.ufl

from dolfinx.mesh import meshtags, locate_entities_boundary

from dolfinx.fem import functionspace, Function, locate_dofs_topological, Constant, locate_dofs_geometrical

from dolfinx.mesh import compute_midpoints

from dolfinx import geometry

import numpy as np

from dolfinx.io import XDMFFile

from dolfinx.cpp.mesh import CellType

from petsc4py import PETSc

from ufl.classes import *

from ufl.algorithms import *

from ufl import Measure

from dolfinx.fem.petsc import NonlinearProblem

from dolfinx.nls.petsc import NewtonSolver

import os

from pathlib import Path

from mpi4py import MPI

from dolfinx import default_real_type, log, plot

from petsc4py.PETSc import ScalarType

import dolfinx_mpc

# -----------------------------

# Time parameters

# -----------------------------

#num_steps = 480000

num_steps= 3600

t = 0.0

t_end = 10*3600

dt = t_end/num_steps

BC_thickness = 100e-6

Oxide_thickness = 1e-6

Lx = 60e-6

Ly = BC_thickness + Oxide_thickness

x_bar = 10e-6

BC_bar = BC_thickness/x_bar

Oxide_bar = Oxide_thickness/x_bar

Lx_bar = Lx/x_bar

Ly_bar = Ly/x_bar

refine_length_bar = 0.01*BC_bar

print(refine_length_bar)

c_characteristic = 0.32e6

#nx, ny = 200,335

nx,ny = 200,335

refined_domain = create_rectangle(MPI.COMM_WORLD, [np.array([0,0]), np.array([Lx_bar,Ly_bar])], [nx,ny])

tdim = refined_domain.topology.dim

def refine_cells(x, tol=1e-16):

return x[1] >= refine_length_bar + tol

with XDMFFile(MPI.COMM_WORLD,"refined_domain_scaled.xdmf","w") as xdmf:

xdmf.write_mesh(refined_domain)

element_degree = 1

variant = basix.LagrangeVariant.equispaced

c_el = basix.ufl.element("Lagrange", refined_domain.topology.cell_name(), element_degree, variant, dtype = default_real_type)

n_el = basix.ufl.element("Lagrange", refined_domain.topology.cell_name(), element_degree, variant, dtype = default_real_type)

u_el = basix.ufl.element("Lagrange", refined_domain.topology.cell_name(), element_degree, variant, dtype = default_real_type, shape = (2,))

s_el = basix.ufl.element("Lagrange", refined_domain.topology.cell_name(), element_degree, variant, dtype = default_real_type, shape = (2,2))

total_elements = basix.ufl.mixed_element([c_el,n_el,u_el])

ME = dolfinx.fem.functionspace(refined_domain, total_elements)

TT = dolfinx.fem.functionspace(refined_domain, s_el)

sig = dolfinx.fem.Function(TT)

w = dolfinx.fem.Function(ME)

w_old = dolfinx.fem.Function(ME)

dw = ufl.TrialFunction(ME)

x= ufl.SpatialCoordinate(refined_domain)

w.sub(1).interpolate(lambda x: np.where(x[1] >= BC_bar, 1.0, 0.0))

with XDMFFile(refined_domain.comm, "phase_sub.xdmf", "w") as xdmf:

xdmf.write_mesh(refined_domain)

xdmf.write_function(w.sub(1))

def top(x):

return np.isclose(x[1], Ly_bar)

def bottom(x):

return np.isclose(x[1], 0.0)

def left(x):

return np.isclose(x[0], 0.0)

def right(x):

return np.isclose(x[0], Lx_bar)

u_space,_ = ME.sub(2).collapse()

u_left = dolfinx.fem.Function(u_space)

u_left.x.array[:] = 0.0

left_facets = dolfinx.mesh.locate_entities_boundary(refined_domain, refined_domain.topology.dim-1,left)

left_u_dofs = dolfinx.fem.locate_dofs_topological((ME.sub(2),u_space),tdim-1,left_facets)

u_left_bc = dolfinx.fem.dirichletbc(u_left,left_u_dofs, ME.sub(2).sub(0))

u_bot = dolfinx.fem.Function(u_space)

u_bot.x.array[:] = 0.0

bot_u_dofs = dolfinx.fem.locate_dofs_topological((ME.sub(2),u_space),tdim-1,bottom_facets)

u_bot_bc = dolfinx.fem.dirichletbc(u_bot,bot_u_dofs, ME.sub(2).sub(1))

right_facets = dolfinx.mesh.locate_entities_boundary(refined_domain, refined_domain.topology.dim-1,right)

right_u_dofs = dolfinx.fem.locate_dofs_topological((ME.sub(2),u_space),tdim-1,right_facets)

ux_dofs = right_u_dofs[0]

master_dof = ux_dofs[0]

MPC_right = dolfinx_mpc.MultiPointConstraint(ME.sub(2).sub(0))

slaves = right_u_dofs[0][1:]

MPC_right.add_constraint(ME.sub(2).sub(0), np.array([slaves], dtype= np.int32), master_dof, np.array([1.0]),0,0)
2

have you had a look at Rigidity constraint · Issue #9 · jorgensd/dolfinx_mpc · GitHub
which I believe implements the constraint you are thinking of

1 Like
3

Massive thanks for your guidance. I applied the procedure in rigidity constrain. Issue #9. Afterwards, the code works but from the results it seems that the multipoint constraint did not function properly because the x-displacement of nodes especially in upper right part of the right edge are not equal. I enclose the updated code and the plot of x displacement of the right edge for further clarification.

ux-animation.0000
ux-animation.00002030×1299 45.9 KB


from mpi4py import MPI

import dolfinx

from dolfinx.mesh import create_rectangle

import ufl

import basix.ufl

from dolfinx.mesh import meshtags, locate_entities_boundary

from dolfinx.fem import functionspace, Function, locate_dofs_topological, Constant, locate_dofs_geometrical

from dolfinx.mesh import compute_midpoints

from dolfinx import geometry

import numpy as np

from dolfinx.io import XDMFFile

from dolfinx.cpp.mesh import CellType

from petsc4py import PETSc

from ufl.classes import *

from ufl.algorithms import *

from ufl import Measure

from dolfinx.fem.petsc import NonlinearProblem

from dolfinx.nls.petsc import NewtonSolver

import os

from pathlib import Path

from mpi4py import MPI

from dolfinx import default_real_type, log, plot

from petsc4py.PETSc import ScalarType

import dolfinx_mpc

# -----------------------------

# Time parameters

# -----------------------------

#num_steps = 480000

num_steps= 3600

t = 0.0

t_end = 10*3600

dt = t_end/num_steps

BC_thickness = 100e-6

Oxide_thickness = 1e-6

Lx = 60e-6

Ly = BC_thickness + Oxide_thickness

x_bar = 10e-6

BC_bar = BC_thickness/x_bar

Oxide_bar = Oxide_thickness/x_bar

Lx_bar = Lx/x_bar

Ly_bar = Ly/x_bar

refine_length_bar = 0.01*BC_bar

print(refine_length_bar)

c_characteristic = 0.32e6

nx, ny = 200,335

refined_domain = create_rectangle(MPI.COMM_WORLD, [np.array([0,0]), np.array([Lx_bar,Ly_bar])], [nx,ny])

tdim = refined_domain.topology.dim

def refine_cells(x, tol=1e-16):

return x[1] >= refine_length_bar + tol

element_degree = 1

variant = basix.LagrangeVariant.equispaced

c_el = basix.ufl.element("Lagrange", refined_domain.topology.cell_name(), element_degree, variant, dtype = default_real_type)

n_el = basix.ufl.element("Lagrange", refined_domain.topology.cell_name(), element_degree, variant, dtype = default_real_type)

u_el = basix.ufl.element("Lagrange", refined_domain.topology.cell_name(), element_degree, variant, dtype = default_real_type, shape = (2,))

s_el = basix.ufl.element("Lagrange", refined_domain.topology.cell_name(), element_degree, variant, dtype = default_real_type, shape = (2,2))

total_elements = basix.ufl.mixed_element([c_el,n_el,u_el])

ME = dolfinx.fem.functionspace(refined_domain, total_elements)

TT = dolfinx.fem.functionspace(refined_domain, s_el)

sig = dolfinx.fem.Function(TT)

w = dolfinx.fem.Function(ME)

w_old = dolfinx.fem.Function(ME)

dw = ufl.TrialFunction(ME)

c_new,n_new,u_new = ufl.split(w)

c_old,n_old,u_old = ufl.split(w_old)

v_c,v_n,v_u = ufl.TestFunctions(ME)

dc,dn,du = ufl.TrialFunctions(ME)

w.sub(0).x.array[:] = 0.0

w.sub(1).x.array[:] = 0.0

w.sub(2).x.array[:] = 0.0

x= ufl.SpatialCoordinate(refined_domain)

w.sub(1).interpolate(lambda x: np.where(x[1] >= BC_bar, 1.0, 0.0))

def top(x):

return np.isclose(x[1], Ly_bar)

def bottom(x):

return np.isclose(x[1], 0.0)

def left(x):

return np.isclose(x[0], 0.0)

def right(x):

return np.isclose(x[0], Lx_bar)

c_space, _ = ME.sub(0).collapse()

c_bc_top = dolfinx.fem.Function(c_space)

c_bc_top.x.array[:] = 5e-5

top_facets = dolfinx.mesh.locate_entities_boundary(refined_domain, tdim-1, top)

top_dofs = dolfinx.fem.locate_dofs_topological((ME.sub(0), c_space), tdim-1, top_facets)

c_space_bot, _ = ME.sub(0).collapse()

c_bc_bot = dolfinx.fem.Function(c_space)

c_bc_bot.x.array[:] = 0.0

bottom_facets = dolfinx.mesh.locate_entities_boundary(refined_domain, tdim-1, bottom)

bot_dofs = dolfinx.fem.locate_dofs_topological((ME.sub(0), c_space_bot), tdim-1, bottom_facets)

bc_c_top = dolfinx.fem.dirichletbc(c_bc_top, top_dofs,ME.sub(0))

bc_c_bot = dolfinx.fem.dirichletbc(c_bc_bot,bot_dofs,ME.sub(0))

u_space,_ = ME.sub(2).collapse()

u_left = dolfinx.fem.Function(u_space)

u_left.x.array[:] = 0.0

left_facets = dolfinx.mesh.locate_entities_boundary(refined_domain, refined_domain.topology.dim-1,left)

left_u_dofs = dolfinx.fem.locate_dofs_topological((ME.sub(2),u_space),tdim-1,left_facets)

u_left_bc = dolfinx.fem.dirichletbc(u_left,left_u_dofs, ME.sub(2).sub(0))

u_bot = dolfinx.fem.Function(u_space)

u_bot.x.array[:] = 0.0

bot_u_dofs = dolfinx.fem.locate_dofs_topological((ME.sub(2),u_space),tdim-1,bottom_facets)

u_bot_bc = dolfinx.fem.dirichletbc(u_bot,bot_u_dofs, ME.sub(2).sub(1))

bcs_problem = [bc_c_top, bc_c_bot,u_left_bc,u_bot_bc]

def mapping(x):

vals = np.zeros(x.shape)

vals[0] = x[0]

vals[1] = 0

vals[2] = 0

return vals

mpc = dolfinx_mpc.MultiPointConstraint(ME.sub(2))

mpc.create_periodic_constraint_geometrical(ME.sub(2).sub(0), right, mapping,bcs_problem)

mpc.finalize()
4

Could you supply your full code, with 3x` encapsulation, as currently the indentation is wrong, and the code ends at finalizing the MPC, not at the end of solve.
To me, it also seems like you are importing things that are not needed, such as

The following is unused:

and I would be careful with:

as those should be from your MPC space, as discussed in:

5

In accordance with your note, I changed the function space of solution function to the function space of mpc and now the code works and it seems that rigidity constraint on the right edge is functioned properly. From this communication, I have recognized that mpc cannot be instantiated in subspaces (MultiPointConstraint(subspace)) but the methods like create_periodic_constraint_geometrical() and create_periodic_constraint_topological() can be defined on space or appropriate subspaces. please verify this intuition.
best regards.

6

The MultiPointConstraint class should be initiated on the patent space, as it is modifying the global (parent) matrix. As you state the constructors can work o subspaces.

7

I am so grateful.
best regards.