I would strongly suggest that you:
ufl.inner(ufl.grad(u), element_tangent)element_tangent is the tangent of the element. As you are doing 1D problems, you could easily compute a vector DG 0 function based on the vertices of each element.Thanks, Dokken, I resolve the issue.
I was confused with the topology in dolfinx, but I am about to find the cells use the code:
num_dofs = mesh.topology.index_map(mesh.topology.dim).size_local
cell_list=[]
for i in range(num_dofs):
cell_list.append(V.dofmap.cell_dofs(i))
To check the points id and their cell connections, where based on the point id and point coordinates I can calculate the unit vector for each cell.
Thanks for your advice and help!!!
Dear Dokken,
Thanks for your amazing MPC library
I am just curious whether dolfinx_MPC can achieve constraints such as:
If the dof D located at [d0,d1] should be constrained to the dofs
E and F at [e0,e1] and [f0,f1] as constraints:
D = alpha E + 1 (constant value)
or nonlinear constraints as:
D^2 = E^2 + F^2
Many Thanks
No, it currently does not support such constraints
Hi Dokken,
I realize that the project(u,v) function:
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()
is taking more time to compute in dolfinx than the function in dolfin. This there any way I can speed up the project(u,v) function a little bit?
Thanks for your time, and happy new year.
Please post a minimal code of what you are projecting (i.e. what kind of mesh/cell type), function space and the expression that you are projecting
For instance running the following with the latest dolfinx/dolfinx container returns:
import time
import dolfinx.mesh
import dolfinx.fem
from mpi4py import MPI
import ufl
N = 15
mesh = dolfinx.mesh.create_unit_cube(MPI.COMM_WORLD, N, N, N)
element = ufl.FiniteElement("Lagrange", mesh.ufl_cell(), 1)
V = dolfinx.fem.FunctionSpace(mesh, element)
v = dolfinx.fem.Function(V)
v.x.array[:] = 1
element_2 = ufl.FiniteElement("Lagrange", mesh.ufl_cell(), 2)
Q = dolfinx.fem.FunctionSpace(mesh, element_2)
def project(function: dolfinx.fem.Function, space: dolfinx.fem.FunctionSpace, jit_parameters: dict = {},
petsc_options: dict = {}):
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, jit_parameters=jit_parameters, petsc_options=petsc_options)
return problem.solve()
# --- No JIT options-----
project(v, Q) # Run to JIT compile before timing
start = time.perf_counter()
q = project(v, Q)
end = time.perf_counter()
print(f"No compile options: {end-start}")
# ---- Optimized jit options -----
jit_parameters = {"cffi_extra_compile_args": [
"-O2", "-march=native"], "cffi_libraries": ["m"]}
# Run to JIT compile before timing
q = project(v, Q, jit_parameters=jit_parameters)
start = time.perf_counter()
q = project(v, Q, jit_parameters=jit_parameters)
end = time.perf_counter()
print(f"JIT options: {end-start}")
# Use linear solver (lu/mumps)
petsc_options = {"ksp_type": "preonly",
"pc_type": "lu", "pc_factor_mat_solver_type": "mumps"}
start = time.perf_counter()
q = project(v, Q, jit_parameters=jit_parameters, petsc_options=petsc_options)
end = time.perf_counter()
print(f"Linear solver (mumps): {end-start}")
# -----Interpolation------
q = dolfinx.fem.Function(Q)
start = time.perf_counter()
q.interpolate(v)
end = time.perf_counter()
print(f"Interpolate: {end-start}")
returns
No compile options: 0.17611332399997082
JIT options: 0.17151404399999137
Linear solver (mumps): 1.1754638479999358
Interpolate: 0.002084933999981331
while running the equivalent in the quay.io/fenicsproject/dev:latest dolfin container
from dolfin import *
import time
N = 15
mesh = UnitCubeMesh(N, N, N)
element = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
V = FunctionSpace(mesh, element)
v = Function(V)
element_2 = FiniteElement("Lagrange", mesh.ufl_cell(), 2)
Q = FunctionSpace(mesh, element_2)
q = project(v, Q) # Run to JIT compile before timing
start = time.perf_counter()
q = project(v, Q, solver_type="mumps")
end = time.perf_counter()
print(f"Dolfin (mumps): {end-start}")
start = time.perf_counter()
q = project(v, Q, solver_type="cg")
end = time.perf_counter()
print(f"Dolfin (cg): {end-start}")
you get:
Dolfin (mumps): 1.9629643200000828
Dolfin (cg): 0.16217956299999514
Hi Dokken,
I tried to run this demo code in the docker image that you provided. (docker pull dokken92/dolfinx_mpc:0.1.0)
However, I am getting an error in
problem = dolfinx_mpc.LinearProblem(a, rhs, mpc, bcs=,
petsc_options={“ksp_type”: “preonly”, “pc_type”: “lu”}),
which says module ‘dolfinx_mpc’ has no attribute ‘LinearProblem’.
Do you know what went wrong?
I am trying to the MPC package but currently, I don’t have a working environment or a working example, so I am using the code you posted.
import dolfinx
import dolfinx.io
import dolfinx_mpc
import dolfinx_mpc.utils
import numpy as np
import pytest
import ufl
from mpi4py import MPI
from petsc4py import PETSc
mesh = dolfinx.IntervalMesh(MPI.COMM_WORLD, 10, [np.array([0]), np.array([2])])
V = dolfinx.FunctionSpace(mesh, (“Lagrange”, 1))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
x = ufl.SpatialCoordinate(mesh)
f = 2 * x[0]**2
a = ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx
rhs = ufl.inner(f, v) * ufl.dx
def l2b(li):
return np.array(li, dtype=np.float64).tobytes()
s_m_c = {l2b([0]): {l2b([2]): 2}}
mpc = dolfinx_mpc.MultiPointConstraint(V)
mpc.create_general_constraint(s_m_c)
mpc.finalize()
problem = dolfinx_mpc.LinearProblem(a, rhs, mpc, bcs=,
petsc_options={“ksp_type”: “preonly”, “pc_type”: “lu”})
uh = problem.solve()
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, “u1d.xdmf”, “w”) as xdmf:
xdmf.write_mesh(mesh)
xdmf.write_function(uh)
You need to use the code matching the particular release. The code you posted above uses the main branch of dolfinx_mpc and dolfinx, and can be installed as described in: GitHub - jorgensd/dolfinx_mpc: Extension for dolfinx to handle multi-point constraints.
If you use the docker image, you need to match it with a suitable release, I.e. 0.1.0 matches Release 0.1.0 · jorgensd/dolfinx_mpc · GitHub
0.3.0 matches Release Release compatible with [DOLFINX v0.3.0](https://github.com/FEniCS/do… · jorgensd/dolfinx_mpc · GitHub
The multi point constraint is vary useful! But there is a similar problem about setting boundary constrains. The function in dolfinx_mpc can handle the slip boundary like “u \dot n = 0”. Then how can I set a constrain like “u \dot n = b”, where n is the given vector and b is the given value which may not equal to zero?
I would suggest enforcing it weakly (using a Nitsche-type method). This type of constraint is on my todo, But it is not of high priority atm.
Thank you for your suggestion! I will try the alternative weak form.