Thank you so much for this example, appreciate your work and code!
Hi Dokken,
Sorry to bother you, I just tested your sample code, and it works very well. May I ask you to teach me about the “subspace_slave” and “subspace_master”?
Say if I have a bifurcation mesh, there are three points at the junction location, I want the outlet of mother branch=First inlet of daughter branch+ Second inlet of daughter branch
How should I define the “subspace_slave” and “subspace_master” here?
Many thanks!
Subspace slave and subspace master is if you use a vector function space for your primary variables. then, if you want to constrain only a component of the solution at a point, you use these variables.
In your problem, it seems to me that you want to use something like:
s_m_c = {l2b([mother_x, mother_y]): {l2b([inlet1_x, inlet1_y]):1, l2b([inlet2_x, inlet2_y]): 1}}
Thank you very much for the explanation, will you make demo instructions about the nonlinear problem and solve with MPC in the future? such as examples at the Algebraic Level, and how to create the correct solver.
Many Thanks
@nate has been working on nonlinear MPC problems. If he has time we can make a demo and maybe a nice interface around his work.
That would be super helpful for a new dolfinx user like me.
Really appreciate you and your team’s work!
Thanks again
Hi Dokken,
Here is the code I wrote to add more than one constraint in the same mesh:
def l2b(li):
return np.array(li, dtype=np.float64).tobytes()
s_m_c = {l2b([0.0279, -0.015132, 0.0]): {l2b([0.0279, 0.0,0.0]): 1}}
s_m_c2 = {l2b([0.0279, 0.015132, 0.0]): {l2b([0.0279, 0.0,0.0]): 1}}
mpc = dolfinx_mpc.MultiPointConstraint(V)
mpc.create_general_constraint(s_m_c)
mpc.create_general_constraint(s_m_c2)
mpc.finalize()
I just want to know is there a better way to write multiple(say 100+) MPCs in one s_m_c instead of assigning each constraint one by one(s_m_c, s_m_c2, s_m_c3, s_m_c4…).
First of all, you can put all the constraints in one dictionary, and only call create_general_constraint once.
You can of course use for loops and append to dictionaries as with any Python program to generate as many constraints as you would like.
I see, many thanks for your explaination!
Sorry to bother again
Based on your instruction: Setting multiple Dirichlet condition
I wrote the code as:
u_R=dolfinx.Function(V)
with u_R.vector.localForm() as loc_bR:
loc_bR.set(uR)
#for p in range(1, len(term_nodes)):
dofs_R1 = dolfinx.fem.locate_dofs_geometrical(V, lambda x: np.isclose(x[0], term_nodes_coord[0]))
bc_R1 = dolfinx.DirichletBC(u_R, dofs_R1)
Is there a way that can define all coordinates ‘x,y,z’ of the node on the mesh, so something like what I wrote in Fenics before:
def u0_boundary(x, on_boundary):
tol = 1e-14;
return abs(x[0] - term_nodes_coord[0][0]) < tol and abs(x[1] - term_nodes_coord[0][1]) < tol and abs(x[2] - term_nodes_coord[0][2]) < tol
bcL = DirichletBC(V, u0, u0_boundary)
I would suggest having a look at my DOLFINx tutorial, as you will find code that illustrates how to do this. Here a slightly modified variant of the code I’ve linked you to
u_D = dolfinx.Function(V)
u_D.interpolate(u_exact)
u_D.x.scatter_forward()
fdim = mesh.topology.dim - 1
def boundary_condition(x):
tol = 1e-14;
return abs(x[0] - term_nodes_coord[0][0]) < tol and abs(x[1] - term_nodes_coord[0][1]) < tol and abs(x[2] - term_nodes_coord[0][2]) < tol
boundary_facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, boundary_condition)
bc = dolfinx.DirichletBC(u_D, dolfinx.fem.locate_dofs_topological(V, fdim, boundary_facets))
Thanks, the BC works now.
So I am still not quite familiar with the code in dolfinx, is there an equivalent command that can find out the value of ‘du/dx’ like what I used in dolfin:
dudx=project(u.dx(0),V)
where u is the solution from the solver and V is function space, where I need the du/dx value only (no need for du/dy or du/dz).
Appreciate all the instructions!
A projection of the function f into V_h is simply solving:
Find u_h\in V_h such that \int_{\Omega} u_h \cdot v_h ~\mathrm{d}x = \int_{\Omega} f \cdot v_h ~\mathrm{d}x \qquad \forall v_h \in V_h.
Therefore, you can achieve a projection with the following:
import dolfinx
from mpi4py import MPI
import ufl
mesh = dolfinx.UnitSquareMesh(MPI.COMM_WORLD, 10, 10)
V = dolfinx.FunctionSpace(mesh, ("CG", 2))
u = dolfinx.Function(V)
u.interpolate(lambda x: 2*x[0]*x[1])
Q = dolfinx.FunctionSpace(mesh, ("DG", 1))
def project(function, space):
p = ufl.TrialFunction(space)
q = ufl.TestFunction(space)
a = ufl.inner(p, q) * ufl.dx
L = ufl.inner(function, q) * ufl.dx
problem = dolfinx.fem.LinearProblem(a, L)
return problem.solve()
udx = project(u.dx(0), Q)
print(udx.vector.array)
I see how it works, Thank you so much
Hi Dokken,
Is there a way in dolfinx can get the direction of norm of grad(u) for 1D problem?
say my code:
V = dolfinx.FunctionSpace(mesh, (“CG”, 1))
u=function(V)
dV = dolfinx.FunctionSpace(mesh, (“DG”, 0))
duds=project((u.dx(0)**2+u.dx(1)**2+u.dx(2)2)(1/2), dV)
for value duds.vector.array[1] for element1 will always be positive
But the value for element1 (u2-u1)/(ds) (where ds=√x^2+y^2+z^2) is negative, u2<u1
So is there a way I can get a correct “du/ds” value in dolfinx?
Many Thanks
One of the problems is that you are projection the norm of the gradient of u into a space it is not in.
Depending on your use-case, you could use dolfinx.Expression to interpolate duds, which should yield better results:
# Copyright (C) 2021 Jørgen S. Dokken
#
# SPDX-License-Identifier: MIT
import ufl
from mpi4py import MPI
import dolfinx
import dolfinx.geometry
import numpy
class Interpolator():
def __init__(self, expr: ufl.core.expr.Expr, u: dolfinx.Function):
i_points = u.function_space.element.interpolation_points()
self.u = u
self._expr = dolfinx.Expression(expr, i_points)
self.value_size = self._expr.value_size
self.dim = self.u.function_space.element.space_dimension() // self.u.function_space.dofmap.index_map_bs
def interpolate(self, cells: numpy.ndarray):
vals = self._expr.eval(cells).reshape(len(cells), self.dim, self.value_size)
for i, cell in enumerate(cells):
blocks = self.u.function_space.dofmap.cell_dofs(cell)
for j, block in enumerate(blocks):
for k in range(self.value_size):
self.u.x.array[block * self.value_size + k] = vals[i, j, k]
mesh = dolfinx.UnitCubeMesh(MPI.COMM_WORLD, 2, 2, 2)
tdim = mesh.topology.dim
num_cells_local = mesh.topology.index_map(tdim).size_local
cells = numpy.arange(num_cells_local, dtype=numpy.int32)
def f(x):
return (0.05 * (x[1] >= 0.75), -0.1 * x[0], 0.2 * x[2] * x[2])
V = dolfinx.VectorFunctionSpace(mesh, ("CG", 1))
u1 = dolfinx.Function(V)
u1.interpolate(f)
u2 = dolfinx.Function(V)
expr = (u1.dx(0)**2 + u1.dx(1)**2 + u1.dx(2)**2)**(1 / 2)
dV = dolfinx.FunctionSpace(mesh, ("DG", 0))
duds = dolfinx.Function(dV)
interpolator = Interpolator(expr, duds)
interpolator.interpolate(cells)
print("Interpolation: ", duds.vector.array)
def project(function, space):
p = ufl.TrialFunction(space)
q = ufl.TestFunction(space)
a = ufl.inner(p, q) * ufl.dx
L = ufl.inner(function, q) * ufl.dx
problem = dolfinx.fem.LinearProblem(a, L)
return problem.solve()
duds_ = project(expr, dV)
print("Projection:", duds_.vector.array)
print("Diff", duds_.vector.array - duds.vector.array)
Thanks for the explanation, but it is still not able to calculate the correct Directional derivative, is there a way can show a duds result with direction such as:
-
(u2-u1)/(ds) (where ds=√x^2+y^2+z^2) if u2 < u1, result will be negative.
-
I was thinking to write something like ufl.grad(u).v where v is the unit vector, and project to element (DG,0) just do not know how to make that happen.
Really appreciate your time!
And my test problem is in 1D so the mesh just lines connect together with x y z coordinates.
As you have not provided a minimal working code example that illustrates what you want to do (and what results you currently get, and explaining what is wrong with them) I cannot give you any further help.
I do not really understand what you are trying to calculate by doing (u2-u1)/ds. If you simply want the direction of the gradient, use my Interpolater class as described above, with the expr: grad(u)/sqrt(u.dx(0)**2+u.dx(1)**2+u.dx(2)**2)
and interpolate it into a suitable space
Sorry for the unclear explanation.
For the code I am testing right now:
- The mesh just in 1D with 4 points(x,y,z coordinates) and 3 lines
V = dolfinx.FunctionSpace(mesh, ("CG", 1))
dV = dolfinx.FunctionSpace(mesh, ("DG", 0))
u = dolfinx.Function(V)
u.vector.array=[1,8,0.5,2]
Gradu=ufl.grad(u)
- what I want is something as:
duds=numpy.dot(ufl.grad(u), unit_vector) #unit_vector is the unit vector of the elements
duds1=project(duds,dV)
I want to make sure
duds1.vector.array[1]=(u.vector.array[2]-u.vector.array[1])/ds #where ds=√x^2+y^2+z^2
As I mentioned (u.dx(0)**2+u.dx(1)**2+u.dx(2)**2)**(1/2) only get the magentiude of grad(u), but not about to present a correct ‘-’ or ‘+’ sign for the value.
Just want to know is that possible for dolfinx to get the unit_vector or just vector for each lines in the mesh?
Many thanks