Coverage for python/src/dolfinx_mpc/multipointconstraint.py: 83%
177 statements
« prev ^ index » next coverage.py v7.15.2, created at 2026-07-21 19:42 +0000
1# Copyright (C) 2020-2023 Jørgen S. Dokken
2#
3# This file is part of DOLFINX_MPC
4#
5# SPDX-License-Identifier: MIT
6from __future__ import annotations
8from typing import Callable, Dict, List, Optional, Tuple, Union
10from petsc4py import PETSc as _PETSc
12import dolfinx.cpp as _cpp
13import dolfinx.fem as _fem
14import dolfinx.mesh as _mesh
15import numpy
16import numpy.typing as npt
17from dolfinx import default_real_type, default_scalar_type
19import dolfinx_mpc.cpp
21from .dictcondition import create_dictionary_constraint
23_mpc_classes = Union[
24 dolfinx_mpc.cpp.mpc.MultiPointConstraint_double,
25 dolfinx_mpc.cpp.mpc.MultiPointConstraint_float,
26 dolfinx_mpc.cpp.mpc.MultiPointConstraint_complex_double,
27 dolfinx_mpc.cpp.mpc.MultiPointConstraint_complex_float,
28]
29_float_classes = Union[numpy.float32, numpy.float64, numpy.complex128, numpy.complex64]
30_float_array_types = Union[
31 npt.NDArray[numpy.float32],
32 npt.NDArray[numpy.float64],
33 npt.NDArray[numpy.complex64],
34 npt.NDArray[numpy.complex128],
35]
36_mpc_data_classes = Union[
37 dolfinx_mpc.cpp.mpc.mpc_data_double,
38 dolfinx_mpc.cpp.mpc.mpc_data_float,
39 dolfinx_mpc.cpp.mpc.mpc_data_complex_double,
40 dolfinx_mpc.cpp.mpc.mpc_data_complex_float,
41]
44class MPCData:
45 _cpp_object: _mpc_data_classes
47 def __init__(
48 self,
49 slaves: npt.NDArray[numpy.int32],
50 masters: npt.NDArray[numpy.int64],
51 coeffs: _float_array_types,
52 owners: npt.NDArray[numpy.int32],
53 offsets: npt.NDArray[numpy.int32],
54 ):
55 if coeffs.dtype.type == numpy.float32:
56 self._cpp_object = dolfinx_mpc.cpp.mpc.mpc_data_float(slaves, masters, coeffs, owners, offsets)
57 elif coeffs.dtype.type == numpy.float64:
58 self._cpp_object = dolfinx_mpc.cpp.mpc.mpc_data_double(slaves, masters, coeffs, owners, offsets)
59 elif coeffs.dtype.type == numpy.complex64:
60 self._cpp_object = dolfinx_mpc.cpp.mpc.mpc_data_complex_float(slaves, masters, coeffs, owners, offsets)
61 elif coeffs.dtype.type == numpy.complex128:
62 self._cpp_object = dolfinx_mpc.cpp.mpc.mpc_data_complex_double(slaves, masters, coeffs, owners, offsets)
63 else:
64 raise ValueError("Unsupported dtype {coeffs.dtype.type} for coefficients")
66 @property
67 def slaves(self):
68 return self._cpp_object.slaves
70 @property
71 def masters(self):
72 return self._cpp_object.masters
74 @property
75 def coeffs(self):
76 return self._cpp_object.coeffs
78 @property
79 def owners(self):
80 return self._cpp_object.owners
82 @property
83 def offsets(self):
84 return self._cpp_object.offsets
87class MultiPointConstraint:
88 """
89 Hold data for multi point constraint relation ships,
90 including new index maps for local assembly of matrices and vectors.
92 Args:
93 V: The function space
94 dtype: The dtype of the underlying functions
95 """
97 _slaves: npt.NDArray[numpy.int32]
98 _masters: npt.NDArray[numpy.int64]
99 _coeffs: _float_array_types
100 _owners: npt.NDArray[numpy.int32]
101 _offsets: npt.NDArray[numpy.int32]
102 V: _fem.FunctionSpace
103 finalized: bool
104 _cpp_object: _mpc_classes
105 _dtype: npt.DTypeLike
106 __slots__ = tuple(__annotations__)
108 def __init__(self, V: _fem.FunctionSpace, dtype: npt.DTypeLike = default_scalar_type):
109 self._slaves = numpy.array([], dtype=numpy.int32)
110 self._masters = numpy.array([], dtype=numpy.int64)
111 self._coeffs = numpy.array([], dtype=dtype) # type: ignore
112 self._owners = numpy.array([], dtype=numpy.int32)
113 self._offsets = numpy.array([0], dtype=numpy.int32)
114 self.V = V
115 self.finalized = False
116 self._dtype = dtype
118 def add_constraint(
119 self,
120 V: _fem.FunctionSpace,
121 slaves: npt.NDArray[numpy.int32],
122 masters: npt.NDArray[numpy.int64],
123 coeffs: _float_array_types,
124 owners: npt.NDArray[numpy.int32],
125 offsets: npt.NDArray[numpy.int32],
126 ):
127 """
128 Add new constraint given by numpy arrays.
130 Args:
131 V: The function space for the constraint
132 slaves: List of all slave dofs (using local dof numbering) on this process
133 masters: List of all master dofs (using global dof numbering) on this process
134 coeffs: The coefficients corresponding to each master.
135 owners: The process each master is owned by.
136 offsets: Array indicating the location in the masters array for the i-th slave
137 in the slaves arrays, i.e.
139 .. highlight:: python
140 .. code-block:: python
142 masters_of_owned_slave[i] = masters[offsets[i]:offsets[i+1]]
144 """
145 assert V == self.V
146 self._already_finalized()
148 if len(slaves) > 0:
149 self._offsets = numpy.append(self._offsets, offsets[1:] + len(self._masters))
150 self._slaves = numpy.append(self._slaves, slaves)
151 self._masters = numpy.append(self._masters, masters)
152 self._coeffs = numpy.array(numpy.append(self._coeffs, coeffs), dtype=self._dtype)
153 self._owners = numpy.append(self._owners, owners)
155 def add_constraint_from_mpc_data(self, V: _fem.FunctionSpace, mpc_data: Union[_mpc_data_classes, MPCData]):
156 """
157 Add new constraint given by an `dolfinc_mpc.cpp.mpc.mpc_data`-object
158 """
159 self._already_finalized()
160 self.add_constraint(
161 V,
162 mpc_data.slaves,
163 mpc_data.masters,
164 mpc_data.coeffs,
165 mpc_data.owners,
166 mpc_data.offsets,
167 )
169 def finalize(self) -> None:
170 """
171 Finializes the multi point constraint. After this function is called, no new constraints can be added
172 to the constraint. This function creates a map from the cells (local to index) to the slave degrees of
173 freedom and builds a new index map and function space where unghosted master dofs are added as ghosts.
174 """
175 self._already_finalized()
176 self._coeffs.astype(numpy.dtype(self._dtype))
177 # Initialize C++ object and create slave->cell maps
178 if self._dtype == numpy.float32:
179 self._cpp_object = dolfinx_mpc.cpp.mpc.MultiPointConstraint_float(
180 self.V._cpp_object,
181 self._slaves,
182 self._masters,
183 self._coeffs.astype(self._dtype),
184 self._owners,
185 self._offsets,
186 )
187 elif self._dtype == numpy.float64:
188 self._cpp_object = dolfinx_mpc.cpp.mpc.MultiPointConstraint_double(
189 self.V._cpp_object,
190 self._slaves,
191 self._masters,
192 self._coeffs.astype(self._dtype),
193 self._owners,
194 self._offsets,
195 )
196 elif self._dtype == numpy.complex64:
197 self._cpp_object = dolfinx_mpc.cpp.mpc.MultiPointConstraint_complex_float(
198 self.V._cpp_object,
199 self._slaves,
200 self._masters,
201 self._coeffs.astype(self._dtype),
202 self._owners,
203 self._offsets,
204 )
206 elif self._dtype == numpy.complex128:
207 self._cpp_object = dolfinx_mpc.cpp.mpc.MultiPointConstraint_complex_double(
208 self.V._cpp_object,
209 self._slaves,
210 self._masters,
211 self._coeffs.astype(self._dtype),
212 self._owners,
213 self._offsets,
214 )
215 else:
216 raise ValueError("Unsupported dtype {coeffs.dtype.type} for coefficients")
218 # Replace function space
219 self.V = _fem.FunctionSpace(self.V.mesh, self.V.ufl_element(), self._cpp_object.function_space)
221 self.finalized = True
222 # Delete variables that are no longer required
223 del (self._slaves, self._masters, self._coeffs, self._owners, self._offsets)
225 def create_periodic_constraint_topological(
226 self,
227 V: _fem.FunctionSpace,
228 meshtag: _mesh.MeshTags,
229 tag: int,
230 relation: Callable[[numpy.ndarray], numpy.ndarray],
231 bcs: List[_fem.DirichletBC],
232 scale: _float_classes = default_scalar_type(1.0), # type: ignore
233 tol: _float_classes = 500 * numpy.finfo(default_real_type).eps,
234 num_threads: Optional[int] = 1,
235 ):
236 """
237 Create periodic condition for all closure dofs of on all entities in `meshtag` with value `tag`.
238 :math:`u(x_i) = scale * u(relation(x_i))` for all of :math:`x_i` on marked entities.
240 Args:
241 V: The function space to assign the condition to. Should either be the space of the MPC or a sub space.
242 meshtag: MeshTag for entity to apply the periodic condition on
243 tag: Tag indicating which entities should be slaves
244 relation: Lambda-function describing the geometrical relation
245 bcs: Dirichlet boundary conditions for the problem (Periodic constraints will be ignored for these dofs)
246 scale: Float for scaling bc
247 tol: Tolerance for adding scaled basis values to MPC. Any contribution that is less than this value
248 is ignored. The tolerance is also added as padding for the bounding box trees and corresponding
249 collision searches to determine periodic degrees of freedom.
250 num_threads: The number of threads to use for certain operations
251 """
252 bcs_ = [bc._cpp_object for bc in bcs]
253 if isinstance(scale, numpy.generic): # nanobind conversion of numpy dtypes to general Python types
254 scale = scale.item() # type: ignore
255 if V is self.V:
256 mpc_data = dolfinx_mpc.cpp.mpc.create_periodic_constraint_topological(
257 self.V._cpp_object,
258 meshtag._cpp_object,
259 tag,
260 relation,
261 bcs_,
262 scale,
263 False,
264 float(tol),
265 num_threads=num_threads,
266 )
267 elif self.V.contains(V):
268 mpc_data = dolfinx_mpc.cpp.mpc.create_periodic_constraint_topological(
269 V._cpp_object,
270 meshtag._cpp_object,
271 tag,
272 relation,
273 bcs_,
274 scale,
275 True,
276 float(tol),
277 num_threads=num_threads,
278 )
279 else:
280 raise RuntimeError("The input space has to be a sub space (or the full space) of the MPC")
281 self.add_constraint_from_mpc_data(self.V, mpc_data=mpc_data)
283 def create_periodic_constraint_geometrical(
284 self,
285 V: _fem.FunctionSpace,
286 indicator: Callable[[numpy.ndarray], numpy.ndarray],
287 relation: Callable[[numpy.ndarray], numpy.ndarray],
288 bcs: List[_fem.DirichletBC],
289 scale: _float_classes = default_scalar_type(1.0), # type: ignore
290 tol: _float_classes = 500 * numpy.finfo(default_real_type).eps,
291 num_threads: Optional[int] = 1,
292 ):
293 """
294 Create a periodic condition for all degrees of freedom whose physical location satisfies
295 :math:`indicator(x_i)==True`, i.e.
296 :math:`u(x_i) = scale * u(relation(x_i))` for all :math:`x_i`
298 Args:
299 V: The function space to assign the condition to. Should either be the space of the MPC or a sub space.
300 indicator: Lambda-function to locate degrees of freedom that should be slaves
301 relation: Lambda-function describing the geometrical relation to master dofs
302 bcs: Dirichlet boundary conditions for the problem
303 (Periodic constraints will be ignored for these dofs)
304 scale: Float for scaling bc
305 tol: Tolerance for adding scaled basis values to MPC. Any contribution that is less than this value
306 is ignored. The tolerance is also added as padding for the bounding box trees and corresponding
307 collision searches to determine periodic degrees of freedom.
308 num_threads: The number of threads to use for certain operations.
309 """
310 if isinstance(scale, numpy.generic): # nanobind conversion of numpy dtypes to general Python types
311 scale = scale.item() # type: ignore
312 bcs = [] if bcs is None else [bc._cpp_object for bc in bcs]
313 if V is self.V:
314 mpc_data = dolfinx_mpc.cpp.mpc.create_periodic_constraint_geometrical(
315 self.V._cpp_object, indicator, relation, bcs, scale, False, float(tol), num_threads
316 )
317 elif self.V.contains(V):
318 mpc_data = dolfinx_mpc.cpp.mpc.create_periodic_constraint_geometrical(
319 V._cpp_object, indicator, relation, bcs, scale, True, float(tol), num_threads
320 )
321 else:
322 raise RuntimeError("The input space has to be a sub space (or the full space) of the MPC")
323 self.add_constraint_from_mpc_data(self.V, mpc_data=mpc_data)
325 def create_slip_constraint(
326 self,
327 space: _fem.FunctionSpace,
328 facet_marker: Tuple[_mesh.MeshTags, int],
329 v: _fem.Function,
330 bcs: List[_fem.DirichletBC] = [],
331 ):
332 """
333 Create a slip constraint :math:`u \\cdot v=0` over the entities defined in `facet_marker` with the given index.
335 Args:
336 space: Function space (possible sub space) for the current constraint
337 facet_marker: Tuple containomg the mesh tag and marker used to locate degrees of freedom
338 v: Function containing the directional vector to dot your slip condition (most commonly a normal vector)
339 bcs: List of Dirichlet BCs (slip conditions will be ignored on these dofs)
341 Examples:
342 Create constaint :math:`u\\cdot n=0` of all indices in `mt` marked with `i`
344 .. highlight:: python
345 .. code-block:: python
347 V = dolfinx.fem.functionspace(mesh, ("CG", 1))
348 mpc = MultiPointConstraint(V)
349 n = dolfinx.fem.Function(V)
350 mpc.create_slip_constaint(V, (mt, i), n)
352 Create slip constaint for a mixed function space:
354 .. highlight:: python
355 .. code-block:: python
357 cellname = mesh.basix_cell()
358 Ve = basix.ufl.element(basix.ElementFamily.P, cellname , 2, shape=(mesh.geometry.dim,))
359 Qe = basix.ufl.element(basix.ElementFamily.P, cellname , 1)
360 me = basix.ufl.mixed_element([Ve, Qe])
361 W = dolfinx.fem.functionspace(mesh, me)
362 mpc = MultiPointConstraint(W)
363 n_space, _ = W.sub(0).collapse()
364 normal = dolfinx.fem.Function(n_space)
365 mpc.create_slip_constraint(W.sub(0), (mt, i), normal, bcs=[])
367 A slip condition cannot be applied on the same degrees of freedom as a Dirichlet BC, and therefore
368 any Dirichlet bc for the space of the multi point constraint should be supplied.
370 .. highlight:: python
371 .. code-block:: python
373 cellname = mesh.basix_cell()
374 Ve = basix.ufl.element(basix.ElementFamily.P, cellname , 2, shape=(mesh.geometry.dim,))
375 Qe = basix.ufl.element(basix.ElementFamily.P, cellname , 1)
376 me = basix.ufl.mixed_element([Ve, Qe])
377 W = dolfinx.fem.functionspace(mesh, me)
378 mpc = MultiPointConstraint(W)
379 n_space, _ = W.sub(0).collapse()
380 normal = Function(n_space)
381 bc = dolfinx.fem.dirichletbc(inlet_velocity, dofs, W.sub(0))
382 mpc.create_slip_constraint(W.sub(0), (mt, i), normal, bcs=[bc])
383 """
384 bcs = [] if bcs is None else [bc._cpp_object for bc in bcs]
385 if space is self.V:
386 sub_space = False
387 elif self.V.contains(space):
388 sub_space = True
389 else:
390 raise ValueError("Input space has to be a sub space of the MPC space")
391 mpc_data = dolfinx_mpc.cpp.mpc.create_slip_condition(
392 space._cpp_object,
393 facet_marker[0]._cpp_object,
394 facet_marker[1],
395 v._cpp_object,
396 bcs,
397 sub_space,
398 )
399 self.add_constraint_from_mpc_data(self.V, mpc_data=mpc_data)
401 def create_general_constraint(
402 self,
403 slave_master_dict: Dict[bytes, Dict[bytes, float]],
404 subspace_slave: Optional[int] = None,
405 subspace_master: Optional[int] = None,
406 ):
407 """
408 Args:
409 V: The function space
410 slave_master_dict: Nested dictionary, where the first key is the bit representing the slave dof's
411 coordinate in the mesh. The item of this key is a dictionary, where each key of this dictionary
412 is the bit representation of the master dof's coordinate, and the item the coefficient for
413 the MPC equation.
414 subspace_slave: If using mixed or vector space, and only want to use dofs from a sub space
415 as slave add index here
416 subspace_master: Subspace index for mixed or vector spaces
418 Example:
419 If the dof `D` located at `[d0, d1]` should be constrained to the dofs
420 `E` and `F` at `[e0, e1]` and `[f0, f1]` as :math:`D = \\alpha E + \\beta F`
421 the dictionary should be:
423 .. highlight:: python
424 .. code-block:: python
426 {numpy.array([d0, d1], dtype=mesh.geometry.x.dtype).tobytes():
427 {numpy.array([e0, e1], dtype=mesh.geometry.x.dtype).tobytes(): alpha,
428 numpy.array([f0, f1], dtype=mesh.geometry.x.dtype).tobytes(): beta}}
429 """
430 slaves, masters, coeffs, owners, offsets = create_dictionary_constraint(
431 self.V, slave_master_dict, subspace_slave, subspace_master
432 )
433 self.add_constraint(self.V, slaves, masters, coeffs, owners, offsets)
435 def create_contact_slip_condition(
436 self,
437 meshtags: _mesh.MeshTags,
438 slave_marker: int,
439 master_marker: int,
440 normal: _fem.Function,
441 eps2: float = 1e-20,
442 num_threads: Optional[int] = 1,
443 ):
444 """
445 Create a slip condition between two sets of facets marker with individual markers.
446 The interfaces should be within machine precision of eachother, but the vertices does not need to align.
447 The condition created is :math:`u_s \\cdot normal_s = u_m \\cdot normal_m` where `s` is the
448 restriction to the slave facets, `m` to the master facets.
450 Args:
451 meshtags: The meshtags of the set of facets to tie together
452 slave_marker: The marker of the slave facets
453 master_marker: The marker of the master facets
454 normal: The function used in the dot-product of the constraint
455 eps2: The tolerance for the squared distance between cells to be considered as a collision
456 num_threads: The number of threads to use for certain operations
457 """
458 if isinstance(eps2, numpy.generic): # nanobind conversion of numpy dtypes to general Python types
459 eps2 = eps2.item() # type: ignore
460 mpc_data = dolfinx_mpc.cpp.mpc.create_contact_slip_condition(
461 self.V._cpp_object, meshtags._cpp_object, slave_marker, master_marker, normal._cpp_object, eps2, num_threads
462 )
463 self.add_constraint_from_mpc_data(self.V, mpc_data)
465 def create_contact_inelastic_condition(
466 self,
467 meshtags: _cpp.mesh.MeshTags_int32,
468 slave_marker: int,
469 master_marker: int,
470 eps2: float = 1e-20,
471 allow_missing_masters: bool = False,
472 num_threads: Optional[int] = 1,
473 ):
474 """
475 Create a contact inelastic condition between two sets of facets marker with individual markers.
476 The interfaces should be within machine precision of eachother, but the vertices does not need to align.
477 The condition created is :math:`u_s = u_m` where `s` is the restriction to the
478 slave facets, `m` to the master facets.
480 Args:
481 meshtags: The meshtags of the set of facets to tie together
482 slave_marker: The marker of the slave facets
483 master_marker: The marker of the master facets
484 eps2: The tolerance for the squared distance between cells to be considered as a collision
485 allow_missing_masters: If true, the function will not throw an error if a degree of freedom
486 in the closure of the master entities does not have a corresponding set of slave degree
487 of freedom.
488 num_threads: The number of threads to use for certain operations
489 """
490 if isinstance(eps2, numpy.generic): # nanobind conversion of numpy dtypes to general Python types
491 eps2 = eps2.item() # type: ignore
492 mpc_data = dolfinx_mpc.cpp.mpc.create_contact_inelastic_condition(
493 self.V._cpp_object,
494 meshtags._cpp_object,
495 slave_marker,
496 master_marker,
497 eps2,
498 allow_missing_masters,
499 num_threads,
500 )
501 self.add_constraint_from_mpc_data(self.V, mpc_data)
503 @property
504 def is_slave(self) -> numpy.ndarray:
505 """
506 Returns a vector of integers where the ith entry indicates if a degree of freedom (local to process) is a slave.
507 """
508 self._not_finalized()
509 return self._cpp_object.is_slave
511 @property
512 def slaves(self):
513 """
514 Returns the degrees of freedom for all slaves local to process
515 """
516 self._not_finalized()
517 return self._cpp_object.slaves
519 @property
520 def masters(self) -> _cpp.graph.AdjacencyList_int32:
521 """
522 Returns an adjacency-list whose ith node corresponds to
523 a degree of freedom (local to process), and links the corresponding master dofs (local to process).
525 Examples:
527 .. highlight:: python
528 .. code-block:: python
530 masters = mpc.masters
531 masters_of_dof_i = masters.links(i)
532 """
533 self._not_finalized()
534 return self._cpp_object.masters
536 def coefficients(self) -> _float_array_types:
537 """
538 Returns a vector containing the coefficients for the constraint, and the corresponding offsets
539 for the ith degree of freedom.
541 Examples:
543 .. highlight:: python
544 .. code-block:: python
546 coeffs, offsets = mpc.coefficients()
547 coeffs_of_slave_i = coeffs[offsets[i]:offsets[i+1]]
548 """
549 self._not_finalized()
550 return self._cpp_object.coefficients()
552 @property
553 def num_local_slaves(self):
554 """
555 Return the number of slaves owned by the current process.
556 """
557 self._not_finalized()
558 return self._cpp_object.num_local_slaves
560 @property
561 def cell_to_slaves(self):
562 """
563 Returns an `dolfinx.cpp.graph.AdjacencyList_int32` whose ith node corresponds to
564 the ith cell (local to process), and links the corresponding slave degrees of
565 freedom in the cell (local to process).
567 Examples:
569 .. highlight:: python
570 .. code-block:: python
572 cell_to_slaves = mpc.cell_to_slaves()
573 slaves_in_cell_i = cell_to_slaves.links(i)
574 """
575 self._not_finalized()
576 return self._cpp_object.cell_to_slaves
578 @property
579 def function_space(self):
580 """
581 Return the function space for the multi-point constraint with the updated index map
582 """
583 self._not_finalized()
584 return self.V
586 def backsubstitution(self, u: Union[_fem.Function, _PETSc.Vec]) -> None: # type: ignore
587 """
588 For a Function, impose the multi-point constraint by backsubstiution.
589 This function is used after solving the reduced problem to obtain the values
590 at the slave degrees of freedom
592 .. note::
593 It is the users responsibility to destroy the PETSc vector
595 Args:
596 u: The input function
597 """
598 try:
599 self._cpp_object.backsubstitution(u.x.array) # type: ignore
600 assert isinstance(u, _fem.Function)
601 u.x.scatter_forward()
602 except AttributeError:
603 assert isinstance(u, _PETSc.Vec)
604 with u.localForm() as vector_local:
605 self._cpp_object.backsubstitution(vector_local.array_w)
606 u.ghostUpdate(addv=_PETSc.InsertMode.INSERT, mode=_PETSc.ScatterMode.FORWARD) # type: ignore
608 def homogenize(self, u: _fem.Function) -> None:
609 """
610 For a vector, homogenize (set to zero) the vector components at the multi-point
611 constraint slave DoF indices. This is particularly useful for nonlinear problems.
613 Args:
614 u: The input vector
615 """
616 self._cpp_object.homogenize(u.x.array)
617 u.x.scatter_forward()
619 def _already_finalized(self):
620 """
621 Check if we have already finalized the multi point constraint
622 """
623 if self.finalized:
624 raise RuntimeError("MultiPointConstraint has already been finalized")
626 def _not_finalized(self):
627 """
628 Check if we have finalized the multi point constraint
629 """
630 if not self.finalized:
631 raise RuntimeError("MultiPointConstraint has not been finalized")